Effects of Geological and Environmental Events on the Diversity and Genetic Divergence of Four Closely Related Pines: Pinus koraiensis, P. armandii, P. griffithii, and P. pumila.

The effects of mountain uplift and environmental oscillations on nucleotide variability and species divergence remain largely unknown in East Asia. In this study, based on multiple nuclear DNA markers, we investigated the levels and patterns of nucleotide diversity and interspecific divergence in four closely related pines in China, i.e., Pinus koraiensis, P. armandii, P. griffithii, and P. pumila. The four pine taxa shared low levels of nucleotide polymorphisms at the species level. P. pumila had the highest silent nucleotide diversity (πsil = 0.00661) whereas P. griffithii had the lowest (πsil = 0.00175), while the levels of genetic polymorphism in P. armandii (πsil = 0.00508) and P. koraiensis (πsil = 0.00652) were intermediate between the other two species. Population genetic structure analysis showed that variations primarily existed within populations of the four pine species, presumably due to habitat fragmentation or the island-like distributions of Pinus species. Population divergence (FST) analysis showed that the genetic divergence between P. griffithii and P. koraiensis was much greater than that between P. koraiensis and the other two pines species. Isolation-with-migration analysis suggested that asymmetric gene flow had occurred between any two pairs of pine species. Phylogenetic analyses indicated that the four allied species split into two groups about 1.37 million years ago, where P. armandii and P. pumila were closer and clustered as sister species, whereas P. koraiensis and P. griffithii were clustered on another branch. Our results and those obtained in previous studies suggest that mountain uplift and geological climate oscillations may have led to the patterns of genetic divergence and nucleotide variations in these four pine species.

Introduction

Nucleotide diversity levels within populations and spatial patterns, as well as species divergence are of great importance in the field of evolutionary biology (Coyne and Orr, 2004; Hao et al., 2015; Ortego et al., 2015). Mountain uplift and past environmental oscillations may have been largely responsible for shaping the spatial patterns of diversity and genetic divergence among species (Coyne and Orr, 2004; Wachowiak et al., 2009). In general, the level and distribution of nucleotide diversity are historical products of the long-term evolution of a species, and they are largely associated with the evolutionary potential or future fate of a species (Wright and Gaut, 2005; Wachowiak et al., 2009; Zhou et al., 2014; Tsuda et al., 2017). In addition, ecological or proximal causes (e.g., mating systems) and various barriers (e.g., geographic and spatio-temporal isolation) due to geological history can cause fragmented of species distributions, which may lead to reduced gene flow between isolated populations and adaptability. This process initiates allopatric divergence, and local adaptation can ultimately drive populations toward speciation and change evolutionary processes (Cutter and Gray, 2016; Ren et al., 2017).

Conifers are anemophilous and outcrossing (Fu et al., 1999). They are mainly characterized by long life cycles, large effective population sizes, incomplete lineage sorting, and extensive introgression/hybridization among populations, which makes their genetic structure and spatial patterns of diversity very different from those found in traditional model plants (Neale and Kremer, 2011; Gao et al., 2012; Li et al., 2012, 2013; Hao et al., 2015). For instance, conifers tend to share haplotypes/genotypes among species, with no distinct genetic divergence across species ranges, and most of the genetic variations are found within populations (Willyard et al., 2009; Chen et al., 2010; Ren et al., 2012; Liu et al., 2014; Zhou et al., 2014). In recent years, many studies have determined nucleotide polymorphisms and speciation history patterns using multiple nuclear loci (Ma et al., 2006; Li et al., 2010; Gao et al., 2012; Wachowiak et al., 2013; Zhou et al., 2014; Zou et al., 2016). These biparentally inherited nuclear genes are functional genes that encode proteins and they are characterized by their orthology, moderate to high rates of evolution, and the presence of many phylogenetically informative sites (Zhou et al., 2014). Therefore, large numbers of nuclear markers can be used to detect the deep evolutionary relationships among closely related species, especially recently diverged taxa (Chen et al., 2010; Zou et al., 2016). In this study, we employed nucleotide polymorphisms as well as the population structure and speciation history to explore the relationships among four Pinus species.

