# Haplotype-resolved genome assembly of a male Dioscorea alata cultivar reveals the structure and evolution of sex chromosomes

Komivi Dossa, Xiaojun Su, Yedomon Ange Bovys Zoclanclounon, Yao Chen, Muhammad Amjad Nawaz, Peng Ni, Ting Qiu, Gang Chen, Erick Malédon, Marie-Claire Gravillon, et al.

► To cite this version:

Komivi Dossa, Xiaojun Su, Yedomon Ange Bovys Zoclanclounon, Yao Chen, Muhammad Amjad Nawaz, et al. Haplotype-resolved genome assembly of a male Dioscorea alata cultivar reveals the structure and evolution of sex chromosomes. 2025. hal-05263805

HAL Id: hal-05263805

https://hal.inrae.fr/hal-05263805v1

Preprint submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Haplotype-resolved genome assembly of a male Dioscorea alata cultivar reveals the structure and evolution of sex chromosomes

**Komivi Dossa**

komivi.dossa@cirad.fr

CIRAD, UMR AGAP Institut

**Xiaojun Su**

Wuhan Benagen Technology Company Limited

**Yedomon Ange Bovys Zoclanclounon**

Rothamsted Research

**Yao Chen**

Wuhan Benagen Technology Company Limited

**Muhammad Amjad Nawaz**

National Research Tomsk State University

**Peng Ni**

Central South University

**Ting Qiu**

Wuhan Benagen Technology Company Limited

**Gang Chen**

Wuhan Benagen Technology Company Limited

**Erick Malédon**

CIRAD, UMR AGAP Institut

**Marie-Claire Gravillon**

CIRAD, UMR AGAP Institut

**Christophe Perrot**

CIRAD, UMR AGAP Institut

**Elie Nudol**

CIRAD, UMR AGAP Institut

**Lévy Laurent**

CIRAD, UMR AGAP Institut

**Moussa Diouf**

Université Cheikh Anta Diop

**Petre I. Dobrev**
**Benjamin Karikari**

University for Development Studies

**Bertrand Fouks**

UMR AGAP Institut, Univ Montpellier, CIRAD, INRAE, Institut Agro

**Gemma Arnau**

UMR AGAP Institut, Univ Montpellier, CIRAD, INRAE, Institut Agro

**Pierre Mournet**

UMR AGAP Institut, Univ Montpellier, CIRAD, INRAE, Institut Agro

**Jianxin Wang**

Central South University

**Sylvain Santoni**

UMR AGAP Institut, Univ Montpellier, CIRAD, INRAE, Institut Agro

**Hana Chair**

UMR AGAP Institut, Univ Montpellier, CIRAD, INRAE, Institut Agro

**Research Article**

**Keywords:** dioecy, root and tuber crops, structural variants, genetic gain, genome

**Posted Date:** September 9th, 2025

**DOI:** https://doi.org/10.21203/rs.3.rs-5627324/v2

**License:** © Ⓟ This work is licensed under a Creative Commons Attribution 4.0 International License.

Read Full License

**Additional Declarations:** The authors declare no competing interests.
1 Haplotype-resolved genome assembly of a male *Dioscorea alata* cultivar
2 reveals the structure and evolution of sex chromosomes

3 Komivi Dossa¹,²*, Xiaojun Su³, Yedomon Ange Bovys Zoclanclounon⁴, Yao Chen³,
4 Muhammad Amjad Nawaz⁵,⁶, Peng Ni⁷,⁸, Ting Qiu³, Gang Chen³, Erick Malédon¹,², Marie-
5 Claire Gravillon¹,², Christophe Perrot¹,², Elie Nudol¹,², Lévy Laurent¹,², Moussa Diouf⁹, Petre
6 I. Dobrev¹⁰, Benjamin Karikari¹¹, Bertrand Fouks²,¹², Gemma Arnau²,¹², Pierre Mournet²,¹²,
7 Jianxin Wang⁷,⁸, Sylvain Santoni², Hanâ Chair²,¹²,*

8

¹CIRAD, UMR AGAP Institut, 97170 Petit Bourg, Guadeloupe, France

10 ²UMR AGAP Institut, Univ Montpellier, CIRAD, INRAE, Institut Agro, F-34398 Montpellier, France

11 ³Wuhan Benagen Technology Company Limited, 430070 Wuhan, Hubei, China

12 ⁴Plant Sciences and the Bioeconomy, Rothamsted Research, Harpenden, AL5 2JQ, UK

13 ⁵Advanced Engineering School (Agrobiotek), National Research Tomsk State University, Lenin Ave, 36,
14 634050, Tomsk, Tomsk Oblast, Russia

15 ⁶Laboratory for Research and Application of Supercritical Fluid Technologies in Agro-Food Biotechnology,

16 National Research Tomsk State University, Lenin Ave, 36, 634050, Tomsk, Tomsk Oblast, Russia

17 ⁷School of Computer Science and Engineering, Central South University, Changsha 410083, China

18 ⁸Hunan Provincial Key Lab on Bioinformatics, Central South University, Changsha 410083, China

19 ⁹Département de Mathématiques et Informatique, Faculté des Sciences et Techniques, Université Cheikh Anta
20 Diop, BP 5005 Dakar-Fann, 10700, Dakar, Senegal

21 ¹⁰Institute of Experimental Botany, Czech Academy of Sciences, Laboratory of hormonal regulations in plants,
22 Rozvojova 263, 165 02 Prague 6, Czech Republic

23 ¹¹Department of Agricultural Biotechnology, University for Development Studies, Nyankpala campus, Tamale,
24 Ghana

25 ¹²CIRAD, UMR AGAP Institut, F-34398 Montpellier, France

26 *Corresponding authors: Komivi Dossa, [komivi.dossa@cirad.fr](mailto:komivi.dossa@cirad.fr), Tel: +590590386161; Hanâ Chair, [hana.chair@cirad.fr](mailto:hana.chair@cirad.fr), Tel: +33(4)32722298.

28

29 Abstract:

30 *Dioscorea alata* L., the most widely cultivated yam species, exhibits dioecy (XY system) and
31 a strong male-biased sex ratio. These two major constraints limit parental combinations and
32 hinder breeding progress. To get insight into the sex chromosome structure and evolution in *D.
33 alata*, we present a haplotype-resolved, near-complete genome assembly of the male cultivar
34 ‘Kabusa’, spanning 958 Mb across 40 chromosomes. DalaChr6A is identified as the Y
35 chromosome which exhibits early signs of heteromorphism, including a slight reduction in size
36 (~400 Kb difference) and a potential centromere shift relative to the X chromosome. Our
37 findings also reveal a sex chromosome turnover between *D. alata* and *D. rotundata*. The sex
38 chromosomes of *D. alata* are evolutionarily young (~4.32 million year ago) and emerged after
39 the divergence from the *D. rotundata* lineage. The sex-determining region (SDR) is refined to
40 ~7.6 Mb, representing ~44% of the Y chromosome. It contains several inversions and a
41 divergence gradient was observed across these inversions. INV4 identified as the oldest
42 pericentric inversion likely marks the early step in *D. alata* SDR evolution. Despite structural
43 divergence, both X and Y chromosomes remain transcriptionally active. Among the 231 genes
annotated in the SDR, 97 are sex-biased and enriched in functions associated with floral organ
formation and hormonal signaling pathways. This study enhances our understanding of sex
chromosome evolution and sex determination in dioecious plants. It provides a gold-standard
reference genome for the *Dioscorea* genus and lays the foundation for accelerated breeding and
genetic improvement in yam.

**Key words:** dioecy, root and tuber crops, structural variants, genetic gain, genome.

50

51

## Main

Dioecy, the condition of separate male and female individuals within a species is present
in only 5–6% of flowering plants species¹. Dioecy poses specific challenges for plant breeding,
particularly when sex expression is stable and strictly unisexual, such as the need for controlled
pollination and limited access to one sex within breeding populations. However, it also offers
unique opportunities, including the development of hybrid breeding strategies, exploitation of
heterosis, and exploration of sex-determining mechanisms². Recent advances have provided
insights into the genetic and molecular mechanisms underlying sex determination across a
variety of plant species. In crops such as papaya³, kiwifruit⁴,⁵, and date palm⁶, a two-gene model
of sex determination has been reported. This model involves two genetically linked genes
located within a non-recombining sex-determining region (SDR), where one gene promotes the
development of one sexual function and the other suppresses the opposite function⁴. In
kiwifruit, the gene *FRIENDLY BOY* (FrBy) has been identified as a crucial male-promoting
factor, while another gene, *SHY GIRL* (SyGI), acts as a suppressor of female organ
development⁴,⁷. The interplay between these two genes is essential for determining the sexual
phenotype of the plant. Their presence within the non-recombining SDR prevents
recombination, thus preserving their functional relationship across generations⁴,⁶. However, in
persimmon⁸ and poplar⁹, a one-gene model appears to control sex differentiation, pointing to
diverse evolutionary pathways. In many species, dioecy is leaky or labile, allowing for a degree
of reproductive flexibility that can sometimes be leveraged in breeding programs².

To date, more than 55 dioecious species genomes have been sequenced, with sex
chromosomes identified in over 35 species¹⁰,¹¹. However, in most cases, only one sex type has
been sequenced, typically the homogametic sex, resulting in incomplete coverage of sex
chromosomes in genomic data¹¹,¹². A near-complete or telomere-to-telomere (T2T) genome
assembly is essential for addressing key biological questions covering genomics-assisted
breeding, the analysis of repetitive elements, genome organization and structure, evolutionary
studies, gene discovery, and sex determination, among other areas of research¹³,¹⁴,¹⁵,¹⁶. For
instance, in maize, the T2T genome assembly has provided insights into the structural variations
(SV) within centromeric regions that were overlooked in previous assemblies¹⁷. This has led to
a better understanding of how these variations influence traits like kernel size and yield,
showcasing the potential for improving crop performance through informed breeding
strategies¹⁷. Moreover, genome assemblies of diploid organisms have collapsed the two
homologous chromosomes into a single mosaic reference genome. While this approach
simplifies the assembly process, it inherently obscures important genetic differences between
the homologous chromosomes, which can be particularly significant in plant species with high
heterozygosity¹⁸ and/or sexual chromosomes¹⁹. This loss of allelic information limits the
understanding of genetic variation and can affect downstream analyses, such as identifying SV,
allele-specific gene expression, sex related genes in plant exhibiting dioecy.

Yams (*Dioscorea* spp.) are economically important tuber crops in tropical regions, with approximately 90% of global production concentrated in West Africa²⁰. *Dioscorea alata*, commonly known as water yam or greater yam, is one of the most widely cultivated yam species due to its high productivity, resilience to environmental stresses, and rich nutritional profile²¹. *Dioscorea alata* is reported to be strictly dioecious and mixed-cytotype, with chromosome numbers of 2n = 40, 60, or 80²². Its haploid genome size is estimated to be ~454 Mb²². Breeding efforts to enhance the yield, disease resistance, and overall tuber quality of *D. alata* have been significantly hampered by its dioecious nature and sex ratio distortion towards male cultivars. This limits the availability of parental combinations for hybridization, hinders the development of improved cultivars and slows genetic gains²³. To overcome these challenges, developing methods to control and manipulate sex phenotype in *D. alata* could facilitate breeding, thereby allowing for more effective selection and improvement of desirable traits. Recent genetic studies have identified Chromosome 6 as the sex chromosome in *D. alata*, supporting an XY (male heterogametic) sex determination system²⁴,²⁵,²⁶. Despite these advances, the development of a reliable sex-specific marker, as well as the elucidation of the molecular mechanisms and identification of key genes governing sex determination in *D. alata*, remain elusive.

The available reference genomes for greater yam is a collapsed assembly of two female
(XX) individuals²⁷,²⁸. The lack of a high-quality Y chromosome sequence restricts
comprehensive sex chromosome analysis in this species. To get insight into the sex
chromosome structure and evolution in greater yam, we employed a combination of long-reads
and chromosome conformation sequencing technologies to assemble a haplotype-resolved,
111 near-complete genome of a male (XY) individual. We also conducted in-depth investigation of
112 the SDR and analyzed the expression profiles of sex-linked genes underpinning sex phenotype
113 expression in *D. alata*. The genome revealed a high degree of similarity between X and Y
114 chromosomes with a 7.6-Mb SDR representing a major pericentric inversion.

115

116 Results

117

118 Genome assembly for the greater yam male cultivar ‘Kabusa’

119 To investigate the genetic mechanisms underlying sex determination in *Dioscorea alata*,
120 we constructed a chromosome-level, haplotype-resolved, near-complete assembly of the male
121 cultivar 'Kabusa', allowing us to generate sequences of both the X and Y sex chromosomes.
122 'Kabusa' is the most popular cultivar in West Indies with a round-shaped tuber, moderate
123 resistance to anthracnose disease and good tuber quality²⁹ (Fig 1a). The haploid genome size of
124 'Kabusa', a diploid cultivar (2n = 40; Supplementary Fig. 1), was estimated to be ~494.13 Mb
125 based on 19-mers analysis, with an unsurprising high heterozygosity level (1.46%)
126 (Supplementary Tables 1,2; Supplementary Fig. 2). We generated 126.5× sequence coverage
127 of raw ultra-long Oxford Nanopore Technology (ONT) data, 55.1× coverage of PacBio HiFi
data and 123.8× coverage of chromatin conformation capture (Hi-C) data for assembling the
129 'Kabusa' genome (Supplementary Tables 3,4,5). Sequencing depths were calculated based on
130 the estimated haploid genome size (1C).

