Interploidy Introgression Shaped Adaptation during the Origin and Domestication History of Brassica napus

Abstract Polyploidy is recurrent across the tree of life and known as an evolutionary driving force in plant diversification and crop domestication. How polyploid plants adapt to various habitats has been a fundamental question that remained largely unanswered. Brassica napus is a major crop cultivated worldwide, resulting from allopolyploidy between unknown accessions of diploid B. rapa and B. oleracea. Here, we used whole-genome resequencing data of accessions representing the majority of morphotypes and ecotypes from the species B. rapa, B. oleracea, and B. napus to investigate the role of polyploidy during domestication. To do so, we first reconstructed the phylogenetic history of B. napus, which supported the hypothesis that the emergence of B. napus derived from the hybridization of European turnip of B. rapa and wild B. oleracea. These analyses also showed that morphotypes of swede and Siberian kale (used as vegetable and fodder) were domesticated before rapeseed (oil crop). We next observed that frequent interploidy introgressions from sympatric diploids were prominent throughout the domestication history of B. napus. Introgressed genomic regions were shown to increase the overall genetic diversity and tend to be localized in regions of high recombination. We detected numerous candidate adaptive introgressed regions and found evidence that some of the genes in these regions contributed to phenotypic diversification and adaptation of different morphotypes. Overall, our results shed light on the origin and domestication of B. napus and demonstrate interploidy introgression as an important mechanism that fuels rapid diversification in polyploid species.


Introduction
Polyploidy, referring to the condition in which cells or organisms possess more than two complete sets of chromosomes, has long been recognized as an important feature in plant history (Stebbins 1971;Hilu 1993;Soltis, et al. 2009).Extensive research over the recent decades has shown that polyploidy is far more prevalent than previously thought in the evolutionary history of plants, with many plant lineages having experienced several rounds of wholegenome duplications over time (Landis, et al. 2018;Leebens-Mack, et al. 2019;Zhao, et al. 2021).This widespread and recurrent status of polyploidy has been considered one of the main driving forces to phenotypic diversification and genome evolution in plants (Soltis, et al. 2015;Van de Peer, et al. 2017, 2020).Intriguingly, many domesticated crops have been demonstrated to be nascent polyploids, and some of which have even undergone multiple rounds of polyploidy events, suggesting that polyploidy can confer preconditions for successful domestication (Udall and Wendel 2006;Cheng, et al. 2014;Renny-Byfield and Wendel 2014;Salman-Minkov, et al. 2016).However, following polyploidization, multiple changes that lead to genomic instabilities can arise, especially when coupled with hybridization as is the case in allopolyploid species.These can include meiotic instability, perturbed gene expression, and epigenetic shock, which can have a direct negative influence on survival and adaptation, depicted as "evolutionary dead end" (Arrigo and Barker 2012;Soltis, et al. 2014;Nieto Feliner, et al. 2020).Understanding the genetic mechanisms of polyploidy that leads to successful domestication is of fundamental importance for a full appreciation of the potential of crops.
There is growing evidence that introgressive hybridization (introgression) across a ploidy barrier can take place Mol.Biol.Evol.40(9):msad199 https://doi.org/10.1093/molbev/msad199Advance Access publication September 14, 2023 1 at the early stages of polyploid adaptation (Chapman and Abbott 2010;Marburger, et al. 2019;Schmickl and Yant 2021).Introgressive hybridization between species is a well-documented process across the tree of life, which can result in the transfer of small amounts of genetic material from one species into another following recurrent backcrossing (Rieseberg and Carney 1998;Baack and Rieseberg 2007).The establishment of a new polyploid species requires a degree of reproductive isolation to remain distinct from its parental species (De Queiroz 2007).Stebbins (1971) has emphasized the potential contributions of introgressions in the early stages of polyploid evolution and pointed out that introgression mainly occurred unidirectionally, from a diploid to a tetraploid species (Stebbins 1971).Such unidirectional introgression has been documented in a handful of studies that show the morphological similarity of tetraploid species to their local diploid species (Kim, et al. 2008;Chapman and Abbott 2010;Whitney, et al. 2010;Han, et al. 2015;Zohren, et al. 2016).Interploidy gene flow can provide a source of novelty for transferring genetic variations into nascent polyploids, thereby permitting polyploids to better adapt to new ecological niches (Schmickl and Yant 2021;Liu, et al. 2022).
Interploidy gene flow is not only relevant during the initial establishment of polyploid species.Even when polyploid species have differentiated to such an extent that they adapt to disparate distribution regions, hybridization might occur after secondary contact, which is especially common postdomestication in crop species (Dempewolf, et al. 2017;Janzen, et al. 2019).Crop domestication is an evolutionary process arising from ancestral species within source centers, followed by expansion to the current distribution ranges (Meyer, et al. 2012).The role of hybridization in crop diversification has been widely reported, and there is evidence for gene flow from wild relatives in domestication, for polyploid crops such as wheat (He, et al. 2019;Zhou, et al. 2020), banana (Cenci, et al. 2020), maize (Wang, et al. 2017), and Brassica (Zhang, et al. 2022).This interploidy gene flow can lead to adaptive introgression, enabling the domesticated polyploid species to be used in new agricultural environments or adapt to new cultural preferences (Schmickl and Yant 2021).Still, despite the prevalence of cross-ploidy introgression in polyploid crop domestication, it is currently not clear to what extent interploidy gene flow affects species morphology and ecology.
The Brassica genus is a relevant target system for explorations of the effects of interploidy gene flow on adaptation and domestication, as it includes both three diploid and three allopolyploid species, resulting from their pairwise hybridization (Nagaharu 1935;Cheng, et al. 2017).Together the relationship of those six species is depicted as the "triangle of U." In addition, the Brassica genus comprises highly diverse morphotypes within species as a result of artificial selection during domestication in different regions of the world (Cheng, et al. 2014(Cheng, et al. , 2017)).Among the Brassica species, Brassica napus (AC genome) provides a well-established study system because of its worldwide cultivation and economic importance and some favorable attributes, including the availability of well-established genomes, genetic transformation, and ease of resynthesis (Heslop-Harrison 2013;Chalhoub, et al. 2014;Song, et al. 2020).Brassica napus is an allopolyploid species, which has diversified into three recognized subspecies, including oil-type B. napus subsp.oleifera (rapeseed or oilseeds), tuber-type B. napus subsp.rapifera (swede or rutabaga), and leafy-type B. napus subsp.pabularia (Siberian kale or leaf rape) (Chalhoub, et al. 2014;Havlickova, et al. 2018).These morphotypes can be further clustered according to their growth habitats into winter ecotypes that require vernalization, semiwinter ecotypes that need a mild winter environment, and spring ecotypes that do not need cold treatment (Leijten, et al. 2018).
Understanding the phylogenetic relationship of a crop species is the fundamental first step in resolving subsequent analyses related to its domestication history.As one of the earlier allopolyploid crops, B. napus originated from hybridization of B. rapa (A genome) and B. oleracea (C genome) followed by polyploidization around 7500-12,500 years ago (Chalhoub, et al. 2014).Recent studies based on large whole-genome resequencing data of oiltype B. napus accessions have suggested that its direct A progenitor is the European turnip (B.rapa) and its C subgenome may have derived from the common ancestor of kohlrabi, cauliflower, broccoli, and Chinese kale (Yang, et al. 2016;An, et al. 2019;Lu, et al. 2019;Wu, et al. 2019).Still, its phylogenetic relationship and domestication history remain elusive due to insufficient sampling of wild B. oleracea accessions, and the leafy and swede subspecies of B. napus.Brassica napus is presumed to have originated geographically in the European-Central Asian regions, where ancient morphotypes of B. rapa and B. oleracea coexisted (Qi, et al. 2017;Mabry, et al. 2021).During the modern domestication process, different B. napus morphotypes have largely been cultivated sympatrically with their diploid progenitor species.Intraspecies introgression has been reported among populations of different wild and diversified crops of diploid species (McAlvay, et al. 2021;Cai, et al. 2022;Saban, et al. 2023).Interploidy introgression can happen through triploid bridges and can result in transfer of desired traits from diploids into allopolyploids resulting in genotypes with improved adaptive value (Mason and Batley 2015).Previous studies have revealed some potential interploidy introgression events in B. napus (Sun, et al. 2017;Zou, et al. 2019).In addition, active intercrossing between B. rapa and B. napus rapeseed were performed in China and Australia to increase B. napus genetic diversity (Udall, et al. 2004;Qian, et al. 2005;Chen, Zou, et al. 2010).This breeding approach has been successful in developing new rapeseed varieties with desirable traits such as high yield, disease resistance, and improved oil quality (Chatterjee, et al. 2016;Mei, et al. 2020;Zhang, et al. 2022).However, the overall extent and genomic location of introgressed regions remain unexplored among different B. napus morphotypes and so do In this study, by generating whole-genome resequencing data together with public data representing major morphotypes of B. rapa, B. oleracea, and B. napus, we analyzed the origin and demographic history of B. napus.We then perform a cross-ploidy comparison of genomic context between different morphotypes and examine the occurrence of interploidy introgression during B. napus domestication.Finally, we qualified the genomic patterns of introgressed regions and clarified how potential interploidy introgression shaped the functional adaptation among the specific morphotypes during their domestication.

