Importance of incomplete lineage sorting and introgression in the origin of shared genetic variation between two closely related pines with overlapping distributions.

Genetic variation shared between closely related species may be due to retention of ancestral polymorphisms because of incomplete lineage sorting (ILS) and/or introgression following secondary contact. It is challenging to distinguish ILS and introgression because they generate similar patterns of shared genetic diversity, but this is nonetheless essential for inferring accurately the history of species with overlapping distributions. To address this issue, we sequenced 33 independent intron loci across the genome of two closely related pine species (Pinus massoniana Lamb. and Pinus hwangshanensis Hisa) from Southeast China. Population structure analyses revealed that the species showed slightly more admixture in parapatric populations than in allopatric populations. Levels of interspecific differentiation were lower in parapatry than in allopatry. Approximate Bayesian computation suggested that the most likely speciation scenario explaining this pattern was a long period of isolation followed by a secondary contact. Ecological niche modeling suggested that a gradual range expansion of P. hwangshanensis during the Pleistocene climatic oscillations could have been the cause of the overlap. Our study therefore suggests that secondary introgression, rather than ILS, explains most of the shared nuclear genomic variation between these two species and demonstrates the complementarity of population genetics and ecological niche modeling in understanding gene flow history. Finally, we discuss the importance of contrasting results from markers with different dynamics of migration, namely nuclear, chloroplast and mitochondrial DNA.

Introduction

Understanding speciation has been a long-standing interest since the publication of ‘On the Origin of Species' (Darwin, 1859). The allopatric speciation model was considered as the null hypothesis for decades (Mayr, 1963; Coyne and Orr, 2004). However, recent investigations have suggested that speciation in the face of gene flow has been more common (Nosil, 2008), especially in plants where gene exchange among closely related species is estimated to occur in at least 25% of species (Mallet, 2005). Scenarios of speciation with gene flow are not all equivalent though. On one hand, gene flow can be continuous from an early phase of speciation to completion of reproductive isolation (sympatric or parapatric speciation), in which case divergent selection is thought to play the major role in promoting reproductive isolation between incipient species (Coyne and Orr, 2004; Feder ). On the other hand, current sympatric or parapatric species distributions may result from secondary contacts occurring during late phases of speciation. This may occur when an allopatric period was not long enough for species to evolve complete reproductive isolation, for example, when they come into secondary contact during climatic cycle oscillations (Barton and Hewitt, 1985; Melo-Ferreira ; Rieseberg ). Introgression following such secondary contact may introduce foreign alleles and reduce the differentiation between species gene pools at loci unlinked to speciation genes (Coyne and Orr, 2004; Petit and Excoffier, 2009).

Patterns of genetic diversity within and among populations carry much information about population demographic history and therefore have been widely used to infer speciation history (Hey and Nielsen, 2004; Won and Hey, 2005; Duvaux ; Roux ; Butlin ). However, retention of ancestral polymorphism because of incomplete lineage sorting (ILS) and secondary gene flow produce very similar patterns of shared genetic diversity, making the investigation of speciation history difficult (Rieseberg ; Charlesworth ; Sousa and Hey, 2013; Huerta-Sánchez ). Indeed, ILS can be the cause of shared genetic diversity a long time after species divergence. Under a simple allopatric speciation model, drift alone requires 9–12 Ne (effective population size) generations to make incipient species reciprocally monophyletic at more than 95% of loci (Hudson and Coyne, 2002). Therefore, species with long generation times, such as coniferous trees, or recently isolated populations, often share genetic variation (see, for example, Du ; Wang ; Li ; Sun ).

The recent availability of new tools in population genetics now allows us to discriminate better between gene flow and ILS. The development of coalescent-based frameworks, for example, isolation with migration (IM) model-based programs (Wakeley and Hey, 1998; Nielsen and Wakeley, 2001; Hey, 2010) and approximate Bayesian computation (ABC) frameworks using coalescent genealogy samplers (Beaumont ; Kuhner, 2009; Csilléry ; Sunnåker ), has allowed direct comparison of different scenarios of population divergence. For example, Qu found, using the program IMa (Hey, 2010), that alleles shared between eastern and western lineages of the parrotbill Paradoxornis webbianus might have resulted from both ILS and secondary gene flow. Ecological niche modeling based on present distribution data can also be applied to infer demographic dynamics and historical variation of species ranges (Elith ). Used jointly with coalescent-based frameworks, this method allows refined interpretations about possible secondary contacts between closely related species during past climate oscillations (see, for example, Levsen ; Sun ).

