Introduction

Accumulating evidence suggests that the Pleistocene climate oscillations, which started around 1.8 million years ago (Ma), have seriously affected the geographic distribution of species and thereby their intraspecific genetic diversity (Hewitt, 2000; Barnosky et al., 2005). However, it is still unclear whether these climatic oscillations also triggered speciation events. For example, fossil evidence indicates a widespread morphological stasis throughout the Pleistocene (Morris et al., 1995; Willis and Niklas, 2004). However, Pleistocene glaciations should have fragmented the continuous distributions of many species into disjunct populations, which might have promoted the production of numerous new species through allopatric speciation (Ribera and Vogler, 2003; Weir and Schluter, 2004). In addition, the new habitats created by climate changes may also have led to parapatric divergence in a few genera through adaptive evolution, especially in some climate-sensitive regions, for example, in desert regions. There are indeed a few well-examined examples of Pleistocene plant speciation (for example, Martín-Bravo et al., 2010; Levsen et al., 2012).

Advances in population genetics inference methods based on multiple loci offer the opportunity to estimate the timing of interspecific divergence and to detect signs of adaptive evolution accompanying speciation (Hey and Nielsen, 2004). The estimated timescales based on coalescent analyses of multilocus DNA sequence data may not be as precise as estimates based on fossils, but, are good enough for the generation of a temporal hierarchy (Wakeley, 2008). In the present study, we analyzed DNA sequence variation to test the influence of Pleistocene climate oscillations on the speciation of two desert poplars from central Asia. The species studied here are Populus euphratica Oliv. and P. pruinosa Schrenk. Both belong to sect. Turanga (a small poplar section with only three species) and occur commonly in central Asia as the dominant species of sparse desert forests (Eckenwalder, 1996). The natural range of P. euphratica extends from west China to south Morocco, whereas P. pruinosa is restricted to west China and adjacent regions. The third species of sect. Turanga, P. ilicifolia, occurs naturally and narrowly in Kenya and is geographically isolated and phylogenetically distant from these two species (Browicz 1966). P. euphratica and P. pruinosa are morphologically well differentiated. For example, P. euphratica has two types of glabrous or almost glabrous leaves. Young individuals and twigs from the basal part of old trees bear lanceolate leaves, whereas leaves from old twigs of mature trees are oval, extravagant oval, broad-triangular or dentate broad-ovate. In contrast, all leaves of P. pruinosa are ovate or kidney-shaped with thick hairs. In addition, the two species have different habitats: P. euphratica is distributed in dry deserts with deep underground water and P. pruinosa occurs in deserts with underground water close to the surface, near ancient or extant rivers. Despite these clear differences, these two species seem to have diverged recently because they share numerous ancestral polymorphisms at chloroplast DNA and Simple Sequence Repeat loci (Wang et al., 2011). Here we used an Isolation-with-Migration (IM) model (Hey and Nielsen, 2004) to estimate the divergence time between P. euphratica and P. pruinosa using 11 unlinked nuclear loci. We also used Approximate Bayesian Computation (De Mita and Siol, 2012) to test for departure from the standard coalescent model in each species.

Materials and methods

Sampled materials

Silica-dried leaves were collected from each of six P. euphratica and six P. pruinosa populations representative of their geographic range in west China (Figure 1; Supplementary Table S1). In each population, we collected five individuals in widely and evenly distributed sample sites and all individuals sampled were at least 500 m apart.

Figure 1
figure 1

Locations of sampling sites in the Northwest of China. White and black circles represent the P. euphratica and P. pruinosa, respectively. (a) Structure analysis based on 42 single nucleotide polymorphisms (SNPs). (b) The number of clusters (K)=2 was showed. Bold lines distinguish populations, and blue and purple sections within each population represent the probability of assigning the individuals to each cluster. The numbers below the bars indicate the population numbers.

Full size image

Loci sampling

