Genome sequencing of three medaka strains

We extracted genomic DNA from three medaka strains (Hd-rR, HNI, and HSOK) and sequenced it using PacBio HiFi and ONT ultralong (ONT-UL) reads, which have complementary strengths. HiFi provides ∼10× lower error rates (∼0.1%) compared with ONT-UL (∼1%), whereas ONT-UL produces ∼10× longer reads (N10 ∼150 kb vs. ∼20 kb). HiFi libraries were prepared from whole-body samples, whereas ONT-UL libraries were prepared from pooled heart, spleen, and testis tissues to obtain DNA of sufficient length and yield. To facilitate the discrimination of sequences divergent between Chr X and Chr Y (which are also known as Chr 1) in medaka, HiFi reads were sequenced from male and female samples separately for each strain. The combined sequencing depth of HiFi reads was 68.4×–75.9×, and that of ONT-UL reads of size >100 kb amounted to 19.7×–32.1×. Additionally, we generated Hi-C long-range data for assembly validation and curation of Hd-rR and HNI. All sequencing statistics are summarized in Supplemental Table S1.

Genome assembly and gene annotation

For each strain, the initial set of chromosomal scaffolds was generated using Hifiasm (Cheng et al. 2024) and RagTag (Alonge et al. 2022). After manually correcting false inversions using Hi-C contact maps (Supplemental Fig. S1), we performed semiautomated curation of the assemblies to obtain a gapless, T2T complete genome for Hd-rR and nearly complete genomes for HNI and HSOK (Methods), which we refer to as the v3 assemblies to distinguish them from the v2 assemblies reported in 2017 (Ichikawa et al. 2017). In particular, we corrected large tandem repeats such as centromeric satellite arrays, whose sequences are prone to collapse during assembly, and extended telomeric sequences at both chromosomal ends (Supplemental Fig. S2). Additionally, sequences of sex-differentiating regions on Chr 1 and the mitochondrial genomes were identified during the curation process.

After the curation, the base-level quality value (QV) of the Hd-rR assembly improved from 49.1 to 57.6 (Fig. 1B; Supplemental Table S2), whereas the QV was limited to 36.3 in the previous v2 assembly. For HNI and HSOK, their assembly graphs were highly tangled (Supplemental Fig. S3), indicating that some repeats remain unresolved, and the assemblies for these two strains remained incomplete. Nevertheless, we anchored 98.8%/97.5% of the sequences to chromosomes with only 15/83 gaps in HNI/HSOK, compared with about 1000 gaps per strain in the v2 assemblies. The comparatively lower assembly contiguity in HSOK may be attributable to the relatively smaller amount of ONT-UL reads (Supplemental Table S1) and the higher copy number of the large transposon Teratorn (see below). In 45S and 5S ribosomal DNA arrays and adjacent megabase-scale tandem repeats, we observed sequence variation and read-depth fluctuations that may reflect sample variation, somatic variation, and/or residual uncertainty (Fig. 1C; Supplemental Figs. S4, S5; Supplemental Discussion). Overall, the present v3 assemblies represent a substantial improvement on reconstructing repetitive elements over the v2 assemblies (Fig. 1C; Supplemental Fig. S6).

We annotated protein-coding genes by integrating multiple gene prediction methods and external evidence (Supplemental Fig. S7). Specifically, gene models of the previous v2 Hd-rR assembly were lifted over to the v3 genomes and merged with de novo predictions with RNA-seq providing homolog hints. In total, 36,040 genes were determined for Hd-rR, and 29,032 and 29,668 genes were predicted for HNI and HSOK, respectively (Fig. 1B; Supplemental Table S3). The higher gene count in Hd-rR than in HNI and HSOK may reflect the larger amount of RNA-seq data available for Hd-rR (387 Gb compared with 69 Gb) and the resulting more comprehensive detection of lowly expressed transcripts. The BUSCO completeness (Manni et al. 2021) ranged from 99.6%–99.8% across the three strains, supporting the high quality of our medaka gene annotations.

Centromere localization and chromosome orientation

