Incongruent range dynamics between co‐occurring Asian temperate tree species facilitated by life history traits

Abstract Postglacial expansion to former range limits varies substantially among species of temperate deciduous forests in eastern Asia. Isolation hypotheses (with or without gene flow) have been proposed to explain this variance, but they ignore detailed population dynamics spanning geological time and neglect the role of life history traits. Using population genetics to uncover these dynamics across their Asian range, we infer processes that formed the disjunct distributions of Ginkgo biloba and the co‐occurring Cercidiphyllum japonicum (published data). Phylogenetic, coalescent, and comparative data suggest that Ginkgo population structure is regional, dichotomous (to west–east refugia), and formed ˜51 kya, resulting from random genetic drift during the last glaciation. This split is far younger than the north–south population structure of Cercidiphyllum (~1.89 Mya). Significant (recent) unidirectional gene flow has not homogenized the two Ginkgo refugia, despite 2Nm > 1. Prior to this split, gene flow was potentially higher, resulting in conflicting support for a priori hypotheses that view isolation as an explanation for the variation in postglacial range limits. Isolation hypotheses (with or without gene flow) are thus not necessarily mutually exclusive due to temporal variation of gene flow and genetic drift. In comparison with Cercidiphyllum, the restricted range of Ginkgo has been facilitated by uncompetitive life history traits associated with seed ecology, highlighting the importance of both demography and lifetime reproductive success when interpreting range shifts.


Introduction
The warm-temperate deciduous forests of eastern Asia harbor a rich plant species diversity including many "living fossils" belonging to the Tertiary flora~65.5-2.6 million years ago (Ma) (Axelrod et al. 1998;Manchester et al. 2009). The greatest concentration of these relicts (e.g., Ginkgo, Metasequoia, Cercidiphyllum, Davidia) presently extends from southwest China eastward across the Yangtze Valley to the east coast and as far as southern Japan (Lu 1999). These floral elements are considered to be relatively evolutionary stable because although affected by climate change, they were not directly affected by the ice sheets from continental glaciation. The region is thus associated with a reduced magnitude of Quaternary environmental change (≤2.6 Ma) compared with other temperate regions in the Northern Hemisphere (Liu 1988). The most recent paleobiome reconstructions, however, suggest rather more complex range dynamics for this forest biome during the Last Glacial Maximum (LGM; 18,000-24,000 ya) as habitat loss (and gain) in different areas of temperate and subtropical China is thought to have occurred (Harrison et al. 2001;Ni et al. 2014). Thus, the genetic potential to detect population structure, population expansion/contraction may be rather higher than the absence of continental glaciation would suggest (reviewed in Qiu et al. 2011).
Consistent with the former view of limited Quaternary environmental change, molecular evidence from Ginkgo biloba suggests that natural populations have persisted in southwestern and eastern China, without exhibiting any strong genetic signatures of range expansion/contraction (Gong et al. 2008). Similar results in understory herbs also implicitly favor a "stable-range" hypothesis, especially in the southern Yangtze Valley (Qiu et al. 2009a,b,c). However, a case study in Cercidiphyllum japonicum, another tree species inhabiting warm-temperate deciduous Asian forests, suggests massive habitat losses in north-central China but, conversely, increases in southwestern/eastern China during the late-Pleistocene .
Two existing competing hypotheses (and a third presented here) may explain the variation in the degree of postglacial expansion to former range limits between temperate deciduous forest plants in eastern Asia. Qian and Ricklefs (2000) hypothesized that during glacial periods, isolated populations admixed/merged at lower elevations (isolation with admixture). Harrison et al. (2001), on the other hand, argued that they remained isolated during glacial as well as interglacial periods (continual isolation). The essential difference between the two hypotheses is one of admixture (negating isolation) during glacial periods. However, as several authors have noted, isolation is relative and the temporal dynamics of population gene flow and random genetic drift during the Pleistocene make the two hypotheses far from mutually exclusive. Thus, the time frame during which population genetic structure may have formed via drift and the pre-/postdivergence dynamics of gene flow are critical in testing both hypotheses as is the use of multiple phylogeographical study systems to form robust independent interpretations of the two hypotheses. A third, albeit less considered, hypothesis is that life history traits can impose restrictions on the range dynamics and population genetic structure of temperate plants by their influence on lifetime reproductive success (Hamrick and Godt 1996).
The range dynamics of Gingko biloba and Cercidiphyllum japonicum, two wind-pollinated dioecious trees, offer a means of testing the same set of hypotheses of genetic isolation and the possible roles of life history traits with data gathered from codistributed species. Both trees are flagship species of warm-temperate forests in eastern Asia, co-occurring in habitats with rocky outcrops in proximity to streams subject to frequent flooding disturbance (Tang et al. 2012), and their co-occurrence has been conserved since the late Cretaceous (c. 66 Ma; Royer et al. 2003). The two species occupy similar climatic niches although the historical distribution of C. japonicum was possibly more northerly, as indicated by introgression with the cool-temperate sister species C. magnificum ). Both G. biloba and C. japonicum possess genetically diverged lineages which may have resulted from vegetation shifts and habitat displacement by the invasion of boreal forests during glacial periods. However, only C. japonicum shows the genetic signal of range expansion . Moreover, the small, flattened, winged seeds of Cercidiphyllum are dispersed by wind (Fu and Endress 2001), while the large fleshy seeds of Ginkgo are mostly dispersed by animals (Del Tredici 2007); a difference in seed dispersal which may lead to different outcomes during recolonization. The comparison of range dynamics between Ginkgo and Cercidiphyllum thus offers an ideal system to test the aforementioned set of hypotheses regarding gene flow and genetic isolation.
Using phylogenetic, coalescent, and comparative methods, we infer the processes responsible for the formation of disjunct distributions in these two species (new Ginkgo data and published Cercidiphyllum data; Qi et al. 2012), simultaneously assessing the degree of congruence among other members of the same ecosystem. Across its complete Asian range, we report the population genetic structure, the split time for divergent populations, and the rates of population migration for Ginkgo. From the split time, we aim to test whether Ginkgo lineages diverged in a congruent time frame as its floral associate, Cercidiphyllum (1.89 Ma). We date the Ginkgo split to a narrow window during glacial periods and compare the degree of genetic isolation experienced by the two species during these episodes. The effects of gene flow and several life history traits on the observed genetic structure are discussed, as are the degree of postglacial expansion to former range limits in eastern Asia and the variance in expansion observed among temperate tree taxa.

