1. Introduction

Polyploidy, or whole-genome duplication (WGD), contributes substantially to speciation and genomic evolution (Otto, 2007; Van de Peer et al., 2017). However, due to high extinction rates (Mayrose et al., 2011), a reduced efficacy of selection (Stebbins, 1971), or the possibility of returning to a diploidized state (Li et al., 2021), few polyploids succeed during evolution (Arrigo and Barker, 2012; Soltis et al., 2014; Levin, 2019). Accordingly, it is not clear how polyploids persist in nature (Fawcett and Van de Peer, 2010; Madlung, 2013). Polyploid success may be evaluated by the longevity, diversification, or adaptation of derived lineages (Fig. 1; Otto, 2007; Hegarty and Hiscock, 2008; Fawcett and Van de Peer, 2010; te Beest et al., 2012; Van de Peer et al., 2017). Success based on longevity emphasizes the survival of polyploidy-derived lineages over tens of millions of years (Fawcett and Van de Peer, 2010). Ancient waves of WGD events during the history of plant evolution highlight the success of polyploid lineages from this perspective (Zhang et al., 2020). Success based on diversification reflects the ability of polyploid lineages to diversify at the species level or above (Stebbins, 1971; Otto, 2007; Kellogg, 2016). Recent work has revealed that polyploidy may promote the diversification of plant lineages, such as Allium or Brassicaceae (Han et al., 2020; Román-Palacios et al., 2020). Finally, success based on adaptation (referred to here as the polyploid adaptation hypothesis) asserts that successful polyploids are capable of consistently colonizing or adapting to new environments after speciation (Hegarty and Hiscock, 2008; te Beest et al., 2012; Van de Peer et al., 2017), irrespective of a short-term advantage in initial establishment. Surprisingly, there is still little evidence for this hypothesis (but see Ramsey, 2011; Selmecki et al., 2015; Wu et al., 2020). The polyploid adaptation hypothesis can be tested at the population level or on contemporary timescales, making it feasible to detect the microevolutionary impacts of polyploid success.

Multiple factors may contribute to adaptation in polyploids, including the initial genomic changes mediated by the speciation mode (e.g., allo- or autopolyploid: polyploidization with or without hybridization) or changes in key traits, ecological amplitudes, or genomic plasticity (Soltis and Soltis, 2000; Comai, 2005; Leitch and Leitch, 2008; te Beest et al., 2012). For example, allopolyploids may attain initial success by fixed-heterozygosity or preadapted ecological tolerances inherited from progenitors (Barker et al., 2016). To achieve establishment and avoid extinction induced by minority cytotype exclusion, polyploids are likely to be asexual, selfing, or perennial (Van Drunen and Husband, 2019; Spoelhof et al., 2020), or to show rapid ecological differentiation (Baniaga et al., 2020; Huynh et al., 2020). However, the speciation mode or associated phenotypic changes may only result in the establishment of a nascent polyploid population (Bomblies, 2020). Genomic modifications, such as unidirectional introgression from diploids during secondary contact (Han et al., 2015; Arnold et al., 2016), are more likely to improve fitness in local populations (Schmickl and Yant, 2021), without influencing allopatric populations. Therefore, to evaluate the evolutionary fates of polyploids, it is important to quantify medium- to long-term adaptation (Van de Peer et al., 2017), either across populations or in broader natural contexts.

Polyploids experiencing extreme climates may provide ideal opportunities to test their evolutionary fates, because the advantages of polyploidy might be only revealed under such environmental stimuli (Van de Peer et al., 2021), for example, extreme salt stress (Chao et al., 2013) or cold and dark conditions (Wu et al., 2020). Polyploids may also benefit from reduced competition in habitats with lower species diversity (Rice et al., 2019), such as denuded or barren areas created by Pleistocene glaciation in high-altitude temperate mountains (Wen et al., 2014).

Located in southwest China, the biodiversity hotspot of the Hengduan Mountains (HDM) provides a natural laboratory for studying the evolution of polyploids (Nie et al., 2005; Wang et al., 2017). The two phytogeographical areas in the HDM (southern and northern) show distinct differences in environments and species assemblages (Zhang et al., 2009). The terrain and vegetation in the southern HDM are more diverse than those in the north (An et al., 2011; Zhao et al., 2020), making the southern HDM a regional biodiversity center or refugium (Yu et al., 2019; Zhang et al., 2021).

Here, we tested the polyploid adaptation hypothesis using natural populations of two polyploid yellowcresses in the HDM, the endemic Rorippa elata and the widespread Rorippa palustris, by speciation modeling, comparative phylogeography, and a transplant experiment. Under the polyploid adaptation hypothesis, we expect to observe adaptation to new environment during the colonization of the HDM. R. elata and R. palustris represent a tractable experimental system to examine adaptive responses of polyploids during colonization in the HDM. First, R. elata and R. palustris are the only two Rorippa species distributed in both parts of the HDM (Al-Shehbaz, 2016). Second, both species are polyploids that have successfully colonized altitudinal ranges from 2000 to 4600 m, significantly higher than the ranges of other Rorippa species in Eurasia (Fig. S1). Third, R. elata plants have evolved a set of putative adaptive traits for high altitudes (Figs. 2 and S2), such as an enlarged plant body, showy flowers, high germination rates, perenniality, clonality, and a mixed mating system (Hu et al., 2021). By contrast, R. palustris plants are relatively weedier, with annual and short-lived life-histories, self-compatibility, and plentiful fruits or seeds (Jonsell, 1968; Les, 2018). In this study, we aim to decipher the polyploid speciation mode and phylogeographic history for R. elata and R. palustris, and then evaluate whether the two polyploids spread via adaptation in the HDM. Our study will better our understanding of the contribution of adaptation to polyploid success.

2. Materials and methods 2.1. Population sampling