Sequencing and Cross-Ploidy Variations Discovery
To fully present genetic variation between and within the morphotypes in B. napus and its progenitors B. rapa and B. oleracea, we collected publicly available resequencing data of large collections of these three species and generated resequencing data of 33 additional accessions of morphotypes with lower representation in these collections.This resulted in a total of 614 accessions from 4 species (B. napus, 283 accessions; B. rapa, 199 accessions; B. oleracea, 130 accessions; and B. nigra, 2 accessions), representing a wide range of ploidy levels, subspecies, and geographic distribution (supplementary table S1, Supplementary Material online).For B. napus, our dataset represents all the morphotypes and ecotypes, including swede, Siberian kale, winter rapeseed, spring rapeseed, and semiwinter rapeseed.Given the different ploidy levels between B. rapa (AA), B. oleracea (CC), and B. napus (AACC), the resequencing data from each accession were mapped to the corresponding B. napus ZS11 reference and further combined to build the A and C lineage SNP datasets based on a cross-ploidy pipeline (supplementary fig.S1, Supplementary Material online).For the resequenced accessions, the effective mapped read depth against the corresponding reference genome averaged 8×, with a range of 2.14× to 26.54× and the mapping rate averaged around 96.8% (supplementary table S2, Supplementary Material online).After variation discovery and filtering, variants were identified with relatively high quality, with 2,731,337 single nucleotide polymorphisms (SNPs) for the A lineage and 4,287,347 SNPs for the C lineage.Those variants were used for subsequent analyses.

Origin and Phylogenetic History of B. napus
We first aim to understand the accurate origin, diploid progenitors, and phylogeny of B. napus as the fundamental first step in studying the adaptation of this polyploid species.SNPs were applied in multiple approaches to obtain robust phylogenetic histories for both the A and C lineages.With two accessions from B. nigra as the outgroup, maximum likelihood (ML) phylogenetic trees of the A and C lineages were constructed respectively, using 4-fold degenerate sites and best-fitting model selected by IQ-TREE (fig.1A and 2A).
To avoid potentially misleading interpretations of different topologies caused by conflicts from different datasets, we also randomly selected 200k SNPs from the A and C lineages for each phylogeny construction, and the results showed similar topology at most nodes (supplementary figs.S2 and S4, Supplementary Material online).
All of the B. napus accessions formed a single clade in both the A and C lineages from nuclear data, indicating B. napus was monophyletic and originated from a single hybridization event (figs.1A and 2A).To further evaluate the hypothesis of the single origin of B. napus, we performed a comprehensive coalescent simulation analysis using fastsim-coal2, in which we compared multiple demographic models.We considered three major scenarios: 1) multiple origin from distinct diploid progenitors, 2) separate origin from the same diploid progenitors, and 3) single origin from the same diploid population.Our coalescent simulation results of both the A and C lineages provided consistent support for the scenario of a single tetraploid origin followed by admixture (supplementary fig.S6 and S7, Supplementary Material online).This result provides evidence in favor of the hypothesis that a single origin of B. napus was followed by extensive interspecies and intraspecies admixture.
The phylogeny of the A lineage revealed European turnip (B.rapa ssp.rapa) as the closest basal group to all B. napus (fig.1A; supplementary fig.S2, Supplementary Material online).In line with this, principal component analysis (PCA) also located European turnips near the B. napus accessions (fig.1B).These results are consistent with previous studies indicating European turnip as the direct progenitor of the B. napus A subgenome (Yang, et al. 2016).Among the B. napus clades, the phylogenetic clustering together with population structure resolved six genetic groups, which is highly congruent with different crop morphotypes and ecotypes (fig.1C; supplementary fig.S3, Supplementary Material online).Interestingly, we found that the two subspecies, swede and Siberian kale, were the basal group to all the other oil-type rapeseeds, indicating these vegetable/feed morphotypes were the first crops to be domesticated.As for the common rapeseeds, the winter ecotype of rapeseeds was the original type, from which spring and semiwinter ecotypes evolved.
The  2B).This suggested that the direct donor of the C subgenome is most likely wildtype B. oleracea.Among the different B. napus morphotypes, Siberian kale and swede were also shown to be the first domesticated forms before oil-type rapeseeds, although the topology between the two subspecies is inconsistent with the A lineage.To mitigate the discordant position of swede and Siberian kale between the A and Interploidy Introgression Shaped Brassica napus Adaptation • https://doi.org/10.1093/molbev/msad199MBE C subgenomes, we utilized the raw sequencing data from European turnip and wild B. oleracea accessions to represent the in silico tetraploid ancestry of B. napus.We randomly selected 150,000 SNPs from each genome to construct a phylogenetic history based on a total of 300,000 SNPs.Using the in silico ancestry accessions as the outgroup, the topology suggested tuber-type swede as the root clade in the phylogeny, which could support the hypothesis that the early domesticated forms of B. napus were used as root vegetables (supplementary fig.S8, Supplementary Material online).
We next estimated historical changes in effective population size (Ne) using the multiple sequentially Markovian coalescent (MSMC).The overall Ne changes for different B. napus populations illustrated similar dynamics (fig.3, C  and D).All populations showed gradual decline in diversity since their divergence with progenitors around 5000-10,000 years ago.This period of declining Ne continued until the recent past, followed by rapid population expansion in recent 1,000 years.

