Evolutionary history of endangered and relict tree species Dipteronia sinensis in response to geological and climatic events in the Qinling Mountains and adjacent areas

Abstract Geological and climatic events are considered to profoundly affect the evolution and lineage divergence of plant species. However, the evolutionary histories of tree species that have responded to past geological and climate oscillations in central China's mountainous areas remain mostly unknown. In this study, we assessed the evolutionary history of the endangered and relict tree species Dipteronia sinensis in the Qinling Mountains (QM) and adjacent areas in East Asia based on variations in the complete chloroplast genomes (cpDNA) and reduced‐genomic scale single nucleotide polymorphisms (SNPs). Population structure and phylogenetic analysis based on the cpDNA variations suggested that D. sinensis could be divided into two intraspecific genetic lineages in the eastern and western sides of the QM (EQM and WQM, respectively) in East Asia. Molecular dating suggested that the intraspecific divergence of D. sinensis occurred approximately 39.2 million years ago during the later Paleogene. It was significantly correlated with the orogeny of the QM, where the formation of this significant geographic barrier in the region may have led to the divergence of independent lineages. Bayesian clustering and demographic analysis showed that intraspecific gene flow was restricted between the EQM and WQM lineages. Isolation‐with‐migration analysis indicated that the two genetic lineages experienced significant demographic expansions after the Pleistocene ice ages. However, the genetic admixture was determined in some populations between the two lineages by the large scale of SNP variations due to DNA incompatibility, the large significant population size, and rapid gene flow of nuclear DNA markers. Our results suggest that two different conservation and management units should be constructed for D. sinensis in the EQM and WQM areas. These findings provide novel insights into the unprecedented effects of tectonic changes and climatic oscillations on lineage divergence and plant population evolution in the QM and adjacent areas in East Asia.


