Medicago super-pangenome reveals adaptive advantages and evolutionary constraints in autotetraploid alfalfa

Domain: biology
Published: 2025-12-15
Editorial outcome: major_revision

Abstract
The genetic basis for the adaptive advantages of polyploids over their diploid relatives remains poorly understood. To address this knowledge gap, we generate a haplotype-resolved autotetraploid alfalfa (Medicago sativa subsp. sativa) genome and construct a super-pangenome from 13 genomes across seven Medicago taxa. We discover substantial gene content variation in alfalfa, with only 20.1% of genes present on all four haplotypes. Within this group, 53.3% are core genes conserved across the Medicago genus, which we term ‘tetra-copy core genes’. We find these genes are significantly enriched in climate-adaptation-associated genes (1.60-fold) and stress-responsive differentially expressed genes (1.61-fold). Paradoxically, they also carry a high genetic burden, with 80.1% of deleterious variants located in coding regions. Indeed, overexpressing a representative tetra-copy core gene, the glycine decarboxylase (MsGDC), improves both biomass and nitrogen use efficiency, despite its high genetic burden. Our study reveals the trade-off between adaptation and evolutionary constraints mediated by tetra-copy core genes, facilitating polyploid genetics and alfalfa breeding.

Introduction
Polyploidy, or whole genome duplication (WGD), has been considered an important evolutionary driving force in plant evolution, impacting from the past to the present and influencing levels from the molecular scale to the ecological level1,2. Polyploidy is estimated to have occurred in 30% to 70% of angiosperms, either as currently existing polyploids or as ancient polyploidies (paleopolyploidies)3–5. This process encompasses gene duplication, gene loss, genetic exchange among sub-genomes, and involves molecular and physiological regulation, as well as gene expression and epigenetic regulation4,6,7. WGD provides genetic resources for evolutionary innovation and adaptation, enhancing fitness and phenotypic diversity. Compared to their diploid progenitors, polyploids tend to have larger organism size, enhanced competitive advantages, and higher fitness levels. Consequently, they often exhibit greater adaptability to environmental changes and occupy broader ecological niches3,4,8–10.

Polyploids exhibit four main potential adaptive strategies. (1) In the short term, polyploids can mask deleterious mutations because recessive harmful variants are effectively hidden in a polyploid population, reducing genetic burden. This provides newly formed polyploids with immediate advantages, allowing them to establish niches. However, these advantages may diminish over time due to the accumulation of mutations3,5. (2) Adaptive processes may become more efficient after WGD due to an increase in effective population size and a higher rate of beneficial mutations, leading to enhanced selective efficiency and fitness4,5,11,12. (3) Duplicated genes offer genetic material for the evolution of new functions (neofunctionalization) or more specific functions (sub-functionalization). This process can lead to significant biological innovations, facilitating adaptation to new niches13–16. (4) Changes in gene expression and epigenetic patterns can alter physiology, metabolism, and morphology, thereby enhancing adaptive abilities4,7,9,15. While polyploidy is associated with an increase in adaptability, it is also regarded as an evolutionary ‘dead-end’. Several theories attempt to explain this fate, including the accumulation of deleterious variants4,5,17, lower diversification rates and higher extinction rates compared to closely related diploids18,19, highly dynamic genomes3,12,20 and abnormal meiosis with multivalent pairing among three or more homologous chromosomes6,13. Yet, the evidence for these hypotheses is largely derived from paleo-polyploids, forcing inferences to be drawn from the complex genomic signatures of ancient events. Furthermore, current research on polyploids predominantly prioritizes allopolyploids, which involve both WGD and hybridization effects. This combination obscures the individual contributions of polyploidy and hybridization, complicating our understanding of how polyploidy alone affects adaptability. Therefore, studying autopolyploids, which result from WGD within a single species, eliminates the confounding effects of hybridization, making them more suitable for investigating how evolutionary adaptability is influenced by WGD3,13.

The construction of a pangenome provides a powerful framework to analyze patterns of gene retention and loss, shedding light on preferences for gene conservation that relate to adaptability. For instance, core gene families were involved in fundamental biological functions in both cicer and arabidopsis21,22. Private genes play a crucial role in adaptation to local environment and climate conditions in populus23. Dispensable and private genes were enriched for abiotic and biotic-related genes in soybean24. Moreover, previous studies have demonstrated that gene retention following WGD is a nonrandom process4,9. The unique patterns of duplicated gene retention are correlated with fitness and adaptability4,9,25. The duplication of genes can create redundancy, which serves as a foundation for genetic and phenotypic innovation and can be utilized for adapting to environmental changes9,25. However, the identification of duplicated genes in polyploids has largely relied on gene family tree1,26 or synonymous substitution rates (Ks values)27,28, which are influenced by divergence time and nucleotide substitution rates. With advancements in sequencing and assembly strategies, the identification of allelic copy number variation, which results from gene duplication and retention following WGD, in autopolyploids can now be achieved with greater accuracy using phased haplotype-resolved genomes29–31. Therefore, a comprehensive approach that integrates super-pangenome analyses with phased autopolyploid genomes, uncovering gene retention over long evolutionary periods and allelic imbalance, provides a powerful framework for an in-depth understanding of adaptive mechanisms in polyploids32,33.

The genus Medicago, a member of the legume family, is highly diverse. It contains 87 species with various ploidy levels (diploid, tetraploid, and hexaploid) that are adapted to a wide range of geographic distributions and climates. This genus is famous for its high protein production, nitrogen-fixing ability, and positive impact on soil ecology. Its members include diploid model plant M. truncatula and the autotetraploid high-quality forage alfalfa (M. sativa subsp. sativa, hereafter subsp. sativa). Compared to the diploid ancestor (M. sativa subsp. caerulea, hereafter subsp. caerulea or diploid alfalfa), the cultivated autotetraploid alfalfa presents higher yield, greater stress tolerance, improved climate adaptability, and a wider distribution range34. Understanding the genomic basis of these adaptive advantages has recently become feasible with the release of several Medicago genomes, including M. truncatula35–37, M. polymorpha38, M. ruthenica39,40, subsp. caerulea41 and alfalfa42–44. However, the absence of a high-quality phased tetraploid alfalfa genome makes it challenging to identify allelic copy number variation following WGD, thereby hindering our understanding of autopolyploid adaptability.

In this study, we generate a high-quality, haplotype-resolved autotetraploid alfalfa genome alongside a comprehensive super-pangenome for the Medicago genus using 13 genomes from seven taxa (six species, including two subspecies of M. sativa). Our analysis of gene retention reveals a distinct set of genes that we term ‘tetra-copy core genes’, which are defined as genes retained on all four allelic copies in alfalfa that are also core members of the Medicago genus. We demonstrate that these tetra-copy core genes are significantly enriched in climate-adaptation-associated genes and stress-related differentially expressed genes. Moreover, tetra-copy core genes exhibit higher overall expression than other genes, with one allelic copy showing significantly higher expression than the other three. Paradoxically, these functionally important genes also carry a disproportionately high genetic burden from deleterious mutations. We functionally validate this trade-off by showing that overexpression of a representative tetra-copy core gene, MsGDC, improves both biomass and nitrogen use efficiency, despite its high genetic load. Our study resolves a central paradox of polyploid adaptation by demonstrating how ‘tetra-copy core genes’ mediate a crucial balance between adaptive advantages and deleterious mutation accumulation. These findings provide a framework for understanding autopolyploid evolution and accelerating the genetic improvement of alfalfa.

Results
Genome assemblies and annotations of Medicago genomes
We sequenced two M. sativa subspecies (autotetraploid subsp. sativa and diploid subsp. caerulea) employing a combination of long reads sequencing platforms (PacBio and Nanopore sequencing) and high-throughput chromosome conformation capture sequencing (Hi-C) approaches. For autotetraploid alfalfa, we selected the high-yielding and salt-tolerant Chinese alfalfa cultivar Zhongmu-4 as the representative material. The Zhongmu-4 genome was sequenced, generating 109.84 Gb (~137×, compared to the ~800 Mb haploid genome of alfalfa) of high-fidelity (HiFi) long reads and 18.25 Gb (~23×) of Oxford Nanopore Technology (ONT) long reads. Combined with our previously sequenced 285 Gb (~356×) Hi-C data44, we conducted the genome assembly (Supplementary Fig. 1). We initially generated contigs using HiFi and ONT reads with Hifiasm45, and phased them into 32 chromosomes using Hi-C data, resulting in an assembly size of 2.54 Gb (ZM4_V1.5, Supplementary Fig. 2). To enhance completeness and contiguity, we then divided the HiFi reads into 32 groups based on the haplotype information of our previously published ZM4_V1.0 and the ZM4_V1.5 genomes44 (Supplementary Fig. 3). The genome was assembled into contigs using haplotype-resolved HiFi reads with Hifiasm, and subsequently improved to scaffold level using Ragtag (see Methods, Supplementary Fig. 1).

Finally, we assembled ZM4_V2.0, a high-quality phased genome with a total size of 3.13 Gb across 32 chromosomes (contig N50=5.45 Mb). This assembly size is consistent with previous estimations for cultivated alfalfa and represents a significant improvement in quality over our previous version, ZM4_V1.0 (Supplementary Fig. 4)44. The raw HiFi and ONT data were mapped to ZM4_V2.0, achieving overall mapping rates of 99.93% and 99.92%, respectively.

To establish a high-confidence set of alignments for assessing haplotype fidelity, we applied a MAPQ≥10 filter, which retained 83.02% of HiFi reads and 88.34% of ONT reads. Critically, over 99% of the mapped HiFi reads aligned correctly to a single corresponding haplotype (Supplementary Fig. 5). Haplotype phasing analysis revealed a switch error rate of 1.97% and a Hamming error rate of 7.12%. This provides strong evidence for the high accuracy of our phasing framework, demonstrating that most haplotypic regions were successfully resolved. The Benchmarking Universal Single-Copy Orthologs (BUSCO) analysis of the genome assembly, using the embryophyta_odb10 database (1614 orthologs), indicated a completeness score of 99.5%. Finally, using Merqury software46, we assessed the base quality and genome completeness, resulting in scores of 61.88 and 97.8%, respectively (Supplementary Table 1). Therefore, the genome quality of ZM4_V2.0 is sufficiently high to serve as a reference genome for alfalfa and to support WGD analysis in autopolyploid. The sizes of four haplotype (hap) genomes were 791.79 Mb (contig N50=6.52 Mb) for hap1, 751.53 Mb (contig N50=4.98 Mb) for hap2, 809.40 Mb (contig N50=5.36 Mb) for hap3 and 775.96 Mb (contig N50=5.45 Mb) for hap4. The four haplotypes showed significant gene collinearity (Supplementary Fig. 6) and strong genomic synteny with certain structural variations (SVs) (Supplementary Fig. 7). Additionally, a high collinearity was observed between the M. truncatula A17 (V5.0) and the ZM4_V2.0 genomes (Supplementary Figs. 8 and 9).

Using a homology-based, short and long-read transcriptomic annotation pipeline, we predicted 202,473 genes, and the number of genes was 51,000, 49,392, 50,958, and 51,123 for hap1, hap2, hap3, and hap4, respectively. The annotation was highly complete, with a BUSCO score of 99.8%, and 187,026 genes (92.4%) were functionally annotated according to the EggNOG database47. Taken together, ZM4_V2.0 represents a substantial improvement over the published alfalfa reference genomes42–44.

