Adaptive responses drive the success of polyploid yellowcresses(Rorippa, Brassicaceae) in the Hengduan Mountains, a temperate biodiversity hotspot

Ting-Shen Han ,Zheng-Yan Hu ,Zhi-Qiang Du ,Quan-Jing Zheng ,Jia Liu ,Thomas Mithell-Olds , Yao-Wu Xing ,**

a CAS Key Laboratory of Tropical Forest Ecology, Xishuangbanna Tropical Botanical Garden, Chinese Academy of Sciences, Mengla, Yunnan 666303, China

b Center of Plant Ecology, Core Botanical Gardens, Chinese Academy of Sciences, Mengla, Yunnan 666303, China

c Department of Biology, Duke University, Box 90338, Durham, NC 27708, USA

d University of Chinese Academy of Sciences, Beijing 100049, China

A B S T R A C T Polyploids contribute substantially to plant evolution and biodiversity; however, the mechanisms by which they succeed are still unclear. According to the polyploid adaptation hypothesis, successful polyploids spread by repeated adaptive responses to new environments. Here, we tested this hypothesis using two tetraploid yellowcresses (Rorippa), the endemic Rorippa elata and the widespread Rorippa palustris, in the temperate biodiversity hotspot of the Hengduan Mountains. Speciation modes were resolved by phylogenetic modeling using 12 low-copy nuclear loci.Phylogeographical patterns were then examined using haplotypes phased from four plastid and ITS markers, coupled with historical niche reconstruction by ecological niche modeling. We inferred the time of hybrid origins for both species as the mid-Pleistocene, with shared glacial refugia within the southern Hengduan Mountains. Phylogeographic and ecological niche reconstruction indicated recurrent northward colonization by both species after speciation,possibly tracking denuded habitats created by glacial retreat during interglacial periods.Common garden experiment involving perennial R. elata conducted over two years revealed significant changes in fitness-related traits across source latitudes or altitudes, including latitudinal increases in survival rate and compactness of plant architecture, suggesting gradual adaptation during range expansion. These findings support the polyploid adaptation hypothesis and suggest that the spread of polyploids was aided by adaptive responses to environmental changes during the Pleistocene.Our results thus provide insight into the evolutionary success of polyploids in high-altitude environments.

Keywords:Adaptation Hengduan mountains Pleistocene Polyploidy Rorippa

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 thelongevity,diversification, oradaptationof 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 onlongevityemphasizes 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 ondiversificationreflects 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 asAlliumor Brassicaceae (Han et al., 2020;Rom′an-Palacios et al., 2020). Finally, success based on adaptation(referred to here as thepolyploid 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). Thepolyploid adaptation hypothesiscan be tested at the population level or on contemporary timescales, making it feasible to detect the microevolutionary impacts of polyploid success.

Fig. 1. Overview of hypotheses for evaluating 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 highaltitude 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 thepolyploid adaptation hypothesisusing natural populations of two polyploid yellowcresses in the HDM,the endemicRorippa elataand the widespreadRorippa palustris,by speciation modeling, comparative phylogeography, and a transplant experiment.Under thepolyploid adaptation hypothesis,we expect to observe adaptation to new environment during the colonization of the HDM.R. elataandR. palustrisrepresent a tractable experimental system to examine adaptive responses of polyploids during colonization in the HDM. First,R. elataandR. palustrisare the only twoRorippaspecies 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 otherRorippaspecies in Eurasia (Fig. S1). Third,R. elataplants 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.palustrisplants 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 forR. elataandR. 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.

Fig.2. Trait variation among Rorippa species in East Asia.

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 346Rorippa elataindividuals 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 ofR. palustrisplus relatedRorippaspecies(R.globosa,R.dubia,R.indica,and the newly identifiedR. hengduanshanensis(Zheng et al., 2021)), as well asNasturtium officinaleandBarbarea intermediaas outgroups(Table S1;Fig.S3).Given the selfing character ofR.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 ofR. elata, and 21 and 17 populations ofR. palustriswere located in the southern and northern HDM, respectively.

To investigate the polyploid speciation histories ofRorippa elataandR. palustris, 10 samples from sevenRorippadiploids(2n= 2x= 16), covering nearly all continents were selected as subgenomic representatives (Bleeker et al., 2002; Table S1),includingRorippa austriaca(2 samples),R. cantoniensis(1),R. globosa(1), andR. islandica(2) from Eurasia; andR. curvisiliqua(1),R. sessiliflora(1), andR. sinuata(2) from North America. Two samples ofNasturtium officinalewere used as outgroups.