| INTRODUC TI ON
Geological events and repeated climate oscillations have affected the genetic architecture of plant species and reshaped the divergence of their lineages (Abbott et al., 2000;Avise, 2000;Hewitt, 2000;Hickerson et al., 2009;Liu et al., 2012;Wen et al., 2014). Previously, it was shown that tectonic uplift could change the environmental conditions to stimulate species divergence and produce novel species lineages and/or adaptive zones. Factors such as mountain uplift and climatic cycles can alter the genetic compositions by limiting gene flow to affect evolutionary processes and result in lineage divergence in mountain populations (Avise, 2000;Hewitt, 2001;Hickerson et al., 2009). However, montane orogenesis led to fluctuations in the climate in high-altitude regions during the Quaternary period to severely impact the vegetation biota and genetic architecture by preventing population expansion and the migration of species with narrow thermal tolerance ranges (Avise, 2000;Knowles, 2001).
Previous highly cited studies broadly explored the effects of tectonic uplifts and climatic oscillations on the diversification patterns of plants species in Europe, North America, and other regions, for example, in the Sierra Madre of Mexico (Bryson et al., 2012), the Andes (Hughes & Eastwood, 2006), Costa Rica (Muñoz et al., 2019), and the Himalayas Yang et al., 2012). In China, many previous studies focused on this phenomenon in the complex and conspicuous topographical regions on the Qinghai-Tibet Plateau and adjacent areas (Ding et al., 2020;Wen et al., 2014;Xing & Ree, 2017). It was recently reported that climate factors and montane orogenesis had strong and direct effects on montane forests, and they shaped the montane plant diversity (Li et al., 2020;Muñoz et al., 2019). However, another critical biodiversity hotspot in central China is the Qinling Mountains (QM), and the detailed evolutionary histories of plant species in this area are generally unknown Wan et al., 2016;Zhang et al., 2015).
The QM comprises a unique range of mountains that run in an east-west direction in China's center, which provide large-scale complex and heterogeneous habitats (Jia et al., 2010;Qin et al., 2014).
The QM is also an important geographical and biodiversity hotspot in Eastern Asia, where it contains more than 337 endangered Chinese plant species. The original development of the QM landforms was formed by the collision between the Yangtze block and north China block during the Mesozoic and Cenozoic periods (Zhang et al., 1996), and the subsequent lateral extrusion of the Tibet Plateau in the early Cenozoic (approximately 40 Mya). In addition, the present high topographic gradients with a height drop of 1,000-2,000 m in the QM and the highest peak, Taibai Mountain, which is 3,767 m, resulted from the latest rapid uplift of the Kun-Huang tectonic changes and Gonghe mountains uplift activated by the Tibetan Plateau during the late Miocene and early Pleistocene Dong et al., 2011;Ying, 1994;Zhang et al., 2006). Moreover, some high mountains in the Qinling region (e.g., Taibai Mountain) were probably covered by ice caps during the Pleistocene ice ages (Li et al., 2004;Rost, 1994).
In combination with complex tectonic changes and harsh climatic oscillations, these mountains influenced the evolutionary history of various organisms (Bennett, 2004;Hua & Wiens, 2013;Zink et al., 2004). However, the populations present on high mountains that expanded during the postglacial periods might have contracted into low elevation refugia during the Pleistocene ice ages.
Moreover, several endemic species inhabit this region, which supports distinct habitats and climatic oscillations on the slopes of both temperate and subtropical ecosystems in the QM area that serve as the boundary of the Palearctic and Oriental realms between north and south China (Chen, 2008;Zhang, 2004). However, the QM could serve as a significant barrier to generate the north-south break for species generally distributed in the low altitude localities of the northern and southern sides, for example, <1,200 m above sea level (Qu et al., 2009;Yan et al., 2010). Some species in the high-altitude localities of the QM were also impacted by Pleistocene climatic fluctuations, which resulted in their diversification in the west-east and west-central-east breaks Wang et al., 2013;Yang et al., 2015). A recent study reported three lineages in the QM, that is, a northwestern QM lineage comprising the Tibet lineage, as well as eastern and western lineages, which separated about 3-4 million years ago (Mya) (Huang et al., 2017). However, our understanding of how geological events and climatic changes during the Pleistocene might be associated with spatial distributions, lineage sorting, and population evolution is restricted due to a lack of experimental observations Wang et al., 2012) in this important region of the QM with high biodiversity.
Dipteronia sinensis (Sapindaceae) is an endangered and relict woody plant species with ecological and economic value, mainly distributed in the QM and adjacent areas (Acevedo-Rodríguez et al., 2011). This species provides a good model for investigating the effects of geological and climatic events on intraspecific divergence and population evolution. Previous molecular phylogenetic studies demonstrated that Dipteronia and Acer are sister genera based on variations in their complete chloroplast genomes (cpDNA) (Renner et al., 2008;Yang et al., 2010) and single gene experiments based on cpDNA/nuclear DNA (Gao et al., 2017). However, other studies based on nuclear ribosomal DNA, cpDNA, and ITS/cpDNA showed that the phylogenetic relationship between Dipteronia and Acer is more distant, where they formed paraphyletic groups Zhou et al., 2016). Moreover, Dipteronia populations may have experienced a genetic bottleneck (Yang et al., 2008) that lasted over a million years. However, the evolutionary history and divergence during the glacial and postglacial periods are not well understood due to data limitations and the lack of management and K E Y W O R D S climatic oscillation, conservation, Dipteronia sinensis, lineage divergence, mountain uplift, population structure conservation for this tree species. Therefore, in this study, we analyzed the cpDNA and nuclear DNA for D. sinensis to determine the diversification, population genetic structure, and lineage divergence in this tree species under the impacts of mountain orogenesis and climatic oscillations in the QM and adjacent regions. In particular, we aimed: (a) to investigate how the divergence of the D. sinensis lineages might have been related to the orogenesis of the QM and past climatic oscillations; (b) to determine whether gene flow among populations was limited due to the geographical barrier and climate oscillations; and (c) to identify suitable conservation and management strategies to maintain the genetic variation in this endangered tree species.

| Plant samples and DNA extraction
We sampled 25 populations of Dipteronia sinensis (one individual from each population) from widespread localities in Sichuan, Hubei, Henan, Gansu, and Shaanxi provinces in its distribution range in the QM and adjacent regions of China. Detailed information for all of the sampled populations is provided in Table 1. All materials and documents were deposited in the College of Life Sciences, Northwest University, Xi'an, China. Genomic DNA was isolated using 5 mg of dried leaves from a sample of each D. sinensis population (preserved in silica gel) before DNA extraction with the CTAB method (Doyle, 1987) or a Tiangen plant DNA extraction kit (Beijing, China).
The DNA quality was checked with a 1% agarose gel based on the high molecular weight bands after gel electrophoresis.