When geographic distribution information is available, ILS and secondary gene flow can also be distinguished by comparing patterns of genetic diversity between pairs of neighboring and distantly located populations of the different species (Muir and Schlötterer, 2005). Gene flow is expected to happen preferentially among neighboring populations that therefore results in higher levels of intraspecific genetic diversity and lower levels of interspecific genetic differentiation than between distantly located populations (Petit and Excoffier, 2009). In contrast, under the ILS scenario, shared polymorphism is expected to be distributed evenly in all populations. Using this approach, Muir and Schlötterer (2005) showed that FCT between two closely related oak species was randomly distributed across distant and nearby populations, suggesting a major role for ILS. However, Lexer showed that FCT was unevenly distributed across the genome and that a limited number of loci showed significantly increased FCT, a pattern they argued was more likely because of recurrent gene flow and selection at differentiated loci.

Closely related coniferous tree species often share genetic variation (Ma ; Li ; Wachowiak ; Ren ). One possible reason lies in their long generation times (>20 years in average) and large Ne. This means that most coniferous species have probably experienced climatic oscillations as rapid successions, in terms of drift units (Tgen/Ne, where Tgen is speciation time in generations), of isolation periods and expansion intervals with secondary contacts. Using this rationale, most coniferous species diverged relatively recently (∼0.1–0.2 million generations, that is, 2–3 Ne generations ago; Leslie ; Mao ), leading them to maintain a high proportion of ancestral diversity. Moreover, interspecific gene exchange promoted by frequent secondary contacts during climate cycle turnover may have contributed to maintaining weak reproductive barriers between closely related species (see, for example, Ma ; Li ).

Pinus massoniana and P. hwangshanensis are two closely related pines occurring in Southeast China with overlapping distributions. They have a different preference for altitude as populations are usually found at low (below 900 m) and high altitudes (above 700 m) for P. massoniana and P. hwangshanensis, respectively (Fu ). The species can occur on the same mountains but at different (partially overlapping) altitudes (parapatric populations) or occur on different mountains (allopatric populations) (Zhou , 2014; Li , 2012). Analysis of cytoplasmic genetic variation revealed that mitochondrial diversity was extensively shared and randomly distributed across populations of the two species (Zhou ), whereas chloroplast variation was clearly species specific. In coniferous trees, chloroplast, mitochondrial and nuclear genomes are paternally, maternally and biparentally inherited respectively (Wagner ; Neale and Sederoff, 1988; Mogensen, 1996). These two pines therefore offer a great opportunity to test for the preponderance of ILS and introgression during speciation and to compare results from markers with differential inheritance.

Here, we analyzed noncoding genetic variation at 33 intron loci across the genome for both parapatric and allopatric populations of the two species. First, we compared patterns of genetic diversity and population admixture between parapatric and allopatric populations of both species. We then conducted an ABC analysis to infer the demography and speciation scenario of both species. Along with the demographic modeling, we used ecological niche modeling to track distribution changes of the two species during the Pleistocene climate changes. Our results suggest that secondary contact(s) rather than ILS explains shared nuclear genetic variation between the two pine species, in contrast to the earlier conclusions based on cytoplasmic markers alone.

Materials and methods

Sampling, sequencing and analyses of genetic diversity

We used a range-wide sample of individuals from eight parapatric as well as seven and three allopatric populations of P. massoniana and P. hwangshanensis, respectively (Figure 1a). In order to reduce the probability of sampling hybrids in parapatric populations, we sampled trees at altitudes at least 100 m lower or higher than the margins of the contact zones, that is, at 600 and 1000 m, respectively. Basic passport information of the sampled populations can be found in Supplementary Table S1, and geographic information is presented in Figure 1a. Two individuals of Pinus koraiensis belonging to subgenus Strobus (Wang ; Gernandt ) from Northeast China were used as outgroup.

Figure 1

(a) Ranges and sampled populations for P. massoniana and P. hwangshanensis. (b) Speciation models tested in the ABC analysis.