Besides autotetraploid alfalfa de novo assembly, we also assembled a diploid subsp. caerulea genome employing 207.00 Gb PacBio CLR long reads, 52.75 Gb Illumina short reads, and 204.28 Gb Hi-C data. The assembled genome was 801.63 Mb (contig N50=1.61 Mb), with 96.4% BUSCO completeness, and the annotated gene number was 50,491. Furthermore, we annotated 47,963 genes for M. lupulina and 46,150 genes for M. arabica using the assembled genome provided by the Darwin Tree of Life Project (https://www.darwintreeoflife.org/)48,49.

Genome evolution and gene conservation revealed by Medicago super-pangenome
Eleven published Medicago genomes were downloaded, and genome annotation was redone. Firstly, the total count, identity, and orientation of chromosomes were revised based on the model plant genome of M. truncatula A17 (V5.0). Specifically, we revised the chromosome number of M. polymorpha from 7 to 8 (Supplementary Fig. 10). Secondly, gene annotation was performed using the same pipeline as that used for ZM4_V2.0. The number of annotated genes and the BUSCO assessment were improved compared to the raw version (Table 1).

Using 13 Medicago genomes (19 haplotypes), we constructed a graph-based super-pangenome (Alfalfapan-V1.0) reference of seven Medicago taxa (Fig. 1A). The Alfalfapan-V1.0 pangenome was constructed using the ZM4_V2.0 hap1 genome as a reference. It encompasses a size of 7.70 Gb across 8 chromosomes, with an average chromosome size of 968.82 Mb. The pangenome is 2.46-fold the size of the alfalfa ZM4_V2.0 genome, 9.72-fold that of the ZM4_V2.0 hap1 genome and 17.91-fold that of the M. truncatula A17 genome (Fig. 1D). We found that genome size in Medicago was highly correlated with transposable elements (TEs) content, indicating that TEs expansion is a primary driver of genome size evolution and divergence. Similarly, a comparable pattern of TEs contribution to genome expansion was observed in Oryza50. The two most abundant types of TEs, LTR Gypsy and LTR Copia, together contributed to ~50% of the total TEs content, underscoring their substantial role in the divergence of Medicago genomes (Fig. 1E).

To investigate the evolution of gene content across the Medicago genus, we constructed a pangenome by identifying orthologous gene groups with OrthoFinder51. The gene based synteny analysis revealed high collinearity among Medicago species (Fig. 1A). Furthermore, we examined a known translocation involving chromosomes 4 and 8 that distinguishes M. truncatula A17 from alfalfa41,43,44. Our multi-genome comparison reveals this translocation to be an accession-specific feature of the M. truncatula A17 genome, rather than a conserved, species-level difference. This highlights how a pangenome approach moves beyond simple pairwise comparison to provide a more accurate and comprehensive view of SVs. Based on 5,106 single-copy orthologous genes of seven taxa, we constructed a phylogenetic tree of Medicago (Fig. 1B). Our primary phylogenetic analysis resolved the overall species relationships among the studied Medicago taxa. In this tree, M. truncatula was identified as the sister group to the clade containing alfalfa and subsp. caerulea, while M. ruthenica was resolved as the earliest-diverging lineage. Interestingly, a more detailed analysis focusing specifically on the alfalfa haplotypes revealed a polyphyletic pattern (Supplementary Fig. 11). This topology suggests multiple independent origins of tetraploid alfalfa from genetically divergent diploid subsp. caerulea ancestors. Divergence time estimation placed the crown age of the studied Medicago clade at ~7.59 million years ago (Mya), with the tetraploid alfalfa lineage originating from its diploid progenitor ~2.13 Mya (Fig. 1B).

A total of 69,767 orthologous gene families were identified from Alfalfapan-V1.0. Based on these gene families, the Alfalfapan-V1.0 super-pangenome is composed of 13,196 core (present in all genomes), 8,454 softcore (present in 11 and 12 genomes), 33,767 dispensable (present in two to ten genomes), and 14,350 private (only present in one genome) gene families (Fig. 2A, Supplementary Table 2). The total number of pan-genes increased as more genomes were added, while the number of core genes steadily decreased (Fig. 2B). Specifically, ZM4_V2.0 contained the highest number of orthologous gene families, totaling 48,151 (69.0% of pangenome). For the other genomes, the count of gene families ranged from 26,997 in Mpoly to 42,646 in XJDY (Fig. 2C, Supplementary Table 3). It suggested that Alfalfapan-V1.0 contains at least 1.45 times more gene families (69,767 vs 48,151) than a single alfalfa genome. Interestingly, core gene patterns tend to cluster at the ends of chromosomes, while TEs are predominantly located in the pericentromeric regions. This suggests that core genes are more stable and resistant to TEs interference (Fig. 2D).

The SVs were called for each genome using the ZM4_V2.0 hap1 as a reference. The results indicated that the M. sativa-falcata complex (including subsp. sativa and subsp. caerulea) exhibited a significantly higher number of SVs compared to other Medicago species (Fig. 2E). This observation is consistent with the significantly lower number of mapped reads observed for the more distantly related species (Supplementary Fig. 12). Consequently, the true number of SVs in these related species was likely underestimated, highlighting the challenges of reliable SV calling at the genus level. Large SVs between species and haplotypes were often found in corresponding chromosomal regions (Supplementary Figs. 13 and 14). Additionally, we generated population-level SVs using the PGGB and vg pipelines52,53. A total of 494,797 SVs were identified among the 13 genomes, with 208,045 SVs found exclusively in the M. sativa-falcata complex (Fig. 2F). Using this comprehensive set of SVs, we further investigated the genetic relationships among these 13 genomes. The PCA results demonstrated clear separation among different species (Fig. 2G). Interestingly, we detected divergence among the four alfalfa haplotypes across both the ZM4_V2.0 and XJDY (Fig. 2H). However, the specific patterns of this divergence varied greatly at the individual chromosome level (Supplementary Fig. 15). Together, these findings demonstrate that SVs are powerful markers for resolving the complex and variable evolutionary histories of different sub-genomic regions in alfalfa.

The retention of core genes correlated with climate adaptation of alfalfa
Compared with the diploid ancestor subsp. caerulea, autotetraploid alfalfa has expanded into a wider range of climate regions and occupied broader ecological niches (Fig. 1C). These results suggested that autopolyploidy plays a major role in the local adaptation of alfalfa. To confirm that the allelic copy number variation originated from WGD without divergence, we carefully identified genes that displayed a 95.0% or higher protein sequence similarity across the four haplotypes in ZM4_V2.0 as allelic copies (Fig. 3A, see Methods). The results indicated that a total of 109,715 (109,715/202,473, 54.2%) genes in ZM4_V2.0 could be considered non-redundant genes (representing each gene only once despite allelic multiplicity). Among these, 49,344 (49,344/109,715, 45.0%) genes had one allelic copy, 23,095 (21.0%) had two allelic copies, 15,264 (13.9%) had three allelic copies, and 22,012 (20.1%) had four allelic copies (Supplementary Table 4). One representative genomic region on Chr1 was selected to illustrate allelic copy number variation (Fig. 3B). Comparing the relationship between gene conservation and allelic copy number variance, we found that genes with four allelic copies constituted a substantial portion of the core genes (proportion: 53.3%; number: 11,735; chi-squared test p < 2.2e-16), whereas genes with one allelic copy represented a substantial portion of the dispensable genes (proportion: 41.7%; number: 20,562; chi-squared test p < 2.2e-16) (Fig. 3C, Supplementary Table 5). It indicated that gene retention and loss following WGD were nonrandom. Analysis of selective pressure indicated that the log(Ka/Ks) value gradually decreased as the allelic copy number increased. Core genes with four allelic copies exhibited the lowest median Ka/Ks value (Fig. 3D).

The core genes that retained all four allelic copies in alfalfa represent a set of evolutionarily conserved duplicates within the Medicago genus; we therefore term them ‘tetra-copy core genes’. In contrast, dispensable genes present as one allelic copy in alfalfa exhibited a high median Ka/Ks ratio, suggesting they are undergoing rapid adaptive evolution, possibly as a mechanism for adaptation to new environmental pressures; we therefore assigned them the name ‘unique-essential genes’. We developed a genome-wide gene classification framework by partitioning all 202,473 genes in ZM4_V2.0 into six groups based on gene conservation and allelic copy number variation. These groups include unique-essential genes (proportion: 10.2%; number: 20,562), one-copy genes (excluding unique-essential, 14.2%, 28,782), two-copy genes (18.4%, 37,344), three-copy genes (17.9%, 36,330), four-copy genes (excluding tetra-copy core genes, 17.7%, 35,873), and tetra-copy core genes (21.5%, 43,582) (Fig. 3E, Supplementary Table 6). These six partitions were established as reference gene groups for various downstream analyses. For instance, they provide a genomic baseline proportion for assessing the enrichment of gene sets such as differentially expressed genes (DEGs) and candidates from genome scans for climate adaptation.

To examine the preferential retention and loss of genes related to adaptation following autopolyploidy, we further analyzed the contribution of each type of genes to local climate adaptation. At first, we examined the proportion of climate adaptation related genes within each of these six groups. The climate adaptation-associated genes were calculated by latent factor mixed model (LFMM) and population branch statistic (PBS) analysis54. A total of 2,429 genes were identified to be associated with climate adaptation. These genes were not evenly distributed among the six gene groups; instead, they were significantly enriched in the tetra-copy core genes. This single group accounted for 836 (34.4%) of the adaptation-associated genes, representing a 1.60-fold enrichment over the genomic baseline proportion of 21.5% (chi-squared test, p < 2.2e-16; Figs. 3E and 4A). Our finding suggested that tetra-copy core genes play a significant role in climate adaptation.

To investigate the role of tetra-copy core genes in environmental adaptation, we analyzed RNA-seq datasets related to biotic/abiotic stress from 23 publications/research groups (covering 124 distinct expression conditions) (Supplementary Data 1 and 2). Our analysis focused on identifying both DEGs and genes ranked within the top percentiles of expression in each condition. Although the total number of DEGs varied significantly across treatments (mean = 689), the proportional contribution of six gene groups to the DEG sets remained remarkably consistent. The group of tetra-copy core genes was significantly overrepresented among the DEGs. The median proportion of these genes was 34.6%, a 1.61-fold enrichment over the genomic baseline proportion of 21.5% (one-sample Wilcoxon test, p = 1.93e-4; Fig. 4B). Furthermore, unique-essential genes and one-copy genes had the highest fold change of DEGs (Supplementary Fig. 16). These results suggest that tetra-copy core genes contribute significantly to stress tolerance, while unique-essential genes are very susceptible to stress.

The level of gene expression is a key determinant of adaptive responses to environmental stress55. After analyzing six gene groups, we observed that genes with more allelic copies tended to have higher expression levels, and similar expression patterns were observed in autotetraploid potato56. Additionally, we found that the tetra-copy core genes had the highest median expression level (Fig. 4C). Analysis of the four allelic copies within each non-redundant gene revealed a characteristic stepwise gradient, in which one allele consistently exhibited significantly higher expression than the other three (Fig. 4D). This suggests that the high expression of a specific allele within the non-redundant gene plays an important role in maintaining the high dosage of gene products. In addition, we examined the proportion of each gene group among the top expressed genes (from top 50% to top 1%) under normal (pre-treatment), biotic, and abiotic stress conditions (Fig. 4E). Across all conditions, we found that gene groups with four allelic copy numbers were disproportionately represented among top expressed genes. Specifically, the tetra-copy core genes were the most prominent, accounting for an average of 38.6% among top expressed genes across the various expression cutoffs (from top 50% to 1%). This represents a significant 1.80-fold enrichment over their 21.5% genomic baseline proportion. The four-copy genes were also consistently overrepresented (contributing an average of 21.1%), while the remaining four gene groups were underrepresented. Interestingly, the proportion of genes classified as unique-essential genes, one-copy and two-copy genes increased among the top 1% of expressed genes, indicating their prominent role in the transcriptional regulation network. In summary, the tetra-copy core genes were consistently the most abundant and enriched group among top expressed genes, suggesting that these highly conserved genes play a key role in basic biological functions and stress adaptation.

To identify genes consistently expressed under various conditions, we counted the occurrences of highly expressed genes. We analyzed a total of 124 different conditions, including variations in accessions, treatments, and tissue types. For each condition, we defined highly expressed genes (HEGs) as the top 10% of expressed genes, with the number of HEGs ranging from 2,106 to 6,225 per condition (Supplementary Data 2). From this dynamic set, we identified a core group of 203 genes that were consistently highly expressed across more than 90% of conditions (Supplementary Fig. 17, Supplementary Data 3). Gene Ontology (GO) analysis revealed that these 203 genes were enriched in responses to water deprivation, cellular response to hypoxia, and reproductive processes. Among the 203 genes, we identified 81 that belong to the tetra-copy core genes (Supplementary Data 3). While a general pattern of unbalanced number of expressed genes was observed among the four haplotypes in both leaf and root tissues (with Hap4 highest and Hap1 lowest), our targeted analysis of this specific 81-gene subset revealed no significant differences in their expression level among the haplotypes (Supplementary Fig. 18A, B, Supplementary Table 7). Most non-redundant genes were expressed in only one haplotype, and none were simultaneously expressed in all four haplotypes (Supplementary Fig. 18C). This suggests that although all four allelic copies are conserved during evolution, their homology-dependent gene silencing mechanisms or genetic regulation prevent their simultaneous expression under stressful conditions7,33,57,58.

Genetic burden of tetra-copy core genes and the long-term evolutionary fate of alfalfa
Our results showed that tetra-copy core genes are important for climate adaptation. However, the genetic burden of these conserved genes may accumulate over evolutionary timescales3. We calculated genomic sequence constraints using GERP score59. The results showed that ~1.5% (45.95 Mb) of the alfalfa genome was predicted as evolutionarily constrained sites (GERP score > 3) (Fig. 5A). We further calculated the number of deleterious variants for each of the six gene groups (Supplementary Data 4). There are 6.41 Mb (13.9%) deleterious variants located within the gene region (promoter + transcribed region). We found that the average number of deleterious variants per 1 kb region increased with the rise in allelic copy number. The mean variants number is 65.25 for unique-essential genes, 53.21 for one-copy gene, 63.50 for two-copy genes, 97.64 for three-copy genes, 142.12 for four-copy genes and 198.89 for tetra-copy core genes. A significant proportion of the variants (63.6%) within the genes are derived from coding sequence (CDS) regions. Moreover, this pattern was consistent across genes with varying allelic copy numbers (Fig. 5B).

Interestingly, the pattern of deleterious variants in unique-essential genes differed from that in other gene groups. Unlike other gene groups where CDS regions harbor the most deleterious variants, in unique-essential genes, the highest accumulation of deleterious variants was found in the 3’ untranslated region (3’ UTR). Tetra-copy core genes exhibited the highest number of deleterious variants among the six gene groups, with 80.1% of these variants located in the CDS regions. This result indicates that for every 1 kb genomic region, an average of 164 deleterious variants are located within CDS, indicating a high deleterious burden in tetra-copy core genes. We observed a consistent spatial distribution of tetra-copy core genes across the four haplotypes, with a significant enrichment clustering near the chromosomal ends (Supplementary Fig. 19), suggesting that deleterious variants are difficult to eliminate through recombination and selection.

To compare the genetic burden among haplotypes, we quantified the number of deleterious variants within each allelic copy of the non-redundant genes. Our analysis of the entire genome revealed an average of 65.24 deleterious variants per gene. Using this value as a threshold, we then classified genes exceeding this average as deleterious genes. Our analysis revealed that 20,203 (46.4%) of the tetra-copy core genes (14,776 non-redundant genes) were identified as deleterious genes. Among them, there are 5,082 in hap1, 5,045 in hap2, 4,999 in hap3, and 5,077 in hap4. Furthermore, we assessed the distribution of deleterious alleles within each of the 14,776 non-redundant genes. Our analysis revealed that 5,070 (34.3%) of these non-redundant genes carried deleterious alleles on all four of their haplotypes (16,414 genes, Supplementary Data 5). In contrast, 4.9%, 5.2%, and 9.5% of the non-redundant genes carried deleterious alleles on only one, two, or three haplotypes, respectively (Fig. 5C). Non-redundant genes carrying four deleterious allelic copies (5,070) were enriched in stress response, reproductive development, and protein ubiquitination (Fig. 5D). This indicates that many key functional genes, including those environmental response genes, have accumulated a high genetic burden across all haplotypes. Specifically, our population-level analysis revealed that autotetraploid alfalfa has accumulated ~1.5-fold more deleterious variants than its diploid ancestor, subsp. caerulea (Fig. 5E). This pattern was most pronounced in the tetra-copy core genes; while this group already had the highest genetic burden in the diploid ancestor, the burden was even greater in the autotetraploid alfalfa subspecies (Fig. 5F). This finding suggests that WGD masked deleterious variants and increased the genetic burden, especially in tetra-copy core genes3.

Overexpressed tetra-copy core gene increased adaptability in alfalfa
We selected the glycine decarboxylase gene (MsGDC) to validate the function of tetra-copy core genes, which is one of the 5,070 non-redundant genes where each of its four allelic copies carries a higher-than-average burden of deleterious variants. The MsGDC gene participates in photorespiration and plays a role in the utilization of CO2 and NH3 60. We aim to determine whether the ability to utilize CO2 or NH3 can be enhanced and lead to phenotypic changes in overexpressed alfalfa. The comparative genomic results revealed that MsGDC is present in seven Medicago taxa, with four allelic copies found in the ZM4_V2.0 alfalfa genome (Fig. 6A). The transcriptome results indicated that only one allelic copy of MsGDC was highly expressed, consistently across different accessions or lines (Fig. 6B). Unexpectedly, only the highly expressed MsGDC (gene Msa125466 in ZM4_V2.0) showed a substantial increase in response to elevated CO2 concentration61. The other three allelic copies exhibited only slight changes in response to shifts in CO2 or NH3 levels. This specific case, combined with our broader expression analyses (Fig. 4D), suggests one allelic copy maintains the major roles in fundamental biological functions and stress tolerance, while the other three allelic copies may have acquired roles in genetic regulation3.

Since the tetra-copy core genes exhibited the highest expression levels among all gene groups (Fig. 4C) and one specific allele of MsGDC showed biased high expression (Fig. 6B), we generated transgenic alfalfa with overexpressed MsGDC (Fig. 6D, Supplementary Table 8). The results showed that the transgenic GS (the acronym for the gene GDC and the promoters ST-LS1) lines exhibited a significant increase in biomass (Supplementary Fig. 20). This indicates that enhancing the expression of even a single tetra-copy core gene can significantly increase alfalfa biomass, suggesting that employing these genes to improve alfalfa yield has substantial potential. The differences in promoter lines (G lines: 35S or GS lines: ST-LS1) significantly affected plant biomass. The GS lines exhibited greater weight and height compared to the wide-type (WT) and G lines (Supplementary Fig. 20), and we also observed significant differences in gene expression levels between the promoter lines (Supplementary Fig. 21). This highlights the critical importance of promoter selection for realizing the biomass-boosting potential of this gene, as demonstrated by the superior performance of the ST-LS1 promoter62.

To further assess stress tolerance, we examined the MsGDC-overexpressing M. truncatula lines. Under normal nutrient conditions, these lines exhibited a phenotypic pattern similar to that of transgenic alfalfa, with the GS lines showing higher biomass than the WT (Fig. 6E). In contrast, when subjected to nitrogen-free treatment, the GS lines maintained good survival rates, whereas the WT plants showed poor growth and survival (Fig. 6F). Moreover, there was a significant difference in chlorophyll content between the GS lines and WT plants (Supplementary Fig. 22). These results suggest that MsGDC can enhance tolerance to nitrogen stress and improve nitrogen use efficiency in M. truncatula. Finally, we calculated the genetic burden of MsGDC. We found that most deleterious variants were present in the CDS regions at a remarkably high density. This observation was consistent for each of the four allelic copies, demonstrating that MsGDC carries a high genetic burden (Fig. 6C). To summarize, MsGDC has the potential to enhance alfalfa’s biomass and improve the nitrogen stress tolerance of M. truncatula. It could be utilized to boost the yield and develop climate-resilient alfalfa varieties. However, the high genetic burdens should be carefully managed in alfalfa breeding programs.

Methods
Genomesequencingandassembly
Genomic DNA was extracted from young leaves of the previously described Zhongmu-4 alfalfa plant44. DNA libraries were constructed using standard protocols for PacBio and Nanopore sequencing. Four PacBio Sequel II HiFi libraries were prepared to generate consensus reads. Additionally, a 100kb Oxford Nanopore Technologies (ONT) library was created to obtain ultra-long reads. The Hi-C data were obtained from our previous sequencing of Zhongmu-444. Full-length transcriptomic sequences were generated using mixed samples from leaves, buds, and flowers.

The haplotype-resolved assembly of Zhongmu-4 was generated in three steps. First, HiFi and ONT reads were used to produce primary contigs using Hifiasm (v0.18.5-r499)45. Hi-C reads were aligned to the primary contigs using the Juicer (v2.0) pipeline77. The resulting contact maps were then used with the 3D-DNA (v.201008) pipeline for an initial round of automated scaffolding78. This draft assembly was loaded into Juicebox (v1.11.08) for manual inspection and correction of large-scale scaffolding errors. To finalize the assembly, the manually curated scaffold list was then used as input to re-run the 3D-DNA pipeline, which generated the haplotype-resolved pseudo-chromosomes (ZM4_V1.5). Second, the genomic continuity and completeness were further improved based on previous version of Zhongmu-4 genome (ZM4_V1.0)44. We mapped all HiFi reads to the ZM4_V1.5 genome, and kept reads with a mapping quality score (MAPQ) of more than 55. Reads that had a mapping quality score below 55 were subsequently remapped to the ZM4_V1.0 genome. Based on the collinearity between ZM4_V1.0 and ZM4_V1.5, the remapped reads were assigned to the corresponding chromosome of ZM4_V1.5 genome. A total of 89.5% of HiFi reads were assigned to 32 pseudo-chromosomes. Third, Hifiasm (v0.19.7-r598) was used to regenerate haplotype-resolved contigs by utilizing HiFi and ONT reads, referencing the strategy of genome assembly of potato C8829. The parameters used were ‘-D 10−5 mapped_high_quality_MQ55_hifi_reads.txt --n-hap 4 --hom-cov 120’ (https://github.com/baozg/Potato_C88). As a result, we generated four groups of phased haplotype-resolved contigs (Supplementary Fig. 1). For each group, the contig reads were oriented and ordered to the chromosome-level genome based on the hap1 of ZM4_V1.5 using RagTag (v2.1.0) with default parameters79,80. The resulting chromosome-level genome was named ZM4_V2.0 and used for subsequent analysis. The genome completeness and base quality were evaluated using Merqury (v1.3)46. Phasing accuracy was evaluated by calculating switch and Hamming error rates using hap1 as the reference with a custom Python script, which compared a FreeBayes (v1.3.8)81 derived variant set against a ground truth variant set established by MUMmer4 (v4.0.1)82.

For subsp. caerulea, young leaves from one accession (PI464714) were used for sequencing. The sequencing libraries for PacBio CLR and Hi-C were constructed according to the manufacturer’s protocols. In total, we generated 207.00 Gb PacBio CLR long reads, 52.75 Gb Illumina short reads, and 204.28 Gb Hi-C data. The PacBio CLR reads were first assembled into contigs using Canu (v1.8)83 and polished with Illumina short reads using Pilon (v1.23)84. The contigs were then aligned and anchored using Juicer (v1.5)77 based on Hi-C data and subsequently assembled into chromosomes using 3D-DNA (v.180419)78.

Annotation of genes and transposable elements
The gene annotation was performed using the GWAP pipeline (https://github.com/unavailable-2374/Genome-Wide-Annotation-Pipeline)85. In summary, the ZM4_V2.0 putative genes were identified using a combination of multiple data sources, including protein annotation information from the ZM4_V1.0 genome, full-length transcriptomic sequences, Illumina short-read RNA data from leaves, stems, and flowers, as well as protein sequences from M. truncatula (Mt_A17(V5.0))35. A preliminary gene model was constructed for the putative genes, followed by additional gene searches utilizing AUGUSTUS (v3.4.0)86. The identified putative gene fragments were then filtered through the following steps: removal of genes with duplicated regions, removal of genes with invalid CDS regions (where the end position is less than the start position), removal of genes with protein sequences containing a point sign or lengths shorter than 55 amino acids, and removal of genes not supported by any evidence. InterProScan (v5.56–89.0)87 was utilized for the functional annotation of the retained genes. Furthermore, inconsistent annotation pipelines can introduce significant biases, leading to inflated estimates of dispensable gene families88. To address this issue, we re-annotated the remaining 12 genomes using a consistent pipeline.

Transposable elements analysis was conducted based on the pangenome of Medicago using our established pipeline85. In brief, the Medicago repeat family database was initially constructed using RepeatModeler2 (open-2.0.3)89 and subsequently identified with RepeatMasker (open-4.1.2)90. TEs families were identified using RepeatModeler2. Single-copy and erroneous annotations of TEs families were removed, and redundant sequences were further eliminated using NCBI-BLAST+ (v2.9.0). Unclassified repeat elements were identified with DeepTE91. Finally, the repetitive sequences were annotated using RepeatMasker.

Super-pangenome construction and gene family analysis
The Medicago graph pangenome was constructed using PGGB (v0.5.4)52 with the ZM4_V2.0 hap1 genome serving as the reference. All 13 genomes, encompassing 19 haplotypes, were included in the construction of the graph pangenome. Single chromosomes of 19 haplotypes were merged into one input fa file. The parameter of PGGB was set as follows: “-s 5000 -l 25000 -p 90 -c 1 -K 19 -F 0.001 -g 30 -k 23 -f 0 -B 10M -n 19 -j 0 -e 0 -G 700,900,1100 -P 1,19,39,3,81,1 -O 0.001 -d 100”. The VCF file was generated based on the gfa file of PGGB using vg with default parameter (v1.58.0)53. Structural variations were identified based on the VCF file results, which only include reference or alternate sequences exceeding 50bp in length. In order to achieve uniform comparisons among the 13 genomes, PCA analysis was conducted using SVs with no more than two alternate variants, and all alternate variants beyond the second were set to NA. The accession-specific SV was identified using SyRI (v1.5.4)92, SVIM-ASM (v1.0.3)93, and CuteSV (v2.1.2)94, with the hap1 of ZM4_V2.0 as the reference, and was subsequently merged with Sniffles (v2.0.7)95.

OrthoFinder (v2.5.4)51 was utilized to identify orthologous gene families across 13 Medicago genomes using default parameters. GENESPACE (v1.2.3)96 was employed to perform synteny analysis among these genomes. Collinear figures were generated using MCscan97 and R (v4.3.2). Phylogenetic trees were constructed based on single-copy orthologs identified by OrthoFinder from seven genomes: M. truncatula (A17_V5.0), M. ruthenica (Mru_zhiwusuo), M. polymorpha, M. lupulina, M. arabica, subsp. caerulea (Mca_long), and ZM4_V2.0 (hap1). The tree inference was performed with raxml-ng (v1.1)98. The divergence time between different species was estimated using PAML (v4.9)99 with default parameters. The calibration points were set at two specific points: the divergence time between the M. sativa-falcata complex and M. truncatula was 6~9 Mya, while the divergence time between the genus Medicago and Trifolium pratense was estimated to be between 13~24 Mya. The final phylogenetic tree was visualized and annotated in FigTree (v1.4.4) (https://tree.bio.ed.ac.uk/software/figtree/).

Gene conservation and expression analysis
To identify allelic copy number variation within haplotype-resolved ZM4_V2.0 genome, we performed an all-to-all BLASTP (v2.2.31)100 analysis using the protein sequences from all four haplotypes. The allelic copy number of a given gene was then defined as the total number of its homologous alleles (including itself) found across the four haplotypes. The filtering parameters were set as follows: “e-value < 1e-10, pident > 95%, (1-mismatch/length) > 95%, and Chr_query = Chr_match”. For each gene, the top-matching sequence with the highest percentage identity was identified as its corresponding allelic copy. These identified pairs were then treated as a single non-redundant gene for all subsequent analyses. To focus on allelic variation across haplotypes, we treated tandemly duplicated genes on the same haplotype as individual entities. This approach allows members of a tandem array to sort into different non-redundant gene groups, with the result that the allelic copy number for any given group is a maximum of four. The allelic copy number of a given gene was defined as one plus the number of other haplotypes on which a homologous allele was identified (via all-to-all BLASTP analysis), resulting in allelic copy numbers ranging from one (no homologs on other haplotypes) to four (homologs on all three other haplotypes). Additionally, we distinguished paralogs from allelic copies by strictly defining the latter as the gene with the highest protein similarity across haplotypes. Furthermore, this protein-based definition differs from identifying genes via presence/absence variation (PAV), which is often linked to structural variation (SV). The key distinction is that PAV is more suited to identifying newly arisen/loss genes, whereas our method offers insights into the detailed evolution of established non-redundant genes.

We determined the conservation status of each orthologous gene family (orthogroup) identified by OrthoFinder. Using the Statistics_Overall.tsv file, families were categorized as core (present in all genomes), softcore (present in 11-12), dispensable (present in 2-10), or private (present in only one). Subsequently, this family-level conservation status was assigned to every individual gene within its respective orthogroup by parsing the Orthogroups.tsv file. The relationship between gene conservation status (i.e., core vs. dispensable) and allelic copy number was analyzed using the R package dplyr. Based on this analysis, core genes with four allelic copies were defined as ‘tetra-copy core genes’, while one allelic copy dispensable genes were classified as ‘unique-essential genes’. The remaining genes were categorized based on their respective allelic copy numbers. For detailed steps, please refer to the code on GitHub (https://github.com/futurefanzhang/genomic-adaptation-of-autotetraploid-alfalfa/tree/main/step3).

The climate adaptation-associated genes were analyzed based on the VCF file generated in our previous study54 (https://ngdc.cncb.ac.cn/gvm/getProjectDetail?Project=GVM000965). The reference transfer was conducted using the following steps: First, the genome sequences of haploid ZM4 were mapped to the ZM4_V2.0 hap1 genome using minimap2 (v2.21-r1071)101. Next, chain files were created using the minimap2chain function of transanno software (v0.4.5) (https://github.com/informationsea/transanno). Finally, the VCF file was generated with ZM4_V2.0 hap1 as the reference through the liftvcf function of transanno. The analyzed sample consisted of 240 accessions exhibiting a clear geographic clustering. The LFMM and PBS analyses were conducted using the same pipeline as that reported in our previous study54. In brief, we performed a genotype-environment association analysis using the univariate LFMM in the R package LEA (v.3.11.4)102. The environmental data for 19 bioclimatic variables (resolution 2.5 arcmin) were downloaded from the WorldClim website103. Local adaptation regions and candidate genes were detected by PBS analysis using PBScan104, which involved grouping samples based on their geographic origin. Additionally, population-related analyses between subsp. caerulea and subsp. sativa were conducted using this VCF file.

The potential distribution regions of alfalfa and subsp. caerulea were inferred using a random forest model54. Briefly, this model was trained with the geographic coordinate data of accessions obtained from the USDA-GRIN database. To model the potential species distribution, we used a set of bioclimatic variables (BIO2, BIO8, BIO9, BIO15, BIO17, BIO18, and BIO19) obtained from the WorldClim version 2.1 database at a spatial resolution of 2.5 arc-minutes103. The final distribution map was generated by plotting the prediction raster from the model in R (v4.3.2) using the ‘raster’ and ‘dismo’ packages105,106.

A total of 23 Bioproject RNA-seq datasets related to biotic/abiotic stress were used for gene expression analysis in this study (22 Bioproject downloaded from NCBI, one generated from this study). The materials, repeat information, treatments, experimental design, and publication details are provided in Supplementary Data 1. For each of 23 datasets, the RNA-seq raw data were initially mapped to the ZM4_V2.0 genome using Tophat2 (v2.1.1)107. The resulting BAM files were sorted using Samtools (v1.9)108, and gene counts were calculated with HTSeq-count (v0.9.1)109. Gene expression FPKM values were calculated using a custom Perl script. Differentially expressed genes (DEGs) were identified using the DESeq2 R package (v1.42.0)110, with thresholds of |log2(FoldChange)| > 1 and an adjusted P-value (FDR) < 0.05. A gene was considered expressed if its FPKM value was greater than zero. For each study, DEGs were identified through pairwise comparisons between samples. Subsequently, only the DEG sets from comparisons yielding more than 100 genes were retained. These selected sets were then used for gene groups enrichment analysis to mitigate potential statistical bias associated with small sample sizes. The corresponding information for DEGs and gene groups was merged using the dplyr package in R (v4.3.2). Gene Ontology analysis was performed using the DAVID website (https://david.ncifcrf.gov/)111 and visualized using the ggplot2 R package112. Statistical analyses were performed to identify significant enrichments of gene groups. To test for associations between gene groups (e.g., tetra-copy genes) and conservation categories (e.g., core genes), we used a chi-squared test of independence based on 2×2 contingency tables. One row of the table contained the counts for the gene list of interest, while the second row contained the corresponding counts from the genomic baseline. For analyzing the proportion of specific gene groups within the DEGs across the 23 datasets, a one-sample Wilcoxon signed-rank test was employed. This non-parametric test was used to assess whether the median of the observed proportions was significantly different from a specific genomic baseline proportion (https://github.com/futurefanzhang/genomic-adaptation-of-autotetraploid-alfalfa/tree/main/step5).

Genetic burden analysis
The constrained sites in the ZM4_V2.0 genome were identified using the GERP++ software59, with the rejected substitutions value considered as the GERP score. Multiple whole-genome sequence alignments were performed using LASTZ (v1.04.15)113. A total of 15 angiosperm genomes were included in the alignment, specifically: grape85, poplar114, Arabidopsis thaliana115, red clover116, soybean117, pea118, chickpea119, Arachis duranensis120, M. truncatula (A17_V5.0)35, M. ruthenica40, M. polymorpha38, M. lupulina49, M. arabica48, subsp. caerulea (Mca_long), and alfalfa (ZM4_V2.0). The alignment files were merged using multiz (v11.2) and converted to FASTA format using maf2fasta (v3)121. A phylogenetic tree was constructed with single copy genes using raxml-ng (v1.1)98, and branch lengths modified based on fourfold degenerate sites using msa_view (v1.4) and phyloFit (v1.4) (http://compgen.cshl.edu/phast/phyloFit-tutorial.php). Finally, GERP scores were calculated for each position using the gerpcol function of GERP++. From this output, we first excluded all uninformative sites (those with a GERP score of zero). We then defined deleterious variants as those occurring at sites with a GERP score > 3 (rejected substitutions > 3). The detailed steps can be found on GitHub (https://github.com/futurefanzhang/genomic-adaptation-of-autotetraploid-alfalfa/blob/main/step4/gerp_score.sh).

Validation of MsGDC in transgenic experiment
The plant materials for treatment were derived from a single plant of the alfalfa cultivar Zhongmu No. 1 (Zhongmu1) and a single plant of M. truncatula (R108). The full-length CDS sequence of the candidate gene MsGDC was cloned from the Zhongmu1 genome43. The stem-leaf-specific expression promoter ST-LS1 sequence was amplified and cloned based on the sequence from NCBI GenBank (accession No. X04753.1)122. The correctly sequenced vectors expressing MsGDC, driven by the 35S and ST-LS1 promoters, were transformed into Agrobacterium tumefaciens competent cells. Alfalfa leaves were transformed via Agrobacterium tumefaciens infection, and the transgenic plantlets were screened on a solid medium containing glyphosate (4 mg·L⁻¹) through tissue culture. The expression level of MsGDC among its four allelic copies was calculated using a publicly available RNA-seq dataset61.

Two-month-old Zhongmu1 plants overexpressing MsGDC were grown in an artificial climate chamber (16 hours light at 24˚C, 8 hours dark at 22˚C, 60% humidity). We quantified leaf MsGDC expression via real-time PCR (7300 Real Time PCR System, ABI) and measured several growth-related phenotypes, including plant height, fresh weight, and dry weight. Three biological replicates were conducted for each test, with one plant selected for each replicate.

One-month-old wild-type (WT) and transgenic M. truncatula seedlings overexpressing MsGDC were selected and subjected to a nitrogen-free nutrient solution for an additional month. Plant phenotypes were examined, and 0.1 g of leaves were collected for chlorophyll extraction. Leaves from the N-free M. truncatula plants were immersed in 10 mL of 95% ethanol and kept in the dark for 48 hours until they were completely discolored. Absorbance values at 665 nm and 649 nm were measured using a NanoPhotometer instrument (Germany). Chlorophyll a, chlorophyll b, and total chlorophyll contents were subsequently calculated123.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Dear Zhang F., Wei C., Shi X et al.,\nThanks for your manuscript: Medicago super-pangenome reveals adaptive advantages and evolutionary constraints in\nautotetraploid alfalfa. The studies were well conducted, and the MS was well drafted. I enjoyed reading your manuscript,\nespecially the results and discussion with respect to gene retention, climate adaption, and genetic burden.\nI have some comments and suggestions that I hope could help you improve your manuscript. Some minor concerns should\nbe addressed, and some potential large structural variations could be cross-checked or validated.\n1. Line 142-164 (and methods): It is appreciated that you have taken care of the polypoid assemblies by walking through\ndetailed mapping, phasing, and HiC assistance. As illustrated in your Fig.1A, 13 genomes across seven species show good\nsyntenic alignments; however, the alignments caught my attention. Taking into account the statics in table 1, some of the\ngenomes might not in a good shape. The large translocations (which look like large chromosome segments), and inverted\nchromosomes will need further investigation. It looks like some 5th chromosomes (dark blue) of Mara, Mlup,\nMru_zhiwusuo/Mru_landa are inverted. Similar observations were made for other chromosomes as well, e.g., the first\nchromosomes for Mlup and Mru_zhiwusuo. PLEASE double-check all genome assemblies. This might be very easily fixed\nby inverting short/long arms. It was noted that 11 of the other genomes were downloaded from the public domain, but it\nshould not stop you from updating them and making them right. Your update should be clarified in your manuscript.\n2. line162-164: Did you look into the large translocation across C4 and C8 (suppl.Fig4). A way to look at it could be long-\nread mapping, as you have both pacbio and nanopore reads. Just curious if that is a real translocation or miss-assemble\nchromosome.\n3. line 180-182: The pan-genome graph using hap1 (ZM4_V2.0) as the reference helped you get around the potential mis-\nassemblies (mentioned above). You should present more regarding the graph, and how it leverages your research. As you\nhave mentioned the variation above 50bp in length, did you look into the large structure variations? I would like to see some\nillustrations for the pangenome graph, and also a closer look, e.g., gene MsGDC. What are the strengths and challenges of\nbuilding the genus/multispecies pangenome graph?\n4. Line 232-250: the beyond analyses are based on your terms, tetra-copy core genes, unique essential genes, etc. It has\nbeen demonstrated that gene retentions, loss and enrichment for biological purposes in some other species, and still worth\ninvestigation in alfalfa. The investigation of genes’ response to climate change and local adaption is valuable.\n5. Line294-296: the expression data supported your analysis and conclusion in your manuscript. You found most genes\nwere expressed in HAP4 (Supple.Fig.7). However, there might be different stories in other species. This triggers me to\nwonder about the evolution of Medicago. Still it was not demonstrated clearly the origin of the sativa autotetraploid, which\nlooks like the most recent ancestor could be caerulea or truncatula. Assume ZM4 has four haplotpyes originated from\ncaerulea, have you analyzed the divergence of four haplotpyes, which might give a clue of why hap4 present most of\nexpressed genes?\n6. Line302-342: I enjoyed reading the genetic burden, constraints, and deleterious variants. I just wondered maybe tetra-\ncopy core genes could be more tolerant to deleterious variants, as they have more copies. In addition to this, is the threshold\nto slow for deleterious variants? It looks full of Deleterious variants in MsGDC genes (Fig 5C), and it functions very well.\n7. Line366: Do the figures show similar expressions? Suppl.Fig10 shows the differences to me. I would suggest conducting\none2one comparison for wild type and edited lines for stem and leaf separately. This led to one of your statements: “despite\nthe similar levels of gene expression.”\n8. Open discussions with the authors. In this manuscript, you deeply investigate the gene retention, number of copies,\nexpressions, and genetic burden and validate the gene’s biological function in Medicago, which includes autotetraploid\nalfalfa. The research is important and valued. However, it is solid for Medicago, and maybe some types of polyploidies, but it\nmight not be able to represent all cases of polyploidy. In some other species, the gene orientation and expression might be\ndifferent cases, e.g., bread wheat, teosinte, or wild rice (polyploids, not refs to wild Asian rice), where the expressions show\ndifferent patterns in different functional groups. This study for sure provides evidence of enhancing genetic analysis and\npotential applied value for autotetraploid and alfalfa species, but I would still be more careful with the overstatements of\n‘facilitating the genetic studies and breeding programs of polyploid crops.’\nBest,", "reviewer_2": "In the manuscript titled \"Medicago Super-Pangenome Reveals Adaptive Advantages and Evolutionary Constraints in\nAutotetraploid Alfalfa,\" the authors present a valuable study by sequencing a diploid Medicago genome and improving a\ntetraploid alfalfa genome. By integrating publicly available Medicago genome sequences, they constructed a pangenome.\nWhile the study provides two high-quality alfalfa genomes and important insights, several major concerns need to be\naddressed:\n1. Insufficient Details in Materials & Methods\nThe methods section lacks sufficient information for replication. Specifically:\n• The definitions of unique/essential genes are unclear.\n• The criteria for classifying genes as core, softcore, dispensable, and private need to be explicitly stated.\n• The parameters used for PGGB should be detailed.\n• The methodology for DEG analysis requires further clarification.\n• Tables should be properly numbered.\n2. Concerns with Differential Gene Expression (DEG) Analysis\nThe approach to DEG analysis is unclear and needs justification:\n• The authors must specify which tissues and conditions were compared.\n• Did the study re-analyze all 23 published RNA-seq experiments using the new reference genome? If so, this could be\nacceptable. However, which 124 conditions were compared? This information is not found in Supplementary Table 4.\n• RNA-seq data is highly specific to plant type, tissue, and environmental conditions. Combining datasets from different\nstudies without proper batch effect correction can lead to misleading conclusions. Publicly available datasets from different\nsources should not be used to infer expression patterns in unrelated plants/populations.\n3. Pangenome Analysis Using OrthoFinder\nBias in gene family classification may arise due to inconsistent genome annotation:\n• Of the 13 genomes analyzed, only four were annotated using the same pipeline, while the rest relied on publicly available\ndata, each annotated differently.\n• Studies (e.g., Lin et al. 2014) have shown that inconsistent annotation pipelines can introduce significant bias, inflating\nestimates of dispensable gene families.\n• To ensure fair comparisons, all genomes should ideally be annotated using the same pipeline. If homology-based\nannotation is used, RNA-seq data should come from the same tissues and conditions.\nTo improve reliability, I strongly recommend re-annotating publicly available genomes using the same pipeline as the other\nfour. Alternatively, apply strict filtering, and include only functionally annotated protein-coding genes across all genomes.\nThis will reduce bias and ensure more accurate gene family classification.\nSpecific Comments\n• Line 54: The last sentence of the abstract is incomplete. Consider deleting “that.”\n• Line 152: Is ZM4_V1.0 from Long et al. 2022? Is ZM4_V1.5 the current version?\n• Line 170: How does this version improve upon the previous one? Please compare the number of gaps and genes between\nZM4_V1.0 and ZM4_V2.0.\n• Lines 222-225: Does this refer to the ZM4_V2.0 genome? Define unique genes clearly. How was the number 109,715\ndetermined? Provide BUSCO analysis results on single-copy and duplicated genes.\n• Lines 225-229: How many four-copy genes are dispensable? How many single-copy genes are core? The statement does\nnot confirm non-random gene retention, duplication, and loss as suggested in the next sentence.\n• Lines 247-248: The sentence is incomplete and unclear.\n• Lines 256-257: The tetra-copy core genes represented the largest proportion of what?\n• Line 283: Which 124 conditions were compared? This is not detailed in Supplementary Table 4. Please refer to my second\nmajor concern.\n• Line 330: How were gene pairs selected for analysis?\n• Line 365: How was plant biomass affected?\n• Lines 492-493: The sentence “Additionally, haploid pseudo-chromosomes were generated using the same pipeline\n(ZM4_V1.5)” appears redundant.\n• Line 539: Specify the version and parameters used for the PGGB pipeline. How was the VCF file obtained and filtered?\n• Line 544: Were publicly available protein sequences used for OrthoFinder in other genomes? Were genes filtered for\nORFs or pseudogenes? See my concern regarding inconsistent annotations.\n• Line 564: Provide more details on how orthologous gene families were aligned with copy number information.\n• Line 571: Clarify the DEG analysis methodology, including specific comparisons made.\nAddressing these concerns will strengthen the manuscript’s validity and reproducibility.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "Comments to the Authors\nThe authors present a study focused on constructing a Medicago super-pangenome using newly assembled autotetraploid and diploid M. sativa genomes alongside published genomic data. They identify tetra-copy core genes in alfalfa, explore their roles in climate adaptation and stress responses, discuss their evolutionary constraints, and validate the functional relevance of one such gene through overexpression. While the topic is highly relevant for understanding polyploid adaptation and breeding applications, the manuscript has several critical limitations in data integration, methodology, and interpretation that require substantial revision. Below are my major and minor concerns.\nMajor points:\n1. The study heavily relies on published genomic data (7 species, 13 genomes), but only 5 genomes were annotated using the same method. This inconsistency in annotation pipelines may introduce significant bias, particularly for identifying tetra-copy core genes, which form the foundation of all downstream analyses. For example, differences in gene prediction tools can inflate or obscure ortholog numbers. The authors should re-annotate all included genomes using the same pipeline to ensure comparability, and maybe also to ddress how annotation heterogeneity might have influenced the tetra-copy core gene set (e.g., sensitivity analysis). Additionally, the term \"comprehensive super-pangenome\" is misleading, as only 7 of ~90 Medicago species are included.\n2. The inclusion of distantly related Medicago species (e.g., M. truncatula and M. polymorpha; genome size ~450 Mb vs. ~850 Mb in diploid alfalfa) introduces evolutionary and genomic divergence that confounds tetra-copy core gene identification. These species differ in life history (annual vs. perennial) and mating systems (selfing vs. outcrossing), which likely influence gene family dynamics. Would it be more reasonable to restrict tetra-copy core gene analysis to closely related species (e.g., M. sativa diploids, tetraploids, and their wild relatives like M. falcata)? Alternatively, it would be possible to perform lineage-specific gene retention/loss analyses to disentangle polyploid-specific signals from broader phylogenetic noise.\n3. The super-pangenome framework is underutilized. Beyond identifying tetra-copy genes, no deeper exploration (e.g., the role of structural variation, presence-absence variation across species) is provided.\n4. The authors directly reuse \"climate-adaptive genes\" identified in their prior work based on LFMM and PBS methods applied to global alfalfa accessions. However, the alfalfa’s global distribution is heavily influenced by human-mediated dispersal and selection. Consequently, how many of these \" climate-adaptive genes\" identified via local climatic factors are actually driven by artificial selection? In their earlier study, the authors also found that nearly 100 of these \" climate-adaptive genes\" were attributed to introgression. Therefore, directly using this gene set to investigate its association with different gene groups in alfalfa is inappropriate, as it conflates natural adaptation signals with human-driven selection and introgression effects. Therefore, the authors should re-analyze the climate-adaptive signals to ensure the reliability of the conclusions.\nMinor points:\nResults\nL184–189: Comparing total genome size (including all four subgenomes) in tetraploid alfalfa against single diploid genomes is methodologically flawed. Use haploid genome equivalents (1C values) for fair comparison.\nL205: The text mentions 11–12 species, but only 7 are included. Clarify this discrepancy.\nL208–209: The claim that \"more homologous gene families indicate higher genome quality\" is unfounded. Gene family counts depend on annotation consistency and evolutionary divergence, not assembly quality.\nL245: The 1671 climate adaptation-related genes were obtained based on the union gene sets of LFMM and PBS analyses, not the results of GWAS.\nL248: Define \"genomic background content\" and clarify how the 19.4% was calculated.\nL270: What is \"normal condition\"?\nL339–341: Figure 4F compares two species but refers to \"highest expression\"—this is ambiguous.\nL348: Should be \"genomic\".\nMethods\nL528–529: Justify why uniform annotation was not applied to all species.\nL531: Clarify whether repeat annotation was performed before or after pangenome construction.\nL548: Detailed information of divergence time estimation (e.g., software, calibration points) are needed.\nL551: The authors should clearly distinguish \"duplicated genes\" (paralogs) from \"allelic genes\" (homologs within a polyploid).\nL581–590: I am completely confused about the genetic burden analysis. For example, how is the phylogenetic tree constructed related to the GERP scores? How are the GERP scores calculated? Are GERP scores calculated for all 15 species, or only for a subset of them?\nFigure 4F\nThe diploid alfalfa genome should have two haplotypes. How are \"three-copy\" and \"four-copy\" genes defined here? Provide a methodological explanation.", "reviewer_5": "Polyploid plants often outperform their diploid sisters, but when they do, the genetic and mechanistic bases for their advantages remain unclear. In this study, researchers generated a haplotype-resolved autotetraploid alfalfa genome and incorporated it into a super-pangenome spanning 13 genomes from seven Medicago species. They report dramatic gene content variation, with only about 20% of genes consistently present across all four haplotypes, a striking and rather surprising finding. A key subset, the tetra-copy core genes, appeared important for climate adaptation and stress responses, though they also carried a heavier genetic burden. While the study suggests that polyploidy can enhance traits like biomass and nitrogen use efficiency, some aspects require better evidence precisey due their importants. In particular, the quality of the tetraploid assembly and its true biological significance remain open questions that require further validation.\nMajor Concern:\nMy single major concern about this paper is that autotetraploid assembly is not sufficiently described in terms of quality. What I see as the major finding of the paper hangs completely on this, so this is important.\nMy concern as someone is that Hifiasm is not built for polyploid (especially autopolyploid) genome assembly: the --n-hap 4 tag is still rather experimental (esp as of v0.18). It is not that I doubt strongly the results, but better documentation (meaning detailed metrics into the supplemental information) describing the process and quality of the assembly is required. Such documentation of correct assembly is crucial for us to believe the key finding figure 3A and therefore this is a major concern.\nThe core issue is this: the finding that half the genes are present only in one haplome is highly striking and thus begs further validation. The contig N50 is not great, but probably not so poor to produce such serious widespread artifacts to invalidate the key finding. That said, greater detail of the assembly process and statistics is needed.\nIn addition, confidence in the finding can also be productively provided by showing gene neighbourhood-level examples of widespread presence/absence variation as well. Ideally this would be shown for several (or better, 8 to 12) regions e.g. in the SI, but at least one large region in the main text figures, demonstrating clearly that presence/absence variation of genes among the four haplomes. To drive this home, a track for each haplome of hifi read mapping depth and gaplessness at least in the regions shown would solidify the finding.\nFinally, relating also to this concern: the wider significance of this finding of what I would term ‘autopolyploid fractionation’ would be very helpful for the reader. Can we expect this to be a common phenomenon, or is this line-specific? The lack of information about the Zhongmu-4 line used is not helpful in this regard. E.g. is it highly heterozygous and/or does it have a history of inbreeding? Could this finding be partly explained by historical inbreeding? Ideally Tajima’s D if the data can support this.\nOther than this single major concern above, my remarks below can be classified as minor, but also hopefully helpful to improve this paper.\n-I’m unsure the title makes sense: in what way is the superpangenome reveal “adaptive advantages?” What is shown is very interesting but really what is shown are correlates of potential adaptive nature, not ‘revealing adaptive advantages’. Such language is overblown. The ‘evolutionary constraints’ also are less clear. Constraint is usually thought of in a different way: can this be better phrased?\nThe introduction is very good, but with some room for improvement. I suggest at line 76 and following that the presentation is overly adaptationist and should cite the dark side of polyploidy. For example Monnahan et al\n(https://doi.org/10.1038/s41559-019-0807-4) describe a decrease in the efficiency of purifying selection and increased increase in load. Mention of this would give a more balanced introduction and presage your interesting findings on burden.\n-Lines 146 and following: I think it would be more effective and intuitive to state the sequencing depth versus the haploid genome size, rather than the four haplotypes together. It took several passes reading this to understand you meant per haplotype (an unusual convention); at first one has the impression that, for example, the nanopore sequencing is only six-fold total, rather than six-fold versus each of the four alleles. That is, you risk misunderstanding for the reader and accidental underselling the excellent data the way it's currently presented.\n-Line 216, the retention of core genes enhanced climate adaptation in alfalfa subheading is indicating a circumstantial finding. It is difficult to make the statement that the retention of the core genes actually did enhance climate adaptation. I would say is correlated with climate adaptation rather than enhanced it.”\n-Line 231 “fully functional”: Nothing here is demonstrated to be ‘fully functional’; this is an inappropriate statement. Indeed, functional in line 233 and unique essential both seem overstated and have not been subject causal testing, only correlation.\nPlease back off on this language. It is unnecessary.\n-Line 249, here, for example, it is fair, the way it is stated, the finding ‘suggested’, that tetra-copy core genes are essential.\n-Lines 262 to 264, the first half is a fragment, join sentences but also simplify.\n-Line 266: “was gradually decreased”: better to say something like, ‘exhibits a quantitative gradation” or so.\n-Line 284: Supplemental Table 4 is simply a list of 23 references. The main text says “we analyzed a total of 124 different conditions”. Please spell these out in the table and make it more clear for the reader what you are talking about.\n-Line 292: “hap2 and hap3 was 25%”. This is completely unclear. Please better show the metrics in question and clarify in the text what this is.\n-Line 318: this on the following line are very unclear to me: can you please rephrase and elaborate this in the text?\n-Line 326: to define genes that accumulate more delirious variants than the genomic average as ‘deleterious genes’ seems extreme.\n-Line 385-386: you do not confirm that tetra-copy core genes ‘are essential for enhancing polypoid adaptability’. Please tone back the language: you have interesting findings but you don’t actually confirm either ‘enhancing adaptability’ or demonstrate that those genes are ‘essential’ for that in particular. The language is overblown unnecessarily.\n-Line 394: ‘past one century’: remove ‘one’\n-Line 412: this does not seem the case at all! You show a clear gradient and not expression so dominated by single haplotype. Is this calling to Figure 3G or am I missing a callout to a panel in the figures? If so, it should be pointed to even in the discussion.\n-Line 419: this is very interesting this part of the discussion.\n-Line 437: ‘TCP1’ is incorrect: it is ‘TPC1’ one. Additionally that paper does not show ‘heterozygote advantage’. It would be more effective to use that paper and others on autopolyploid population genomics, such as the Hamala or Monnahan (mentioned above) to indicate that increased load also allows masking of potentially useful alleles for adaptation to, for example here, serpentine.\n-Line 447: this distribution pattern and supplemental figure 8 is completely opaque. Please it would be better to show the relative enrichment of tetra-copy core genes versus all genes in across the genome in supplemental figure eight: not just the pattern of the tetra-copy core genes. I expect this pattern that you show closely reflects all genes not just tetra-copy core genes. That is, of course tetra-copy core genes very crudely as shown tend to cluster at the ends of chromosomes relative to the centromeres: this is because all genes do! What would help your case is to show the relative density of tetra-copy genes, normalized vs overall gene density in the figure. This would then allow the hypothesis to be tested that core genes are, in fact, differently distributed from the rest of the genes (e.g. harder to purge or so if this is what you hypothesise).\n-Line 463: here again it would be very good to relate your findings to Monnahan and Hamala MSs, indicating that this paper is congruent with findings that load is increased in autopolyploids, and indeed that purifying selection is perhaps less efficacious.\n-Line 470: this is very interesting.\n-In table one, the BUSCO scores are incomplete. I am intetested to see the C, D and S metrics.\nOverall, this work provides exciting insights, highlighting extensive gene presence/absence variation and the potential role of tetra-copy core genes in stress adaptation. However, the reliability of the key findings hinges on the quality of the tetraploid genome assembly, which requires more detailed documentation and validation. Addressing this, along with refining overstatements in the text and improving clarity in data presentation, will strengthen the study's impact."}, "author_response": "Response to Review Comments\nContents\nPart I: Summary of responses to the major questions: ............................................................. 2\nPart II: Point by point response to specific comments: ............................................................. 3\nReviewer #1 (Remarks to the Author): ....................................................................................... 3\nReviewer #2 (Remarks to the Author): ....................................................................................... 7\nReviewer #3 (Remarks to the Author): ..................................................................................... 12\nReviewer #4 (Remarks to the Author): ..................................................................................... 12\nReviewer #5 (Remarks to the Author): ..................................................................................... 17\nWe sincerely appreciate the interest shown by the editor and the five reviewers. Before presenting\na detailed point-by-point response in Part II, we first address four major issues collectively:\nFirst, we have reannotated all the genomes using the same pipeline. The number of annotated genes\nand the BUSCO assessment have shown improvement compared to the raw version.\nSecond, the chromosome numbers, order, and orientation were modified based on the model plant\ngenome of M. truncatula A17 (5.0). The genomic collinearity in Medicago has been enhanced\ncompared to the raw version.\nThird, we have provided more detailed information regarding genome assembly, gene copy\nnumber variation, identification of tetra-copy genes, and differentially expressed genes in the\nmethods section, the supplemental file, and on GitHub\n(https://github.com/futurefanzhang/genomic-adaptation-of-autotetraploid-alfalfa/tree/main).\nFourth, we have added further analyses to support our conclusions, including SV analysis of the\npangenome, conservation studies of tetra-copy genes, and a reanalysis of climate-related genes.\nIn conclusion, we have provided detailed explanations in response to the reviewers' comments. All\nupdates in the manuscript are highlighted with tracked changes. Additionally, we have revised the\nmanuscript, supplemental figures and tables.\n(The reviewer’s comments are shown in black, and the authors’ responses are displayed in blue.)\nReviewer #1 (Remarks to the Author):\nDear Zhang F., Wei C., Shi X et al.,\nThanks for your manuscript: Medicago super-pangenome reveals adaptive advantages and\nevolutionary constraints in autotetraploid alfalfa. The studies were well conducted, and the MS\nwas well drafted. I enjoyed reading your manuscript, especially the results and discussion with\nrespect to gene retention, climate adaption, and genetic burden.\nI have some comments and suggestions that I hope could help you improve your manuscript. Some\nminor concerns should be addressed, and some potential large structural variations could be cross-\nchecked or validated.\nResponse: We thank the positive comments from the reviewer and appreciate the reviewer’s time,\npatience and constructive comments. We had revised the manuscript accordingly and we hope the\ncurrent version is more suitable for publication.\n1. Line 142-164 (and methods): It is appreciated that you have taken care of the polypoid\nassemblies by walking through detailed mapping, phasing, and HiC assistance. As illustrated in\nyour Fig.1A, 13 genomes across seven species show good syntenic alignments; however, the\nalignments caught my attention. Taking into account the statics in table 1, some of the genomes\nmight not in a good shape. The large translocations (which look like large chromosome\nsegments), and inverted chromosomes will need further investigation. It looks like some 5th\nchromosomes (dark blue) of Mara, Mlup, Mru_zhiwusuo/Mru_landa are inverted. Similar\nobservations were made for other chromosomes as well, e.g., the first chromosomes for Mlup\nand Mru_zhiwusuo. PLEASE double-check all genome assemblies. This might be very easily\nfixed by inverting short/long arms. It was noted that 11 of the other genomes were downloaded\nfrom the public domain, but it should not stop you from updating them and making them right.\nYour update should be clarified in your manuscript.\nResponse: Thank the reviewer for raising this important point. We reanalyzed and modified the\ninverted chromosome. Furthermore, we reannotated all genomes using the same pipeline (Table 1,\nLine 210-215). Good syntenic alignments can be observed in Medicago species (Figure 1A).\n2. line162-164: Did you look into the large translocation across C4 and C8 (suppl.Fig4). A way to\nlook at it could be long-read mapping, as you have both pacbio and nanopore reads. Just curious\nif that is a real translocation or miss-assemble chromosome.\nResponse: Thank you for raising this question. The large translocation is observed exclusively in\nM. truncatula A17, measuring over 10 Mb (Figure 1A). Tracking this translocation using PacBio\nor Nanopore reads is challenging; therefore, we utilized Hi-C data to verify its existence (the raw\nusing Juicer, we observed strong Hi-C interaction signals within each chromosome (Response Fig.\n1). This suggests that the genome quality is sufficiently high and that the translocation is indeed\nreal.\nResponse Fig. 1 The hic interaction matrices of M. truncatula A17\nReference:\nLi A, Liu A, Wu S, et al. Comparison of structural variants in the whole genome sequences of two\nMedicago truncatula ecotypes: Jemalong A17 and R108[J]. BMC Plant Biology, 2022, 22(1): 77.\n3. line 180-182: The pan-genome graph using hap1 (ZM4_V2.0) as the reference helped you get\naround the potential mis-assemblies (mentioned above). You should present more regarding the\ngraph, and how it leverages your research. As you have mentioned the variation above 50bp in\nlength, did you look into the large structure variations? I would like to see some illustrations for\nthe pangenome graph, and also a closer look, e.g., gene MsGDC. What are the strengths and\nchallenges of building the genus/multispecies pangenome graph?\nResponse: We appreciate the reviewer’s suggestion. We have conducted additional analyses related\nto SVs. More than 70,000 SVs have been identified between ZM4_V2.0_hap1 and other genomes\nwithin the M. sativa-falcata complex. However, we detected only a small number of SVs across\nspecies (Figure 2E, F, Line 264-277). We took a closer look at Chr1, as the MsGDS is too small\nfor effective SV detection. Our analysis revealed several large SVs among different species\n(Supplemental Fig. 11 and 12). We believe that the strength of the genus pangenome lies in its\nabundant gene diversity, while the challenge remains that mapping SVs is quite difficult.\n4. Line 232-250: the beyond analyses are based on your terms, tetra-copy core genes, unique\nessential genes, etc. It has been demonstrated that gene retentions, loss and enrichment for\n4\nbiological purposes in some other species, and still worth investigation in alfalfa. The investigation\nof genes’ response to climate change and local adaption is valuable.\nResponse: We thank the reviewer for this suggestion. We agree that gene retention, loss, and\nenrichment play crucial roles in alfalfa’s adaptability. This study aims to identify key factors\ncontributing to climate adaptation, and the results showed tetra-copy core genes play the most\nimportant role for climate adaptation (Figure 4A, B, and E). Your suggestion provides a good\ndirection for our future study. To enhance understanding, we have added relevant content in the\ndiscussion section to clarify the effect of gene retention, loss and enrichment on local adaptation\n(Line 523-525).\n5. Line294-296: the expression data supported your analysis and conclusion in your manuscript.\nYou found most genes were expressed in HAP4 (Supple.Fig.7). However, there might be different\nstories in other species. This triggers me to wonder about the evolution of Medicago. Still it was\nnot demonstrated clearly the origin of the sativa autotetraploid, which looks like the most recent\nancestor could be caerulea or truncatula. Assume ZM4 has four haplotypes originated from\ncaerulea, have you analyzed the divergence of four haplotypes, which might give a clue of why\nhap4 present most of expressed genes?\nResponse: We thank the reviewer for this constructive suggestion. We analyzed the divergence of\nfour haplotypes by phylogenetic tree (Supplemental Fig.10). Our analysis revealed that the four\nalfalfa haplotypes are polyphyletic, which strongly suggests they have multiple independent\norigins from different diploid subsp. caerulea ancestors. Furthermore, the distinct phylogenetic\nposition of hap4 suggests it might possess unique genetic attributes, potentially leading to specific\nexpression effects. Recognizing the importance of this result, we added related contents in the\nmanuscript (Line 242-245).\n6. Line302-342: I enjoyed reading the genetic burden, constraints, and deleterious variants. I just\nwondered maybe tetra-copy core genes could be more tolerant to deleterious variants, as they have\nmore copies. In addition to this, is the threshold to slow for deleterious variants? It looks full of\nDeleterious variants in MsGDC genes (Fig 5C), and it functions very well.\nResponse: Thanks for rising this important question. We agree that tetra-copy core genes should\nexhibit greater tolerance to deleterious variants due to the presence of four copies (Figure 5F). The\nthreshold for the GERP value was set to greater than three (Lines 772-773), with the expected\nsubstitution under neutrality being 4. This means that an observed substitution at a given position\nwith fewer than one substitution was classified as a deleterious variant (4-3=1). So the threshold\nwas appropriate for identification of deleterious variants. We would like to clarify that the MsGDC\ngene itself is not saturated with deleterious variants. We apologize for the previous figure display.\nA new figure has been redrawn (Figure 6C).\n5\n7. Line366: Do the figures show similar expressions? Suppl.Fig10 shows the differences to me.\nI would suggest conducting one2one comparison for wild type and edited lines for stem and leaf\nseparately. This led to one of your statements: “despite the similar levels of gene expression.”\nResponse: Thank you for highlighting this important issue. We have redrawn the figure and found\nthat there are different expression levels among the wild type, G lines, and GS lines (Response Fig.\n2). We have revised our statements and figure (Line 477-478, Supplementary Fig. 19).\nResponse Fig. 2 Relative MsGDC expression levels in WT, G lines, and GS lines.\n8. Open discussions with the authors. In this manuscript, you deeply investigate the gene retention,\nnumber of copies, expressions, and genetic burden and validate the gene’s biological function in\nMedicago, which includes autotetraploid alfalfa. The research is important and valued. However,\nit is solid for Medicago, and maybe some types of polyploidies, but it might not be able to represent\nall cases of polyploidy. In some other species, the gene orientation and expression might be\ndifferent cases, e.g., bread wheat, teosinte, or wild rice (polyploids, not refs to wild Asian rice),\nwhere the expressions show different patterns in different functional groups. This study for sure\nprovides evidence of enhancing genetic analysis and potential applied value for autotetraploid and\nalfalfa species, but I would still be more careful with the overstatements of ‘facilitating the genetic\nstudies and breeding programs of polyploid crops.’\nallopolyploids may have different stories. We have modified our abstract and discussion\naccordingly (Line 66-67, 504-506).\nReviewer #2 (Remarks to the Author):\nIn the manuscript titled \"Medicago Super-Pangenome Reveals Adaptive Advantages and\nEvolutionary Constraints in Autotetraploid Alfalfa,\" the authors present a valuable study by\nsequencing a diploid Medicago genome and improving a tetraploid alfalfa genome. By integrating\npublicly available Medicago genome sequences, they constructed a pangenome. While the study\nprovides two high-quality alfalfa genomes and important insights, several major concerns need to\nbe addressed.\nResponse: We appreciate the reviewer’s time, patience and constructive comments that help us to\nimprove this manuscript. We have revised the manuscript accordingly and we hope the revised\nmanuscript is more suitable for publication.\n1. Insufficient Details in Materials & Methods\nThe methods section lacks sufficient information for replication. Specifically:\n1.1 The definitions of unique/essential genes are unclear.\nResponse: We apologize for the lack of clarity. We have added more detailed information in the\nMethods section (lines 720-727). Additionally, to assist readers in replicating our analysis, we have\nincluded the detailed code on GitHub (https://github.com/futurefanzhang/genomic-adaptation-of-\nautotetraploid-alfalfa/tree/main/step3).\n1.2 The criteria for classifying genes as core, softcore, dispensable, and private need to be\nexplicitly stated.\nResponse: Thank you for your suggestion. We have added more explicit statements (Line 717-\n719).\n1.3 The parameters used for PGGB should be detailed.\nResponse: Thanks for your suggestion. We have added more detailed information (Line 676-678).\n1.4 The methodology for DEG analysis requires further clarification.\nResponse: Thanks for your suggestion. We have added more detailed information related to DEG\nanalysis (Line 752-756).\n1.5 Tables should be properly numbered.\nResponse: Thanks for your suggestion. We have revised table numbers in the manuscript.\n2. Concerns with Differential Gene Expression (DEG) Analysis. The approach to DEG analysis is\nunclear and needs justification:\nin the Methods section (lines 752-756). Additionally, the detailed code has been released on\nGitHub (https://github.com/futurefanzhang/genomic-adaptation-of-autotetraploid-\nalfalfa/tree/main/step5).\n2.1 The authors must specify which tissues and conditions were compared.\nResponse: Thanks for your suggestion. We have added the related information in Supplementary\nTable 7 to show the tissue, stress condition, DEG numbers and gene group information.\n2.2 Did the study re-analyze all 23 published RNA-seq experiments using the new reference\ngenome? If so, this could be acceptable. However, which 124 conditions were compared? This\ninformation is not found in Supplementary Table 4.\nResponse: Thanks for raising this question. The 23 published RNA-seq experiments were\nreanalyzed using our new reference ZM4_V2.0 (added more sentences in Line 747). The 124\nconditions were supplied in Supplementary Table 6.\n2.3 RNA-seq data is highly specific to plant type, tissue, and environmental conditions. Combining\ndatasets from different studies without proper batch effect correction can lead to misleading\nconclusions. Publicly available datasets from different sources should not be used to infer\nexpression patterns in unrelated plants/populations.\nResponse: Thanks for raising this import issue. We acknowledge that combining different datasets\nmay lead to biased expression patterns. Our objective was to determine which gene group exhibited\nthe highest expression levels, with each dataset containing the six gene groups we aimed to\ncompare. Therefore, DEGs and FPKM values were analyzed within each dataset, allowing us to\ncalculate the percentage of each gene group (Figure 4B, E). To validate the expression patterns\npresented in Figures 4C and 4D, we analyzed the gene expression profiles from our sequenced\ndataset (Response Fig. 3) and found that they exhibited similar patterns, thus confirming the\nreliability of the original results.\nexpression pattern of six groups genes. B, Expression levels of the four alleles in tetra-copy core\ngenes.\n3. Pangenome Analysis Using OrthoFinder\nBias in gene family classification may arise due to inconsistent genome annotation:\n3.1 Of the 13 genomes analyzed, only four were annotated using the same pipeline, while the rest\nrelied on publicly available data, each annotated differently.\nResponse: Thank you for highlighting this important question. We have re-annotated the remaining\ngenomes (Table 1). The number of annotated genes and the BUSCO assessment have shown\nimprovement compared to the raw version.\n3.2 Studies (e.g., Lin et al. 2014) have shown that inconsistent annotation pipelines can introduce\nsignificant bias, inflating estimates of dispensable gene families.\nResponse: Thanks for this suggestion. We have added related contents in the manuscript (Line\n660-663).\n3.3 To ensure fair comparisons, all genomes should ideally be annotated using the same pipeline.\nIf homology-based annotation is used, RNA-seq data should come from the same tissues and\nconditions.\nResponse: Thanks for this suggestion. All genomes have been re-annotated using the same pipeline.\nThe RNA-seq related question have been answered in Question 2.3\n4. To improve reliability, I strongly recommend re-annotating publicly available genomes using\nthe same pipeline as the other four. Alternatively, apply strict filtering, and include only\nfunctionally annotated protein-coding genes across all genomes. This will reduce bias and ensure\nmore accurate gene family classification.\npipeline (Line 210-215, Table 1).\n5. Specific Comments\n5.01 Line 54: The last sentence of the abstract is incomplete. Consider deleting “that.”\nResponse: we modified it in manuscript (Line 65).\n5.02 Line 152: Is ZM4_V1.0 from Long et al. 2022? Is ZM4_V1.5 the current version?\nResponse: Thanks for rising this key question. The ZM4_V1.0 is from Long et al. 2022. The\nZM4_V1.5 was the primary version, and the updated genome is ZM4_V2.0. We have added a\nfigure to illustrate our assembly pipeline (Supplementary Fig. 1).\nReference:\nLong R, Zhang F, Zhang Z, et al. Genome assembly of alfalfa cultivar zhongmu-4 and\nidentification of SNPs associated with agronomic traits[J]. Genomics, proteomics &\nbioinformatics, 2022, 20(1): 14-28.\n5.03 Line 170: How does this version improve upon the previous one? Please compare the number\nof gaps and genes between ZM4_V1.0 and ZM4_V2.0.\nResponse: thanks for your suggestion. We compared the genome size, contig length, number of\ngenes, and BUSCO results (Response Fig. 4 and Supplementary Fig. 4). The genome has shown\nsignificant improvement compared to the previous version. However, despite the contig size being\n5.45 Mb, there are still numerous gaps present, which makes it challenging to evaluate this\nparameter effectively.\nResponse Fig. 4 The genome quality comparation between ZM4_V1.0 and ZM4_V2.0. The y-axis\nrepresents the improved percentage of ZM4_V2.0.\n10\n5.04 Lines 222-225: Does this refer to the ZM4_V2.0 genome? Define unique genes clearly. How\nwas the number 109,715 determined? Provide BUSCO analysis results on single-copy and\nduplicated genes.\nResponse: Thanks for your suggestion. We rewrite the related sentence and made clear explanation\nfor the ZM4_V2.0 genome and number 109,715 (Line 286-288). The BUSCO analysis of single-\ncopy and duplicated genes of ZM4_V2.0 has been provided in Supplemental Table 1.\n5.05 Lines 225-229: How many four-copy genes are dispensable? How many single-copy genes\nare core? The statement does not confirm non-random gene retention, duplication, and loss as\nsuggested in the next sentence.\nResponse: Thanks for rising this question. There are 2,867 four-copy genes are dispensable. There\nare 7,782 genes single-copy genes are core (added detailed information in Supplementary Table\n3). We performed a chi-square test to examine gene retention patterns. The results yielded a p-\nvalue of less than 2.2e-16, indicating significant evidence of non-random patterns in gene retention,\nduplication, and loss (related contents were added in Line 293-296).\n5.06 Lines 247-248: The sentence is incomplete and unclear.\nResponse: we modified it in manuscript (Line 330-332).\n5.07 Lines 256-257: The tetra-copy core genes represented the largest proportion of what?\nResponse: The largest proportion of DEGs. We modified it in manuscript (Line 341).\n5.08 Line 283: Which 124 conditions were compared? This is not detailed in Supplementary Table\n4. Please refer to my second major concern.\nResponse: The 124 conditions were supplied in Supplementary Table 6.\n5.09 Line 330: How were gene pairs selected for analysis?\nResponse: The gene pair was selected based on the blastp results. The filter parameter was as\nfollow: e-value < 1e-10, pident > 95%, (1 - mismatch/length) > 95%, and Chr_query = Chr_match.\n(Methods, Line 700-701).\n5.10 Line 365: How was plant biomass affected?\nResponse: The biomass was increased. We added related sentence in manuscript (Line 476).\n5.11 Lines 492-493: The sentence “Additionally, haploid pseudo-chromosomes were generated\nusing the same pipeline (ZM4_V1.5)” appears redundant.\nResponse: The sentence has been removed.\n5.12 Line 539: Specify the version and parameters used for the PGGB pipeline. How was the\nVCF file obtained and filtered?\nResponse: The detailed information of PGGB has been added in manuscript (Line 676-678). The\nVCF file was generated using vg and filtered with default parameter (Line 678-679).\n5.13 Line 544: Were publicly available protein sequences used for OrthoFinder in other genomes?\nWere genes filtered for ORFs or pseudogenes? See my concern regarding inconsistent annotations.\n210-215, Table 1). The ORFs and pseudogenes were filtered by transcript support, exon and intron\ndepth, and protein length (https://github.com/futurefanzhang/genomic-adaptation-of-\nautotetraploid-alfalfa/blob/main/step1/genomic_annotation.sh). The genome and protein\nsequences will be released following the publication of our paper.\n5.14 Line 564: Provide more details on how orthologous gene families were aligned with copy\nnumber information.\nResponse: The detailed information has been added in manuscript (Line 720-727).\n5.15 Line 571: Clarify the DEG analysis methodology, including specific comparisons made.\nAddressing these concerns will strengthen the manuscript’s validity and reproducibility.\nResponse: We added more detail information in manuscript (Line 752-756). Additionally, the step-\nby-step code has been made available on GitHub (https://github.com/futurefanzhang/genomic-adaptation-of-autotetraploid-alfalfa/tree/main/step5).\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We appreciate the reviewer’s time, patience and constructive comments that help us to\nimprove this manuscript. We hope the revised manuscript is more suitable for publication.\nReviewer #4 (Remarks to the Author):\nComments to the Authors\nThe authors present a study focused on constructing a Medicago super-pangenome using newly\nassembled autotetraploid and diploid M. sativa genomes alongside published genomic data. They\nidentify tetra-copy core genes in alfalfa, explore their roles in climate adaptation and stress\nresponses, discuss their evolutionary constraints, and validate the functional relevance of one such\ngene through overexpression. While the topic is highly relevant for understanding polyploid\nadaptation and breeding applications, the manuscript has several critical limitations in data\nintegration, methodology, and interpretation that require substantial revision. Below are my major\nand minor concerns.\nimprove this manuscript. We have revised the manuscript accordingly and we hope the revised\nmanuscript is more suitable for publication.\nMajor points:\n1. The study heavily relies on published genomic data (7 species, 13 genomes), but only 5 genomes\nwere annotated using the same method. This inconsistency in annotation pipelines may introduce\nsignificant bias, particularly for identifying tetra-copy core genes, which form the foundation of\nall downstream analyses. For example, differences in gene prediction tools can inflate or obscure\northolog numbers. The authors should re-annotate all included genomes using the same pipeline\nto ensure comparability, and maybe also to address how annotation heterogeneity might have\ninfluenced the tetra-copy core gene set (e.g., sensitivity analysis). Additionally, the term\n\"comprehensive super-pangenome\" is misleading, as only 7 of ~90 Medicago species are included.\nResponse: Thanks for pointing out this important problem. We have re-annotated all remaining\ngenomes using the same pipeline. The number of genes and BUSCO assessment showed\nimprovement compared to the raw versions (Table 1). The number of tetra-copy core genes\nincreased relative to the initial version, suggesting that a uniform annotation pipeline can reduce\nbias in gene identification (Lines 660-663). Additionally, the term 'comprehensive super-\npangenome' has been revised to 'super-pangenome' (Line 497).\n2. The inclusion of distantly related Medicago species (e.g., M. truncatula and M. polymorpha;\ngenome size ~450 Mb vs. ~850 Mb in diploid alfalfa) introduces evolutionary and genomic\ndivergence that confounds tetra-copy core gene identification. These species differ in life history\n(annual vs. perennial) and mating systems (selfing vs. outcrossing), which likely influence gene\nfamily dynamics. Would it be more reasonable to restrict tetra-copy core gene analysis to closely\nrelated species (e.g., M. sativa diploids, tetraploids, and their wild relatives like M. falcata)?\nAlternatively, it would be possible to perform lineage-specific gene retention/loss analyses to\ndisentangle polyploid-specific signals from broader phylogenetic noise.\nResponse: Thank you for this insightful suggestion. We agree that focusing solely on the M.\nfalcata-sativa complex (which includes diploid M. sativa subsp. caerulea and tetraploid M. sativa\nsubsp. sativa, as well as diploid and tetraploid M. sativa subsp. falcata) could provide more\ndetailed and accurate information regarding gene retention and loss related to whole genome\nduplication. Our study has two main objectives: first, to identify which genes have been retained\nover long evolutionary periods (core genes), particularly in relation to changes in life history and\nmating systems; and second, to determine the role of these duplicated core genes in alfalfa's climate\nadaptation. Therefore, we conducted relevant analyses using Medicago species. We appreciate this\nvaluable suggestion, as it opens up a new avenue for our future studies. Accordingly, we have\nincorporated this idea into the discussion (Line 527-530).\n13\n3. The super-pangenome framework is underutilized. Beyond identifying tetra-copy genes, no\ndeeper exploration (e.g., the role of structural variation, presence-absence variation across species)\nis provided.\nResponse: Thanks for this suggestion. We have incorporated additional analyses related to SV in\nthe manuscript. These analyses include accession-specific SV, various types of SV, an assessment\nof SV bias at the genus level, and the PCA analysis based on SV (Line 264-277, Figure 2 E-H).\n4. The authors directly reuse \"climate-adaptive genes\" identified in their prior work based on\nLFMM and PBS methods applied to global alfalfa accessions. However, the alfalfa’s global\ndistribution is heavily influenced by human-mediated dispersal and selection. Consequently, how\nmany of these \" climate-adaptive genes\" identified via local climatic factors are actually driven by\nartificial selection? In their earlier study, the authors also found that nearly 100 of these \" climate-\nadaptive genes\" were attributed to introgression. Therefore, directly using this gene set to\ninvestigate its association with different gene groups in alfalfa is inappropriate, as it conflates\nnatural adaptation signals with human-driven selection and introgression effects. Therefore, the\nauthors should re-analyze the climate-adaptive signals to ensure the reliability of the conclusions.\nResponse: Thanks for rising this key question. We agree that human-mediated dispersal and\nselection could mislead the understanding of alfalfa’s global distribution. It is challenging for\nalfalfa to identify one accession that has been grown in a specific climatic region for thousands of\nyears. We believe that one accession can grow and adapt to specific climatic conditions, suggesting\nthat it must possess the relevant genes or alleles that confer this adaptability. Regarding the LFMM\nand PBS analyses, we recognize the difficulty in balancing sample size and statistical power\n(Spencer et al, 2009). To mitigate potential human influence, we focused on 240 samples that\nexhibited a clear geographic clustering pattern for the LFMM and PBS analyses. The results\ndemonstrated patterns consistent with previous findings using 397 samples (Response Fig. 5). The\nrelated contents have been modified (Figure 4A, Line 326-328).\ngroup. The red line indicates the enrichment fold.\nReference:\nSpencer C C A, Su Z, Donnelly P, et al. Designing genome-wide association studies: sample size,\npower, imputation, and the choice of genotyping chip[J]. PLoS genetics, 2009, 5(5): e1000477.\nMinor points:\nResults\n5. L184–189: Comparing total genome size (including all four subgenomes) in tetraploid alfalfa\nagainst single diploid genomes is methodologically flawed. Use haploid genome equivalents (1C\nvalues) for fair comparison.\nResponse: Thank you for your suggestion. We have added the comparison between the super\npangenome and the ZM4_V2.0 hap1 genome (Line 221). Our pangenome was constructed using\nthe four alfalfa sub-genomes along with 12 other genomes. To understand the specific contribution\nof the alfalfa sub-genomes, we compared them against the overall pangenome (Figure 1D).\n6. L205: The text mentions 11–12 species, but only 7 are included. Clarify this discrepancy.\nResponse: Thanks for your suggestion. We have modified the species to genomes (Line 251-253).\n7. L208–209: The claim that \"more homologous gene families indicate higher genome quality\"\nis unfounded. Gene family counts depend on annotation consistency and evolutionary divergence,\nnot assembly quality.\nResponse: Thanks for this suggestion. We have removed the sentence (Line 256).\n8. L245: The 1671 climate adaptation-related genes were obtained based on the union gene sets of\nLFMM and PBS analyses, not the results of GWAS.\nResponse: Thanks for this suggestion. We have modified this sentence (Line 326).\n9. L248: Define \"genomic background content\" and clarify how the 19.4% was calculated.\nResponse: we added more detailed explanation and defined “genomic background content” in the\nmanuscript (Line 310-316). The percentage has been revised from 19.4% to 21.5% due to re-\nannotated gene data.\n10. L270: What is \"normal condition\"?\nResponse: the normal condition means the condition before treatment (added explanation in the\nmanuscript Line 362-363).\n11. L339–341: Figure 4F compares two species but refers to \"highest expression\"—this is\nambiguous.\nto “highest deleterious variants accumulation” in the initial version. We have modified the Figure\nlegend (Line 450).\n12. L348: Should be \"genomic\".\nResponse: The word has been revised (Line 457).\nMethods\n13. L528–529: Justify why uniform annotation was not applied to all species.\nResponse: The uniform re-annotation has been done for all genomes (Line 210-215).\n14. L531: Clarify whether repeat annotation was performed before or after pangenome\nconstruction.\nResponse: The de novo annotation was conducted for each genome before pangenome\nconstruction (Line 660-663, and Line 675-678).\n15. L548: Detailed information of divergence time estimation (e.g., software, calibration points)\nare needed.\nResponse: The software detailed information and calibration points was supplied in the manuscript\n(Line 691-694).\n16. L551: The authors should clearly distinguish \"duplicated genes\" (paralogs) from \"allelic genes\"\n(homologs within a polyploid).\nResponse: The paralogs were filtered out during the filtering of the blastp results. Related details\nhave been added to the manuscript (Line 710-711).\n17. L581–590: I am completely confused about the genetic burden analysis. For example, how is\nthe phylogenetic tree constructed related to the GERP scores? How are the GERP scores calculated?\nAre GERP scores calculated for all 15 species, or only for a subset of them?\nResponse: The initial phylogenetic tree was constructed based on single-copy genes. Subsequently,\na neutral tree was created by combining four-fold sites (those with A/T/C/G genotypes across the\n15 genomes). The GERP score was calculated using the multiple genome alignment file of the 15\nspecies in conjunction with the neutral tree. To facilitate understanding, we have provided detailed\ncode on GitHub (https://github.com/futurefanzhang/genomic-adaptation-of-autotetraploid-\nalfalfa/blob/main/step4/gerp_score.sh). (Line 763-773, A more detailed explanation about GEPR\nscore was provided in response to question 6 from Reviewer #1.)\n18. Figure 4F. The diploid alfalfa genome should have two haplotypes. How are \"three-copy\"\nand \"four-copy\" genes defined here? Provide a methodological explanation.\nResponse: This part was analysis based on population resequencing data. The raw data of diploid\nalfalfa was mapping to the ZM4_V2.0 hap1 genome. The \"three-copy\" and \"four-copy\" gene was\n16\nidentified based on the reference of ZM4_V2.0 hap1 genome (related contents has been added in\nFigure 5F and line 731-734).\nReviewer #5 (Remarks to the Author):\nPolyploid plants often outperform their diploid sisters, but when they do, the genetic and\nmechanistic bases for their advantages remain unclear. In this study, researchers generated a\nhaplotype-resolved autotetraploid alfalfa genome and incorporated it into a super-pangenome\nspanning 13 genomes from seven Medicago species. They report dramatic gene content variation,\nwith only about 20% of genes consistently present across all four haplotypes, a striking and rather\nsurprising finding. A key subset, the tetra-copy core genes, appeared important for climate\nadaptation and stress responses, though they also carried a heavier genetic burden. While the study\nsuggests that polyploidy can enhance traits like biomass and nitrogen use efficiency, some aspects\nrequire better evidence precisey due their importants. In particular, the quality of the tetraploid\nassembly and its true biological significance remain open questions that require further validation.\nResponse: We are grateful for the reviewer’s constructive comments. We have carefully revised\nthe manuscript accordingly based on all these suggestions. We hope the revised manuscript is more\nsuitable for publication.\nMajor Concern:\n1. My single major concern about this paper is that autotetraploid assembly is not sufficiently\ndescribed in terms of quality. What I see as the major finding of the paper hangs completely on\nthis, so this is important.\nResponse: Thank you for highlighting this important point. We compared several key parameters\nto illustrate the quality of the genome. Firstly, the estimated genome size of alfalfa was reported\nto be between 3.0 and 3.3 Gb based on previous studies (Long et al., 2022; Shen et al., 2020; Chen\net al., 2020). The assembled genome size of ZM4_V2.0 was 3.13 Gb. Furthermore, we mapped\nthe raw HiFi and ONT data, achieving mapping rates of 99.93% and 99.92%, respectively. The\nBUSCO analysis of the genome yielded a completeness score of 99.5%. Finally, using Merqury\nsoftware, we assessed the base quality and genome completeness, resulting in scores of 61.88 and\n97.78%, respectively (Line 180-185). The quality of the genome is comparable to that of\nautotetraploid potato (Sun et al., 2025). Therefore, we believe the genome quality of ZM4_V2.0\nis high enough to serve as a reference genome for alfalfa and to be used in the analysis of\nautopolyploid.\nReference:\nLong R, Zhang F, Zhang Z, et al. Genome assembly of alfalfa cultivar zhongmu-4 and\nidentification of SNPs associated with agronomic traits[J]. Genomics, proteomics &\nbioinformatics, 2022, 20(1): 14-28.\nand resequencing of core germplasms provide genomic resources for alfalfa research[J]. Molecular\nPlant, 2020, 13(9): 1250-1261.\nChen H, Zeng Y, Yang Y, et al. Allele-aware chromosome-level genome assembly and efficient\ntransgene-free genome editing for the autotetraploid cultivated alfalfa[J]. Nature communications,\n2020, 11(1): 2494.\nSun H, Tusso S, Dent C I, et al. The phased pan-genome of tetraploid European potato[J]. Nature,\n2025: 1-9.\n2. My concern as someone is that Hifiasm is not built for polyploid (especially autopolyploid)\ngenome assembly: the --n-hap 4 tag is still rather experimental (esp as of v0.18). It is not that I\ndoubt strongly the results, but better documentation (meaning detailed metrics into the\nsupplemental information) describing the process and quality of the assembly is required. Such\ndocumentation of correct assembly is crucial for us to believe the key finding figure 3A and\ntherefore this is a major concern.\nResponse: Thank you for rising this important question. We have included a detailed assembly\npipeline in Supplemental Fig. 1. The haplotype-resolved reads were separated based on Hi-C\ninformation and ZM4_V1.0 (Supplemental Figs. 2-3). The parameter '--n-hap 4' was applied using\nhifiasm V0.19.7 (we apologize for any lack of clarity). Additionally, the detailed steps of the\nassembly process have been uploaded to GitHub (https://github.com/futurefanzhang/genomic-adaptation-of-autotetraploid-alfalfa/blob/main/step1/genome_assembly.sh).\n3. The core issue is this: the finding that half the genes are present only in one haplome is highly\nstriking and thus begs further validation. The contig N50 is not great, but probably not so poor to\nproduce such serious widespread artifacts to invalidate the key finding. That said, greater detail of\nthe assembly process and statistics is needed.\nResponse: We apologize for the lack of clarity. We define a gene as having multiple allelic copies\nif they exhibit at least 95% similarity, meaning that only highly conserved genes are considered as\nallelic copies. Furthermore, the 45.0% genes belong to one-copy genes was calculated by unigenes\n(where multiple copies are counted only once). Among the total genes, the proportion of one-copy\ngenes is 24.4%. We have added more detailed information in the manuscript (Line 284-290, and\nLine 310-314, Supplemental table 3).\n4. In addition, confidence in the finding can also be productively provided by showing gene\nneighbourhood-level examples of widespread presence/absence variation as well. Ideally this\nwould be shown for several (or better, 8 to 12) regions e.g. in the SI, but at least one large region\nin the main text figures, demonstrating clearly that presence/absence variation of genes among the\nfour haplomes. To drive this home, a track for each haplome of hifi read mapping depth and\ngaplessness at least in the regions shown would solidify the finding.\nResponse: Thank you for your valuable suggestion. We have added one region (chr1: 8.4-8.8 Mb)\nto illustrate the allelic copy number variation (Figure 3B). Regarding mapping depth and gapless\nassembly, we would like to clarify that we did not identify genes based on SV presence or absence;\ninstead, we identified duplicate genes based on protein conservation. The key distinction is that\nSV-based PAV is more likely to identify newly arisen/loss genes, whereas protein conservation\nprovides possibility to identify allelic copy number variation (We have revised the manuscript to\ninclude further details that clarify this concept, Line 282-285, 711-715). Nevertheless, your\nsuggestion opens up a new direction for us to further explore the polyploid adaptation of alfalfa.\n5. Finally, relating also to this concern: the wider significance of this finding of what I would term\n‘autopolyploid fractionation’ would be very helpful for the reader. Can we expect this to be a\ncommon phenomenon, or is this line-specific? The lack of information about the Zhongmu-4 line\nused is not helpful in this regard. E.g. is it highly heterozygous and/or does it have a history of\ninbreeding? Could this finding be partly explained by historical inbreeding? Ideally Tajima’s D if\nthe data can support this.\nResponse: Thanks for this insightful suggestion. We believe that autopolyploid fractionation is a\ncommon phenomenon. This has been observed in autopolyploid potatoes, where 54% of genes are\npresent across all four haplotypes (Sun et al., 2022), as well as in sugarcane, where 10.3% of genes\nexhibit four alleles (Bao et al., 2024). Therefore, we consider that autopolyploid fractionation is\nlikely a prevalent occurrence in autopolyploid species.\nRegarding the exploration of line-specific effects in Zhongmu-4, we believe that the phased\nhaplotypes of the pangenome in the M. sativa-falcata complex will be very beneficial. Performing\nTajima’s D analysis using haplotype-resolved variants should yield valuable insights into the\neffects of heterozygosity and inbreeding. It will be a significant direction for our future research\n(added discussion in Line 527-530).\nReference:\nSun H, Jiao W B, Krause K, et al. Chromosome-scale and haplotype-resolved genome assembly\nof a tetraploid potato cultivar[J]. Nature genetics, 2022, 54(3): 342-348.\nBao Y, Zhang Q, Huang J, et al. A chromosomal-scale genome assembly of modern cultivated\nhybrid sugarcane provides insights into origination and evolution[J]. Nature Communications,\n2024, 15(1): 3041.\n6. Other than this single major concern above, my remarks below can be classified as minor, but\nalso hopefully helpful to improve this paper.\nResponse: We sincerely appreciate the reviewer’s time, patience, and constructive feedback. The\nsuggestions have proven invaluable to us and have significantly enhanced the quality of the\nmanuscript.\n7. I’m unsure the title makes sense: in what way is the superpangenome reveal “adaptive\nadvantages?” What is shown is very interesting but really what is shown are correlates of potential\nadaptive nature, not ‘revealing adaptive advantages’. Such language is overblown. The\n‘evolutionary constraints’ also are less clear. Constraint is usually thought of in a different way:\ncan this be better phrased?\nResponse: We thank the reviewer for this suggestion. Our aim is to uncover the adaptive\nadvantages and evolutionary constraints following the WGD in alfalfa. The title reflects potential\ncomparisons made before and after WGD. Furthermore, similar titles have been used in related\nstudies on polyploidy (Comai et al., 2005; Weiss et al., 2013 and Fawcett et al., 2010).\nReference:\nComai L. The advantages and disadvantages of being polyploid[J]. Nature reviews genetics, 2005,\n6(11): 836-846.\nWeiss-Schneeweiss H, Emadzade K, Jang T S, et al. Evolutionary consequences, constraints and\npotential of polyploidy in plants[J]. Cytogenetic and genome research, 2013, 140(2-4): 137-150.\nFawcett J A, Van de Peer Y. Angiosperm polyploids and their road to evolutionary success[J].\nTrends in Evolutionary Biology, 2010, 2(1): e3-e3.\n8.The introduction is very good, but with some room for improvement. I suggest at line 76 and\nfollowing that the presentation is overly adaptationist and should cite the dark side of polyploidy.\nFor example Monnahan et al (https://doi.org/10.1038/s41559-019-0807-4) describe a decrease in\nthe efficiency of purifying selection and increased increase in load. Mention of this would give a\nmore balanced introduction and presage your interesting findings on burden.\nResponse: Thanks for this suggestion. We have added the dark side of polyploidy and cited related\npaper (Line 98-103).\n9. Lines 146 and following: I think it would be more effective and intuitive to state the sequencing\ndepth versus the haploid genome size, rather than the four haplotypes together. It took several\npasses reading this to understand you meant per haplotype (an unusual convention); at first one\nhas the impression that, for example, the nanopore sequencing is only six-fold total, rather than\nsix-fold versus each of the four alleles. That is, you risk misunderstanding for the reader and\naccidental underselling the excellent data the way it's currently presented.\nResponse: Thanks for giving us this suggestion. We have revised our expression strategy (Line\n166-168).\n10. Line 216, the retention of core genes enhanced climate adaptation in alfalfa subheading is\nindicating a circumstantial finding. It is difficult to make the statement that the retention of the\ncore genes actually did enhance climate adaptation. I would say is correlated with climate\nadaptation rather than enhanced it.”\nResponse: We have modified the subheading as your suggestion (Line 279).\n11. Line 231 “fully functional”: Nothing here is demonstrated to be ‘fully functional’; this is an\ninappropriate statement. Indeed, functional in line 233 and unique essential both seem overstated\nand have not been subject causal testing, only correlation.\nPlease back off on this language. It is unnecessary.\nResponse: We have modified the related expression contents (Line 302-304).\n12. Line 249, here, for example, it is fair, the way it is stated, the finding ‘suggested’, that tetra-\ncopy core genes are essential.\nResponse: We have modified the related sentence (Line 332-333).\n13. Lines 262 to 264, the first half is a fragment, join sentences but also simplify.\nResponse: We have modified the related sentence (Line 348-349).\n14. Line 266: “was gradually decreased”: better to say something like, ‘exhibits a quantitative\ngradation” or so.\nResponse: We have modified the related sentence (Line 357-359).\n15. Line 284: Supplemental Table 4 is simply a list of 23 references. The main text says “we\nanalyzed a total of 124 different conditions”. Please spell these out in the table and make it more\nclear for the reader what you are talking about.\nResponse: We have added one table to explain 124 conditions (Supplementary Table 6).\n16. Line 292: “hap2 and hap3 was 25%”. This is completely unclear. Please better show the metrics\nin question and clarify in the text what this is.\nResponse: We have removed this expression (Line 388).\n17. Line 318: this on the following line are very unclear to me: can you please rephrase and\nelaborate this in the text?\nResponse: We apologize for the unclear description. We have modified the related sentence (Line\n416-417).\n18. Line 326: to define genes that accumulate more delirious variants than the genomic average as\n‘deleterious genes’ seems extreme.\nResponse: We agree that define deleterious genes is hard. Our objective is to investigate the\nevolutionary properties of genes that accumulate a disproportionately high load of deleterious\nvariants. Therefore, we believe that utilizing this threshold is a reasonable approach.\n19. Line 385-386: you do not confirm that tetra-copy core genes ‘are essential for enhancing\npolypoid adaptability’. Please tone back the language: you have interesting findings but you don’t\nactually confirm either ‘enhancing adaptability’ or demonstrate that those genes are ‘essential’ for\nthat in particular. The language is overblown unnecessarily.\nResponse: Thanks for this suggestion. We have modified this sentence (Line 499-500).\n20. Line 394: ‘past one century’: remove ‘one’\nResponse: We have modified this sentence (Line 509).\n21. Line 412: this does not seem the case at all! You show a clear gradient and not expression so\ndominated by single haplotype. Is this calling to Figure 3G or am I missing a callout to a panel in\nthe figures? If so, it should be pointed to even in the discussion.\nResponse: We apologize for the inappropriate description. We have modified this sentence (Line\n532-535).\n22. Line 419: this is very interesting this part of the discussion.\nResponse: Thanks for your commons.\n23. Line 437: ‘TCP1’ is incorrect: it is ‘TPC1’ one. Additionally that paper does not show\n‘heterozygote advantage’. It would be more effective to use that paper and others on autopolyploid\npopulation genomics, such as the Hamala or Monnahan (mentioned above) to indicate that\nincreased load also allows masking of potentially useful alleles for adaptation to, for example here,\nserpentine.\nResponse: Thanks for your suggestion. We have modified this sentence (Line 560-561). Since we\nhave utilized the paper you mentioned to illustrate the effects of increased load on long-term fate\n(Line 100-101), we want to use this paper to illustrate the relationship between deleterious\nmutations and adaptation.\n24. Line 447: this distribution pattern and supplemental figure 8 is completely opaque. Please it\nwould be better to show the relative enrichment of tetra-copy core genes versus all genes in across\nthe genome in supplemental figure eight: not just the pattern of the tetra-copy core genes. I expect\nthis pattern that you show closely reflects all genes not just tetra-copy core genes. That is, of course\ntetra-copy core genes very crudely as shown tend to cluster at the ends of chromosomes relative\nto the centromeres: this is because all genes do! What would help your case is to show the relative\ndensity of tetra-copy genes, normalized vs overall gene density in the figure. This would then allow\nthe hypothesis to be tested that core genes are, in fact, differently distributed from the rest of the\ngenes (e.g. harder to purge or so if this is what your hypothesis).\nResponse: Thanks for your valuable suggestion. We have redrawn this figure and modified the\nsentence (Line 420-423 and Supplemental Figure 17). The enrichment of tetra-copy core genes\n(the density of tetra-copy core genes / the density of all genes) indicates that tetra-copy core genes\nare preferentially located at both ends of the chromosomes and are situated in similar positions\nacross the corresponding haplotypes (Response Fig. 6).\nrepresents the enrichment fold of tetra-copy core genes compared to all genes present on the\nchromosome\n25. Line 463: here again it would be very good to relate your findings to Monnahan and Hamala\nMSs, indicating that this paper is congruent with findings that load is increased in autopolyploids,\nand indeed that purifying selection is perhaps less efficacious.\nResponse: Thanks for your suggestion. We have added related contents in discussion (Line 587-\n589).\n26. Line 470: this is very interesting.\nResponse: Thanks for your commons.\n27. In table one, the BUSCO scores are incomplete. I am intetested to see the C, D and S metrics.\nResponse: we have added it in Supplemental Table1.\n28. Overall, this work provides exciting insights, highlighting extensive gene presence/absence\nvariation and the potential role of tetra-copy core genes in stress adaptation. However, the\nreliability of the key findings hinges on the quality of the tetraploid genome assembly, which\nrequires more detailed documentation and validation. Addressing this, along with refining\noverstatements in the text and improving clarity in data presentation, will strengthen the study's\nimpact.\nResponse: We sincerely appreciate the reviewer’s time, patience, and constructive feedback, which\nhave greatly contributed to improving this manuscript. We have revised the manuscript in\naccordance with your valuable suggestions.\nContents\nPart I: Summary of responses to the major questions: ............................................................. 2\nPart II: Point by point response to specific comments: ............................................................. 3\nReviewer #1 (Remarks to the Author): ....................................................................................... 3\nReviewer #2 (Remarks to the Author): ....................................................................................... 4\nReviewer #3 (Remarks to the Author): ....................................................................................... 6\nReviewer #4 (Remarks to the Author): ....................................................................................... 6\nReviewer #5 (Remarks to the Author): ....................................................................................... 6\nWe sincerely appreciate the constructive feedback and valuable suggestions provided by the editor\nand the five reviewers. Before presenting a detailed point-by-point response in Part II, we first\naddress two major points:\nFirst, we have validated the high accuracy of our four-haplotype phasing. This is supported by a\nlow switch error rate (1.97%) and by the fact that over 99% of high-confidence mapped HiFi reads\naligned to a single, corresponding haplotype (83.02% total mapping rate, MAPQ ≥ 10). These\nfindings suggest that the majority of genome regions were successfully phased into four haplotypes.\nSecond, we have carefully revised the manuscript to improve its clarity. Specifically, we have\nrevised the term 'genomic background' to 'genomic baseline proportion', and clarified the term\nwherever else it is used.\nIn summary, we have provided detailed explanations in response to the reviewers' comments. All\nmodifications in the manuscript are highlighted with tracked changes. Additionally, we have\nrevised the manuscript, supplemental figures and tables.\n(The reviewer’s comments are shown in black, and the authors’ responses are displayed in blue.)\nReviewer #1 (Remarks to the Author):\nThank you for the revision. The updated version has been significantly improved. Most of my\nconcerns have been addressed, except for the following:\nFig. 2E is still a concern. The SVs of species following R108 show approximately tenfold\ndifferences compared to ZM4/XJDY/Mca. The small number of SVs across species might be due\nto sequence similarity. I suggest surveying sequence similarity across these species. This should\nnot conflict with the gene-based synteny analysis (Fig. 1A), as most genes should show their\npresence.\nResponse: We thank the reviewer for this valuable suggestion. To assess the sequence similarity\nprior to SV calling, we calculated the mapped reads of each species. We found that a species from\noutside the M. sativa-falcata complex produced only ~35% of the number of mapped reads\ncompared to a subspecies within the complex (Response Fig. 1). Furthermore, this divergence in\nSVs was magnified after applying standard filters for mapping quality, sequence length, and\nvariant types. This suggests that robust cross-species SV calling at the species level remains a\nsignificant challenge (Lines 277-280, Supplementary Fig. 12).\nResponse Fig. 1 Relationship between mapped reads and the number of SVs, based on\nalignments generated by (a) mummer4 and (b) minimap2.\nGeneral Remarks: Thank you to the authors for providing clarification on my previous concerns.\nWith the major issues regarding uniform annotation and clarification of the DEG analyses now\nresolved, I have reread the manuscript with a fresh perspective. The overall structure and flow is\nimproved. However, there are still some minor issues that need to be resolved before final\nacceptance:\nResponse: We thank the reviewer for the positive feedback and helpful suggestions. We have\nrevised the manuscript accordingly and believe it is now much improved and suitable for\npublication.\nLine 47: Please spell out “DEGs” (differentially expressed genes) the first time it appears in the\nmanuscript.\nResponse: Thanks for your suggestion. We have revised the sentence (Line 43).\nLine 164: Change “high-yield and salt-tolerance” to “high-yielding and salt-tolerant”\nResponse: We have revised the sentence (Line 164).\nLine 181: Indicate which database was used for the BUSCO analysis (e.g.,embryophyta_odb10?).\nResponse: We have added the embryophyta_odb10 database information (Line 192).\nLine 198: Please provide a citation for the EggNOG database.\nResponse: We thank the reviewer for pointing out this omission. We apologize for the oversight\nand have now added the citation (Line 207).\nLine 286-288: If I understand correctly, the term unigenes here refers to a non-redundant set of\nannotated genes, where each gene is counted only once, regardless of allelic variation or copy\nnumber. I recommend using a clearer term, such as “non-redundant genes” instead of unigene. The\nterm unigene has a specific historical meaning in biology, particularly referring to the NCBI\nUniGene database, which clustered transcripts from the same gene. Using it here could mislead\nreaders unless clearly defined.\nResponse: Thanks for your kindly suggestion. We have replaced all “unigenes” to “non-redundant\ngenes” in the manuscript.\nLine 315: Similarly, I suggest replacing “genomic background” with a clearer term such as\n“background gene set” or specifying whether it refers to expression-based background or copy-\ngroup background. Clarification is needed.\nResponse: Thanks for your suggestion. We have used “genomic baseline proportion” for the\ncomparation in DEGs analysis and climate adaption candidate gene analysis (Lines 327-330), and\nclarified the term wherever else it is used.\nLine 330: In the sentence starting with “Compared to its genomic background content,” what does\n“its” refer to? The focal gene group? The whole genome? Please clarify what is being enriched\n4\n(climate adaptation-associated genes?) and against what baseline. The sentence is currently\nambiguous.\nResponse: We thank the reviewer for highlighting this point. We agree that the original sentence\nwas unclear and have revised it for greater clarity (Lines 338-342).\nLine 365: “Higher percentage” of what? Is this referring to differentially expressed genes? Please\nspecify.\nResponse: It referred to top expressed genes (from top 50% to top 1%). We have revised the related\nsentence (Lines 380-386).\nLine 367: The term “genomic background” is used somewhat loosely and becomes increasingly\nunclear later in the manuscript. It’s not obvious whether it refers to the entire genome or a group-\nspecific background set. If I understand correctly, this term is used in two contexts:\n1. To describe gene enrichment among six groups in relation to climate adaptation.\n2. To describe DEG comparisons from RNA-seq data.\nPlease clearly define genomic background in each context to prevent confusion.\nResponse: We thank the reviewer for highlighting this important point. To improve clarity and\nconsistency, we have revised the term 'genomic background' to 'genomic baseline proportion'\nthroughout the manuscript. This term is defined as the proportion of a given gene group within the\nentire genome and serves as the reference for comparison in both the climate adaptation-associated\ngenes and DEGs analyses (Lines 327-330, Figure 3E).\nLine 377-379: How many total genes (from all conditions) were present in the top 10% expression\nquantile? Please specify.\nResponse: The total genes present in the top 10% expression quantile ranged from 2,106 to 6,225\nin 124 expressed conditions (Lines 401-402, Supplemental Table 9).\nLine 395: Genetic burden section is very interesting.\nResponse: We thank the reviewer for this encouraging comment.\nLine 735-736: Define the acronyms LFMM (Latent Factor Mixed Model?) and PBS (Population\nBranch Statistic?) upon first use in the manuscript.\nResponse: The acronyms have now been defined in the manuscript (Lines 336-337).\nAlso, earlier in the Results section (Lines 318–333), there is a discussion of a custom enrichment\nanalysis for climate adaptation–related genes across six gene groups. Please provide a clear\nmethodology for this analysis in the Methods section e.g. what was considered the observed versus\nexpected values in the chi-squared test?\nResponse: We have now provided a detailed description of the chi-squared test in the Methods\nsection (Lines 853-863). Furthermore, to ensure reproducibility, the code implementing this\nanalysis is publicly available on GitHub (https://github.com/futurefanzhang/genomic-adaptation-\nof-autotetraploid-alfalfa/blob/main/step5/p_value_and_chi_squared.R).\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We thank the reviewer for the valuable time and effort. We have carefully addressed all\nthe points raised, and believe these revisions have significantly strengthened the manuscript. We\nbelieve the manuscript has been substantially improved with these revisions and hope it is now\nsuitable for publication.\nReviewer #4 (Remarks to the Author):\nAll major concerns have been effectively addressed. However, one minor issue remains\nunresolved- the explanation of genomic background content is still unclear (Point 9). This term\nappears to be used in two different contexts (climate adapation gene retention and DEG analyses),\nplease clarify whether it carries the same meaning in both cases or if it refers to distinct concepts.\nA clear definition in each context would improve the manuscript’s clarity.\nResponse: We thank the reviewer for this helpful suggestion. To improve clarity, we have now\nstandardized the terminology and consistently use 'genomic baseline proportion' throughout the\nmanuscript for both the adaptation-associated genes and DEGs. This makes it explicit that the same\nreference dataset was used for both analyses. A clear definition has been added to the manuscript\n(Lines 327-330, Figure 3E).\nReviewer #5 (Remarks to the Author):\nOverall this revision is very much improved and a pleasure to read. Broadly, the authors have made\nsome progress in alleviating my major concerns in this revision.\nThat said, I still do have a substantial concern about the autotetraploid assembly that goes beyond\nregular QC metrics. I appreciate very much that the authors present a much clearer exposition of\nthe assembly strategy and QC now, which is very helpful. I see that in terms of completeness and\nbase accuracy, as well as mapping, quality is high. But this is only normal nowadays and does not\nspeak to the key point: the accurate representation of all four haplotypes. The authors provide more\ncomplete BUSCOs (Supp table 1), where D busco scores are 10-12% in each haplome. It seems\nvery much in conflict with the overall pattern of widespread allele fractionation across the genome.\nCan the authors resolve this appearance, or provide some convincing explanation here, based on\nanalysis? In our work, we see very much inflated D buscos in autotetraploid long read assemblies\nusing hifiasm, which are hard to reduce, and almost certainly artifactual due to the fact that the\nassembler is simply not designed for autopolyploids. What can the authors say about these high D\nscores in the face of widespread (and dramatic) signal of autopolyploid *fractionation* in their\nbroader analysis? I note that while the authors rebut that their autotetraploid assembly is of similar\nquality to the Sun et al potato, the latter have only 2% D scores.\nResponse: We thank the reviewer for the insightful question, which touches upon the core\nchallenge of assembling a complex autotetraploid genome.\nAs for the accurate representation of all four haplotypes, we calculated the switch error\nand obtained a rate of 1.97%. Furthermore, we remapped the raw reads back to our phased\nhaplotypes. Among the high-confidence mapped HiFi reads (83.02% mapping rate, MAPQ ≥\n10), over 99% align to a single, corresponding haplotype (Response Fig. 2). This demonstrates that\nfor the majority of the resolvable genome, which can also be considered fractionation regions,\nour phasing framework is highly accurate and the four haplotypes are clearly separated (Lines 184-\n190, Supplementary Fig. 5).\nWe suggest that the inflated D-score is a realistic reflection of alfalfa's complex genome and\ninherent limitations of current technology (sequencing and assembly). We got this conclusion\nbased on two key observations. First, within the context of our 19-haplotype pangenome, the D-\nscore of ZM4_V2.0 is not an outlier and follows a positive correlation with genome size,\nshowing it is comparable to assemblies of similar scale (Response Fig. 3). Second, the genome contains\nextensive \"difficult-to-resolve\" regions, a fact directly quantified by our read mapping results and\nHamming error of 7.12%. Despite their >99.9% accuracy, a substantial portion (~17%) of HiFi\nreads failed to map uniquely (83.02% mapping rate, MAPQ ≥ 10). The 14% of cross-haplotype\nmappings observed with ONT reads (88.34% mapping rate, MAPQ≥10) further underscores the\nchallenge of analyzing this genome, as their intrinsic error rate makes it difficult to distinguish\nbetween highly similar haplotypes. These complex regions not only explain the assembly\nchallenges but also directly contribute to the inflated D-score.\nIn comparison to potato, we propose that alfalfa's outcrossing reproductive strategy leads to more\nfrequent and complex recombination among its four haplotypes. This results in extensive\nhomologous regions with high sequence similarity, which likely inflates the D score.\nIn summary, autopolyploid fractionation and inflated D score represent different regions of the\ngenome. The majority of autopolyploid fractionation region was successful phased to four\nhaplotypes, while inflated D score arises primarily from the \"difficult-to-resolve\" regions where\nphasing is intrinsically ambiguous. The high rate of unmappable HiFi reads quantifies the large\nscale of these regions, while the ambiguous ONT mappings illustrate the mechanism by which\nphasing errors can occur within them.\ncolors are used to distinguish among three read mapping categories: (i) single-haplotype mapping,\nwhere a read maps uniquely to one haplotype; (ii) intra-haplotype multi-mapping, where a read\nmaps to multiple locations within the same haplotype; and (iii) cross-haplotype mapping, where a\nread maps to different haplotypes.\nResponse Fig. 3 The relationship between genome size and duplicated BUSCOs score.\n164 'tolerance' to 'tolerant'\nResponse: we have corrected the word from 'tolerance' to 'tolerant' (Line 164).\n282 'essential' is overplayed language and not something one can test this way. I suggest not\noverplaying it.\nResponse: we have revised the word from 'essential' to 'major' (Line 297).\n295 suggest change 'significant' to 'substantial'\nResponse: we have revised the word from ' significant' to ' substantial' (Line 309).\n360 language like 'crucial' risks unnecessarily overplaying results\nResponse: we have corrected the word from ' crucial' to 'important' (Line 375).\n504 'autotetraploid' to 'autotetraploids'\nResponse: we have revised the word from ' autotetraploid' to ' autotetraploids' (Line 545).\n507 'genes' to 'gene'; delete 'the'\nResponse: we have revised the sentence (Line 549).\n545 'one haplotype gene': rephrase?\nResponse: we have revised the sentence (Lines 596-599).\n562-3 'railways': more accurate to simply say 'adaptive early flowering' or so, to clearly specify\nthe adaptation itself: rapid cycling in dry environments.\nResponse: Thanks for this suggestion. We have revised the sentence (Line 624).\nAgain, I appreciate the effort put into this improved revision.\nResponse: We are very grateful to the reviewer’s time, effort, and constructive suggestions, which\nhave helped us significantly improve the manuscript. We have carefully addressed all the points\nraised and trust that the revised manuscript is now suitable for publication in Nature\nCommunications.\nReviewer #1 (Remarks to the Author):\nThank you for your update.\nCross-checking the correlations between mapping reads and the number of SVs helped clarify the\nvariation for SVs (Fig. 2E and Supplementary Fig. 12). This has explained my concern clearly and\naddressed my question from the previous version.\nResponse: We thank the reviewer for the positive comments and support of our manuscript.\nReviewer #2 (Remarks to the Author):\nThe revised version is acceptable.\nResponse: We are grateful for the reviewer's positive assessment and support for publication.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We thank the reviewer for the constructive feedback. Having addressed all points raised,\nwe believe the revised manuscript now meets the standards for publication in Nature\nCommunications.\nReviewer #5 (Remarks to the Author):\nThe revision and replies have satisfied my queries.\nI congratulate the authors on a very well-presented manuscript that will be of broad interest.\nResponse: We are very grateful for the reviewer's positive feedback and constructive comments.\nAfter addressing the concerns raised, we believe the manuscript is now suitable for publication."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thank you for your update.\nCross-checking the correlations between mapping reads and the number of SVs helped clarify the variation for SVs (Fig. 2E and Supplementary Fig. 12). This has explained my concern clearly and addressed my question from the previous version.", "reviewer_2": "The revised version is acceptable", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "The revision and replies have satisfied my queries.\nI congratulate the authors on a very well-presented manuscript that will be of broad interest."}, "author_response": "Response to Review Comments\nReviewer #1 (Remarks to the Author):\nThank you for your update.\nCross-checking the correlations between mapping reads and the number of SVs helped clarify the\nvariation for SVs (Fig. 2E and Supplementary Fig. 12). This has explained my concern clearly and\naddressed my question from the previous version.\nResponse: We thank the reviewer for the positive comments and support of our manuscript.\nReviewer #2 (Remarks to the Author):\nThe revised version is acceptable.\nResponse: We are grateful for the reviewer's positive assessment and support for publication.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We thank the reviewer for the constructive feedback. Having addressed all points raised,\nwe believe the revised manuscript now meets the standards for publication in Nature\nCommunications.\nReviewer #5 (Remarks to the Author):\nThe revision and replies have satisfied my queries.\nI congratulate the authors on a very well-presented manuscript that will be of broad interest.\nResponse: We are very grateful for the reviewer's positive feedback and constructive comments.\nAfter addressing the concerns raised, we believe the manuscript is now suitable for publication."}]

Decision letter