| Assembly and annotation of cpDNA
High-quality DNA was sequenced using a 350-bp insert fragment by Novogene Company (Beijing, China). We then used the default parameters in NGSQC Toolkit v 2.3.3 to extract the clean data (Patel & Jain, 2012). After trimming, the clean reads were assembled using MITObim v 1.8 with GenBank Accession No. NC_029338 as the reference sequence (Zhou et al., 2017). Finally, the consensus sequences for each population data set were aligned and visually inspected in Geneious software V.7.1.5 and compared with a previously published Dipteronia sinensis genome (NC_029338) (Zhou et al., 2017). The newly sequenced complete chloroplast genomes were deposited in the NCBI GenBank under accession numbers MK193760-MK193784.

| Single nucleotide polymorphism (SNP) detection
In order to detect SNPs in nuclear DNA, 23 qualitative genomic DNA samples were selected from 23 populations (one sample from each population). We used quantitative real time-PCR to assess the quality of the DNA fragments before nuclear DNA library construction. Restriction endonuclease enzymes were applied to amplify each genome sample. All of the samples were mixed after PCR before selecting good quality fragments with the desired length based on the results of capture arrays, which were then used to construct the DNA library. The genotyping-by-sequencing technique was employed to sequence the selected nDNA library fragments. After sequencing, the raw reads were filtered based on the paired-end reads.
High-throughput sequencing obtained 23.15 GB of raw reads data from 23 samples, with an average of 0.99 GB per sample.
The clean reads comprised 23.14 GB and the average amount per sample was 0.99 GB. After high-quality sequencing (Q20 ≥ 92.8%, Q30 ≥ 82.85%), we also obtained the GC contents of the 23 Dipteronia sinensis samples. Burrows-Wheeler aligner software was used to align the clean data , and the DJ population was used as the reference genome. High-quality SNP loci were filtered with SAMtools and used for subsequent analyses ). GENALEX V.6.5 (Peakall & Smouse, 2012) was employed to detect the physical linkages between any pairs of loci.

| Phylogenetic analysis
Gene genealogies based on complete chloroplast genomes (cpDNA) were reconstructed by using MEGA v 7.0 and MrBayes v 3.2.3 (Ronquist & Huelsenbeck, 2003) to obtain the maximum likelihood (ML), maximum parsimony (MP), and Bayesian inference (BI) cladograms. To examine the phylogenetic relationships among the 25 populations of Dipteronia sinensis, we used the complete chloroplast genomes for Acer buergerianum (NC_034744.1), A.