We amplified and Sanger sequenced a set of 37 loci distributed across the genome. Genomic DNA was extracted from haploid megagametophytes of germinated seeds using the QIAGEN DNeasy Plant Mini Kit (QIAGEN, Inc., Valencia, CA, USA). The primers were designed based on previous published sequences of closely related pines (Zhou ). PCR and sequencing information have been described in detail in Zhou . We only used intron sequences in order to reduce the confounding effect of natural selection. Despite this, four outlier loci showed Tajima's D strongly deviating from our sample mean, suggesting a close genetic proximity with selected sites. As a precaution, we therefore excluded these loci from analyses. The remaining 33 loci were aligned using MUSCLE (Edgar, 2004) as implemented in Mega 5 (Tamura ).

In order to compare genetic diversity and differentiation in parapatric and allopatric populations of both species, we computed the following statistics for each locus: nucleotide diversity (πs, Nei, 1987; θs, Watterson, 1975), the number of haplotypes (K), Tajima's D (Tajima, 1989), average number of nucleotide substitutions (Dxy, Nei and Kumar, 2000), average number of net nucleotide substitutions (Da, Nei and Kumar, 2000), differentiation among populations within species (FST) and differentiation between the two species (FCT). All the summary statistics (Supplementary Table S2) were computed using DnaSP v 5.10 (Librado and Rozas, 2009).

Bayesian admixture analysis

To examine genetic structure in each species and to compare genetic admixture between parapatric and allopatric groups, we used the Bayesian clustering approach using the software STRUCTURE v.2.3.3 (Hubisz ). We ran 6 replicates for each value of K from 1 to 8 (100 000 burn-in cycles followed by 1 000 000 cycles of data collection, admixture allowed) and identified the best K for our data set according to the statistic ΔK (Evanno ). Clustering analyses were first conducted for each species to examine intraspecific population structure. As within-species population structure was weak (Figure 2), we merged populations into allopatric and parapatric groups for each species to compare levels of admixture at different biogeographic settings.

Figure 2

Population structure and admixture in parapatric and allopatric populations of the two species (STRUCTURE results).

Approximate Bayesian computation

ABC allows inference on complex models by using data simulated under a model to replace its intractable likelihood function (Csilléry ).

ABC data set

In order to reproduce the structure of the observed data set closely during simulations, we used msnsam, a modified version of the program ms, allowing different sample sizes across loci (Hudson, 2002; Ross-Ibarra ). Each simulated data set consisted of 33 loci with 48–60 and 40–44 individuals for P. massoniana and P. hwangshanensis, respectively. As ms uses an infinite site model and only simulates ancestral-derived states, we removed all positions with indels, missing data and triallelic single-nucleotide polymorphisms before computing statistics.

Model description

We inferred the best speciation scenario among four models (Figure 1b), all derived from the ‘classical' IM model (Nielsen and Wakeley, 2001). In these models, an ancestral population of size Na split at time Ts in the past to form two populations of size N1 (P. massoniana) and N2 (P. hwangshanensis). These models are mainly characterized by the following features: in model I (for Isolation, Figure 1b), there is no gene flow between species; in model F (for gene Flow, Figure 1b), there is continuous gene flow (potentially at asymmetrical rates m12 and m21); in model FI, gene flow ceases at time Tm; and in model IF, gene flow starts at time Tm (Figure 1b).

Priors on demographic parameters

All demographic parameters had uniform prior distributions and some were scaled by a factor N=100 000 for convenience. Prior ranges were defined by trial and error to fully cover the posterior distribution ranges of each parameter for each model, thereby allowing us to properly estimate model posterior probabilities and parameter posterior distributions. Prior distributions of population parameters N/N0, N2/N0 and N/N0 ranged over 0.001–6. Prior distributions of migration parameters 4N0m12 and 4N0m21 ranged over 0.0001–8 and those of time parameters Tm/Ts and Ts/4N0 over 0.0001–1 and 0.0001–20, respectively.

Priors on locus-specific parameters