Genome-Wide Comparison of Genetic Diversity and Linkage Disequilibrium across A, C Lineages
To estimate and compare the genetic diversity across the A and C lineages, we conducted a comprehensive analysis of nucleotide diversity (π) and fixation statistics (F st ) among different groups.The overall nucleotide diversity value for B. napus (π A = 1.55 × 10 −3 , π C = 1.37 × 10 −3 ) was much lower compared with that of diploid Brassica vegetable crops B. rapa (2.17 × 10 −3 ) and B. oleracea (2.98 × 10 −3 ) (fig.3A).Significant reduction of genetic diversity between B. napus and progenitors was observed, confirming the strong domestication and founder effect during recent polyploidization.Taking the two early diversified morphotypes swede and Siberian kale as example, the nucleotide diversity of their A subgenome captured only 50.8% from that of European turnip (one-tailed t-test, P < 2.2e−16) and for their C subgenome only 28.1% of that of wild B. oleracea (one-tailed t-test, P < 2.2e−16).Among B. napus populations, asymmetric distribution of nucleotide diversity between their A and C subgenomes with π A > π C was revealed in most groups, except for Siberian kale, which showed slightly higher π C than π A .Based on the A subgenome, semiwinter rapeseed and Siberian kale had the highest and lowest genetic diversity, respectively, whereas based on the C subgenome, winter rapeseed and swede showed the most and least genetic diversity.Interestingly, nucleotide diversity of the A subgenome in swede and semiwinter rapeseed groups was considerably higher than that of the C subgenome, being 1.88-fold and 1.33-fold higher, respectively.
The fixation index (F st ) was further calculated on subgenome level to determine genetic differentiation between  3B).Among B. napus populations, swede and semiwinter rapeseed showed the two lowest F st values with European turnip, compared with other groups.This corroborates with the finding of higher nucleotide diversity on A subgenome, indicating a closer genetic relationship of the A genomes of swede and semiwinter rapeseed with their A progenitors.The linkage disequilibrium (LD) decay (indicated by r 2 ) of different groups across the A and C subgenomes generally showed that LD decay in progenitors was faster than in B. napus morphotypes, consistent with a bottleneck in both subgenomes (supplementary fig.S11, Supplementary Material online).The overall LD decay was stronger in the A subgenome than in the C subgenome.Among the B. napus populations, the level of LD varied, with Siberian kale showing the highest LD value and swede and spring rapeseed showing the lowest.