131 Four assembly strategies including Hifiasm (HiFi)³⁰, Hifiasm (HiFi+ONT), Hifiasm
132 (HiFi+Hi-C), and Nextdenovo (ONT)³¹ were evaluated to identify the most effective
133 preliminary backbone assembly to prioritize phasing accuracy. The combination of HiFi and
134 Hi-C reads assembled using Hifiasm produced the most robust foundation for the genome
135 assembly with contig N50 of 12.23 Mb and 17.46 Mb for HaplotypeA and HaplotypeB,
136 respectively (Supplementary Table 6). It is worth noting that, although the Hifiasm
137 (HiFi+ONT) assembly demonstrated good contiguity metrics, the relatively low base-level
138 accuracy of ONT reads may compromise phasing precision. Incorporating Hi-C data in the
139 assembly process theoretically enhances the accuracy of haplotype phasing results and enables
140 chromosome scaffolding. After manual adjustment, a total of 958.71 Mb of contig sequences
141 (contig N50 = 23.83 Mb for HaplotypeA and contig N50 = 24 Mb for HaplotypeB) were
142 mapped to 40 chromosomes in the hybrid assembly version (Supplementary Table 7; Figs 1b,
143 c). The assembly was iteratively polished by HiFi and Illumina data. HiFi reads were aligned
144 to the genome assembly for telomere sequence extension and repair. The telomere sequence
identification results, based on the repetitive units CCCTAAA for the 5' end and TTTAGGG
for the 3' end, are shown in **Supplementary Table 8** and **Supplementary Fig. 4**. Gap filling was
performed using ONT and HiFi contigs as well as raw reads. After gap filling, we assessed the
reliability of the filled regions by aligning both preliminary alternative assemblies and ONT
reads to the gap-filled genome. A gap was considered reliably filled if either data type provided
continuous coverage spanning both flanking regions (**Supplementary Fig. 5**). In total, out of the
54 gaps detected in the assembly, 51 were successfully filled (**Supplementary Table 9**). Only
DalaChr4A, DalaChr4B and DalaChr13B contained single gaps and 37 out of the 40
chromosomes were covered by single contigs (**Supplementary Table 10**). A total of 99.88% of
the contig sequences was anchored on the chromosomes. Altogether, statistics shown in **Table
1** confirmed that we generated a highly accurate and contiguous chromosome-level haplotype-
resolved assembly for 'Kabusa'. The genome assembly is available via Yam Genome Hub
(https://yam-genome-hub.southgreen.fr/).

158

159 **Figure 1. Chromosome-level haplotype-resolved near-complete assembly of the male *D. alata* cultivar 'Kabusa'**. (a) Pictures of 'Kabusa' cultivar displaying the aerial biomass in the field, the tuber shape and flesh color, and the male flowers. (b) Genome-wide Hi-C contact heatmap, visualized with Juicebox³². (c) Circle diagram of basic genome information. Concentric circles from outside to inside represent chromosomes, gene density, repeat density, and GC density calculated in a 50 Kb window. Lines in the interior indicate collinearity. The scale bar = 10 cm. (d) Syntenic mapping between the two haplotypes of the 'Kabusa' genome using SyRI³³. Collinear regions between HaplotypeA and HaplotypeB are indicated with gray lines, while inversion, translocation, and duplication regions are indicated with orange, green, and blue lines, respectively.

169

170 Annotation of the ‘Kabusa’ genome and analysis of sequence variation between
171 haplotypes

172

173 The annotation of the two haplotype genome assemblies revealed remarkably similar
174 features. About 63% of each haplotype genome sequences was annotated as repetitive elements
175 (Supplementary Table 11). LTR retrotransposons were identified as the most abundant
176 sequences (43% and 44%), of which Gypsy (20% and 21%) was the top superfamily in
177 HaplotypeA and HaplotypeB, respectively. Moreover, 8% of each haplotype genome sequence
178 was made of tandem repeat sequences (Supplementary Table 11).

179 ONT-based full-length transcriptome data generated from eight tissues collected at
180 developmental stages were used for genome annotation (Supplementary Table 12). We
181 employed an automated annotation pipeline supported by ONT-based full-length transcriptome
182 data, *ab initio* gene prediction, and protein homology, successfully annotating 27,767 and
27,609 protein-coding genes in HaplotypeA and HaplotypeB, respectively (Table 1; Supplementary Table 13). In total, 21,652 allelic gene pairs were identified between the two haplotypes. Functional annotations were assigned to more than 94% of genes in each haplotype using multiple reference databases.

BUSCO assessment of gene annotation completeness demonstrated 97.5% and 97.8% completeness in the two haplotypes, confirming the high quality of the gene annotations. (Supplementary Table 14). In addition, 139 miRNAs, 376 tRNAs, 2642 rRNAs, and 8609 snRNAs (HaplotypeA) and 140 miRNAs, 374 tRNAs, 3521 rRNAs, and 8738 snRNAs (HaplotypeB) were annotated (Supplementary Table 13). We also performed putative centromere detection by resolving concentration of transposable elements (TE) and genes along all the 40 chromosomes, with size ranging from 122,139 to 5,871,010 bp (Supplementary Table 15; Supplementary Figs. 4,6). For example, the homologous copies of the chromosomes DalaChr6 and DalaChr20 exhibit significant differences in putative centromere position. To enhance the gene annotation of the Dala_V2 female genome assembly²⁷, we projected gene models from the high-quality haplotype-resolved assembly of ‘Kabusa’ (HaplotypeA) onto the Dala_V2 assembly using Liftoff v1.6.3³⁴ with minimap2 v2.24-r1122³⁵. A total of 27,240 genes were successfully mapped to Dala_V2, with only 527 genes failing to align, likely representing Kabusa-specific sequences. Among the mapped genes, 99.9% were placed on chromosomal sequences. The updated annotation (Dala_V3, https://zenodo.org/records/15700131) represents an increase of 2051 genes compared to the previous version (Dala_V2, 25,189 genes), highlighting the improved completeness and accuracy of gene discovery using the ‘Kabusa’ HaplotypeA as a reference (Supplementary Table 16).

We compared the two haplotype assemblies to highlight their sequence variation. A total of 4.2 M single nucleotide polymorphisms (SNPs) including 1.17 M genic SNPs and 741,026 small insertions and deletions (indels), of which half located in genic regions, were identified between the two haplotypes (Supplementary Table 17). In addition, 49,694 structural variants (SVs) were found, including 5075 interchromosomal translocations, 202 inversions, 24,563 duplications, 9824 large insertions and 10,030 large deletions (Fig 1d). Nearly half of these SVs were located in genic regions (Supplementary Fig. 7A). Functional enrichment analysis of SV-influenced genes revealed a wide range of biological pathways, including nucleic acid binding, DNA integration, defense responses, post-embryonic development, regulation of cell shape, and chloroplast organization (Supplementary Fig. 7B)

<caption>**Table 1.** Summary statistics of genome assembly and gene annotation of ‘Kabusa’</caption>
<table><thead><tr><th>Assembly</th><th>'Kabusa' HaplotypeA</th><th>'Kabusa' HaplotypeB</th><th>Combined</th></tr></thead><tbody><tr><td>Number of chromosomes</td><td>20</td><td>20</td><td>40</td></tr><tr><td>Number of contigs</td><td>21</td><td>22</td><td></td></tr><tr><td>Contig N50 (bp)</td><td>23,828,568</td><td>24,039,800</td><td></td></tr><tr><td>Number of gaps</td><td>1</td><td>2</td><td></td></tr><tr><td>Assembly length (bp)</td><td>479,656,039</td><td>480,287,338</td><td>958,714,780</td></tr><tr><td>Illumina read-mapping rate (%)</td><td>99.73</td><td>99.72</td><td></td></tr><tr><td>Unanchored scaffolds</td><td>5</td><td>5</td><td>10</td></tr><tr><td>HiFi read-mapping rate (%)</td><td>100</td><td>100</td><td></td></tr><tr><td>ONT read-mapping rate (%)</td><td>99.65</td><td>99.62</td><td></td></tr><tr><td>Genome assembly BUSCO (%)</td><td>98.1</td><td>97.9</td><td></td></tr><tr><td>Genome assembly QV</td><td>53.1006</td><td>52.3324</td><td></td></tr><tr><td>Annotated genes</td><td>27,767</td><td>27,609</td><td></td></tr><tr><td>Genome annotation BUSCO (%)</td><td>97.5</td><td>97.8</td><td></td></tr><tr><td>Mean mRNA_length (bp)</td><td>5,029.52</td><td>5,077.03</td><td></td></tr><tr><td>Mean cds_length of per gene (bp)</td><td>1,117.66</td><td>1,117.43</td><td></td></tr><tr><td>Mean exon_number of per gene</td><td>5.27</td><td>5.28</td><td></td></tr><tr><td>Mean exon_length (bp)</td><td>302.35</td><td>300.83</td><td></td></tr><tr><td>Mean intron_length (bp)</td><td>801.46</td><td>811.85</td><td></td></tr></tbody></table>

## Detecting 'Kabusa' Y chromosome and sex chromosome evolution in yams

The haplotype-resolved genome of 'Kabusa' offers the opportunity to determine the Y sequence that is absent in the previous female genomes<sup>27,28</sup>. *Dioscorea alata* Chromosome 6 has been identified as the sex chromosome via multiple studies<sup>24,25,26</sup>. Although discrepancies exist in the precise location of the sex-determination region (SDR), evidence points to the presence of large SVs on the sex chromosome. The assembled DalaChr6A was slightly smaller than DalaChr6B (~400 Kb difference) and a large pericentric inversion (from 1,619,245 bp to 9,202,430 bp) was found between them (Fig 1d). A total of 1087 and 1102 genes were annotated in DalaChr6A and DalaChr6B, respectively, with 840 shared genes. Aside from the single pericentric inversion, the two chromosomes exhibit strong collinearity in their distal regions, which likely correspond to the pseudoautosomal regions (PARs).

We compared the sequences of DalaChr6A and DalaChr6B to Chromosome 6 of the female genome Dala_V2 (DalaChr6F)<sup>27</sup> and found that DalaChr6A was the least collinear to DalaChr6F (Supplementary Fig. 8). Remarkably, the region from 1,619,245 bp to 9,202,430 bp on DalaChr6A displayed high variation with DalaChr6F. Moreover, intergenomic comparisons revealed a pericentric inversion between DalaChr6A and DalaChr6F, which was absent between DalaChr6B and DalaChr6F (Supplementary Fig. 9). We confirmed the existence of
this inversion with HiFi long reads spanning the two breakpoints of the inverted regions (Supplementary Fig. 10). We therefore concluded that HaplotypeA genome in ‘Kabusa’ contains the Y chromosome (DalaChr6A) and the pericentric inversion represents the SDR.

We aligned the ‘Kabusa’ assemblies to the sequenced genomes of related *Dioscorea* species, focusing on chromosomes previously reported to be associated with sex determination (Supplementary Fig. 11). *Dioscorea alata* and *D. rotundata*, both from section Enantiophyllum²³, diverged approximately 17.4 million years ago (Mya). We observed strong collinearity between DalaChr6 and DrotChr6. Notably, the sex chromosome of *D. rotundata* has been identified as DrotChr11³⁶, suggesting a sex chromosome turnover between these two species, which also differ in sex determination systems (XY vs. ZW).

In more distantly related species of section Stenophora (2n = 20), such as *D. tokoro* and *D. zingiberensis*, the SDRs have not yet been fully characterized. However, previous work suggests that DtokChr3 harbors sex-determining loci in *D. tokoro*³⁷. In our analysis, part of DtokChr3 aligns with DalaChr6 (Supplementary Fig. 11), and a similar pattern is observed between DzinChr2 and DalaChr6. While this may reflect conservation of ancestral chromosomal segments, we cannot conclude that the SDR itself is shared. Moreover, differences in ploidy could influence genome architecture and synteny. With an increasing number of yam genome sequences becoming available³⁸,³⁹, further studies are needed to identify the SDRs in these species and to determine whether the observed synteny reflects functional conservation or independent evolutionary trajectories.

## Structure of the sex-determination region in *D. alata*

To gain insights into the structure of the SDR, we generated Hi-C contact map for the X and Y chromosomes and further conducted a detailed comparison of the X and Y haplotypes of *D. alata* ‘Kabusa’ using SyRI³³ (Fig. 2a,b). We could observe four different inversions: INV1- INV4 (Table 2). INV1 and INV4, located at the 5' and 3' boundaries respectively, were also detected in comparisons between Y chromosome and two independently assembled X chromosomes from female genomes (Dala_V2²⁷ and Ziyushenshu²⁸). INV2 was also detected in the Y vs. Ziyushenshu comparison, but not in Dala_V2, likely due to the lower resolution of Dala_V2 genome. These observations support a model in which inversions have structured the SDR during the evolution of sex chromosomes in *D. alata*. Between the four inversions, blocks enriched in repeated sequences were observed. These regions are non-collinear with the X chromosome assemblies and may represent putative male-specific region of the Y chromosome (MSY), contributing to the structural differentiation of the Y chromosome. Reciprocally, the
corresponding regions on the X chromosome, which lack homology with the Y, may define X-specific regions (XSR). Interestingly, we observed that the SDR is approximately 400 Kb shorter than the homologous X region (HXR). While such size differences can occur in phased assemblies, the reduction here corresponds to specific deletions and rearrangements detected in the Y chromosome, supporting the hypothesis of a deletion event during sex chromosome differentiation (Table 2).

To further delineate the SDR in the *D. alata* genome, we employed five complementary analytical approaches. First, we constructed female and male read pools using resequencing data from 37 *D. alata* accessions (17 females and 20 males) from Mota et al.⁴⁰. Reads were aligned to the HaplotypeA genome of ‘Kabusa’ to extract variants. This yielded 15.5 million SNPs, from which we calculated the fixation index (*Fst*) between the two sex groups. The highest *Fst* values were observed between 2 and 8 Mb on DalaChr6A (Fig. 2c). Second, analysis of male-specific SNPs, defined as variants homozygous in the female pool and heterozygous in the male pool, identified 4,948 such SNPs. Approximately one-quarter of these variants were located between 2.6 and 9.1 Mb on DalaChr6A (Fig. 2c). Third, we conducted genome-wide association studies (GWAS) using a high-density variant set (>2 M SNPs and 658,000 InDels) from an independent population comprising 68 males and 40 females (Figs. 2d,e). Sex phenotypes were scored over two years at two sites in Guadeloupe (West Indies). GWAS confirmed DalaChr6 as the sole sex chromosome, with 2152 significant SNPs (1,589,777–9,195,813 bp) and 383 significant InDels (1,638,603–9,195,052 bp) associated with sex (*p* < 1e-6; Supplementary Table 18). Notably, most sex-associated SNPs were located within the four inversions identified between Y and X (Figs. 2f,g). Previously reported sex-linked markers²⁵,²⁵,²⁶ also localized within INV1 and INV4. The colocalization of INVs with sex-associated markers suggests that these structural rearrangements have played a central role in the evolution and differentiation of the SDR in *D. alata*. Fourth, we examined differential read coverage (10 Kb window) between male and female genotypes along the Y chromosome (Fig. 2h). Several regions within the SDR displayed markedly higher coverage in males, corresponding to the boundaries of the putative MSY. These regions exhibited low variant density and contained few significant GWAS markers, suggesting suppressed recombination and high sequence divergence. Finally, we analyzed male-specific 27-mers (> 8.2 M) along the Y chromosome, which showed a distribution pattern that precisely mirrored the male-biased coverage (Fig. 2i), further supporting the presence of Y-specific sequences.
The SDR contains 231 annotated genes, while the HXR includes 248 genes (Supplementary Table 19). Upstream, the SDR is delimited by *DalaChr6AG081680*, encoding a 3-isopropylmalate dehydrogenase, and downstream by *DalaChr6AG083980*, encoding a transmembrane protein 50. Only 108 genes had homologous counterparts in the HXR and were considered putative gametologues. An additional five genes within the SDR showed homology with PAR-linked genes in the X chromosome. Conversely, 118 genes within the SDR lacked identifiable X-linked homologues and are considered Y-specific while 140 genes are considered X-specific. We identified 43 and 41 genes on the Y and X chromosomes, respectively, that have functional homologues on autosomal chromosomes in the *D. alata* 'Kabusa' genome, including *DalaChr18*, *DalaChr19*, *DalaChr2*, *DalaChr7*, and *DalaChr16A* (Supplementary Table 20). Notably, only four autosomal genes (from *DalaChr19*) have homologues located within *INV1* of the SDR and HXR, while another four genes (from *DalaChr16A*) have homologues in the XSR. The majority of these syntenic gene pairs are found within the PARs of both chromosomes, suggesting that these relationships likely predate sex chromosome differentiation.

Together, these lines of evidence point to a region of pronounced structural and sequence divergence between the X and Y chromosomes, consistent with an evolving non-recombining SDR (Fig. 2j).

<table><caption>**Table 2:** Structure and gene content of the SDR and HXR in the 'Kabusa' genome</caption><thead><tr><th>Structure</th><th colspan="3">SDR (DalaChr6A)</th><th colspan="3">HXR (DalaChr6B)</th></tr><tr><th></th><th>Location/<br>(Mb)</th><th>size</th><th>Gene content</th><th>Location:<br>(Mb)</th><th>size</th><th>Gene content</th></tr></thead><tbody><tr><td>INV1</td><td>1.6-3.8/2.2</td><td>144</td><td></td><td>7.5-9.6/2.1</td><td>141</td><td></td></tr><tr><td>INV2</td><td>6.15-6.3/0.15</td><td>15</td><td></td><td>4.6-4.75/0.15</td><td>11</td><td></td></tr><tr><td>INV3</td><td>6.65-7.1/0.45</td><td>3</td><td></td><td>6.55-7/0.45</td><td>9</td><td></td></tr></tbody></table>
<table><tr><td>INV4</td><td>8.86-9.2/0.34</td><td>15</td><td>1.6-1.79: 0.19</td><td>9</td></tr><tr><td>MSY/XSR</td><td>-/4.46</td><td>54</td><td>-/5.11</td><td>78</td></tr><tr><td>Total</td><td>1.6-9.2/7.6</td><td>231</td><td>1.6-9.6/8</td><td>248</td></tr></table>

**Figure 2. Structure of the sex determination region (SDR) in *D. alata*.** (a) Hi-C contact heatmap between DalaChr6A (Y) and DalaChr6B (X), visualized with Juicebox [31]. (b) Syntenic mapping between ‘Kabusa’ Y chromosome and X chromosomes from ‘Kabusa’, Dala_V2<sup>27</sup> and Ziyushenshu<sup>28</sup> using SyRI<sup>33</sup>. Four inversions (INV1-INV4) were detected. (c) The genomic landscape of genetic variants between male and female pools. The outer track represents the 20 chromosomes (chr) of HaplotypeA in the ‘Kabusa’ genome. The inner tracks from outside to inside represent *Fst* between females and males and male-specific SNP density, respectively (window size 100 Kb). (d) Manhattan plot for SNP-sex phenotype genome-wide association study (GWAS), with the peaks indicating significant GWAS signals, and the black and red dotted horizontal lines indicating the suggestive and genome-wide significance thresholds, respectively. (e) Manhattan plot for InDel-sex phenotype GWAS, with the peaks indicating significant GWAS signals, and the black and red dotted horizontal lines indicating the suggestive and genome-wide significance thresholds, respectively. Regional Manhattan plots on DalaChr6A showing the dispersion of the significant SNPs (f) and Indels (g). (h) Differential read coverage (10 Kb window) between male and female groups along the Y chromosome. (i) Exact matches of male specific 27-mers along the Y chromosome. (j) Representation of the X (DalaChr6B) and Y (DalaChr6A) chromosomes along with the gene density (heatmap), centromere and telomere locations, pseudoautosomal regions (PAR), and the SDR containing the 4 inversions (INV1-INV4) detected.

## Evolutionary divergence of the sex chromosomes in *D. alata*

We computed synonymous site divergence (*Ks*) between X- and Y-linked allelic gene pairs. In total, we analyzed 87 gametologues within the SDR and 272 one-to-one homologues genes in the PARs. After removing outliers, *Ks* values varied across the SDR-HXR segment of the sex chromosomes (range: 0.00486-0.137; mean = 0.0345; s.d. = 0.0259) and were significantly higher than those in the PARs (range: 0.00205-0.127; mean = 0.0292; s.d. = 0.0245; Kruskal-Wallis *p* = 0.0019) (Fig. 3a), consistent with ongoing divergence between the SDR and the HXR. Notably, the gene *DalaChr6AG083010*, located within INV1, showed the highest *Ks* value and encodes a PLATZ transcription factor previously reported as essential for female flower morphology determination in grapevine<sup>41</sup>.
354

355

356

357

358

359

360

361

362

To gain further insight into the evolution of the SDR, we examined *Ks* distributions within the three chromosomal inversions (INV1, INV2, and INV4) that distinguish the Kabusa Y chromosome from the three assembled X chromosomes (Fig. 3b). We analyzed 74, 7, and 6 gametologues within INV1, INV2, and INV4, respectively. As expected, INV4 (range: 0.0202–0.0724; mean = 0.0438; s.d. = 0.0234) and INV1 (range: 0.00535–0.137; mean = 0.0367; s.d. = 0.0275) exhibited significantly higher *Ks* values than INV2 (range: 0.00486–0.0274; mean = 0.0145; s.d. = 0.00743; Kruskal–Wallis *p* = 0.0072). These results suggest that INV1 and INV4 represent ancestral chromosomal rearrangements associated with the origin of the SDR.

363

364

365

366

367

368

369

370

371

372

373

Using a synonymous substitution rate of r = 4 × 10⁻⁹ substitutions per site per year previously estimated for *Dioscorea*³⁹, we estimated that the divergence between the SDR and HXR in *D. alata* occurred approximately 4.32 Mya. The estimated divergence times of the three inversions suggest that INV4 (~5.48 Mya), followed by INV1 (~4.58 Mya), predate the SDR. These two inversions may therefore have played a role in the initial expansion of the SDR and the evolution of the sex chromosomes, similarly as reported in *Populus euphratica*⁴². In contrast, INV2 exhibits a considerably more recent divergence (~1.81 Mya). Altogether, these findings suggest that the sex chromosomes of *D. alata* are evolutionarily young, with their emergence occurring well after the split from the *D. rotundata* lineage (~17.4 Mya). The data reveal a gradient of divergence across the inverted regions, with INV4 emerging as the oldest pericentric inversion likely involved in the early steps of SDR evolution.

374

375

376

377

378

379

380

381

382

383

384

385

386

387

TEs can influence the dynamics of the SDR by promoting recombination suppression under selective pressure, thereby shaping the organization and structural evolution of sex chromosomes<sup>43,44</sup>. We further analyzed TE abundance and distribution pattern along chromosomes in HaplotypeA and HaplotypeB with an emphasis on LTR/Gypsy and LTR/Copia as the most abundant TEs (Supplementary Fig. 12). Detailed analysis of LTR/Gypsy and LTR/Copia densities across X and Y in ‘Kabusa’ revealed region-specific accumulation patterns. Significant differences in both LTR/Gypsy and LTR/Copia densities were observed when comparing the SDR and HXR to their respective flanking PARs (Figs. 3c, d), indicating a localized enrichment of LTR/Gypsy in sex-linked regions. The comparison between the SDR and HXR did not yield a statistically significant difference in TE density (Fig. 3e), suggesting that the overall levels of TE accumulation in the X and Y are broadly similar. However, smoothed density plots along X and Y highlighted marked peaks of LTR/Gypsy accumulation in regions that correspond to non-colinear segments between haplotypes (Fig. 3f). These regions coincide with the putative MSY, as well as XSR. This finding suggests that both the X and Y
chromosomes harbor structurally TE-rich segments as reservoir of sex-specific divergence. Such patterns align with observations in other dioecious plant species such as papaya⁴⁵ or *Silene latifolia*⁴⁶ where sex chromosome differentiation is associated with the accumulation/specific localization of LTR elements and chromosomal rearrangements. Although global TE density may not differ significantly between the SDR and HXR in ‘Kabusa’, the presence of localized LTR/Gypsy-rich regions likely contributes to the disruption of collinearity and may play a structural role in X–Y divergence. In conclusion, our findings support a model where localized accumulation of LTR/Gypsy elements contributes to sex chromosome differentiation in *D. alata*.

**Figure 3: Sex chromosome divergence in *D. alata*.** (a) Comparison of the synonymous substitution rates (*Ks*) of genes between SDR and PAR on the Y chromosome. (b) *Ks* comparison between Y chromosome components. ** = significant difference at *p* <0.001. (c) Comparison of the LTR/Copia and LTR/Gypsy densities between PAR and SDR in the Y chromosome. (d) Comparison of the LTR/Copia and LTR/Gypsy densities between PAR and HXR in the X chromosome. (e) Comparison of the LTR/Copia and LTR/Gypsy densities between SDR and HXR. (f) Smoothed density plots along X and Y chromosomes. Densities are calculated per 1 Mb overlapping windows.

### Sex-biased gene expression along sex chromosomes and SDR in *D. alata*

In *D. alata*, ontogenetic differences between female and male flowers become morphologically evident between Day 8 and Day 12 after blooming (Fig. 4a), indicating that key sex-determination processes likely occur during this window. To capture relevant transcriptional activity during the differentiation of sexual organs, we collected flower buds at Day 5, when a clear developmental divergence is not yet observed. Transcriptome sequencing was conducted on six female and six male individuals (Supplementary Table 21). Libraries of poor quality were excluded from further analyses (Supplementary Table 22). More than 918 million 150-bp paired-end clean reads were mapped to the haploid consensus genome (HapA + HXR) of Kabusa. While this approach enables consistent gene expression analysis, we acknowledge that some degree of cross-mapping between homologous X and Y regions may occur, potentially introducing minor biases in quantification. These limitations are inherent to mapping in sex chromosome systems with closely related haplotypes. Transcript quantification and differential expression analysis were performed to detect the sex-biased genes (SBG). Male and female groups were clearly separated at the transcriptional level with globally higher gene expression levels within the male group (Fig 4b; Supplementary Fig. 13).
Globally, male-biased genes (4905) were twice as abundant as female-biased genes (2012) (Fig. 4c), a trend frequently reported in dioecious plant species⁴⁷. Gene ontology enrichment analysis of the SBGs revealed several biological pathways including photosynthesis, cellular morphogenesis, pollen wall assembly, oxidoreductase activity, and carbohydrate metabolism. While not all enriched categories are directly involved in sex determination, they may reflect sex-specific physiological processes or developmental demands during floral organ formation in *D. alata* (Supplementary Fig. 14). The most strongly-induced SBGs were located within or in close proximity to the SDR, further supporting the central role of this region in sex differentiation (Fig. 4c). A total of 751 genes on the Y chromosome (69% of total annotated Y genes) and 762 on the X chromosome (69% of total annotated X genes) were transcriptionally active in the sampled tissue types, indicating that both sex chromosomes retain substantial expression capacity and that Y chromosome of *D. alata* has not yet undergone widespread loss of gene expression typically associated with advanced stages of sex chromosome differentiation. A detailed examination of gene expression patterns along the sex chromosomes revealed a clear concentration of male-biased expression along the Y chromosome (Fig. 4d). In total, 230 SBGs were identified on the Y chromosome, compared to 154 SBGs on the X chromosome. As illustrated in Fig. 4e, SBG located in the Y chromosome exhibited significantly higher log₂ fold changes in expression between male and female groups, compared to their counterparts on the X chromosome (*p* = 4.994e-08). Along the sex chromosomes, the magnitude of differential expression was more pronounced within the SDR and HXR, highlighting a potential regulatory role of these regions in sex-specific transcriptional activity (Fig. 4d). Within the SDR, 97 SBGs were detected (Supplementary Table 23), with all but one—*DalaChr6AG081990*, a GERMIN-LIKE PROTEIN SUBFAMILY 3, showing higher expression in male flowers. This supports the idea that the SDR is a transcriptional hub for male-specific functions. Notably, genes located in the MSY displayed markedly stronger expression in males, as would be expected for loci involved in male reproductive development and function (Fig. 4f). In contrast, the HXR contained 45 SBGs, but only a single gene *DalaChr6BG081880*, annotated as *Cytochrome P450* was upregulated in male samples (Supplementary Table 23). Interestingly, these X- and Y-linked SBGs identified in floral tissues also exhibited differential expression in mature tubers (three males *D. alata* accessions compared to three female *D. alata* accessions), suggesting that sex-associated transcriptional divergence may extend beyond reproductive tissues in *D. alata*⁴²,⁴⁸ (Fig. 4g).

Notably, we identified 10 genes within the SDR that were only expressed in males (Supplementary Table 23). These genes may contribute to sex differentiation by promoting
male organ development and/or repressing female organ development. For instance, *DalaChr6AG082390*, annotated as a style cell-cycle inhibitor 1-A (SCI1) homolog, showed strong expression in all male flower samples but was undetectable in female flowers. In *Nicotiana tabacum* and *Arabidopsis thaliana*, SCI1 has been shown to control cell proliferation in pistil development⁴⁹,⁵⁰. While its function in male reproductive tissues has not been reported in these model species, the male-specific expression pattern in *D. alata* raises the possibility of a lineage-specific functional divergence or a novel regulatory role in male floral development. However, functional validation is required to confirm its role in sex differentiation in yam. Several other SBGs within the SDR have known functions in flower sexual organ development, including *DalaChr6AG081680* (3-isopropylmalate dehydrogenase) essential for pollen development and proper embryo sac development⁵¹; *DalaChr6AG082570* (SWR1) involved in female germline development⁵²; *DalaChr6AG083470* (NYFC) known for its role in flowering and organ morphogenesis⁵³; and *DalaChr6AG082810* (ARR) involved in cytokinin (CK) hormone signaling⁵⁴. Since INV4 is the most ancestral inversion associated with sex chromosome evolution in *D. alata*, we hypothesize that SBGs located within this region may represent strong candidates for primary sex-determining factors. Among them, four have homologs in model species with functions related to reproductive development. *DalaChr6AG083890* encodes a rhomboid-like protein, homologous to *Arabidopsis thaliana* genes involved in pollen development and pollen wall formation⁵⁵. *DalaChr6AG083930* encodes a PHD-finger protein with strong similarity to genes essential for male fertility in *A. thaliana*, soybean, and rice⁵⁶⁻⁵⁸. *DalaChr6AG083880* is a homolog of CORONATINE INSENSITIVE 1 (COI1), an F-box protein known to play a central role in jasmonic acid signaling and reproductive organ development across multiple plant species⁵⁹⁻⁶⁰. Finally, *DalaChr6AG083840* is similar to the MATERNAL EFFECT EMBRYO ARREST 9 gene in *A. thaliana*, which is involved in female gametophyte development and function⁶¹. While functional validation remains to be conducted in *D. alata*, the presence of these homologs highlights candidate genes potentially involved in sex determination and reproductive processes in this species.

To explore whether epigenetic regulation may contribute to the observed biased gene expression, we analyzed DNA methylation levels of 5-methylcytosine (5mC) at CpG sites. Globally, quite similar patterns of methylation levels and distribution across homologous chromosomes of HaplotypeA and HaplotypeB could be observed (Supplementary Fig. 15, Fig. 4h). However, a finer-scale comparison of X and Y chromosomes revealed a significantly higher methylation level within the SDR and HXR compared to the PARs (Figs. 4i,j). No
significant difference was observed between the SDR and HXR, suggesting that the epigenetic
landscape of this region is largely conserved despite major structural rearrangements (Fig. 4k).
It is important to note that our analysis focused solely on CpG methylation and did not account
for other types of DNA methylation (e.g., CHG or CHH contexts), which may also play a role
in gene regulation⁶².

496

**Figure 4: Sex-biased gene expression in *D. alata*.** (a) Development of *D. alata* flower at Day5, Day8 and Day12 after blooming in male ‘Kabusa’ and female ‘74F’. Ontogenetic differences between female and male flowers become morphologically evident between Day8 and Day12 after blooming. Scale = 1 cm. (b) Principal component analysis based on the gene expression levels for flower samples (Day5) expressed as fragments per kilobase of transcript per million fragments mapped in the male and female groups. The first and the second principal components explaining 55% and 20% of the total variation, respectively. (c) Volcano plot displaying the differential expression analysis in flowers between male and female groups. Dots in green, red, and gray represent genes downregulated, upregulated and not differentially expressed in male flowers compared to female flowers, respectively. (d) Gene expression patterns in flower buds (log2fold change between male and female groups) along the sex chromosomes. Each dot represents a sex-biased gene (SBG). Statistical comparison of SBG located within the sex chromosomes (e) and within the Y chromosome components (f). (g) Statistical comparison of SGB in mature tuber samples located within the sex chromosomes. Samples represent mature tuber from three female *D. alata* accessions (‘Dou’, ‘Kinabayo’, and ‘74F’) and three male accessions (‘Divin’, ‘Plimbite’, and ‘Toufi-Tetea’). \* = significant difference at *p* <0.05. (h) Comparison of the methylation levels (5-methylcytosine at CpG sites) across sex chromosomes and statistical comparison between sex chromosomes components: SDR vs PAR (i), HXR vs PAR (j) and SDR vs HXR (k). Methylation levels are calculated per 100 Kb overlapping windows.

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

Development of a reliable sex-prediction molecular marker

519

520

521

522

523

524

525

526

527

528

529

530

A sex-prediction PCR allele competitive extension (PACE) marker for *D. alata* was previously developed²⁶. However, it does not achieve 100% accuracy in distinguishing male and female individuals. This suggests that the marker is not tightly linked to the sex determination loci on DalaChr6. In this study, we developed a new PACE marker, *sda6*, based on allelic variation within the coding sequence of one candidate sex-controlling gene *DalaChr6AG083080* (psbP domain-containing protein 6) and its homolog *DalaChr6BG082060*. The *sda6* marker is located within the INV1. From the transcriptome data, all female genotypes were homozygous T/T (DalaChr6B: 7,643,820 bp) while all male genotypes were heterozygous G/T (DalaChr6A: 3,714,267 bp) (Supplementary Table 24). Genotyping 146 *D. alata* accessions using the *sda6* marker revealed two distinct clusters corresponding to male and female groups (Supplementary Fig. 16). Among the 74 accessions with known sex phenotypes, *sda6* accurately distinguished male and female plants with 100%
accuracy (Supplementary Table 25). This result underscores *sda6* as a reliable molecular
marker for predicting the sex phenotype in *D. alata*, addressing an essential requirement for
advancing breeding efforts²⁴.

534

535 Discussion

536
537 The haplotype-resolved genome assembly of the 'Kabusa' cultivar surpasses the quality
538 of existing *D. alata* assemblies<sup>27,28</sup>, establishing a new reference standard for the *Dioscorea*
539
genus. It features a significantly higher contig N50, improved genome annotation, and greater
540 accuracy. The assembly is nearly gap-free and haplotype-resolved, containing a large repertoire
541
of novel genes. Importantly, it includes the previously missing full Y chromosome sequence.
542
As demonstrated in other plant species, such as pear<sup>63</sup> and walnut<sup>64</sup>, haplotype-resolved genome
543
assemblies offer significant advantages over mosaic genomes. These include enhanced
544
resolution for studying complex traits and evolutionary dynamics, especially in species with
545
high heterozygosity or sex-specific chromosomes, such as *D. alata*. The haplotype assemblies
546
of the 'Kabusa' genome provide valuable resources for a wide range of applications, including
547
variant calling, gene expression analysis, genome-wide association studies, genomic prediction,
548
genome editing, and comparative genomics. Additionally, they enable precise evaluation of the
549
specific contributions of each haplotype to phenotypic traits, supporting both functional and
550
applied research efforts<sup>14,65</sup>.

The identification of telomeric and centromeric sequences is crucial for ensuring a
550
complete, high-quality genome assembly. Such completeness enhances the precision of trait
551
mapping (e.g., GWAS, QTL), enabling more accurate association between genetic loci and
552
agronomic traits, which in turn facilitates targeted breeding strategies¹³. In the 'Kabusa'
553
genome, we identified the telomeres and putative centromeric positions for all 40 chromosomes,
554
revealing notable putative shifts in centromere positions between some homologous
555
chromosomes. These shifts can profoundly affect chromosome behavior during meiosis and
556
mitosis⁶⁶. The intricate relationship between centromere dynamics, genomic rearrangements,
557
and their influence on trait variation has been well-documented in various species⁶⁷,⁶⁸,⁶⁹. In the
558
'Kabusa' genome, we also detected several para- and pericentric SVs. Consistent with our
559
findings, a paracentric inversion on linkage group 16 (corresponding to DalaChr16) in 'Kabusa'
560
has been previously identified²⁹. Their study further demonstrated that this SV reduces fertility
561
and recombination rates in 'Kabusa'. The prevalence and impacts of such SVs on *D. alata*
562
breeding programs warrant further investigation. Overall, we anticipate that the high-quality
563
564

‘Kabusa’ genome will play a pivotal role in accelerating greater yam improvement by providing a robust foundation for functional genomics and targeted breeding strategies.

565

566

The sex chromosomes of *D. alata* are young and at the early stage of differentiation, features observed in other plant species such as willows, poplars, strawberry and silene (*Silene otites* and *Silene pseudotites*)<sup>70,71,72,73</sup>. In addition to the pericentric inversion, the Y chromosome exhibits a slight size reduction and a shift in putative centromere position compared to the X chromosome. Cormier et al.<sup>25</sup> demonstrated that *D. alata* sex chromosomes still undergo recombination, albeit at a significantly reduced rate. This reduced recombination can impact the detection of quantitative trait loci (QTL) on these chromosomes. Genes located on sex chromosomes, particularly those linked to important agronomic traits, present both challenges and opportunities for crop breeding<sup>74</sup>. To date, only a single minor QTL associated with tuber oxidative ratio has been identified on DalaChr6<sup>29,75</sup>. It is likely that the reduced recombination rate on DalaChr6 obscures the detection of additional QTL present on this chromosome.

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

The collinearity of sex chromosomes in related species offers valuable insights into their evolutionary dynamics<sup>76</sup>. Our analysis suggests that sex chromosome formation in *D. alata* occurred after its divergence from the *D. rotundata* lineage. In *D. rotundata* (ZW system), a relatively small sex-determining region (SDR) of ~161 Kb containing only three genes has been identified<sup>36</sup>. In contrast, the SDR in *D. alata* spans 7.6 Mb and harbors 231 genes, reflecting a markedly different genomic architecture associated with sex determination. These differences in sex determination systems, as well as in SDR size and structure, indicate a clear case of sex chromosome turnover between the two species. Notably, *D. rotundata* frequently displays sex lability and includes monoecious cultivars<sup>77,78</sup>, whereas such reproductive flexibility has not been reported in *D. alata* to date. Given that yams display diverse sex determination systems, including XY, ZW, and XO (*D. reticulata* and *D. sinuata*)<sup>23</sup>, the *Dioscorea* genus represents an excellent model for studying sex chromosome evolution and turnovers. However, high-quality genome resources of heterogametic sex, like the ‘Kabusa’ genome, are essential to enable accurate detection of sex chromosomes, SDR, and sex-determining genes.

592

593

594

595

596

597

We specifically investigated the SDR of *D. alata*, which corresponds to a large pericentric inversion on the Y chromosome and is composed of multiple inversion-derived evolutionary strata. Our results perfectly match reports in a wide range of animals and plants, where suppression of recombination between sex chromosomes is thought to occur progressively, often mediated by chromosomal inversions. This leads to the formation of evolutionary strata, marked by distinct levels of sequence divergence between the sex
chromosomes<sup>79,80,81,82,83</sup>. The *D. alata* SDR harbors sex-limited, sex-specific and sex-biased genes critical for sex phenotype expression. Many of these SDR-associated genes are predicted to have functions related to flower organ development. This aligns with previous findings that male/female sterility and sex determination are closely tied to processes such as growth and development, pollen viability and mobility, ovule and embryo sac formation, and other gametogenic processes<sup>84</sup>, as well as hormone signaling pathways<sup>85</sup>. The convergent model of dioecy in plants involves a cascade of regulatory pathways, with sex-regulatory genes playing central roles in determining sex<sup>10</sup>. Overall, the candidate sex-determining genes identified in *D. alata* provide a valuable foundation for advancing functional genomics and offer promising targets for future efforts in sex manipulation and breeding in this species.

608

## 609 Methods

610 Plant material

*Dioscorea alata* cultivars ‘Kabusa’ and ‘74F’ were obtained from our germplasm collection of the Centre de coopération Internationale en Recherche Agronomique pour le Développement (CIRAD), Guadeloupe. Tubers of ‘Kabusa’ and ‘74F’ were planted at the CIRAD research station located at Roujol (16°10′ 56″ N, 61° 35′ 24″ W, 10 meters above sea level), Guadeloupe in May 2022. We also planted ten offspring (‘A36’, ‘A94’, ‘A113’, ‘A74’, ‘A96’, ‘A117’, ‘A88’, ‘A111’, ‘A64’, ‘A89’) from the cross ‘Kabusa’ x ‘74F’<sup>25</sup>. Genotyping data from 37 and 108 *D. alata* accessions were retrieved from the population studied by Mota et al.<sup>40</sup>.

619 Library construction and sequencing

During vegetative stage, young leaves of ‘Kabusa’ were collected for high-quality genomic DNA extraction using the cetyl trimethyl ammonium bromide method<sup>86</sup>. DNA quality and concentration were tested by 0.75% agarose gel electrophoresis, NanoDrop One spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and Qubit 3.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). After the DNA quality and integrity were tested, it was randomly sheared by Covaris ultrasonic disruptor. DNBSEQ sequencing pair-end libraries were prepared using VAHTS® Universal Plus DNA Library Prep Kit for MGI V2 (Vazyme, Nanjing, China). For Oxford Nanopore sequencing, the libraries were prepared using the SQK-LSK110 ligation kit according to the standard protocol. The purified library was loaded onto primed R9.4.1 Spot-On Flow Cells and sequenced using a PromethION sequencer

620

621

622

623

624

625

626

627

628

629
630 (Oxford Nanopore Technologies (ONT), Oxford, UK) with 72-hour runs. For PacBio 631 sequencing, the libraries were prepared using the SMRTbell Prep Kit 3.0 (Pacific Biosciences, 632 Menlo Park, CA, USA) according to the standard protocol. The purified library was loaded and 633 sequenced using REVIO Polymerase Kit (Pacific Biosciences, Menlo Park, CA, USA) and 634 REVIO sequencing plate. For Hi-C sequencing, nuclear DNA collected from young leaves was 635 first cross-linked to capture DNA fragment interaction and then digested with a restriction 636 enzyme, Dpn II. Next, sequencing libraries were constructed based on the biotinylated and 637 ligated DNA fragments and subsequently sequenced on Illumina NovaSeq 6000 (San Diego, 638 CA, USA). Mixed tissues (flower, young leaf, old leaf, vine, tuber flesh, tuber skin, root and 639 petiole) of 'Kabusa' were collected at different developmental stages and used for ONT full- 640 length Iso-seq sequencing. The preparation of long-read cDNA libraries was initiated using the 641 ONT PCR-amplified cDNA kit (SQK-PCS109). Specifically, full-length cDNA libraries were 642 constructed from poly(A)+ mRNA using the cDNA-PCR Sequencing kit (SQK-PCS109). 643 Subsequently, the cDNA was PCR-amplified for 13-14 cycles with barcoded adapters from the 644 Oxford Nanopore PCR Barcoding kit (SQKPBK004). Finally, the 1D sequencing adapter was 645 ligated to the cDNA prior to loading onto an R9.4.1 FLO PRO002 flow cell for sequencing on 646 a PromethION platform. Flower buds were chosen at Day5 because they are RNA-rich and at 647 a stage when a clear developmental divergence between male and female is not yet observed 648 (Fig. 4a). Samples from ten offspring and their parents ('Kabusa' x '74F') were collected in 649 triplicate and used for RNA-seq. Total RNA was isolated from young flower buds using the 650 E.Z.N.A.® Plant RNA Kit (Omega Biotek, Norcross, GA, USA) according to the 651 manufacturer's instructions. Quality, purity and quantity of the isolated RNA was checked on 652 Agarose electrophoresis, Nanodrop (NanoDrop Technologies, Wilmington, DE, USA) and 653 Qubit 3.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), respectively. The libraries 654 were prepared using the NEBNext® UltraTM Kit of Illumina® according to the manufacturer's 655 instructions. RNA sequencing was conducted on Illumina Novaseq 6000 (San Diego, CA, 656 USA).

657 Genome assembly

658 The k-mer depth distribution was calculated with Jellyfish (v2.2.10)⁸⁷ based on the 659 DNBSEQ short reads, and GCE (v1.0.0)⁸⁸ was used to estimate the genome size and 660 heterozygosity. Four preliminary genome assemblies were generated:

661 1- with nextDenovo (v2.5.0, parameters: read_cutoff = 1k, blocksize = 1g, 662 nextgraph_options = -a 1)³¹, we performed the preliminary assembly of ONT
663 sequencing data

664 2- HiFi data was assembled with hifiasm, v0.16.1-r375³⁰

665 3- A mixed assembly of HiFi and ONT ultra-long data using hifiasm, v0.18.2-r467

666 4- A mixed assembly of HiFi and Hi-C data using hifiasm, v0.19.5-r587

667 The quality of the four assemblies were evaluated and the 'HiFi + Hi-C' preliminary assembly version was selected as the foundation for the subsequent near-complete genome assembly.

669 Contig sequences from the assembly were clustered into 40 chromosomes by the agglomerative hierarchical clustering method and then oriented using the ALLHiC (v0.9.8)⁸⁹. We convert the pairwise interactions between contigs into specified binary files (.hic files) using 3D-DNA (v180419)⁹⁰ and juicer (v1.6)³². We further refined the ordered and oriented contigs (.assembly) using Juicebox (v1.11.08)³² to manually adjust into a reviewed assembly (.review.assembly). Gap regions were represented with 100 Ns to obtain chromosome-level genome sequences. The Hi-C heatmaps were generated based on contig interaction strength and genomic positions using HiCExplorer (v3.6)⁹¹.

670 Contig sequences from the assembly were clustered into 40 chromosomes by the agglomerative hierarchical clustering method and then oriented using the ALLHiC (v0.9.8)⁸⁹. We convert the pairwise interactions between contigs into specified binary files (.hic files) using 3D-DNA (v180419)⁹⁰ and juicer (v1.6)³². We further refined the ordered and oriented contigs (.assembly) using Juicebox (v1.11.08)³² to manually adjust into a reviewed assembly (.review.assembly). Gap regions were represented with 100 Ns to obtain chromosome-level genome sequences. The Hi-C heatmaps were generated based on contig interaction strength and genomic positions using HiCExplorer (v3.6)⁹¹.

675 Gap regions were represented with 100 Ns to obtain chromosome-level genome sequences. The Hi-C heatmaps were generated based on contig interaction strength and genomic positions using HiCExplorer (v3.6)⁹¹.
676 Gap regions were represented with 100 Ns to obtain chromosome-level genome sequences. The Hi-C heatmaps were generated based on contig interaction strength and genomic positions using HiCExplorer (v3.6)⁹¹.
677 HiCExplorer (v3.6)⁹¹.
678 To proceed to the extension and repair of telomere sequences, we employed Winnowmap (v2.03, parameters: k=15, -MD)⁹². All ONT reads were aligned to the genome and reads aligning once within 50 bp from the chromosome ends were collected. Occurrence of telomere repeat motifs (CCCTAAA/TTTAGGG) was calculated in each read, and we defined the read with the highest occurrence as the reference and others as query. With medaka_consensus (v1.2.1, parameters: -m r941_min_high_g360), we reassembled the reference telomere reads and query telomere reads to obtain the consensus sequence. Thereafter, the consensus telomere sequences were aligned to each chromosome with nucmer (v3.1)⁹³, and we replaced the telomere sequence at the chromosome ends with the best alignment result.
679 Winnowmap (v2.03, parameters: k=15, -MD)⁹². All ONT reads were aligned to the genome and reads aligning once within 50 bp from the chromosome ends were collected. Occurrence of telomere repeat motifs (CCCTAAA/TTTAGGG) was calculated in each read, and we defined the read with the highest occurrence as the reference and others as query. With medaka_consensus (v1.2.1, parameters: -m r941_min_high_g360), we reassembled the reference telomere reads and query telomere reads to obtain the consensus sequence. Thereafter, the consensus telomere sequences were aligned to each chromosome with nucmer (v3.1)⁹³, and we replaced the telomere sequence at the chromosome ends with the best alignment result.
680
681
682
683
684
685
686
687
The next step consists of gap-filling. Winnowmap (v2.03, parameters: k=15, -MD)⁹² was used to align gap-filling data (excluding Ns) to the genome gap regions (including Ns). The priority for gap-filling data was as follows: preliminary assembly genome versions > ONT data (consensus.fa) > HiFi data (ccs.fasta). When the alignment spans exactly across both ends of the gap, the best alignment region with the longest length were chosen to replace the gap-containing sequence on the genome using gap-filling data. After the gap-filling is completed, the replaced regions in the gap-filled genome were aligned with other assembly versions/ONT reads/HiFi reads to verify the reliability of the gap-filling results.
688
689
690
691
692
693
694
695
Genome evaluation
Genome assembly continuity was assessed by locating the gaps on the genome. For genome assembly consistency, the clean short reads were aligned to the genome using bwa (v0.7.17-r1188) while ONT and HiFi reads were mapped to the genome using minimap2 (v 2.17-r941)³⁵. Alignment and coverage rates were calculated. Genome completeness was assessed using BUSCO protein mode with embryophyta_odb10 data set (v 5.3.0; parameters: -evalue 1e-05 --augustus)⁹⁴. Likewise, genome quality and accuracy were evaluated using Meryl (v1.3)⁹⁵. We employed MUMmer (v4.0.0)⁹⁶ to align the ‘Kabusa’ haplotype genome assemblies to the female genome Dala_V2²⁷ to determine the chromosome numbering and orientation. Then, we used SyRI (v1.6)³³ to identify collinear regions and variations between genomes.

706

## 707 Repeat sequence annotation

The same pipeline was applied for HaplotypeA and HaplotypeB. First, the repeat sequences were identified and annotated. Model sequences were predicted based on the genome sequence using RepeatModeler (v2.0.4, parameters: -database mydb -threads 16, https://github.com/Dfam-consortium/RepeatModeler). Additionally, LTR sequences were predicted using LTR_FINDER (version: Official release of LTR_FINDER_parallel, parameters: -threads 16 -harvest_out -size 1000000 -time 300)⁹⁷ and LTRharvest (v1.62, parameters: -minlenltr 100 -maxlenltr 7000 -mintsd 4 -maxtsd 6 -motif TGCA -motifmis 1 -similar 85 -vic 10 -seed 20 -seqids yes)⁹⁸. We used LTR_retriever (v2.9.0, parameters: -threads 16 -noanno)⁹⁹ to remove redundancy from the sequences predicted by the two methods to obtain non-redundant LTR sequences and merged these two sets of sequences to create a *de novo* repeat sequence library. TEclass (v2.1.3, default parameters)¹⁰⁰ was employed to redefine the content labeled as "Unknown". Then, we merged the RepBase library (v20181026) with the *de novo* library. We also used the RepeatMasker subprogram RepeatProteinMask (v4.1.5, parameters: -noLowSimple -pvalue 0.0001) to predict TE proteins. Finally, all repetitive prediction results were merged, and redundancy was removed to obtain the final collection of repetitive sequences in the ‘Kabusa’ haplotype assemblies. Tandem repeat sequences were predicted using TRF software (v4.09, parameters: 2 7 7 80 10 50 2000 -d -h)¹⁰¹ and MISA software (v2.1, default parameters)¹⁰².

708

709

710

711

712

713

714

715

716

717

718

719

720

721

722

723

724

725

726

## 727 Gene annotation

728

729

Protein-coding gene models were predicted by integrating three strategies: *ab initio* prediction, homology-based search, and gene expression-based prediction.
For gene expression-based prediction, (a) NGS RNA-seq clean reads were aligned to the genome using hisat2 (v2.2.1, parameters: -dta -p 16), and the resulting BAM files were used by stringtie (v2.2.1, parameters: -p 16)¹⁰³ to reconstruct the transcripts; (b) ONT Iso-seq data were filtered using NanoFilt (v2.8.0, parameters: -l 50 -q 7) and full-length sequence identification performed with pychopper (v2.7.5, parameters: -t 16). The obtained full-length sequences were aligned to the 'Kabusa' haplotype genomes using minimap2 (v2.26-r1175, parameters: minimap2 -t 16 -ax splice -uf --secondary=no)³⁵. The alignment results in BAM files were used by stringtie (v2.2.1, parameters: -p 16 -R -L)¹⁰³ to reconstruct the transcripts. (c) The transcripts predicted by the two methods were merged using tama (v1.0, parameters: -f filelist.txt -p merge, https://github.com/GenomeRIK/tama). The merged transcript results were then used with TransDecoder (v5.7.0, parameters: TransDecoder.Predict - retain_long_orfs_length 150, TransDecoder.LongOrfs -m 50, https://github.com/TransDecoder/TransDecoder) to predict open reading frames and finally obtain predicted coding genes.

For homology-based search, homologous protein sequences from *D. alata* female genome²⁷, *D. rotundata*³⁶, *D. dumetorum*¹⁰⁴, *D. zingiberensis*¹⁰⁵ and *T. zeylanicus*¹⁰⁶ were aligned to the 'Kabusa' haplotype genomes with tblastn (v2.13.0, parameters: default)¹⁰⁷ to obtain candidate regions. Then, Exonerate (v2.4.0, parameters: -model protein2genome -showtargetgff 1)¹⁰⁸ was used to predict transcripts and coding regions based on the alignment results.

The *ab initio* prediction was conducted using Augustus (v3.5.0, parameters: -uniqueGeneId=true -noInFrameStop=true -gff3=on -strand=both)¹⁰⁹ and glimmerhmm (v3.0.4, parameters: -f -g)¹¹⁰.

All the predicted gene sets were integrated using MAKER2 (v3.01.03, default parameters)¹¹¹. The completeness of the predicted annotation sets was assessed using BUSCO protein mode with the Embryophyta_odb10 data set⁹⁴. Functional annotation of protein-coding genes was performed based on existing databases such as KEGG, GO and Uniprot. Concerning non-coding RNAs prediction, tRNAscan-SE (v2.0.12, parameters: -E -j tRNA.gff -o tRNA.result -f tRNA struct -thread 16)¹¹² was used to find tRNA sequences in the genome, RNAmmer (v1.2, parameters: -S euk -m tsu,lsu,ssu)¹¹³ was used to predict rRNA, and INFERNAL (v1.1.4)¹¹⁴ was employed to search for ncRNA in the genome based on Rfam database (v14)¹¹⁵.
762 Centromere detection analyses

763 We used two complementary approaches to predict centromere regions.

764 1- Tandem repeat density and gene coverage approach: This method leveraged tandem repeat (TR) density and low gene density, as described in previous studies¹¹⁶,¹¹⁷. Using BEDTools (v2.30.0)¹¹⁸, we calculated TR and gene coverage in 50 Kb windows across the ‘Kabusa’ genome, visualizing continuous regions of high TR and low gene density as potential centromeric regions.

769 2- TR cluster (TRC)-based approach: Following the methods used in T2T kiwifruit and *Rhodomyrtus tomentosa* genomes¹¹⁹,¹²⁰, TR sequences were identified using TRF (v4.09.1, parameters: 2 7 7 80 10 50 500 -d -h)¹⁰¹ and clustered into TRCs using cd-hit (v4.8.1, parameters: -c 0.9 -n 5 -d 0 -g 1 -aS 0.9, https://github.com/weizhongli/cdhit/releases). The most abundant TRCs were plotted to aid centromere identification.

775 3- Final prediction and visualization: We used the longest TRC as the representative sequence and minimized overlaps using RepeatMasker. BEDTools (v2.30.0)¹¹⁸ was then used to merge and define centromere regions, which were visualized with RIdeogram¹²¹ and LINKVIEW2 (v1.0.5, https://github.com/YangJianshun/LINKVIEW2), alongside features like rDNA and telomeres.

780 Sequence variation between haplotype genomes of ‘Kabusa’ and evolutionary analysis
781 with other *Dioscorea* species

782 Haplotype genomes of ‘Kabusa’ were aligned using MUMmer (v4.0.0rc1; default parameters)⁹⁶, with HaplotypeA as the reference. Then, variants were detected using SyRI (v1.5; default parameters)³³. For synteny analysis between haplotype genomes of ‘Kabusa’ and the female *D. alata* genome²⁷,²⁸, gene sequences were aligned using LAST (v1170)¹²² to identify similar gene pairs. These gene pairs were then analyzed using JCVI (v0.9.13)¹²³ to determine if they are adjacent on the chromosomes based on annotation files (gff3), thereby identifying all genes in syntenic blocks. The synteny plots were generated using SyRI (v1.5; default parameters)³³.

790 To identify shared and unique genes between the two haplotype genomes (HaplotypeA and HaplotypeB), the genomes were clustered into gene families using Orthofinder (v2.3.12)¹²⁴. Protein sequences from HaplotypeA were extracted and compared to the HaplotypeB genome using tblastn (v2.13.0)¹⁰⁷. Structural integrity of matching genes was verified using gffread (v0.12.7). Protein sequences of structurally complete genes were further compared to HaplotypeA’s proteins using blastp (v2.6.0)¹⁰⁷. Genes with alignment scores >150 were deemed
homologous and removed from HaplotypeA's unique gene list. The remaining genes were considered unique to HaplotypeA. From the unique gene list of HaplotypeA, genes located on Chr6 were identified as DalaChr6A-specific unique genes. The same process was applied to identify unique genes in HaplotypeB and DalaChr6B-specific unique genes.

For phylogenetic tree construction between *Dioscorea* species, a multiple sequence alignment of protein sequences for each single-copy gene family was performed using MUSCLE (v3.8.31)¹²⁵. The alignment results were then filtered using *trimAl* (v1.2rev59; parameters: -gt 0.2)¹²⁶. The filtered alignments were concatenated into the supergene. A maximum likelihood phylogenetic tree was constructed based on the supergene using RAxML (v8.2.10)¹²⁷ with the PROTGAMMAWAG model. Divergence times for the selected species (*D. alata*, *D. rotundata*, *D. tokoro* and *D. zingiberensis*) were estimated using the mcmctree program from the PAML package (v4.9; parameters: nsample=3000000; burnin=8000000; seqtype=0; model=4)¹²⁸.

## Transcriptome analysis

Newly generated raw reads from flower buds, along with short reads from mature tubers of three female *D. alata* accessions ('Dou', 'Kinabayo', and '74F') and three male accessions ('Divin', 'Plimbite', and 'Toufi-Tetea'), previously generated by Mota et al.⁴⁰ (PRJNA918625), were used to perform a comparative gene expression analysis between male and female groups. The package `fastp` (v0.21.0)¹²⁹ was used to filter raw reads (150 bp paired-end) and discard low-quality reads. Several libraries were found to be of poor quality or contaminated after sequencing; therefore, we excluded them from further analyses. The cleaned reads from each sample were mapped to the haploid consensus genome of the 'Kabusa' genome (specifically the HapA assembly, which includes the Y chromosome + the HXR) using STAR (v2.7.9, parameters:default)¹³⁰. This choice was made to enable a focused analysis of male-specific and shared gene expression. The RSEM package¹³¹ was used to calculate gene expression level for each sample expressed as fragments per kilobase of transcript per million fragments mapped. The edgeR (v3.32.1) software¹³² was used for differential expression analysis. Selection thresholds of FDR < 0.05 and |log₂ FC| > 1. Functional enrichment analysis was performed using the R package `clusterProfiler` v4.0¹³³ based on the functional annotations of all genes in this species, including Gene Ontology¹³⁴ and Kyoto Encyclopedia of Genes and Genomes¹³⁵.
Refining sex-determination region

Whole genome resequencing raw data of 20 male and 17 female accessions were downloaded from NCBI (PRJNA918625). We constructed a male-pool by merging male data and a female-pool by merging female data. The two pools were mapped to the HaplotypeA of the ‘Kabusa’ genome. GATK (v4.6.1.0)¹³⁶ was used for SNP calling and *Fst* was estimated between the two pools with VCFtools (0.1.17)¹³⁷. Next, SNPs that are homozygous in the female pool but heterozygous in the two pools were filtered. The variant density in a 100 Kb window was calculated and plotted with the circlize package¹³⁸ in R (v4.4.2).

Genome wide association studies

A total of 200 *D. alata* accessions were planted at Roujol and Godet research stations in Guadeloupe during seasons 2023 and 2024. Sex phenotypes were recorded during flowering period (September to December) each year. Flowering was recorded in 108 accessions including 68 males and 40 females. Genotyping data from the 108 *D. alata* accessions were retrieved from the population studied by Mota et al.⁴⁰ (PRJNA918625). SNPs and InDels calls were performed as described above using GATK (v4.6.1.0)¹³⁶. GWAS was conducted independently for SNPs and InDels genotype matrices using PLINK v1.90b6.16¹³⁹ with the logistic regression model that included the top 3 principal components from SNPs and InDels genotype matrices. The Manhattan plots were generated in R (v4.4.2). Markers having a significant association with traits were determined by the adjusted *p* value.

K-mer analysis

To identify sex-specific sequences based on k-mer composition, we performed a k-mer counting analysis using KMC3 (v3.2.1)¹⁴⁰. Cleaned Illumina paired-end reads from 20 male and 20 female individuals were used for this analysis. For each sample, canonical 27-mers (k = 27) were counted independently using the following parameters: -k27 -ci5 -fq -m64 -t4. The minimum count threshold (-ci5) was applied to filter out low-frequency k-mers likely resulting from sequencing errors. Reads were processed in paired-end mode (-fq) using compressed FASTQ files (\*\_clean_R1.fq.gz and \*\_clean_R2.fq.gz). The k-mer databases generated for each individual were stored separately to allow downstream comparison between male and female groups. (1) Canonical k-mer databases previously generated using KMC3 were compared using kmc_tools. Specifically, k-mers unique to males were identified by subtracting male using kmers_subtract producing male-specific k-mers. (2) The unique k-mers were dumped into text files and filtered to exclude sequences containing ambiguous nucleotides
858 ('N'). Cleaned lists were then converted into FASTA format using awk, with each k-mer
859 assigned a unique identifier. (3) Male-specific k-mers were mapped to the *D. alata* HaplotypeA
860 using bwa mem, with relaxed parameters (-k 21 -T 21 -a) to accommodate short reads and allow
861 multi-mapping. The resulting SAM file was converted to sorted BAM format using SAMtools,
862 (v1.18)¹⁴¹ followed by indexing for downstream visualization and coverage analysis.

863 Read coverage analysis

864 To quantify sequencing read coverage, we used Mosdepth (v0.3.3)¹⁴². A 10 Kb sliding
865 window BED file was generated to cover DalaChr6A. BAM files were grouped by sex based
866 on metadata and processed independently for the 68 male and 40 female samples. Coverage
867 depth was calculated across each window using Mosdepth with the options -n -x to exclude
868 duplicate and secondary alignments. For each window, the mean depth was calculated across
869 all individuals of each sex using custom AWK scripts. The final dataset consisted of a merged
870 table reporting average male and female coverage per window across the DalaChr6A.

871 Sex chromosome divergence analysis

872 We estimated nonsynonymous (Ka) and synonymous (Ks) substitution rates for gene
873 pairs located on the X and Y chromosomes. First, we extracted gene coordinates from specific
874 regions of interest SDR, INV1-4 and XSR, PAR, MSY—using bedtools intersect and GFF
875 annotations from the ‘Kabusa’ haplotype-resolved genome assembly. The total number of
876 genes in each region was quantified based on gene annotations from chromosomes DalaChr6A
877 (Y) and DalaChr6B (X). Next, their corresponding coding DNA sequences and protein
878 sequences were retrievd using seqtk-1.5 (r133; https://github.com/lh3/seqtk). Orthologous gene
879 pairs between Haplotypes A and B were determined using a anchor file from whole-genome
880 synteny analysis. For each orthologous gene pair, protein sequences were aligned using MAFFT
881 (v7.453)¹⁴³. Codon alignments were then generated using PAL2NAL (v14)¹⁴⁴. doi:
882 10.1093/nar/gkl315.). The resulting codon alignments were cleaned, concatenated, and used as
883 input to KaKs_Calcula tor (v2.0)¹⁴⁵ with the Yang-Nielsen model to estimate pairwise Ka, Ks,
884 and Ka/Ks ratios.

885 PACE marker development and experimental procedure

886 SNP genotyping was performed as described by von Maydell¹⁴⁶. DNA was extracted
887 from leaf tissue 146 *D. alata* accessions using mixed alkyl trimethylammonium bromide buffer
and NucleoMag Plant Kit (Macherey-Nagel, Düren, Germany). We used the PACE™ master mix (3CR Bioscience, UK) (standard protocol), and runs were performed on LightCycler480 (Roche, Basel, Switzerland). Five negative controls (water) were included in the experiment. From the raw fluorescence data, allele calling was performed by defining fluorescence clusters. Plotting was conducted with the R package snpclust (v1.1; https://github.com/jframi/snpclust).

## 893 Methylation analysis

HiFi CCS reads were used for estimating the methylation frequencies of the CpGs according to the ccsmethphase pipeline developed by Ni et al.¹⁴⁷. To phase the reads, we used Meryl (v1.4.1)⁹⁵ to find unique k-mers in the reference genome of each haplotype as described by Nie et al.¹⁴⁸ and classify the reads using the unique k-mers. The methylation frequencies were estimated for binned regions (100 Kb) in each haplotype assembly.

## 899 Statistical analysis

All statistical analyses were performed in R (v4.4.2). Depending on the data distribution and experimental design, both parametric tests (t-tests and one-way ANOVA) and non-parametric tests (Wilcoxon rank-sum test and Kruskal-Wallis test) were applied. A significance threshold of *p* < 0.05 was used throughout.

## 904 Supplementary information

906
**Supplementary Table 1:** Statistics of the Illumina short reads data
908
**Supplementary Table 2:** Statistics of k-mer analysis based on Illumina short reads data
909
**Supplementary Table 3:** Oxford Nanopore Technology ultra-long data statistics
910
**Supplementary Table 4:** Pacbio HiFi ultra-long data statistics
911
**Supplementary Table 5:** Statistics of the Hi-C data
912
**Supplementary Table 6:** Four assembly strategies evaluated on the quality of the contig-level genomes
913
**Supplementary Table 7:** Detailed statistical results of each chromosome and its length in the hybrid assembly version (HiFi + Hi-C)
914
**Supplementary Table 8:** Position of telomeric repeat sequences along the chromosomes in the 'Kabusa' genome
915
**Supplementary Table 9:** Gap filling in the genome assembly
916
**Supplementary Table 10:** Detailed statistical results of each chromosome and its length in the final near-complete 'Kabusa' genome assembly
917
**Supplementary Table 11:** Statistics of repetitive element annotation in the genome of 'Kabusa'
918
**Supplementary Table 12:** Statistics of ONT based full-length transcriptome data for genome annotation
919
**Supplementary Table 13:** Statistics of prediction, functional annotation of coding genes and Non-coding RNA annotation in the HaplotypeA and HaplotypeB genome assemblies of 'Kabusa'
927 **Supplementary Table 14:** Statistics of gene set BUSCO evaluation
928 (busco_lineage:embryophyta_odb10) in the HaplotypeA and HaplotypeB genome assemblies
929 of 'Kabusa'

930 **Supplementary Table 15:** Position of centromeric repeat sequences along the chromosomes
931 in the 'Kabusa' genome

932 **Supplementary Table 16:** Report of the Liftoff annotation of the female genome *D. alata_v2*
933 based on HaplotypeA of the 'Kabusa' genome

934 **Supplementary Table 17:** Statistics of variants detected between the two haplotype genome
935 assemblies of 'Kabusa' and their impact on gene function

936 **Supplementary Table 18:** PLINK logistic GWAS results based on SNPs and Indels markers.
937 102 *D. alata* accessions were evaluated for sex phenotype in the field over 2 years at 2 sites in
938 Guadeloupe.

939 **Supplementary Table 19:** Genes located within the SDR and HXR and their functional
940 annotation. INV1-INV4 = Inversion 1-Inversion 4. XSR = X-specific region; MSY = Male-
941 Specific region of the Y chromosome

942 **Supplementary Table 20:** Sex-linked genes that have functional homologues on autosomal
943 chromosomes in the *D. alata* 'Kabusa' genome

944 **Supplementary Table 21:** Samples used for RNA-seq transcriptome

945 **Supplementary Table 22:** Statistics of RNA-seq transcriptome data from male and female
946 flower buds. -1, -2, -3 represent the three biological replicates. If a genotype does not have 3
947 replicates, it means a replicate was contaminated and was excluded.

948 **Supplementary Table 23:** Sex-biased genes located within the SDR and HXR in *D. alata*.
949 Up/down expression in male compared to female flower buds

950 **Supplementary Table 24:** Development of a PACE sda6 marker linked to sex determination
951 in Dioscorea alata. A sex-linked locus was detected via GWAS in the gene DalaChr6AG083080
952 at the position 91 bp. All female offspring are homozygous TT (XX) while male are
953 heterozygous GT (XY).

954 **Supplementary Table 25:** Details on the 146 *D. alata* germplasm genotyped with the sda6
955 PACE marker

956 **Supplementary Table 26:** Details on the *D. alata* germplasm used for genetic analyses

957 **Supplementary Fig. 1:** Flow cytometry-based determination of the ploidy level in 'Kabusa'.
958 The ploidy level was assessed by comparing the position of the main peak with the references
959 indicated by the arrows. According to the results, 'Kabusa' is identified as a diploid cultivar
960 (2x).

961 **Supplementary Fig. 2:** k-mer analysis for estimation of the 'Kabusa' genome size. The genome
962 size was estimated to approximately 494.13 Mb.

963 **Supplementary Fig. 3:** Distribution of the length of gene and gene components in the 'Kabusa'
964 genome assembly compared to related species. A. Gene length distribution. B. Coding sequence
965 (CDS) length distribution. C. Exon length distribution. D. Intron length distribution.

967 **Supplementary Fig. 4:** Centromere and telomere detection map in the two haplotype genome
968 assemblies of 'Kabusa'.
969 Orange squares and green triangles represent centromere and telomeres within the assembled
970 genome. Blue indicates low gene density; red indicates
971 high gene density. Deep blue and light blue circles represent the distribution along the
972 chromosomes of 45S_rDNA and 5S_rDNA, respectively.

973 **Supplementary Fig. 5:** Illustrations of alignments of preliminary assemblies (A and B) or ONT
974 reads (C) that clearly spanned gap regions providing strong evidence of the accuracy of the gap
975 filling stage.
976 **Supplementary Fig. 6:** Centromere detection in the two haplotype genome assemblies of
977 'Kabusa'. The first layer representing LTR/Copia density (blue color), the second layer
978 representing LTR/Gypsy density (orange color), the third layer representing tandem repeat
979 density (ruby red color) and the fourth layer representing the gene density (turquoise color)
980 along the chromosomes.

981 **Supplementary Fig. 7:** Structural variation (SV) between the two haplotype assemblies of the
982 'Kabusa' genome. A. Quantity of SVs mapped to different genomic regions. B. Gene ontology
983 enrichment analysis of SV-influenced genes.

984 **Supplementary Fig. 8:** Collinearity analysis between the HaplotypeA, HaplotypeB of the
985 'Kabusa' genome and the female genome assembly from Bredeson et al. [27].

986 **Supplementary Fig. 9:** Intergenomic comparison of the HaplotypeA, HaplotypeB of the
987 'Kabusa' genome and the female genome assembly from Bredeson et al. [27]. Red circle shows
988 the location of the pericentric inversion in DalaChr6A vs DalaChr6F but absent in DalaChr6B
989 vs DalaChr6F.

990 **Supplementary Fig. 10:** IGV snapshots showing the alignment of HiFi reads spanning the
991 breakpoints of the inverted region bewteen the Y chromosome (DalaChr6A) and X
992 chromosome (DalaChr6B) of the 'Kabusa' genome

993 **Supplementary Fig. 11:** Sex chromosome evolution in yams. (A) Synteny analysis between
994 genomes of four yam species with emphasis on sex chromosomes reported to date. Dala_V2 =
995 female genome of *Dioscorea alata*, DalA and DalB = Two haplotype assemblies of the
996 'Kabusa' genome (*D. alata*), Drot = *D. rotundata*, Dtok = *D. tokoro* and Dzin = *D.
997 zingiberensis*. *Dioscorea alata* and *D. rotundata* belong to the section Enantiophyllum while
998 *D. tokoro* and *D. zingiberensis* belong to the section Stenophora. XY and ZW denote the sex
999 determination systems reported in the yam species. Numbers before and within brackets
1000 represent divergence time million years ago

1001 **Supplementary Fig. 12:** Distribution of transposable elements in the two haplotype assemblies
1002 of the 'Kabusa' genome.

1003 **Supplementary Fig. 13:** Gene expression profiling and heatmap clustering of male and female
groups.

1005 **Supplementary Fig. 14:** Gene ontology enrichment analysis of the 6917 sex-biased expressed
genes detected between male and female groups in *Dioscorea alata*.

1007 **Supplementary Fig. 15:** Patterns of methylation level (DNA methylation 5mC at CpG sites)
1008 along homologous chromosomes in the two haplotype assemblies of the 'Kabusa' genome. Plot
1009 was generated using a window of 100 Kb.

1010 **Supplementary Fig. 16:** Genotyping of *D. alata* diversity panel (146 accessions from different
countries) with the *sda6* PACE marker. Fluorescence signals are plotted by accession and
observed sex. In x-axis, the "T" fluorescence allele (X), and in y-axis, the "G" fluorescence
allele (Y).

1014
1015 **Code**
1016 No new code has been developed in this project.

1017
1018 **Data availability**

1019 All raw high-throughput sequencing data comprising Illumina, Oxford Nanopore, PacBio, and
1020 Hi-C reads are deposited in NCBI under project numbers PRJNA1175864 and PRJNA880983.

1021 Whole genome assemblies are available under the Umbrella BioProject number
1022 PRJNA1182344 with accession numbers PRJNA1181019 (HaplotypeA) and PRJNA1181020
1023 (HaplotypeB). The genome assemblies and annotation files are also available at https://yam-genome-hub.southgreen.fr.

1024
1025
1026 Funding
1027 This work was supported by the grants from CGIAR Research Program on Roots, Tubers and Bananas (CRP-RTB), Bill and Melinda Gates Foundation (BMGF) through the RTB breeding project (grant opportunity INV-003446, formerly OPP1112307), the French National Research Agency (ANR) project EVOSEXYAM (ANR-24-CE20-6201-01), the ERDF program Johannes Amos Comenius project TowArds Next GENeration Crops (CZ.02.01.01/00/22_008/0004581) and European Union and Guadeloupe Region (Programme Opérationnel FEDER – Guadeloupe - Conseil Régional 2021-2027) project: TAESA.

1034
1035 Authors' contributions

1036 Conceptualization: Komivi Dossa and Hanâ Chair; Methodology: Komivi Dossa, Ting Qiu, Yao Chen, Peng Ni, Xiaojun Su, Gang Chen, Petre I. Dobrev, Bertrand Fouks, Gemma Arnau, Pierre Mournet, Jianxin Wang, Sylvain Santoni, Hanâ Chair; Software: Komivi Dossa, Ting Qiu, Yao Chen, Peng Ni, Xiaojun Su, Gang Chen, Jianxin Wang; Formal Analysis: Yedomon Ange Bovys Zoclanclounon, Muhammad Amjad Nawaz, Erick Malédon, Moussa Diouf, Petre I. Dobrev, Benjamin Karikari, Bertrand Fouks, Pierre Mournet; Investigation: Erick Malédon, Marie-Claire Gravillon, Christophe Perrot, Elie Nudol, Lévy Laurent, Pierre Mournet, Sylvain Santoni; Resources: Komivi Dossa and Hanâ Chair; Data Curation: Komivi Dossa, Ting Qiu, Yedomon Ange Bovys Zoclanclounon, Pierre Mournet, Yao Chen, Moussa Diouf, Peng Ni, Xiaojun Su, Gang Chen, Petre I. Dobrev, Bertrand Fouks, Jianxin Wang; Writing—Original Draft Preparation: Komivi Dossa, Ting Qiu, Yedomon Ange Bovys Zoclanclounon, Muhammad Amjad Nawaz; Writing—Review and Editing: Hanâ Chair, Gemma Arnau; Visualization, Komivi Dossa, Ting Qiu, Gang Chen, Muhammad Amjad Nawaz; Supervision: Hanâ Chair, Benjamin Karikari; Project Administration: Komivi Dossa and Hanâ Chair; Funding Acquisition: Komivi Dossa and Hanâ Chair. All authors have read and agreed to the published version of the manuscript.

1052
1053 Acknowledgment
1054 We thank Dr Song Chi for providing access to his computing platform.

1055
1056 Declarations
1057
1058 **Ethics approval and consent to participate**

1059 Not applicable.

1060

1061 **Consent for publication**

1062 Not applicable.

1063

1064 **Competing interests**

1065 The authors declare that they have no competing interests.

1066

1067 **References**

1068 1. Renner, S.S. The relative and absolute frequencies of angiosperm sexual systems: dioecy, monoecy, gynodioecy, and an updated online database. *Am. J. Bot.* **101**, 1588-96 (2014).

1069 2. Barrett, S.C. The evolution of plant reproductive systems: how often are transitions irreversible? *Proc. Biol. Sci.* **280**, 20130913 (2013).

1070 3. Liu, Z. et al. A primitive Y chromosome in papaya marks incipient sex chromosome evolution. *Nature* **427**, 348-52 (2004).

1071 4. Akagi, T. et al. Two Y-chromosome-encoded genes determine sex in kiwifruit. *Nat. Plants* **5**, 801-809 (2019).

1072 5. Akagi, T. et al. Recurrent neo-sex chromosome evolution in kiwifruit. *Nat. Plants* **9**, 393-402 (2023).

1073 6. Torres, M.F. et al. Genus-wide sequencing supports a two-locus model for sex-determination in Phoenix. *Nat. Commun.* **9**, 3969 (2018).

1074 7. Akagi, T. et al. A Y-Encoded Suppressor of Feminization Arose via Lineage-Specific Duplication of a Cytokinin Response Regulator in Kiwifruit. *Plant Cell* **30**, 780-795 (2018).

1075 8. Akagi, T. et al. The persimmon genome reveals clues to the evolution of a lineage-specific sex determination system in plants. *PLoS Genet.* **16**, e1008566 (2020).

1076 9. Müller, N.A. et al. A single gene underlies the dynamic evolution of poplar sex determination. *Nat. Plants* **6**, 630-637 (2020).

1077 10. Zhang, X. et al. A convergent mechanism of sex determination in dioecious plants: Distinct sex-determining genes display converged regulation on floral B-class genes. *Front. Plant Sci.* **13**, 953445 (2022).

1078 11. Garcia, S. et al. Sex-chrom v. 2.0: a database of green plant species with sex chromosomes. *Chromosoma* **132**, 55-58 (2023).

1079 12. Baránková, S. et al. Sex-chrom, a database on plant sex chromosomes. *New Phytol.* **227**, 1594-1604 (2020).

1080 13. Garg, V. et al. Unlocking plant genetics with telomere-to-telomere genome assemblies. *Nat. Genet.* **56**, 1788-1799 (2024).

1081 14. Liao, Z. et al. A telomere-to-telomere reference genome of ficus (*Ficus hispida*) provides new insights into sex determination. *Hortic. Res.* **11**, uhad257 (2024).

1082 15. Liu, A.W. et al. H3S28P Antibody Staining of Okinawan Oikopleura dioica Suggests the Presence of Three Chromosomes. *F1000Res.* **9**, 780 (2020).

1083 16. Massonnet, M. et al. The genetic basis of sex determination in grapes. *Nat. Commun.* **11**, 2902 (2020).

1084 17. Chen, J. et al. A complete telomere-to-telomere assembly of the maize genome. *Nat. Genet.* **55**, 1221-1231 (2023).

1085 18. Zhang, X. et al. Unzipping haplotypes in diploid and polyploid genomes. *Comput. Struct. Biotechnol. J.* **18**, 66-72 (2020).

1086 19. Li, H. and Durbin R. Genome assembly in the telomere-to-telomere era. *Nat. Rev. Genet.* **25**, 658-670 (2024).

1087 20. Owusu Danquah, E. et al. Sustainable Intensification and Climate-Smart Yam Production for Improved Food Security in West Africa: A Review. *Front. Agron.* **4** (2022).

1088 21. Neina, D. Ecological and Edaphic Drivers of Yam Production in West Africa. *Appl. Environ. Soil. Sci.* **2021**, 5019481 (2021).

1089 22. Arnau, G. et al. Revision of ploidy status of *Dioscorea alata* L. (*Dioscoreaceae*) by cytogenetic and microsatellite segregation analysis. *Theor. Appl. Genet.* **118**, 1239-49 (2009).

1090 23. Sugihara, Y. et al. Population Genomics of Yams: Evolution and Domestication of *Dioscorea* Species, in Population Genomics: Crop Plants, O. Rajora, Editor. Springer International Publishing: Cham. 837-864 (2024).

1091 24. Mondo, J.M. et al. Genome-Wide Association Studies for Sex Determination and Cross-Compatible in Water Yam (*Dioscorea alata* L.). *Plants* **10** (2021).
<table>
  <tr>
    <td>1116</td>
    <td>25.</td>
    <td>Cormier, F. et al. A reference high-density genetic map of greater yam (Dioscorea alata L.). Theor. Appl. Genet. 132, 1733-1744. (2019)</td>
  </tr>
  <tr>
    <td>1117</td>
    <td>26.</td>
    <td>Cormier, F. et al. Genetic control of flowering in greater yam (Dioscorea alata L.). BMC Plant Biol. 21, 163 (2021).</td>
  </tr>
  <tr>
    <td>1118</td>
    <td>27.</td>
    <td>Bredeson, J.V. et al. Chromosome evolution and the genetic basis of agronomically important traits in greater yam. Nat. Commun. 13, 2001 (2022).</td>
  </tr>
  <tr>
    <td>1119</td>
    <td>28.</td>
    <td>Zhang YM, et al. A telomere-to-telomere genome assembly for greater yam (Dioscorea alata). Plant Commun. 2, 101326 (2025).</td>
  </tr>
  <tr>
    <td>1120</td>
    <td>29.</td>
    <td>Ehounou, A.E. et al. Identification and validation of QTLs for tuber quality related traits in greater yam (Dioscorea alata L.). Sci. Rep. 12, 8423 (2022).</td>
  </tr>
  <tr>
    <td>1121</td>
    <td>30.</td>
    <td>Cheng, H. et al. Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nat. Method 18, 170-175 (2021).</td>
  </tr>
  <tr>
    <td>1122</td>
    <td>31.</td>
    <td>Hu, J. et al. NextDenovo: an efficient error correction and accurate assembly tool for noisy long reads. Genome Biol. 25, 107 (2024).</td>
  </tr>
  <tr>
    <td>1123</td>
    <td>32.</td>
    <td>Durand, N.C. et al. Juicebox Provides a Visualization System for Hi-C Contact Maps with Unlimited Zoom. Cell Syst. 3, 99-101 (2016).</td>
  </tr>
  <tr>
    <td>1124</td>
    <td>33.</td>
    <td>Goel, M. et al. SyRI: finding genomic rearrangements and local sequence differences from whole-genome assemblies. Genome Biol. 20, 277 (2019).</td>
  </tr>
  <tr>
    <td>1125</td>
    <td>34.</td>
    <td>Shumate A, Salzberg SL. Liftoff: accurate mapping of gene annotations. Bioinformatics 19, 1639-1643 (2021).</td>
  </tr>
  <tr>
    <td>1126</td>
    <td>35.</td>
    <td>Li, H. Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences. Bioinformatics 32, 2103-10 (2016).</td>
  </tr>
  <tr>
    <td>1127</td>
    <td>36.</td>
    <td>Tamiru, M. et al. Genome sequencing of the staple food crop white Guinea yam enables the development of a molecular marker for sex determination. BMC Biol. 15, 86 (2017).</td>
  </tr>
  <tr>
    <td>1128</td>
    <td>37.</td>
    <td>Natsume, S. et al. Whole genome sequencing of a wild yam species Dioscorea tokoro reveals a genomic region associated with sex. BioRxiv 2022.06.11.495741 (2022).</td>
  </tr>
  <tr>
    <td>1129</td>
    <td>38.</td>
    <td>Ried, T. et al. Genome Sequencing of Chinese yam (Dioscorea polystachya), analysis of PEBP gene family diversity and identification of a potential tuber inducing factor. BioRxiv 2025.03.10.641988 (2025).</td>
  </tr>
  <tr>
    <td>1130</td>
    <td>39.</td>
    <td>Hu K, Feng Y, Li P, Chen M, Shen ZJ, Sun XQ, Lu RS. Haplotype-resolved genome and population genomics provide insights into dioscin biosynthesis and evolutionary history of the medicinal species Dioscorea nipponica. Plant J. 121, e17237 (2025).</td>
  </tr>
  <tr>
    <td>1131</td>
    <td>40.</td>
    <td>Mota, A.Z. et al. Whole-genome sequencing and comparative genomics reveal candidate genes associated with quality traits in Dioscorea alata. BMC Genomics 25, 248 (2024).</td>
  </tr>
  <tr>
    <td>1132</td>
    <td>41.</td>
    <td>Iocco-Corena, et al. VviPLATZ1 is a major factor that controls female flower morphology determination in grapevine. Nat. Commun. 12, 6995 (2021).</td>
  </tr>
  <tr>
    <td>1133</td>
    <td>42.</td>
    <td>Zhang, S. et al. Chromosome-scale assemblies of the male and female Populus euphratica genomes reveal the molecular basis of sex determination and sexual dimorphism. Commun. Biol. 5, 1186 (2022).</td>
  </tr>
  <tr>
    <td>1134</td>
    <td>43.</td>
    <td>Srikulnath, K., Ahmad, S.F., Singchat, W. & Panthum, T. Do Ty3/Gypsy Transposable Elements Play Preferential Roles in Sex Chromosome Differentiation? Life 12, 522 (2022).</td>
  </tr>
  <tr>
    <td>1135</td>
    <td>44.</td>
    <td>Hobza, R. et al. Impact of Repetitive Elements on the Y Chromosome Formation in Plants. Genes 8, (2017).</td>
  </tr>
  <tr>
    <td>1136</td>
    <td>45.</td>
    <td>Wang, J. et al. Sequencing papaya X and Yh chromosomes reveals molecular basis of incipient sex chromosome evolution. Proc. Natl. Acad. Sci. U.S.A. 109, 13710-13715 (2012).</td>
  </tr>
  <tr>
    <td>1137</td>
    <td>46.</td>
    <td>Moraga, C. et al. The Silene latifolia genome and its giant Y chromosome. Science 7, 630-636 (2025).</td>
  </tr>
  <tr>
    <td>1138</td>
    <td>47.</td>
    <td>Prentout, D. et al. An efficient RNA-seq-based segregation analysis identifies the sex chromosomes of Cannabis sativa. Genome Res. 30, 164-172 (2020).</td>
  </tr>
  <tr>
    <td>1139</td>
    <td>48.</td>
    <td>Cossard G.G., Toups M.A. & Pannell J.R. Sexual dimorphism and rapid turnover in gene expression in pre-reproductive seedlings of a dioecious herb. Ann. Bot. 8, 1119-1131. (2019).</td>
  </tr>
  <tr>
    <td>1140</td>
    <td>49.</td>
    <td>DePaoli, H.C. et al. Stigma/style cell cycle inhibitor 1 (SCI1), a tissue-specific cell cycle regulator that controls upper pistil development. New Phytol. 190, 882-895 (2011).</td>
  </tr>
  <tr>
    <td>1141</td>
    <td>50.</td>
    <td>DePaoli, H.C., Dornelas, M.C. & Goldman, M.H.S. SCI1 is a component of the auxin-dependent control of cell proliferation in Arabidopsis upper pistil. Plant Sci. 229, 122-130 (2014).</td>
  </tr>
  <tr>
    <td>1142</td>
    <td>51.</td>
    <td>He, Y. et al. Functional characterization of Arabidopsis thaliana isopropylmalate dehydrogenases reveals their important roles in gametophyte development. New Phytol. 189, 160-75 (2011).</td>
  </tr>
  <tr>
    <td>1143</td>
    <td>52.</td>
    <td>Huang, Y. et al. Epigenetic regulation of female germline development through ERECTA signaling pathway. New Phytol. 240, 1015-1033 (2023).</td>
  </tr>
  <tr>
    <td>1144</td>
    <td>53.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1145</td>
    <td>54.</td>
    <td>Luo, Y. et al. Developmental basis for flower sex determination and effects of cytokinin on sex determination in Plukenetia volubilis (Euphorbiaceae). Plant Reprod. 33, 21-34 (2020).</td>
  </tr>
  <tr>
    <td>1146</td>
    <td>55.</td>
    <td>He, Y. et al. Functional characterization of Arabidopsis thaliana isopropylmalate dehydrogenases reveals their important roles in gametophyte development. New Phytol. 189, 160-75 (2011).</td>
  </tr>
  <tr>
    <td>1147</td>
    <td>56.</td>
    <td>Huang, Y. et al. Epigenetic regulation of female germline development through ERECTA signaling pathway. New Phytol. 240, 1015-1033 (2023).</td>
  </tr>
  <tr>
    <td>1148</td>
    <td>57.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1149</td>
    <td>58.</td>
    <td>Luo, Y. et al. Developmental basis for flower sex determination and effects of cytokinin on sex determination in Plukenetia volubilis (Euphorbiaceae). Plant Reprod. 33, 21-34 (2020).</td>
  </tr>
  <tr>
    <td>1150</td>
    <td>59.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1151</td>
    <td>60.</td>
    <td>Luo, Y. et al. Developmental basis for flower sex determination and effects of cytokinin on sex determination in Plukenetia volubilis (Euphorbiaceae). Plant Reprod. 33, 21-34 (2020).</td>
  </tr>
  <tr>
    <td>1152</td>
    <td>61.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1153</td>
    <td>62.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1154</td>
    <td>63.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1155</td>
    <td>64.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1156</td>
    <td>65.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1157</td>
    <td>66.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1158</td>
    <td>67.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1159</td>
    <td>68.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1160</td>
    <td>69.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1161</td>
    <td>70.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1162</td>
    <td>71.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1163</td>
    <td>72.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1164</td>
    <td>73.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1165</td>
    <td>74.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1166</td>
    <td>75.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1167</td>
    <td>76.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1168</td>
    <td>77.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1169</td>
    <td>78.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1170</td>
    <td>79.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1171</td>
    <td>80.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1172</td>
    <td>81.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1173</td>
    <td>82.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1174</td>
    <td>83.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
  <tr>
    <td>1175</td>
    <td>84.</td>
    <td>Siriwardana, C.L. et al. NUCLEAR FACTOR Y, Subunit A (NF-YA) Proteins Positively Regulate Flowering and Act Through FLOWERING LOCUS T. PLoS Genet. 12, e1006496 (2016).</td>
  </tr>
</table>
1176 55. Kanaoka, M.M. et al. KOMPEITO, an Atypical *Arabidopsis* Rhomboid-Related Gene, Is Required for Callose Accumulation and Pollen Wall Development. *Int. J. Mol. Sci.* **23**, (2022).
1177
1178 56. Thu, S.W. et al. Mutation in a PHD-finger protein MS4 causes male sterility in soybean. *BMC Plant Biol.* **19**, 378 (2019).
1179
1180 57. Yang X. Makaroff C. A. & Ma H. The Arabidopsis MALE MEIOCYTE DEATH1 gene encodes a PHD-finger protein that is required for male meiosis. *Plant Cell* **15**, 1281-1295 (2003).
1181
1182 58. Yang Z. F. et al. TDR INTERACTING PROTEIN 3, encoding a PHD-finger transcription factor, regulates Ubisch bodies and pollen wall formation in rice. *Plant J.* **99**, 844-861 (2019).
1183
1184 59. Wasternack, C. Determination of sex by jasmonate. *J. Integr. Plant Biol.* **62**, 162-164 (2020).
1185 60. Xie, D.X. et al. COII: an *Arabidopsis* gene required for jasmonate-regulated defense and fertility. *Science* **1280**, 1091-4 (1998).
1186
1187 61. Pagnussat, G.C. et al. Genetic and molecular identification of genes required for female gametophyte development and function in Arabidopsis. *Development* **132**, 603-614 (2005).
1188
1189 62. Fernandes, J.B. et al. Structural variation and DNA methylation shape the centromere-proximal meiotic crossover landscape in Arabidopsis. *Genome Biol.* **22**, 30 (2024).
1190
1191 63. Li, Q. et al. Haplotype-resolved T2T genome assemblies and pangenome graph of pear reveal diverse patterns of allele-specific expression and the genomic basis of fruit quality traits. *Plant Commun.* **5**, 101000 (2024).
1192
1193
1194 64. Han, L. et al. A haplotype-resolved genome provides insight into allele-specific expression in wild walnut (* Juglans regia* L.). *Sci. Data* **11**, 278 (2024).
1195
1196 65. Raiyemo, D.A. et al. A phased chromosome-level genome assembly provides insights into the evolution of sex chromosomes in *Amaranthus tuberculatus*. *BioRxiv* 2024.05.30.596720 (2024).
1197
1198 66. Lu, M. & He, X. Centromere repositioning causes inversion of meiosis and generates a reproductive barrier. *Proc. Natl. Acad. Sci. U S A* **116**, 21580-21591 (2019).
1199
1200 67. Han, Y. et al. Centromere repositioning in cucurbit species: implication of the genomic impact from centromere activation and inactivation. *Proc. Natl. Acad. Sci. U S A* **106**, 14937-41 (2009).
1201
1202 68. Gent, J.I. et al. Stable Patterns of CENH3 Occupancy Through Maize Lineages Containing Genetically Similar Centromeres. *Genetics* **200**, 1105-16 (2015).
1203
1204 69. Taagen, E., Bogdanove, A.J. & Sorrells, M.E. Counting on Crossovers: Controlled Recombination for Plant Breeding. *Trends Plant Sci.* **25**, 455-465 (2020).
1205
1206 70. Hou, J. et al. Different autosomes evolved into sex chromosomes in the sister genera of Salix and Populus. *Sci. Rep.* **5**, 9076 (2015).
1207
1208 71. Cauret, C.M.S. et al. Chromosome-scale assembly with a phased sex-determining region resolves features of early Z and W chromosome differentiation in a wild octoploid strawberry. *G3* **12** (2022).
1209
1210 72. Martin, H. et al. Evolution of Young Sex Chromosomes in Two Dioecious Sister Plant Species with Distinct Sex Determination Systems. *Genome Biol. Evol.* **11**, 350-361 (2019).
1211
1212 73. Pucholt, P. et al. Recent Sex Chromosome Divergence despite Ancient Dioecy in the Willow *Salix viminalis*. *Mol. Biol. Evol.* **34**, 1991-2001 (2017).
1213
1214 74. Charlesworth, D. and A. Harkess, Why should we study plant sex chromosomes? *Plant Cell*, 2024. 36(5, 1242-1256.
1215
1216 75. Dossa, K. et al. Genome-wide association studies reveal novel loci controlling tuber flesh color and oxidative browning in *Dioscorea alata*. *J. Sci. Food Agric.* **104**, 4895-4906 (2024).
1217
1218 76. Zhou, Q. et al. Complex evolutionary trajectories of sex chromosomes across bird taxa. *Science* **346**, 1246338.
1219
1220 77. Denadi, N. Déterminisme du sexe chez les ignames cultivées *Dioscorea rotundata* (Poir.) du Bénin. PhD Thesis, UCLouvain, Belgium. pp : 1-364. available at http://hdl.handle.net/2078.1/270468. (2022).
1221
1222 78. Denadi, N. et al. Plant Sex Prediction Using Genetic Markers in Cultivated Yams (*Dioscorea rotundata* Poir.) in Benin. *Agronomy* **10** (2020).
1223
1224 79. Charlesworth, D. et al. Steps in the evolution of heteromorphic sex chromosomes. *Heredity* **95**, 118-128 (2005).
1225
1226 80. Hoffmann, A.A. & Rieseberg L.H. Revisiting the Impact of Inversions in Evolution: From Population Genetic Markers to Drivers of Adaptive Shifts and Speciation? *Annu. Rev. Ecol. Evol. Syst.* **39**, 21-42 (2008).
1227
1228
1229 81. She, H. et al. Evolution of the spinach sex-linked region within a rarely recombining pericentromeric region. *Plant Physiol.* **193**, 1263-1280 (2023).
1230
1231 82. Wang, Y. et al. Evolution of Sex-linked Genes and the Role of Pericentromeric Regions in Sex Chromosomes: Insights from Diploid Willows. *Mol. Biol. Evol.* **41** (2024).
1232
1233 83. Olito, C. and J.K. Abbott, The evolution of suppressed recombination between sex chromosomes and the lengths of evolutionary strata. *Evolution* **77**, 1077-1090 (2023).
1234
1235 84. Renner, S.S. & Müller N.A. Plant sex chromosomes defy evolutionary models of expanding
<table>
  <tr>
    <td>1236</td>
    <td>recombination suppression and genetic degeneration. Nat. Plants 7, 392-402 (2021).</td>
  </tr>
  <tr>
    <td>1237</td>
    <td>85. Golenberg, E.M. & West N.W. Hormonal interactions and gene regulation can link monoecy and environmental plasticity to the evolution of dioecy in plants. Am. J. Bot. 100, 1022-37 (2013).</td>
  </tr>
  <tr>
    <td>1238</td>
    <td>86. Porebski, S. L.G. Bailey, & Baum, B.R. Modification of a CTAB DNA extraction protocol for plants containing high polysaccharide and polyphenol components. Plant Mol. Biol. Rep. 15, 8-15 (1997).</td>
  </tr>
  <tr>
    <td>1239</td>
    <td>87. Marçais, G. and C. Kingsford, A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics 27, 764-70 (2011).</td>
  </tr>
  <tr>
    <td>1240</td>
    <td>88. Liu, B. et al. Estimation of genomic characteristics by analyzing k-mer frequency in de novo genome projects. arXiv arXiv:1308.2012 (2013).</td>
  </tr>
  <tr>
    <td>1241</td>
    <td>89. Zhang, X. et al. Assembly of allele-aware, chromosomal-scale autopolyploid genomes based on Hi-C data. Nat. Plants 5, 833-845 (2019).</td>
  </tr>
  <tr>
    <td>1242</td>
    <td>90. Dudchenko, O. et al. De novo assembly of the Aedes aegypti genome using Hi-C yields chromosome-length scaffolds. Science 356, 92-95 (2017).</td>
  </tr>
  <tr>
    <td>1243</td>
    <td>91. Wolff, J. et al. Galaxy HiCExplorer 3: a web server for reproducible Hi-C, capture Hi-C and single-cell Hi-C data analysis, quality control and visualization. Nucleic Acids Res. 48, W177-W184 (2020).</td>
  </tr>
  <tr>
    <td>1244</td>
    <td>92. Jain, C. et al. Weighted minimizer sampling improves long read mapping. Bioinformatics 36, i111-i118 (2020).</td>
  </tr>
  <tr>
    <td>1245</td>
    <td>93. Kurtz, S. et al. Versatile and open software for comparing large genomes. Genome Biol. 5, R12 (2004).</td>
  </tr>
  <tr>
    <td>1246</td>
    <td>94. Manni, M. et al. BUSCO: Assessing Genomic Data Quality and Beyond. Curr. Protoc. 1, e323 (2021).</td>
  </tr>
  <tr>
    <td>1247</td>
    <td>95. Rhie, A. et al. Merqury: reference-free quality, completeness, and phasing assessment for genome assemblies. Genome Biol. 21, 245 (2020).</td>
  </tr>
  <tr>
    <td>1248</td>
    <td>96. Marçais, G. et al. MUMmer4: A fast and versatile genome alignment system. PLoS Comput. Biol. 14, e1005944 (2018).</td>
  </tr>
  <tr>
    <td>1249</td>
    <td>97. Xu, Z. & Wang, H. LTR_FINDER: an efficient tool for the prediction of full-length LTR retrotransposons. Nucleic Acids Res. 35, W265-8 (2007).</td>
  </tr>
  <tr>
    <td>1250</td>
    <td>98. Ellinghaus, D. S. Kurtz, & Willhoeft U. LTRharvest, an efficient and flexible software for de novo detection of LTR retrotransposons. BMC Bioinformatics 9, 18 (2008).</td>
  </tr>
  <tr>
    <td>1251</td>
    <td>99. Ou, S. & Jiang N. LTR_retriever: A Highly Accurate and Sensitive Program for Identification of Long Terminal Repeat Retrotransposons. Plant Physiol. 176, 1410-1422 (2018).</td>
  </tr>
  <tr>
    <td>1252</td>
    <td>100. Abrusán, G. et al. TEclass--a tool for automated classification of unknown eukaryotic transposable elements. Bioinformatics 25, 1329-30 (2009).</td>
  </tr>
  <tr>
    <td>1253</td>
    <td>101. Benson, G. Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Res. 27, 573-80 (1999).</td>
  </tr>
  <tr>
    <td>1254</td>
    <td>102. Beier, S. et al. MISA-web: a web server for microsatellite prediction. Bioinformatics 33, 2583-2585 (2017).</td>
  </tr>
  <tr>
    <td>1255</td>
    <td>103. Kovaka, S. et al. Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol. 20, 278 (2019).</td>
  </tr>
  <tr>
    <td>1256</td>
    <td>104. Siadjeu, C. E. Mayland-Quellhorst, & Albach D.C. Genetic diversity and population structure of trifoliate yam (Dioscorea dumetorum Kunth) in Cameroon revealed by genotyping-by-sequencing (GBS). BMC Plant Biol. 18, 359 (2018).</td>
  </tr>
  <tr>
    <td>1257</td>
    <td>105. Li, Y. et al. The genome of Dioscorea zingiberensis sheds light on the biosynthesis, origin and evolution of the medicinally important diosgenin saponins. Hortic. Res. 9, uhac165 (2022).</td>
  </tr>
  <tr>
    <td>1258</td>
    <td>106. Vadakkemukadiyil Chellappan, B. et al. High Quality Draft Genome of Arogyapacha (Trichopus zeylanicus), an Important Medicinal Plant Endemic to Western Ghats of India. G3 9, 2395-2404 (2019).</td>
  </tr>
  <tr>
    <td>1259</td>
    <td>107. Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinformatics 10, 421 (2009).</td>
  </tr>
  <tr>
    <td>1260</td>
    <td>108. Slater, G.S. & Birney E. Automated generation of heuristics for biological sequence comparison. BMC Bioinformatics 6, 31 (2005).</td>
  </tr>
  <tr>
    <td>1261</td>
    <td>109. Stanke, M. et al. Using native and syntenically mapped cDNA alignments to improve de novo gene finding. Bioinformatics 24, 637-44 (2008).</td>
  </tr>
  <tr>
    <td>1262</td>
    <td>110. Delcher, A.L. et al. Identifying bacterial genes and endosymbiont DNA with Glimmer. Bioinformatics 23, 673-9 (2007).</td>
  </tr>
  <tr>
    <td>1263</td>
    <td>111. Holt, C. & Yandell, M. MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects. BMC Bioinformatics 12, 491 (2011).</td>
  </tr>
  <tr>
    <td>1264</td>
    <td>112. Chan, et al. tRNAscan-SE 2.0: improved detection and functional classification of transfer RNA genes. Nucleic Acids Res. 49, 9077-9096 (2021).</td>
  </tr>
  <tr>
    <td>1265</td>
    <td>113. Lagesen, K. et al. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. Nucleic Acids Res. 35, 3100-8 (2007).</td>
  </tr>
  <tr>
    <td>1266</td>
    <td>114. Nawrocki, E. and S.R. Eddy, Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics 29, 2933-5 (2013).</td>
  </tr>
  <tr>
    <td>1267</td>
    <td>115. Kalvari, I. et al. Rfam 14: expanded coverage of metagenomic, viral and microRNA families. Nucleic</td>
  </tr>
</table>
1296 Acids Res. 49, D192-D200 (2021).
1297 116. Deng, Y. et al. A telomere-to-telomere gap-free reference genome of watermelon and its mutation library provide important resources for gene discovery and breeding. Mol. Plant 15, 1268-1284 (2022).
1298 117. Song, J.M. et al. Two gap-free reference genomes and a global view of the centromere architecture in rice. Mol. Plant 14, 1757-1767 (2021).
1299 118. Quinlan, A.R. & Hall I.M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841-2 (2010).
1300 119. Han, X. et al. Two haplotype-resolved, gap-free genome assemblies for Actinidia latifolia and Actinidia chinensis shed light on the regulatory mechanisms of vitamin C and sucrose metabolism in kiwifruit. Mol. Plant 16, 452-470 (2023).
1301 120. Li, F. et al. Gap-free genome assembly and comparative analysis reveal the evolution and anthocyanin accumulation mechanism of Rhodomyrtus tomentosa. Hortic. Res. 10, uhad005 (2023).
1302 121. Hao, Z. et al. RIdeogram: drawing SVG graphics to visualize and map genome-wide data on the idiograms. PeerJ Comput. Sci. 6, e251 (2020).
1303 122. Frith, M.C., Hamada, M. & Horton, P. Parameters for accurate genome alignment. BMC Bioinformatics 11, 80 (2010).
1304 123. Tang, H. et al. JCVI: A versatile toolkit for comparative genomics analysis. Imeta 3, e211 (2024).
1305 124. Emms, D.M. & Kelly, S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 20, 238 (2019).
1306 125. Edgar, R.C. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792-7 (2004).
1307 126. Capella-Gutiérrez, S., Silla-Martínez, J.M. & Gabaldón T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25, 1972-3 (2009).
1308 127. Stamatakis, A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 20, 1312-3 (2014).
1309 128. Yang, Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586-91 (2007).
1310 129. Chen, S. Ultrafast one-pass FASTQ data preprocessing, quality control, and deduplication using fast Imeta 2, e107 (2023).
1311 130. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15-21 (2013).
1312 131. Li, B. & Dewey C.N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323 (2011).
1313 132. Robinson, M.D., McCarthy, D.J. & Smyth, G.K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139-40 (2010).
1314 133. Wu, T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2, 100141 (2021).
1315 134. Ashburner, M. et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 25, 25-9 (2000).
1316 135. Kanehisa, M. & Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27-30 (2000).
1317 136. McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297-303 (2010).
1318 137. Danecek, et al. The variant call format and VCFtools. Bioinformatics 27, 2156-8 (2011).
1319 138. Gu, Z. et al. circlize Implements and enhances circular visualization in R. Bioinformatics 30, 2811-2 (2014).
1320 139. Purcell, S. et al. PLINK: a toolset for whole-genome association and population-based linkage analysis. American Journal of Human Genetics 81 (2007).
1321 140. Kokot M., Długosz, M., Deorowicz, S. KMC 3: counting and manipulating k-mer statistics. Bioinformatics 33, 2759-2761 (2017).
1322 141. Danecek, P. et al. Twelve years of SAMtools and BCFtools. GigaScience 10, giab008 (2021).
1323 142. Pedersen, B. S. & Quinlan, A. R. Mosdepth: quick coverage calculation for genomes and exomes. Bioinformatics 34, 867-868 (2018).
1324 143. Rozewicki, J. et al. MAFFT-DASH: integrated protein sequence and structural alignment. Nucleic Acids Research 47, W5-W10 (2019).
1325 144. Suyama, M., Torrents, D. & Bork P. PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments. Nucleic Acids Res. 1, W609-12 (2006).
1326 145. Wang, D. et al. KaKs_Calculator 2.0: a toolkit incorporating gamma-series methods and sliding window strategies. Genomics Proteomics Bioinformatics 8, 77-80 (2010).
1327 146. von Maydell, D. PCR Allele Competitive Extension (PACE). Methods Mol. Biol. 2638, 263-271 (2023).
1328 147. Ni, P. et al. DNA 5-methylcytosine detection and methylation phasing using PacBio circular consensus sequencing. Nat. Commun. 14, 4054 (2023).
1329 138. Gu, Z. et al. circlize Implements and enhances circular visualization in R. Bioinformatics 30, 2811-2 (2014).
1330 139. Purcell, S. et al. PLINK: a toolset for whole-genome association and population-based linkage analysis. American Journal of Human Genetics 81 (2007).
1331 140. Kokot M., Długosz, M., Deorowicz, S. KMC 3: counting and manipulating k-mer statistics. Bioinformatics 33, 2759-2761 (2017).
1332 141. Danecek, P. et al. Twelve years of SAMtools and BCFtools. GigaScience 10, giab008 (2021).
1333 142. Pedersen, B. S. & Quinlan, A. R. Mosdepth: quick coverage calculation for genomes and exomes. Bioinformatics 34, 867-868 (2018).
1334 143. Rozewicki, J. et al. MAFFT-DASH: integrated protein sequence and structural alignment. Nucleic Acids Research 47, W5-W10 (2019).
1335 144. Suyama, M., Torrents, D. & Bork P. PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments. Nucleic Acids Res. 1, W609-12 (2006).
1336 145. Wang, D. et al. KaKs_Calculator 2.0: a toolkit incorporating gamma-series methods and sliding window strategies. Genomics Proteomics Bioinformatics 8, 77-80 (2010).
1337 146. von Maydell, D. PCR Allele Competitive Extension (PACE). Methods Mol. Biol. 2638, 263-271 (2023).
1338 147. Ni, P. et al. DNA 5-methylcytosine detection and methylation phasing using PacBio circular consensus sequencing. Nat. Commun. 14, 4054 (2023).
1356 148. Nie, F. et al. De novo diploid genome assembly using long noisy reads. *Nat. Commun.* **15**, 2964 (2024).

1357

1358
## Figures

a

<figure><img src="image_2.png" /></figure>

b

<figure><img src="image_3.png" /></figure>

c

<figure><img src="image_4.png" /></figure>

d

<figure><img src="image_5.png" /><figcaption>**Figure 1**</figcaption></figure>

**Chromosome-level haplotype-resolved near-complete assembly of the male *D. alata* cultivar 'Kabusa'.** (a) Pictures of 'Kabusa' cultivar displaying the aerial biomass in the field, the tuber shape and flesh color, and the male flowers. (b) Genome-wide Hi-C contact heatmap, visualized with Juicebox³². (c) Circle diagram of basic genome information. Concentric circles from outside to inside represent chromosomes, gene density, repeat density, and GC density calculated in a 50 Kb window. Lines in the interior indicate collinearity. The scale bar = 10 cm. (d) Syntenic mapping between the two haplotypes of the 'Kabusa' genome using SyRI³³. Collinear regions between HaplotypeA and HaplotypeB are indicated with gray lines, while inversion, translocation, and duplication regions are indicated with orange, green, and blue lines, respectively.
<figure><img src="image_6.png" /><figcaption>Figure 2</figcaption></figure>

**Structure of the sex determination region (SDR) in *D. alata*.** (a) Hi-C contact heatmap between DalaChr6A (Y) and DalaChr6B (X), visualized with Juicebox [31]. (b) Syntenic mapping between 'Kabusa' Y chromosome and X chromosomes from 'Kabusa', Dala_V2<sup>27</sup> and Ziyushenshu<sup>28</sup> using SyRI<sup>33</sup>. Four inversions (INV1-INV4) were detected. (c) The genomic landscape of genetic variants between male and female pools. The outer track represents the 20 chromosomes (chr) of HaplotypeA in the 'Kabusa' genome. The inner tracks from outside to inside represent *Fst* between females and males and male-specific SNP density, respectively (window size 100 Kb). (d) Manhattan plot for SNP-sex phenotype genome-wide association study (GWAS), with the peaks indicating significant GWAS signals, and the black and red dotted horizontal lines indicating the suggestive and genome-wide significance thresholds, respectively. (e) Manhattan plot for InDel-sex phenotype GWAS, with the peaks indicating significant GWAS signals, and the black and red dotted horizontal lines indicating the suggestive and genome-wide significance thresholds, respectively. Regional Manhattan plots on DalaChr6A showing the dispersion of the significant SNPs (f) and Indels (g). (h) Differential read coverage (10 Kb window) between male and female groups along the Y chromosome. (i) Exact matches of male specific 27-mers along the Y chromosome. (j) Representation of the X (DalaChr6B) and Y (DalaChr6A) chromosomes along with the gene density (heatmap), centromere and telomere locations, pseudoautosomal regions (PAR), and the SDR containing the 4 inversions (INV1-INV4) detected.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 3</figcaption></figure>

**Sex chromosome divergence in *D. alata*.** (a) Comparison of the synonymous substitution rates (*Ks*) of genes between SDR and PAR on the Y chromosome. (b) *Ks* comparison between Y chromosome components. \*\* = significant difference at *p*<0.001. (c) Comparison of the LTR/Copia and LTR/Gypsy densities between PAR and SDR in the Y chromosome. (d) Comparison of the LTR/Copia and LTR/Gypsy densities between PAR and HXR in the X chromosome. (e) Comparison of the LTR/Copia and LTR/Gypsy densities between SDR and HXR. (f) Smoothed density plots along X and Y chromosomes. Densities are calculated per 1 Mb overlapping windows.
<figure><img src="image_13.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>i</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>j</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>k</figcaption></figure>

<figure><img src="image_23.png" /></figure>

Figure 4

**Sex-biased gene expression in *D. alata*.** (a) Development of *D. alata* flower at Day5, Day8 and Day12 after blooming in male 'Kabusa' and female '74F'. Ontogenetic differences between female and male flowers become morphologically evident between Day8 and Day12 after blooming. Scale = 1 cm. (b) Principal component analysis based on the gene expression levels for flower samples (Day5) expressed as fragments per kilobase of transcript per million fragments mapped in the male and female groups. The first and the second principal components explaining 55% and 20% of the total variation, respectively. (c) Volcano plot displaying the differential expression analysis in flowers between male and female groups. Dots in green, red, and gray represent genes downregulated, upregulated and not differentially expressed in male flowers compared to female flowers, respectively. (d) Gene expression patterns in flower buds (log₂fold change between male and female groups) along the sex chromosomes. Each dot represents a sex-biased gene (SBG). Statistical comparison of SBG located within the sex chromosomes (e) and within the Y chromosome components (f). (g) Statistical comparison of SGB in mature tuber samples located within the sex chromosomes. Samples represent mature tuber from three female *D. alata* accessions ('Dou', 'Kinabayo', and '74F') and three male accessions ('Divin', 'Plimbite', and 'Toufi-Tetea'). * = significant difference at *p* <0.05. (h) Comparison of the methylation levels (5-methylcytosine at CpG
sites) across sex chromosomes and statistical comparison between sex chromosomes components: SDR vs PAR (i), HXR vs PAR (j) and SDR vs HXR (k). Methylation levels are calculated per 100 Kb overlapping windows.

## Supplementary Files

This is a list of supplementary files associated with this preprint. Click to download.

- SupplementaryTables.xlsx

- SupplementaryFigures.pdf