In the prior v2 medaka assemblies, only ∼10% of centromeric sequences were inferred, whereas the v3 assemblies resolved all 24 centromeres in Hd-rR and yielded 21 and 14 gapless centromeres in HNI and HSOK, respectively. Compared with v2, centromere positions were newly determined or corrected for 42 chromosomes (Supplemental Table S4). The updated centromeric satellite locations were largely shared across the three strains and consistent with previous FISH images in Hd-rR (Ichikawa et al. 2017). Specifically, of the 24 chromosomes, 12 were acrocentric whose centromeres are located near telomere, and the other 12 were nonacrocentric, except that five HNI chromosomes (Chr 5, 6, 7, 10, 16) contained two distant centromeric satellite arrays at submetacentric and acrocentric positions, respectively. These putative pseudodicentric chromosomes were likely generated by large inversions, which were supported by multiple ONT-UL reads (Supplemental Figs. S8, S9; Supplemental Discussion).

Precise identification of centromeric locations enabled p/q-arm assignment and standardization of chromosome orientation across all strains by defining position 1 at the telomeric end of the p-arm. Compared with v2, 12 chromosomes including the sex chromosomes were reverse-complemented in each v3 assembly.

Comparative analysis of centromeric satellite monomers

Leveraging fully assembled centromeric sequences, for each chromosome of the three strains we clustered ∼160 nt centromeric satellite monomers into representative monomers based on sequence similarity. As similarly defined in primates (Alexandrov et al. 2001), we further categorized these chromosomal representative monomers into suprachromosomal families (SFs), in which each SF is a group of monomers that share sequence similarity and span multiple chromosomes. We newly identified five SFs (SF1–SF5), whereas only four SFs were inferred from the fragmented v2 assemblies (Supplemental Fig. S10A–D).

In the phylogenetic tree of the 350 representative monomers belonging to SF1–SF5 (Fig. 2A; Supplemental Data File 1), monomers from different strains were intermixed, indicating that all three strains share the same set of SFs. We then tested whether representative monomer sequence similarity varies with strain relationship (intra- vs. interstrain) and chromosome type (acrocentric vs. nonacrocentric) and detected a strong effect of chromosome type on sequence divergence (Freedman–Lane permutation ANOVA, 20,000 permutations, P = 5.0 × 10⁻⁵), with no significant effect of strain relationship (P = 0.190) and no interaction (P = 0.683). Consistently, pairwise comparisons confirmed that monomers were significantly less diverged within acrocentric chromosomes than within nonacrocentric ones (two-sided Wilcoxon test, Bonferroni-corrected for all possible comparisons, P < 1.0 × 10⁻⁵⁰) (Fig. 2B) as previously suggested (Ichikawa et al. 2017). We further identified 20 groups of highly similar (>97%) representative monomers shared between different acrocentric chromosomes in a genome, which we refer to as interchromosomal ultraconserved (ICUC) monomers (Fig. 2A,B). All ICUC monomers were found within individual strains and not across strains, as well as exclusively among acrocentric chromosomes (Supplemental Fig. S10E). These results support preferential sequence homogenization among acrocentric chromosomes, likely through interchromosomal gene conversions, as previously suggested in humans (Choo et al. 1989; Guarracino et al. 2023).

Figure 2.

Compositional and phylogenetic landscape of centromeric satellites. (A) Phylogenetic tree of 350 representative monomers excluding divergent monomers. The outermost track is the monomer count per representative, and the “hypo” track indicates ONT-inferred CpG hypomethylation level (i.e., the proportion of hypomethylated monomers within each representative). Dark gray shades on the dotted lines connecting the tree leaves and the strain track indicate interchromosomal ultraconserved (ICUC) monomers. (B) Pairwise interchromosomal monomer sequence similarities for four categories by two aspects: intra- versus interstrain, and within acrocentric versus nonacrocentric chromosomes. Wilcoxon test; (*) P < 1.0 × 10⁻⁵⁰, (n.s.) P > 0.01. (C) Centromeric satellite array lengths by chromosome and SF group. Acrocentric chromosomes are placed on the upper side and nonacrocentric chromosomes on the lower side.

Open in new tabDownload PowerPointLink to figure