Four related Pinus species in subsection Strobus occur in East Asia: P. armandii Franch., P. koraiensis Sieb. et Zucc., P. griffithii McClelland, and P. pumila (Pall.) Regel. These species share some common features, such as possessing five needle leaves in a bundle. There are obvious differences among these species in terms of their ecological niche, natural geographic distribution, morphology, wood anatomy, and cytology (Wu and Feng, 1995). P. armandii, P. pumila, and P. koraiensis occur according to the changes in the hydrothermal conditions, and P. griffithii is distributed on the China–Nepal and China–Bhutan borders (Supplementary Figure ; Wu and Feng, 1995). The distributions of these pines also increase successively from low to high altitudes, where P. koraiensis occurs at altitudes of 150–1,800 m and P. pumila always forms copses with other coniferous trees on mountain tops at altitudes of 1,000–2,300 m. P. armandii usually occurs in pure forest or mixed forest at altitudes of 1,000–3,300 m, and P. griffithii is distributed in the same manner at altitudes of 1,600–3,300 m on the Qinghai–Tibet Plateau and Mount Everest. According to the classic categorization of Pinus sect. Strobus, P. pumila is categorized into the P. koraiensis taxon and P. armandii into the P. griffithii group (Fu et al., 1999). A study of the divergence of the resin ducts in Pinus sect. Strobus suggested that P. armandii is the most primitive species and its southward spread gave rise to P. griffithii, whereas its northward spread gave rise to P. pumila and P. griffithii (Peng, 1999). In recent years, several studies based on plastid molecular markers have shown that P. koraiensis and P. pumila are most closely related to each other (Peng, 1999; Wang et al., 2016). However, the accurate phylogenetic relationships and interspecific divergence among these related pine species is still controversial due to the limited availability of morphological and molecular biological evidence (Peng, 1999; Liu et al., 2014; Hao et al., 2015).

In addition, some studies found low variations in DNA barcodes, such as rbcL and matK, among related species due to low levels of cpDNA diversity and genetic divergence (Syring et al., 2007; Li et al., 2015). Moreover, frequent interspecific introgression and hybridization among species have important effects on their genetic diversity levels, especially in parapatric or allopatric of species in China. Thus, in the current study, we used multiple nuclear genes to investigate the genetic diversity and divergence in four closely related pine species comprising P. koraiensis, P. armandii, P. griffithii, and P. pumila. We specifically addressed the following two questions. (1) How is the level and pattern of population divergence among the four pines species? (2) How is the pattern of gene flow and interspecific introgression between these closely related species in East Asia?

Materials and Methods

Population Sampling

To accurately determine the nucleotide diversity and interspecific relationships among pines, we sampled 216 individuals from 16 allopatric populations of the four pine species (Figure ). The distance between any two trees of the same species was at least 50 m (Supplementary Table ). We isolated the haploid megagametophyte from each of the sampled trees.

(A) Geographical distribution of the sampled populations of the four closely related species: Pinus pumila (red), P. griffithii (green), P. koraiensis (yellow), and P. armandii (blue). Color scales indicate different altitudes. (B) Bayesian clustering analysis to determine the population structure of the four pine species. Red, green, and blue represent the dominant clusters (K = 3) identified by STRUCTURE in each population.

DNA Extraction, PCR Amplification, and Sequencing

Total genomic DNA was extracted from the megagametophyte samples for each individual using the modified CTAB method (Doyle and Doyle, 1987). In preliminary studies, about 40 nuclear gene loci were screened for cross-amplification in the four pine species (Ma et al., 2006; Eckert et al., 2013). Finally, six polymorphic loci (1_1609_01, CL1694, PTIFG2009, 0_12929_02, 0_14221_01, and 0_1688_02) associated with protein kinase family protein, serine-tRNA ligase, and leucine-rich repeat family protein were selected for subsequent sequence amplification and analysis (Supplementary Table ). PCR amplification was conducted in a volume of 25 μL with a DNA concentration of 10–40 ng/μL, 50 mM of Tris-HCl, 0.5 mM of each dNTP, l.5 mM of MgCl2, 2 μM of each primer, and 0.75 U of Ex Taq DNA polymerase (Runde, Xi’an, China). The PCR program comprised initial denaturation at 94°C for 5 min, followed by 35 cycles for 1 min at 94°C, at a specific annealing temperature (53–60°C, see Supplementary Table for details) for 1.5 min and extension for 1 min at 72°C, and a final extension for 10 min at 72°C. Primer synthesis and sequencing of the PCR products were performed by Shanghai Biological Engineering Co. Ltd. (Shanghai, China). Sequencing was conducted using both forward and reverse primers for each gene (Supplementary Table ) on an ABI Prism 3730xl sequencer (Applied Biosystems, Foster City, CA, United States).

Data Analysis

Data Reconciliation

Sequences were aligned and manually adjusted with Chromas and MEGA5.0 (Librado and Rozas, 2009; Tamura et al., 2011) to correct random errors generated by sequencing.

Nucleotide Diversity and Neutral Tests

