Yihui Wang
Yuan Chen3†
Zeljana Prijic
Tatjana Markovic
Yingmin Lyu
Caihuan Tian
Xiuxin Zhang- 1Beijing Key Laboratory of Ornamental Germplasm Innovation and Molecular Breeding, China National Engineering Research Center for Floriculture, College of Landscape Architecture, Beijing Forestry University, Beijing, China
- 2State Key Laboratory of Vegetable Biobreeding, Key Laboratory of Biology and Genetic Improvement of Flower Crops (North China), Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Ministry of Agriculture and Rural Affairs, Beijing, China
- 3Chongqing Academy of Agricultural Sciences, Chongqing, China
- 4Institute for Medicinal Plants Research “Dr Josif Pančić”, Belgrade, Serbia
Ecological changes have been observed to promote rates of lineage diversification, yet the precise roles of ecological factors, species evolution, and environmental variability in driving species diversity remain research hot spots. The association between ecological change and lineage diversification, particularly with regard to the size of the time scale, remains poorly understood. To explore whether ecological change facilitates species evolution, we focused on the unique family Paeoniaceae, which encompasses both herbaceous and woody taxa, to investigate the evolutionary rates. As a unique family characterized by a single genus of angiosperms and comprising various climatic types, the ecological niche changes of Paeoniaceae are closely associated with the evolution, making it an ideal model for conducting association analysis. In this study, we integrated the molecular fragments and ecological factors to explore the relationship between species evolution and niche changes in Paeoniaceae. The phylogenetic tree revealed that Paeoniaceae forms a sister relationship with Penthoraceae, Haloragidaceae, Iteaceae, Crassulaceae, and Saxifragaceae, constituting an independent clade based on the positive selection of molecular fragments including two protein-coding genes and eight non-coding regions. The divergence time was estimated to be between 102 and 116 Mya (Million years ago). The phylogenetic tree within Paeonia revealed a clear division into three groups: sections of Paeonia, Moutan, and Onaepia with high support values for each branch based on the ten positive selection of molecular fragments. The rapid rate of evolution observed in Paeonia, about 0-5 Mya. In addition, ecological niche modeling showed that the potential distributions for Paeonia expanded from middle Asia to eastern Asia, and from central North America to the Northern part of North America during the Last Glacial Maximum (LGM) to Mid Holocene (MID) period. This suggests that Paeonia continuously adapted to changing ecological environments over time. Compared to the rate of climatic niche divergence and lineage diversification, the ecological niche of Paeonia underwent significant changes during the period of 3-11 Mya, occurring 5 Mya earlier than the period of evolutionary rate change. These findings offer comprehensive insights into the relationship between niche change and the evolution of species, providing valuable perspectives for further ecological cultivation efforts.
1 Introduction
Evolution is a long-term process influencing species diversity, regulated at a macroscopic level by climate change and at a microscopic level by nucleotide variation (Pyron and Burbrink, 2013). For microevolution, the inferences regarding adaptive evolution within DNA sequences serve as research hot spots for biologists (Chen et al., 2021). The divergence time of species, as determined by plastid super-barcodes such as ndhH, rpl23, and rbcL, is extensively utilized (Krawczyk et al., 2018; Wu et al., 2021a). In recent years, with advancements in sequencing technology, the utilization of the entire chloroplast (cp) genome has emerged as a versatile tool for phylogenetics (Dong et al., 2018; Zhai et al., 2019; Raman and Park, 2016). Indeed, some studies have shown that the recognition of pervasive adaptive evolution, based on high levels of genetic diversity, relies on positive selection of molecular fragments (Shapiro et al., 2007; Yang et al., 2005). Signals of positive selection could accurately assess the divergence time between species (Tajima, 1989; Fay and Wu, 2000; Sabeti et al., 2002; Chen et al., 2021). It has previously been reported that CemA, rps7, and rpl23 genes exhibit relatively higher nucleotide substitution rates than other genes within Poaceae and Typhaceae (Guisinger et al., 2010). Similar patterns were observed within Saxifragales. In addition, ycf1 and ycf2 emerged as the predominant genes with high nucleotide substitution rates in many plants, including Caryophyllaceae (Sloan et al., 2012), Geraniaceae (Guisinger et al., 2011) and Prunus (Xue et al., 2019). Given their extensive coding lengths, ycf1 and ycf2 are prone to frequent substitution rate variations (Neubig et al., 2009; Huang et al., 2010; Dong et al., 2015).
For macroevolution, aspects of climate change influence the distribution of angiosperms globally and the clade of species for evolutionary processes (Haffer, 1997; Barnagaud et al., 2012; Ali and Aitchison, 2014). Studies indicate that climate and niche change had a profound impact on the diversity of plant species in Eastern Eurasia since the Last Glacial Maximum (Song et al., 2024). Additionally, the global diversity of angiosperms has undergone significant changes due to Quaternary climate and niche fluctuations (Xu et al., 2023). Clearly, climate changes are closely related to species distribution. However, the impact of natural selection on species unfolds over extended periods, leaving the question of whether ecological niches change during species evolution still open for debate (Laland et al., 1996; Laland and Feldman, 1999; Wiens and Graham, 2005; Losos, 2008). In recent years, studies on niche changes and species evolution have been increased. Broadly speaking, the periods of rapid diversification in angiosperms are similar to the periods of rapid changes in average annual temperature and precipitation (Liu et al., 2020). The average diversification rate of angiosperm niches is comparable to the rapid divergence time of the Crassulaceae based on 301 protein-coding loci (Folk et al., 2019). Therefore, there is a universal need for studies investigating the relationship between changes in ecological niches and the evolutionary history of species.
Paeoniaceae Raf. is a family of angiosperms, comprising a single genus and 34 species, including both herbaceous and woody plants found across Asia, Europe, and North America (Hong, 2010, 2021). This taxon exhibits diverse habitats in temperate conditions, including Monsoon Climate, Temperate Continental Climate, Temperate Marine Climate, and Plateau Climate. Its distribution spans three continents, North America, Asia, and Europe (Hong, 2010, 2021), indicating a rich diversity of ecological environments. The Paeoniaceae has been classified within Saxifragales based on molecular data (Dong et al., 2018; Folk et al., 2019; Li et al., 2019, 2021), diverging from its previous placement within Ranunculidae based on morphological characteristics (Guan, 1979). The bootstrap and posterior probability values of the Paeoniaceae clade are relatively low due to the limited diversity of genetic fragments (Li et al., 2019, 2021). In comparison with the findings in the literature, there exists a significant disparity in the estimated divergence times of Paeoniaceae based on chloroplast genomes or molecular fragments (Li et al., 2019; Zhou et al., 2021). Therefore, the selection of appropriate molecular fragments is vital for determining both the systematic placement of Paeoniaceae within Saxifragales and the accurate estimation of its divergence time.
Among the 34 species of Paeoniaceae, 15 species are distributed in Asia and North America (Hong, 2010, 2021). Despite efforts to estimate the divergence times of Paeonia in Asia based on fossils from other families within Saxifragales (Zhou et al., 2021; Chen et al., 2023), the phylogenetic position of Paeoniaceae within Saxifragales remains uncertain. To explore whether ecological changes drive species evolution, it is urgent to study the relationship between ecological niche change and divergence times of Paeonia. However, only the ecological niche of Paeonia rockii has been studied thus far (Dong et al., 2022).
In this study, we adopted an approach that utilizes a positive selection of molecular fragments from whole chloroplast genome sequences to comprehensively infer phylogenetic relationships and divergence times within Paeoniaceae and other related families. Besides, we explored the change in ecological niches and the evolutionary history of Paeonia in Asia to elucidate whether niche changes facilitate species evolution. Our study offers novel insights into niche dynamics within a family encompassing species evolution.
2 Materials and methods
2.1 Sequence data
We obtained sequence data for 294 species within Saxifragales and outgroups from the National Center for Biotechnology Information (NCBI: https://www.ncbi.nlm.nih.gov/). These data are detailed in Supplementary Table S1. The dataset included species from the following families: Paeoniaceae (29), Hamamelidaceae (38), Altingiaceae (10), Cercidiphyllaceae (2), Daphniphyllaceae (5), Saxifragaceae (108), Grossulariaceae (5), Iteaceae (1), Crassulaceae (79), Haloragidaceae (7), Penthoraceae (1), and outgroups (9). The Genbank numbers sourced from NCBI are listed in Supplementary Table S1.
2.2 Positive-selection signals screening and genome comparison
To estimate rates of nucleotide substitution, we calculated the ratio of the number of nonsynonymous substitutions per nonsynonymous site (Ka) to the number of synonymous substitutions per synonymous site (Ks) using the YN model in KaKs_Calculator (Zhang, 2022) for the protein-coding genes (PCGs). We employed the genetic code with 11 bacterial, archaeal, and plant plastid codes (Yang and Nielsen, 2000) for the species obtained from NCBI (Supplementary Table S2). The ratio of non-coding nucleotide substitution rate (Kn) to neutral substitution rate (Ks) was determined for the non-coding sequences by comparing aligned adjacent coding sequences. Notably, the genetic code used remains the 11-bacterial, archaeal, and plant plastid codes (Yang and Nielsen, 2000). Subsequently, we generated plots of Ka/Ks and Kn/Ks using the R-packages ‘tidyverse’ and ‘readxl’ (Wickham and Lionel, 2019; Hadley and Bryan, 2019). To investigate plastid genome divergence, we computed the nucleotide diversity (Pi) for the 74 common protein-coding genes and 72 common non-coding regions using DnaSP v5 (Librado and Rozas, 2009). Next, we identified hotspots characterized by high numbers of nucleotide substitutions and sequence distance (Pi) and then generated Pi plots using R-packages ‘tidyverse’ and ‘readxl’ (Wickham and Lionel, 2019; Hadley and Bryan, 2019).
2.3 Phylogenetic analysis within Saxifragales
In order to infer phylogenetic relationships within Paeoniaceae and other related families, we obtained 294 high-quality chloroplast genome data and comprehensive annotation files from NCBI (Supplementary Table S1). Given that the important lineages within Saxifragales consist of the clade of Saxifragaceae, Penthoraceae, Haloragidaceae, Crassulaceae, Iteaceae, and Grossulariaceae, we selected 201 representative species. As the focus of this study, we selected all 29 published data within Paeoniaceae, and selected 64 species to represent the other clade within Saxifragales and outgroups. Based on the 294 sequence data, we utilized MAFFT version 7.49 to compare the sequence fragments matrix (including ycf1, ycf2, ndhD-psaC, trnC-petN, trnK-rps16, ycf4-cemA, trnH-psbA, petG-trnW, petN-psbM, and petA-psbJ) using the ‘–auto’ strategy and normal alignment mode (Katoh and Standley, 2013). Ambiguously aligned fragments were subsequently removed using Gblocks version 0.91b with the following parameter settings: minimum number of sequences for a conserved/flank position (10/10), maximum number of contiguous non-conserved positions (8), minimum length of a block (10), and allowed half of the gap positions (Talavera and Castresana, 2007). For the species tree, a maximum likelihood (ML) phylogeny for the genera was inferred using RAxML v8.0.0 (Stamatakis et al., 2008). The best-fitting DNA substitution model for the ML tree was determined to be GTR (General Time Reversible) + F (Felsenstein) + I (proportion of Invariable sites) using ModelFinder (Kalyaanamoorthy et al., 2017). A standard bootstrap with 5000 replicates is used to infer the ML tree. In addition, Bayesian inference (BI) analysis was conducted using MrBayes version 3.2.6 (Ronquist and Huelsenbeck, 2003), with the best-fitting DNA substitution model being GTR + F + I. Markov chain Monte Carlo simulations (MCMC) were run for 10,000,000 generations. The BI analysis commenced with a random tree and sampled trees every 1,000 generations. The first 25% of the trees were discarded as burn-in, and the remaining trees were used to generate a majority-rule consensus tree. Subsequently, the trees were edited using ITOL to display bootstrap (BS) and posterior probability (PP) values (https://itol.embl.de/; Letunic and Bork, 2021).
2.4 Molecular dating within Saxifragales
In BEAST 1.10.4, a lognormal distribution with an uncorrelated relaxed clock model was run using the GTR + F + I site model, with a random starting tree and a Yule speciation process tree prior (Yule, 1927; Suchard et al., 2018). MCMC was performed for 500 million generations, with samples taken every 50,000 generations. The effective sample size (ESS) values were verified to exceed 200 for all parameters. The initial 25% of the trees were designated as burn-in and discarded. The remaining trees were then used to generate an all-compatible consensus tree using TreeAnnotator software (Drummond and Rambaut, 2007). Moreover, 95% highest posterior density intervals (HPD) for each node are shown on the tree. Besides, the phylogeny was calibrated using five fossils. The first calibration point was set to the age of Saxifragaceae (81 Mya with HPD about 72.1-89.8 Mya; Bell, 1957). The second calibration point was the age of Hamamelidaceae (77 Mya with HPD of about 70.6-83.5 Mya; Itterbeeck et al., 2007). The third calibration point was the age of Altingiaceae (91 Mya with HPD of about 89.3-93.5 Mya; Zhou et al., 2001). The fourth calibration point was the age of Cercidiphyllaceae (68.3 Mya with HPD of about 66-70.6 Mya; Brown, 1962). The last calibration point was the age of Vitales (123 Mya with an HPD of about 101.8-144.1 Mya; Li et al., 2019). The tree was viewed and edited using ITOL and FigTree version 1.4.0 (https://itol.embl.de/; http://tree.bio.ed.ac.uk/software/Figtree/; Letunic and Bork, 2021).
2.5 Phylogenetic analysis and molecular dating within Paeoniaceae
The phylogenetic analysis within the Paeoniaceae family was conducted using 15 wild species as listed in Supplementary Table S1. The analysis was based on a matrix comprising 10 sequence fragments and was performed using both the ML and BI methods utilizing RAxML v8.0.0 and MrBayes version 3.2.6 (Ronquist and Huelsenbeck, 2003; Stamatakis et al., 2008). Among the 15 species, Chrysosplenium album (NC067019), Heuchera abramsii (MN496062), and Saxifraga sikkimensis (NC070509) were designated as outgroups. The parameters used for MAFFT version 7.49, Gblocks version 0.91b, RAxML v8.0.0, and MrBayes version 3.2.6 were consistent with those employed for the phylogenetic analysis within Saxifragales (Ronquist and Huelsenbeck, 2003; Talavera and Castresana, 2007; Stamatakis et al., 2008; Katoh and Standley, 2013). The best-fitting DNA substitution model for both methods was determined to be GTR + F + I. Given the absence of reliable fossils within Paeoniaceae, we initially used the estimated time frame of 106-112 Mya, as previously described. Next, we used the 70.6-83.5 Mya divergence as the calibration point for Saxifragaceae (Bell, 1957). The parameters employed in BEAST 1.10.4 and TreeAnnotator software were consistent with those utilized in the phylogenetic analysis within Saxifragales (Drummond and Rambaut, 2007; Suchard et al., 2018). Finally, the ‘ape’ R package was used to calculate and plot the evolutionary rate of Paeonia (Paradis and Schliep, 2019).
2.6 Geographic distributions of species
We acquired distribution data of the 15 species distributed in Asia and North America listed in Supplementary Table S1 of Paeoniaceae from various sources, including the Chinese Virtual Herbarium (http://www.cvh.org.cn/), Global Biodiversity Information Facility database (http://www.gbif.org/), Flora of China (Guan, 1979), and records from field surveys conducted by our research group. Each data about latitude and longitude information for each data point was meticulously verified, and specimens lacking such coordinates were reconstructed based on microclimate data. A total of 9,702 record points were screened and included in the analysis (Supplementary Table S3). Subsequently, we filtered and duplicated the data with similar latitude and longitude, ensuring a minimum distance of less than 5 kilometers between each point, and visualized the distribution of the dataset using the R-packages ‘dplyr’ and ‘geosphere’ (Karney, 2013; Morgan and Cheng, 2023).
2.7 Climate data screening and filtering
We obtained climate data for two historical periods comprising 19 bioclimatic and elevation layers for the last glacial maximum (LGM, about 22,000 years) and Mid Holocene (MID, about 6000 years ago) in the WorldClim global climate database (https://www.worldclim.org/). The LGM and MID data were obtained from the CCSM4 climate model. To ensure greater accuracy about the variables, we performed the correlation analysis between individual environmental factors and species with the method of spearman and remove environment variables with low relevance (Murdoch and Chow, 1996).
2.8 Analysis of suitable distribution areas
The collected data, comprising distribution points and the precise variables without bio8, bio9, bio18 and bio19, were imported into MaxEnt 3.4.1 software for modeling analysis with the best model using the R package “kuenm” (Phillips et al., 2017; Muscarella et al., 2014). The 75% of the distribution points were treated as the training data, while the remaining 25% were designated as random test data. We optimized the 29 feature combinations (FC) and regularization multipliers (RM) from the “kuenm” package (ranging from 0 to 4.0 with an interval of 0.1), resulting in the generation of 2,480 candidate models. The model with the smallest delta AICc was selected as the optimal solution and used to establish the final model. Next, the final model (linear features and quadratic features) was replicated 10 times using cross-validation to assess accuracy and the output data format was logistic (Moreno et al., 2011). A regularization multiplier of 0.1 was applied, and the maximum number of background points was set to 10,000. Besides, the number of iterations for each dataset in MaxEnt maximum was set to 500. The convergence threshold for the predictions was established at 0.0001, and the probability of presence at ordinary occurrence points was set to 0.5 (Elith et al., 2011). The criterion for assessing the accuracy of the model’s predictions is based on the area under the curve (AUC) value of the Receiver Operator Characteristic (ROC) (Hebbar et al., 2022). AUC >0.7 is typically considered indicative of a good model (Rebelo et al., 2010). The determination of the main bioclimatic factors influencing the suitable distribution areas of Paeonia in Asia was conducted using the method of Jackknifing. Finally, we classified suitable areas into the following four categories using natural breaks (Jenks): non-suitability (p < 0.08), low suitability (0.08 ≤ p < 0.26), moderate suitability (0.26 ≤p < 0.43), and high suitability (0.43 ≤ p < 1) during the LGM period (Saleh, 2020). During the MID period, suitable areas were classified into the following four categories: non-suitability (p <0.09), low suitability (0.09≤p < 0.27), moderate suitability (0.27 ≤p < 0.46), and high suitability (0.46≤p < 1).
To explore the change in distributions between LGM and MID periods, we used the SDM_Toolbox_v2.4 package to visualize the potential variations in suitable areas (Brown et al., 2017). Initially, the potential suitable areas were converted into binary vectors, with suitability probabilities of species set at 0.43 and 0.46, respectively.
2.9 Rate of climatic niche divergence
To estimate the rate of climatic niche divergence in Paeonia, we analyzed the similarity of climatic niche spaces from each clade using the R package ‘ecospat’ (Broennimann et al., 2011; Cola et al., 2017). For each clade and the overall environmental niches from the phylogenetic tree, we compared the similarity of the environmental niches using Schoener’s D niche metric ranges, where D-Values range from 0 to 1 (Warren et al., 2008). Principal component analysis (PCA) was used to determine the proportion of environmental factors contributing to the ecological niche (Broennimann et al., 2011). Finally, we compared the rate of evolution with the rate of climatic niche divergence in Paeonia using the R package ‘ggplot2’.
3 Results
3.1 Nucleotide substitution rates
To trace the evolutionary rates along Paeonia, we analyzed nucleotide substitution rates relative to neutral substitution rates on sequences, along with Pi analysis. We screened for common high nucleotide substitution rates of sequence fragments based on methods established in previous research (Zhou et al., 2021; Chen et al., 2023). First, Nonsynonymous (Ka) and Synonymous (Ks) substitution rates, along with Ka/Ks, were estimated for the 74 common protein-coding genes to detect evolutionary rate heterogeneity. Among the 74 genes analyzed, large ribosomal protein (rpl23), hypothetical chloroplast reading frames (ycf2), hypothetical chloroplast reading frames (ycf1), large ribosomal protein (rpl22), and small ribosomal protein (rps18) exhibited relatively higher Ka/Ks values (Ka/Ks > 0.6), suggesting they have been subjected to positive selection pressures. The remaining genes exhibited relatively lower Ka/Ks values (Ka/Ks < 0.6), suggesting their exposure to purifying selection pressure (Figure 1A; Supplementary Figure S1). We also compared Kn/Ks across non-coding regions to identify variations in evolutionary rate. Among 72 non-coding regions, we selected those with high Kn/Ks values across 11 species. Elevated Kn/Ks values (Kn/Ks > 1) had been identified within the following non-coding regions: ndhD-psaC, trnC-petN, trnK-rps16, ycf4-cemA, trnH-psbA, petG-trnW, petN-psbM, petA-psbJ, psaI-ycf4, psbE-petL, psbI-trnS, and psbJ-psbL (Figure 1B; Supplementary Figure S1).
Figure 1. The map of Ka (Kn)/Ks and Pi analysis within 11 species. (A) The estimations of Ka to Ks of plastid protein-coding genes (PCG). The estimations of nonsynonymous (dN), synonymous (dS) substitution rates, and dN/dS of plastid protein-coding genes (PCG). (B) The estimations of Kn/Ks of plastid genes. The estimations of nonsynonymous (dN), synonymous (dS) substitution rates, and dN/dS of plastid noncoding regions.
The Pi analysis identified several chloroplast regions (atpF-atpH, atpI-rps2, ccsA-ndhD, ndhD-psaC, petA-psbJ, petG-trnW, petN-psbM, psbZ-trnG, rbcL-accD, rpl33-rpl18, rps15-ycf1, rps16-trnQ, trnC-petN, trnD-trnY, trnG-trnR, trnH-psbA, trnK-rps16, trnL-trnF, trnS-trnG, trnY-trnE, and ycf4-cemA) and genes (ccsA, matK. ndhE, ndhF, psaI, rps15, ycf1, ycf2, and ycf3) to be divergent hotspot regions (Figure 1; Supplementary Figure S2). In comparison to Pi analysis, Ka/Ks and Kn/Ks, the commonly identified potential molecular markers for species identification and phylogenetic evolution within Saxifragales include the following 10 regions: ndhD-psaC, trnC-petN, trnK-rps16, ycf4-cemA, trnH-psbA, petG-trnW, petN-psbM, petA-psbJ, ycf1, and ycf2.
3.2 Phylogenetic relationships and divergence time estimation within Saxifragales
We conducted a phylogenetic analysis to determine the systematic classification of Paeoniaceae within Saxifragales. Based on data from 10 regions extracted from chloroplast genomes, our phylogenetic reconstructions showed a similar topological tree between the ML tree and BI tree within the family clades (Figure 2; Supplementary Figure S3). The nodes within each family exhibited high bootstrap values and posterior probabilities, indicating a reliable inference of phylogenetic relationships. The tree comprehensively resolved the phylogenetic relationships within Saxifragales among the major clades. Saxifragales underwent a divergence event at approximately111 Mya (HPD:102-117 Mya), giving rise to two large sub-clades: Clade I (comprising Hamamelidaceae, Altingiaceae, Cercidiphyllaceae, and Daphniphyllaceae) and Clade II (encompassing Paeoniaceae, Penthoraceae, Haloragidaceae, Iteaceae, Crassulaceae, and Saxifragaceae) with a PP of 1 and a BS of 100. Paeoniaceae diverged from Clade II at approximately 110 Mya (HPD:102-116 Mya) forming an independent clade with a PP of 0.95 and a BS of 71. The lineages, including Penthoraceae, Haloragidaceae, and Crassulaceae, initiated diversification at 109 Mya (HPD:102-116 Mya) from the Iteaceae, Grossulariaceae, and Saxifragaceae groups, with a PP of 1 and a BS of 96. Crassulaceae, as the sister group to Haloragidaceae and Penthoraceae, underwent diversification at 108 Mya (HPD:100-117 Mya) with a PP of 1 and a BS of 100. Haloragidaceae began diversifying from Penthoraceae at 90 Mya (HPD:56-115 Mya), with a PP of 1 and a BS of 100. Iteaceae commenced diversification from Grossulariaceae and Saxifragaceae groups at 107 Mya (HPD: 100-117 Mya), with a PP of 1 and a BS of 100. Furthermore, Saxifragaceae initiated diversification from Grossulariaceae at 106 Mya (HPD: 97-117 Mya), with a PP of 1 and a BS of 100. Within Clade I, Daphniphyllaceae diverged approximately 107 Mya (HPD: 99-117 Mya) as an independent clade from Cercidiphyllaceae, Altingiaceae, and Hamamelidaceae groups, with a PP of 1 and a BS of 100. Hamamelidaceae diverged from Altingiaceae and Hamamelidaceae groups at around 104 Mya (HPD: 87-115 Mya), with a PP of 0.95 and a BS of 70. Finally, it was found that Altingiaceae diverged from Hamamelidaceae at 77 Mya (HPD:70-83 Mya), with a PP of 0.95 and a BS of 100. It was evident that the majority of species within Saxifragales underwent divergence events between 100 and 115 Mya, during the lower Cretaceous period.
Figure 2. Phylogenetic and divergence time analysis within Saxifragales (A) Results of phylogenetic and divergence time analysis within Saxifragales. The 95% highest posterior density (HPD) credibility intervals for node ages are labelled above the line. The numbers on the clade represent the support value and posterior probability. (B) ML tree within Saxifragales. The numbers on the clade represent the support value.
3.3 Phylogenetic relationships and divergence time estimation within Paeoniaceae
To elucidate the evolutionary relationships within Paeoniaceae, we performed phylogenetic analyses of the genus Paeonia in Asia based on data from 10 regions extracted from chloroplast genomes. Notably, Heuchera abramsii, Chrysasplenium album, and Saxifraga sikkimensis were employed as outgroups (Figure 3A). The topological structures of the trees remained consistent across both ML and BI analysis methods. Paeonia species exhibited subdivision into three large sub-clades, corresponding to the taxonomic classification of sections of Paeonia, Moutan, and Onaepia (Hong, 2010, 2021). The Moutan section diverged from the herbaceous clades (sections of Paeonia and Onaepia) at 24 Mya (HPD:7-53 Mya), with a PP of 1 and a BS of 100. P. brownie, representing section Onaepia, formed a sister relationship with section Paeonia with a PP of 1 and BS of 100. This divergence occurred at approximately 18 Mya (HPD:3-36 Mya). The species in the Moutan section can clustered into two distinct clades: group one (P. delavayi and P. ludlowii) and group two (P. ostii, P. rockii, P. qiui, P. jishanensis, and P. decomposita), with a PP of 1 and a BS of 99. The divergence between the two clades occurred around 11 Mya (HPD:1-28 Mya). In section Paeonia, P. emodi was positioned at the base with a PP of 1 and BS of 100. Within this section, a subclade comprising P. lactiflora, P. anomala, and P. veitchii formed a sister relationship with another subclade containing P. daurica, P. obovata, and P. intermedia with a PP of 1 and a BS of 100. The divergence between these two subclades occurred approximately 7 Mya (HPD:1-21 Mya). Collectively, these results suggest that the rate of evolution of Paeonia rose rapidly from 0 to 5 Mya.
Figure 3. Phylogenetic and species evolution rate analysis within Paeonia. (A) Results of comparative analysis of niche changes and phylogenetic trees based on 10 regions of Paeonia. The 95% highest posterior density (HPD). Credibility intervals for node ages are labelled above the line. The numbers on the clade represent the support value and posterior probability. (B) Rate of species evolution (C) Rate of niche change. (D) Niche similarity tests.
3.4 Potential distributions and comparisons of Paeonia in Asia between LGM and MID
The single-factor correlation analysis reveals that the all variables without bio8, bio9, bio18 and bio19 exhibit significant correlations with the distribution of the species (Supplementary Figure S4). Based on the distribution patterns of Paeonia in Asia and North America (Supplementary Figure S5) and analysis of the potential distributions during the LGM and MID periods, the mean AUC for test replicates runs 0.838 for LGM and 0.843 for MID (Supplementary Figure S6). Plots illustrating the contributions of respective variables were generated from a pool of 16 factors influencing distribution patterns across the two periods (Supplementary Figure S7). The contributions of bio11, bio12 and bio5 to the LGM period were determined to be 27.2%, 22.2% and 19.5% respectively, with a cumulative contribution of 70%. On the other hand, the contributions of bio11, bio1, bio12 and elevation, to the MID period were found to be 34.7%, 17.1%, 12.6%, and 8.5%, respectively, with a cumulative contribution of 72.9%.
Based on the natural breaks method, the potential distribution of Paeonia in Asia and North America during the LGM period was classified into four grades: unsuitable area, lowly suitable area, moderately suitable area, and highly suitable area. The regions identified as highly suitable for Paeonia distribution include central China, Japan, the Korean Peninsula, regions in middle Asia such as Russia and Kazakhstan, and North America (Figure 4A). The respective areas for highly, moderately, and lowly suitable distribution areas were 1079.491×104 km2, 1658.386×104 km2, and 1792.825×104 km2, whereas the unsuitable area measured 10908.483×104 km2 (Supplementary Table S4). For the MID period, the potential distribution of Paeonia in Asia and North America was also categorized into the four grades. Specifically, the regions identified as highly suitable for Paeonia distribution were southern China, Japan, the Korean Peninsula, middle Asia such as Russia and Kazakhstan, and western North America. The respective areas for the highly, moderately, and lowly suitable distribution zones were 1118.525×104 km2, 2141.561×104 km2, and 2363.145×104 km2, while the unsuitable area measured 9811.784×104 km2 (Supplementary Table S4). A comparison of the LGM and MID periods revealed that the distribution of Paeonia expanded from central China to southern China while contracting from northern Russia to southern Russia. Moreover, the core area of Paeonia extended from North America to Canada (Figure 4B).
Figure 4. Distributions and comparisons of Paeonia in Asia between LGM and MID periods. (A) Potential distribution for Paeonia during LGM and MID periods using the MaxEnt model. Deep blue, blue, yellow, and red areas represent unsuitable, marginally suitable, moderately suitable, and highly suitable areas, respectively. (B) Spatial changes of Paeonia during LGM and MID periods. Gray, green, green, red, and purple areas represent unsuitable, unchanged suitable, expansion suitable, and contraction suitable areas, respectively.
3.5 Niche comparisons
PCA analysis was performed for niche comparisons. The findings revealed that the identified factors accounted for 99.31% of the 20 environmental factors in the distribution of Paeonia. PC1 accounted for 98.42%, with bio10 emerging as the main impact factor, whereas elevation and bio4 were found to be associated with PC2 (Supplementary Figure S8). The niche similarity, expressed as D-Value, was calculated for each main clade of the Paeonia tree, and the obtained results are displayed in Supplementary Table S5 and Figure 3B. The plotted change in D-Value for each branch showed a significant shift in the ecological niche of Paeonia around 3-11 Mya, predating the period of evolutionary rate change by 5 Mya. This suggests that alterations in the ecological environment promote the evolution of species, and the time scale is 5 Mya.
4 Discussion
4.1 The positive-selection signals from plastid genomes
By estimating nucleotide substitution rates across various genes and non-coding regions, a more comprehensive understanding of species evolution can be attained (Weng et al., 2016; Schwarz et al., 2017). First, we identified two high nucleotide substitution rates of genes, ycf1 and ycf2, in the coding regions, which is consistent with the findings of Dong et al. (2015). Therefore, leveraging the information from these ycf genes enables more detailed profiling of the angiosperm tree of life. Moreover, in comparison with the non-coding regions, a super-barcode region (trnH-psbA) was identified, as analogous to findings in Oryza (Zhang et al., 2021) and Dracaena (Zhang et al., 2019). In this study, additional non-coding regions (ndhD-psaC, trnC-petN, trnK-rps16, ycf4-cemA, petG-trnW, petN-psbM, and petA-psbJ) were newly identified as super-barcodes for Saxifragales. Specifically, trnK-rps16 was used to determine the systematic status of Paeonia in Pan-Himalayan regions (Zhou et al., 2021). Even though the regions exhibiting high nucleotide substitution rates varied when analyzed by KaKs, KnKs, and Pi methods, the commonly shared regions served as DNA super-barcodes, facilitating species delineation and the estimation of evolutionary history within Saxifragales.
4.2 Phylogenetic relationships and species divergence within Saxifragales
Phylogenetic analysis revealed that Paeoniaceae forms a sister relationship with the groups Penthoraceae, Haloragidaceae, Iteaceae, Crassulaceae, and Saxifragaceae. This finding is consistent with previous studies based on both coding and non-coding gene fragments (Dong et al., 2018; Zhou et al., 2021), but contradicts the results of a study focused on atpB, matK, and rbcL genes (Fishbein et al., 2001). Previous studies constructing phylogenetic trees to determine the position of Paeoniaceae using molecular fragments reported lower PP and BS. Our work defined the status of Paeoniaceae and improved the value of PP and BS. Divergence time estimates for Paeoniaceae vary significantly across studies due to the utilization of different molecular fragments. For instance, reported estimates include 80-100 Mya in the Upper Cretaceous (Tank et al., 2015; Folk et al., 2019; Dong et al., 2022), 65-90 Mya in the Upper Cretaceous (Zhou et al., 2021), and 88-130 Mya in the lower Cretaceous (Li et al., 2019). The divergence of Paeoniaceae is estimated to have occurred during the lower Cretaceous at approximately 110 Mya (HPD:102-116 Mya) based on the positive-selection regions, which represents a more accurate assessment compared to previous studies. In Saxiragales, the status of other families was consistent with the majority of studies (Dong et al., 2018; Li et al., 2019; Zhou et al., 2021) and the 95% HPD interval was narrower.
4.3 Phylogenetic relationships and species divergence within Paeoniaceae
The inclusion of three sections of Paeoniaceae (sections Paeonia, Moutan, and Onaepia) within a broader evolutionary lineage garnered high BS and PP. These findings were consistent with the previous phylogenetic reconstructions (Wu et al., 2021b; Zhou et al., 2021). Furthermore, our study enhanced the BS and PP of the three main clades. We also provided evidence supporting the proximity of section Paeonia to section Onaepia, a finding corroborated by Wu et al., 2021b and Zhou et al., 2021. In addition, we subdivided section Moutan and section Paeonia into two subsections, aligning with findings from previous studies (Zhou et al., 2021; Dong et al., 2022). Regarding the divergence time, the split between woody and herbaceous Paeonia species was approximately 24 Mya (HPD: 7-53 Mya). This estimate falls within the range previously determined based on whole chloroplast genome sequences, which was approximately between 20.78 Mya (HPD: 12-29 Mya; Dong et al., 2022) and 28 Mya (Zhou et al., 2021). The divergence of sections Paeonia and Onaepia occurred during the Miocene period, approximately 18 Mya, consistent with previous estimates (Zhou et al., 2021). Furthermore, the rapid rate of evolution observed in Paeonia, about 0-5 Mya, aligns with recent research (Chen et al., 2023).
4.4 The change of spatial distributions of Paeonia in Asia and North America
Our study suggests that Paeoniaceae originated during the lower Cretaceous, potentially influenced by the warm and humid conditions prevalent in the peripheral tropical environment during the Cretaceous era (Klages et al., 2020). During the LGM period, the regions identified as highly suitable for Paeonia were mainly concentrated in central China, Japan, the Korean Peninsula, middle Asia including Russia and Kazakhstan, and North America. The prediction results were consistent with the actual distributions (Dong et al., 2022). Given that during the LGM period much of the Earth was characterized by extensive glaciation, resulting in a reduction of warm-tropical and subtropical habitats compared to present conditions (Clark et al., 2009), the suitable distribution areas were more concentrated. After the LGM period, as temperatures gradually increased during the post-ice age, the Earth transitioned to being predominantly covered in warm-tropical and subtropical habitats. Cool-temperate biomes were primarily confined to polar regions (Pound et al., 2012). Consequently, the suitable distribution areas expanded from Northern Asia to Southern Asia. Similar predictions have been made for other rare and endangered plants, such as Cedrus and Firmiana kwangsiensis (Xiao et al., 2022; Gao et al., 2022). Besides, the dispersal of anthers by frugivorous birds, coupled with rapid radiation evolution occurring in the early Miocene, facilitated the swift expansion of the distribution range of Paeonia (Moyle et al., 2016). The observed spatial distribution patterns of Paeonia, coinciding with climate change, are consistent with the radiation evolution of birds.
The expansion of Paeonia into eastern Asia can be attributed to the intensification of the East Asian summer monsoon, which occurred around the Oligocene-Miocene boundary, and was triggered by the uplift of the northern Tibetan Plateau at around 20-25 Mya (Tada et al., 2016). This resulted in the enlargement of the potential highly suitable area of Paeonia in eastern Asia. The presence of endemic Paeonia species in Southwest China can be verified by evidence of population isolation and formation of endemic species in the region largely due to tectonic shifts and river course dynamics throughout the Quaternary period (Qiu et al., 2011). The widespread spatial distribution of Paeonia can be attributed to a significant transition from a zonal distribution of arid climates to warmer and wetter climates across central and eastern Asia since the early Miocene (Guo et al., 2008). Thus, it is plausible that the intensification of the East Asian summer monsoon during the late Miocene played a crucial role in triggering species or population-level diversification within tropical regions.
4.5 Niche differentiation and species divergence
Utilizing 20 ecological factors, our analysis revealed compelling evidence indicating a strong positive correlation between ongoing climatic variability and species divergence within Paeonia. This relationship is consistent with simulation-based theoretical findings suggesting a synergistic effect between phenotypic evolution and species differentiation (Hua and Wiens, 2013), where macroevolutionary dynamics and rates are driven by ecological and microevolutionary factors, with decreasing niche width accelerating speciation diversification (Aguilée et al., 2018). Both ecological factors, specifically precipitation of the coldest quarter and temperature seasonality, exert influence on the evolutionary history of Paeonia, with a notable contribution to niche conservatism. Similar studies have indicated that decreasing temperature and precipitation can facilitate higher evolutionary rates (Herbert et al., 2016). Furthermore, our study suggests that the diversification rates of Paeonia in Asia estimated were highest during the Quaternary period, specifically around 0-5 Mya, which is consistent with findings in Theaceae (Yu et al., 2017) and evergreen oak (Jiang et al., 2016). The rates of niche change within the clade were highest at 5-11 Mya, which precedes the period of highest diversification rates. However, this result contradicts the findings of a previous study (Folk et al., 2019).
5 Conclusions and prospects
In this study, we used positive selection analysis of molecular fragments, including ndhD-psaC, trnC-petN, trnK-rps16, ycf4-cemA, trnH-psbA, petG-trnW, petN-psbM, petA-psbJ, ycf1, and ycf2 to investigate the phylogenetic relationships within Saxifragales and Paeoniaceae. The phylogenetic results showed that Paeoniaceae formed an independent clade, and was a sister to Penthoraceae, Haloragidaceae, Iteaceae, Crassulaceae, and Saxifragaceae, with an estimated divergence time of 102-116 Mya. In addition, we inferred the ecological niche and spatial distributions of Paeonia in Asia and America, revealing that the potential habitat of species expanded to the southwestern edge of Asia and northern America. Moreover, the rates of niche change within the clade were observed to precede the diversification rates, occurring about 5 Mya. This study offers valuable insights into the evolution of Paeoniaceae and enhances our understanding of the relationships between evolutionary history and ecological niche change.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author contributions
YW: Writing – original draft. YC: Writing – original draft. ZP: Writing – review & editing. TM: Writing – original draft. YL: Writing – review & editing. CT: Writing – review & editing. XZ: Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by grants from the National Key R & D Program of China (2022YFD1200205), China Agricultural Research System (CARS-21), and Fundamental Research Funds for Chinese Academy of Agricultural Sciences (IVF-BRF2022016).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2024.1413707/full#supplementary-material
Supplementary Figure 1 | The estimations of Ka/Ks and Kn/Ks analysis within 11 species. (A) The estimations of Ka/Ks of all plastid protein-coding genes (PCG). (B) The estimations of Kn/Ks of all plastid noncoding regions.
Supplementary Figure 2 | The estimations of Pi analysis within 11 species. (A) The estimations of nonsynonymous (dN), synonymous (dS) substitution rates and dN/dS of all plastid protein-coding genes (PCG). (B) The estimations of nonsynonymous (dN), synonymous (dS) substitution rates and dN/dS of all plastid noncoding regions.
Supplementary Figure 3 | BI tree within Saxifragales. The numbers on the clade represent the Posterior probability.
Supplementary Figure 4 | The single-factor correlation analysis of 20 variables. The number of * symbols represents the magnitude of the P-value (* for P-value <= 0.05, ** for P-value <= 0.01, and *** for P-value <= 0.001).
Supplementary Figure 5 | The map of the collection and distribution of the Paeonia genus in Asian and North American.
Supplementary Figure 6 | The average test of AUC.
Supplementary Figure 7 | The contributions of the corresponding variable were plotted from 16 factors that affected the distributions of two periods.
Supplementary Figure 8 | Principal component analysis (PCA) in the ecological niche of 20 environmental factors.
References
Aguilée, R., Gascuel, F., Lambert, A., Ferriere, R. (2018). Clade diversification dynamics and the biotic and abiotic controls of speciation and extinction rates. Nat. Commun. 9, 3013. doi: 10.1038/s41467-018-05419-7
Ali, J. ,. R., Aitchison, J. C. (2014). Exploring the combined role of eustasy and oceanic island thermal subsidence in shaping biodiversity on the Galápagos. J. Biogeogr. 41, 1227–1241. doi: 10.1111/jbi.12313
Barnagaud, J. Y., Devictor, V., Jiguet, F., Barbet-Massin, M., Viol, I. L., Archaux, F. (2012). Relating habitat and climatic niches in birds. PloS One 7, e32819. doi: 10.1371/journal.pone.0032819
Bell, W. A. (1957). Flora of the upper cretaceous nanaimo group of vancouver island, british columbia. Geol. Surv. Can. 293, 1–84. doi: 10.4095/101457
Broennimann, O., Fitzpatrick, M. C., Pearman, P. B., Petitpierre, B., Pellissier, L., Yoccoz N, G., et al. (2011). Measuring ecological niche overlap from occurrence and spatial environmental. Glob. Ecol. Biogeogr. 21, 481–497. doi: 10.1111/j.1466-8238.2011.00698.x
Brown, J. L., Bennett, J. R., French, C. M. (2017). SDMtoolbox 2.0: the next generation Python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. PeerJ 5, e4095. doi: 10.7717/peerj.4095
Brown, R. W. (1962). Paleocene flora of the rocky mountains and great plains. U.S. Geol. Surv. Prof. Pap. 375, 1–119. doi: 10.3133/pp375
Chen, Q., Chen, L., Teixeira da Silva, J. A., Yu, X. (2023). The plastome reveals new insights into the evolutionary and domestication history of peonies in East Asia. BMC Plant Biol. 23, 243. doi: 10.1186/s12870-023-04246-3
Chen, Q. P., Yang, H., Feng, X., Chen, Q. J., Shi, S. H., Wu, C. I., et al. (2021). Two decades of suspect evidence for adaptive molecular evolution - Negative selection confounding positive selection signals. Natl. Sci. Rev. 9, nwab217. doi: 10.1093/nsr/nwab217
Clark, P. U., Dyke, A. S., Shakun, J. D., Carlson, A. E., Clark, J., Wohlfarth, B., et al. (2009). The last glacial maximum. Science 325, 710–714. doi: 10.1126/science.1172873
Cola, V. D., Broennimann, O., Petitpierre, B., Breiner, F. T., D’Amen, M., Randin, C., et al. (2017). ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography 40, 774–787. doi: 10.1111/ecog.02671
Dong, P. B., Wang, L. J., Jia, Y., Li, Z. H., Wang, H. Y., Guo, F. X., et al. (2022). Niche divergence at the intraspecific level in an endemic rare peony (Paeonia rockii): A phylogenetic, climatic and environmental survey. Front. Plant Sci. 13. doi: 10.3389/fpls.2022.978011
Dong, W. P., Xu, C., Li, C. H., Sun, J. H., Zuo, Y. J., Shi, S., et al. (2015). ycf1, the most promising plastid DNA barcode of land plants. Sci. Rep. 5, 8348. doi: 10.1038/srep08348
Dong, W. P., Xu, C., Wu, P., Cheng, T., Yu, J., Zhou, S. L., et al. (2018). Resolving the systematic positions of enigmatic taxa: manipulating the chloroplast genome data of Saxifragales. Mol. Phylogenet. Evol. 126, 321–330. doi: 10.1016/j.ympev.2018.04.033
Drummond, A. J., Rambaut, A. (2007). BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 8, 214–222. doi: 10.1186/1471-2148-7-214
Elith, J., Phillips, S. J., Hastie, T., Dudík, M., Chee, Y. E., Yates, C. J. (2011). A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 17, 43–57. doi: 10.1111/j.1472-4642.2010.00725.x
Fay, J. C., Wu, C. I. (2000). Hitchhiking under positive Darwinian selection. Genetics 155, 1405–1413. doi: 10.1093/genetics/155.3.1405
Fishbein, M., Hibsch-Jetter, C., Soltis, D. E., Hufford, L. (2001). Phylogeny of Saxifragales (angiosperms, eudicots): analysis of a rapid, ancient radiation. Syst. Biol. 50, 817–847. doi: 10.1080/106351501753462821
Folk, R. A., Stubbs, R. L., Mort, M. E., Cellinese, N., Allen, J. M., Soltis, P. S., et al. (2019). Rates of niche and phenotype evolution lag behind diversification in a temperate radiation. Proc. Natl. Acad. Sci. U.S.A. 116, 10874–10882. doi: 10.1073/pnas.1817999116
Gao, X. X., Liu, J., Huang, Z. H. (2022). The impact of climate change on the distribution of rare and endangered tree Firmiana kwangsiensis using the Maxent modeling. Ecol. Evol. 12, e9165. doi: 10.1002/ece3.9165
Guan, K. J. (1979). Flora reipublicae Popularis Sinicae Vol. 27 (Beijing, China: Science Press), 37–41.
Guisinger, M. M., Chumley, T. W., Kuehl, J. V., Boore, J. L., Jansen, R. K. (2010). Implications of the plastid genome sequence of Typha (Typhaceae, Poales) for understanding genome evolution in Poaceae. J. Mol. Evol. 70, 149–166. doi: 10.1007/s00239-009-9317-3
Guisinger, M. M., Kuehl, J. V., Boore, J. L., Jansen, R. K. (2011). Extreme reconfiguration of plastid genomes in the angiosperm family Geraniaceae: rearrangements, repeats, and codon usage. Mol. Biol. Evol. 28, 583–600. doi: 10.1093/molbev/msq229
Guo, Z. T., Sun, B., Zhang, Z. S., Peng, S. Z., Xiao, G. Q., Ge, J. Y., et al. (2008). A major reorganization of Asian climate by the early Miocene. Clim. Past. 4, 153–174. doi: 10.5194/cp-4-153-2008
Hadley, W., Bryan, J. (2019). readxl: Read Excel Files. Available online at: https://CRAN.R-project.org/package=readxl. (accessed February 16, 2019).
Haffer, J. (1997). Alternative models of vertebrate speciation in Amazonia: an overview. Biodivers. Conserv. 6, 451–476. doi: 10.1023/A:1018320925954
Hebbar, K. B., Abhin, P. S., Sanjo Jose, V., Neethu, P., Santhosh, A., Shil, S., et al. (2022). Predicting the potential suitable climate for coconut (Cocos nucifera l.) cultivation in India under climate change scenarios using the MaxEnt model. Plants 11, 731–754. doi: 10.3390/plants11060731
Herbert, T. D., Lawrence, K. T., Tzanova, A., Peterson, L. C., CaballeroGill, R., Kelly, C. S. (2016). Late Miocene global cooling and the rise of modern ecosystems. Nat. Geosci. 9, 843–847. doi: 10.1038/ngeo2813
Hong, D. Y. (2010). Peonies of the world: taxonomy and phytogeography (Richmond, UK: Royal Botanic Gardens).
Hong, D. Y. (2021). Peonies of the world: Phylogeny and evolution (London: Royal Botanic Gardens, Kew). In press.
Hua, X., Wiens, J. J. (2013). How does climate influence speciation? Am. Naturalist. 182, 1–12. doi: 10.1086/670690
Huang, J. L., Sun, G. L., Zhang, D. M. (2010). Molecular evolution and phylogeny of the angiosperm ycf2 gene. J. Syst. Evol. 48, 240–248. doi: 10.1111/j.1759-6831.2010.00080.x
Itterbeeck, J. V., Missiaen, P., Folie, A., Markevich, V. S., Damme, D. V., Guo, D. Y., et al. (2007). Woodland in a fluvio-lacustrine environment on the dry Mongolian Plateau during the late Paleocene: Evidence from the mammal bearing Subeng section (Inner Mongolia, P.R. China). Palaeogeogr. Palaeoclimatol. Palaeoecol. 243, 55–78. doi: 10.1016/j.palaeo.2006.07.005
Jiang, X. L., Deng, M., Li, Y. (2016). Evolutionary history of subtropical evergreen broad-leaved forest in Yunnan Plateau and adjacent areas: an insight from Quercus schottkyana (Fagaceae). Tree Genet. Genomes. 12, 104. doi: 10.1007/s11295-016-1063-2
Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A., Jermiin, L. S. (2017). ModelFinder: fast model selection for accurate phylogenetic estimates. Nat. Methods 14, 587–589. doi: 10.1038/nmeth.4285
Karney, C. F. F. (2013). Algorithms for geodesics. J. Geod. 87, 43–55. doi: 10.1007/s00190-012-0578-z
Katoh, K., Standley, D. M. (2013). MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772–780. doi: 10.1093/molbev/mst010
Klages, J. P., Salzmann, U., Bickert, T., Hillenbrand, C. D., Gohl, K., Kuhn, G., et al. (2020). Temperate rainforests near the South Pole during peak Cretaceous warmth. Nature 580, 81–86. doi: 10.1038/s41586-020-2148-5
Krawczyk, K., Nobis, M., Myszczyński, K., Klichowska, E., Sawicki, J. (2018). Plastid super-barcodes as a tool for species discrimination in feather grasses (Poaceae: Stipa). Sci. Rep. 8, 1924. doi: 10.1038/s41598-018-20399-w
Laland, K. N., Feldman, O. (1999). Evolutionary consequences of niche construction and their implications for ecology. Proc. Natl. Acad. Sci. U.S.A. 96, 10242–10247. doi: 10.1073/pnas.96.18.10242
Laland, K. N., Odling-Smee, F. J., Feldman, M. W. (1996). The evolutionary consequences of niche construction: A theoretical investigation using two-locus theory. J. Evol. Biol. 9, 293–316. doi: 10.1046/j.1420-9101.1996.9030293.x
Letunic, I., Bork, P. (2021). Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. NAR 49, W293–W296. doi: 10.1093/nar/gkab301
Li, H. T., Luo, Y., Gan, L., Ma, P. F., Gao, L. M., Yang, J. B., et al. (2021). Plastid phylogenomic insights into relationships of all flowering plant families. BMC Biol. 19, 232. doi: 10.1186/s12915-021-01166-2
Li, H. T., Yi, T. S., Gao, L. M., Ma, P. F., Zhang, T., Yang, J. B., et al. (2019). Origin of angiosperms and the puzzle of the Jurassic gap. Nat. Plants. 5, 461–470. doi: 10.1038/s41477-019-0421-0
Librado, P., Rozas, J. (2009). DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. doi: 10.1093/bioinformatics/btp187
Liu, H., Ye, Q., Wiens, J. ,. J. (2020). Climatic-niche evolution follows similar rules in plants and animals. Nat. Ecol. Evol. 4, 753–763. doi: 10.1038/s41559-020-1158-x
Losos, J. B. (2008). Phylogenetic niche conservatism, phylogenetic signal and the relationship between phylogenetic relatedness and ecological similarity among species. Ecol. Lett. 11, 995–1003. doi: 10.1111/j.1461-0248.2008.01229.x
Moreno, R., Zamora, R., Molina, J. R., Vasquez, A., Herrera, M. Á. (2011). Predictive modeling of microhabitats for endemic birds in south Chilean temperate forests using maximum entropy (Maxent). Ecol. Inform. 6, 364–370. doi: 10.1016/j.ecoinf.2011.07.003
Morgan, M., Cheng, Y. (2023). Organism.dplyr: dplyr-based access to bioconductor annotation resources. Bioconducter. Available online at: https://bioconductor.org/packages/Organism.dplyr. (accessed February 24, 2023).
Moyle, R. G., Oliveros, C. H., Andersen, M. J., Hosner, P. A., Benz, B. W., Manthey, J. D., et al. (2016). Tectonic collision and uplift of Wallacea triggered the global songbird radiation. Nat. Commun. 7, 12709. doi: 10.1038/ncomms12709
Murdoch, D. J., Chow, E. D. (1996). A graphical display of large correlation matrices. Am. Stat. 50, 178–180. doi: 10.1080/00031305.1996.10474371
Muscarella, R., Galante, P. J., Soley-Guardia, M., Boria, R. A., Kass, J. M., Uriarte, M., et al. (2014). ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods Ecol. Evol. 5, 1198–1205. doi: 10.1111/2041-210X.12261
Neubig, K. M., Whitten, W. M., Carlsward, B. S., Blanco, M. A., Endara, L., Williams, N. H., et al. (2009). Phylogenetic utility of ycf1 in orchids: a plastid gene more variable than matK. Plant Syst. Evol. 277, 75–84. doi: 10.1007/s00606-008-0105-0
Paradis, E., Schliep, K. (2019). ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528. doi: 10.1093/bioinformatics/bty633
Phillips, S. J., Anderson, R. P., Dudık, M., Schapire, R. E., Blair, M. E. (2017). [amp]]#x0301;Opening the black box: An open-source release of maxent. Ecography 40, 887–893. doi: 10.1111/ECOG.03049
Pound, M. J., Haywood, A. M., Salzmann, U., Riding, J. B. (2012). Global vegetation dynamics and latitudinal temperature gradients during the Mid to Late Miocene (15.97-5.33 Mya). Earth Sci. Rev. 112, 1–22. doi: 10.1016/j.earscirev.2012.02.005
Pyron, R. A., Burbrink, F. T. (2013). Phylogenetic estimates of speciation and extinction rates for testing ecological and evolutionary hypotheses. Trends Ecol. Evol. 28, 729–736. doi: 10.1016/j.tree.2013.09.007
Qiu, Y. X., Fu, C. X., Comes, H. P. (2011). Plant molecular phylogeography in China and adjacent regions: Tracing the genetic imprints of Quaternary climate and environmental change in the world’s most diverse temperate flora. Mol. Phylogenet. Evol. 59, 225–244. doi: 10.1016/j.ympev.2011.01.012
Raman, G., Park, S. J. (2016). The complete chloroplast genome sequence of Ampelopsis: gene organization, comparative analysis, and phylogenetic relationships to other angiosperms. Front. Plant Sci. 7. doi: 10.3389/fpls.2016.00341
Rebelo, H., Tarroso, P., Jones, G. (2010). Predicted impact of climate change on European bats in relation to their biogeographic patterns. Glob. Chang. Biol. 16, 561–576. doi: 10.1111/j.1365-2486.2009.02021.x
Ronquist, F., Huelsenbeck, J. P. (2003). MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. doi: 10.1093/bioinformatics/btg180
Sabeti, P. C., Reich, D. E., Higgins, J. M., Levine, H. Z., Richter, D. J., Schaffner, S. F., et al. (2002). Detecting recent positive selection in the human genome from haplotype structure. Nature 419, 832–837. doi: 10.1038/nature01140
Saleh, M. (2020). Clustering via Jenks Natural Breaks. Available online at: https://github.com/MSH19/Clustering-via-Jenks-Natural-Breaks. (accessed October, 2020).
Schwarz, E. N., Ruhlman, T. A., Weng, M. L., Khiyami, M. A., Sabir, J. S. M., Hajarah, N. H., et al. (2017). Plastome-wide nucleotide substitution rates reveal accelerated rates in Papilionoideae and correlations with genome features across legume subfamilies. J. Mol. Evol. 84, 187–203. doi: 10.1007/s00239-017-9792-x
Shapiro, J. A., Huang, W., Zhang, C., Hubisz, M. J., Lu, J., Turissini, D. A., et al. (2007). Adaptive genic evolution in the Drosophila genomes. Proc. Natl. Acad. Sci. U.S.A. 104, 2271–2276. doi: 10.1073/pnas.0610385104
Sloan, D. B., Alverson, A. J., Wu, M., Palmer, J. D., Taylor, D. R. (2012). Recent acceleration of plastid sequence and structural evolution coincides with extreme mitochondrial divergence in the angiosperm genus Silene. Genome Biol. Evol. 4, 294–306. doi: 10.1093/gbe/evs006
Song, W., Li, Y., Luo, A., Su, X., Wang, Q., Liu, Y., et al. (2024). Historical and contemporary climate jointly determine angiosperm plant diversity patterns across east Eurasia. Ecography 2024, e07062. doi: 10.1111/ecog.07062
Stamatakis, A., Hoover, P., Rougemont, J. A. (2008). Rapid bootstrap algorithm for the RAxML web servers. Syst. Biol. 57, 758–771. doi: 10.1080/10635150802429642
Suchard, M. A., Lemey, P., Baele, G., Ayres, D. L., Drummond, A. J., Rambaut, A. (2018). Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 4, vey016–vey021. doi: 10.1093/ve/vey016
Tada, R., Zheng, H. B., Clift, P. D. (2016). Evolution and variability of the Asian monsoon and its potential linkage with uplift of the Himalaya and Tibetan Plateau. Prog. Earth Planet. Sc. 3, 4. doi: 10.1186/s40645-016-0080-y
Tajima, F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595. doi: 10.1093/genetics/123.3.585
Talavera, G., Castresana, J. (2007). Improvement of Phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments. Syst. Biol. 56, 564–577. doi: 10.1080/10635150701472164
Tank, D. C., Eastman, J. M., Pennell, M. W., Soltis, P. S., Soltis, D. E., Hinchliff, C. E., et al. (2015). Nested radiations and the pulse of angiosperm diversification: increased diversification rates often follow whole genome duplications. New Phytol. 207, 454–467. doi: 10.1111/nph.13491
Warren, D. L., Glor, R. E., Turelli, M. (2008). Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62, 2868–2883. doi: 10.1111/j.1558-5646.2008.00482.x
Weng, M. L., Ruhlman, T. A., Jansen, R. K. (2016). Plastid-nuclear interaction and accelerated coevolution in plastid ribosomal genes in Geraniaceae. Genome Biol. Evol. 8, 1824–1838. doi: 10.1093/gbe/evw115
Wickham, H., Lionel, H. (2019). tidyr: Tidy Messy Data. Available online at: https://CRAN.R-project.org/package=tidyr.
Wiens, J. J., Graham, C. H. (2005). Niche conservatism: integrating evolution, ecology, and conservation biology. Annu. Rev. Ecol. Evol. Syst. 36, 519–539. doi: 10.1146/annurev.ecolsys.36.102803.095431
Wu, L. W., Nie, L. P., Wang, Q., Xu, Z. C., Wang, Y., He, C. N., et al. (2021b). Comparative and phylogenetic analyses of the chloroplast genomes of species of Paeoniaceae. Sci. Rep. 11, 1–16. doi: 10.1038/s41598-021-94137-0
Wu, L., Wu, M. L., Cui, N., Xiang, L., Li, Y., Li, X. W., et al. (2021a). Plant super-barcode: a case study on genome-based identification for closely related species of Fritillaria. Chin. Med. 16, 52. doi: 10.1186/s13020-021-00460-z
Xiao, S. M., Li, S. F., Wang, X. J., Chen, L. L., Su, T. (2022). Cedrus distribution change: past, present, and future. Ecol. Indic. 142, 159. doi: 10.1016/j.ecolind.2022.109159
Xu, W. B., Guo, W. Y., Serra-Diaz, J. M., Schrodt, F., Eiserhardt, W. L., Enquist, B. J., et al. (2023). Global beta-diversity of angiosperm trees is shaped by Quaternary climate change. Sci. Adv. 9, eadd8553. doi: 10.1126/sciadv.add8553
Xue, S., Shi, T., Luo, W., Ni, X., Iqbal, S., Ni, Z., et al. (2019). Comparative analysis of the complete chloroplast genome among Prunus mume, P. Armeniaca, and P. salicina. Hort. Res. 21, 89. doi: 10.1038/s41438-019-0171-1
Yang, Z., Nielsen, R. (2000). Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models. Mol. Biol. Evol. 17, 32–43. doi: 10.1093/oxfordjournals.molbev.a026236
Yang, Z. H., Wong, W. S. W., Nielsen, R. (2005). Bayes empirical Bayes inference of amino acid sites under positive selection. Mol. Biol. Evol. 22, 1107–1118. doi: 10.1093/molbev/msi097
Yu, X. Q., Gao, L. M., Soltis, D. E., Soltis, P. S., Yang, J. B., Fang, L., et al. (2017). Insights into the historical assembly of East Asian subtropical evergreen broadleaved forests revealed by the temporal history of the tea family. New Phytol. 215, 1235–1248. doi: 10.1111/nph.14683
Yule, G. U. (1927). On a method of investigating periodicities in disturbed series. Philos. Trans. R. Soc. London. 226, 267–298. doi: 10.1098/rsta.1927.0007
Zhai, W., Duan, X. S., Zhang, R., Guo, C. C., Li, L., Xu, G. X., et al. (2019). Chloroplast genomic data provide new and robust insights into the phylogeny and evolution of the Ranunculaceae. Mol. Phylogenet. Evol. 135, 12–21. doi: 10.1016/j.ympev.2019.02.024
Zhang, Z. (2022). KaKs_Calculator 3.0: calculating selective pressure on coding and non-coding sequences. Genomics Proteomics Bioinf. 20, 536–540. doi: 10.1016/j.gpb.2021.12.002
Zhang, W., Sun, Y. Z., Liu, J., Xu, C., Zou, X. H., Chen, X., et al. (2021). DNA barcoding of Oryza: conventional, specific, and super barcodes. Plant Mol. Biol. 105, 215–228. doi: 10.1007/s11103-020-01054-3
Zhang, Z. L., Zhang, Y., Song, M. F., Guan, Y. H., Ma, X. J. (2019). Species identification of dracaena using the complete chloroplast genome as a super-barcode. Front. Pharmacol. 10. doi: 10.3389/fphar.2019.01441
Zhou, Z. K., Crepet, W. L., Nixon, K. C. (2001). The earliest fossil evidence of the Hamamelidaceae: Late Cretaceous (Turonian) inflorescences and fruits of Altingioideae. Am. J. Bot. 88, 753–766. doi: 10.2307/2657028
Keywords: phylogenetic relationship, Paeoniaceae, Paeonia, niche change, niche
Citation: Wang Y, Chen Y, Prijic Z, Markovic T, Lyu Y, Tian C and Zhang X (2024) Whether niche changes promote the evolution of species: a case study of Paeonia in Asia and North America. Front. Plant Sci. 15:1413707. doi: 10.3389/fpls.2024.1413707
Received: 07 April 2024; Accepted: 28 October 2024;
Published: 15 November 2024.
Edited by:
Xiaohua Jin, Chinese Academy of Sciences (CAS), ChinaReviewed by:
Robinson J. Herrera Feijoo, Universidad Técnica Estatal de Quevedo, EcuadorClaudia Bita-Nicolae, Institute of Biology Bucharest of the Romanian Academy, Romania
Copyright © 2024 Wang, Chen, Prijic, Markovic, Lyu, Tian and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yingmin Lyu, bHV5aW5nbWluQGJqZnUuZWR1LmNu; Caihuan Tian, dGlhbmNhaWh1YW5AY2Fhcy5jbg==; Xiuxin Zhang, emhhbmd4aXV4aW5AY2Fhcy5jbg==
†These authors have contributed equally to this work