The composition of SFs as a set was generally conserved for each chromosome across the three strains, although individual array sizes varied (e.g., from 190 kb in Chr 22 to 1.4 Mb in Chr 4 in Hd-rR) (Fig. 2C), probably owing to unequal crossover and replication slippage during meiosis (Smith 1976; Levinson and Gutman 1987). Although nonacrocentric chromosomes showed greater variability in SF composition, acrocentric chromosomes were consistently dominated by SF2 monomers with minor fractions of SF1 and SF3.

Positional distributions and methylation patterns in centromeric satellite arrays

Despite their low abundance in acrocentric chromosomes (Fig. 2C), SF1 and SF3 monomers displayed a unique positional and putative epigenomic organization distinct from other reported eukaryotic centromeres (Naish et al. 2021; Altemose et al. 2022; Chen et al. 2023; Chaudhry et al. 2025). Specifically, in the centromeric sequences of acrocentric chromosomes in Hd-rR, short composite arrays containing both SF1 and SF3 (hereinafter referred to as SF1+3 arrays) were consistently interspersed at multiple positions within large SF2 arrays, except on Chr 22 (Fig. 3A). In line with this structural organization, these SF1+3 arrays corresponded predominantly to hypomethylated regions inferred by ONT-UL reads, whereas the large SF2 arrays were almost entirely hypermethylated (Fig. 3A). The assembly of these arrays was validated using ONT-UL reads (Supplemental Fig. S11). In contrast, nonacrocentric centromeres were largely homogenized by a single SF or a composite of SF1+3, supporting independent and rapid concerted evolution of satellite arrays (Fig. 3B). These observations were common in HNI and HSOK as well (Supplemental Figs. S12, S13).

Figure 3.

Positional landscape of sequence and CpG methylation status of centromeric satellites in the complete Hd-rR genome. (A,B) Positional distributions of centromeric monomers colored by SF group, along with average ONT-inferred CpG methylation levels from zero to 100, for acrocentric (A) and nonacrocentric (B) chromosomes of the T2T Hd-rR assembly. Thick lines indicate ICUC monomers. For HNI and HSOK, see Supplemental Figures S12 and S13. (C) Graph representations of representative monomers for SF1 and SF3, for which monomer sequences were progressively aligned and sequence variations appear as branching paths, shown around the sequence motifs enriched in putatively hypomethylated monomers. Edges are colored by relative enrichment in hypo- and hypermethylated monomers.

Open in new tabDownload PowerPointLink to figure

We detected putative CpG hypomethylated regions in centromeric sequences, termed centromere dip regions (CDRs), although future experimental validation will be required. Each chromosome had at least one CDR, whereas some nonacrocentric chromosomes exhibited relatively lower dip levels (Fig. 3A,B). Of note, in acrocentric chromosomes, 93.1% of monomers in the interspersed SF1+3 arrays were hypomethylated at CpG sites, whereas 99.7% of monomers in the large SF2 arrays were hypermethylated (Figs. 2A, 3A; Supplemental Fig. S14). Among the acrocentric SF1 and SF3 monomers, we further identified sequence motifs enriched in hypomethylated monomers (Fig. 3C; Supplemental Fig. S15; Supplemental Table S5). Specifically, GCATATTGTACATAAAAA (18 nt) was the significant motif enriched in hypomethylated SF1 monomers (86.8%, n = 81, χ² test, Bonferroni-corrected, P = 4.8 × 10⁻⁶). In SF3, ATAGCGT (7 nt) containing a CpG site was the significant motif, although its enrichment was smaller than that of SF1 (62.9%, n = 72, χ² test, Bonferroni-corrected, P = 3.9 × 10⁻⁴). These results suggest that putative CpG hypomethylated regions in acrocentric chromosomes are characterized by SF1+SF3 composite arrays within SF2 arrays and by distinct sequence motifs.

Structural rearrangements among three medaka strains

We detected large-scale structural variation between the three nearly complete genomes (Fig. 4). In addition to the previously reported large (>10 Mb) inversion on Chr 11 in Hd-rR, we identified large subtelomeric inversions in HNI on five chromosomes causing division of centromeric satellite arrays. We also identified large (∼20 Mb) highly nonsyntenic regions in the distal part of Chr 2 (Fig. 4; Supplemental Fig. S8), most of which were unassembled in previous assemblies.

Figure 4.