A total of 11 primer sets (NHX1, NHX6, NDPK, GA20-oxl, GO, PhyB2, pC, HRGP, CPI, GOl and TR) were designed based on expressed sequence tag sequences of P. euphratica (Brosche et al., 2005). These loci are located on nine different chromosomes and those from the same chromosome are far apart on the Populus trichocarpa genome (Supplementary Tables S2, S3). In addition, we had already sequenced two chloroplast DNA fragments (trnV and trnL) and the nuclear internal transcribed spacer (ITS) region (Wang et al., 2011).

Polymerase chain reaction amplification and sequencing

PCR amplifications were performed in a 25 μl volume containing 10–50 ng of diploid DNA, 50 mM Tris-HCI, 1.5 mM MgCl2, 0.5 mM dNTPs, 2 μM of each primer and 1.0 units of Taq polymerase (TAKARA, Dalian, China), using a Gene-Amp PCR system 9700 DNA Thermal Cycler (Applied Biosystems, Foster City, CA, USA). The thermal profile consisted of the primary denaturing at 94 °C for 5 min, followed by 36 cycles, each of which involved 55 s at 94 °C, 1 min at the primer-specific annealing temperature and 55 s at 72 °C, and then a final extension phase of 7 min at 72 °C. Amplification products were then purified using a TIAN quick MidiPurification Kit (TIANGEN, Beijing, China) following the manufacturer’s protocol. Sequencing reactions were performed with the PCR primers to cover the whole PCR segment using an ABI Prism BigDye Terminator Cycle Version 3.1 Sequencing Kit and an ABI3130 xlor3730 xl Genetic Analyzer (Applied Biosystems). The heterozygous PCR DNA fragments were cloned into pGEM T-easy vector and randomly selected clones were sequenced (TAKARA). Clone sequences (5–10 per allele) were compared with sequences obtained through direct sequencing and with other clone sequences for the same individual to help identify polymerase errors and PCR-mediated recombination. All sequences of the same locus were aligned with Mega 4.1 (Kumar et al., 2008). All new sequences obtained in this article have been deposited in the GenBank Databases under accession numbers KC962665-KC962856.

Nucleotide diversity

The average number of nucleotide differences per site between two sequences in a sample (π) and Watterson’s estimate of the population mutation rate (θw), the number of segregating sites (S), the number of haplotypes (Nh) and the haplotype diversity (He) were calculated using DnaSP 5.00.04 (Librado and Rozas, 2009).

Linkage disequilibrium and recombination