In order to account for heterogeneity in local mutation and recombination rates (μ and r respectively), we computed per generation θ (4Nμl, where l is the locus length) and ρ (4Nrl) for each locus i using outgroup sequence information. We first estimated μ as μ=D/2T, where D is the mean of the Ds (net substitution rates per site, equation 12.67 in Nei and Kumar, 2000) between P. massoniana or P. hwangshanensis and the outgroup and T is the divergence time with the outgroup in generations (that is, 2.25–4.25 M generations assuming a generation time of 20 years; Willyard ). Second, we estimated r/μ by calculating the average of γ/π across species where γ is an estimator of ρ (Hey and Wakeley, 1997). Because the divergence time T is a range estimate rather than a point estimate, we obtained two estimates of μ and r corresponding to bounding values of T. In order to account for this uncertainty, we used log-normal distributions as hyperpriors for each θ and ρ. For any θ or ρ distribution, the log mean and log s.d. were set equal to the mean and variance computed from the corresponding two estimates of μ or r. Although simple, this approach is conservative to estimate θ and ρ uncertainties. For loci with no outgroup sequences available, we used the mean and s.d. of μ and r across loci instead. Computation of population genetics estimators were made using SITES (Hey and Wakeley, 1997; https://bio.cst.temple.edu/~hey/software/software.htm#SITES).

Summary statistics for ABC analysis

We used the means and s.d. values across loci of the following statistics: per population π and Tajima's D, the number of sites with fixed derived alleles in one population yet polymorphic in the other (sxf and sxf, where i and j refer to P. massoniana and P. hwangshanensis, respectively), the number of sites with a derived allele fixed in one population and the ancestral allele fixed in the other (sf and sf) and the number of sites with a shared derived allele (ss), FST, DA and DXY (Nei and Kumar, 2000, equations 12.66 and 12.67). Although this set of 24 summary statistics is probably not optimal, it is only slightly different from the one Duvaux have shown to be informative enough to distinguish between speciation models. We therefore chose not to implement recent methods optimizing the choice of summary statistics (Aeschbacher ; Fearnhead and Prangle, 2012) in order to keep the same pipeline during the whole study. Summary statistics were computed using a modified version of the program used in Roux (http://www.abcgwh.sitew.ch).

Model selections and parameter estimations

We used tools from the ‘abc' R package (Csilléry ; cran.r-project.org/web/packages/abc/vignettes/abcvignette.pdf). We first determined the best model with the function ‘postpr'. We set the number of units in the hidden layer to 9, and loaded 2.5 million simulations per model. We retained the best 5000 simulations during the rejection step and fitted a neural network regression on them (100 neural networks). The false positive rate for model selection was estimated by using 500 pseudo-observed data sets from the second best model and by performing model choice using each pseudo-observed data set in turn in place of the real data set. The false discovery rate is thus the proportion of times we observe a posterior probability equal or superior to the real posterior probability of the best model. To do so, we used a modified version of the function ‘cv4postpr' that uses as pseudo-observed data sets simulations retained by the rejection step only (that is, not drawn from the whole parameter space as in the original function). We used the function ‘abc' with the best 5000 simulations to estimate the posterior distributions of the model parameters. Statistics were log transformed for the regression step and simulations with parameter values outside the prior ranges were removed before estimation of the posterior distributions (Duvaux ).

Goodness of fit (GoF) of the best model

To assess the GoF of the best model, we used posterior predictive simulations (Gelman, 2003). We ran 10 000 simulations for the 33 loci using parameter values sampled from the joint posterior distribution and computed the resulting summary statistics. We used two sets of statistics for the GoF: those used to perform the ABC analysis and, in order to avoid overestimation of the fit, a set of new statistics. The second set was made of 10 statistics including mean and variance across loci for S the expected number of segregating sites across both species, an estimator of θw for each species ( and ) and the number of sites polymorphic in one population yet fixed for the ancestral allele in the other (sx and sx). For both data sets, an initial check of the fit was obtained by performing principal component analysis on a subset of 2000 posterior predictive simulations along with 2000 prior (original) simulations and the observed data set. Finally, using posterior predictive distributions, we computed a two-sided P-value for each statistic (that is, two times the probability of obtaining the observed statistic or a more extreme value under the estimated model).

The full ms command lines for model simulations and other R code developed (R Core Team, 2013) for the ABC analysis are available in Dryad.

Species distribution modeling

Along with ABC inferences of species demography, we modeled the dynamic changes in species distributions during the Pleistocene climatic oscillations (for example, the Last Glacial Maximum (LGM; ~21 000 years before present) and the Last Inter-Glacial (LIG; ~120 000–140 000 years before present)). The knowledge of current species geographical distribution along with series of historical climate data allow ecological niche modeling that can also be used to track down the distribution changes of closely related species that can potentially be used to indicate the possibility of gene flow during speciation (Elith ). We used geographical information system-based methods with layers of historical and current environmental variables under the maximum entropy model implemented in the program MaxEnt 3.3.3k (Phillips ; Phillips and Dudík, 2008). GPS data from 123 and 109 localities (our sampling sites and museum records, Supplementary Table S1) were used to generate species distribution models for P. massoniana and P. hwangshanensis, respectively. We used 19 environmental variables (Supplementary Table S4) from the WorldClim database (Hijmans ) with spatial resolution of 2.5′ (arc-minutes) and 30′ as environment layers. The methods assume that the probability of occurrence of a species has a uniform probability distribution considering present distributional data (that is, maximum entropy; Phillips ; Phillips and Dudík, 2008). Assuming niche conservation over time (Wiens and Graham, 2005), we applied this model to LGM and LIG climatic layers to predict historical species distributions. Both the CCSM (community climate system model) (Collins ) and MIROC (model for interdisciplinary research on climate) (Hasumi and Emori, 2004) were used to predict species distributions during the LGM. We generated 10 jackknife replicates with deletion of 20% of species occurrence localities for each analysis to account for error in the predictive modeling. Replicate was run for 500 interactions with a convergence threshold 10−5. We configured the machine-learning algorithm to use 75% of species records for training and 25% for testing the model. The program DIVA-GIS v7.5 (Hijmans ) was used to create the graphics of species distributional ranges based on raster cells from original species distribution modeling with a logistic probability of presence >0.1.

To measure the levels of climate niche overlap between P. massoniana and P. hwangshanensis, we calculated the niche overlap statistic Schoener's D (Schoener, 1968) for the two species using the program ENMTools 1.3 (Warren ). We also computed the overlap between species ranges at LGM and present with a threshold 0.5 in the ENMTools.

Results

Comparisons of genetic diversity and differentiation between parapatric and allopatric populations

In total, we aligned 11 588 bp noncoding nuclear sequences. Summary statistics on genetic diversity and divergence at each locus are listed in Supplementary Table S2. As interspecific gene flow should bring foreign genetic diversity and reduce genetic differentiation, we expected to observe higher levels of genetic diversity and lower differentiation in parapatric populations than in allopatric populations. Genetic diversity was found to be similar in allopatric and parapatric populations of both species (θallo=0.0053±0.0044, θpara=0.0056±0.0050 for P. massoniana and θallo=0.0083±0.0064, θpara=0.0090±0.0061 for P. hwangshanensis). Although 24 out of the 33 examined loci had higher FCT values in allopatric than parapatric populations (Supplementary Table S1), the average interspecific differentiation was not significantly different between allopatric and parapatric populations (0.41±0.26 and 0.37±0.26, respectively; Supplementary Table S2).

Genetic structure and admixture

STRUCTURE (Hubisz ) analyses showed that our data were best described by K=2 genetic clusters. The inferred ancestry of each of the 104 pine samples was calculated and reported as the fraction assigned to each of the two clusters (Figure 2). We observed that each species was made up of one specific cluster with some individuals showing admixed patterns. These individuals were slightly more frequent in parapatric groups (2.9% and 4.9% for P. massoniana and P. hwangshanensis, respectively) than in allopatric groups (2.0 and 2.6% for P. massoniana and P. hwangshanensis, respectively; Figure 2).

Demographic histories and speciation models

In order to infer the speciation history of the two species, we performed an ABC analysis (Beaumont ; Csilléry ) on the set of 33 noncoding loci. Among our four speciation scenarios, those that did not incorporate any current gene exchange, namely scenarios I and FI, had posterior probability equal to zero, consistent with the present distribution of the two species (see Table 1). Furthermore, the secondary contact scenario was the most compatible with observed data (posterior probability=0.79). The false discovery rate, computed by including the F and IF models only, was very low (0.015), meaning that the model selection was robust to our choice of statistics (Csilléry ). The GoF further demonstrates the relatively good ability of the IF scenario to generate data close to the observed data, highlighting its overall appropriateness to model the speciation history of our two species. Apart from the second axis of the first principal component analysis (statistics of the ABC analysis, Supplementary Figure S1), the observed data always lie in regions of the principal component analysis showing a reasonable density of prior and posterior predictive simulations (Supplementary Figures S1 and S2). The P-values of most summary statistics, including those not used for the ABC, were usually far above 0.05 (Supplementary Figure S3 and Supplementary Table S2). Where not, the posterior predictive simulations show better P-values than the prior simulations (for example, FST). The only aspect of the data poorly predicted by our model was the variation of genetic diversity among loci that was overestimated in the simulations (see θ_std, π_std and D_std in Supplementary Figure S3 and Supplementary Table S3, and the Discussion section).

Table 1

Results of ABC model selection

Model	I	F	FI	IF
PP	0	0.2066	0.0003	0.7932

Abbreviations: ABC, approximate Bayesian computation; F, there is continuous gene flow; FI, gene flow ceases at time Tm; I, isolation model; IF, gene flow starts at time Tm; PP, posterior probability.

One striking result of the IF scenario was the length of the isolation period between the two species, estimated to have lasted 90% of the divergence time (with Ts=591 605 and Tm=45 984 generations (11.83 and 0.92 Mya, respectively), Table 2 and Supplementary Figure S4 modal values). The effective population size of P. hwangshanensis (N2) was estimated to be ∼3 times larger than that of P. massoniana (N1) (253 000 (167 000–351 044) and 80 000 (48 712–143 850) respectively, Table 2) and this difference was significant (P=0.003 estimated from the posterior distribution of the difference), in line with higher genetic diversity found in P. hwangshanensis. However, this did not seem to have a strong impact on the direction of migration as there was no significant gene flow asymmetry between species (2N1m12=1.05, 2N2m21=1.44, P=0.32, Table 2).

Table 2

Parameters estimates and 90% higher probability density (HPD) intervals under the IF (secondary contact) models, with times in generations (generation of 20 years)

Parameter	Mode	Median	HPD90low	HPD90High
N1	79 737	86 639	48 712	143 850
N2	253 189	253 281	167 302	351 044
Na	352 588	393 784	69 280	801 875
m12	4.31E−06	5.92E−06	1.65E−06	1.53E−05
2N1m12	0.871	1.008	0.367	2.26
m21	1.80E−06	2.95E−06	6.26E−07	9.76E−06
2N2m21	1.015	1.452	0.34	4.339
Tm	45 985	79 159	6630	369 906
Ts	591 605	1 033 626	202 235	2 315 909
(Ts−Tm)/Ts	0.94	0.91	0.76	0.98

N and 2Nm parameters are expressed as number of individuals and m are migration rates. Population 1 and 2 refers to P. massoniana and P. hwangshanensis, respectively. Abbreviation: IF, gene flow starts at time Tm.

Species distribution change

In order to track historical distribution changes of both species, we modeled the potential distribution for the LIG, the LGM and present (Figure 3). The accuracy of species distribution modeling was relatively high for both species (area under the curve=0.972 and s.d.=0.022 for P. massoniana; area under the curve=0.952 and s.d.=0.039 for P. hwangshanensis). The simulated distributions (Figures 3e and f) based on present climate data were mostly congruent with current ranges (Figure 1) of the two species. The climate variables that contributed the most percentage in predictions of species distributions were the precipitation of warmest quarter (56%) and the minimum temperature of coldest month (38%) for P. massoniana and P. hwangshanensis, respectively.

Figure 3

Species distribution modeling for P. massoniana (a, c and e) and P. hwangshanensis (b, d and f) from the last inter-glacial to present. The different colors represent the probabilities of occurrence of the species.

Paleodistribution models suggested that P. massioniana had a nearly constant distributional range but P. hwangshanensis experienced a gradual expansion during the Pleistocene climatic oscillations (Figure 3). We compared the modeled distributions at the LGM and the present using ENMTools with a threshold 0.5. The two species showed moderate reduction of ranges sizes at LGM (0.81 and 0.79 for P. massoniana and P. hwangshanensis respectively). The two closely related species showed overall overlapping distributions since the LGM (21 Kya) (Figures 3c and d). The current range overlap (0.69, threshold 0.5) between the two species was similar with the overlap at the LGM (0.46, threshold 0.5).

Discussion

We investigated two alternative hypotheses to explain the extent of shared genetic diversity between two closely related pine species with overlapping distributions, namely (1) retention of ancestral polymorphisms because of ILS and (2) introgression through secondary contact. By contrasting genetic variation of allopatric and parapatric populations, we observed that patterns of between-species genetic differentiation and admixed ancestries, although not significant, were more compatible with the introgression hypothesis. Comparison of speciation scenarios using ABC also favored the secondary contact hypothesis (Figure 1b). Species distribution modeling supported a gradual expansion of P. hwangshanensis during the Pleistocene climatic oscillations leading to the current distribution overlap with P. massoniana. In contrast, results of a previous study showed that mitotypes were commonly shared and randomly distributed across allopatric and parapatric populations (Zhou ). We discuss how these apparently contradictory results can be synthesized. We argue that the discordant information observed from the chloroplast DNA (cpDNA), mitochondrial DNA (mtDNA) and nuclear DNA (nuDNA) genomes actually reflects the differential dynamics in conifers of their main dispersal vectors, that is, seeds or pollen.

Successes and limitations of the ABC approach

As shown by the GoF analysis, the secondary contact (IF) model simulates data reasonably well, although there is still room for improvement. First, the results of the ecological niche modeling suggest that incorporating an exponential expansion for P. hwangshanensis could have improved the model fit. However, we do not expect such improvements to change our main result, namely the secondary contact scenario being largely favored. Indeed, the absence of real expansions from models is expected to have little effect on gene flow estimations (see, for example, table 6 in Strasburg and Rieseberg, 2010). The poorest predicted aspect of the data is the variance of diversity and the distribution of fixed alleles among loci that was over- and under-estimated by the posterior predictive simulations, respectively (for example, θi_std, πi_std and D_std, sfX in Supplementary Figure S3, where std stands for s.d.). A supplementary analysis based on the IF model (not shown) indicated that the genetic variance among loci within a population/species is positively correlated with the effective population size Ne, meaning that if Ne was overestimated by the ABC analysis, so was the variance in the posterior predictive simulations. A well-known cause of Ne overestimation is the exclusion of population structure from models where some actually does exist in reality (Chikhi ; Strasburg and Rieseberg, 2010). Including structure into models of speciation with gene flow is not trivial, however, as it is difficult to implement and, most importantly, increase the ‘curse of dimensionality' (Sunnåker ). Because our data showed only moderate within-species structure (Supplementary Figure S2), we decided to include none in our models. As another possible explanation, Roux showed that modeling the heterogeneity of migration rates across loci considerably improves the fit of speciation models to the data in the context of ancient divergence. Because the divergence history of our pine species appears ancient, it seems likely that modeling the heterogeneity of gene flow across loci would improve the inference of the parameter posterior distributions (for example, by inferring smaller effective population sizes) and thus the quality of posterior predictive simulation (especially for the sfX statistics). It is important to note though that properly fitting the hyperparameters of the β-distribution used to model heterogeneity of gene flow requires a huge amount of data (for example, Roux used 852 loci for a total alignment length of 270 kb) that are not yet available for these two pine species. We therefore restrain ourselves to investigate basic but essential models of divergence. Finally, although 33 loci may appear to be a very small sample for ABC analyses, our false discovery rate analysis showed that this number of loci was sufficient to distinguish among models of speciation with gene flow.

Both ILS and secondary contact contributed to current shared genetic diversity between species

Shared polymorphisms across closely related species might play very important roles in processes of adaptation and speciation. They can be generated in various ways, including hybridization (Rieseberg, 2009; Abbott ), adaptive introgression (Song ; Jones ; Pardo-Diaz ; Staubach ; Hedrick, 2013; Huerta-Sánchez ), balancing selection favoring specific ancestral alleles (Leffer ) and incomplete lineage sorting (Hobolth ). In this study, using a pair of closely related parapatric pine species, we assessed the importance of ILS and secondary introgression in explaining shared polymorphisms. We estimated the species divergence time for nuclear loci to be 2.62 Ne generations (Table 2) corresponding to an early stage of divergence. As the divergence date is much lower than 9 Ne (Hudson and Coyne, 2002), we expect that lineage sorting should be limited, meaning that the two species still share much of their ancestral polymorphism.

Although ILS cannot be ignored in our case, climate oscillations strongly influence species distributions and thus the possibility of secondary contacts (Hewitt, 1996, 2000; Taberlet and Cheddadi, 2002). If the duration of an allopatric episode is not long enough to build complete reproductive isolation between two incipient species, secondary contacts may promote interspecific gene flow resulting in shared genetic polymorphism (Melo-Ferreira ; Rieseberg ). For these two pine species, we showed that our results were more compatible with a secondary introgression scenario rather than with a strict isolation model. Our STRUCTURE analysis suggested the presence of admixture in parapatric populations and the ABC analysis showed that this may be due to an ancient secondary contact dating back to 920 Kya (Table 2). This value is unlikely to be strongly inflated as the modal estimation of the divergence time (~12 Mya) falls within the range of a conservative estimation based on observed Ks (9.17–17.33 Mya—calibrated using the divergence time between Pinus and Strobus and 20 years per generation; Willyard ). Accordingly, our ecological niche modeling suggested that the two species probably had overlapping distribution during the LIG (~140 KYA; Figure 3). Most importantly, it also showed that the range of P. hwangshanensis was more sensitive to climatic oscillations (Figure 3). Together, the ABC and niche modeling therefore suggest that the species demographic history has been complex and probably involves multiple range contractions and expansions leading to successive periods of gene flow that were initiated before the LIG. The possibility of current and past secondary gene flow is also supported by recent studies showing, along elevation gradients, that the two species still hybridize at a low rate in intermediate altitude contact zones (Luo and Zou, 2001; Zhang ). This low introgression rate is more likely to be associated with postzygotic rather than prezygotic isolation. First, the long period of initial isolation (~550 000 generations) should have permitted fixation, through selection or drift, of many new mutations in both species and therefore generation of Bateson–Dobzhansky–Muller incompatibilities before any secondary contact (Orr and Turelli, 2001). Second, field observations suggest that postzygotic barriers may be more prevalent than prezygotic ones. Indeed, seeds from hybrid zones have been shown to have lower germination rate than other seeds, suggesting hybrids may carry Bateson–Dobzhansky–Muller incompatibilities (Li ). In parallel, prezygotic isolation appears relatively weak as the two species still share a synchronous flowering period in April–May (Fu ).

The differential evolutionary dynamics of mtDNA, cpDNA and nuDNA reveal different aspect of conifer histories

In a previous study, we investigated cytoplasmic genetic differentiation between the two pine species and found that chlorotypes were highly species specific, whereas mitotypes were extensively shared and were randomly distributed among parapatric and allopatric populations. The cause of shared mtDNA variation was thought to be the retention of ancestral polymorphism rather than recurrent interspecific gene flow (Zhou ). In contrast, we show here that the nuDNA diversity presents a different picture, that is, 0.6% of fixed and 38% of shared polymorphism sites between species, respectively (Supplementary Table S2). Such conflicting results from markers transmitted differentially through males and females have been observed in other species (Petit ; Petit and Excoffier, 2009; Toews and Brelsford, 2012; Poznik ). Most of them have been explained by sex-specific introgression patterns resulting from the interaction of two distinct effects (Petit and Excoffier, 2009; Toews and Brelsford, 2012). First, loci associated with different sex-specific dispersal vectors have different rates of gene flow, pollen being more efficient to disperse than seeds (Liepelt ). Second, interspecific introgression is negatively correlated with intraspecific gene flow (Currat ). Indeed, high level of intraspecific gene flow efficiently homogenizes within species gene pools. Coupled with partial reproductive isolation (that is, low interspecific introgression rates), this leads to reduce the probability of fixation of locally introgressed alleles. Inversely, low intraspecific gene flow increases the probability of local fixation of introgressed alleles (Rieseberg ; Petit and Excoffier, 2009; Zhou ). In conifers interestingly, mitochondrial, chloroplast and nuclear genomes are maternally, paternally or biparentally inherited through seeds, pollen or both, respectively (Wagner ; Neale and Sederoff, 1988; Mogensen, 1996). In agreement with the above expectations, cpDNA, which is always efficiently dispersed by pollen, shows weak within-species structure and high between-species differentiation, whereas mtDNA, which is only dispersed locally by seeds, shows high intraspecific structure and more commonly shared polymorphisms between species (Du ; Zhou ; Jiang ). Finally, we showed in this study that nuDNA, which is dispersed both by pollen and seeds, shows intermediate patterns.

The occurrence of gene flow is also strongly consistent with cpDNA and mtDNA having an expected effective population size four times smaller than nuDNA (Ne equal to N and 2N, respectively). Given our ABC results, this implies that the cpDNA and mtDNA lineages should have diverged ∼10.5 Ne generations ago and would have a very high probability to be reciprocally monophyletic currently. Only long-lasting or recurrent period gene flow seems to be able to explain why it is clearly not the case for mtDNA. In addition, it is worth noting that the nucleotide substitution rate varies greatly among seed plants genomes, with the mtDNA and nuDNA having the lowest and highest rates of nucleotide substitution, respectively (Wolfe ). In practice, mitochondrial genetic variation has been used very efficiently to infer population structure and has been found to be shared between closely related coniferous species, whereas chloroplast markers have shown weak population structure but clear species boundaries (Petit ; Du ; Zhou ). By reflecting evolutionary dynamics at different timescales and in different sexes, the three kinds of coniferous DNAs provide complementary views of the evolutionary processes.

Data archiving

The sequences of each locus reported in this study were deposited in GenBank under accession numbers: KJ921127–KJ921496.