The genetic diversity parameters for the four pine species were calculated using DnaSP v. 5.10 (Librado and Rozas, 2009), including the number of segregation sites S, Watterson parameters θw (Watterson, 1975), total nucleotide polymorphisms πt (Li and Nei, 1975), nucleotide diversities of non-synonymous sites and silent loci (synonymous sites and non-coding positions), πa and πsil, number of haplotypes Nh and haplotype diversity Hd (Nei and Tajima, 1981; Fu, 1997; Depaulis and Veuille, 1998; Depaulis et al., 2001), and intragenic minimum recombination events (RM) (Hudson and Kaplan, 1985). In addition, in order to accurately detect departure from the neutral model of molecular evolution at each locus, the neutral equilibrium was tested for various parameters using Tajima’s D (Tajima, 1989), Fu and Li’s D and F (Fu and Li, 1993), and the standardized Fay and Wu’s H (Fay and Wu, 2000). Tajima’s D measures the standardized difference between π and θW, whereas Fay and Wu’s H measures the difference between π and θH. The former is more sensitive to an excess of rare variants whereas the latter is more sensitive to an excess of high-frequency-derived variants. Both D and H are expected to be zero under the standard neutral model (Zou et al., 2013). We also conducted maximum frequency of derived mutations (MFDM) tests to examine the likelihood of natural selection acting on individual loci at species levels. The MFDM tests exclude the confounding effects of demography completely when detecting recent positive selection (Li, 2011). In practice, a single DNA fragment (i.e., a locus) may have a short length and only contain a few RM. The MFDM v. 1.1 test always depends on the estimate of RM (Li, 2011).

Genetic Divergence and Population Structure

The sources of genetic variation among the four species (group), populations and individuals were analyzed by analysis of molecular variance (AMOVA) with ARLEQUIN v.3.11 (Excoffier et al., 2005). We estimated F statistics hierarchically, both among species (FCT) and among populations within species (FST). FST (Wright, 1949) is a coancestry statistic that provides the variance within populations relative to the total population. We used NETWORK v. 4.6.1.3 (Bandelt et al., 1999) to construct phylogenetic relationships based on the haplotypes of each species at the six loci (gaps were excluded). In addition, genetic clustering based on individuals was estimated by Bayesian clustering using the STRUCTURE V.2.3 program (Hubisz et al., 2009). To estimate the number of clusters (K) in the data, K values from 1 to 16 were explored using 10 independent runs per K and an admixture model. As described in previous studies, in order to generate a reliable estimate of the optimal K, the burn-in was set to at least 200,000 and the run length was at least 500,000 (e.g., Zou et al., 2013; Zhou et al., 2014; Tsuda et al., 2015). We also utilized the STRUCTURE HARVESTER program to estimate the most likely number (K) of genetic clusters (Evanno et al., 2005; Earl and vonHoldt, 2012).

Reconstruction of Historical Dynamics and Species Relationship

The migration rates, effective population sizes, and population split times were calculated based on the isolation-with-migration (IM) model using the IMa2 program to infer the population history dynamics of the four pine species (Nielsen and Wakeley, 2001; Hey, 2006, 2010; Kuhner, 2009). We analyzed sibling species in a pairwise manner using a basic two-population model. We extracted the largest region with no recombination for each of the six nuclear loci. Functions of the model parameters were estimated in the M-mode based on 1 × 106 Monte Carlo Markov chain (MCMC) steps following 5 × 105 burn-in periods in order to obtain reliable estimates (i.e., similar posterior distributions for the parameter), and the effective sample size for each parameter was at least 200. The divergence time between species was estimated based on a mean mutation rate of μ = 4.875 × 10−9 (per site per generation), and the generation time for pines was assumed to be 25 years (Ma et al., 2006). In addition, we constructed the phylogenetic relationships among the four pine species using ∗BEASTv1.8.0 (Heled and Drummond, 2010). The species tree was computed using the six nuclear genes sequenced for the sampled species. We selected a Yule model as the species tree prior, a constant population size, and relaxed lognormal clock models for all nuclear loci (Heled, 2012). Pinus bungeana was used as outgroup. We ran the MCMC analysis for one billion generations with sampling every 50,000 generations. Two independent runs were conducted. Tracer[1] v1.5 (Rambaut and Drummond, 2009) was used to assess the convergence of chains to the stationary distribution (effective sample size >200). After discarding the first 2,500 trees as a burn-in, the remaining trees were summarized in a maximum clade credibility tree with the TreeAnnotator v1.8.0 program (Drummond and Rambaut, 2007). Joint Bayesian species delimitation and species tree estimation were also analyzed using the BPP v3.4 program based on the multispecies coalescent model (Yang, 2015). In addition, phylogenetic relationships based on nuclear haplotypes were reconstructed with the maximum likelihood (ML) model using PAUP 4.0 (Swofford, 2002). The ML analysis employed the HKY substitution model, where support values for the nodes were estimated based on 1,000 bootstrap replicates.