Structural rearrangements and transposons. Structural rearrangements among the three medaka strains with annotations of major repetitive elements. Strains are ordered as HNI, Hd-rR, and HSOK to clearly illustrate that megabase-scale inversions associated with divided centromeric arrays are specific to HNI. For Chr 11, however, the order is Hd-rR, HNI, and HSOK to highlight a large inversion specific to Hd-rR.

Open in new tabDownload PowerPointLink to figure

Pairwise comparisons identified 21,495 SVs (insertions, deletions, and inversions >50 bp) between Hd-rR and HNI compared with 18,713 SVs between HNI and HSOK and 18,229 SVs between Hd-rR and HSOK. These SV counts did not simply scale with strain divergence time. Interpreting this pattern is limited by technical factors such as assembly quality and alignment sensitivity that vary with genome divergence. Moreover, pairwise genome comparisons cannot unambiguously assign individual SVs to specific evolutionary lineages. Additional high-quality outgroup genomes will be needed to place each SV on a particular branch of the medaka phylogeny.

Evolution of Teratorn copies

The giant (180 kb) DNA transposon Teratorn has two subtypes, which share ∼88% sequence identity within their coding regions, with an ∼80 kb inversion in the middle of subtype 2 (Supplemental Fig. S16A). Our new assemblies revealed the complete architecture of Teratorn, showing that the copy number and composition of each subtype varied among strains, ranging from 35–84 copies for subtype 1 and one to 15 copies for subtype 2 (Supplemental Fig. S16B; Supplemental Table S6). Overall sequence structures were conserved among copies except for small indels and long-repeat expansion in some copies (Fig. 5; Supplemental Figs. S16C,D, S17).

Figure 5.

Evolutionary relationships of Teratorn copies. Maximum-likelihood trees based on six concatenated herpesvirus genes (17,865 nt in total) and 17 host genes (13,803 nt). Branch values indicate support from 100 replicates. Structural categories of Teratorn copies are shown next to Teratorn copy names.

Open in new tabDownload PowerPointLink to figure

Phylogenetic analyses demonstrated that Teratorn copies generally clustered by strain and followed the topology of host genes (Fig. 5; Supplemental Fig. S18; Supplemental Data File 2), supporting independent expansion and vertical inheritance of Teratorn within each strain, although herpesviruses are generally transmitted horizontally without being integrated into the host genome. Teratorn could produce virus particles as it contains the whole genome of the alloherpesvirus species (Inoue et al. 2017). Consistent with this, horizontal transfer, although possibly less frequent, was also suggested by the shorter branch length of the phylogenetic trees between Hd-rR and HNI for subtype 2 compared with that of host genes (Fig. 5; Supplemental Fig. S18).

Most Teratorn copies retained coding capacity; namely, the genes essential for viral replication and propagation are intact (Supplemental Fig. S19A). In addition, the nonsynonymous-to-synonymous mutation rate (d_N/d_S) of herpesvirus genes was much lower than one for almost all copies (Supplemental Fig. S19B), suggesting that herpesvirus genes have undergone purifying selection. Indeed, there seemed to be subclusters with higher proportion of full-length Teratorn copies (especially for subtype 1 in Hd-rR) (Fig. 5), implying that full-length copies exhibit higher propagation capacity. In contrast, medaka Tol2, one of the conventional transposons, showed extremely high sequence identity between the copies of each strain and also between the three strains (less than two mutations among 4678 nt aligned regions) (Supplemental Fig. S20) as reported previously (Ichikawa et al. 2017). Tol2 is thus likely to have recently invaded into medaka populations widely distributed in southeast Asia (Koga et al. 2000). Taken together, although the underlying mechanism remains elusive, the transposase gene and the herpesvirus genes may be coordinately utilized for the propagation of Teratorn, presumably through virus-like particles in medaka.

Evolution of the Y Chromosome