Frequent Interploidy Introgression during B. napus Domestication
Since different populations of B. napus have largely been cultivated sympatrically with diploid progenitor populations during the species expansion and domestication, our prior hypothesis was that interploidy gene flow from diploids might play a significant role in the adaptation and worldwide expansion of different tetraploid B. napus populations.To test this hypothesis, we first systematically tested for the overall scale of introgression signatures in the A and C lineages by calculating ABBA-BABA statistics (Patterson's D) and f 4 admixture ratio ( f 4 -ratio) statistics (Patterson, et al. 2012;Martin, et al. 2015).The D and f 4 -ratio statistics are commonly used to assess the evidence of gene flow and proportion of introgression between populations in genomic datasets.Both of them are based on examining excess patterns of shared alleles between potential P3 and either sister species P1 or P2 on a four-species tree model as ([P1, P2], P3, O), where O is the outgroup.The estimates of D and f 4 -ratio for the same P2-P3 species pairs displayed variation contingent on the distinct P1 populations, but overall trends were rather consistent.Hence, we used their maximal values to summarize the data and focus on the overall support for introgression between P2 and P3.We used B. nigra as the outgroup for all trios (combinations of three different populations) and estimated significance using a block jackknife approach.Interploidy Introgression Shaped Brassica napus Adaptation • https://doi.org/10.1093/molbev/msad199MBE Overall, we tested a total of 140 trios in the A lineage, including 5 morphotypes of B. napus and 14 morphotypes of B. rapa.Of these trios, 48 had a significant D value at Z-score > 4, representing a total of 13 P2-P3 groups pairs.In the C lineage with 9 morphotypes of B. oleracea, 90 trios were tested, resulting in 15 P2-P3 groups pairs of D value to be significant.This provides strong evidence for pervasive historical introgression during B. napus expansion and domestication (fig.4  and f 4 -ratio showed significant values for the introgression into semiwinter rapeseed from pak choi, caixin, taicai, komastsuna, and other Asian-derived B. rapa populations, with f 4 -ratio explaining 5-21% admixture proportion (fig.4C).These introgression events reveal that when rapeseed was adopted from Europe into East Asia, intercrossing with cultivated diploid B. rapa accessions may have been applied to adapt rapeseed to the local climate with mild winters.In the C lineage, pervasive introgression events were also found from European-derived kale (B.oleracea) into Siberian kale and rapeseed ecotypes as they were largely distributed sympatrically (fig.4D; supplementary table S4, Supplementary Material online).The proportion of introgression estimated by f 4 -ratio ranged from 6% to 22% (fig.4, E and F).
We also conducted Treemix analysis to investigate potential gene flow events in both the A and C lineages (Pickrell and Pritchard 2012).The results indicate that intraspecies gene flow is widespread in both B. rapa and B. oleracea, consistent with previous studies (Mabry, et al. 2021;McAlvay, et al. 2021).Importantly, interploidy gene flow events were also detected in both lineages.In the A lineage, we observed gene flow from European turnip to swede and from different Asian cultivated morphotypes to semiwinter rapeseed (supplementary fig.S14, Supplementary Material online).Similarly, in the C lineage, we observed gene flow into different B. napus morphotypes from wild B. oleracea, kale, cabbage, and brussels sprouts, along with extensive intraspecies and interploidy introgressions (supplementary fig.S15, Supplementary Material online).Treemix analysis can differ compared with the D statistics (Patterson, et al. 2012), but these interploidy introgression signals were detected by both methods.In addition, given that homoeologous exchange (HE) has been demonstrated to be a widespread phenomenon in B. napus, we also examined the potential impact of HE on detecting introgression signals using the pipeline from a previous study (He, et al. 2017).Our findings show the high correlation coefficient between the D values obtained with and without excluding HE regions, suggesting that HE did not impact the identification of introgression events (supplementary fig.S16, Supplementary Material online).In summary, these results show that interploidy introgression was prevalent during B. napus domestication, with a likely role in adaptation to local environments and domestication for favorable traits.

Identification of Interploidy Introgressions across the Genome and Their Characteristics
Given the prevalence of interploidy introgression from both progenitor populations into B. napus, we then investigated the potential adaptive importance of these interploidy introgression events during domestication.Putative introgression to chromosome regions for the trios with significant D values were localized through genome-wide calculation of both f d and f dM statistics (supplementary fig.S17, Supplementary Material online) (Patterson, et al. 2012;Martin, et al. 2015).Putative introgressed genomic regions were defined as the top f dM windows that summed to the genomic proportion estimated from the f 4 -ratio (supplementary figs.S18 and S19, Supplementary Material online).

Functional Significance of Interploidy Introgression during Domestication
We next selected representative interploidy introgression events and sought to determine their functional value during B. napus expansion and domestication.The first example is the introgression events from Asian-cultivated B. rapa morphotypes into semiwinter rapeseed (fig.4, A  and C).This ecotype has been recorded to be developed in East Asia from winter rapeseed over the last 200 years, where diploid B. rapa had already been domesticated into different vegetable morphotypes for thousands of years (Qi, et al. 2017).After localizing interploidy introgressions from B. rapa morphotypes into semiwinter rapeseed, 54-251 genomic regions were identified, according to different donor morphotypes, of which 2,201-6,950 genes were found to be potentially affected (supplementary fig.S18, Supplementary Material online).These genes were further examined using selective sweep analysis to determine their potential adaptive significance.We used both haplotype-based methods XP-EHH and XP-CLR to scan for candidate selective sweeps in semiwinter rapeseed with windows of the top 5% of maximum values considered as selection regions (supplementary fig.S24, Supplementary Material online).Two hundred fifty-three selective signals spanning 13.48 Mb were identified in both selective sweep methods, among which 67 signals overlapped with interploidy regions (significant according to Fisher's exact test, P < 0.001) (supplementary tables S5 and S6, Supplementary Material online).To gain further insight into the potential functions of overlapping selective introgression genes, we performed gene ontology (GO) enrichment analyses (supplementary fig.S25, Supplementary Material online).GO analyses for genes in those introgressed regions showed enriched terms involving regulation of developmental growth (GO: 0048638), regulation of flower development (GO: 0009909), seed growth (GO: 0080112), lipid storage (GO: 0019915), glucosinolate catabolic process (GO: 0019759), and response to gibberellin (GO: 0009739), revealing the potential adaptive significance for interploidy introgression.
We then conducted more detailed investigations of the potential connection between adaptive introgression and flowering time diversification, as it was an indispensable trait for semiwinter rapeseed to be successfully adapted from Europe into Asian regions.Compared with the rapeseed winter ecotype, semiwinter ecotypes have a weaker vernalization requirement as they are mainly grown in Asian regions with mild winters.We identified 631 orthologs of 306 Arabidopsis thaliana genes related to flowering time (FLOR-ID; Flowering Interactive Database) from the B. napus A subgenome (Bouché, et al. 2015).Pak choi morphotype was chosen as representative of the Asian B. rapa populations for localizing candidate regions because it showed one of the highest f 4 -ratio values in semiwinter rapeseed.We identified 105 flowering time-related genes that were among the 6,423 introgressed genes affected by pak choi (supplementary table S7, Supplementary Material online).Enrichment analysis further supported that introgressed regions were functionally enriched in flowering time genes (Fisher's exact test, P = 0.021, gene ratio: 105/6423, background ratio: 631/47,274), suggesting interploidy introgression played an important role for flowering time adaptation to the local environment.Those flowering genes in the introgressed regions included important vernalization-related genes (SOC1, FRI, SVP), circadian clock-related genes (TOC4, PRR9, CKB4), as well as photoperiod-related genes (COL9, TSF) and autonomousrelated genes (PRMT4A, FPF1, HUB1).A specific example is given by SOC1 (SUPPRESSOR OF OVEREXPRESSION OF CO 1, BnaA04G0287900ZS) on chromosome A04 and PRMT4A (PROTEIN ARGININE METHYLTRANSFERASE 4A, BnaA06G0359200ZS) on chromosome A06 (fig.5, F and G).Both genes were found in introgressed regions, which showed lower genetic divergence between the donor pak choi and receptor semiwinter rapeseed.In addition, selection analysis on those shared regions revealed higher selective sweep signals in both receptor and donor populations, indicating a potentially important role for these two regions in domestication.Taken together, those results provide support for the role of interploidy introgression from locally adapted B. rapa morphotypes in flowering time adaptation to cultivation in Asian regions during the domestication of semiwinter rapeseed.
Similar interploidy selection significance was also found in European-cultivated B. napus morphotypes swede from diploid turnip.We identified 1,284 potential introgressed regions, spanning 32.4 Mb of the whole genome (fig.4, A  and B; supplementary fig.S18, Supplementary Material online).Selective sweep analysis further showed that 239 introgressed regions overlapped with selection regions, containing 1,204 genes (significant according to Fisher's exact test, P < 0.001) (supplementary fig.S26 and supplementary tables S8 and S9, Supplementary Material online).We next performed GO enrichment analyses of those potential adaptive introgressed genes, which identified significant overrepresentations in categories related to regulation of immune response (GO:0050776), photosynthesis (GO:0015979), sucrose metabolic process (GO:0005985), and other biological pathways (supplementary fig.S27, Supplementary Material online).Among the GO categories enriched, there were many relevant genes that may have functions during swede domestication.For example, sucrose transporter genes have been regarded as important meditators to transfer sugars from mature leaves to the swelling root.(Braun 2022).We found two representative genes SWEET7 (BnaA09T0254100ZS) and SUS3 (BnaA09T0024200ZS) on chromosome A09 Interploidy Introgression Shaped Brassica napus Adaptation • https://doi.org/10.1093/molbev/msad199MBE that were affected by adaptive introgression in swede (supplementary fig.S28, Supplementary Material online).In regulation of immune response category, we identified the gene WRKY33 (BnaA05T0071200ZS), whose orthologs are involved in multiple abiotic stresses (Wang, et al. 2018).Overall, our results suggest that introgression from diploid progenitors played an important role during domestication of the different B. napus morphotypes (supplementary fig.S29, Supplementary Material online).

Discussion
Origin and Evolutionary History of B. napus Populations Species in the genus Brassica are known for their extreme diversity (Cheng, et al. 2014).Although several previous studies have reported the phylogeny and relationships of B. napus, its direct progenitors and intraspecies relationships have remained elusive (Li, et al. 2017;An, et al. 2019;Lu, et al. 2019;Wu, et al. 2019).In this study, to fully understand the phylogenetic relationship among the diverse morphotypes in B. napus (AACC) and its progenitors B. rapa (AA) and B. oleracea (CC), we generated wholegenome resequencing data together with public data, which represents the majority of morphotypes and ecotypes among these species.With B. nigra as the outgroup and multiple sources of variant datasets being used, our study has provided results about the phylogenetic relationship and origin of B. napus, as well as the adaptive importance of interploidy introgression from diploid species in the successful domestication of different B. napus morphotypes.
Understanding the origin and phylogeny of a species is always the fundamental first step in resolving other analyses and its better utilization in breeding.Brassica napus was thought to have originated from hybridization between B. rapa and B. oleracea less than 10,000 years ago (Chalhoub, et al. 2014;Lu, et al. 2019).Still, identification of specific A and C progenitor genotypes has remained elusive.In the present study, we used two different variant datasets (4-fold degenerate sites and randomly selected whole-genome variants) to build the ML phylogenetic tree for the A and C lineages, and both results showed similar topology (figs.1A and 2A; supplementary figs.S2 and S4, Supplementary Material online).In combination with the Structure and PCA analysis, our phylogenetic results revealed that the A progenitor of B. napus was closest to the European turnip (B.rapa) and the C progenitor was closest to wild type B. oleracea.Our result of A progenitor is in line with several previous studies, which confirm that European turnip is the direct progenitor of the A subgenome of B. napus (Yang, et al. 2016;Lu, et al. 2019).As for the C progenitor, Lu et al. (2019) proposed that it evolved from the common ancestor of kohlrabi, cauliflower, broccoli, and Chinese kale, although they also discussed that broader sampling of B. oleracea accessions would be helpful to better understand the complex origin.In our study, owing to the extended sampling including the nonoil type subspecies of B. napus and wild B. oleracea accessions that were underrepresented in previous research, we identified wild B. oleracea as direct progenitor of the C subgenome.This result is in line with a prior finding based on AFLP data, which also put the wild-type B. oleracea at the closest position to all B. napus accessions (Allender and King 2010).Besides, given the domesticated status of the European turnip progenitor and the absence of wild populations of B. napus, it seems reasonable to assume that the original hybridization event probably took place in a cultivated rather than a natural environment, where wild B. oleracea plants were distributed sympatrically with European turnip populations.As European turnip was largely cultivated in European regions, and wild B. oleracea was found along the coastal regions in Europe (Mabry, et al. 2021;McAlvay, et al. 2021;Cai, et al. 2022), their vicinity enabled interspecific hybridization that finally led to the speciation of B. napus.
It is debated whether the allopolyploid species B. napus is the result of a single hybridization event or from multiple events.In our research, based on multiple datasets from nuclear genomes, we found that all B. napus morphotypes formed a single clade, supporting a monophyletic origin from a single hybridization event.In addition, coalescent simulations of population quartets from distinct diploid and tetraploid populations in both the A and C lineages consistently favor scenario of a single tetraploid origin, followed by interploidy and intraspecies admixture (supplementary figs.S6 and S7, Supplementary Material online).Our result is congruent with several previous findings, which also support a monophyletic origin using nuclear genome data (An, et al. 2019;Lu, et al. 2019).On the other hand, several other studies based on chloroplast genome sequences concluded a potential multiorigin ancestry for the maternal B. rapa parent, with B. napus accessions grouping with different B. rapa morphotypes (Allender and King 2010;Li, et al. 2017).Given the prevalence of interploidy introgressions from different parental populations into B. napus, we think it is likely that those introgression events may have contributed to the diverse chloroplast genomes in the cases where B. rapa was the maternal parent.
The inconsistent topology of swede and Siberian kale we found in the A and C lineages leads to the question which morphotype evolved first during the intraspecific diversification in B. napus.To mitigate potential conflicts, we used datasets from combinations of European turnips and wild B. oleracea as the pseudo ancestry B. napus.We found that the root-vegetable swede was placed at the base position, supporting the hypothesis that swede was the first domesticated morphotype following B. napus speciation (supplementary fig.S8, Supplementary Material online).Furthermore, we mitigate the potential impact of introgression on topology inference by filtering SNPs from introgressed genomic regions.The discrepancy observed in the phylogenetic tree of the C lineage indicates that interploidy introgression from the C subgenome may have blurred inference on tree topology (supplementary figs.S9 and S10, Supplementary Material online).
Wang et al. • https://doi.org/10.1093/molbev/msad199MBE Intriguingly, we realized that the swollen hypocotyl-root vegetables in the Brassica genus including European turnip from B. rapa (AA), swede from B. napus (AACC), and root mustard from B. juncea (AABB) were all earliest morphotypes to emerge during their independent domestications (Yang, et al. 2018;McAlvay, et al. 2021).These morphotypes with enlarged root-hypocotyl tubers were cultivated before the oilseed morphotypes, possibly as the tubers provide starch and sugars as energy for human needs, whereas other crops may have been used as oil source.The swollen root-hypocotyl morphotypes of the two allopolyploid species B. napus and B. juncea represent interesting models for studying role of polyploidization and convergent evolution in the future.

Contribution of Adaptive Introgressions from Diploids to the Success of Tetraploid B. napus Domestication
Hybridization across ploidy barriers can also break down species barriers, promote genetic variations in polyploids, and bring about adaptive traits (Marburger, et al. 2019;Novikova, et al. 2020;Edelman and Mallet 2021).Intriguingly, we found evidence for interploidy introgressions from various diploids into different B. napus morphotypes, which likely played an important role in adaptation to different environments and trait domestication by fueling adaptive genetic variations.In our study, both D statistics and f 4 -ratio statistics reveal that introgression is prevalent in B. napus and has shaped appreciable proportion of extant genomes of different B. napus morphotypes (fig.4).After comparing various genomic landscapes, introgressed regions were shown to have significantly higher genetic diversity and were associated with regions with higher recombination frequency (supplementary fig.S23, Supplementary Material online).This was in accordance with previous studies, which suggested that high recombinant regions tend to be more permissive to introgression because selection is more effective in separating neutral or beneficial alleles from deleterious alleles in areas with high recombination rates (Schumer, et al. 2018;Suarez-Gonzalez, et al. 2018).Although our dataset includes the majority of morphotypes and ecotypes in the three Brassica species, it is possible that some unsampled "ghost lineages," including extinct ones, may have had direct interploidy introgression between the two diploid species and the allopolyploid species (Dagilis, et al. 2022;Tricou, et al. 2022).
We investigated whether those interploidy introgression events played a role in diversification and adaptation of semiwinter rapeseed and swede morphotypes to local environments.The semiwinter rapeseed was diversified from winter rapeseed based on our phylogeny and was recorded to be brought to Asian region to replace original oil-type B. rapa as new oilseed crop less than 200 years ago (Sun, et al. 2017).Swede represents a relatively old crop, which was largely cultivated in European regions (Gowers 2010).Both morphotypes have been cultivated sympatrically with various B. rapa and B. oleracea populations.Our introgression analysis reveals that there was indeed consistent interploidy gene flow from Asian cultivated B. rapa morphotypes (like pak choi, caixin, Chinese cabbage, etc.) into semiwinter rapeseed and from European turnips into swede.This could possibly explain the higher nucleotide diversity of their A subgenomes compared with the A genomes of the other morphotypes like spring and winter oilseed (fig.3A).Such interploidy introgression can mitigate the severe bottleneck effect of allotetraploid B. napus.We identified several functional important genes, such as those related to flowering time diversification, in the introgressed regions, which colocated with strong selective sweep signals.Interestingly, in swede, the genes selected in introgressed regions related to sucrose transport, important for its fast-growing tubers, and immune response, likely related to their long growth and presence in soil.This study provides evidence that frequent interploidy introgressions from the congeneric diploids are important factors for successful domestication and adaptation of different B. napus morphotypes to a wide range of local environments.

Plant Materials and Sequencing
Seventeen accessions of swede/rutabaga (B.napus subsp.rapifera), 5 accessions of Siberian kale (B.napus subsp.pabularia), and 11 accessions of kale (B.oleracea subsp.acephala) and wild B. oleracea were collected for whole genome resequencing.Seeds were planted in the greenhouse during spring 2020 in order to confirm their morphotypes.Together with public data from previous studies (Cheng, et al. 2016;An, et al. 2019;Lu, et al. 2019;Wu, et al. 2019), the diversity panel used in this study contains 283 accessions of B. napus, 199 accessions of B. rapa, and 130 accessions of B. oleracea, representing most cultivar morphotypes and ecotypes of AA, CC, and AACC genomes.Besides, we also collected two accessions of B. nigra from public resources as the outgroup (Perumal, et al. 2020).Detailed information about morphotypes and ecotypes and geographic origins for each accession are shown in supplementary table S1, Supplementary Material online.

DNA Extraction and Sequencing
Genomic DNA was extracted from leaf tissue for each accession using Qiagen DNeasy plant kit.Libraries with an insert size of 350 bp were constructed according to the standard manufacturer's protocol.Paired-end reads (2 × 150 bp) were generated using an Illumina NovaSeq 6,000 platform at Novogene-Tianjin.Samples were sequenced with an average depth of 15×.Reads with more than 5% of "N" bases or with more than half of bases having quality value less than Q20 were removed from the raw data.

Variant Calling and Quality Control
The raw reads data from all accessions were filtered to remove adapters and low-quality bases using Trimmomatic (Bolger, Interploidy Introgression (Zhang, Cai, et al. 2018) (AA, version 3), B. oleracea JZC (Cai, et al. 2020) (CC, version 2), and B. napus ZS11 (AACC) were selected (Song, et al. 2020).Filtered B. napus sequencing reads were aligned to AACC genome, and B. rapa and B. oleracea data were mapped to the corresponding ZS11 AA and CC subgenomes using BWA-MEM (Li and Durbin 2009) (version 0.7.12) with default parameters.The alignment results were sorted, and PCR duplicates were marked by Sambamba (Tarasov, et al. 2015) (version 0.7).Following Genome Analysis Toolkit (GATK) Best Practices (McKenna, et al. 2010), Variant calling for each accession was then carried out with the GATK HaplotypeCaller module and consolidated into a single GVCF file model, from which SNPs and InDels were finally identified using joint calling approach.

Detection of Syntenic Sites for A and C Lineages
Given the different ploidy level across species, we used a cross-ploidy pipeline, and only syntenic or conserved SNPs from the A and C subgenomes were retained to build A and C variants lineages (supplementary fig.S1, Supplementary Material online).NUCmer module from MUMmer4 (Marçais, et al. 2018) (version 4.0.2) was used to align the two subgenomes in the A and C lineages, respectively, with default parameters, following which the delta-filter module was used to obtain one-to-one syntenic blocks in the alignment results with parameters "-r -q."Finally, SNPs retrieved from the A subgenome of B. rapa and B. napus were combined and filtered by syntenic blocks.These SNPs were defined as the A lineage variants.SNPs retrieved from the B. oleracea and B. napus C subgenome were similarly combined and filtered to represent the C lineage variants.

Phylogenetic Inference and Population Structure
To construct a ML phylogenetic tree, we used 4-fold degenerate sites SNPs to reduce the potential influence of natural or artificial selection.Besides, we also randomly select 200k SNPs from corresponding subgenomes for each phylogeny construction.Phylogenetic tree for each lineage was reconstructed using IQ-TREE (Minh, et al. 2020) (version 2.0.3), based on the best fitting model (TVM + R10), determined by the Bayesian information criterion.Bootstrap values were calculated using the ultrafast bootstrap method (UFboot) with 1,000 replicates.The output tree was then plotted and visualized by R package ggtree (Yu 2020) (version 3.0.2) with two accessions from B. nigra as the outgroup.
PCA was performed by PLINK (Purcell, et al. 2007) (version 1.90b4).The top three principal components were used and plotted in R. We then used fastSTRUCTURE (Raj, et al. 2014) to perform population structure analysis for each lineage, with the number of clusters (K) been set from 2 to 20.The optimal K, which maximizes the marginal likelihood, was evaluated by the script chooseK.py.The output structure results were further visualized using pophelper (Francis 2017) package in R.

Estimation of Demographic History
We used MSMC2 software (Malaspinas, et al. 2016) based on MSMC approach to estimate historical patterns of effective population sizes for each lineage.We selected 10 accessions from each population with high mapping depth.The mask files for the genome were calculated using the bamCaller.pyscript of msmc-tools, where genomic sites having less than average coverage were filtered out.In addition, the software SNPable (http://lh3lh3.users.sourceforge.net/snpable.shtml)was used to create a mappability mask for the reference genome.Sites with the majority of overlapping 100-mers without mismatch were defined as the SNPable sites and used for the following MSMC analysis.We then estimated effective population size using MSMC2 with the pattern parameters "1 × 2 + 25 × 1 + 1 × 2 + 1 × 3." To convert the coalescent scaled time to absolute time in years, we used a mutation rate of upper 1.5e−8 and lower 9e−9 from previous research and a generation time of 1 year.The output files were then plotted and visualized in R software.
We also performed a comparative analysis of demographic models and estimated parameters utilizing the coalescent simulation software, fastsimcoal2, to determine whether B. napus had a single or multiple origins, as established in the previous pipeline (Monnahan, et al. 2019).Population quartets consisting of representatives from both parental diploid populations and tetraploids were utilized for this analysis.The models were constructed with varying topologies and accounted for the presence or absence of migration (admixture) events.Each model was then fitted to a multidimensional site frequency spectrum calculated from the 4-fold degenerate SNP data.By using these methods, we aimed to discern the most plausible demographic scenario for the origin of B. napus.For each demographic scenario and population quartet, we conducted 30 independent fastsimcoal2 runs.We then extracted the best likelihood partition for each fastsimcoal2 run, calculated the Akaike Information Criterion (AIC), and aggregated the AIC values across the 30 independent fastsimcoal2 runs over the scenarios tested within each population quartet.The scenario exhibiting a consistently Wang et al. • https://doi.org/10.1093/molbev/msad199MBE lower AIC value (ΔAIC > 2) within a particular population quartet was considered to be the preferred model.

Genetic Diversity and Divergence Analysis
We calculated the nucleotide diversity for each subgenome in a window of nonoverlapping 10 kb using VCFtools.Weir and Cockerham's F ST value and D xy value between two populations were calculated using VCFtools with a window size 10 kb.Average F ST values were then calculated to represent the mean F ST value between two populations.Linkage disequilibrium (r 2 ) was calculated by PopLDdecay (Zhang, Dong, et al. 2018) (version 3.41) with parameters "-MaxDist 1000 -Het 0.1 -Miss 0.1."

Historical Recombination Rate Estimation
SNPs for each lineage were firstly phased using Beagle (Browning, et al. 2018) (version 4.1) default parameter settings and 30 iterations in a 50 kb sliding window.Then, FastEPRR (Gao, et al. 2016) (version 2.0) was employed for estimating population recombination rate, which is denoted as ρ = 4N e r (N is the effective population size and r is the recombination rate of the window), with a 50 kb nonoverlapping window size.

Detection of Past Introgression
To estimate the signals of past introgression among different morphotypes in the A and C lineages, the overall Patterson's D statistics (ABBA/BABA) (Patterson, et al. 2012) was used to examine introgression sites with a defined tree topology for the four groups as ([{P1, P2}, P3], O).For each lineage, two accessions from B. nigra were used as the outgroup (O) to evaluate whether P1 or P2 shared more alleles with a potential introgression P3 than with outgroup.D statistics for all trios of population from the A and C lineages were calculated using Dtrios module in Dsuite (Malinsky, et al. 2021) (version 0.4) with default parameters.D statistics significantly differing from 0 indicates introgression between P1 and P3 (D < 0) or between P2 and P3 (D > 0).The overall f 4 admixture ratio for all trios, which estimates the proportion of introgressed genome from donor, was calculated using admixr (Petr, et al. 2019) (version 0.9.1).The significance of the D statistics and standard error were calculated using a block jackknifing approach.To avoid repeating data and focus on the main evidence of gene flow between P2 and P3, we used the highest estimates of D and f4-ratio for each P2-P3 pair for the following analysis.
To further locate the introgressed regions across the genome for each trio, we used the f d and its modified statistics f dM using a sliding window of 10 kb with steps of 1 kb throughout the genome (Martin, et al. 2015).We defined the putative introgressed regions as windows with highest x% of f dM values, where x was determined for each trio by the proportion of introgression estimated by the f 4 ratio statistics, following the approach used by (Morales-Cruz, et al. 2021).Filtered windows were further merged to represent the final introgressed regions.
We also applied Treemix (version 1.3) to infer migration events and population relatedness in both the A and C lineages (Pickrell and Pritchard 2012).To reduce the impact of highly linked genomic regions, we employed a 50 kb sliding window (with steps of 10 SNPs) to scan the entire genome and then eliminated any SNPs that had a strong association (r 2 > 0.2).We used the LD-pruned sites to construct a tree without any migration events and then used this tree as the basis for migration models.We built admixture trees with 1-15 migration events and evaluated the model fit for each migration event by estimating the proportion of variance explained by each migration model among all the subgroups.The numbers of migration edges were estimated aided by the optM (v0.1.5)package (Fitak 2021).The resulting gene flow and migration events were further visualized in R.

Relationship between Introgression and Genetic Characteristics
To characterize the relationship between genomic features and introgression, we compared the nucleotide diversity, F st , D xy , and recombination frequency between introgressed regions and other genomic regions.We estimated those genomic features in the same window of f dM statistics using the software described before.The boxplot of the comparison for each genomic pattern was plotted using ggplot2 in R.

Detection of Selective Sweeps
A cross-population composite likelihood ratio test (XP-CLR) was used to identify selective sweeps regions in different morphotypes based on the SNPs with less than 10% missing data (Chen, Patterson, et al. 2010).The XP-CLR score between two populations was calculated using parameters of -w1 0.005 500 10,000 -p1 0.95 for each chromosome.Genetic distances were estimated according to physical distances in a previous high-density genetic map (Yang, et al. 2017).The mean XP-CLR score was calculated using 100 kb sliding windows with a 10 kb step size.The R package GenWin was used for the normalization and detecting the boundary of genomic regions with smoothness = 1,000 and method = 4.Moreover, we also estimated the cross-population extended haplotype homozygosity (XP-EHH) for each population using Selscan, after filtering all missing data, with 50 kb sliding windows and 25 kb step size (Szpiech and Hernandez 2014).Sliding windows with average XP-CLR scores and XP-EHH scores higher than 95th percentile were selected as significant windows.

GO Enrichment Analysis
GO terms for the B. napus genes were assigned based on syntenic relationship with Arabidopsis genes detected by Synorths (Cheng, et al. 2012).To identify the biological processes of genes in introgressed regions and selective sweeps regions, GO analysis was performed with ClusterProfiler (Wu, et al. 2021) (version 4.0).Enrichment Interploidy Introgression Shaped Brassica napus Adaptation • https://doi.org/10.1093/molbev/msad199MBE significance was analyzed with Fisher's exact test.P values were further corrected for multiple comparisons using the method of Bonferroni.

FIG. 1 .
FIG. 1. Phylogenetic relationship and population structure of the A lineage. (A) The phylogeny of the A lineage with B. nigra as the outgroup.The tip colors of the phylogeny indicate subspecies/morphotypes, whereas the branch colors denote the ploidy level.Branches with reliable bootstrap value (>70) are labeled with black point at the corresponding nodes.(B) PCA of the B. rapa and B. napus accessions.The proportions of variance explained by the top three principal components are presented in the axis labels.Colored points represent different morphotypes and are the same as the tip colors in the phylogeny.(C) Model-based Bayesian clustering performed with the number of ancestry kinship (K ) set to 11.The different colors of vertical bar represent contributions to the K-groups.The colored segments of each horizontal bar indicate morphotypes.

FIG. 2 .
FIG. 2. Phylogenetic relationship and population structure of the C lineage.(A) The phylogeny of the C lineage with B. nigra as the outgroup.The tip colors of the phylogeny indicate subspecies/morphotypes, whereas the branch colors denote the ploidy level.Branches with reliable bootstrap value (>70) are labeled with black point at the corresponding nodes.(B) PCA of the B. oleracea and B. napus accessions.The proportions of variance explained by the top three principal components are presented in the axis labels.Colored points represent different morphotypes and are the same as the tip colors in the phylogeny.(C) Model-based Bayesian clustering performed with the number of ancestry kinship (K ) set to 10.The different colors of vertical bar represent contributions to the K-groups.The colored segments of each horizontal bar indicate morphotypes.

FIG. 3 .
FIG. 3. Genomic features of the A and C lineages.(A) The nucleotide diversity over 100-kb nonoverlapping windows in the A and C lineages.The middle line indicates the median value.The top and bottom of whisker denote the maximum and minimum value or the third quartile plus 1.5× the interquartile range (IQR).(B) Comparison of F st and nucleotide diversity between B. napus morphotypes and its direct progenitors.The value in each circle indicates the nucleotide diversity of the group, whereas the value on each line represents F st value between groups.(C and D) Demographic history of the A and C lineages inferred by MSMC model.Generation estimates were inferred by assuming that mutation rates were 1.5 × 10 −8 per synonymous site per generation, respectively, and that the generation time was 1 year.

FIG. 4 .
FIG. 4. Frequent past interploidy introgression during B. napus domestication.(A and D) Heatmaps indicate maximum pairwise Patterson's D statistics measurements between pairs of morphotypes across all combinations in the A and C lineages.Asterisks indicate a significant value (|Z| > 4).(B, C, E, and F ) f 4 -ratio statistics to test the proportion of interploidy introgression in specific groups.Colors of filled circles indicate a significant value threshold (|Z| > 4).The top and bottom whiskers correspond to 1 SE calculated across the A and C subgenomes using a weighted blocked jackknife.(Complete data sets are available in supplementary tables S3 and S4, Supplementary Material online.) Wang et al. • https://doi.org/10.1093/molbev/msad199MBE the potential functional and adaptive values that may influence successful domestication.
, A and D; supplementary figs.S12 and S13 and supplementary tables S3 and S4, Supplementary Material online).In the A lineage, two interploidy introgression events were found in B. napus populations (fig.4A; supplementary table S3, Supplementary Material online).The first one happened between tetraploid swede and diploid European turnip, which were largely cocultivated in European regions.The D statistics revealed a strong introgression signal in swede from European turnip compared with other B. rapa populations including Asian turnip.The f 4 -ratio further estimated around 30% proportion of introgressed genome between European turnip and swede, which is the highest value for introgression in B. na- pus (fig.4B).Another large-scale introgression event happened between the semiwinter rapeseed and most Asian-cultivated B. rapa populations.Both D statistics Wang et al. • https://doi.org/10.1093/molbev/msad199MBE recombination regions tend to be more susceptible to introgression than regions with low recombination rates (fig.5E; supplementary fig.S23, Supplementary Material online).In general, these genomic characteristics suggest some potential role of interploidy introgression for B. napus, possibly helping the various populations better adapt to local environment during domestication.