2.2. Flow cytometry

Genome sizes of germinated seedling samples were estimated by flow cytometry (FCM) as described previously (Doleˇzel et al.,2007), using a BD FACSVerse Flow Cytometer (USA) withCapsella 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 diploidRorippa 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 186R.elatasamples from 92 populations and 62R. palustrissamples 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, andtrnL-trnF, were tested and filtered for amplification efficiency and nucleotide polymorphism across 72 randomly selectedRorippa elatasamples.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 lowcopy nuclear loci among 14 candidates, includingACO1,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)forRorippa elataandR.palustrisaccording 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 singlegene trees, subgenomic origins of PURC-inferred homoeologs inR. elataorR. palustriswere classified as g1 and g2 for each locus,according to their phylogenetic positions within the backbone built for seven diploidRorippaspecies (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 polyploidR.elataorR.palustriswere 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 amongRorippahaplotypes, 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 72R.elatasamples used during marker screening tests(Fig. S7).A dated phylogenetic tree was generated forRorippaplastid haplotypes and outgroups(Barbarea intermediaandNasturtium officinalein 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 diversityHdand 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 ofR. elataandR. palustrisusing MaxEnt v.3.4.3(Phillips et al.,2006).Five bioclimatic variables with low correlations (absolute Pearson's correlation coefficientR＜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)forR.elataand the HDMR.palustriswere 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 ofR. elataandR. palustrisacross 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 DIVAGIS 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 forRorippa 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).

ForRorippa 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°C and 60% humidity. Plants were randomized into 5 × 5 plots (75 × 75 cm2) 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 Augustof 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 primarysprouts(PS)orfruiting primarysprouts(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 localRorippapolyploid species (R. elataandR. palustris) were better adapted to the environment of the HDM,we quantified germination rates forR. elata,R. palustris, and other lowland diploids or polyploids in the SABG common garden.EightRorippaspecies were sampled, including three diploids,Rorippa barbareifolia,R. cantoniensis, andR. globosa; three tetraploids,R. dubia,R. elata, andR. palustris; and two hexaploids,R. hengduanshanensisandR. indica. For each species,10 accessions(except seven forR. globosaand nine forR. 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-1gibberellin in the dark at 4°C 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 selectedRorippa palustrisaccessions from the HDM. BecauseR. elatacannot flower even under long-term vernalization treatment (e.g., 2 months), we did not include these plants in the growth chamber experiment. Seeds were surfacesterilized and germinated on half-strength MS media and then stratified in the dark at 4°C for 7 days.After 10 days of germination under LD conditions at 20°C,seedlings were transplanted into soil(nutrient soil to vermiculite = 1:3) and grown in the growth chamber with LD conditions, at 20°C and 60% humidity. For the vernalization treatment, 7-day-old seedlings were moved to 4°C and short-day(SD)conditions(8 hours light and 16 hours dark)for 4 weeks and then grown under 20°C 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 coefficientRbetween 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 amongRorippaspecies, with Bonferroni adjustment ofP-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 theRBook (Crawley,2013).

3. Results

3.1. Speciation mode

The estimated mean (±SD) genome sizes were 562.431±11.219 Mb forRorippa elata(Fig.S4;range,541-587 Mb)and 545.750 ± 10.779 Mb (514-568 Mb) forR. palustris. In comparison with diploidR.islandica(245 Mb;2n=2x=16),samples of bothR. elataandR. palustriswere identified as tetraploids with 2n= 4x= 32, with further verification by chromosome counting(Fig. S4).

We inferred the hybrid speciation history for tetraploidRorippa elataandR. 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 ofR. elataorR. palustrisclustered with diploidR.globosaacross all gene trees and the other homoeolog ofR.elatadiverged from all otherRorippaspecies, except atMLH1.The other homoeolog ofR. palustrisclustered with diploidR. 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,bothR.elataandR. palustrisshowed hybrid origins within the Eurasian clade(Fig. 3c-d). Based on the phylogenetic relationships, one ancestral progenitor ofR.elataorR.palustriswas inferred as closely related toR. globosaand the other progenitor was identified as an unknown species, possibly extinct or unsampled forR. elataand closely related toR. islandicaforR. palustris.

Fig.3. Speciation modeling.

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 forRorippaand related plants when including indels(Figs.4,S9 and S10).Seven plastid haplotypes and two ITS genotypes were unique toR. elata, and one plastid haplotype and one ITS genotype were unique toR. palustrisin 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. palustrishad a higher overall plastid haplotypic and nucleotide diversity (Hd= 0.257,π=0.00080)than didR.elata(0.158 and 0.00052;Table S6).At the regional scale, populations ofR. palustrisfrom the southern HDM harbored more genetic variation(Hd=0.325,π=0.00131)than did populations from the northern HDM (0.166 and 0.00012). ForR.elata,genetic diversity was similar in the southern and northern HDM populations (Hd= 0.162 vs.0.155,π = 0.00054 vs. 0.00051).

Fig.4. Haplotype network,phylogeny,and distribution of Rorippa species in the Hengduan Mountains.

3.3. Haplotype network and phylogeny

Patterns of genetic variation were inferred through the construction of a haplotypic or genotypic network and phylogeny acrossRorippapopulations in the HDM. The plastid haplotypic network showed two main genetic clades, including one clade dominated byR. elata(named theR. elataclade) and another dominated byR. palustrisor otherRorippaspecies (named theR. palustrisclade; Fig. 4a). The core haplotype of theR. elataclade,labeled Hap_2, was shared by 84% ofR. elataindividuals. The core haplotype of theR. palustrisclade, Hap_6, was shared by 78% ofR.palustrisindividuals and 16R.elataindividuals.Six haplotypes ofR. elata(Hap_3, 6, 11, 13, 14, and 16) and two haplotypes ofR. palustris(Hap_2 and 12) were not within these main clades(estimated by pairwisep-distances; Fig. S9). These haplotypes accounted for only 6% ofR. elataindividuals or 17% ofR. palustrisindividuals,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 forR. elataor Hap_2 forR. palustris.

The two genetic clades were also well-supported(posterior = 0.99) in the dated BEAST tree (Fig. 4a). TheR. elataclade andR. palustrisclade originatedc. 0.55 Ma (95% HPD:1.30-0.14 Ma) andc. 0.86 Ma (1.94-0.22 Ma), respectively, with a crown age ofc.1.55 Ma (3.25-0.51 Ma). The core clade ofR. elata(excluding Hap_5,7,and 8)covered 92%ofR.elatasamples,with a crown age ofc.0.37 Ma (0.89-0.04 Ma).

3.4. Phylogeography and niche dynamics

Ancestral geographic ranges were reconstructed forRorippaplants in the HDM along the well-supported plastid haplotypic tree(Fig. 4a). The southern HDM was the most likely ancestral area for both theR. elataclade andR. palustrisclade in a BioGeoBEARS analysis.Divergence between theR.elataandR.palustrisclades was mediated by an ancient dispersal event from the southern to southern or northern HDM.Multiple dispersal(two forR.elataand four forR. palustris) or vicariance events (one forR. elataand two forR. 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). ForR. 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. ForR. palustris, the core haplotype (Hap_6) was mainly found in the northern HDM (Fig. 4b), whereas most of the relict haplotypes ofR.palustriswere located in a relict region that fully overlapped with the relict region forR. elata. Based on the above phylogeographic results, the southern HDM would be a stable ancestral area for eitherR. elataor the HDMR. 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 forR.elataorR.palustrispopulations 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 ofR. elataandR. palustrisexpanded 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 eitherR.elataorR.palustris(Figs.4b and S10b).Geologically,niche ranges during the LGM were mainly limited by the reconstructed fallensnowline 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 forR. elataorR. palustrisplants. 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 ofR. elataorR. palustrisin response to glacial dynamics in the HDM.

3.5. Adaptive response

AmongRorippaspecies in Eurasia,plants ofR.elatainhabited areas with the highest altitudes (mean elevation = 3589 ± 612.096 m,P＜ 0.0001; Fig. S1), together with the HDMR. palustris(3192±823.859 m in HDM vs.387±67.171 m in other places).AmongRorippaspecies found in sympatry in East Asia, tetraploids had significantly higher 18DASgermination rates(LSDtest,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, HDMR. palustrisshowed the highest 18DAS germination rate(lsmean=4.73;Fig.2f),followed byR.elata(3.90),R.dubia(3.25),R.indica(2.50),R.globosa(2.45),R.barbareifolia(2.03),R. cantoniensis(1.72), andR. hengduanshanensis(1.17). The germination rate of HDMR.palustrisdiffered significantly from those of other species(df=66,P≤0.030)exceptR.elata(df=66,P=1.000).Accessions ofR.elatahad significantly higher germination rates than those of otherRorippaspecies(df=66,P≤0.002)exceptR.palustris,R. dubia,R. indica, andR. globosa(df = 66,P≥0.069). The 2-year common garden experiment further revealed clinal variation in multiple FRTs inR.elataalong source latitudes or altitudes(Table 1;Fig.S15).

Table 1 Estimates of the effects of source latitude or altitude on trait variation.

3.5.1. Overallfitness

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 ofRorippa 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 coefficientR=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.

Fig. 5. Niche evolution of Rorippa elata and R. palustris in the Hengduan Mountains.

3.5.2. Survival

More than 91% ofRorippa elataplants 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. elatapopulations in the highlatitude 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 ofRorippa elatadiffered 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 ofR≤0.10,P≥0.610) or the total number of secondary inflorescences (TN_SI; absolute value ofR≤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 ofR. elataat 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 forRorippa 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 ofR≤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 ofR. elatashowed early flowering and produced fewer ovules compared with populations at low latitudes or altitudes.

3.5.5. Flowering time

ForRorippa 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 nonmutually 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 thepolyploid 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 postglacial colonization of HDM from the south to north by two allotetraploids,R. elataandR. 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 bothRorippa elataand HDMR.palustris(Figs.3,S5 and S6) during the mid-Pleistocene aroundc. 0.55 Ma andc.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 asPicea purpurea(1.30 Ma)(Sun et al.,2014)orOstryopsis 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 allopolyploidAegilopswheats 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 allotetraploidCapsella bursa-pastoris(Han et al., 2015)andCardamineflexuosa(Mand′akov′a et al., 2014), as well as the allohexaploidEchinochloa 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 ofRorippa elataandR. 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 ofRorippapolyploids in the HDM, with recurrent northward expansion during interglacial periods and repeated southward withdrawals during glacial periods (Figs. 4,5 and S11-S14).

The initial expansion of the core clade ofRorippa elata(covered 92%ofR.elatasamples,with the crown age ofc.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 ofR. elatamay 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 ofR. elataorR. palustriswere 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 ofR.elataandR.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 ofRorippa elataorR. palustrisin 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′e and Dangles, 2019), as well as the reduced summer precipitation driven by weak monsoons sincec. 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 likeR. elataorR. palustristo 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 supportingthisconjectureare scarce (Ramsey, 2011; but seeertner etal., 2019).In the HDM, speciesrichness 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 ofRorippa elataandR. palustris(Figs. 2, S1 and S2). How theseRorippaallotetraploids have adapted better than diploids across the HDM is an intriguing question that awaits further investigation. Our results agree with thepolyploid adaptation hypothesis(Fig. 1) (Hegarty and Hiscock,2008; Van de Peer et al., 2017), which stipulates that successfulRorippapolyploids may spread via adaptive responses to the environments in the HDM.

We tested the hypothesis inRorippapolyploids through common garden and growth chamber experiments. Seedling survival rates were higher forRorippapolyploids than for diploids (Fig. 2),implying an evolutionary advantage under the high-altitude environments of the HDM.ForR.elataandR.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.elataplants 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 annualR.palustrisfrom the northern HDM highlands showed early flowering and low vernalization sensitivity(Figs.S17 and S18), consistent with findings in the annual model plantArabidopsis 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¨orner, 2020). In contrast to lowland plants with more diffuse architectures, the highlandRorippa elataevolved 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 inR. elatapopulations.Along latitudes or altitudes,survival rates increased as ovule numbers decreased(Fig.S15),perhaps reflecting a generalKselection 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 ofRorippapolyploids 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 inRorippapolyploids 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 forRorippapolyploids in the HDM.

5. Conclusion

We found thatRorippa elata,which is endemic to the HDM,and the widely distributedR. 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 thepolyploid 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 fitnessrelated 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 providingCapsella rubellaseeds 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.