In our new assemblies, we clearly defined, for the first time, the Y-differentiating (Y-dif) region as a genomic segment largely homologous to Chr X but showing sequence or structural divergence. The X/Y-dif regions were identified by reduced male read depth in X-dif (Fig. 6A) and by almost no female reads in Y-dif (Supplemental Fig. S21). We focused on Hd-rR and HNI because the HSOK sex chromosomes were less complete, although the sequence accuracy around the Dmy locus in HSOK was supported by raw reads (Supplemental Fig. S21E). In HNI, the Y-dif region spans ∼23.6 Mb and contains 858 genes, extending beyond the centromere. In Hd-rR, the primary Y-dif region spans ∼11.6 Mb and contains 596 genes; it does not extend to the centromere, but instead, two additional small Y-dif regions are present on opposite sides of the centromere (Supplemental Tables S7, S8). The expansion of Y-dif regions along the chromosome is consistent with theoretical predictions (Kirkpatrick and Guerrero 2014). Within the Y-dif region, the portions flanking the Y-sp region exhibit extensive structural divergence, especially on the telomeric side, including large inversions, indels, and long stretches of repetitive sequences (Fig. 6B). Beyond these flanking regions, a certain degree of diversification was broadly detected over the Y-dif region (Fig. 6A), indicating that the structural differentiation extends far beyond the Y-sp region.

Figure 6.

Structural divergence of sex chromosomes and conservation of the Dmy locus. (A) Dot plots between Chr X and Chr Y for Hd-rR and HNI, with male HiFi read depth. Blue and purple dots indicate forward and reverse complement matches, respectively. Yellow lines indicate global average depth. Reduced male read depth marks Y-dif regions. Centromeres showing reduced HiFi depth are supported by much longer ONT-UL reads (Fig. 1C). (B) Structural comparison of X/Y-dif regions (dark gray bars) and Y-sp regions (orange bars) between Hd-rR and HNI. (C) Schematic representation of the structural evolution of the Y-sp region and its corresponding region on Chr 9. Blue and purple boxes indicate genes; arrows, transcriptional direction. A 43 kb segment containing Dmrt1 was duplicated and inserted into the proto-Y Chromosome, giving rise to Dmy.

Open in new tabDownload PowerPointLink to figure

A previous study estimated the Y-sp region of HNI to be ∼258 kb, having expanded from the original 43 kb segment derived from Chr 9 (Kondo et al. 2006). In our assemblies, we precisely determined the Y-sp region in the three strains (250 kb in HNI, 507 kb in Hd-rR, and 351 kb in HSOK) and identified strain-specific patterns of internal repeats (Supplemental Fig. S22). Despite extensive divergence, comparison of the Y-sp regions revealed a small, structurally conserved region surrounding Dmy, spanning 53 kb (from 0.2 kb upstream of exon1 to 3.1 kb downstream of exon6) in HNI, 31 kb in Hd-rR, and 24 kb in HSOK. We refer to these conserved intervals as the candidate “minimum Dmy cassette” (Fig. 6C). Although further experimental validation will be required to confirm its regulatory activity, this conserved cassette could be sufficient for Dmy to function as the sex-determining gene.

Sample preparation and sequencing

