Late Pleistocene climate change promoted divergence between Picea asperata and P. crassifolia on the Qinghai–Tibet Plateau through recent bottlenecks

Abstract Divergence during the early stage of speciation can be driven by a population bottleneck via reduced gene flow and enhanced lineage sorting. In this study, we aimed to examine whether such bottlenecks occurred during the initial speciation of two closely related spruce species Picea asperata and P. crassifolia occurring on the Qinghai–Tibet Plateau (QTP). We analyzed sequences of three chloroplast, two mitochondrial DNA fragments and a further 13 nuclear loci from 216 individuals of the two species. Both species showed a low level of genetic diversity in contrast to other congeners occurring in the QTP and adjacent regions. The estimated population sizes of P. asperata and P. crassifolia are less than the ancestral population size before splitting. These results together with multiple statistical tests (Tajima's D, Fu and Li's D* and F*) suggest that these two species underwent recent bottlenecks. Based on approximate Bayesian computation (ABC), we also determined that the period of the population shrinkage was consistent with the interspecific divergence during the late Pleistocene. The reduced population sizes and the divergent selection may together have triggered the initial divergence under high gene flow between these two species. Our results therefore highlight the importance of climatic oscillations during the late Pleistocene in promoting speciation through changing demographic sizes of the ancestral species on the QTP and in adjacent regions.


Introduction
Divergence can be promoted by population bottlenecks because of reduced quantity of gene exchange and that small populations with increased rates of genetic drift can fix more derived mutations than large populations (Templeton 1980(Templeton , 2008Barton and Charlesworth 1984). Recently, many studies have focused on the earlier stage of speciation Liu et al. 2014), and some of them have shown that divergence among species can occur even in the presence of gene flow and that species boundaries are maintained by divergent selection (Hewitt 1996;Petit and Excoffier 2009;Schluter 2009;Nosil 2012;Sousa and Hey 2013). However, whether population bottlenecks have played a role in promoting divergence between closely related taxa has rarely been studied (Butlin et al. 2012).
When the strength of random genetic drift increases following reduced population size during a bottleneck period or during the colonization of a new environment, most of the genetic variation within a species is likely to be lost (Young et al. 1996;Excoffier et al. 2009). In addition, when the population size decreases below a critical level, effective mating with conspecifics tends to be difficult to maintain (May and McLean 2007). Populations during a bottleneck period are always at a high risk of extinction due to the lack of adapted variations or mating partners (May and McLean 2007;Excoffier et al. 2009). The lack of genetic variations makes it difficult, experimentally, to model and examine historical bottlenecks of target species (Hein et al. 2004). Therefore, it is challenging to find empirical evidence for the hypotheses that historical bottlenecks promoted or accompanied the initial speciation stage. However, a high level of genetic variation is frequently found in woody species pollinated over long distances by the wind, for example, spruces (Heuertz et al. 2006;Li et al. 2015;Sun et al. 2015). Thus, enough but fluctuated levels of genetic variation in these species may allow us to detect the signals of bottlenecks that occurred during the evolutionary divergence of closely related species. In this study, we aimed to examine whether bottlenecks promoted divergence between Picea asperata Masters and P. crassifolia Komarov.
Picea asperata and P. crassifolia have wide distributions across the eastern Qinghai-Tibet Plateau (QTP), but P. crassifolia extends to a more northerly latitude. Phenotypic differences between P. asperata and P. crassifolia are obvious. The leaf apices of P. asperata are always acute; in contrast, the leaf apex of P. crassifolia is obtuse and the length/width ratio of the leaf is always lower than that of P. asperata (Fig. 1). The winter buds of P. asperata are resinous, unlike those of P. crassifolia (Fu et al. 1999). Previous work based on chloroplast and mitochondrial DNA variations showed that these two species were closely related to each other and none of the derived mutations were species specific (Du et al. 2009), indicating that they may have diverged very recently. We used sequence data from 13 nuclear loci and three chloroplast and two mitochondrial DNA fragments to investigate the evolutionary histories of P. asperata and P. crassifolia, in combination with planting experiments in common garden. We addressed the following questions. Did historical population bottlenecks and gene flow occur when they diverged? If population shrinkages did occur, when did these bottlenecks take place? Were they correlated with the late Pleistocene climate oscillations?

