Molecular Phylogeography and Population Genetic Structure of O. longilobus and O. taihangensis (Opisthopappus) on the Taihang Mountains

Historic events such as the uplift of mountains and climatic oscillations in the Quaternary periods greatly affected the evolution and modern distribution of the flora. We sequenced the trnL–trnF, ndhJ-trnL and ITS from populations throughout the known distributions of O. longilobus and O. taihangensis to understand the evolutionary history and the divergence related to the past shifts of habitats in the Taihang Mountains regions. The results showed high genetic diversity and pronounced genetic differentiation among the populations of the two species with a significant phylogeographical pattern (N ST>G ST, P<0.05), which imply restricted gene flow among the populations and significant geographical or environmental isolation. Ten chloroplast DNA (cpDNA) and eighteen nucleus ribosome DNA (nrDNA) haplotypes were identified and clustered into two lineages. Two corresponding refuge areas were revealed across the entire distribution ranges of O. longilobus and at least three refuge areas for O. taihangensis. O. longilobus underwent an evolutionary historical process of long-distance dispersal and colonization, whereas O. taihangensis underwent a population expansion before the main uplift of Taihang Mountains. The differentiation time between O. longilobus and O. taihangensis is estimated to have occurred at the early Pleistocene. Physiographic complexity and paleovegetation transition of Taihang Mountains mainly shaped the specific formation and effected the present distribution of these two species. The results therefore support the inference that Quaternary refugial isolation promoted allopatric speciation in Taihang Mountains. This may help to explain the existence of high diversity and endemism of plant species in central/northern China.