Three inbred medaka strains, Hd-rR (southern Japan), HNI (northern Japan), and HSOK (Korea), were used in this study. All three strains were obtained from the Japan NBRP Medaka resource center (https://shigen.nig.ac.jp/medaka/).

For PacBio HiFi sequencing, we used whole-body samples excluding the head, skin, and gut. Male and female samples were prepared separately for each strain. The numbers of fish used were two males and two females for Hd-rR, two males and three females for HNI, and one male and one female for HSOK.

For ONT-UL read sequencing, we did not use whole-body samples because pilot extractions from whole bodies did not yield DNA of sufficient length and quantity for ultralong library preparation. We therefore sought tissues suitable for ultra-high-molecular-weight DNA extraction. As a result, a combination of heart, spleen, and testis yielded the best balance between DNA length and total yield, whereas single tissues alone did not provide sufficient amount of DNA. Hd-rR and HNI tissues from both males and females were mixed for the UL DNA library, whereas HSOK tissues were collected only from males. The numbers of fish used were 15 for Hd-rR, 17 for HNI, and 15 for HSOK. Details of sequencing are provided in the Supplemental Methods.

Genome assembly

For each strain, we assembled HiFi and ONT-UL reads of >50 kb jointly using Hifiasm v0.19.6 (Cheng et al. 2024) and Verkko v2.0 (Antipov et al. 2025). Because the strains are inbred, their genomes were treated as effectively haploid except for sex-specific regions. Therefore, we specified the “-l0” option for Hifiasm and the “‐‐haploid” option for Verkko. Because Hifiasm contigs showed higher contiguity, we primarily proceeded with the Hifiasm assembly, curating it using both Verkko contigs and raw reads. Initial chromosomal scaffolds, named v3.0.0, were obtained by scaffolding the Hifiasm contigs using RagTag v2.1.0 (Alonge et al. 2022) based on the v2.2.4 assembly (NCBI Genomes database [https://www.ncbi.nlm.nih.gov/home/genomes/] accession number GCF_002234675.1).

Sex-differentiating sequences were identified based on female and male HiFi read depths. Specifically, sequences whose male read depths were approximately half of the global average depth were classified as sex-differentiating sequences. Among these, sex-differentiating sequences with normal female read depths were categorized into X-dif, whereas those with almost no female read depths were designated as Y-dif. We searched for sex-differentiating sequences from not only primary contigs but also alternative contigs Hifiasm generated. Details of assembly are provided in the Supplemental Methods.

Semiautomated curation of potentially misassembled regions

We defined potentially misassembled regions (PMRs) in an assembly as intervals flagged by sequence gaps, abnormal ONT-UL read depth, or local accumulation of homozygous variants detected from same-strain reads. ONT-UL reads were used preferentially for these metrics because they provide more unambiguous mappings. Additionally, we assessed assembly quality using k-mer QVs and the annotations of major important elements such as telomeric and centromeric repeats, rDNA arrays, Teratorn, and sex-specific sequences.

We fixed PMRs particularly for completing the Hd-rR genome. First, Verkko contigs were aligned to the scaffolds, and PMRs were replaced with the corresponding Verkko sequence when the number of variants detected becomes smaller. Then we curated most of the remaining PMRs using ONT-UL reads spanning each PMR. For each PMR, we selected the spanning ONT-UL read with the longest flanking sequence as a backbone, aligned other spanning reads to it, calculated a consensus sequence, and substituted the consensus for the PMR. Telomeric sequence extension was performed similarly, except that only one flanking side of each PMR was available. In the end, we curated remaining small homozygous variants detected by either HiFi or ONT-UL reads by taking sequence consensus among mapped reads. We iterated this small variant correction step twice. Details of assembly curation are provided in the Supplemental Methods.

Orienting chromosomes

We categorized each chromosome into one of the metacentric, submetacentric, subtelocentric, acrocentric, and (pseudo-)dicentric. Twelve chromosomes (Chr 1, 2, 7, 8, 9, 11, 12, 13, 18, 19, 22, and 24) were initially oriented from the q-arm to the p-arm and were therefore reverse-complemented. We named these final scaffolds as v3.1 assemblies or simply as v3 assemblies. For compatibility, the same 12 chromosomes in the previous v2.2.4 assemblies were also reverse-complemented, and we named these assemblies as v2.3. In this study, when it is clear whether v2.2.4 or v2.3 is being referred to, it is simply called v2.

Inference of protein-coding regions

Protein-coding genes were predicted from the assemblies by integrating three different programs. LiftOff v1.6.3 (Shumate and Salzberg 2021) was run using the existing protein-coding gene models in GCF_002234675.1_ASM223467v1_genomic.gff associated with the v2 Hd-rR assembly. BRAKER3 v3.0.8 (Gabriel et al. 2024) was run using odb11_actinopterygii_fasta from OrthoDB (Zdobnov et al. 2021) as hints, as well as short-read RNA-seq alignments (NCBI Sequence Read Archive [SRA; https://www.ncbi.nlm.nih.gov/sra] accession number DRP003809) as evidence. Helixer v0.3.3 (Holst et al. 2026) was run using the developer-provided model vertebrate_v0.3_m_0080.h5 with the parameters “‐‐lineage vertebrate ‐‐batch-size 16 ‐‐subsequence-length 213840.”

These three gene sets were merged, and genes encoding peptides shorter than 20 amino acids were removed. Completeness of the annotation was assessed with BUSCO v5.1.2 (Manni et al. 2021) using the actinopterygian BUSCO genes in OrthoDB v10. Details of gene annotation are provided in the Supplemental Methods.

Analysis of centromeric monomers

Centromeric arrays in the v3 assemblies were initially decomposed into monomers with StringDecomposer v1.1.2 (Dvorkina et al. 2020) using the four SF consensus sequences previously defined in the v2 assemblies. For each chromosome, representative monomers were selected by grouping highly similar (>80%) monomers. Representative monomers accounting for >1% of all centromeric monomers on a chromosome were retained.

These chromosome-level representative monomers were compared by hierarchical clustering and principal component analysis. Based on the correspondence with the v2 SFs and the clustering patterns in v3, we redefined medaka centromeric monomers into five SFs. Three v3 SFs corresponded clearly to v2 SF1, SF2, and SF4 and retained these names. The two additional groups represented intermediate monomer families between v2 SF1/SF2 and v2 SF2/SF3 and were designated v3 SF3 and v3 SF5, respectively. The resulting five SF consensus sequences were then used to reannotate centromeric arrays. Monomers with <65% identity to any SF were treated as divergent monomers.

ICUC centromeric monomers were defined as clusters of chromosome-level representative monomers with >97% sequence similarity. Individual monomers were assigned to an ultraconserved monomer class when their similarity to the corresponding representative monomer was >99%.

CpG methylation levels derived from ONT-UL reads were summarized for each centromeric monomer. A raw monomer was classified as hypomethylated or hypermethylated when the average CpG methylation level across its CpG sites was <30% or >70%, respectively. A representative monomer was classified as hypomethylated or hypermethylated when the corresponding class was more frequent among the raw monomers assigned to it. CDRs were detected from 10 kb sliding-window methylation profiles within each array. Partial order alignment graphs were constructed for SF1 and SF3 from each SF consensus sequence and its representative monomers. Details of centromere analysis are provided in the Supplemental Methods.

Teratorn analysis

Teratorn insertions were screened in the three medaka genomes using 1000 nt terminal sequences from both ends of a Teratorn copy as queries. Each Teratorn copy was defined as the interval between the left and right terminal inverted repeats (TIRs). Copies overlapping PMRs were removed in HSOK. Subtype 1 and subtype 2 Teratorn copies were analyzed separately because their terminal sequences share little sequence identity. Structural categorization was performed based on dot plots between each copy and a representative full-length copy. ORFs were identified using ORFFinder (https://github.com/Chokyotager/ORFFinder).

Sequences of six herpesvirus core genes, piggyBac-like transposase, ORF56, Tol2, and 17 host genes (Betancur-R. et al. 2013) were obtained from the three medaka genomes and the genome of Oryzias curvinotus. Multiple alignments were built up for each gene, and alignments were concatenated for copies containing at least five of the six herpesvirus core genes. Maximum-likelihood trees were constructed from these alignments, with branch support estimated from 100 replicates.

Nonsynonymous-to-synonymous ratios (d_N/d_S) for the six concatenated herpesvirus genes were calculated for each medaka inbred strain and each Teratorn subtype. Consensus sequences were generated from multiple alignments, and d_N/d_S values were calculated by comparing each copy with the corresponding consensus sequence using the Nei–Gojobori model. Details of Teratornanalysis are provided in the Supplemental Methods.

Comparative analysis of genomes and sex chromosomes

Genome-wide structural rearrangements were computed using SyRI v1.6.3 (Goel et al. 2019) with the parameter “-x asm5” for genome alignments and were visualized using pyGenomeViz (https://github.com/moshi4/pyGenomeViz; v1.5.0).

Y-dif sequences were identified in each strain based on sex-specific HiFi read depths. Y-sp sequences were identified from the insertion signature in the alignment plot between the X- and Y-dif sequences. Internal repeats were identified with RepeatMasker (http://www.repeatmasker.org; v4.1.5). The Dmy gene was identified from the gene annotation by mapping the previously known Dmy transcript sequence of HNI (NCBI: AB071534.1). Dmy amino acid sequences were compared among the v3 assemblies and the previously registered HNI sequence (NCBI: BAB92012.1). Synteny of regions containing X- and Y-dif sequences was calculated with SyRI and visualized with pyGenomeViz. Details of sex chromosome analysis are provided in the Supplemental Methods.