The minimum number of historical recombination events (Rm) was estimated using the four-gamete test of Hudson and Kaplan’s (1985). We further investigated the decay of linkage disequilibrium (LD; McVean et al., 2002) in terms of the Zns statistic (Kelly, 1997) over physical pairwise distance following Remington et al., (2001). Pairwise LD was estimated using the average of squared allele-frequency correlation (r2) between pairs of polymorphic sites. The nonlinear regression is simulated using the R statistical package (http://www.r-project.org/).

Departure from the standard neutral model

The site frequency spectrum of each locus was summarized with Tajima’s D (Tajima, 1989), Fu and Li’s D* and F* (Fu and Li, 1993), Fay and Wu’s H (Fay and Wu, 2000) as well as the E and DH test (Zeng et al., 2006). Departure from the standard neutral model was tested by comparing the observed values of the summary statistics with their expected distributions. The latter were obtained through 10 000 coalescent simulations. We also used the multilocus Hudson-Kreitman-Aguadé test (HKA; Hudson et al., 1987) implemented in Jody Hey’s program HKA (http://lifesci.rutgers.edu/~heylab) to assess the fit of the data to the neutral equilibrium model. Besides, we repeated the HKA tests for P. euphratica and P. pruinosa separately using one sequence of P. trichocarpa as outgroup.

The likelihood that natural selection has occurred at individual loci was estimated using the recently developed Maximum Frequency of Derived Mutations (MFDM) test (Li, 2011), which is based on the fact that selection leads to an unbalanced gene genealogy topology. As migration may also cause unbalanced trees, we used individual sequence from the other species as a migration detector Li, 2011 to analyze this possibility. It is worth pointing out that the MFDM test is insensitive to the presence of population bottleneck or population expansion if a single Wright-Fisher population is assumed.

Population structure

In order to quantify the extent of population genetic structure, we estimated Wright’s fixation index FST among populations using an analysis of molecular variance as implemented in Arlequin 3.0 program (Excoffier et al., 2005). The significance of FST was tested by comparing the observed value with the distribution of FST after 10 000 permutations of sequences among populations. Furthermore, we used the program STRUCTURE version 2.3.2 (Hubisz et al., 2009) to assess the correspondence between geographical grouping and genotypic clustering. An admixture model was used that includes linkage among polymorphisms for the ancestry model in which polymorphic sites in each locus were treated as a single linkage group and the independence of allele frequencies among populations was assumed. The number of clusters, K (from 1 to 12, corresponding to the number of populations), was assessed and allowance was made for the correlation of allele frequencies between clusters. Twenty runs were performed with a burn-in of 100 000 and then 500 000 iterations. The most likely number of clusters was estimated using the original method from Pritchard et al. (2000) and the ΔK statistics given in Evanno et al. (2005).

Within species demography

To compare the fit of the two species to different demographic models we used Approximate Bayesian Computation implemented in the egglib software package (De Mita and Siol, 2012). Three models were evaluated: (i) a standard neutral model (SNMR); (ii) a population expansion model (PEMR); and (iii) a population bottleneck model (BNMR). The first model is determined by two parameters, the population mutation rate (θ) and the population recombination rate (ρ). The PEMR model has one additional parameter, the exponential growth rate (α). The BNMR is a more complex model that is determined by θ, ρ, T (time until end of bottleneck), d (the length of the bottleneck) and f (the size of the population during the bottleneck). The three models and associated parameters are described in more details in St Onge et al. (2011). We used the following wide uniform priors: θ: 0–0.05, ρ: 0–0.5, α: 0–50, T: 0–10, d: 0–10, f: 0–1 and the number of segregating sites, Tajima’s D (Tajima, 1989) and the standardized version of Fay and Wu’s H (Fay and Wu, 2002; Zeng et al., 2006) as summary statistics. For each of the models we simulated 1 million values by randomly sampling from the prior and retained 1% of these values for linear regression of parameter values. The three models were compared using a simple rejection algorithm essentially as described in Fagundes et al. (2007).

Haplotype structure

We constructed relationships between all recovered haplotypes at each nuclear locus, chloroplast DNA fragments and ITS region by the Median-Joining network. Its implementation followed the default parsimony connection limit of 95%, with indels considered as single mutation events and coded as a fifth state.

Between-species divergence

Three approaches were used to obtain the estimate of divergence between the two species.

First, we used the log-likelihood statistics of nested models, which can be assessed using a χ2 test to detect whether the data are consistent with a strict allopatric speciation model (Hey and Nielsen, 2007). Once runs had converged, we tested for the fit of the data to simpler demographic models by the ‘Load-Trees’ mode of IMa2. We created three replicate data sets that each contained a randomly chosen subset of five loci and of two populations of different distances. Then, parameters including current and ancestral effective population sizes, divergence time, long-term rates of effective migration and the ancestral population splitting were estimated using the program IMa2 based on the IM model. In cases when recombination within loci was detected, we used the IMgc program to extract maximally informative blocks of non-recombining sequences for individual markers (Woerner et al., 2007).

Based on the HKY mutation model, three independent (different random number seeds) runs of at least 2 to 3 million steps following a 100 000 steps burn-in were used to ensure convergence. All runs involved 10 independent chains, and the lowest effective sample size among the six parameters was at least 200 in each case, as recommended in the IMa2 documentation. Long, well-mixing runs were repeated at least twice with different random seed. If these independent runs generated similar posterior distributions, we assumed that the chains had converged to a stationary distribution. Upper bounds of the prior distributions for each parameter were set based on the results of a preliminary run. Maxima of posterior distributions for each parameter were well under the upper bounds of the prior distributions for all three runs. Maximum-likelihood estimates and 95% highest posterior density ranges for each parameter were converted to meaningful biological quantities using a mutation rate estimate of 2.5 × 10−9 synonymous substitutions/site/year (Tuskan et al., 2006; Invarsson 2008) and the generation time of 10 (8–15) years (Li et al., 2003). Because our data include nonsynonymous and noncoding sites as well as synonymous ones, the ratio of the average total genetic divergence to the average synonymous genetic divergence (Ktotal/KS) was used to obtain the mutation rate (μ=Ktotal/KS × 2.5 × 10−9) (Ikeda et al., 2009). The mean per-locus mutation rate was estimated to be 2.6 × 10−6 substitutions/locus/year, which was then used to estimate divergence time between the two species in years.

Second, we also estimated the same parameters using MIMAR (MCMC estimation of the Isolation-Migration model Allowing for Recombination; Becquet and Przeworski, 2009). Compared with IM, MIMAR takes recombination into account and is based on Hey and Wakeley (1997) summary statistics of the joint frequency spectrum. For each locus, S1 and S2 are the number of polymorphic sites unique to P. euphratica and P. pruinosa, respectively, Ss is the number of sites with shared alleles between the two species, and Sf is the number of fixed alleles in either species. The prior range for θ1 is Uniform (0.001, 0.01), the prior range for T is uniform (0, 100 000) generations, the expected number of migrants, M, will range between (0.135 and 7.389). To test the validity of the MIMAR runs, a goodness-of-fit test was performed using the package MIMARGOF, with simulations being carried out using the posterior distributions of the parameters estimated in MIMAR.

Third, divergence time was estimated using ITS variation between the two species, which was calculated as the average DNA sequence distance divided by twice the sequence mutation rate (μ). The woody perennial substitution rates ranged from 0.38 × 10−9 to 7.83 × 10−9 substitutions/site/year (mean=2.15 × 10−9 subs/site/year; Kay et al., 2006).

Results

Nucleotide diversity and neutrality tests

Two haplotypes per individual were obtained for each locus, representing 60 alleles per species. The aligned sizes of the eleven loci ranged from 621 to 1863 bp, with an average size of 868 bp (9546 bp in total, Supplementary Table S2). The variation was on the average higher in P. pruinosa (πt=0.00407; θwt=0.00380) than in P. euphratica (πt=0.00327; θwt=0.00297; Table 1). Our average estimates of silent site sequence diversity are 0.00398 and 0.00506 for P. euphratica and P. pruinosa, with the highest estimates within species being 0.01442 and 0.01529, respectively.

Table 1 Nucleotide variation and haplotype diversity at 11 loci in 12 populations of P. euphratica and P. pruinosa
Full size table

Most tests of the standard neutral model at individual loci were nonsignificant, as summarized by the average values of the test statistics, except for PhyB2 in P. euphratica that showed significant values for all tests (Table 2). We further conducted a multilocus HKA test. We did not detect any significant departure from the standard neutral model with the HKA test when considering P. euphratica and P. pruinosa (χ2=15.6774; P=0.736). We compared the test statistics with a distribution generated from 10 000 coalescent simulations and did not detect significant values for both pairs of contrasts (P. euphratica/P. trichocarpa, χ2=2.297, P>0.05; P. pruinosa/P. trichocarpa, χ2=0.663, P>0.05).

Table 2 Neutrality tests at 11 loci in 12 populations of P. euphratica and P. pruinosa. Significance was estimated with coalescent simulations under the standard neutral model
Full size table

The MFDM test detected evidence of selection at one locus: GO, in P. euphratica (P=0.0339; at a significance level of P<0.05). In addition, a low P-value was also obtained for PhyB2 in P. euphratica (P=0.0678; 0.05<P<0.07). For each locus in P. euphratica, we arbitrarily picked one individual from P. pruinosa as the migration detector and one individual from P. trichocarpa as outgroup to build a UPGMA (Unweighted Pair Group Method with Arithmatic Mean) tree (Supplementary Figure S1). The analyses indicated that migration is not responsible for these unbalanced trees.

Recombination and interspecific LD

Four-gamete tests indicated that the minimum number of recombination events (Rm) in all loci for each species ranged from 0 to 3 in P. euphratica and from 0 to 5 in P. pruinosa, with the estimates being generally higher in P. pruinosa than in P. euphratica at all nuclear loci (Table 1). Supplementary Figure S2 illustrated the decline of LD with physical distance, using pooled data of all 11 loci. LD decays more rapidly to lower levels (r2<0.2) in P. pruinosa than in P. euphratica.

Within species demographic history

In both species, the average Tajima’s D values were positive (P. euphratica=0.72 and P. pruinosa=0.29) indicating an excess of intermediate frequency variants compared with neutral expectations, even though single loci only rarely deviated significantly. In order to identify a demographic scenario compatible with the collected sequence data three simple demographic models, a standard neutral model (SNMR), a population expansion model (PEMR) and a bottleneck scenario (BNMR), were compared. The precision in parameter estimates was in general not very high and under the BNMR, θ was the only parameter showing a narrow, unimodal distribution (Supplementary Figure S3). The model comparisons lend support for either the SNMR or the BNMR and clearly show (as expected with positive Tajima’s D values) that the data are not compatible with the simple growth model (Table 3).

Table 3 The proportion of accepted values for each of the three tested demographic models at four different tolerance levels
Full size table

Population structure and differentiation

A hierarchical analysis of molecular variance approach was used to additionally quantify differentiation between the species. We observed 36.8% of the total variation among species and only 18.7% among populations within species, whereas 44.5% of the total variation was found within populations (Table 4). The average FST estimate is 0.212 for P. pruinosa, which was lower than the value for P. euphratica (FST=0.367). There were substantial variation across the eleven loci in levels of genetic differentiation, with the lowest estimate of FST=0.322 for GO and the highest estimate FST=0.868 for NDPK (Table 4). Structure indicated that the most likely number of clusters for the entire data set was K=2 (Figure 1; Supplementary Figure S4) with a clear delineation of the two species, except for one population (Pop10) of P. pruinosa showing a strong signal of admixture with P. euphratica.

Table 4 Summary of population differentiation (FST) and hierarchical analysis of molecular variance (AMOVA) at 11 loci in P. euphratica and P. pruinosa
Full size table

Divergences between the two species

Genealogies of the recovered haplotypes observed at each locus were constructed (Figure 2). Most of the sampled loci did not show complete lineage sorting except for GA20-oxl and NDPK. Similarly, networks for the chloroplast DNA showed no obvious divergence between the two species. In contrast, the ITS network was composed of two clusters with species-specific lineages.

Figure 2
figure 2

Gene genealogies of the eleven nuclear loci, nrITS and two region aligned chloroplast DNA (cpDNA). Colors in the pie chart indicate the haplotype origin. The size of the pie is proportional to the haplotype frequency found in the two poplars. Branch lengths longer than one mutation step are marked on each branch.

Full size image

The strict allopatric speciation model was rejected for the complete (all loci, all individuals), and random five loci by the ‘Load-Trees’ mode (P<0.001 for all tests; Supplementary Table S4). We repeated runs of simulations with the IMa2 program and revealed unambiguous marginal posterior probability distributions of the demographic parameters for the two species. Because the loci, GO and PhyB2, are putatively under selection, we removed them from the demographic analysis. The effective sample size values for the time parameter are above 200 in the different analyses. Results are shown in Table 5 and Figure 3. Estimates of effective population sizes for each species are 97 907 for P. euphratica and 93 365 for P. pruinosa. These sizes are about two times smaller than the effective size of the ancestral population that gave rise to these species, indicating that both P. euphratica and P. pruinosa effective population sizes have decreased since the split of the two species. The marginal posterior probability distribution of the divergence parameter, t, showed a sharp peak at 2.149, which was converted into a divergence time of 0.83 Ma (0.66–1.23 Ma) with 95% highest posterior density ranged from 0.24 to 5.17Ma (Table 5). MIMAR estimate of the divergence of the two species was of the same order of magnitude and around 1.09 Ma (0.93–1.28 Ma). Estimate of divergence time based on ITS genetic distance between the two species, and the woody perennial substitution rates (2.15 × 10−9 substitutions/site/year; Kay et al., 2006), was also of the same order of magnitude with a value of 1.37 Ma. Estimates of the long-term effective migration rates between P. euphratica and P. pruinosa are asymmetric. Maximum-likelihood estimates of 2Nm from P. pruinosa into P. euphratica and from P. euphratica into P. pruinosa are 0.009 and 0.293, respectively.

Table 5 MLE and the 95% HPD intervals of demographic parameters from IMa2 analyses
Full size table
Figure 3
figure 3

Marginal distribution of the posterior probability of six demographic parameters estimated by the IM model. (a) Effective population sizes for P. euphratica and P. pruinosa. (b) Effective population sizes for their common ancestor at the time of divergence. (c) Effective migration rates between the two species. (d) Divergence time.

Full size image

Discussion

We conducted population genetic analyses of two desert poplar species, P. euphratica and P. pruinosa, based on sequence data from 11 nuclear loci. These loci are evenly distributed across the genome (Supplementary Table S2), and should therefore provide a reasonable estimation of species-level divergence. Coalescent analyses indicated that these two poplar species diverged within the Pleistocene, at a stage when the largest glaciation occurred in Asia. These findings suggest that the increased aridity owing to the Pleistocene climate change might have triggered the recent speciation of these two poplar species in response to continuous differentiation of desert habitats. We also found that these two species diverged in the presence of gene flow and that two loci, GO and PhyB2, showed evidence of selection in P. euphratica.

Gene flow and selection signature

Our IM analysis rejected a strict allopatric speciation hypothesis (Supplementary Table S4) but favored the hypothesis that the two species may have diverged in the presence of asymmetrical gene flow (Table 5). The asymmetrical gene flow from P. euphratica to P. pruinosa is consistent with our previous study based on nuclear Simple Sequence Repeat loci (Wang et al., 2011), even though it is worth pointing out that estimates of effective population sizes were vastly different between the two studies. The low effective population sizes reported in our previous study might simply reflect the inadequacy of the model assumed by Migrate (Beerli, 2009). In contrast to IM models that enforce a split between the two descendant species at a certain time, Migrate assumes an equilibrium island model. If migration is limited, the time to observe the last coalescent event may therefore be very long and smaller effective population size would therefore be required to obtain the same level of polymorphism in the two descendant species. In addition, our samples did not cover the total natural ranges of the two species and this sampling limitation may bias our estimates of effective population sizes, genetic structure and gene flow. However, the asymmetrical gene flow between the two species also agrees well with the habitat distribution of two species: P. euphratica occurs in the drier desert associated with deep underground water levels, whereas P. pruinosa occurs in damper environments within general a more limited underground water level. Therefore, P. euphratica can survive in the habitats of P. pruinosa, but the reverse seems to be more difficult, which may have accelerated the interspecific gene flow from P. euphratica to P. pruinosa.

In this study, we used various neutrality tests and the MFDM test to examine the selection signatures of each locus. The advantage of the latter is to avoid effects of demographic processes by focusing on the topology of the gene genealogy, which theoretically should not be affected by demographic events such as population expansion or bottlenecks (Li, 2011). In addition, this method considers possible gene flow and helps to minimize the number of false positives. The MFDM tests indicated that the gene, GO, was under selection in P. euphratica and also suggested the presence of selection for PhyB2 in the same species although the signal was weaker. On the other hand, all values from Tajima’s D, Fu and Li’s D*, Fu and Li’s F*, DH and E tests suggested positive selection for PhyB2 (Table 2). The glycolate oxidase edited by GO has an important role in the oxidation of glycolate to glyoxylate in plants, whereas PhyB2 is associated with phenology (Ingvarsson et al., 2006). Hence, the observed deviations at these genes do not come as a surprise, and reinforce their status as major candidate genes for further studies on adaptation in poplar species.

Interspecific divergence

Both IM analyses of multiple nuclear loci and ITS mutations indicate that the two species diverged during the early Pleistocene. This conclusion depends on the mutation rate and generation times adopted here. We used the same mutation rate of 2.5 × 10−9 per site per year for multiple nuclear loci as in previous poplar studies and the generation time of 10 (8–15) years. For ITS, we also used previous mutation calibrated for woody perennials (Kay et al., 2006). The mutation rates adopted here are not direct estimates and could therefore be under- or overestimates of the true mutation rate. Similarly, estimating the average generation time remains difficult in most tree species with long lifespan (Petit and Hampe, 2006). We chose a generation time of 10 years but as the trees have a longer lifespan a higher value would also be reasonable. In a recent study of P. euphratica, reproductive characteristics in 20 stands from Xinjiang the mean age of the trees was 15.4 years with an average s.d. of 2.6 (Cao et al., 2012). Assuming that the stands were a random sample that would imply that our estimates of generation time may not be strongly off the mark. In addition, we only sampled the co-occurring regions of two species, whereas P. euphratica extends further westward to Europe and southwestward to northern Africa, where P. pruinosa is absent. Allopatric populations may increase the overall interspecific divergence. All these uncertainties may cause some substantial error, but we believe that the estimated timescales would likely remain within the Pleistocene or the late Pliocene. These estimated stages are in line with the largest glaciation and coldest climate between 0.5 and 1.2 million years again recorded in Asian mountains. For example, the ice sheet at this stage in the Qinghai-Tibetan Plateau was five to seven times larger than it was at later stages (Owen et al., 2005). The decreased temperature caused the enlargements of the deserts in central Asia (Zheng et al., 2002). This increase in desertification and diversity of desert habitats during the Pleistocene may have acted as an effective stimulus to promote the divergence of the two species.

The estimated divergence between P. euphratica and P. pruinosa is older than that between Populus balsamifera and P. trichocarpa in North America (Levsen et al., 2012). However, all these estimates fall within the Pleistocene, providing further evidence for Pleistocene tree speciation in Populus. In contrast to the macroevolutionary stasis observed in trees since the Pliocene, these two examples, as well as those previously reported (Martín-Bravo et al., 2010) highlight the importance of the Pleistocene climate oscillations for plant species diversification.

Data archiving

Data deposited in the Dryad repository: doi:10.5061/dryad.k70h6. Sequence data have been submitted to GenBank: accession numbers KC962665-KC962856, including the following sequence data:

CPI_Pe: KC962665—KC962669

CPI_Pp: KC962670—KC962678

GO_Pe: KC962679—KC962693

GO_Pp: KC962694—KC962709

NHX1_Pe: KC962710—KC962717

NHX1_Pp: KC962718—KC962728

NHX6_Pe: KC962729—KC962730

NHX6_Pp: KC962731—KC962732

PhyB2_Pe: KC962733—KC962745

PhyB2_Pp: KC962746—KC962756

TR_Pe: KC962757—KC962762

TR_Pp: KC962763—KC962770

GA-20oxl_Pp: KC962771—KC962772

GA-20oxl_Pe: KC962773 - KC962777

NDPK_Pp: KC962778—KC962787

NDPK_Pe: KC962788—KC962798

pC_Pp: KC962799—KC962813

pC_Pe: KC962814—KC962823

GOl_Pe: KC962824—KC962829

GOl_Pp: KC962830—KC962834

HRGP_Pe: KC962835—KC962849

HRGP_Pp: KC962850—KC962856