Results

Nucleotide Polymorphisms

For all nuclear loci, P. griffithii had the lowest estimates for the total average segregating sites and average values of the segregating sites in silent sites compared with the other three species (Table ). P. griffithii had two singleton mutation sites in PTIFG2009. The numbers of shared polymorphisms (SS) were similar among the four closely related species and the numbers of fixed differences (Sf) were low. The differences in the polymorphisms between P. pumila and P. griffithii were mainly due to the 0_14221_01 gene locus, with 18 polymorphic sites in P. pumila but only three in P. griffithii (Supplementary Table ). Similarly, the differences between P. koraiensis and P. pumila were mainly due to the CL1694 locus, with 20 polymorphic sites in the former but only five in the latter. The 0_14221_01 and 0_12929_02 loci accounted for the observed differences between P. armandii and P. pumila.

Neutrality Tests

Positive Fu and Li’s D∗ and Fu and Li’s F∗ values were estimated for most loci, although most of these values were not significant in each species (Table ). The mean Tajima’s D (D) values were negative for P. pumila (−0.236) and P. griffithii (−0.030), but positive for P. koraiensis (0.274) and P. armandii (0.254) (Table ). In addition, the mean Fay and Wu’s H (H) values were negative for P. koraiensis and P. griffithii but positive for P. pumila and P. armandii (Table ). However, with the exception of locus PtIFG2009 (P = 0.04674) in P. pumila, no significant deviation from neutrality were detected for the six loci using the MFDM test (Supplementary Table ). The MFDM test detected slight deviation from the standard neutral model at the PtIFG2009 locus in the four pine species by considering genetic recombination (P < 0.05).

Population Genetic Structure

Within the four species, the 0_12929_02 locus had the highest genetic divergence among populations (FST = 0.714, P < 0.001), whereas the 1_1609_01 locus had the lowest genetic divergence among populations (FST = 0.085, P < 0.001) (Supplementary Table ). The population genetic divergence was also significant across all loci (FST = 0.624, P < 0.001). It should be noted that FST was much higher for the overall loci than interspecific genetic differentiation (FCT) except for the 0_12929_02 locus (Supplementary Table ). Between pairs of species, FST varied from 0.008 to 1.000 (Supplementary Table ).

The divergence among the four pine species at the six nuclear loci was also supported by Bayesian clustering analysis (Figure ). The most likely number of clusters for the entire dataset was K = 3 (Supplementary Figure ). P. pumila and P armandii individuals were separated into two groups that corresponded to their respective species, whereas the majority of P. koraiensis and P. griffithii individuals were assigned to another cluster. Remarkable levels of gene flow and gene introgression were apparent between P. koraiensis and P. armandii (Figure and Supplementary Figure ).

Genealogy of Each Locus

The average number of haplotypes (Nh) and haplotype diversity (Hd) were much higher in P. pumila (Nh = 12.167, Hd = 0.822) and P. armandii (Nh = 10, Hd = 0.781) than P. koraiensis (Nh = 7.167, Hd = 0.655) and P. griffithii (Nh = 4.167, Hd = 0.420) (Table ). The haplotypes in the center of the network were shared (Figure ), but most haplotypes were exclusive to specific species at the five loci (1_1609_01, 0_1688_02, PTIFG2009, 0_14221_01, and CL1694). In addition, there was no shared haplotype at the 0_12929_02 locus. According to the ΦST values for all loci among the species in Supplementary Table , each two Pinus species exhibited significant genetic divergence, where the highest ΦST value (0.62006) was found between P. koraiensis and P. griffithii, whereas the divergence between P. pumila and P. griffithii was lowest (ΦST = 0.18941).

Networks obtained for the six nuclear genes in the four species comprising Pinus pumila (red), P. griffithii (green), P. koraiensis (yellow), and P. armandii (blue). Each sector of a circle corresponds to the frequency of the haplotype for each species.

Evolutionary Relationships Among the Four Species