Introduction
During the last three million years, repeated glacial and interglacial cycles have greatly affected the landform and fauna of the total earth [1][2][3]. Historical ecological and biogeographic factors are often considered to have played significant roles in shaping global biodiversity by influencing regional differences in speciation, extinction, and migration [4][5][6]. The genetic structures of the current species record the simultaneous consequences of two fundamental processes: population dynamics in response to past geological or climatic changes, and lineage sorting within a species under natural selection. Thus, knowledge of the evolutionary history of many plant species is central to the identification of divergence and speciation processes [3,[7][8].
Phylogeographical analyses can provide an understanding of how paleo-environmental changes in landscape and climate have influenced species distributions and population demography [9]. Such analyses play an important role in understanding the evolutionary history of species with changes [3]. Over the last decade, molecular research on the evolutionary history of plant species, both in China and throughout East Asia, has been conducted with reference to past climatic oscillations [10][11][12][13][14][15][16][17][18]. One of the results of such research is the proposal by some authors that glacial refugia were maintained in both the northern and the southern regions, or at different spatial-temporal scales in China, during these glacial periods. These refugia are thus suggested to have acted as sites for subsequent range expansion during the interglacial (or postglacial) periods [12][13][14][19][20]. Some plant species have extended their distribution from southwestern China to central/northern China [21]. Reduction of genetic diversity is expected under a scenario of rapid postglacial expansion, as has been found in northern Europe and America [2]. Because central/ northern China is believed not to have been glaciated during the Quaternary [2,[22][23], it is an open question whether species experienced past northern expansion and southern retreat during the Quaternary climatic oscillations.
Taihang Mountains locates on 35u199-40u519N and 113u109-115u489E. As a natural boundary mountain of the east, southeast of Shanxi province and Hebei and Henan provinces, it is the important mountain range and geographical boundary for eastern China. The northern part of Taihang Mountains is higher than its southern part. Most areas of Taihang Mountains are higher than 1200 meters. The uplift of Taihang Mountains, followed by the formation of high mountains and deep valleys within the plateau [24][25], was one of the most important geological events. Several lines of evidence suggest that the rapid uplift of Taihang Mountains took place after the late Pleistocene [26][27]. Rich sources of species, most of which are endemic, are found in this region [24]. The high species richness in Taihang Mountains has led to the hypothesis that this region is a distribution and diversity center for many plant genera [24].
Molecular techniques have provided powerful tools for studying the phylogeography or migratory footprints of species [7][8][9][28][29]. Maternally inherited chloroplast DNA lineages in natural populations often display geographical structure [9], which is useful for deciphering the evolutionary history of species. The cpDNA markers are thought to be the most appropriate candidates because of their slow evolution and lack of recombination [30]. Nevertheless, the joint use of molecular markers derived from different genomes provides a more complete description of population structure and insights into population history and dynamics, particularly for comparisons of maternally inherited organelle and bi-parentally inherited nuclear markers [31]. Several studies about inter-specific and intra-specific divergence from China, using different molecular markers, have played an important role in discovering phylogeographical patterns in East Asian flora. This is particularly true for the mechanisms of plant speciation, along with the production of the high plant biodiversity and endemism found in East Asian flora [12][13]15,20,32]. However, most of the studied species from China have been trees [13,[33][34]. Herbaceous vegetation has received much less attention [16]. Therefore, it would be of great interest to study the phylogeography of an herbaceous species to understand the evolution and modern distribution of the vegetation, especially in Taihang Mountains of central/northern China.
Opisthopappus longilobus Shih and Opisthopappus taihangensis (Ling) Shih belong to Opisthopappus (Asteraceae). The genus Opisthopappus is endemic to China, and its wild distribution is mainly restricted to Taihang Mountains across the provinces of Shanxi, Hebei, and Henan, occurring on the slopes at an elevation of about 1000 m or in the cracks of the steep cliffs [35]. O. longilobus is mainly distributed in the province of Hebei, whereas O. taihangensis is found in the provinces of Shanxi and Henan [35][36]. Both species are diploid (2n = 18) [25]. A strict morphological differentiation occurred between O. longilobus and O. taihangensis. The leaf blade of O. longilobus is smooth and subcylindrical, and most stem and leaves are pinnatifid, with one pair of bracteal leaf. O. taihangensis, apprised pubescent on both surfaces of leaf blade, stem, and leaves, is bipinnatifid, with no bracteal leaf [35][36]

Ethics statement
This study was conducted in accordance with all People's Republic of China laws. No specific permits were required for the described field studies because all researchers collecting the samples had introduction letters from the College of Life Science, Shanxi Normal University, Linfen.  Figure 1). Each population included 10 to 25 individuals that were collected at least 5 m apart. Healthy leaves were collected and dried in silica gel until total DNA was extracted. In each of the 13 populations, the ecological factors of each location were determined and recorded (Table 1).

DNA extraction, amplification and sequencing
Total DNA was extracted from the silica gel-dried leaves using the modified 26CTAB procedure [37]. DNA quality was checked by electrophoresis on 0.8% agarose gels, and DNA concentration was determined using an Eppendorf biophotometer protein nucleotide analyzer (Eppendorf China Ltd., Beijing, Germany). The DNA samples were diluted to 10 ng?mL 21 and stored at 20uC for subsequent use.
The trnL-trnF sequences were amplified using trnL and trnF [38].The primers ndhJ and trnL [39] were used to amplify the ndhJ-trnL sequences. ITS sequences were then amplified using the primers ITS4 and ITS5 [40]. A polymerase chain reaction (PCR) was then performed in a 25 mL volume, with 50 ng plant DNA, 26MasterMix (0.2 mM dNTPs, 3 mM MgCl 2 , 16PCR buffer, and 0.1 unit Taq DNA polymerase), and 0.6 mM of both forward and reverse primers. The PCR parameters for all amplification programs of ITS and cpDNA were as follows: 4 min of predenaturation at 94uC, followed by 34 cycles of 30 s of denaturation, 40 s of annealing (49.2uC for ITS or 58.8uC for trnL-trnF or 59.8uC ndhJ-trnL), 1 min 20 s of elongation at 72uC, and a final elongation step of 7 min at 72uC. PCR products were purified using the Wizard PCR Preps DNA purification system (Promega, Madison, WI, USA) following the manufacturer's instructions. Cycle sequencing reactions were conducted using the purified PCR product, AmpliTaq DNA polymerase, and fluorescent BigDye terminators. The sequencing products were analyzed using an ABI Prism 310 DNA sequencer (Applied Biosystems Inc., Foster City, CA, USA).

Data analysis
The sequences were edited manually based on the chromatograms and aligned by CLUSTAL X [41] and then adjusted manually. Inserts and indels within all cpDNA and nrDNA sequences were firstly treated as a single character resulting from  Table 1, and locality numbers correspond to those in Table 1  one mutation. Haplotype diversity (h) and nucleotide diversity (p) were calculated for each population (h, p) and at the species level (h d , p d ) using DNAsp [42]. The program PERMUT [44] was used to calculate the within-population diversity (h S ), total diversity (h T ), geographical total haplotype diversity (V T ), geographical average haplotype diversity (V S ), level of population differentiation at the species level (G ST ), and an estimate of population subdivisions for phylogenetically ordered alleles (N ST ). We further tested the phylogeographical structure at the species range between (N ST ) and (G ST ) by using the U-statistic. An N ST value higher than the G ST value indicates that closely related haplotypes are observed more often in a given geographical area than would be expected by chance [43].
To quantify the genetic differentiation partitioned among different groups and total genetic variance, analyses of molecular variance (AMOVA) were carried out using ARLEQUIN [44], with 1000 random permutations to test for significance of partitions. The spatial genetic structure of haplotypes was analyzed by spatial analysis of molecular variance using SAMOVA [45]. This program uses a simulated annealing approach to define groups of populations (K) that are geographically homogenous and maximally differentiated from each other. In this analysis of haplotypes, K varied from 2 to 13, with each simulation starting from 100 random initial conditions. An F CT index of genetic differentiation among initial K groups was computed, followed by an iterative simulated annealing process to obtain the optimal configuration of groups and final F CT values. The simulated annealing process was repeated 1000 times. The configuration with the largest F CT value among the 100 tested was retained as the best grouping of populations.
Genealogical relationships between haplotypes were inferred from a maximum parsimony median-joining network calculated in NETWORK 4.5.0.2 [46]. To complement the results of NETWORK, TCS 1.21 [47] was also used to construct haplotype relationships. Phylogenetic analyses were conducted for the nrDNA and cpDNA sequence data using maximum likelihood (ML) and Bayesian inference (BI), respectively.
The time of the most recent common ancestor (TMRCA) of all haplotypes was estimated via a Bayesian approach implemented in *BEAST (Star BEAST) [48][49] using a GTR+I substitution model. The best substitution model was determined according to the Akaike Information Criterion (AIC) in jModeltest [50]. These models were applied in ML, BI, and subsequent *BEAST analyses. Bayesian inference was performed using MrBayes 3.1.2 [51]. Two independent runs of Metropolis-coupled Markov chain (MCMC) analysis were executed, each including one cold chain and three incrementally heated chains that started randomly in the parameter space. Two independent runs of 10 8 generations were carried out, with sampling at every 1000 generations. The first 25% of sampled trees were discarded as burn-in, and the remaining trees were used to construct a Bayesian consensus tree. The convergence of chains was checked using Tracer 1.5 [52]. The remaining trees were pooled to estimate the posterior probabilities (PPs).
No fossil records or biogeographic events isolating distinct populations are available to calibrate a cpDNA-IGS substitution rate for O. longilobus and O. taihangensis. Therefore, for our combined chloroplast non-coding regions, we assumed minimum and maximum values of a range of average mutation rates reported for synonymous sites of plant chloroplast genes [i.e., 1.2 and 1.7610 29 substitutions per site per year (s/s/y)] [53]. These rates were then used for estimating TMCRA in *BEAST under a relaxed molecular clock assumption.
To detect whether population groups experienced recent population expansion and satisfied the assumption of neutrality, a mismatch distribution analysis was performed using the DNAsp program. Tajima's D and Fu and Li's D* were calculated for the entire genus and groups of populations. The statistical significance of D and D* was estimated with coalescent simulations as implemented in this program [54][55]. To further infer demographic processes, the null hypotheses of a spatial expansion and of a pure demographic expansion were tested in ARLEQUIN by comparing observed and expected distributions of pairwise sequence differences (mismatch distributions). For each expansion model, goodness-of-fit was tested using the sum of squared differences (SSD) and Harpending's [56] raggedness index (HRag).

Patterns of variability in cpDNA and ITS sequences
The trnL-trnF and ndhJ-trnL intergenic spacers varied from 800 to 866 bp and from 790 to 835 bp, respectively. The trnL-trnF intergenic spacer, which exhibits considerable nucleotide polymorphism, is characterized by 10 substitutions and a 4 bp insertions and deletions (Table S1). The 4 bp insertions and deletions region was mainly observed in O. longilobus populations. The ndhJ-trnL spacer was characterized by one indel, three substitutions, and a 10 bp insertions and deletions (Table S1). The total cpDNA combined matrix comprised 1701 sites, of which 13 positions were variable and 18 harbored gaps. Ten haplotypes ( (Table 1).
The full ITS sequences, including ITS1+5.8S+ITS2, were 682 bp in length, comprising of 254 bp and 224 bp for ITS1 and ITS2, respectively. Eighteen haplotypes ( Figure 2B) and 17 polymorphic sites, the latter consisting of 15 parsimony informative sites and 2 singleton variable sites, were detected at ITS sequences. Two populations, SBY and XTS, were monomorphic, whereas the remaining populations were polymorphic (  Based on ITS sequences, phylogenetic analyses were carried ( Figure S1). All studied populations made up a single monophyletic lineage, indicating that all of these divergent populations were derived within Opisthopappus itself, rather than the result of introgression by other species. Spatial genetic analyses of cpDNA and nrDNA haplotypes in all 13 populations using SAMOVA indicated that F CT increased to a maximal value of 0.921 when K = 2 (K, the number of groups). Thus, the division of all the 13 sampled populations approximately into two groups is appropriate.

Genetic diversity and population structure
AMOVA analysis revealed that 60.58% and 67.70% (P,0.01) of the total variation was due to differences among the cpDNA and nrDNA populations for all studied populations, respectively ( Table 2). When the populations were grouped by taxonomic variety, the cpDNA showed 27.28% variation among populations within species (P,0.01; Table 2). Similarly, for the nrDNA, AMOVA revealed that 26.68% of the variation was attributed to differences among populations within the species (P,0.01; Table 2). There was a significant (P,0.01) variation between species, suggesting that the two separately distributed groups have a genetic differentiation. In addition, 60.71% of the total cpDNA variation and 68.24% of the total nrDNA variation existed among O. longilobus populations ( Table 2). The genetic structure of O. taihangensis populations showed a similar trend with that of O. longilobus populations, with most of the genetic differentiation occurring among populations ( Table 2).
Results of the Mantel test showed a significant correlation between nrDNA data differentiation among populations and the natural logarithm of the geographical distances throughout the sampled range of O. longilobus (r = 0.442; P = 0.001) and O. taihangensis (r = 0.385; P = 0.001). Likewise, a significant correlation was detected among populations from cpDNA data of O. longilobus (r = 0.401; P = 0.001) and O. taihangensis (r = 0.399; P = 0.001). This indicates a significant correlation between genetic differentiation and geographical distance, implying strong IBD.

Haplotype and phylogeographical structure
In this study, similar topologies were obtained from the two different network approaches used to infer the relationships among the cpDNA and nrDNA haplotypes. Only the median-joining network is shown. The haplotype H8 was the most frequent haplotype in cpDNA data, occurring in 6 of the 13 populations. Three haplotypes, H1, H3, and H8, were shared by O. longilobus and O. taihangensis populations (Table 1 and Figure 3). In nrDNA data, Hn1 was the most frequent haplotype, occurring in 10 of the 13 populations (Figure 3 and Figure 4). The haplotypes Hn1, Hn4, and Hn5 were shared by O. longilobus and O. taihangensis populations (Table 1 and Figure 4). The distributions of haplotypes that were restricted to a single population were extremely skewed. Four cpDNA and nine nrDNA haplotypes were found in O. longilobus single populations. Only one cpDNA haplotype was identified for O. taihangensis populations.
Based on cpDNA data, we estimated a G ST value (0.649, P, 0.05), which is significantly smaller than the N ST value (0.761). When the nrDNA was examined, the G ST value (0.253) also   The phylogeographical relationships of the 10 cpDNA haplotypes and 18 nrDNA haplotypes were assessed under Neighborjoining analyses (NJ) and Bayesian inferences drawn using Hippolytia and Nipponanthemum as outgroups (Figure 3 and Figure 4). Similar topology to phylogeographical relationships indicated a single phylogenetic split. All haplotypes were clustered into two lineages: I and II (Figure 3 and Figure 4). Lineage I contained some populations of O. longilobus with a high bootstrap support. The remaining haplotypes formed lineage II, containing some O. longilobus populations and all O. taihangensis populations, also with a high bootstrap support. This unrooted network of cpDNA haplotypes is consistent with the strict consensus tree produced by NJ and Bayesian inference and furthermore displays the relationship of interior (ancestral) and tip (derived) haplotypes ( Figure 3). For lineage I (Figure 3), the haplotype H3 located in the center and fixed in three populations, being an ancestral haplotype. In clade II, H1 was shared by five populations. It appears to be an ancestral haplotype distributed in the center of clade II. The haplotype H8 was located in the tip of the network despite though contained a higher frequency of haplotypes. In the nrDNA haplotype network (Figure 4), Hn8 was ancestral for the lineage I; Hn1 was the predominant and widespread one in the clade II.
To test whether the current expansion had occurred in different populations, we calculated the frequency distribution of pairwise nucleotide differences among O. longilobus populations and O. taihangensis populations respectively. The resultant mismatch distribution consisted of multiple multimodal curves ( Figure S2). Neutrality tests revealed nonsignificant positive values in considering all populations ( Table 3). The observed mismatch distributions of cpDNA haplotypes were multimodal for all groups except for O. taihangensis. The mismatch distributions based on nrDNA haplotypes were present a similar pattern with cpDNA data (Table 3)

Taxonomic implications
At the morphological character levels, the evidences that illustrated O. taihangensis as a separate species from O. longilobus were only leaf form and involucres. So, these two species were once regarded as a species. In this study, the analyses of the combined cpDNA and nrDNA matrix support the monophyly of two clades ( Figure 3A, Figure 4A, and Figure S1). The phylogenetic relationships among all populations show an overall congruence with the result of the haplotype network ( Figure 3B and Figure 4B). All haplotypes are split into two major groups according to cpDNA and nrDNA data. The populations from two distinct clades approximately exhibited specific geographical distributions in Taihang Table 2). By used SRAP and cpSSR markers [58][59], AMOVA indicated that a significant genetic differentiation between O. longilobus and O. taihangensis. Accordingly, our molecular results support that O. longilobus and O. taihangensis were regarded as two species rather than a species. That result is consistent with the point of Hu [25]. The divergence times between O. longilobus and O. taihangensis ( Figure 3A) indicate that genetic divergence between two species occurred in early Pleistocene. This split was related to the ecological and geographical habitat changes resulting from climate oscillations during the Quaternary glacial periods and complicated topography with the uplift of Taihang Mountains. Geographic patterns of genetic diversity and structure The high levels of genetic diversity were observed in O. longilobus and O. taihangensis populations (Table 1). These results are similar to previous research on herbaceous plants in East Asia and significantly exceed the average value of 0.67 for the 170 plant species documented [8,12,19,31,60]. In O. longilobus, the WDS population obtained the most haplotypes and highest genetic diversity, followed by the SHS population. This suggests the presence of two genetic diversity centers for O. longilobus, i.e., the regions where WDS and SHS populations are located. For O. taihangensis, the BJY and WML populations, which obtained the high genetic diversity, were genetic diversity centers for O. taihangensis. Significant genetic divergence and a highly structure were illustrated in all studied populations (Figure 3, Figure 4, and Table 2). For O. longilobus (F STcpDNA = 0.607, F STnrDNA = 0.682) and O. taihangensis (F STcpDNA = 0.686, F STnrDNA = 0.805), the results of our analysis suggested that significant variation was detected from among populations rather than from variation within groups (Table 2). Moreover, very high genetic differentiation in O. longilobus and O. taihangensis greatly exceeding the mean value (G ST = 0.22) reported by Nybom [61].
Species genetic diversity and structure can be affected by life history traits (e.g., life cycle, breeding system, pollination mechanism) and environmental effects (e.g., geographical range, climate, topography) [61][62]. O. longilobus and O. taihangensis possess a sexual reproductive mode and a relatively long life cycle as perennial herbs, all of which are traits associated with low total genetic diversity [63][64]. However, high genetic diversity was observed in our study, and this diversity appears to result from high among-population genetic variation. It has been widely documented that the dispersal of pollen and seeds in plant species is strongly linked to the development of a population genetic structure in bi-parentally inherited (nuclear-encoded) and maternally inherited (cytoplasmic) genetic markers [12,31,61,62,64]. Generally, outcrossing taxa with high seed dispersal capacity retain the majority of both types of genetic marker variation within populations, whereas selfing (and/or asexual) taxa with restricted seed dispersal allocate the majority of such variation among populations. O. longilobus and O. taihangensis are self-compatible plant species [25]. These two species, therefore, might be expected to have a heterogeneous distribution of both cpDNA and nrDNA variation among populations. Taken at face value, the above results may suggest that restricted gene flow via pollen and seed has resulted in significant population genetic differentiation in O. longilobus (N mcpDNA = 0.12; N mnrDNA = 0.10) and O. taihangensis (N mcpDNA = 0.19; N mnrDNA = 0.14). The physiographic complexity of Taihang Mountains and the deeply carved valleys and ravines between the inhabit areas probably would imposes significant barriers to gene flow among populations of O. longilobus and O. taihangensis. The genetic divergence between the two lineages occurred at 1.28 Ma, corresponding to the early Pleistocene. This divergence falls within the glaciations during the Taku Glaciation (1.05-1.20 Ma). The glaciers and/or extremely low temperature in the mountains might also have created barriers to gene flow between geographically isolated populations, which therefore promoted the divergence. For both marker systems, we generally observed a good fit of these populations to an IBD model. When coupled with high levels of population subdivision, such conditions strongly suggest that geographical isolation had a larger historical role in extant population structure compared with limited pollen and seed flow alone. The habitats of O. longilobus and O. taihangensis are either rocky or massifs, and these species grow discontinuously in different habitats along the altitudinal gradient from 700 to 1500 m, in which the variable micro-surrounding, Table 3. Parameters of mismatch distribution analysis. complex topography, and great altitudinal variability of the region might promote the high degree of genetic population differentiation [65][66]. Limited migration coupled with a heterogeneous environment could promote local adaptation and fixation of different alleles in O. longilobus and O. taihangensis populations.
Phylogeographical structure and inference of demographic history Both O. longilobus and O. taihangensis populations have a significant phylogeographic structure. The current distribution of haplotypes could be a result of past climatic changes related to the advance/retreat of glaciations [11]. During the Quaternary periods, climatic oscillations resulted in repeated drastic environmental changes, which further caused massive range shifts of most plants and animals, leading to accumulated genetic differences and particular phylogeographic patterns [2][3].
Our analyses suggested that Opisthopappus originated during the Pliocene. Before the Pleistocene, the current range of Taihang Mountains would have been covered by grassland with roughly similar geographic conditions [67]. As herbs, O. longilobus and O. taihangensis can disperse following grassland migration. On this basis, we suspect that O. longilobus and O. taihangensis were once distributed throughout Taihang Mountains and its adjacent region. With uplift of Taihang Mountains [27] and climate fluctuation of Quaternary glacial periods, O. longilobus and O. taihangensis were fragmented and forced into the refuge areas. The flora would migrate and expand during interglacial and glacial periods. This migration and expansion resulted in ecological displacement, which might have resulted in populations situated at different spatial-temporal scales. Different habitats might enhance the isolation of plant populations and prevent gene flow, which in turn could lead to incipient allopatric speciation or further allopatric speciation [12,18,19,68,69]. During the Quaternary periods, the paleovegetation in Taihang Mountains repeatedly appear the replacement of grasslands and forests cycles [67]. Those multiple replacements likely fragmented and isolated the habitat of herbs. Therefore, the refuge areas for O. longilobus and O. taihangensis would have remained separated and fragmented at different spatial-temporal scales for a long time. Moreover, dramatic geological changes induced by the Taihang Mountains uplift contributed to the extremely complicated landscape, which likely acted as a major driving force for the separation of the O. longilobus and O. taihangensis. Divergent selection between populations in contrasting environments, longer temporal and spatial separation and fragmentation would accelerate their differentiation, and eventually promote allopatric speciation.
Following  Table 2). The 10 bp insertions of ndhJ-trnL occurring in population LLS also supported this viewpoint (Table  S1). Subsequently, with the advent of inter-glaciations or postglaciations, the climate changed, which led to the paleovegetation of Taihang Mountains experiencing a transition from grasslands to forests. The corridor for O. longilobus migration gradually disappeared, which would have resulted in isolation of the population [11,12,68]. Additionally, we did not detect any population expansion due to the clear multimodel mismatch distribution ( Figure S2). The complex geological conditions of Taihang Mountains would limit O. longilobus population expansion on a macro scale and thus restricting O. longilobus to sustain in situ during the Quaternary periods, perhaps by moving upwards and downwards in their mountain ranges. For O. taihangensis, a pattern of population expansion was observed, because the observed mismatch distribution for cpDNA and nrDNA haplotypes was unimodal (a pattern consistent with expansion) ( Figure S2). This expansion was dated to be 0.78 Ma, which corresponding the interglacial stage between the Taku Glaciation and Lushan Glaciation. The above results suggest an expansion consistent with a range expansion and re-colonization before the mainly uplift of Taihang Mountains. With the uplift of Taihang Mountains, large geographical distances and significant topographical barriers could restrict gene flow and limit any large-scale population expansion [12,70]. The cpDNA sequences characteristic of the trnL-trnF and ndhJ-trnL intergenic spacers further confirmed this hypothesis (Table S1).
Crandall and Templeton [71] suggest that ancestral haplotypes are located in a central position within the haplotype phylogeographical network and that potential refugia are usually characterized by widely dispersed ancestral haplotypes and other tip and unique haplotypes. Geographic areas displaying increased levels of genetic diversity are thus good candidates in the search for past refuges [8,34]. These regions are characterized by relatively stable ecological conditions during environmental fluctuations, which foster the accumulation of genetic diversity [72]. Most O. longilobus individuals or populations gathered in lineage I (Figure 3 and Figure 4). Haplotypes H3 and Hn8 are located in the center of the haplotype network. And Haplotype H3 is detected in three populations (WDS, SHS and WML), Hn8 in two populations (WDS and SHS) (Figure 2). WDS and SHS are located in the southern part of Hebei province and have higher haplotype diversity (Table 1). Thus, we suggest that the WDS and SHS populations' regions were two major potential refugia for O. taihangensis. These regions may therefore have been characterized by a relatively mild climate, potentially containing microclimatic environments able to impart relative stability to a range of habitats. Lineage II consisted of all populations of O. taihangensis and some individuals of O. longilobus. Haplotypes H1 and Hn1, which are located in the center of the haplotype network, were the most dominant haplotypes found in five and ten populations, respectively. Moreover, some populations in these haplotypes (H1 and Hn1) had higher haplotype diversity (Table 1, e.g., BJY, WML, and YTS). Thus, we suggest that there were multiple small refuges located in Shanxi and Henan provinces during the last glacial maximum or earlier cold periods for O. taihangensis. The above populations of O. taihangensis, distributed in southern Taihang Mountains, would have been provided with a complex and stable environment. Thus, the potential refugia area may have been located in those regions.
However, our comprehensive sampling showed that assessing the species-level relationships in our study is complicated by the fact that not all cpDNA and nrDNA haplotypes are speciesspecific. Incongruence between the species boundaries and the genealogy of their cpDNA and nrDNA sequences is illustrated by the sharing of haplotypes (e.g., H1 and Hn8) between the two species. On the other hand, several haplotypes were found within a single species (Table 1 and Figure 2). Haplotype distribution does not strictly follow species circumscription. Moreover, an association of haplotypes with geographically circumscribed regions rather than with taxonomic boundaries is a phenomenon observed in O. longilobus and O. taihangensis. The sharing of haplotypes and lack of reciprocal monophyly might be explained by the persistence of ancestral polymorphisms during speciation events and/or exchange of genes by inter-specific hybridization (Table S1 and Figure S1). These two species have similar habitats and the same period of flowering [25], and possible visitation by the same pollinators in the absence of physiological or genetic reproductive barriers might have enabled hybridization events. Nonetheless, further work on chromosomes and the sequencing of additional genomic regions are needed to make more in-depth conclusions about the hybrid status of these taxa.

Conclusions
In summary, O. longilobus and O. taihangensis have a high level of genetic diversity. The following phylogeographical structure and historical scenario of O. longilobus and O. taihangensis can also be drawn. The vicariance following the uplift of Taihang Mountains and a transition of paleovegetation of Taihang Mountains mainly shaped the present distribution of these two species. The above findings imply that multiple small refugia could have been maintained in the Taihang Mountains regions, which allowed O. longilobus and O. taihangensis to persist in situ and maintain sizable haplotype and nucleotide diversity at the lineagewide scale during the glaciations. The current distribution pattern of O. longilobus and O. taihangensis refuges is similar to that of other East Asian plants (e.g., Cathaya argyrophylla, Picea crassifolia, Sinopodophyllum hexandrum, Lagochilus Bunge ex Bentham) [8,20,73], with multiple refuges sustained across their distribution ranges. Furthermore, complex topography rendered these refuge areas both isolated and fragmented in different geographic units. Geographical and ecological isolation likely restricted seed and pollen migration ability. The isolation and fragmentation further promoted the intra-specific split. These results support the conclusion that the high plant diversity and endemism found in central/northern China and throughout eastern Asia has mainly resulted from allopatric speciation due to the complex topography of Mountains and allopatric fragmentation during the late Tertiary and Quaternary periods.