Fieldwork was performed in southwest China in Yunnan, Sichuan, Qinghai, Gansu, and Tibet (Xizang) between 2017 and 2019. A total of 346 Rorippa elata individuals were sampled from 46 populations, covering most of the native range (Table S1; Fig. S3). For each population, we sampled 4–17 individuals (mean = 7.5) that were separated from each other by at least 10 m. In addition, samples were analyzed from 46 populations of R. palustris plus related Rorippa species (R. globosa, R. dubia, R. indica, and the newly identified R. hengduanshanensis (Zheng et al., 2021)), as well as Nasturtium officinale and Barbarea intermedia as outgroups (Table S1; Fig. S3). Given the selfing character of R. palustris (Jonsell, 1968), only one individual was sampled from most populations. Seeds or healthy leaves were collected from each sampled individual. In total, 27 and 18 populations of R. elata, and 21 and 17 populations of R. palustris were located in the southern and northern HDM, respectively.

To investigate the polyploid speciation histories of Rorippa elata and R. palustris, 10 samples from seven Rorippa diploids (2n = 2x = 16), covering nearly all continents were selected as subgenomic representatives (Bleeker et al., 2002; Table S1), including Rorippa austriaca (2 samples), R. cantoniensis (1), R. globosa (1), and R. islandica (2) from Eurasia; and R. curvisiliqua (1), R. sessiliflora (1), and R. sinuata (2) from North America. Two samples of Nasturtium officinale were used as outgroups.

2.2. Flow cytometry

Genome sizes of germinated seedling samples were estimated by flow cytometry (FCM) as described previously (Doležel et al., 2007), using a BD FACSVerse Flow Cytometer (USA) with Capsella rubella (haploid genome size = 219 Mb) as an internal standard (Slotte et al., 2013). For each sample, at least 5000 nuclei were collected, with a CV threshold for gated FCM peaks of < 5%. The ploidy level was calibrated using the ratio of genome sizes between the tested sample and diploid Rorippa islandica (2n = 2x = 16; haploid genome size = 245 Mb) (Koenig and Weigel, 2015). Chromosome numbers for representative samples with small, middle, or large genome sizes were counted according to the method of Han et al. (2015). In total, the genome sizes of 186 R. elata samples from 92 populations and 62 R. palustris samples from 38 populations were estimated by FCM (Fig. S4).

2.3. Sequence data

DNA was extracted using a modified CTAB method (Doyle and Doyle, 1987) from either fresh leaves of seedlings germinated in a growth chamber or silica gel-dried materials obtained through fieldwork. For phylogeographic reconstruction, four of ten plastid markers, psbC-trnS, trnG-trnM, trnL, and trnL-trnF, were tested and filtered for amplification efficiency and nucleotide polymorphism across 72 randomly selected Rorippa elata samples. In addition, one ITS marker was included. PCR primers were designed based on previous work (White et al., 1990; Taberlet et al., 1991; Stanford et al., 2000; Nakayama et al., 2014) (Table S2). PCR products were sequenced directly using an ABI 3730 automated sequencer (Applied Biosystems, Foster City, CA, USA).

To infer the speciation mode, we cloned and sequenced 12 low-copy nuclear loci among 14 candidates, including ACO1, CHS, CIP7, CRD1, FTSZ1-1, HY4, MCM5, MLH1, SMC2, SYP61, AT3G50910, and AT5G52810 (Table S2). PCR primers and procedures were designed and performed according to previous work (Koch et al., 2000; Stockenhuber et al., 2015; Cai and Ma, 2016). PCR products were cloned into the pEASY-T1 Cloning Vector System (TransGen, Beijing, China) and proceeded to Sanger sequencing. Sequences for 8–16 clones were collected for each locus and sample. For error-free sequences, raw clonal sequencing data were processed, and haplotypes or genotypes were then inferred using PURC and Fluidigm2PURC scripts in Python v.2.7 (Rothfels et al., 2017; Blischak et al., 2018). A total of 12 regimes with different clustering thresholds or UCHIME settings were implemented in two sets (Fig. S5). Nucleotide sequences were aligned and edited using Geneious Prime v.2019.0.3 (Kearse et al., 2012).

The sequences can be accessed by GenBank IDs: MW031315-MW031737, and MW120371-MW122062.

2.4. Speciation modeling