The mean divergence time between P. pumila and P. armandii was estimated at 1.13 million years ago (Mya). A younger divergence time (0.319 Mya) was estimated between P. griffithii and P. koraiensis (Table ). In addition, we found asymmetric historical gene flow between pairs of species. In particular, the migration rate from P. pumila to P. griffithii was 2.0450, with 0.0005 in the reverse direction (Table ). Pinus griffithii and P. koraiensis had smaller population sizes (0.0918–0.1846 and 0.3264–0.5369, respectively; Table ) than P. pumila and P. armandii (0.3738–0.8326 and 0.4087–0.7019, respectively; Table ). The species tree analyses demonstrated that the relationship was closer between P. armandii and P. pumila where they clustered as sister groups, whereas P. koraiensis and P. griffithii were located in another clade. The four species split into two groups about 1.37 Mya (Figure ). Moreover, we obtained the best species-tree model using BPP v3.4 software, where the posterior probability of the species tree was one and the acceptance proportion was near to zero (0.025) based on multiple runs (Supplementary Figure ). The topology of the tree was consistent with the results obtained by ∗BEAST (Figure and Supplementary Figure ).

Divergence times and phylogenetic relationships among four pine species. The tree was constructed based on six nuclear genes using ∗BEAST. The tree was rooted with Pinus bungeana. The numbers on the branches indicate the corresponding posterior probabilities values, mean divergence dates, and 95% credibility interval.

Discussion

Nucleotide Diversity

The nucleotide diversity at silent sites basically agreed with the neutral model of molecular evolution (Ma et al., 2006; Wachowiak et al., 2016). We detected low levels of silent polymorphisms in the four closely related pines species, because the average values were much lower than the average polymorphism for most conifers at multiple nuclear genes (πsil = 0.0029–0.0122) (Ma et al., 2006; Li et al., 2012). Among the four species, P. pumila had the highest silent nucleotide diversity (πsil = 0.00661), whereas P. griffithii had the lowest (πsil = 0.00175). However, these diversity values were much lower than those found in other Pinus species, such as P. densata and P. yunnanensis (Ma et al., 2006). Factors such as the nuclear gene loci selected in the study, sample size variations, mutation rates within species, demographic effects, and natural selection can influence the nucleotide diversity levels and patterns in species or populations (Loveless and Hamrick, 1984; Lande, 1988; Li et al., 2012; Zou et al., 2013). The four related species and six nuclear loci investigated in our research have also been studied previously, and thus these nuclear loci were not the main cause of the low levels of nucleotide polymorphisms. However, the unequal sample sizes for different genes and species may have caused differences in the nucleotide variability among species. To verify this bias, we detected the nucleotide diversity parameters based on the same sample sizes for each gene from each Pinus species. The results showed that there were significant differences in the levels of diversity among different species compared with the previous estimates (Table and Supplementary Table ). We concluded that the levels of nucleotide variability among species were significantly associated with the samples sizes of the pine species. Similar differences in the patterns of diversity have also been detected in some other gymnosperm species (Ma et al., 2006; Du et al., 2009; Zou et al., 2013; Zhou et al., 2014). In addition, the long life cycles and low mutation rates in conifers may explain the low levels of nucleotide polymorphisms in P. pumila, P. griffithii, P. koraiensis, and P. armandii. Moreover, high levels of linkage disequilibrium were detected and some species deviated from neutrality according to the tests conducted in our study. This was particularly evident at one locus, i.e., PtIFG2009, which suggests that this locus might have undergone selection or population shrinkage according to the results obtained from the Tajima’s D and MFDM tests. In particular, for the populations of P. griffithii and P. koraiensis, Tajima’s D was positive, and Fay and Wu’s H was negative, which is a pattern that is consistent with a recent bottleneck. In addition, P. griffithii and P. pumila descendant populations had a somewhat smaller size than the ancestral population (Table ), and thus it is possible that the populations have experienced from genetic bottlenecks. The population dynamics due to geological isolation and climatic oscillations probably contributed to their relatively lower diversity (Chen et al., 2017). In addition, the mean Tajima’s D (D) and mean Fay and Wu’s H (H) values were negative but close to zero for P. armandii at the PtIFG2009 locus, which may indicate a neutral equilibrium (Holliday et al., 2010). The numbers of nucleotide polymorphism were higher in P. pumila, P. koraiensis, and P. armandii than P. griffithii. These results can partly be explained by the fact that their seeds are food for nutcrackers and rodents such as squirrels. These animals may screen the seeds and transport them over long or short distances for secondary storage and dispersal, thereby also enhancing the spread of the seeds, and this may affect their genetic differences (Li et al., 2007; Fan and Jin, 2011). The low level of nucleotide polymorphism in P. griffithii may be explained by its small geographic distribution compared with more common and widespread species, because of drift, founder events, and other stochastic processes (Cole, 2003; Chen et al., 2017). In addition, P. pumila is well known because of its larger island-like distribution and it rarely develops into pure forest, thereby accounting for its higher diversity compared with the other three species (Qiu et al., 2007; Chen et al., 2017).

Interspecific Gene Flow and Species Divergence