Materials and Methods
Sampling, planting, and DNA sequencing We sampled needles from 11 and 15 populations, respectively, through the main natural distribution range of Picea asperata Masters and P. crassifolia Komarov between 2007 and 2012 ( Fig. S1; Table 1). The latitude, longitude, and altitude of each sampling location were recorded using an eTrex GIS monitor (Germany); detailed information is presented in Table 1. About 3~11 individuals were sampled from each population. The distance between sample trees was at least 100 m. Needles were collected in the field and dried immediately with silica gel. We compared the width of needles with length of~1.7 cm, and five needles for each of ten trees per species were measured. The significance level of differences was calculated using R version 3.1.1 (two-sample Student's t-test, https://www.r-project.org/). Seeds collected from 20 trees (10 trees per species from different populations) were planted in a common garden at Lanzhou University (Xiaguanying), in order to observe whether their morphological traits were stable under common environmental conditions. A total of 13 nuclear loci (4CL, EBS, FT3, GI, M002, M007D1, PCH, Sb16, Sb29, Sb62, se1364, se1390, and xy1420) were sequenced using an ABI 3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA) following the methods in Li et al. (2010). Sequence data were basecalled using PHRED with CodonCode Aligner software (CodonCode Corporation) and checked using MEGA version 5.0 (Kumar et al. 2008). Sequences were aligned by CLUSTAL X (Thompson et al. 1997) with MEGA5. Sequences with heterozygous sites were rephased and separated into two allelic sequences by PHASE (Stephens et al. 2001;Stephens and Donnelly 2003) in the software package DnaSP V5 (Librado and Rozas 2009) with default parameters. All sequences have been deposited in Gen-Bank (Accession Numbers: KX212991-KX213062). Chloroplast and mitochondrial DNA variation data (trnL-trnF, trnS-trnG, ndhK/C, nad1 intron b/c, and nad5 intron1) were obtained from previous studies (Meng et al. 2007;Du et al. 2009).

Analyses of population genetic diversity and structure
Population genetic parameters, including the number of segregating sites (S), Watterson's parameter h w (Watterson 1975), and nucleotide diversity p (Tajima 1983), were computed using DnaSP V5.10.01 (Librado and Rozas 2009), after excluding insertions/deletions (Table S2). For each locus, Tajima's D statistic (Tajima 1989), Fu and Li's D* and F* (Fu and Li 1993), and the number of segregating sites were also calculated using DnaSP V5.10.01 (Librado and Rozas 2009). Sequences of P. breweriana were downloaded from GenBank and were used as an out-group when performing D* and F* tests.
To understand interspecific genetic differentiation at neutral nuclear loci, we computed Φ ST using Arlequin V3 (Excoffier et al. 2005) and the significance level was based on 1000 permutations. To check whether data from these two species could be delimitated into two groups, NET-WORK version 4.6.1.3 was used to determine genealogical relationships among the mtDNA, cpDNA, and nuclear haplotypes at each locus (Bandelt et al. 1999). Indels at each locus were coded as characters using simple indel coding as implemented in SeqState (Muller 2006). In addition, we used the Bayesian algorithm implemented in STRUCTURE ver. 2.3.3 (Pritchard et al. 2000) to infer the possible genetic structure with the admixture model. All 26 populations and 13 loci were used in STRUCTURE analysis. Sites that showed significant statistical association after Bonferroni correction were removed. Linkage disequilibrium among SNPs within genes was limited as the number of SNPs per gene was small. For each K-value (1 ≤ K ≤ 10), we performed 20 independent runs with a burn-in of 200,000 and 500,000 iterations, respectively. The program Distruct version 1.1 (Rosenberg 2004) was used to generate graphical representations of the data. The most likely number of clusters was estimated using the original method from Pritchard et al. (2000) and with the DK statistics (Evanno et al. 2005).

Analyses of the isolation and migration (IM) model and bottleneck model
We estimated current and ancient population size, divergence time, and interspecific migration rate under the isolation and migration model using an approximate Bayesian computation (ABC) approach implemented in ABCtoolbox (Wegmann et al. 2010). This divergence model under gene flow assumed that an ancestral lineage with population size N diverged into two lineages T generations ago. The population sizes of P. crassifolia and P. asperata are referred to as Nc and Na. The migration rate from P. crassifolia to P. asperata is referred to as Mcra->Masp and that in the reverse direction as Masp->Mcra. We used 13 statistics to summarize the pattern of molecular variation at all nuclear loci, including the number of haplotypes within each species, the number of polymorphic sites within each species, the number of private polymorphic sites for each species, the number of pairwise differences, Tajima's D (Tajima 1989), Fu's F S (Fu 1997), pairwise F ST , the number of pairwise differences between species, and the total number of polymorphic sites. These statistics were computed by Arlequin version 3.0 for both observed and simulated datasets (Excoffier et al. 2005). Simulated datasets were generated using a standard algorithm in the ABCtoolbox software package, and fastsimcoal was used to perform the simulation (Excoffier and Foll 2011). The gene information is listed in Table S1. For the simulations, the mutation rate was assumed to be 1 9 10 -8 per site per generation (Li et al. 2010;Sun et al. 2014). The prior of each parameter was assumed to follow the uniform distribution. A total of 500,000 datasets were simulated by fastsimcoal, and 5000 of them were retained. We applied the regression adjustment general linear model to compute marginal density and to generate posterior distributions of all parameters (Wegmann et al. 2010). The low genetic diversity in these two species and positive D, D* and F* values at many loci (Table 2) along with the estimated population parameters indicated that a population bottleneck could have occurred during the initial divergence of these two species. So, we dated the population bottlenecks for P. crassifolia and P. asperata using ABCtoolbox. In addition, due to the low differentiation between these two species, we also analyzed the pooled sequence data from the two species to infer the likely historical population bottleneck. The bottleneck model assumed that an ancestral population with population size N underwent population contraction Tb generations ago, then began to expand Te generations ago. During the bottleneck period (from Tb to Te, forward in time), the population size was Nb. Four statistics, including the segregating sites, Tajima's D, nucleotide diversity, and Fu's Fs, were used to summarize the genetic information in the sequence dataset. We used the standard algorithm in ABCsampler (Wegmann et al. 2010) and the program fastsimcoal (Excoffier and Foll 2011) to simulate samples using the gene information given in Table S1. The mutation rate was set to 1 9 10 À8 substitutions per site per generation (Li et al. 2010;Sun et al. 2014), and the prior of each parameter was assumed to follow the uniform distribution. We excluded simulated samples with exceptional values because their inclusion would lead to failure when processing them. A total of 1,500,000 samples (500,000 samples for each species and also for pooled data) were simulated. We retained the 5000 simulated samples which were closest to the observed values for each model to calculate marginal density and to generate posterior distributions for Tb and Te with the regression adjustment general linear model.

Results
Stable morphological difference between two species in a common garden Needles of the trees that were five or more years old were collected from the upper side of branchlets (five needles per tree) and compared. The differences in leaf apex between species were obvious, and the needles of P. crassifolia were always wider than those of P. asperata, when comparing needles of the same length (means, 1.44 mm < 1.69 mm, Student's t-test, P < 0.001; Fig. 1). These observations suggested that the morphological differences between these two species are stable and are likely to be caused by genetic changes during their divergent histories.

Nucleotide diversity and population structure
Genetic diversity was detected at 12 loci (Tables S1, S2). No variation was detected at nuclear locus 4CL (Table S2). A total of 71 and 78 polymorphic sites were detected in P. asperata and P. crassifolia, respectively (Table S1). The silent nucleotide diversities within P. asperata and P. crassifolia were 0.00334 and 0.00393, respectively. Similarly,  the total nucleotide diversity in P. asperata (0.00264) was slightly lower than in P. crassifolia (0.00304). The fixed index Φ ST between P. asperata and P. crassifolia (0.053) was significantly higher than in the random scenario (P < 0.001), indicating clear population differentiation between the species. However, genealogies constructed based on the variations in cpDNA, mtDNA, and at each nuclear locus revealed that P. asperata and P. crassifolia shared most of their genetic variations (Figs. 2, 3). Like the NETWORK analyses which revealed high levels of sharing allele diversity (Figs. 2, 3), the result of STRUCTURE showed that the most likely number of clusters (K) is 2 (Fig. S2). When K had a value of 2-4, we found that these two species could not be grouped into distinct clusters (Fig. 4).

Demographic inference and the parameterized IM model
To examine the demographic dynamics of P. asperata and P. crassifolia, we performed three tests: Tajima's D, Fu and Li's D* and F* tests. Positive average values of D, D* and F* were found (Table 2), indicating P. asperata and P. crassifolia may have experienced recent population bottlenecks.
Next, using ABCtoolbox (Wegmann et al. 2010), we estimated the time of the bottleneck based on data from each species and the pooled data. The bottlenecks that these two species experienced occurred~196,000-37,000 years (Fig. 6, 95% HPDI, a scale of 50 years per generation) ago, and the divergence between P. asperata and P. crassifolia accompanied with the bottleneck. Thus, the bottleneck likely began earlier than the divergence between P. asperata and P. crassifolia, and the reduced population size should have restricted gene flow between P. asperata and P. crassifolia when splitting occurred.

Discussion
In this study, we integrated two approachestransplant experiments and population genetic analysesto examine the phenotypic difference and whether the population demography promoted the initial divergence between two conifer species, Picea asperata and P. crassifolia. The population size of the ancestral lineage before splitting is larger than the current population sizes of P. asperata and P. crassifolia (Fig. 5). The interlineage divergence likely occurred following the onset of population shrinkage ( Fig. 6), indicating that serious bottlenecks during the late Pleistocene period might have promoted the divergence between P. asperata and P. crassifolia. The distinct phenotypic differences between P. asperata and P. crassifolia and extensive genetic sharing indicate that some undetected divergent loci may have contributed to these species' morphological divergences. Our results highlight the importance of Pleistocene climate change in promoting species divergence through serious bottlenecks.
Population shrinkages in response to climate warming Silent and total nucleotide diversities for P. asperata and P. crassifolia are lower than other related spruce species found in the east of the QTP and adjacent regions (Li et al. 2010). For example, the total nucleotide diversities of P. likiangensis (0.0093) and P. wilsonii (0.0087) are much higher than either P. crassifolia (0.00304) or  . Posterior distribution of parameters within the isolation and migration model based on the sequences from nuclear loci. Nasp, Ncra, and N indicate the effective population size of P. asperata, P. crassifolia, and their ancestral lineage, respectively. Masp->cra denotes the migration rate from P. asperata to P. crassifolia, and the reverse migration rate is denoted by Mcra->asp. T denotes the divergence time between the two species (logarithmically transformed). Table 3. Posterior mode estimates and 95% highest posterior density (HPD) intervals for demographic parameters in the IM model based on the nuclear multilocus nuclear sequence data. Nasp and Ncra, current population size of Picea asperata and P. crassifolia, respectively; N, effective population size of the common ancestor; T, divergence time between P. asperata and P. crassifolia; Mcra->asp, the effective migration rate from P. crassifolia to P. asperata; Masp ->cra, the effective migration rate from P. asperata to P. crassifolia.

Parameter
Nasp  . Thus, the lower levels of nucleotide diversity in both P. asperata and P. crassifolia suggest that these two species may have experienced more severe reductions in their population sizes. This suggestion is also supported by the positive values of Tajima's D, Fu and Li's D* and F* at 12 nuclear loci (Table 2). Furthermore, the estimated current population sizes of P. asperata and P. crassifolia are much smaller than their ancestral population size ( Fig. 5 and Table 3). Therefore, we conclude that the population bottlenecks might occur during the evolution of P. asperata and P. crassifolia. During the late Pleistocene period, the average temperature in the eastern QTP increased following the retreat of the largest glaciation, which occurred around 800,000 years ago (Shi et al. 1998). However, the low temperature continued until the late Ionian stage, around 300,000-126,000 years ago (Shi et al. 1998). In response to the climate warming, the ranges of cold-adapted species, such as spruce and yew, contracted, in contrast to the range expansions found for other plants during the inter-or postglacial periods Sun et al. 2015). In this study, according to the estimate based on the sequences from all nuclear loci, the bottlenecks of P. asperata and P. crassifolia began 196,000 years ago (3918 generations, Fig. 6). This estimate is similar to the onset of population size reductions (around 200,000 years ago) in the evolutionary histories of P. likiangensis and P. wilsonii . Therefore, the climate rewarming during the late Pleistocene period appears to have reduced the ranges and population sizes of many cold-adapted species distributed in the eastern QTP.
Divergence between Picea crassifolia and P. asperata For spruce species occurring on the QTP and in adjacent regions, most divergences between related species have been found to predate the largest glacial event (Li et al. 2010Zou et al. 2013). A recent divergence (around 127,000 years ago) between P. asperata and P. crassifolia was revealed here based on the analysis of the IM model (Table 3), indicating that these two species are at the initial stage of recent speciation. The lack of fixed variation (Table S2) also supported this conclusion of recent divergence. Simulations have shown that at such a stage of speciation, there is insufficient time to accumulate genetic differentiation between the newly diverged lineages (Nielsen and Wakeley 2001). Incomplete lineage sorting results in extensive genetic sharing between lineages, as revealed by our STRUCTURE and NETWORK analyses (Figs. 2, 4). An alternative explanation for these shared variations is the gene flow between the two species (Table 3). Primary/second contacts of the distributional ranges of P. asperata and P. crassifolia after the initial divergence could also have provided opportunities for gene exchange (Fig. S1). Picea asperata and P. crassifolia bear stable phenotypic differences despite of high level of genetic sharing (Fig. 1), suggesting that the species boundary might be mainly maintained by the divergent selection although the underlying genetic mutations are unknown.
Environmental differences can promote adaptive divergence through selection. Similarly, the reduced gene flow caused by limited dispersal ability and decreased population size can also promote interspecific divergence by accelerating the fixing of different adaptive alleles in the diverging lineages (R€ as€ anen and Hendry 2008). In this study, we found that when P. asperata and P. crassifolia diverged, their populations were much smaller than their ancestral lineage (Figs. 5, 6). Therefore, we tentatively conclude that the reduced population sizes restricted interspecific gene flow and further promoted divergence between P. asperata and P. crassifolia, although some interspecific gene flow continued (Table 3).

Supporting Information
Additional Supporting Information may be found online in the supporting information tab for this article: Figure S1. Locations of each of the sampled populations for the two spruce species studied here, P. asperata and P. crassifolia. Figure S2. Estimated number of clusters (K) obtained with Structure. Table S1. Nucleotide variation at 13 loci in P. asperata and P. crassifolia. Table S2. The number of segregating sites at 13 loci for P. asperata and P. crassifolia.