Phylogenetic methods were used to infer the speciation mode (autopolyploidy vs. allopolyploidy) for Rorippa elata and R. palustris according to a proposed pipeline for polyploid phylogenetic studies (Rothfels, 2021). First, single-gene trees for these two species and seven diploid progenitors were constructed for all 12 nuclear loci using BEAST v.1.8.4 (Drummond et al., 2012), as described previously (Han et al., 2020). Based on the topologies of these single-gene trees, subgenomic origins of PURC-inferred homoeologs in R. elata or R. palustris were classified as g1 and g2 for each locus, according to their phylogenetic positions within the backbone built for seven diploid Rorippa species (Fig. S6). Second, a concatenated sequence matrix was generated for each species based on the classified homoeologs and was used to construct a Bayesian species tree in BEAST v.1.8.4. Third, using the gene-tree reconciliation algorithm implemented in GRAMPA (Thomas et al., 2017), speciation modes for polyploid R. elata or R. palustris were distinguished (e.g., autopolyploidy, allopolyploidy, or no polyploidy). The species tree was used to generate two multi-labeled trees (MUL-trees) with g1 or g2 tips dropped. Based on each single-gene tree, 1000 new gene trees were generated with simulated gene gains or losses in JPrIME (https://github.com/arvestad/jprime). A matrix of 12, 000 simulated gene trees was used in the two sets of GRAMPA runs for g1 or g2. The MUL-tree with the minimum parsimony score was selected as the best-fitted species tree. The speciation mode was identified according to the topology of the selected MUL-tree.

2.5. Phylogeography and ecological niche modeling

To evaluate evolutionary relationships among Rorippa haplotypes, both reticulate and bifurcating trees were generated. Haplotypes of concatenated plastid sequences or ITS genotypes were generated separately using DnaSP v.6.10.04 (Rozas et al., 2017), with consideration of indels. The geographic distribution of plastid haplotypes or ITS genotypes was checked for the 72 R. elata samples used during marker screening tests (Fig. S7). A dated phylogenetic tree was generated for Rorippa plastid haplotypes and outgroups (Barbarea intermedia and Nasturtium officinale in the tribe Cardamineae) using BEAST v.1.8.4 (Drummond et al., 2012), with secondary calibration at the crown of the tribe Cardamineae according to the 95% highest posterior density (95% HPD: 10.2–18.6 Myr, million years), based on the previous estimation in Guo et al. (2017). Additionally, a phylogenetic network was constructed using TCS v.1.21 and visualized using tcsBU (Clement et al., 2000; Santos et al., 2015). The maximum connection unit for TCS inference was defined as the 95% limit (= 20 mutation steps), and indels were treated as the fifth state. Estimators of genetic diversity, including the haplotypic diversity H_d and the nucleotide diversity π, were evaluated for sites without indels using DnaSP v.6.10.04.

To evaluate the geographic pattern of haplotypic variation, the distribution and frequency of haplotypes were mapped for each population. The ancestral geographic ranges were reconstructed along the dated BEAST tree without outgroups using BioGeoBEARS v.1.1.2 implemented in R v.4.0.2 (Matzke, 2018). The best-fitted biogeographic model was tested under six scenarios (Table S3). Two bioregions, southern HDM (S) and northern HDM (N), were assigned to haplotypes according to the Flora of the Pan-Himalaya (Al-Shehbaz, 2016).

Bioclimatic variables obtained from WorldClim (Fick and Hijmans, 2017) and Oscillayers (Gamisch, 2019) were used to model potential ecological niches of R. elata and R. palustris using MaxEnt v.3.4.3 (Phillips et al., 2006). Five bioclimatic variables with low correlations (absolute Pearson's correlation coefficient R < 0.70; Fig. S8) within each species were selected, including mean diurnal temperature range (bio2), isothermality (bio3), mean temperature of coldest quarter (bio11), precipitation of wettest month (bio13), and precipitation seasonality (bio15). Based on multiple field collections, 119 and 102 geo-referenced occurrence records (pairwise geographic distance > 300 m) for R. elata and the HDM R. palustris were used, respectively (Table S4). The ecological niches were first modeled across three periods, the last interglacial (LIG: 0.14–0.12 Ma, million years ago), the last glacial maxima (LGM: 0.021–0.018 Ma), and the present, with bioclimatic data obtained from WorldClim, at a spatial resolution of 30 arc-seconds. To accurately reflect niche evolution of R. elata and R. palustris across historical periods in the HDM, ecological niche modeling (ENM) was then performed using bioclimatic data extracted from Oscillayers, ranging from 0.46 to 0.01 Ma with 10 Kyr (thousand years) intervals, at a spatial resolution of 2.5 arc-minutes. For each ENM panel per species per time period, the range size was estimated as the number of grids at different probabilities of modeled niches (e.g., ≥0.5, medium; ≥ 0.7, high; and ≥0.9, full) using DIVA-GIS v.7.5.0.0 (Hijmans et al., 2012). The latitudinal borders were recorded at the southern margin (SM) or northern margin (NM) for each species across the HDM.

2.6. Common garden experiment

A field common garden was constructed in Shangri-La, Yunnan, a potential site of glacial refugium for Rorippa elata. The site was located in the Shangri-La Alpine Botanical Garden (SABG), at an latitude of 27.908°N, longitude of 99.64°E, and altitude of 3304 m (Hu et al., 2021).

For Rorippa elata, 633 seedlings were germinated from seeds of 150 accessions and 77 populations. They were then transplanted into the SABG on April 29, 2019, after 8 weeks of propagation in a growth chamber under long-day (LD) conditions (16 hours light and 8 hours dark) at 20 ℃ and 60% humidity. Plants were randomized into 5 × 5 plots (75 × 75 cm²) using custom R code (https://github.com/Ting-Shen/Rorippa-polyploids-in-HDM), with 15 cm as the interval between individuals. At the end of August 2019, approximately 91% of plants (576) survived, excluding 17 living individuals that were virtually destroyed by rodents in one plot. The main data sets with enough replicates (≥10) per population were collected from 29 populations across the HDM in 2020 (with no observed rodent activity). In each population, the mean number of individuals was 13.9, ranging from 10 to 22 plants per population. Several potential fitness-related traits (FRTs) were measured in 2020 (Table S5), including the survival rate, flowering time, number and length of primary sprouts or secondary inflorescences, total number of fruits, fruit length, and number of ovules per fruit.

All traits were categorized into three groups: survival, plant architecture, and reproduction (Table S5). First, two survival rates (SR) were recorded, the overwinter survival rate based on the survival or mortality of plants in May 2020 (SR_01) and the juvenile survival rate at the end of the growing season, collected in late August of 2020 (SR_02). Second, multiple traits related to plant architecture were measured at the end of the growing season, including the total number (TN_), mean length (ML_), or mean basal diameter (MD_) of primary sprouts (PS) or fruiting primary sprouts (PSF) and secondary inflorescences (SI). Plant height (PH) was quantified based on the longest primary sprout. Third, reproductive characters were recorded, including the flowering time (FT), total number of fruits (TN_F), mean fruit length (ML_F), and mean number of ovules per fruit (MN_O). In particular, for MN_O, about 2160 fruits were counted for 44 populations in the common garden, with an average of 50 fruits (30–105 fruits) per population. Fitness for each accession was estimated based on the rate of survival to reproduction age and the total number of seeds produced, including three components (SR_02) × (TN_F) × (MN_O).

To test whether local Rorippa polyploid species (R. elata and R. palustris) were better adapted to the environment of the HDM, we quantified germination rates for R. elata, R. palustris, and other lowland diploids or polyploids in the SABG common garden. Eight Rorippa species were sampled, including three diploids, Rorippa barbareifolia, R. cantoniensis, and R. globosa; three tetraploids, R. dubia, R. elata, and R. palustris; and two hexaploids, R. hengduanshanensis and R. indica. For each species, 10 accessions (except seven for R. globosa and nine for R. indica) were selected from collections in or around the HDM, with six replicates for each accession and 10 seeds for each replicate. All replicates were randomized into 5 × 10 plots. Seeds were stratified in 0.1% agar and 3.3 mg L⁻¹ gibberellin in the dark at 4 ℃ for 10 days and then sowed into soil (nutrient soil to vermiculite = 1:3) directly. The number of germinated seedlings in each replicate plot were counted at 18 days after sowing (18DAS).

2.7. Growth chamber experiment

We investigated flowering responses to vernalization treatment for 27 randomly selected Rorippa palustris accessions from the HDM. Because R. elata cannot flower even under long-term vernalization treatment (e.g., 2 months), we did not include these plants in the growth chamber experiment. Seeds were surface-sterilized and germinated on half-strength MS media and then stratified in the dark at 4 ℃ for 7 days. After 10 days of germination under LD conditions at 20 ℃, seedlings were transplanted into soil (nutrient soil to vermiculite = 1:3) and grown in the growth chamber with LD conditions, at 20 ℃ and 60% humidity. For the vernalization treatment, 7-day-old seedlings were moved to 4 ℃ and short-day (SD) conditions (8 hours light and 16 hours dark) for 4 weeks and then grown under 20 ℃ and LD conditions. Three to six biological replicates were used for each accession. Individual plants were randomized to control for position effects. Bolting time (BT) was recorded as the time to the first sprout or elongated inflorescence from rosettes. Flowering time (FT) was based on the date of the first opened flower. To evaluate vernalization sensitivity, the differences in BT or FT between non-vernalization and vernalization treatments were calculated for each accession.

2.8. Statistical analysis

The effects of source latitude or altitude or their interaction on traits, as well as the relationship between traits and geographic parameters or fitness components, were evaluated by analysis of variance (ANOVA) and a correlation analysis. Variables were treated as putative proxies for selective pressures (e.g., latitude) or variables under selection (e.g., FRTs). To test combined effects, traits belonging to the same functional category were treated as one multivariate response variable (Table S5); the joint influence of source latitude or altitude was then analyzed by a multivariate analysis of variance (MANOVA). The relationships between the single trait variable and source latitude or altitude or estimated fitness were evaluated by Pearson's correlation coefficients (R) across populations. Traits significantly related to fitness components (e.g., SR_02, TN_F, or MN_O) or fitness itself were classified as FRTs. The fitness function was quantified based on the Pearson's correlation coefficient R between the trait and fitness component, and was used to quantify how selection would act on quantitative FRTs. In addition, to test whether the above relationships (e.g., between source latitudes and FRTs) are obscured by the effects of other conditioning processes (e.g., source altitudes), coplots conditioned on additional variables were generated when exploring the focal relationship. Multiple comparisons with least square means (lsmean) were performed to test differences in germination rates among Rorippa species, with Bonferroni adjustment of P-values using the lsmeans function in the 'lsmeans' package. All statistical analyses were conducted in R v.4.0.2 (R Core Team, 2018) according to descriptions in the R Book (Crawley, 2013).

3. Results 3.1. Speciation mode

The estimated mean (±SD) genome sizes were 562.431 ± 11.219 Mb for Rorippa elata (Fig. S4; range, 541–587 Mb) and 545.750 ± 10.779 Mb (514–568 Mb) for R. palustris. In comparison with diploid R. islandica (245 Mb; 2n = 2x = 16), samples of both R. elata and R. palustris were identified as tetraploids with 2n = 4x = 32, with further verification by chromosome counting (Fig. S4).

We inferred the hybrid speciation history for tetraploid Rorippa elata and R. palustris. First, after two sets of PURC runs under six analysis regimes, the number of inferred homoeologs at each locus decreased to two across regimes (Fig. S5). Second, according to phylogenetic topologies of 12 single-gene trees (Fig. S6), one homoeolog of R. elata or R. palustris clustered with diploid R. globosa across all gene trees and the other homoeolog of R. elata diverged from all other Rorippa species, except at MLH1. The other homoeolog of R. palustris clustered with diploid R. islandica. These two kinds of homoeologs were distinguished as subgenome g1 or g2. Third, based on two GRAMPA runs, the MUL-trees with the lowest parsimony score supported a mode of allotetraploid speciation for both tetraploids (Fig. 3a-b). According to the Bayesian species tree built on concatenated sequence data, both R. elata and R. palustris showed hybrid origins within the Eurasian clade (Fig. 3c-d). Based on the phylogenetic relationships, one ancestral progenitor of R. elata or R. palustris was inferred as closely related to R. globosa and the other progenitor was identified as an unknown species, possibly extinct or unsampled for R. elata and closely related to R. islandica for R. palustris.

3.2. Genetic diversity

Haplotypic or genetic diversity varied among species as well as geographic regions. Based on the investigated plastid and ITS nucleotides, 17 plastid haplotypes and 14 ITS genotypes were phased for Rorippa and related plants when including indels (Figs. 4, S9 and S10). Seven plastid haplotypes and two ITS genotypes were unique to R. elata, and one plastid haplotype and one ITS genotype were unique to R. palustris in the HDM, with four plastid haplotypes (Hap_2, 3, 6, and 14) and three ITS genotypes (Gen_3, 4, and 11) shared between plants in the HDM and Himalayas (HIM; Fig. 4, S10). For populations in the HDM, R. palustris had a higher overall plastid haplotypic and nucleotide diversity (H_d = 0.257, π = 0.00080) than did R. elata (0.158 and 0.00052; Table S6). At the regional scale, populations of R. palustris from the southern HDM harbored more genetic variation (H_d = 0.325, π = 0.00131) than did populations from the northern HDM (0.166 and 0.00012). For R. elata, genetic diversity was similar in the southern and northern HDM populations (H_d = 0.162 vs. 0.155, π = 0.00054 vs. 0.00051).

3.3. Haplotype network and phylogeny

Patterns of genetic variation were inferred through the construction of a haplotypic or genotypic network and phylogeny across Rorippa populations in the HDM. The plastid haplotypic network showed two main genetic clades, including one clade dominated by R. elata (named the R. elata clade) and another dominated by R. palustris or other Rorippa species (named the R. palustris clade; Fig. 4a). The core haplotype of the R. elata clade, labeled Hap_2, was shared by 84% of R. elata individuals. The core haplotype of the R. palustris clade, Hap_6, was shared by 78% of R. palustris individuals and 16 R. elata individuals. Six haplotypes of R. elata (Hap_3, 6, 11, 13, 14, and 16) and two haplotypes of R. palustris (Hap_2 and 12) were not within these main clades (estimated by pairwise p-distances; Fig. S9). These haplotypes accounted for only 6% of R. elata individuals or 17% of R. palustris individuals, indicating that they are rare haplotypes or relicts with extremely high pairwise divergence by the criteria of Alonso-Blanco et al. (2016). However, some of these haplotypes may also result from hybridization with related species, especially haplotypes shared among multiple species, such as Hap_3, 6, and 14 for R. elata or Hap_2 for R. palustris.

The two genetic clades were also well-supported (posterior = 0.99) in the dated BEAST tree (Fig. 4a). The R. elata clade and R. palustris clade originated c. 0.55 Ma (95% HPD: 1.30–0.14 Ma) and c. 0.86 Ma (1.94–0.22 Ma), respectively, with a crown age of c. 1.55 Ma (3.25–0.51 Ma). The core clade of R. elata (excluding Hap_5, 7, and 8) covered 92% of R. elata samples, with a crown age of c. 0.37 Ma (0.89–0.04 Ma).

3.4. Phylogeography and niche dynamics

Ancestral geographic ranges were reconstructed for Rorippa plants in the HDM along the well-supported plastid haplotypic tree (Fig. 4a). The southern HDM was the most likely ancestral area for both the R. elata clade and R. palustris clade in a BioGeoBEARS analysis. Divergence between the R. elata and R. palustris clades was mediated by an ancient dispersal event from the southern to southern or northern HDM. Multiple dispersal (two for R. elata and four for R. palustris) or vicariance events (one for R. elata and two for R. palustris) were also identified within each clade (Fig. 4a).

The geographic distribution of haplotypic variation revealed different population structures across species or regions (Figs. 4 and S10). For R. elata, no clear phylogeographic structure was observed among populations. Plastid haplotype 2 (Hap_2) was the most frequent type across populations, with a mean frequency per population of up to 91% (Fig. 4b). Only four populations were dominated by other haplotypes (i.e., the summed percentage of other haplotypes exceeded 50% per population), including two populations located in the southern HDM (Re3 and Re8), one in the northern HDM (Re44), and one in HIM (Re46). However, populations with relictual haplotypes were unevenly distributed. Among 10 populations harboring relicts or more than three haplotypes (Mao et al., 2021), three were found in the southern HDM, particularly in a region below the 29°N latitudinal line (referred to as the relict region, hereafter) along the Muli, Jianchuan, Shangri-La, and Daocheng counties, clockwise from east to west. For R. palustris, the core haplotype (Hap_6) was mainly found in the northern HDM (Fig. 4b), whereas most of the relict haplotypes of R. palustris were located in a relict region that fully overlapped with the relict region for R. elata. Based on the above phylogeographic results, the southern HDM would be a stable ancestral area for either R. elata or the HDM R. palustris.

Temporal dynamics of suitable niche ranges were reconstructed from 0.46 Ma to the present by ENM (Figs. 5 and S11-S14). Irrespective of the probability levels of modeled niches (e.g., 0.5, 0.7, or 0.9), suitable niches for R. elata or R. palustris populations expanded northward recurrently, approximately along the main rivers during interglacial periods (e.g., 0.14–0.12 Ma, LIG), and withdrew to south of the 29°N latitudinal line during glacial periods (e.g., 0.021–0.018 Ma, LGM; Figs. 5 and S11-S14). During the past 0.46 myr, the main ranges of R. elata and R. palustris expanded northward at least six times (counted as the number of peaks above the 29°N latitudinal line), with the largest expansion around the LIG (c. 0.12 Ma) and the latest one at present (Figs. 5 and S13-S14). The ranges with high habitat suitability (≥0.9) were stably localized in the southern HDM (e.g., below the 29°N latitudinal line; Fig. S14), overlapping with areas inhabited by relict haplotypes of either R. elata or R. palustris (Figs. 4b and S10b). Geologically, niche ranges during the LGM were mainly limited by the reconstructed fallen-snowline and the expansion of glaciers in the northern HDM (Fig. S12), in which areas below the lower snowline altitudes (≤4600 m) and glacier intervals were suitable for R. elata or R. palustris plants. By contrast, present niche ranges shifted to a higher altitude snowline (≤5200 m) and were restricted to the wet glacier forefields at the northern rim of the maritime-type glaciers (type I; Fig. S12). Overall, we demonstrated a clear pattern of recurrent northward movement by populations of R. elata or R. palustris in response to glacial dynamics in the HDM.

3.5. Adaptive response

Among Rorippa species in Eurasia, plants of R. elata inhabited areas with the highest altitudes (mean elevation = 3589 ± 612.096 m, P < 0.0001; Fig. S1), together with the HDM R. palustris (3192 ± 823.859 m in HDM vs. 387 ± 67.171 m in other places). Among Rorippa species found in sympatry in East Asia, tetraploids had significantly higher 18DAS germination rates (LSD test, lsmean = 3.96, df = 71) than did those of either diploids (2.02, P < 0.0001) or hexaploids (1.87, P < 0.0001; Fig. 2f) in the field common garden of SABG. Of the eight species tested, HDM R. palustris showed the highest 18DAS germination rate (lsmean = 4.73; Fig. 2f), followed by R. elata (3.90), R. dubia (3.25), R. indica (2.50), R. globosa (2.45), R. barbareifolia (2.03), R. cantoniensis (1.72), and R. hengduanshanensis (1.17). The germination rate of HDM R. palustris differed significantly from those of other species (df = 66, P ≤ 0.030) except R. elata (df = 66, P = 1.000). Accessions of R. elata had significantly higher germination rates than those of other Rorippa species (df = 66, P ≤ 0.002) except R. palustris, R. dubia, R. indica, and R. globosa (df = 66, P ≥ 0.069). The 2-year common garden experiment further revealed clinal variation in multiple FRTs in R. elata along source latitudes or altitudes (Table 1; Fig. S15).

3.5.1. Overall fitness

The estimated fitness varied from 279 ± 68.525 to 6922 ± 84.709 (estimated as the survival rate to reproduction age plus the total seed produced) across populations of Rorippa elata (Fig. S15), with the lowest value for a population from the southern HDM (27.853°N, 102.478°E) and the highest value for a population from the northern HDM (31.541°N, 99.978°E). Fitness function analysis suggested that plant architecture traits, such as plant height (Pearson's correlation coefficient R = 0.45, P = 0.032; Fig. S15) and the total number of primary sprouts or secondary inflorescence (e.g., TN_PS/PSF/SI; R ≥ 0.59, P ≤ 0.007), and reproductive traits, such as flowering time (R = 0.60, P = 0.006) and mean length of fruit (ML_F, R = 0.40, P = 0.037), were non-randomly related to fitness or its components. Therefore, these traits are candidate FRTs, functioning as indicators of the strength of the phenotypic response to selection.

3.5.2. Survival

More than 91% of Rorippa elata plants survived to the end of August 2019 in the common garden of SABG; however, only 63% resprouted in May 2020. Significant effects were detected for both source latitudes and altitudes on survival rate variation (SR_01/02; F ≥ 5.049, P ≤ 0.034; Table 1) but not for the interaction (P ≥ 0.795), particularly at the lowest (2047–3200 m, mainly located in the southern HDM) or highest levels (3560–4309 m, mainly located in the northern HDM) of the conditioned altitudinal variable (Fig. S16). The calculated mean overwinter germination rates (SR_01) per population were positively correlated with source latitudes (Fig. S15; R = 0.47, P = 0.011) but not altitudes (P = 0.96). A significant positive relationship was also observed between source latitudes and survival rates recorded in August 2020 (SR_02; R = 0.53, P = 0.003). Collectively, R. elata populations in the high-latitude northern HDM had overwintering germination rates of 68% ± 19% and juvenile survival rate of 58% ± 13%, compared with values of 55% ± 18% and 49% ± 19% for populations in the southern HDM.

3.5.3. Plant architecture

Variation in primary sprout or secondary inflorescence traits of Rorippa elata differed with respect to geographic variables. We detected significant effects of source latitude on most FRTs of primary sprouts (PS, including PH, TN_PS/PSF, and Fre_PSF; F ≥ 5.553, df = 1, P ≤ 0.027) other than mean basal diameter (MD_PS; P = 0.882), and a significant interaction effect of source latitude and altitude on combined variation in secondary inflorescences (SI; F = 5.235, P = 0.005; Table 1). There were no significant effects of source altitude on primary sprout or secondary inflorescence traits or their combined multivariate trait (P ≥ 0.09; Table 1). In particular, the total number of fruiting sprouts (TN_PSF) decreased across both source latitudes and altitudes (R = −0.45, P = 0.015; R = −0.51, P = 0.004; Fig. S15). Also, plant height was negatively correlated with source latitudes and altitudes (PH; R = −0.43, P = 0.022; R = −0.38, P = 0.044; Fig. S15). No significant latitudinal or altitudinal correlations were detected with the mean basal diameter of sprouts (MD_PS; absolute value of R ≤ 0.10, P ≥ 0.610) or the total number of secondary inflorescences (TN_SI; absolute value of R ≤ 0.23, P ≥ 0.170). However, the mean length of secondary inflorescences (ML_SI) clearly increased across source latitudes and altitudes (R = 0.47, P = 0.014; R = 0.46, P = 0.019). Taken together, the plant architecture of R. elata at high latitudes or altitudes was more compact than the more diffuse growth form at low latitudes or altitudes, characterized by shorter stems, fewer sprouts, and longer branches.

3.5.4. Reproduction

Both source latitude and altitude significantly influenced flowering time variation for Rorippa elata (FT; F = 9.913, P = 0.006; F = 19.184, P < 0.001; Table 1) but not fruit number or length or ovule number per fruit. Genotypes from higher latitudes or altitudes flowered earlier than those from low latitudes or altitudes (R = −0.42, P = 0.024; R = −0.57, P = 0.001; Fig. S15). The total number of fruits or mean fruit length did not differ across source latitudes or altitudes (TN_F and ML_F; absolute value of R ≤ 0.13, P ≥ 0.410). Statistically significant negative correlations were found between the mean number of ovules per fruit and latitude or altitude (MN_O; R = −0.42, P = 0.004; R = −0.34, P = 0.025; Fig. S15). Overall, high-latitude or altitude populations of R. elata showed early flowering and produced fewer ovules compared with populations at low latitudes or altitudes.

3.5.5. Flowering time

For Rorippa palustris, variation in both bolting time and flowering time was significantly influenced by or correlated with source latitude and altitude under non-vernalizing conditions (Table 1). The observed changes in bolting time were negatively related to altitude (R = −0.49, P = 0.009; Fig. S17). Similarly, flowering time decreased significantly from lower to higher latitudes (R = −0.49, P = 0.010) or altitudes (R = −0.45, P = 0.020). However, these patterns were not observed for plants in the vernalization treatment (P ≥ 0.100). Furthermore, the calculated vernalization sensitivity was also negatively related to source latitudes (R = −0.52, P = 0.005) or altitudes (R = −0.42, P = 0.031; Fig. S18).

4. Discussion

The success of polyploids can be attributed to a number of non-mutually exclusive factors, including their contributions to the longevity or diversification of descendent lineages or to long-term adaptability after speciation (Fig. 1). Here, we tested the polyploid adaptation hypothesis, which predicts that successful polyploids have an advantage with respect to persistence and adaptation to new environments. We found evidence for the recurrent post-glacial colonization of HDM from the south to north by two allotetraploids, R. elata and R. palustris, with adaptive alterations in fitness-related traits across latitudinal gradients.

4.1. Hybrid origins of tetraploid Rorippa elata and Rorippa palustris in the mid-Pleistocene

Polyploid success is thought to be influenced by the mode of speciation, under the assumption that allopolyploids have an evolutionary advantage (Otto, 2007; Barker et al., 2016). We inferred hybrid origins for both Rorippa elata and HDM R. palustris (Figs. 3, S5 and S6) during the mid-Pleistocene around c. 0.55 Ma and c. 0.86 Ma, respectively, in the southern HDM (Fig. 4). The estimated time of origin is younger than those for other hybrid species in the HDM, such as Picea purpurea (1.30 Ma) (Sun et al., 2014) or Ostryopsis intermedia (1.80 Ma) (Wang et al., 2021). This period corresponded with the Mid–Pleistocene Transition (1.25–0.70 Ma), with ice–age cycles shifting from symmetric 41 Kyr to asymmetric 100 Kyr, as well as an increased severity of glaciations (Chalk et al., 2017; Sun et al., 2021). These climatic fluctuations resulted in stepwise alterations in habitat connectivity, which may in turn have influenced gene flow among populations and facilitated hybrid speciation (Muellner-Riehl, 2019; Rahbek et al., 2019).

A hybrid origin could facilitate the establishment of polyploids by the co-option of preadapted genes or traits from progenitors (Comai, 2005; Otto, 2007; Wu et al., 2022). For example, it has been proposed that the establishment and spread of wild allopolyploid Aegilops wheats was favored by the eco-genetic additivity of their progenitors (Huynh et al., 2020). Similarly, such eco-genetic additivity may have promoted the success of several widespread weeds, such as the allotetraploid Capsella bursa-pastoris (Han et al., 2015) and Cardamine flexuosa (Mandáková et al., 2014), as well as the allohexaploid Echinochloa crus-galli (Ye et al., 2020), thus supporting the role of allopolyploidization in plant colonization (Stebbins, 1985). Our results indicate that allopolyploidy may have facilitated the initial establishment and later success of Rorippa elata and R. palustris, contrary to previous arguments that polyploidy had a minor role in promoting species diversity in the HDM (Nie et al., 2005; Wang et al., 2017; but see Mao et al., 2021).

4.2. Northward colonization shaped by glacial dynamics in the late-Pleistocene

Neopolyploids are common in deglaciated areas (Brochmann et al., 2004; Novikova et al., 2018), reflecting a biogeographical response of plants to the Quaternary climatic oscillations. Here, we found that Pleistocene glaciation shaped the phylogeography of Rorippa polyploids in the HDM, with recurrent northward expansion during interglacial periods and repeated southward withdrawals during glacial periods (Fig. 4, Fig. 5 and S11-S14).

The initial expansion of the core clade of Rorippa elata (covered 92% of R. elata samples, with the crown age of c. 0.37 Ma) followed weak summer monsoons in the HDM (c. 0.41 Ma) (An et al., 2011; Zhao et al., 2020) and was influenced by glacial dynamics (Figs. 5 and S11-S14). The larger interglacial expansions of R. elata may have been driven by the reduced precipitation and drought stresses associated with weakening summer monsoons in either the southern (c. 0.40 to 0.13 Ma, in Heqing Basin, 26.562°N, 100.171°E; An et al., 2011) or northern (c. 0.62 or 0.41 Ma to present, in Zoige Basin, 33.969°N, 102.331°E; Zhao et al., 2020) HDM. In agreement with paleoclimatic or ancient DNA evidence (Li et al., 1996; Liu et al., 2021), the niche ranges of R. elata or R. palustris were largely shaped by shifts in snowline altitudes (Fig. S12). For example, the northern niche margins for both species were limited by the snowline altitudes ranging from 4600 m during the LGM (c. 0.021 Ma) to 5200 m at present (Fig. S12), accompanied by the availability of ice-free and newly-opened lands in the northern HDM (c. 0.014 Ma, in Lake Naleng, 31.100°N, 99.750°E; Liu et al., 2021). The climatic region defined by maritime-type glaciers also may have played an important role in limiting the northern margins of R. elata and R. palustris (Fig. S12), supporting the influence of glacial dynamics in shaping the geography of alpine plants in the HDM (Sun et al., 2017).

In summary, we propose that the current distributions of Rorippa elata or R. palustris in the HDM was molded by the increasing altitude of the snowline and the expansion of glacier forefields during glacial retreats (Wallis et al., 2016; Cauvy-Fraunié and Dangles, 2019), as well as the reduced summer precipitation driven by weak monsoons since c. 0.41 Ma (An et al., 2011; Zhao et al., 2020). The denuded habitats created during the late-Pleistocene were stressful but lacking interspecific competition (Liu et al., 2021), which may have provided a chance for polyploids like R. elata or R. palustris to conquer. However, the causal relationship and the underlying adaptive basis for those inferences have yet to be clearly demonstrated (Mao et al., 2021).

4.3. Clinal adaptation across latitudinal gradients

The stress response is an important determinant for the success of polyploidy (Van de Peer et al., 2021) in recent or geologic eras (Chao et al., 2013; Wu et al., 2020), although unambiguous data supporting this conjecture are scarce (Ramsey, 2011; but see Čertner et al., 2019). In the HDM, species richness decreases significantly from south to north (Zhang et al., 2009), corresponding with the higher prevalence of species with hybrid origins in the south (Yang et al., 2019; Wang et al., 2021). Only a few hybrids are known to have spread northwards or expanded their ranges across the HDM (Wu et al., 2022), such as populations of Rorippa elata and R. palustris (Figs. 2, S1 and S2). How these Rorippa allotetraploids have adapted better than diploids across the HDM is an intriguing question that awaits further investigation. Our results agree with the polyploid adaptation hypothesis (Fig. 1) (Hegarty and Hiscock, 2008; Van de Peer et al., 2017), which stipulates that successful Rorippa polyploids may spread via adaptive responses to the environments in the HDM.

We tested the hypothesis in Rorippa polyploids through common garden and growth chamber experiments. Seedling survival rates were higher for Rorippa polyploids than for diploids (Fig. 2), implying an evolutionary advantage under the high-altitude environments of the HDM. For R. elata and R. palustris, clinal variation in fitness-related traits was found along latitudinal and altitudinal gradients, suggesting that adaptive changes accompanied their range expansions (Table 1; Figs. S15-S18). In general, R. elata plants at high latitudes or altitudes are biennials or short-lived perennials (e.g., high overwinter resprouting rate) with early flowering, which could confer an increased tolerance to seasonal climatic changes, such as a short growing period (Molau, 1993; Folk et al., 2020). Similarly, the annual R. palustris from the northern HDM highlands showed early flowering and low vernalization sensitivity (Figs. S17 and S18), consistent with findings in the annual model plant Arabidopsis thaliana (Stinchcombe et al., 2005).

For plants in cold or dry conditions, compact architectures may be advantageous to conserve warmth, save water, or resist wind (Körner, 2020). In contrast to lowland plants with more diffuse architectures, the highland Rorippa elata evolved a relatively compact architecture, characterized by a shorter height, reduced sprout numbers, and longer branches (Table 1; Fig. S15). A trade-off between survival and reproduction was also found in R. elata populations. Along latitudes or altitudes, survival rates increased as ovule numbers decreased (Fig. S15), perhaps reflecting a general K-selection strategy for plants inhabiting cold, density-dependent environments (Gadgil and Solbrig, 1972).

Latitudinal or altitudinal gradients may covary with other climatic or geographic factors, such as diurnal amplitude, seasonal precipitation, or rate of snow melt. In addition to shaping the ecological niches of Rorippa polyploids in the HDM, these factors may modulate the sensitivity of plants to local environments undergoing climate change. Plant responses to these variables require further investigation. In addition, the correlation-based adaptive potential estimated at one site (e.g., SABG) does not fully capture the performance of populations in other places (e.g., the northern HDM). Additional work is needed to evaluate adaptation in Rorippa polyploids across their ranges in the HDM. Further studies should also examine the specific genomic regions or genes responsible for the observed adaptive changes, along with more exact demographic inferences for Rorippa polyploids in the HDM.

5. Conclusion

We found that Rorippa elata, which is endemic to the HDM, and the widely distributed R. palustris, are allotetraploids with mid-Pleistocene origins. Both species have repeatedly colonized the HDM from south to north by taking advantage of the availability of new niches created by glacial retreats. Following colonization of the northern HDM, clinal variation in several fitness-related traits was found along a latitudinal gradient. These findings for both species support the polyploid adaptation hypothesis (Fig. 1) and suggest that polyploids may have spread via adaptation in response to environmental changes during the Pleistocene (Novikova et al., 2018). In conclusion, our results link the evolutionary success of polyploids to both intrinsic (e.g., adaptive clinal responses in fitness-related traits) and extrinsic (e.g., glacial dynamics) factors during medium-to long-term evolution, contributing to our understanding of the significance of polyploidy from the perspective of adaptation.

Author contributions

T-S.H. and Y–W.X. conceived the project. T-S.H. designed the experiments, collected the data, and performed analyses. Z-Y.H. and Z-Q.D. performed the common garden work. Q-J.Z. did the cytological experiments and phylogenetic analysis. J.L. explained the ecological results. T-S.H. drafted the manuscript. T.M-O. revised the manuscript and edited the language. All authors reviewed the final manuscript.

Declaration of competing interest

There are no conflicts of interest to declare.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (31800177, 32170224, and U1802242), and the Strategic Priority Research Program of Chinese Academy of Sciences (XDB31000000). T-S.H. is also supported by the Youth Innovation Promotion Association CAS (2020391), China Scholarship Council (201804910061), and CAS Light of West China Program. We acknowledge Ya-Long Guo at the Institute of Botany CAS for providing Capsella rubella seeds and revising the manuscript. We also thank members of the Paleoecology group, and Biogeography and Ecology group at the Xishuangbanna Tropical Botanical Garden CAS, as well as Zhendong Fang and Lin Liu at the Shangri-La Alpine Botanical Garden who helped with valuable discussions on the transplanting experiment. We thank the anonymous reviewers for helpful comments which greatly improved the manuscript.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2022.02.002.