Analysis of molecular variance detected remarkable divergence in the four Pinus species where the variations were mainly within populations, which agreed with the small differences among populations of wind-pollinated gymnosperms (Supplementary Table ). However, we also found a high level of genetic divergence within groups, possibly due to habitat fragmentation or the island-like distribution of the four species, particularly when considering that the habitats of P. pumila and P. griffithii are harsher than those of the other two species. The former often grows in barren soil on bare rocky peaks. This type of habitat is vulnerable to fragmentation but mountains and ravines may partly hinder the gene flow between populations, thereby leading to the isolation of groups. Genetic divergence was found among the four species, although some degree of gene flow and introgression was detected. In particular, populations of P. koraiensis had a mosaic-like pattern and they were further subdivided into independent sub-clusters when K = 4, which suggests that a high level of introgression in this species. The pairwise migration rate between P. koraiensis and P. pumila was relatively high compared with that between P. koraiensis and P. armandii. In addition, P. pumila and P. koraiensis had a relatively limited distribution in the northwest and northeast of North China, and there was no clear phylogenetic resolution among P. pumila, P. koraiensis, and P. armandii based on DNA fragments from the chloroplast, mitochondrial, and nuclear genomes according to previous phylogenetic studies (Liu et al., 2014; Wang and Wang, 2014; Hao et al., 2015). These results suggest that migration may have led to a sympatric distribution in addition to the existing incomplete reproductive isolation. The phylogenetic relationships determined based on ML and NETWORK analysis also showed that the shared haplotypes were located in the center of the topological structure (Figure and Supplementary Figure ), and thus incomplete genealogical screening based on a large effective population of Pinus may have led to the sharing of ancestral polymorphisms. However, there were no shared haplotypes based on the 0_12929_02 locus, and the interspecies variation was similar to the genetic differentiation among populations (FCT = 0.718, FST = 0.714; P < 0.001; Supplementary Table ). In general, different nuclear loci have different evolutionary rates and molecular functions (Li et al., 2012; Eckert et al., 2013; Zhou et al., 2014). Previous studies have shown that the 0_12929_02 locus is associated with the protein kinase family and that it has been under selection (Eckert et al., 2013). The rapid fixation of genetic variation in this locus may have led to greater species divergence (Nei et al., 1975; Nei and Tajima, 1981; Ellstrand and Elam, 1993). In addition, nuclear DNA introgression in ancestral populations among species may have also affected the topology of the phylogenetic trees. The significant topological incongruence among the nuclear gene trees (Supplementary Figure ) indicates a complex evolutionary history, thereby providing novel insights into the evolution of Pinus. The four species split into two clades about 1.37 Mya (Figure ). However, we should be cautious when inferring divergence times based on the assumption that the generation time is 25 years in the four Pinus species because of the longevity of Pinus, the long overlaps between generations, the variable age of maturity and the replacement speed of forests. In addition, our multilocus analysis determined that the genetic divergence among the four pine species, was consistent with geological events and climatic oscillations in the mid- to late Tertiary period about 5 Mya. The uplift of the Tibetan Plateau caused by Himalayan orogeny had a great impact on the climate in China, with decreases of in temperature in some areas, but increases of 4–8°C in the region east of 100°E (across the Inner Mongolia, Gansu, Qinghai, Sichuan, and Yunnan regions of China) and of 1–4°C to the west (Jiang, 2009). These climatic conditions may have changed the geographic distributions of plants, and thus we suggest that P. armandii and its ancestral population spread eastward to the northeast of China and westward to Tibet. However, the intensities of the winter and summer monsoons were reduced greatly during the middle-late Pliocene, and the dispersal of Pinaceae pollen by the wind might have been affected (Zhou et al., 2014). Moreover, after gradual changes in the microhabitats and variations in the directions and amounts of gene flow, as well as the accumulation of mutations, new relatives may have emerged by gradual divergence. Effective migration, hybridization, and introgression among species can increase genetic diversity (Wachowiak et al., 2016), and other factors such as selection, isolation, and genetic drift among different microhabitats can promote divergence and speciation. Indeed, significant and asymmetric gene flow and introgression were detected in the four closely related Pinus species. Gene flow and genetic introgression among different pines could have led to changes in genetic variability (Hao et al., 2015; Wachowiak et al., 2016).

Data Availability

Sequence data obtained in this study were deposited in GenBank (KF286539 – KF286739).

Author Contributions

Z-HL conceived the study. YJ and JZ performed the experiments. Z-HL, YJ, YW, W-BF, and G-FZ contributed materials and analysis tools. Z-HL, YJ, and JZ wrote the manuscript. YJ and ZL revised the manuscript. All authors approved the final version of the manuscript.