and
Dimocarpus longan (NC_037447.1) as outgroup taxa. jModelTest and PAUP* v 4.0 software were used to obtain the general time-reversible model and gamma (G) distribution for the rate variation among sites based on Akaike's information criterion (Swofford, 2002). For the ML and MP phylogenetic cladogram analyses, we conducted 1,000 bootstrap replicates. To construct the best phylogenetic tree in MrBayes, we kept the burn-in as 2,500 and retained every 1,000 generations from 30,000,000 random tree rotations. The FigTree program was used to visualize the output (Rambaut, 2009).
We used 23 Dipteronia sinensis populations to perform phylogenetic analysis with high-quality nDNA SNP data. TreeBest software (http://trees oft.sourc eforge.net/treeb est.shtml) was used to construct a neighbor-joining phylogenetic tree based on 1,000 bootstrap values.

| Estimation of divergence time (cpDNA)
The BI method in BEAST v 1.8.4 was used to estimate the divergence times for the intraspecific Dipteronia sinensis lineages (Drummond & Rambaut, 2007). We used a fossil age of 56 Mya as the calibration point for divergence between Dipteronia and Acer (Feng et al., 2018).
We modeled two most recent coalescent ancestral (MRCA) points to obtain the best divergence time results (Guo et al., 2014) by using the log-normal distribution of the first calibration point with mean = 0, SD = 1, and offset = 60 Mya for the complete chloroplast genome data sets including outgroups. The crown age for the sister genera Acer and Dipteronia was used as one calibration point with a normal distribution, mean = 56 Mya, and SD = 1.
The Yule speciation tree prior and an "uncorrelated relaxed clock" model were employed with a normal prior distribution for the branch lengths. We set 50,000,000 Markov chain Monte Carlo (MCMC) algorithm steps, with 5,000,000 steps as the burn-in, and collected the parameters after every 5,000 steps. We used Tracer v 1.5 software (Rambaut & Drummond, 2009)

| STRUCTURE analysis
The model-based Bayesian algorithm in STRUCTURE v 2.3.4 software was used to analyze the genetic structure of Dipteronia sinensis. We implemented an admixture model with the additional correlated allele frequency between populations, as recommended in a previous study (Falush et al., 2007). Unlinked SNPs (

| Population dynamics history
To analyze the coalescent-based isolation-with-migration (IM) model for Dipteronia sinensis, we implemented six parameters in IMa software (Hey, 2009). All of the parameters were scaled by   theta (θ) (Nei & Li, 1979).
Principal component analysis (PCA) was conducted based on nDNA SNP data to visualize the significant axes of variation in the population genetics using the adegenet package (Jombart, 2008) in R software. We conducted PCA using the following formula.
PCA calculates the number of individuals "n" = XX autosomal data, while ignoring more than two alleles and inconsistent data. In the equation above, "i" is the number of individuals, "k" is the position of an SNP, d ik indicates that if an individual "i" and the reference allele are homozygous, then d ik = 0, whereas if they are heterozygous, then d ik = 1, and if the "i" individual is homozygous for the non-reference allele, then d ik = 2. In this study, the eigenvectors and eigenvalues were calculated using GCTA (http://cnsge nomics. com/softw are/gcta/pca.html) software and the PCA plot was drawn using R software.

| Lineage relationships and divergence
The cpDNA data set comprised 25 population samples, where each genome had similar a typical quadripartite structure to those of most land plants. The newly sequenced chloroplast genomes were identical to the previously published plastomes for Dipteronia sinensis in terms of their structure, tRNAs, mRNAs, and gene contents (Table 2) (Zhou et al., 2017). The chloroplast genome data set for each population had different assembly reads and aligned base-pair length ranges. The numbers of tRNAs, mRNA, and protein-coding genes were the same as those in the previously published genome (Table 2).
Phylogenetic analysis based on the complete chloroplast genome data set obtained a similar topology using the three methods employed (ML, MP, and BI cladograms) (Figure 2). The phylogenetic tree formed a larger bi-phyletic clade with high bootstrap support values, but two clades clearly diverged between the eastern and western sides of the QM (EQM and WQM, respectively) (Figures 1 and 2).
An identical topology was obtained with the highest posterior probability of 1.0 by BI analysis. We also estimated the divergence times for the intraspecific Dipteronia sinensis lineages based on the cpDNA data. Our molecular dating results suggested that the species D. sinensis diverged from its relatives during the early Paleogene at 62.7

| STRUCTURE analyses
The results obtained by STRUCTURE analyses based on the cpDNA variations were the same as those produced using the phylogenetic tree. They suggested the clear intraspecific divergence of Dipteronia sinensis into EQM and WQM lineages (Figures 2 and 4a). However, no gene flow was detected between the EQM and WQM lineage populations (Figure 4a). respectively. Bootstrap support values are shown above the branches (>50% in each case). The first bootstrap value represents that for maximum likelihood (ML), the middle value represents that for maximum parsimony (MP), and the last value represents that for Bayesian inference (BI) also detected three major populations groups, where two groups were clearly the EQM and WQM populations, and the other was a mixture of both sides.

| Demographic dynamics
We calculated the nucleotide diversity Pi (π) and Theta (θ) by using DnaSP software based on the SNP loci (cpDNA). We estimated that the nucleotide diversity was much higher for the EQM lineage (π = 0.00026 and θ = 0.00031) than the WQM lineage (π = 0.00009 and θ = 0.00011).
IMa software was used to analyze the gene flow migration be- PCA (PC1, PC2, and PC3) identified the same genetic pattern as STRUCTURE analysis based on the nDNA data. Three different genetic groups were determined by PCA, which contained the EQM and WQM populations, and a mixture of the populations from both sides, as shown in Figure 5. The three axes comprising PC1, PC2, and PC3 estimated very low genetic variation among each group, that is, 7.14%, 6.07%, and 5.73%, respectively ( Figure 5). The PC1 geographic axis indicated the high correspondence between the genetic data for the EQM and WQM population groups. Figure 5 shows that the EQM populations (TY, TBS, HPS, and HH) and WQM populations (LD, MJS, ZA, NBL, HHG, LY, and DC) were on the two extreme sides of the distribution map. By contrast, the remaining populations were isolated from the groups and formed a "scattered populations group" on the PC1 axis. The geographic distance on the PC3 axis was a substructure of PC1 representing the increase in the geographic range from the EQM to WQM population groups ( Figure 5). The large gap between the two groups of populations suggested a significant difference in those genetic heredity groups. Ideally, individuals with a similar genetic background will be clustered together by PCA.

F I G U R E 3
Chronogram obtained forDipteronia sinensisbased on the whole chloroplast genome sequences using BEAST. Red and green indicate the populations from the east and west sides of the Qinling Mountains denoted as "EQM and WQM," respectively. The mean divergence time and 95% HPDs are labeled above the line in Mya

| Intraspecific divergence of D. sinensis
We examined the divergence of Dipteronia sinensis lineages based on the complete chloroplast genomes (cpDNA) and large-scale SNP variations in the QM and adjacent areas. Phylogenetic analysis based on cpDNA clearly determined the bi-phyletic clades within the species (EQM versus. WQM) (Figures 1 and 2) which is consistent with a previous phylogenetic analysis based on transcriptome data sets (Feng et al., 2018). Therefore, our phylogenomic dating strongly suggests that D. sinensis and D. dyeriana qualify as "living fossils" that have experienced long-term morphological stasis, possibly since the Paleocene/Eocene (Figure 3). They survived the late Tertiary/Quaternary periods with adverse climatic conditions in subtropical China's mountainous areas (Qiu et al., 2011).
The orogenesis of the QM was estimated to have begun in the Late Triassic (~230 Mya) (Zhang et al., 1996). However, the current topographical features of the QM mainly formed starting with a small-magnitude uplift at 40-70 Mya and a rapid large-magnitude uplift beginning at 1.2-2.4 Mya (or 0.7 Mya), especially for the higher mountain features of the Taibai Mountains (Dong et al., 2011;Liu et al., 2010;Ying, 1994). Given that the main phylogeographical patterns occur along the west-east axis, we suggest that histori- indicated a smaller effect on the population size due to low mutation and evolutionary rates compared with the nDNA data. By contrast, the nDNA data revealed a more rapid evolutionary rate and larger significant population size than the cpDNA data (Wolfe et al., 1987). Therefore, the rate of lineage sorting was more rapid for cpDNA than nDNA (Toumi & Lumaret, 2010).

F I G U R E 5
Principal component analysis results obtained forDipteronia sinensisbased on SNPs (nDNA) data. Red, green, and blue dots represent the east, west, and mixed populations, respectively, and red, green, and blue circles represent the east, west, and mixed population groups, respectively. Each axis shows the percentage of genetic variation in brackets Moreover, the low evolutionary rate of the cpDNA markers restricted the gene flow (Philippe et al., 2011), and the fixed lineage sorting of the ancestral polymorphisms may have contributed to the incongruence between the results obtained based on cpDNA and nDNA data . We found that when K = 2, the STRUCTURE results (K = 2) based on nDNA SNP data comprised two  (Meng, 2017;Vitelli et al., 2017). Another factor related to the incongruent results is nDNA incompatibility, as described by Abe's rules of nDNA and cpDNA dispersal (as in Haldane's rule), or sex differences (Abe et al., 2005;Di-Candia & Routman, 2007). For example, cpDNA is maternally inherited (e.g., seeds) (Bartish et al., 2002) and nDNA is biparentally inherited (e.g., seeds and pollens) (Sun et al., 2002).

| Demographic history
In the present study, STRUCTURE analysis obtained a clear genetic pattern for Dipteronia sinensis, and the results indicated no gene flow between two lineages in the EQM and WQM areas. The results obtained by IM analysis (with IMa) also suggested little gene flow between the two lineages. Previously, it was reported that divergence due to geographical tectonics or climatic isolation could produce a strong signature because of highly restricted gene flow (Avise, 2012;Parchman et al., 2006). Several studies have suggested that climatic change during the ice ages hindered the gene flow among populations, thereby resulting in the divergence of species lineages and the formation of new population lineages/taxa (Jiang et al., 2007;Ye et al., 2016). Studies have also shown that the gene flow might have been limited due to the QM forming a robust geographical barrier between northern and southern China Yuan et al., 2012;Zhang et al., 2014).
Similarly, evidence indicates that the formation of the Yangtze and Yellow Rivers due to the QM uplift during the Pliocene affected the genetic structure and limited the gene flow (Duan et al., 2009;Zhao et al., 2011). A study based on Paeonia rockii showed that the gene flow was limited between populations due to geographic isolation in the QM area (Yuan et al., 2010(Yuan et al., , 2011. Thus, these results suggest that tectonic changes and climatic adaptation during the Pliocene might have counteracted interspecific gene flow, thereby promoting the divergence of lineages (Feder et al., 2012).
Interestingly, our IMa estimates indicated that the daughter populations (q1) of Dipteronia sinensis were much larger than those of the ancestral population (qA) ( with climatic changes (Avise, 2000;Hewitt, 2001). Therefore, our IM results support the hypothesis that the Pleistocene climatic oscillations were related to genetic variation and the expansion of D.
sinensis in the QM and adjacent areas (Table 3). Previously, it was reported that climatic changes could drive range expansions of species (Avise, 2000) to promote phenological radiation and the generation of various morphotypes or species diversification (Milá et al., 2007).
Phylogeographical studies suggest that populations of alpine species exhibited evident range retreat patterns and expansion during the Pleistocene climate oscillations (Zhang et al., 2005). Based on the discussion above, we may conclude that the divergence and expansion of species populations triggered by Pleistocene changes promoted species differentiation in plants.

| Conservation strategies for D. sinensis
The ecological and evolutionary units that maintain biodiversity should be protected and conserved. Dipteronia sinensis is a naturally endangered and relict tree species included in the Red List, and thus, there is a critical need to design effective conservation strategies for this plant (Cowling & Pressey, 2001;Davis et al., 2008). This is the first study to provide details of the evolutionary history and genetic structure of D. sinensis based on cpDNA and nDNA SNP data.
In addition, our analysis of the demographic history of D. sinensis indicated low nucleotide diversity and restricted gene flow between the WQM and EQM clades of the widely distributed populations of D. sinensis, thereby supporting our conclusion that the WQM and EQM populations are not interbreeding (Figure 4a). However, during field investigations, we found that illegal human activities (e.g., continuous harvesting, development of agriculture, the rise of tourism activities, low reproductive capacity, hydroelectric development, and other industrial actions) and genetic fragmentation have significantly depleted several natural populations of D. sinensis (Hamilton et al., 2012;Zhou et al., 2010). Due to the continuous decreases in the sizes of wild populations, it is necessary to protect the populations affected by human activities in order to save the natural ecosystem. The direct effects of climate change are even more significant on alpine ecosystems, and they will be a much greater challenge for nature conservation in the 21st century. Details of factors such as topographic features, climatic changes, the growth of alpine ecology, and the population genetic structure and phylogeography of D. sinensis are necessary to facilitate the conservation of valuable plant habitat resources (Schmitt, 2007). Therefore, two different conservation and management units should be constructed for D.
sinensis in the EQM and WQM areas. In addition, the conservation of natural populations, habitat diversity, and enhanced gene exchange among populations in different locations can be improved by collecting mature seeds from various populations and artificially planting them in botanical gardens (Cabrera-Toledo et al., 2012).

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare. Zhong-Hu Li https://orcid.org/0000-0001-6763-5586