Tree samples
Young leaves from a total of 761 individuals of G. biloba were collected from 25 populations in China ( Fig. 1A; Table S1 in Appendix S1) and preserved in silica gel. This total includes 126 individuals from 13 populations analyzed previously (Gong et al. 2008). To exclude recent planting by humans, we sampled trees only with a diameter at breast height (DBH) of larger than 50 cm which corresponds to a minimum age of~125 years (Tang et al. 2012). As far as possible, the populations sampled were located in or near natural habitats, confirmed using additional field investigations or reference to historical litera-ture . Twelve presumably cultivated populations (XA, SX,  ZH, TC, TX, JR, GX, CX, WY, YX, SG, HS) were included as controls, many of which have been genetically characterized with chloroplast markers (Shen et al. 2005;Gong et al. 2008). The latter studies suggest that any human-mediated gene flow (assisted planting) is usually local and would not result in dichotomous genetic structure at a countrywide level.

Molecular protocols
Genomic DNA was extracted from~100 mg of silica geldried leaf material using a modified CTAB protocol from Doyle and Doyle (1987). Three chloroplast DNA (cpDNA) regions were amplified in 354 individuals using the following primers: trnK1 and trnK2 for the trnK intron (Demensure et al. 1995), trnS and trnG for the trnS-trnG intergenic spacer (IGS) (Hamilton 1999), and atpH and atpI for the atpH-atpI IGS (Grivet et al. 2001). We included previously published trnK and trnS-trnG regions for 124 individuals (Gong et al. 2008 Table S5 in Appendix S1. Asterisks indicate cultivated populations. Bayesian reconstruction of phylogenetic relationships (B) among cpDNA haplotypes using Cycas taitungensis as the outgroup (posterior probabilities are indicated for branches with >50% support). A statistical parsimony network is shown immediately below (C). In the latter, pie size is proportional to haplotype frequency, open circles represent missing haplotypes, and the asterisk indicates the most ancestral haplotype (with the highest probability). The dotted line indicates a potentially closed haplotype "loop." Note, however, that haplotypes H12 and H8, respectively, connect to H1 and H6 with high certainty rather than to the missing haplotypes (see text). tems, Foster City, CA) using a total volume of 50 lL containing 3 lL 109 PCR buffer, 2 mmol/L MgCl 2 , 0.2 lmol/L of each dNTP, 0.2 lmol/L of each primer, 1.25 U Taq polymerase (Takara, Dalian, Liaoning, China), and 50 ng genomic DNA. Cycling conditions were 94°C for 5 min, followed by 35 cycles of 94°C for 1 min, annealing for 1 min (62°C for trnK, 58°C for trnS-trnG, and 52°C for atpH-atpI), and 72°C for 1.5 min (with a final extension step of 72°C for 10 min). Amplicons were purified using PCR Clean-Up Kit (Bio-engineering, Shanghai, China). Both strands were sequenced using Taq dye deoxy terminator cycle sequencing kit on an ABI 377XL DNA Sequencer (Applied Biosystems). Sequence alignment was performed using GENEIOUS 4.8 ) with visual inspection and manual adjustment. Unique haplotype sequences were deposited in GenBank (accession numbers of newly sequenced samples: KF220583-KF220586; and retrieved: EF468632-EF468642) ( Table S2 in Appendix S1). All 761 individuals were genotyped at 14 nuclear microsatellite (nSSR) loci (Yan et al. 2006(Yan et al. , 2009 (Table S3 in Appendix S1). PCR amplification was performed in a reaction volume of 10 lL containing 1 lL of 109 PCR buffer, 1 mmol/L dNTPs, 0.5 lmol/L each primer, 0.5 U of Taq polymerase, and 30-50 ng genomic DNA. The forward primers were labeled with a fluorescent dye (6-FAM or HEX; Applied Biosystems). Amplification conditions were as follows: 94°C for 5 min followed by 35 cycles of 94°C for 1 min, the optimized annealing temperature (Table S3 in Appendix S1) for 30 sec, and 72°C for 45 sec (with a final extension step of 72°C for 7 min). Fragment analyses were performed on a MegaBACE 1000 autosequencer (GE Healthcare Biosciences, Pittsburgh, PA). To reduce scoring errors, allele sizes were manually scored twice with reference to an internal size standard (ET-ROX 400) using Genetic Profiler 2.2 (GE Healthcare Biosciences).

Plastid DNA sequence analyses
Each indel and inversion was treated as a single mutation, and coded as a nucleotide substitution, after Caicedo and Schaal (2004). Plastid DNA haplotypes were identified using DnaSP 5.10 (Librado and Rozas 2009) and unique haplotypes from the three cpDNA regions differed by at least one mutation from other sequence variants (Table S4 in Appendix S1). Haplotype diversity (h) and nucleotide diversity (p) were estimated for each population (h S , p S ) and at the species levels (h T , p T ) using the same program.
To quantify the proportion of total genetic variance explained by differences between regional groups of populations, we applied analyses of molecular variance (AMOVAs) based on Tamura and Nei's distance (1993) with 10,000 permutations using ARLEQUIN 3.5 (Excoffier and Lischer 2010).
Phylogenetic relationships among cpDNA haplotypes were reconstructed using Bayesian inference implemented in MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003) with Cycas taitungensis as the outgroup (Wu et al. 2013; Gen-Bank accession No. NC_009618). The best fitting substitution model, HKY + G, was determined according to Akaike Information Criterion (AIC) using jModelTest 0.1 (Posada 2008). Three independent runs of Metropoliscoupled Markov chain Monte Carlo (MCMC) analysis were executed, each including one cold chain and three incrementally heated chains that started randomly in the parameter space. Each run included ten million generations with tree sampling every 1000 generations and a burn-in of the first 20% of sampled trees. The convergence of chains was checked using Tracer 1.5 (Rambaut and Drummond 2009). The remaining trees were pooled to estimate the posterior probabilities. Genealogical relationships were also inferred by constructing a statistical parsimony network using TCS 1.21 (Clement et al. 2000). Because indels may be more prone to homoplasy (Kelchner 2000), we performed separate analyses by including and excluding indel sequences. The resulting network patterns were congruent.

Microsatellite analyses
Linkage disequilibrium and Hardy-Weinberg equilibrium (10,000 dememorization steps, 100 batches, and 5000 iterations) were tested for each locus-population combination using GENEPOP 4.2 (http://www.genepop.curtin.edu.au/). Significance levels were adjusted using the sequential Bonferroni correction for multiple comparisons (Rice 1989). Genetic diversity, as measured by the observed number of alleles (N A ), allele richness (A R ), expected heterozygosity (H E ), and observed heterozygosity (H O ), was estimated per locus and per population using FSTAT 2.9.3 (Goudet 2001). Rare allelic richness (R A ) was manually calculated following Hartl and Clark (1997).
Genetic structure was examined in both cpDNA and nSSR data sets using AMOVA based on the sum of squared size differences (Michalakis and Excoffier 1996) performed in ARLEQUIN 3.5. Genetic structure was further identified in the nSSRs using the Bayesian MCMC programs STRUCTURE 2.3.3 (Pritchard et al. 2000) and BAPS 6.0 (Cheng et al. 2013). The most likely number of clusters (K) was inferred using Evanno's DK method (Evanno et al. 2005)  local optima, ten independent runs for each K value were performed, assuming the admixture model with correlated allele frequencies and using burn-in and MCMC lengths of 100,000 and 1,000,000 iterations, respectively. The BAPS analysis was conducted to incorporate geographical information using prior information on spatial coordinates via the spatial clustering module. STRUCTURE and BAPS vary in their approaches to estimating distinct gene pools and how much ancestry individuals draw from each pool. STRUCTURE infers the highest likelihood for each gene pool and the admixture of genotypes using allele frequency and linkage disequilibrium information from the data set directly. BAPS, on the other hand, first infers the most likely number of gene pools in the sample population and then performs the most likely admixture of genotypes (Corander et al. 2004). This approach has more power in identifying hidden structure within populations (Corander and Marttinen 2006).
Demographic parameters associated with regional populations Two regional populations were defined according to (1) the genetic clusters revealed in BAPS and (2) the phylogeography of cpDNA haplotypes (see results). The twelve known cultivated populations were excluded from this analysis. F ST between these two groups was calculated for each cpDNA locus using DnaSP and for each nSSR locus using Microsatellite Analyzer (MSA) 4.0.5 (Dieringer and Schl€ otterer 2003).
We used an isolation-with-migration (IM) model (Nielsen and Wakeley 2001), implemented in IMa2 (Hey 2010) to estimate the population parameters associated with a split model in which two daughter populations (of size h W and h E ) diverge from a common ancestor (of size h A ) at time t and continue to exchange migrants at rates m 0-1 and m 1-0 . As the complete data set was too large to analyze using an MCMC approach, we generated a smaller subset consisting of 25 individuals from each regional population (by random sampling of both the nuclear and the chloroplast data sets). IMa2 was also used to conduct likelihood ratio tests of nested models to test whether models with unidirectional gene flow are a significantly better fit to the data than models with/without bidirectional gene flow by the method of Nielsen and Wakeley (2001).
Note that to accommodate the uncertainty in mutation rates, we allowed for relaxed intervals, 2.2 9 10 À10 -1.7 9 10 À9 (Graur and Li 2000). Mutation rate priors for nSSR were derived from Thuja plicata (6.3 9 10 À4 , range 3.0 9 10 À5 -4.0 9 10 À3 mutations per locus per generation; O'Connell and Ritland 2004). One hundred geometrically heated chains (using the heating parameters h a = 0.99, h b = 0.75) were run for 2,500,000 steps beyond a 500,000-step burn-in, and trees were sampled every 100 steps. Four runs with different seed numbers were implemented independently.
Because results from each replicate were similar, 10 5 trees were concatenated into a single run in load-tree mode and the "test nested models" option was activated. This option evaluates the likelihood of 24 models simpler than the full IM model by constraining parameters (other than divergence time) and rejecting those that are significantly lower scored than the full model based on a likelihood ratio test. To convert the parameter estimates to demographic units, we assumed a generation time of 160 years for Ginkgo. This is in line with the average age (as a proxy for generation time) of one refugial population (from which our population JF was sampled) based on a direct measurement of increment cores of 46 Ginkgo trees with DBH greater than 5 cm (Tang et al. 2012). Note that simple demographic calculations of generation time might bias the timing estimation and we apply it here with some caution. This generation time is comparable to the mean age (~200 years, DBH > 10 cm) of another temperate tree, Quercus robur, again using age as a surrogate for generation time (Ranius et al. 2009). Increasing the generation time to 210 years (the modal age of population JF) had negligible effects on parameter estimation (data not shown).

Tests for departures from neutrality
ARLEQUIN was used to perform two neutrality tests: Tajima's D (Tajima 1989) and Fu's F s (Fu 1997). Critical values for these statistics were obtained using 100,000 coalescent simulations. Mismatch distribution analyses (MDAs) were conducted using ARLEQUIN, and a goodness-of-fit based on the sum of squared deviations (SSD) and Harpending's (1994) raggedness index (H Rag ) was tested using parametric bootstrapping (Schneider and Excoffier 1999) with 1000 replicates.

Results
Rooting all chloroplast polymorphisms reveals H1 as the ancestral haplotype The reconstructed Bayesian tree of 13 cpDNA haplotype sequences indicates the monophyly of all the observed haplotypes from central and western China (except H1 and H3) (86% Bayesian posterior support, Fig. 1B). The remaining haplotypes are not grouped into one single clade. The parsimony network (Fig. 1C) provided better resolution, although the network contained a potentially closed loop. This ambiguity can be easily resolved, however, as the haplotype in question (H12) connects with high certainty to the co-occurring H11 haplotype rather than to the alternative H8 haplotype from a remote population (Crandall and Templeton 1993). The widespread haplotype H1 was inferred to be the most ancestral, giving rise to two groups of haplotypes (Fig. 1C).
Chloroplast sequences also indicate that more than half of the analyzed populations (14/25) are polymorphic (Table S5; Fig. 1A) with the majority occurring in the west (9/10 populations). Central and eastern populations were significantly less polymorphic (P < 0.05, Wilcoxon sign test) compared to western ones with only 3/8 and 2/ 7 populations being polymorphic, respectively. Nine of the twelve known cultivated populations were monomorphic (except XA, GX, and CX) and fixed with the widespread haplotype H1.

Maternally and biparentally inherited markers exhibit strong west-east population structure
The majority of derived cpDNA haplotypes show westeast geographical divisions ( Fig. 1; Table 1). Haplotypes H3-H9, for example, are specific to central and western China, while H2 and H10-H13 are distributed in eastern China (one exception; H12 does occur in the western cultivated population XA, however). The widespread ancestral haplotype H1 is present in nearly all sampled populations across Ginkgo's natural range.
Bayesian-based mixture analysis of all populations (based on the nSSR data) using prior information on spatial coordinates (BAPS) revealed that the best genetic partition within the data had a clear west-east spatial component (Fig. 2B), supporting the dichotomous west-east division inferred from both the haplotype network and the AMOVA.
Similarly, the STRUCTURE analysis revealed that DK peaked at K = 2 (Fig. S1 in Appendix S1). Individuals from western China (QC, JF, PX, SP, WC, ES, HP) were assigned to one cluster ("green"), while those from the eastern population (TM) were assigned to a different cluster ("red") ( Fig. 2A,C). Individuals from central China were variously assigned to one or other of these two genetic clusters.
AMOVA suggests that genetic variation is significantly structured at the regional (west/east) level for both the maternally inherited cpDNA (Ф CT = 0.166, P < 0.001) and biparentally inherited nSSR (R CT = 0.025, P < 0.001; Table 1). Within either the "west" or "east" group, however, most of the variation is partitioned within populations (~53% cpDNA and~88% nSSR, respectively), with a significant component of variation among populations in these two regions.

Recent population split and asymmetrical gene flow
We fitted an isolation-with-migration model to the west/ east population structure apparent within Ginkgo in order to estimate the split time from a common ancestor for these two regional groups. The estimated divergence time for the split between the two groups identified in the BAPS and STRUCTURE analyses was as recent as 51,355 years ago (95% HPD range: 17,103-95,000) (Fig. 3). Interestingly, following the split, the two regions have exchanged genes asymmetrically. Gene flow from west to east was 2.502 population migrants per generation (0.456-4.525), while gene flow in the converse direction was not statistically significant, 0.889 per generation (0-1.889). A model which assumes unidirectional gene flow between the two regions (from west to east) was a significantly better fit than models without gene flow or gene flow in the opposite direction (from east to west, Table S6 in Appendix S1). Φ is a fixation index similar to Wright's F-statistics. It reflects the correlation of random pairs of haplotypes drawn from a group relative to the correlation of pairs of random haplotypes drawn from the whole data set (Excoffier et al. 1992). Significance is indicated as follows: ***P < 0.001.

No footprint of a genetic bottleneck
From the IMa2 analysis, we also inferred extremely small effective population sizes for the derived (west and east) daughter populations (in the order of 100's), which is an order of magnitude smaller than their census sizes. The magnitude of these effective population sizes may be indicative of demographic changes and thus genetic bottlenecks. The latter, for example, may have resulted in the loss of (rare) genetic variants; detectable in a skew in allele frequency spectra. Nevertheless, G. biloba shows comparatively high levels of both cpDNA and nSSR diversity, which vary minimally among natural populations (Table S5 in Appendix S1). Analysis of the nucleotide frequency spectra in cpDNA indicated no significant signal of population expansion; Tajima's D was slightly negative but not significantly so (Table S7 in Appendix S1).  Table S5 in Appendix S1. Asterisks indicate cultivated populations. Effective population size

Discussion
By comparing population genetic data from Ginkgo biloba with our (in-house) published results of Cercidiphyllum japonicum ) together with published Asian case studies, we provide some robust interpretations below of the hypotheses concerning variation in the degree of postglacial expansion of warm-temperate deciduous forests in eastern Asia. The detailed dynamic history gleaned from the population genetic analyses allows the investigation of a process that includes isolation-withmigration models, random genetic drift, and the split times of population genetic structure in the context of isolation hypotheses proposed for the range dynamics of flora during glacial and interglacial periods.
Regional Ginkgo population structure results from genetic drift in the late-Pleistocene The IMa2 analysis suggests significant evidence for unidirectional gene flow from west to east. Nevertheless, significant population differentiation between west and east is observed for each chloroplast locus (F ST = 0.148-0.301, P < 0.001) and for each nSSR locus (0.002-0.081, P < 0.01). Gene flow is expected to reduce or eliminate this differentiation. Low levels of migration (i.e., > one migrant per generation) are sufficient to prevent population differentiation by drift (Slatkin 1987). However, west-east population differentiation is apparent despite asymmetrical gene flow acting as a homogenizing force between these regions. Trees have long generation times, and together with the recent split, it may take many thousands of generations of gene flow to diminish this population structure. To illustrate this point, we ran coalescent simulations (in ms, Hudson 2002) for a two-deme model with asymmetrical migration, using parameters estimated for G. biloba from our IM analysis (Fig. 3). From these simulations, at least 0.5Ne generations of gene flow are required for F ST to reach equilibrium (Fig. 4). The magnitude of simulated F ST apparent between populations at 0.4N e generations ago is comparable to the average F ST (0.14) we observe from our cpDNA haplotype and nSSR data. This timing (0.4N e generations) in turn corresponds to approximately 39,000 years ago (using an average estimate for N e from western and eastern populations and a generation time of 160 years). The simulations suggest that the relatively high F ST we observe between the regional populations dates from~39 k years ago and may mirror the lag before F ST reaches equilibrium between the lineage sorting (random genetic drift) of population divergence and the homogenizing effects of asymmetrical gene flow.

Conflicting support for isolation hypotheses during glacial periods
While the present study confirms the presence of two genetically distinct populations associated with independent glacial refugia identified by Gong et al. (2008), one intention here was to shed light on the demography behind these west-east regional refugia and to identify where the population split (and associated genetic drift) lies in relation to geological time during the last glacial period. The recent split between the two populations at 51 kya (95% HPD: 17,103-95,000; Fig. 3), fits firmly within the last glaciation (10-73 kya) in China (Yi et al. 2005). This suggests that the reduced range (and effective population size) of G. biloba dates, in part, from genetic drift in the late-Pleistocene. Moreover, prior to the split, gene flow among populations across the range might have been higher, and/or the effective population size might have been larger in the historical migrant gene pool than observed after the split. These inferences are supported by Tajima's D values and the mismatch distribution, both of which did not reveal any signals of population expansion (Table S7 in Appendix S1). We may thus infer that the demography of Ginkgo includes recent fragmentation of a formerly larger population with probably long-term in situ persistence (see also Eurycorymbus cavaleriei; Wang et al. 2009). Polymorphism data were generated by coalescent simulations in the ms program (Hudson 2002). The simulations assumed that all demes were of equal size with h = 1.7, and 50 sequences were sampled from each subpopulation, respectively. One subpopulation received migrants at the same rate from the other subpopulation. A migration parameter (4N 0 m) of 5.0 and 1.86 was assumed for subpopulation WEST and EAST, respectively. For each run, gene flow was introduced at time T (measured in N e ) generations ago at intervals of 0.05N e . Dots represent the average differentiation (F ST ; error bars represent standard deviations from 10,000 replicates) measured T generations after gene flow was introduced.
Assuming this inference of gene flow and/or a large effective population size prior to the split is correct, together with the asymmetrical gene flow observed following the split (Fig. 3B), it appears that west-east regional Ginkgo populations have not remained genetically isolated during the last glaciation. There is thus both evidence to support the Harrison et al. (2001) isolation hypothesis during glacial periods (i.e., sufficient isolation for drift to create the regional split) and logical inference to reject it, as there is genetic connectivity prior to the split, and recent asymmetrical gene flow between regions following the split.
This illustrates that isolation hypotheses (with or without gene flow) fail to account for the temporal dynamics of both isolation and gene flow among populations during glacial and interglacial periods. A model which includes isolation with migration is a good fit to the Ginkgo system and yet the coalescent simulations (for a two-deme model with asymmetrical migration) confirm that genetic population structure (and apparent isolation) can persist in the presence of (recent) gene flow. Harrison et al. (2001) and Qian and Ricklefs' (2000) hypotheses suffer by not recognizing the potential for this temporal complexity. Only by unraveling population genetic processes such as those revealed in this study, are these subtleties revealed.

Incongruence of demography between cooccurring species may be facilitated by life history traits
The west-east dichotomous regional genetic structure that we observe in Ginkgo is also apparent in numerous other temperate woody plants from eastern Asia, for example, Taxus wallichiana (Gao et al. 2007), Quercus variabilis , and Kalopanax septemlobus . The genetic structure observed in Cercidiphyllum is approximately north-south ( Fig. 1 in Qi et al. 2012). However, for most documented Asian temperate tree species, the division either lies between central-west and eastern China/Asia (e.g., Ginkgo, Q. variabilis) or shifted clockwise along a northwest-southeast axis between Sino-Himalayan and Sino-Japanese forests (e.g., T. wallichiana, K. septemlobus). Among all of the taxa, the ages of these regional divisions varies substantially, from the youngest in Ginkgo~51 kya to the oldest in Cercidiphyllum 1.89 Ma. In between these ranges, the genetic structure has been dated to 1.45 Ma in Q. variabilis and to 0.28-0.74 Ma in K. septemlobus. Interestingly, particularly for inferences regarding in situ persistence versus expansion hypotheses, postglacial expansions to former range limits differ among these species. Temporal and spatial incongruences suggest that these four trees did not respond to climate-driven vegetation shifts in the same manner. One explanation for the incongruences may be differences in life history traits which we examine in detail below for two of the more ecologically associated taxa, Gingko and Cercidiphyllum.
Given the maternal inheritance of cpDNA in both Ginkgo and Cercidiphyllum, the range expansion observed in Cercidiphyllum following glaciation may imply comparatively poorer seed dispersal in Ginkgo. Indeed, the fleshy (heavy) seeds of Ginkgo are probably dispersed less efficiently than the wind-dispersed seeds of Cercidiphyllum. Knowledge of the seed dispersers in Ginkgo has been speculative (Del Tredici 2007;Del Tredici 1992;Jiang et al., 1990) but some information about the effectiveness of seed dispersers can be inferred from the distribution of maternally inherited chloroplast haplotypes. With the exception of the widespread ancestral H1, no derived cpDNA haplotypes were shared between closely located populations, for example, TM and CX (<100 km apart), suggesting that Ginkgo seeds were not dispersed efficiently. Nevertheless, the widespread H1 and its direct relatives were shared between distantly located populations, suggesting that the commonality among geographically dispersed haplotypes results from shared ancestry rather than animal (or human) mediated gene flow. The split time (~51 ka) for the regional population structure far predates that of human-mediated introductions (<2 ka), and is too recent to have allowed new mutations to generate five newly derived haplotypes specific to eastern China. We thus consider animal mediated gene flow via seeds to be limited.
Ginkgo is also at a comparative disadvantage to Cercidiphyllum in disturbed riparian rocky environments because its large seeds with slowly developing embryos require a level of stability less often available in habitats with frequent disturbance and flooding (Del Tredici 2007). Furthermore, the slow and complex sexual reproductive cycle in Ginkgo is highly constrained by temperature (Hatano and Kano 1952;Del Tredici 2007). As a result, sexual reproduction may have been inhibited during the last glaciation, resulting in range contraction. Limited seed dispersal, poor competition to recolonize riparian habitats, and inefficient sexual reproduction would have severely constrained the plant's ability to spread beyond refugia following range contraction. Of the two species, Gingko may have been more sensitive to climatic cooling and aridification at the beginning of the Pleistocene. Therefore, the influence of climatic niche differentiation in the demographic histories of these two species cannot be discounted.
On balance, seed dispersal coupled with the long generation time in Ginkgo, are important factors in explaining the presence of (restricted) regional Ginkgo populations in China. Life history traits can impose restrictions on the range dynamics and population genetic structure of temperate plants by their influence on lifetime reproductive success (Hamrick and Godt 1996) and should be considered as a valid alternative to previously proposed hypotheses to explain variation in the degree of postglacial expansion to former range limits between cooccurring species of warm-temperate deciduous forests.

Conclusions
Two competing isolation hypotheses (with or without gene flow) have been proposed to explain the variance of postglacial expansion to former range limits among warm-temperate deciduous Asian tree species. During glacial periods, isolation with admixture at lower altitudes (Qian and Ricklefs 2000), or continual isolation without admixture (Harrison et al. 2001; negating any isolation, if 2Nm > 1), may explain this variance. However, proving that isolation occurred continuously throughout a glacial period is challenging. Gene flow and random genetic drift can vary both temporally and spatially throughout glacial periods. Both isolation hypotheses are thus not necessarily mutually exclusive, and they also neglect the effect of life history traits on range dynamics and population genetic structure via their influence on lifetime reproductive success. Using population genetics to uncover the dynamism of postglacial history, this study illustrates the inadequacies of "isolation" hypotheses that span geological time where population dynamics are the norm. The comparison between Ginkgo and other taxa, particularly Cercidiphyllum, sampled across the same ecosystem indicates that co-occurring species may not have congruent range dynamics. This highlights the importance of both demography and life history traits when interpreting range shifts and the evolution of ecological communities.

Supporting Information
Additional Supporting Information may be found in the online version of this article: Appendix S1. Sampling information, GenBank accession numbers, population diversity of nSSRs and cpDNA haplotypes, tests of IM nested models, neutrality tests, and likelihood values of the cluster number in Bayesian assignment analyses. Table S1. Details of sample locations, voucher number, latitude, longitude and altitude in 25 G. biloba populations. Table S2. GenBank accession numbers of cpDNA variants. Table S3. Locus, repeat motifs, size ranges, number of alleles (N A ), observed (H O ) and expected heterozygosity (H E ) for fourteen microsatellite loci analyzed in 25 populations of G. biloba. detected in the three regions of Ginkgo biloba identify 13 haplotypes (H1-H13). Table S5. Plastid and nuclear diversity in 25 G. biloba populations sampled across its complete natural range in China. Table S6. Tests of nested models for the WEST versus EAST group for the three cpDNA loci and 14 nSSR loci data set. Table S7. Summary of neutrality tests and mismatch distribution parameters for groups WEST and EAST of G. biloba. Estimates of mismatch distribution were obtained under models of spatial expansion using ARLE-QUIN. Fig. S1 Likelihood values across multiple values of the cluster number in Bayesian assignment analyses, K (1-10) using Evanno's methods (Evanno et al. 2005).