Conflict of Interest Statement

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.

Table 1

Nucleotide variations in four Pinus species: Pinus pumila, P. griffithii, P. koraiensis, and P. armandii.

Species	Locus	Total					Nonsynonymous sites				Silent Sites				Rm
		N	L	S (Singl.)	θ	πt	L	S	θW	πa	L	S	θW	πSil
P. pumila	1_1609_01	80	387	14 (0)	0.00850	0.00827	101	3	0.00700	0.00916	282	11	0.00916	0.00807	1
	PTIFG2009	88	505	13 (1)	0.00592	0.00648	194	4	0.00474	0.00600	311	9	0.00666	0.00678	7
	0_1688_02	80	613	6 (0)	0.00233	0.00258	217	0	0.00000	0.00000	387	6	0.00364	0.00404	1
	CL1694	56	305	5 (2)	0.00421	0.00330	59	0	0.00000	0.00000	245	5	0.00524	0.00410	0
	0_12929_02	100	735	14 (1)	0.00425	0.00255	288	2	0.00155	0.00078	443	12	0.00604	0.00372	2
	0_14221_01	68	585	18 (1)	0.00801	0.00779	289	4	0.00391	0.00359	285	15	0.01315	0.01297	4
Average		78		11.67	0.00554	0.00516	191.3	2.2	0.00287	0.00326	325.5	9.7	0.00732	0.00661	15
P. griffithii	1_1609_01	56	387	3 (0)	0.00200	0.00162	101	2	0.00509	0.00483	280	1	0.00092	0.00049	0
	PTIFG2009	52	505	9 (2)	0.00467	0.00568	194	3	0.00405	0.00541	311	6	0.00506	0.00585	4
	0_1688_02	52	613	10 (0)	0.00432	0.00291	217	6	0.00726	0.00409	387	4	0.00271	0.00226	0
	CL1694	16	305	0 (0)	0.00000	0.00000	59	0	0.00000	0.00000	245	0	0.00000	0.00000	–
	0_12929_02	40	735	0 (0)	0.00000	0.00000	290	0	0.00000	0.00000	441	0	0.00000	0.00000	–
	0_14221_01	60	585	3 (0)	0.00140	0.00180	289	1	0.00101	0.00174	285	2	0.00178	0.00188	1
Average		46		4.17	0.00207	0.00200	191.7	2	0.00290	0.00268	324.8	2.2	0.00175	0.00175	5
P. koraiensis	1_1609_01	24	387	7 (0)	0.00599	0.00775	101	2	0.00659	0.01011	282	5	0.00586	0.00702	0
	PTIFG2009	52	505	11 (2)	0.00571	0.00684	194	2	0.00271	0.00394	311	9	0.00758	0.00865	1
	0_1688_02	20	613	5 (0)	0.00292	0.00374	217	0	0.00000	0.00000	387	5	0.00456	0.00585	1
	CL1694	172	305	19 (0)	0.01305	0.01144	59	2	0.00676	0.00156	245	18	0.01461	0.01386	5
	0_12929_02	20	735	0 (0)	0.00000	0.00000	287	0	0.00000	0.00000	444	0	0.00000	0.00000	–
	0_14221_01	24	539	6 (3)	0.00372	0.00281	289	2	0.00264	0.00187	285	4	0.00477	0.00371	1
Average		52		8.00	0.00523	0.00543	191.2	1.3	0.00312	0.00291	325.7	6.8	0.00623	0.00652	8
P. armandii	1_1609_01	80	387	9 (1)	0.00547	0.00519	101	2	0.00467	0.00622	282	7	0.00583	0.00490	0
	PTIFG2009	56	505	12 (4)	0.00611	0.00485	194	3	0.00397	0.00326	311	9	0.00744	0.00584	4
	0_1688_02	96	613	10 (0)	0.00372	0.00661	217	0	0.00000	0.00000	387	10	0.00582	0.01034	2
	CL1694	36	305	6 (0)	0.00572	0.00514	59	1	0.00494	0.00355	245	5	0.00593	0.00555	0
	0_12929_02	52	735	4 (1)	0.00143	0.00212	288	4	0.00364	0.00540	440	0	0.00000	0.00000	0
	0_14221_01	108	539	14 (1)	0.00525	0.00378	293	7	0.00524	0.00380	287	7	0.00536	0.00382	4
Average		71		8.00	0.00462	0.00462	192	2.8	0.00374	0.00371	325.3	6.3	0.00506	0.00508	8

Table 2

