Refining polyploid breeding in sweet potato through allele dosage enhancement

Refining polyploid breeding in sweetpotato through allele dosage enhancement Xiangbo Zhang1,9, Chaochen Tang1,9, Bingzhi Jiang1,9, Rong Zhang1,9, Ming Li2,3,9, Yaoyao Wu4, Zhufang Yao1, Lifei Huang1, Zhongxia Luo1, Hongda Zou1, Yiling Yang1, Minyi Wu5, Ao Chen5, Shan Wu6, Xingliang Hou5, Xu Liu5,*, Zhangjun Fei6,7,*, Junjie Fu8,*, Zhangying Wang1,* 1Crops Research Institute, Guangdong Academy of Agricultural Sciences & Key Laboratory of Crop Genetic Improvement of Guangdong Province, Guangzhou, China 2College of life Sciences, Chongqing Normal University, Chongqing, China 3Institute of Biotechnology and Nuclear Technology, Sichuan Academy of Agricultural Sciences, Chengdu, Sichuan, China 4College of Horticulture, Nanjing Agricultural University, Nanjing, China 5Guangdong Provincial Key Laboratory of Applied Botany, and State Key Laboratory of Plant Diversity and Specialty Crops, South China Botanical Garden, Chinese Academy of Sciences, Guangzhou, China 6Boyce Thompson Institute, Cornell University, Ithaca, NY, USA 7USDA-ARS, Robert W. Holley Center for Agriculture and Health, Ithaca, NY, USA 8State Key Laboratory of Crop Gene Resources and Breeding, Institute of Crop Sciences, Chinese Academy of Agricultural Sciences, Beijing, China 9These authors contributed equally *Email: [email protected]; [email protected]; [email protected]; [email protected] Abstract Allele dosage plays a key role in phenotypic variation of polyploids. Here we present a genomewide variation map of hexaploid sweetpotato that captures allele dosage information, constructed from deep sequencing of 294 hexaploid accessions. Genome-wide association studies identified quantitative trait loci with dosage effects on 23 agronomic traits. Our analyses reveal that sweetpotato breeding has progressively increased the dosage of favorable alleles to enhance trait performance. Notably, the Mesoamerican gene pool has evolved toward higher dosages of favorable alleles at multiple loci, which have been increasingly introgressed into modern Chinese cultivars. We substantiated the breeding-driven dosage accumulation through transgenic validation of IbEXPA4, an expansin gene influencing tuberous root weight. Additionally, we explored causative sequence variations that alter the expression of the Orange gene, which regulates flesh color. Our findings illuminate the breeding history of sweetpotato and establish a foundation for leveraging allele dosages in polyploid breeding practices. Introduction Approximately 40% of cultivated plant species are polyploid1, a characteristic driven by wholegenome duplications (WGDs) and central to the origin and evolution of most crop species2-5. Compared to their diploid counterparts, polyploids typically exhibit enhanced growth vigor and a broader geographic distribution, reflecting the increased phenotypic variation enabled by multiple alleles in their homoeologous chromosomes6. Genome-wide gene expression studies in polyploids have revealed both additive and nonadditive dosage effects on gene expression7-9. For instance, the proportion of expressed genes with additive effect ranges from 65% to 95% in Arabidopsis allotetraploids, while that with nonadditive effect varies from 5% to 38%10,11. In hexaploid Camelina sativa, the additive effect of allele dosage at the FAD2 loci on seed oleic acid profile has been demonstrated in CRISPR-Cas9 edited lines12. To date, dosage effects have primarily been studied at the gene expression level and on several cloned genes. However, the connection between allele dosage and agronomic traits in large populations remains poorly understood and presents a promising avenue for further trait improvement. Recent research has advanced crop breeding by identifying superior alleles, as seen in maize (Zea mays ssp. mays L), rice (Oryza sativa), common bean (Phaseolus vulgaris L.), and wheat (Triticum aestivum)13-18. However, previous population genetic studies in polyploids have employed a simplified approach, treating polyploids as diploids for variant calling and population genetic analyses19-21, which has led to the omission of information on heterozygous loci and inaccurate estimation of allele dosage effects. With the advances of sequencing technologies and computational algorithms, it is now possible to predict allele dosage and assess its effects on phenotypic variation in polyploids. Furthermore, a comprehensive genome-wide scan of breeding footprints in polyploids to identify and evaluate allele dosage in the breeding process would provide crucial insights for crop improvement. Sweetpotato [Ipomoea batatas (L.) Lam.] is a globally important food crop, with an annual production of approximately 86.4 million tons in 2022 (FAOSTAT; http://faostat.fao.org/). China is the world’s largest producer and consumer of sweetpotato, accounting for ~54% of the global production. Given the hexaploid nature of sweetpotato, exploring genetic variants with dosage effects and understanding their roles in sweetpotato breeding is critical. In this study, we deepsequenced 294 sweetpotato accessions and systematically investigated the selection of agronomic traits during sweetpotato breeding, the population structure of these accessions, and the involvement of allele dosage in sweetpotato breeding. We identified genome regions with allele dosage effects on 23 agronomic traits and conducted transgenic confirmation for candidate genes/variations related to tuberous root weight and flesh color. This study highlights the importance of allele dosage in influencing agronomic traits and provides valuable resources of superior alleles to accelerate sweetpotato breeding efforts. Results Phenotypic evolution during sweetpotato breeding To analyze the impact of allele dosage on phenotypic variation and polyploid sweetpotato breeding, we collected 294 accessions spanning different stages along the sweetpotato breeding history (Supplementary Table 1). This collection comprised 250 accessions from China, as well as 44 from other countries that have been frequently used as donors in sweetpotato breeding in China. Of the 250 accessions from China, 135 were landraces and 115 were cultivars (Fig. 1a and Supplementary Table 1). The 115 Chinese cultivars were further divided into two groups based on distinct breeding targets22, i.e., 51 early-stage cultivars released before 2000 with white/lightyellow flesh color (denoted as B2000) and 64 recent cultivars released after 2000 with purple/orange flesh color (denoted as A2000). To elucidate the improvement characteristics of sweetpotato, 23 agronomic traits were investigated in these 294 accessions grown in Guangzhou, China, in 2020 and 2021 (Extended Data Fig. 1 and Supplementary Table 2). Principal component analysis (PCA) of the phenotypic data for these traits largely separated the B2000 and A2000 cultivars (Extended Data Fig. 2), further supporting the classification of Chinese cultivars into these two groups. Of the 23 traits, 17 (73.91%) showed improvement during sweetpotato breeding, including five tuberous root weight traits, four tuberous root appearance traits, one flesh color trait, six high-density planting traits, and the leaf shape trait (Fig. 1b-e and Extended Data Fig. 3). Tuberous root production, the primary target for sweetpotato improvement, is determined by a combination of planting density and tuberous root weight. We observed that plant architecture traits were preferentially selected from landrace to B2000 (early stage), including more erect branches, shorter internode length, lower ground fresh weight, shorter vine length, and more compact tuberous roots, which together contribute to increased plant density tolerance (Fig. 1b). Tuberous root weight (per individual) was preferentially selected from B2000 to A2000 (later stage) (Fig. 1c). Consequently, the increase in total tuberous root production reflected a shift from high-density planting to a combination of highdensity planting and increased tuberous root weight (Fig. 1f). Notably, tuberous root appearance was preferentially selected after the 2000s due to consumer preference, with favored traits including purple skin color, smooth tuberous root surface, spindle-shaped tuberous roots, and orange flesh color (Fig. 1d,e). In summary, sweetpotato breeding has resulted in increased tuberous root weight, improved tuberous root appearance, higher flesh quality, and better adaptation to high planting density. Population structure and genetic diversity To construct a genome-wide variation map of hexaploid sweetpotato that captures allele dosage information, we deep-sequenced these 294 accessions at an average depth of 26.2× per haplotype, equivalent to 157× relative to the monoploid genome of hexaploid sweetpotato (Supplementary Table 1). The high-quality reads were mapped to the genome assembly of the wild diploid I. trifida, which has been reported to be a robust reference for hexaploid sweetpotato23 (available at https://ngdc.cncb.ac.cn; accession number PRJCA015454). The mapping rates ranged from 98.1% to 99.5%, allowing us to accurately identify full hexaploid genotypes with allele dosage, as described in a previous study24. Using GATK25, we generated a final set of 6,828,068 singlenucleotide polymorphisms (SNPs) and 1,679,020 small insertions or deletions (indels) (Supplementary Table 3 and Supplementary Fig. 1). Of these SNPs, approximately 7.38% were located in coding regions (Supplementary Fig. 2), with 3.14% and 0.06% annotated as missense and stop-gain variations, respectively (Supplementary Table 4). SNP distribution varied across chromosomes, exhibiting a profile similar to that of gene density (Supplementary Fig. 1). Within this SNP set, we calculated the percentages of SNPs with three different dosages: simplex (single dose), duplex and triplex. Simplex SNPs were the predominant type, with an average of 54.1% in each accession, while duplex and triplex SNPs represented 33.0% and 12.9%, respectively (Supplementary Fig. 3). We performed phylogenetic analysis to provide a better understanding of the history of sweetpotato breeding, which revealed that accessions from outside China sporadically clustered with landraces and cultivars from China, consistent with the frequent usage of germplasms from other regions as donors in sweetpotato breeding in China (Fig. 2a). Furthermore, population structure and PCA uncovered three distinct ancestry components (optimal K=3), corresponding to the Peruvian-Ecuadorian, Colombian-Venezuelan, and Mesoamerican gene pools identified in a previous study that analyzed 5,792 sweetpotato accessions in the CIP sweetpotato germplasm collection using simple sequence repeat (SSR) markers26 (Fig. 2b, Extended Data Fig. 4 and Supplementary Figs. 4,5). The Peruvian-Ecuadorian and Mesoamerican gene pools contained representative accessions Okinanwa100 from Japan and Nancy Hall from USA, respectively, consistent with the previous finding that these two accessions clustered with those from Peru and Ecuador and Mesoamerica in the CIP sweetpotato germplasm collection26. Accessions in the Colombian-Venezuelan gene pool exhibited genetic admixture with the Peruvian-Ecuadorian and Mesoamerican gene pools at K=2, both in this study and in the CIP collection26. We found that the Mesoamerican gene pool showed a lower genetic differentiation from the Colombian-Venezuelan gene pool (FST=0.04181) than from the Peruvian-Ecuadorian gene pool (FST=0.06399) (Fig. 2c), and that the Colombian-Venezuelan gene pool displayed intermediate phenotypic traits that were further improved in the Mesoamerican gene pool (Extended Data Fig. 5). Therefore, we hypothesize that the Colombian-Venezuelan gene pool may have originated from the PeruvianEcuadorian gene pool (primitive gene pool) and evolved into the Mesoamerican gene pool (improved gene pool). Admixture due to inter-gene pool introgression was observed in the ‘Mixed’ subpopulation, which contained ancestry components from all three gene pools (Fig. 2b). This subpopulation included 66 landraces and 74 cultivars, suggesting that these landraces could be close ancestors used for sweetpotato breeding in China. As expected, the improved gene pool (Mesoamerican) represented the predominant ancestry component in the cultivars (52.7%), while the primitive gene pool (Peruvian-Ecuadorian) and the intermediate gene pool (Colombian-Venezuelan) represented 16.0% and 31.3%, respectively (Fig. 2d). Moreover, the component of the improved Mesoamerican gene pool in cultivars gradually expanded during sweetpotato breeding (Fig. 2e). These results suggest that the Mesoamerican gene pool played a major role in sweetpotato breeding in China. Finally, we observed a modest reduction in genetic diversity during sweetpotato breeding. The genetic diversity (expected heterozygosity; HS) of landraces (8.471×10-3) was significantly higher than that of B2000 (8.407×10-3), while no significant difference was observed between B2000 and A2000 (8.402×10-3). Overall, cultivars retained 99.2% of the genetic diversity compared to landraces, suggesting a weak bottleneck during sweetpotato improvement. Capturing QTLs with dosage effects We next conducted genome-wide association studies (GWAS) using SNP dosage data with GWASpoly27, to identify QTLs with allele dosage effects on 23 agronomic traits and to evaluate their contributions to phenotypic variation and breeding. We identified a total of 650 unique QTLs significantly associated with traits measured over two independent years (2020 and 2021), along with their best linear unbiased predictor (BLUP) values (Supplementary Table 5 and Supplementary Fig. 6). Of these, 353 QTLs were associated with the BLUP phenotypic values. These 353 QTLs typically spanned less than 200 kb and encompassed 6,938 genes (Supplementary Table 6 and Supplementary Fig. 7). Among these genes, 1,488 harbored SNPs annotated as missense or stop-gain variations, as well as those located within 5 kb upstream and downstream regions of the genes. The lead SNPs within these QTLs exhibited dosage effects on phenotypic variation, indicating that changes in allele dosage were associated with proportional changes in phenotypes (Fig. 3a and Supplementary Table 7). Notably, accessions with heterozygous genotypes exhibited intermediate phenotypes compared to those with either of the two homozygous genotypes (Extended Data Fig. 6). The additive effects of these QTLs substantially contributed to trait variation in sweetpotato, with the phenotypic variance explained (PVE) by all lead SNPs for each trait ranging from 59.60% to 75.23%, e.g., 72.58% for tuberous root weight and 72.64% for flesh color (Supplementary Fig. 8). Among the 353 QTLs, 68 (19.2%) had PVE values greater than 10%, highlighting their potential as valuable genetic resources for molecular-assisted breeding (Supplementary Table 8). We revisited several previously reported QTLs and cloned genes and found that they appeared to function in a dosage-dependent manner. Four QTLs identified in our study have been previously reported28, including those associated with both overground traits (internode length and leaf shape) and underground traits (flesh color and tuberous root weight) (Supplementary Table 9). Previously identified genes, such as IbMYB129 associated with skin color and IbFBOX28,30 associated with leaf shape, were found within GWAS signals (Fig. 3b,c). Additionally, we identified several plausible candidate genes. The QTLs for flesh color on chromosomes 1, 7, and 11 contained IbCCD4, IbPSY, and IbOr, which are homologous to Arabidopsis CCD431 and PSY32, and cauliflower (Brassica oleracea var. botrytis) Or33, respectively (Fig. 3d and Supplementary Table 10). The QTL associated with tuberous root weight on chromosome 2 included an expansin gene, IbEXPA4, potentially involved in cell wall loosening and cell enlargement34 (Fig. 3e). A correlation between different allele dosages and their corresponding phenotypes was consistent with these genes functioning in a dosage-dependent manner (Fig. 3f). Contribution of gene pools to favorable allele accumulation The Mesoamerican gene pool exhibited improved agronomic traits compared to the ColombianVenezuelan and Peruvian-Ecuadorian gene pools (Extended Data Fig. 5). To investigate whether this improvement involves dosage accumulation of favorable alleles (defined as those positively associated with desirable traits such as higher tuberous root weight, increased planting density, and improved tuberous root appearance, based on lead SNPs in QTL intervals), we focused on nine major genomic loci associated with high planting density traits, tuberous root weight, and flesh color and analyzed the distribution of different favorable allele dosages across accessions with varying compositions of the three gene pools. This analysis revealed that high-dosage favorable alleles for these nine QTLs were mainly contributed by the Colombian-Venezuelan and/or the Mesoamerican gene pool (Fig. 4a). Specifically, high-dosage favorable alleles for five of the nine QTLs were primarily contributed by the Mesoamerican gene pool, three by the ColombianVenezuelan gene pool, and one by both. These findings support the hypothesis that sweetpotato improvement involved dosage accumulation of favorable alleles primarily contributed by the Mesoamerican gene pool. Tuberous root weight is a key target in sweetpotato breeding. We identified two significantly associated genomic loci carrying high-dosage favorable alleles: one primarily contributed by the Mesoamerican gene pool (Chr2_3326986) and another by the ColombianVenezuelan gene pool (Chr10_20914146) (Fig. 4a). This suggests the frequent utilization of these two gene pools in sweetpotato breeding programs to increase tuberous root weight. The Mesoamerican gene pool, represented by Nancy Hall, also carried valuable genetic potential for flesh color improvement, as it contained three major genomic loci with high-dosage favorable alleles for this trait (Fig. 4a). This finding aligns with that Nancy Hall possessed a relatively high number of favorable alleles associated with flesh color among the 294 accessions (Supplementary Fig. 9). Finally, we found that the enhancement of cultivars with more erect branches could be achieved by introducing the high-dosage favorable allele (Chr7_4547410) from the ColombianVenezuelan gene pool (Fig. 4a). Our findings highlight the contributions of the Mesoamerican and Colombian-Venezuelan gene pools, each providing specific high-dosage favorable alleles to the improved traits in modern sweetpotato (Fig. 4b). Leveraging dosage of favorable alleles to enhance breeding and those released after 2000 (A2000), reflecting the shift in the focus of sweetpotato breeding in China from processing to fresh-eating around the year 2000. A two-year field experiment was conducted for this core collection at the Baiyun Experimental Station (23º23’N, 113º26’E) of the Guangdong Academy of Agricultural Sciences in 2020 and 2021, using a randomized complete block design with three replicates. Each plot (4.4 m2) consisted of two rows (2 m long in the NorthSouth direction), with a plant spacing of 20 cm and a row spacing of 110 cm. A total of 23 agronomic traits were investigated, and the phenotype descriptors were adopted from the “Specifications and Data Standards for the Description of Sweetpotato Germplasm Resources” and “Morphological identification of duplicates in collections of Ipomoea batatas”49,50, which provides illustrative graphs of both aboveand under-ground traits to assist in classifying these traits into different categorical values. Detailed descriptions of these traits are provided in the Supplementary Note. To eliminate environmental deviations and obtain the true genetic effect of phenotypic values, phenotypic values for the 23 traits across all trials were adjusted using the best linear unbiased prediction (BLUP) method. The lme4 package (v1.1-35.5) in R was used to calculate the BLUP values for each trait (phenotypes ~ (1|year) + (1|lines) + (1|year:lines) + (1|year:replicate)). Pairwise comparisons of the trait values were conducted using Student’s t-test with the ‘stat_compare_means’ function in the R package ggplot2 (v3.5.1). Genome sequencing and SNP calling Young fresh leaves from the 294 sweetpotato accessions were sampled, and genomic DNA was extracted using the cetyltrimethylammonium bromide (CTAB) method. The extracted DNA was then fragmented by ultrasound on a Covaris E220 (Covaris, Brighton, UK), and fragments ranging from 300 to 500 bp were selected using magnetic bead size selection. For each of the 294 accessions, a sequencing library was constructed using the MGIeasy Kit (BGI) following the manufacturer’s recommendations. Briefly, the selected DNA fragments were repaired to create blunt ends and modified at the 3’-end to generate dATP sticky ends. The dTTP-tailed adaptor was ligated to both ends of the DNA fragments, which were then amplified by PCR and circularized to obtain a single-stranded circular (ssCir) library. The ssCir library was then amplified through rolling circle amplification to obtain DNA nanoballs (DNBs). These DNBs were loaded onto a flow cell and sequenced using the DNBSEQ Platform. Raw reads were processed to remove adaptor and low-quality sequences using Trimmomatic51 (v0.39). The cleaned reads were subsequently mapped to the I. trifida genome (https://ngdc.cncb.ac.cn; accession number PRJCA015454) using BWA-MEM (v0.7.17) with default parameters52. The resulting SAM files were converted to BAM files using SAMtools53 (v1.9), which were sorted according to the physical position of read alignments on the genome. Duplicated reads in the sorted BAM files were marked using the MarkDuplicates command in GATK25 (v4.2.2.0). The HaplotypeCaller function in GATK was employed to call SNPs for each accession with the parameters ‘--sample-ploidy 6 -ERC GVCF’. The GenomicsDBImport command was used to integrate the gvcf files of individual accessions, and the GenotypeGVCFs tool was used to call SNPs. Raw SNPs were filtered using the VariantFiltration tool with the following parameters: “-filter ‘QD < 2.0’ --filter-name ‘QD2’ -filter ‘QUAL < 30.0’ --filter-name ‘QUAL30’ -filter ‘FS > 60.0’ --filter-name ‘FS60’ -filter ‘MQ < 40.0’ --filter-name ‘MQ40’”. Additional filtering was performed as follows: 1) SNPs with an average read depth > the mean plus two standard deviations were filtered out, as these sites were likely located in repetitive genome regions; 2) Excessive heterozygosity sites (more than 80% of accessions with a heterozygosity rate between 0.33 and 0.67), which could represent potential misassembled regions, were also filtered out; 3) SNPs with lower read coverage (<20´) were defined as missing genotypes; 4) SNPs with a minor allele frequency < 0.05 or a missing rate > 0.5 were filtered out. Finally, a total of 6,828,068 SNPs and 1,679,029 indels were obtained. Population structure, PCA, and phylogenetic analyses To infer the ancestry components in the 294 sweetpotato accessions, we conducted a population structure analysis based on spatially explicit Bayesian clustering models for clustering, as implemented in the LEA R package54 (v3.8.3). We performed 100 repetitions for K values ranging from 2 to 6 and determined the optimal K value (3) with the lowest cross-entropy value. Initially, genotypes 0/0/0/0/0/0, 0/0/0/0/0/1, 0/0/0/0/1/1, 0/0/0/1/1/1, 0/0/1/1/1/1, 0/1/1/1/1/1, 1/1/1/1/1/1, and missing, where ‘0’ represented the reference allele and ‘1’ represented the alternative allele, were converted to 0, 1, 2, 3, 4, 5, 6, and 9, respectively. The snmf function in LEA was then used to generate a project associated with the genotype file using the parameter ‘ploidy = 6’. We assessed the robustness of the results by running snmf with four values of the α regularization parameter: 10, 50, 100, and 500. Subsequently, the Q function was used to perform the population structure analysis. PCA was performed using the glPCA function implemented in the adegenet package55 (v2.2.10). For phylogenetic tree construction, we randomly selected SNPs from the whole SNP dataset with replacement to generate 1,000 genotype files, which were used to calculate the genetic distance for each pair of accessions using the formula described in a previous study56, as implemented in the StAMPP package57 (v1.6.3). Genotypes in the input files were coded as dosage values: AAAAAA, AAAAAB, AAAABB, AAABBB, AABBBB, ABBBBB, and BBBBBB. Phylogenetic trees for the 1,000 replicates were constructed using Phylip (v3.697; https://evolution.genetics.washington.edu/phylip.html), employing the neighbor-joining method. The consensus module in Phylip was then used to generate the final consensus phylogenetic tree with bootstrap confidence levels. The R package ggtree (v1.4.11; https://bioconductor.org/packages/release/bioc/html/ggtree.html) was used to visualize the phylogenetic tree. Genetic diversity and population fixation index (FST) Genetic diversity can be estimated using the expected heterozygosity (𝐻!)58. 𝐻! was calculated using the Hs function in the adegenet package55 (v2.2.10) with the parameter of ploidy=6. Genotypes in the input file were coded as dosage values ranging from 0 to 6. FST was calculated using the “stamppFst” function in the StAMPP package57 (v1.6.3) with parameters of nboots=100, percent=95, and nclusters=50. Genotypes in the input file were coded as dosage values: AAAAAA, AAAAAB, AAAABB, AAABBB, AABBBB, ABBBBB, and BBBBBB. Genome-wide association studies GWAS were conducted using SNP dosage data (a minor allele frequency >=0.1 and a missing rate <= 0.5) with a linear mixed model that accounted for population structure and relative kinship, as implemented in the R package GWASpoly27 (v2.10). First, the set.K function was used to compute the covariance matrix for the polygenic effect. Subsequently, the GWASpoly function was used to perform GWAS using an ‘additive’ model, where the SNP effect was proportional to the dosage of the minor allele. We used the permutation test to determine the P-value threshold. A total of 1000 permutations and a genome-wide significance level of α = 0.05 were used, resulting in a threshold of P <10-5 or -log10(P) >5. Significant signals were identified by grouping adjacent significant SNPs (within 500 kb) that had a pairwise linkage disequilibrium (LD) r2 ≥0.5, as calculated using the ldseq59 package (v2.1.5) in R. The resulting genomic intervals were further assessed, and those containing fewer than two SNPs with P <10-5 or fewer than 10 SNPs with P <10-3 were filtered out. Candidate genes were identified as those located within the genomic intervals and within 50 kb upstream or downstream of the genomic intervals. Phenotypic variance explained (PVE) by QTLs We calculated the PVE for each individual QTL, as well as the total PVE for all QTLs associated with a particular trait, using a previously reported cross-validation approach60. A training set consisting of 75% of the sweetpotato accessions and a test set comprising the remaining 25% were randomly selected. The training set was fitted using the mixed linear model61 implemented in the rrBLUP package (v4.6.3; http://cran.r-project.org/web/packages/rrBLUP), and was subsequently used to predict the phenotypic values of the test set. PVE was then calculated using the following formula: PVE=1 − 𝑣𝑎𝑟(𝑦"#!" − 𝑦)/𝑣𝑎𝑟(𝑦"#!") where 𝑣𝑎𝑟(+) represents the variance, 𝑦"#!" is the phenotypic value of the test set, and y is the predicted value of the test set. RNA-Seq analysis Tuberous root tissues from six-week-old sweetpotato cultivar Guangshu87 and the IbEXPA4overexpression lines, as well as stem tissues from Guangshu87 at 11 different time points, were sampled. Total RNA was extracted from these samples using the Plant RNA kit (Omega Bio-tek) and subsequently used for RNA-Seq library construction and sequencing on the DNBSEQ Platform with the paired-end mode. Two to three biological replicates were conducted for each sample. Raw RNA-Seq reads were processed using Trimmomatic51 (v0.93) for quality filtering, and the cleaned reads were then mapped to the I. trifida genome using STAR62 (v2.7.7a). The resulting BAM alignment files were analyzed using featureCounts63 (v1.6.3) to calculate the raw counts of protein-coding genes, which were subsequently normalized to fragments per kilobase of transcript per million mapped fragments (FPKM) values. Genes differentially expressed between the IbEXPA4-overexpressing lines and wild-type Guangshu87 plants were identified using DEseq264 (v1.36.0) with a cutoff of false discovery rate (FDR) < 0.05 and |log2 (OE/WT)| ³ 1. GO terms enriched in the differentially expressed genes were identified using the clusterProfiler package65 (v4.4.4) with a cutoff of FDR < 0.05. Gene cloning and transformation Plasmid construction and transformation were conducted as previously described66,67. Briefly, for p35S:IbEXPA4-FLAG, the coding region of IbEXPA4 was amplified from Guangshu87 and cloned into the modified binary vector p35S:FLAG with the pCambia1300 backbone (Novagen). For the p35S:IbEXPA4-RNAi (IbEXPA4-RNAi) constructs, the complementary specific sequence of IbEXPA4 was amplified and cloned into the modified binary vector p35S-RNAi with the pCambia1300 backbone (Novagen). Primers used are listed in Supplementary Table 14. For sweetpotato transformation, a single colony of Agrobacterium tumefaciens strain AGL1 containing the target binary plasmid was selected and cultured at 28℃ for 2 days. The cultures were then collected by centrifugation at 5000 g for 3 min, washed with a wash buffer, and diluted to an optical density (OD) of 0.5. One ml of the Agrobacterium culture was injected into the fresh stem section of Guangshu87 using a syringe, and the stem was then transplanted into soil and selected for hygromycin resistance until tuberization. Transient expression assay For the pIbOR:LUC construct, promoter sequences (2 kb) of IbOr were amplified from Nancy Hall and Pengwei, respectively, and cloned into the Dual-LUC vector (pGreenII 0800-LUC), which contains a Firefly luciferase (LUCFly) reporter. The Renilla luciferase (LUCRen), driven by the CMV35S promoter, was used as an internal control. Primers used are listed in Supplementary Table 14. Tobacco leaves were transfected and cultured as previously described68. The activities of LUCFly and LUCRen were quantified using the Dual-LUC Reporter Assay Kit (Promega) according to the manufacturer’s instructions. The LUC activity was calculated as the relative activity of LUCFly / LUCRen. Quantitative real-time PCR analysis Total RNA was extracted using the Plant RNA kit (Omega Bio-tek), and complementary DNA was synthesized using M-MLV reverse transcriptase (Promega). qRT-PCR was performed in triplicate on a Roche LightCycler480 real-time system with the SYBR qPCR mix (Q711-02, Vazyme). IbTUA was used as the internal control, and the relative expression levels were normalized to that of IbTUA. Primers used for qRT–PCR are listed in Supplementary Table 14. Data availability Raw genomic sequencing and RNA-Seq reads have been deposited in the National Center for Biotechnology Information BioProject database under the accession number PRJNA1089346 and PRJNA1151231, respectively. The I. trifida genome is available at the National Genomics Data Center under accession number PRJCA015454. Code availability Custom codes used in this study are available on GitHub (https://github.com/xiangboabc/population-genetics-analysis-of-autopolyploid-species.git). Acknowledgements We thank Xueying Mo and Meixian Zhi for their assistance in growing plants and phenotyping, and Baobao Wang (Chinese Academy of Agricultural Sciences) and Junpeng Shi (Sun Yat-sen University) for critical reading of the manuscript. This work was supported by grants from the National Key R&D Program of China (2019YFD1000701) to Z.W., the Protection and Utilization for Crop Germplasm Resources (NWB014) to Z.W., the earmarked fund of China Agriculture Research System (CARS-10, Sweetpotato) to Z.W., the Special Fund for Scientific Innovation Strategy-construction of High-Level Academy of Agriculture Science (R2022YJ-YB1004) to X.Z., the Innovation Program of the Chinese Academy of Agricultural Sciences to J.F., and the Bill & Melinda Gates Foundation through SweetGAINS (OPP1213329) to Z.F. and RTB Breeding (CGIAR Investment ID1523-BMGF) Projects under a subcontract with the International Potato Center to Z.F. Author contributions Z.W., X.Z., C.T., J.F., and Z.F. designed and managed the project. X.Z. conducted bioinformatics analyses. J.F. participated in data analysis. C.T., B.J. and R.Z. performed phenotyping and generated DNA sequencing data. X.H., X.L. M.W., and A.C. carried out the experimental validation of candidate genes. M.L. provided the reference genome of I. trifida. Z.Y., H.Z., L.H., Y.Y., and Z.L. participated in material preparation. X.Z. wrote the manuscript. Z.W., J.F., Z.F., C.T., Y.W. and S.W. revised the manuscript. Competing interests The authors declared no competing interests. Figure legends Figure 1 Phenotypic changes along sweetpotato breeding. a Summary of the sequenced sweetpotato accessions. “Introduced” refers to accessions from outside China. Landrace, B2000 (cultivars released before 2000), and A2000 (cultivars released after 2000) represent accessions from the three distinct breeding stages. b-e Boxplots showing phenotypic changes during sweetpotato breeding: high planting density traits (tuberous root compactness, vine length, overground flesh weight, internode length, branching angle, and stem diameter) (b), tuberous root weight (number and weight of small and medium tuberous roots) (c), tuberous root appearance traits (tuberous root shape, surface defect, skin smoothness, and skin color) (d), and tuberous root quality traits (orange flesh color) (e). Significance tests for continuous trait data were performed using two-sided Student’s t-tests. *, ** and *** indicate significant differences at P < 0.05, 0.01, and 0.001, respectively. ns, not significant. In each boxplot, the lower and upper bounds indicate the first and third quartiles, respectively, the center line indicates the median, and whiskers represent 1.5× the interquartile range. Sample sizes are displayed above boxplots. f Model for tuberous root production improvement during sweetpotato breeding. Figure 2 Population structure of diverse sweetpotato accessions. a,b Phylogenetic tree (a) and population structure (K = 2 and 3) (b) of sweetpotato accessions. The phylogenetic tree was rooted using hexaploid I. trifida, a close wild relative of sweetpotato. Bootstrap confidence levels are shown along the tree from the root to the branches, continuing until they reach 100%. Accessions with an ancestry component ≥ 0.65 for a single gene pool at K = 3 were assigned to that gene pool, otherwise to the ‘Mixed’ subpopulation. c Population differentiation (FST) between and expected heterozygosity (HS) of the three gene pools used in sweetpotato breeding. d Percentage of the three gene pools in modern cultivars. e Percentage of the three gene pools in landraces, B2000, and A2000 cultivars. Statistical significances of differences between groups were assessed using twosided Student’s t-tests. In each boxplot, the lower and upper bounds indicate the first and third quartiles, respectively, the center line indicates the median, and whiskers represent 1.5× the interquartile range. Sample sizes are displayed above boxplots. Figure 3 Genome-wide association studies (GWAS) of sweetpotato. a Allele dosage effect on phenotypic variation. For each significant GWAS locus, the average phenotypic value was calculated for each of the seven different dosage genotypes of the lead SNP. In the heatmap, gray, purple and orange colors represent missing, low, and high phenotypic values, respectively. QTLs are listed along the x-axis, while the corresponding dosage genotype classes are shown on the yaxis. “a” denotes the allele associated with low phenotype value and “A” denotes the allele associated with high phenotype value. b-e Manhattan plots of GWAS for skin color (b), leaf shape (c), flesh color (d), and tuberous root weight (e). Dashed horizontal lines indicate the significance threshold of GWAS at P = 10-5. Plausible candidate genes are indicated. f Distribution of phenotypes in accessions with different genotypes. Letters above the box plots represent significant differences from multiple comparisons at P < 0.05 (analysis of variance followed by Duncan’s multiple-comparison test). In each boxplot, the lower and upper bounds indicate the first and third quartiles, respectively, the center line indicates the median, and whiskers represent 1.5× the interquartile range. Sample sizes are displayed above boxplots. Figure 4 Contribution of the three gene pools to the accumulation of favorable alleles for sweetpotato trait improvement. a Compositions of the three gene pools in accessions with different favorable allele dosages. The x-axis represents genotypes with favorable allele dosages ranging from D0 to D6, corresponding to 0 to 6 favorable alleles. The y-axis indicates the proportion of each gene pool within the accessions, based on the population structure analysis. In each boxplot, the lower and upper bounds indicate the first and third quartiles, respectively, the center line indicates the median, and the whiskers represent 1.5× the interquartile range. Sample sizes are displayed above boxplots. b Summary of the contribution of the three gene pools to sweetpotato trait improvement. Figure 5 Accumulation of favorable alleles along sweetpotato breeding. a Accumulation of favorable alleles in landraces, and B2000 and A2000 cultivars. *, ** and *** indicate significant differences at P < 0.05, 0.01, and 0.001, respectively (two-sided paired Student’s t-test). ns, not significant. In each boxplot, the lower and upper bounds indicate the first and third quartiles, respectively, the center line indicates the median, and whiskers represent 1.5× the interquartile range. Sample sizes are displayed above boxplots. b Pedigree information for the sweetpotato cultivar Guangshu87. Red circles indicate sequenced accessions in this study. c Accumulation of favorable alleles in the Guangshu87 pedigree. d Pedigree information for the sweetpotato cultivar Longshu9. e Accumulation of favorable alleles in the Longshu9 pedigree. f Number of favorable alleles accumulated in representative cultivars from China. The number of accumulated favorable alleles for each accession was calculated and then normalized to Z-score using the “scale” function in the R package. “*” denotes a normalized Z-score value <-1.2, indicating undesirable traits in the cultivars. Figure 6 Functional validation of candidate genes/variants for tuberous root weight and flesh color. a,g Local Manhattan plots of loci associated with tuberous root weight (a) and flesh color (g). Dashed horizontal lines indicate the significance threshold of GWAS at P = 10-5. b,c Allele dosage effect of the locus around IbEXPA4 on tuberous root weight (b) and the expression of IbEXPA4 (c). In each boxplot, the lower and upper bounds indicate the first and third quartiles, respectively, the center line indicates the median, and whiskers represent 1.5× the interquartile range. Sample sizes are displayed above boxplots. d Expression levels of IbEXPA4 in wild-type (Guangshu87) and three independent IbEXPA4-overexpressing (OE) lines. e Phenotype of wildtype and IbEXPA4-OE lines grown for 4 months under standard laboratory conditions. Arrows indicate lignified storage roots in IbEXPA4-OE lines. Scale bar = 5 cm; f Root lignification ratios of wild-type and IbEXPA4-OE lines, shown as the percentage of lignified roots to total tuberous roots. h Significant GWAS signal around the IbOr gene. i Sequence variations in the IbOr promoter between the orange-fleshed sweetpotato variety T044 (Nancy Hall) and the light-yellow-fleshed Figure 2 Figure 3 Figure 4 Figure 5 Figure 6