Haplotype diversity and neutrality tests for Pinus pumila, P. griffithii, P. koraiensis, and P. armandii: number of haplotypes (Nh), haplotype diversity (Hd), Tajima’s D (D), Fu and Li’s D∗ (D∗), Fu and Li’s F∗ (F∗), and Fay and Wu’s H (H).

		Haplotype diversity			Neutrality tests
Species	locus	Nh	Hd	H	D	D∗	F∗
P. pumila	1_1609_01	11	0.730	2.503	−0.088	1.535∗	1.183
	PTIFG2009	20	0.936	1.801	0.292	1.015	0.916
	0_1688_02	8	0.796	−0.305	0.293	1.192	1.069
	CL1694	6	0.735	0.825	−0.604	−0.576	−0.679
	0_12929_02	12	0.809	1.584	−1.216	1.059	0.346
	0_14221_01	16	0.925	2.624	−0.095	1.288	0.988
Average		12.167	0.822	1.505	−0.236	0.919	0.637
P. griffithii	1_1609_01	4	0.550	0.508	−0.456	0.959	0.642
	PTIFG2009	10	0.778	−0.345	0.694	0.253	0.449
	0_1688_02	4	0.551	−0.074	−1.071	1.407∗	0.772
	CL1694	1	0.000	−	−	−	−
	0_12929_02	1	0.000	−	−	−	−
	0_14221_01	5	0.641	−0.202	0.654	0.950	1.001
Average		4.167	0.420	−0.019	−0.030	0.432	0.477
P. koraiensis	1_1609_01	6	0.818	1.636	1.151	1.351	1.474
	PTIFG2009	7	0.689	−7.409	0.660	0.479	0.623
	0_1688_02	6	0.844	0.533	1.133	1.300	1.411
	CL1694	17	0.820	1.378	−0.363	1.730∗∗	1.141
	0_12929_02	1	0.000	−	−	−	−
	0_14221_01	6	0.758	−4.030	−0.934	−0.504	−0.696
Average		7.167	0.655	−1.315	0.274	0.726	0.659
P. armandii	1_1609_01	9	0.762	1.554	−0.149	0.733	0.532
	PTIFG2009	13	0.902	−0.138	−0.681	−0.368	−0.543
	0_1688_02	11	0.861	0.051	2.239∗	1.399	1.972∗∗
	CL1694	6	0.837	0.078	−0.332	1.259	0.943
	0_12929_02	3	0.520	0.019	1.289	0.089	0.501
	0_14221_01	18	0.802	0.122	−0.842	1.051	0.482
Average		10	0.781	0.281	0.254	0.694	0.648

Table 3

Maximum-likelihood estimates and 90% highest posterior density (HPD) intervals for demographic parameters obtained from pairwise IM multilocus analyses.

Comparison	q1	q2	qa	m1	m2	t	T (Mya)
P. pumila and P. griffithii	0.8326	0.0918	2.8650	0.0005	2.0450	0.9050	0.356
HPD90Lo	0.5481	0.0407	0.7525	0.0005	0.7850	0.4850	0.191
HPD90Hi	1.1916	0.1898	12.9304	0.4605	5.1450	9.1150	3.582
P. pumila and P. koraiensis	0.8142	0.3264	4.8591	0.2275	1.8865	1.9190	0.754
HPD90Lo	0.5571	0.1800	1.3648	0.0025	0.0035	1.0450	0.411
HPD90Hi	1.1505	0.5637	13.1795	1.6375	3.8465	33.9150	13.327
P. pumila and P. armandii	0.3738	0.4087	0.5456	0.4550	3.1950	2.6200	1.030
HPD90Lo	0.1840	0.2111	0.0032	0.0050	0.5650	1.3400	0.527
HPD90Hi	0.6528	0.7457	5.8013	3.9750	6.9750	39.9800	15.711
P. griffithii and P. koraiensis	0.1470	0.5369	1.5049	1.0750	0.2950	0.8125	0.319
HPD90Lo	0.0624	0.2989	0.6903	0.0350	0.0050	0.3175	0.125
HPD90Hi	0.3057	0.9019	5.2869	3.5450	1.7350	4.5675	1.795
P. griffithii and P. armandii	0.1846	0.7019	1.4730	1.1525	0.3050	5.6875	2.235
HPD90Lo	0.0726	0.4046	0.0346	0.0575	0.0050	1.7875	0.702
HPD90Hi	0.4087	1.1445	12.4550	3.5175	1.8650	64.8375	25.479
P. koraiensis and P. armandii	0.4364	0.6272	0.2650	0.4305	0.0100	5.7900	2.275
HPD90Lo	0.2209	0.3327	0.0059	0.0035	0.0100	1.8700	0.735
HPD90Hi	0.8251	1.0748	10.5714	5.0505	4.4100	18.7300	7.360