Spatial and temporal population genetic variation and structure of Nothotsuga longibracteata (Pinaceae), a relic conifer species endemic to subtropical China.

Nothotsuga longibracteata, a relic and endangered conifer species endemic to subtropical China, was studied for examining the spatial-temporal population genetic variation and structure to understand the historical biogeographical processes underlying the present geographical distribution. Ten populations were sampled over the entire natural range of the species for spatial analysis, while three key populations with large population sizes and varied age structure were selected for temporal analyses using both nuclear microsatellites (nSSR) and chloroplast microsatellites (cpSSR). A recent bottleneck was detected in the natural populations of N. longibracteata. The spatial genetic analysis showed significant population genetic differentiation across its total geographical range. Notwithstanding, the temporal genetic analysis revealed that the level of genetic diversity between different age class subpopulations remained constant over time. Eleven refugia of the Last Glacial Maximum were identified, which deserve particular attention for conservation management.

Introduction

Genetic structure of a plant natural populations is mainly influenced by the species breeding system, gene flow, genetic drift and natural selection (Epperson, 1993; Austerlitz ). Furthermore, climatic and/or geological changes can also result in shifts in the genetic composition of a population through creating spatial and/or temporal heterogeneity in the environments (Jump and Penuelas, 2005). Although spatial genetic structure has been investigated in many plant groups at different geographical scales (Erickson and Hamrick, 2003; Geng ; Soldati ), most previous studies did not distinguish among age classes in a population, and the temporal dynamics of gene flow and selection are not well understood, possibly due to divergence among groups of individuals of different ages in their responses to ecological and evolutionary factors (Erickson and Hamrick, 2003; Soldati ). Obviously, the patterns of population genetic structure of adult plants reflect the accumulation of the effects of ecological and evolutionary processes acting in the past and present, while the genetic structure of juvenile plants represents the processes of a few recent years (Kalisz ; Chung ; Ozbek ). Dissecting both the spatial and temporal genetic structure of a population therefore permits a better understanding of the evolutionary dynamics and biogeographic history in space and time.

Severe climatic oscillations associated with glacial cycles during the late Tertiary and throughout the Quaternary resulted in profound changes in species distribution and population structure (Hewitt, 2000; Petit ; Wang and Ge, 2006). Subtropical mainland China, where no ice sheets were formed during the climatic oscillations, had served as an important refugium for lineages that evolved prior to the late Tertiary and Quaternary glaciations (Gao ). In particular, mountain regions such as the Nanling Mountains and the adjacent Wuyi and Daiyun mountains have been proposed as major centers of the great floristic diversity in subtropical China, and an active center for speciation and evolution (Ying, 2001; Wang ). The complex physiogeographic and climatic features in these mountain regions have allowed the survival of many relict species. The current geographical distribution and population genetic structure of these relict species represent the results of the long-term impacts of climatic and geological events together with the long or short-term effects of demographical processes related to anthropogenic disturbances. Thus, detailed assessment of the population genetic structure of these relict species in space and time is crucial to a comprehensive understanding of the genetic and biogeographical patterns of species endemic in the mountain regions of subtropical China.

Nothotsuga longibracteata (W. C. Cheng) Hu ex C. N. Page is a relict tree species belonging to the genus Nothotsuga in the family Pinaceae. Although it has first been placed in the genus Tsuga, the species is notably distinguished from Tsuga by its radially spreading leaves with stomatal lines on both upper and lower surfaces and the upright seed cones. Most natural populations of this species are currently centered in several highly fragmented and isolated mountain areas of south China, including the Nanling mountains locating in the cross-border region of the Guangdong, Jiangxi, Hunan, Guangxi and Guizhou provinces and the Daiyun mountains in the Fujian province. Notwithstanding, fossil records of N. longibracteata found in Russia, Germany and Japan suggest a wide geographical distribution of this species during the Pliocene and Miocene (Karavaev, 1958; Miki, 1954). The climatic and environmental changes during the Quaternary glaciations and the increasing anthropogenic disturbance and habitat fragmentation during the past 100 years may have severely decreased the natural populations of N. longibracteata. A recent investigation of the demographic composition of several undisturbed populations in the Fujian and Guangdong provinces showed that the age spectrum of this species is large, ranging from one year old seedlings to > 500 years old trees. Thus, N. longibracteata provides an opportunity to investigate both the spatial and temporal genetic structure of populations in response to climatic and geological changes in the past together with recent anthropogenic disturbance and local habitat fragmentation.

In the present study, we used biparentally inherited nuclear microsatellites (nSSRs) and paternally inherited chloroplast microsatellites (cpSSRs) to screen the genetic variation in natural populations of N. longibracteata. Ten populations were sampled across their entire natural range for spatial genetic structure analyses, and three key populations with large population sizes and age classes were selected for temporal genetic analyses. We were particularly interested in the following questions: (1) the degree of spatial genetic differentiation among the populations of N. longibracteata throughout the entire natural distribution in China; (2) whether there is a significant difference in gene diversity between different age class subpopulations; and (3) the possible location of N. longibracteata during the Last Glacial Maximum.

Materials and Methods

Sample collections

We conducted an exhaustive field survey for the entire natural range of N. longibracteata in 2005. The geographic distribution of N. longibracteata is scattered and fragmental (Qiu ), and only 10 natural populations are relatively large (with more than 30 individuals). A total of 631 individuals were sampled from these 10 populations across their entire geographical distribution (Figure 1). The geographical distances between populations varied from 1.3 km to 700 km. About 30 individuals per populations were collected depending on the accessibility of samples and the original population size (Table 1). Moreover, three populations, FJM (149 individuals), FJT (160 individuals) and HNZ (106 individuals), with large population sizes and broad age classes were exhaustively sampled for temporal genetic structure analyses by grouping individuals of different age class into subpopulations (Table 2). The age of individuals in the three large populations was determined using the methods described by Wu . Needle leaves from each sampled individual were collected and immediately buried in silica gel until DNA extraction.

Figure 1

Population location and geographical distribution and frequencies of cpDNA SSR haplotypes of Nothotsuga longibracteata in subtropical China.

Table 1

Genetic diversity within populations in Nothotsuga longibracteata.

Region (abbr.)	Population	N	A	AR	HO	HE	FIS	NP
Daiyun mountain(DY)	FJM	32	11.3	10.6	0.748	0.849	0.109	9
Daiyun mountain(DY)	FJT	32	10.3	9.8	0.708	0.845	0.153***	4
Nanling mountain(NL)	HNZ	32	9.8	9.1	0.807	0.797	−0.029	10
Nanling mountain(NL)	HNG	32	8.5	8.1	0.743	0.770	0.036	4
Nanling mountain(NL)	GDL	31	10.7	10.2	0.767	0.851	0.084	9
Jiangxi province(JXS)	JXS	32	5.0	5.0	0.846	0.761	−0.119	0
Hunan province(HNQ)	HNQ	32	2.0	2.0	0.451	0.413	−0.105	1
Hunan province(HNH)	HNH	26	1.7	1.7	0.310	0.291	−0.022	1
Fanjing mountain(GZF)	GZF	31	9.2	8.8	0.545	0.760	0.274***	11
Maoer mountain(GXM)	GXM	32	8.3	7.9	0.697	0.710	0.018	8
Average		31	7.6	7.3	0.658	0.704		5.7

N, number of individuals sampled for each population; A, average number of alleles per locus; AR, allele richness; HO, observed heterozygosity; HE, expected heterozygosity; FIS, within-population coefficient of inbreeding; NP, number of private alleles. Levels of significance for FIS of each population were determined after 1000 permutations:

, p < 0.001.

Table 2

Summary of genetic diversity measures for each age classes of Nothotsuga longibracteata at FJM, FJT and HNZ populations.

Stand	Age class	Age classes analysis						Hierarchical age class analysis
		Age range (yr)	N	A	AR	HO	HE	Age range (yr)	N	A	AR	HO	HE
FJM	I	seedling	38	11.0	8.3	0.775	0.840	> 0	149	15.7	9.1	0.755	0.844
	II	50–100	46	10.5	7.7	0.726	0.815	> 50	111	14.7	8.9	0.749	0.835
	III	100–150	14	7.5	7.5	0.774	0.835	> 100	65	14.0	9.3	0.765	0.845
	IV	150–200	18	9.2	8.6	0.795	0.818	> 150	51	13.0	9.3	0.764	0.847
	V	200–500	18	10.0	9.2	0.741	0.857	> 200	33	11.5	9.3	0.747	0.851
	VI	> 500	15	9.2	8.9	0.756	0.839	> 500	15	9.2	9.2	0.756	0.839
FJT	I	1–50	25	10.2	7.0	0.719	0.850	> 10	160	14.5	9.6	0.708	0.851
	II	50–100	93	13.8	7.3	0.696	0.846	> 50	135	14.3	9.5	0.706	0.846
	III	100–150	27	10.2	6.9	0.752	0.845	> 100	42	10.5	8.8	0.729	0.843
	IV	> 150	15	7.7	6.2	0.687	0.817	> 150	15	7.7	7.7	0.687	0.817
HNZ	I	< 50	60	9.5	8.9	0.744	0.757	> 0	106	10.8	9.3	0.753	0.770
	II	> 50	46	9.3	9.3	0.764	0.780	> 50	46	9.3	9.3	0.764	0.780

N, sample size; A, average number of alleles per locus; AR, allele richness; HO, observed heterozygosity; HE, expected heterozygosity.

DNA extraction and microsatellite analysis

Genomic DNA was extracted from needle leaf tissue using the CTAB method (Doyle and Doyle, 1987). DNA quality and quantity were determined by electrophoresis in 1% agarose gels with λDNA markers. Microsatellite genotyping of the nuclear genome of N. longibracteata was performed according to the methods described by Qiu using the loci NT01, NT02, NT03, NT04, NT06 and NT07. In addition, three cpDNA microsatellite primer pairs (Pt15169, Pt63718, Pt71936) derived from Pinus thunbergii (Vendramin ) were used due to their high genetic polymorphism found in a preliminary screen for genetic variation of cpDNA in N. longibracteata. Polymerase chain reaction (PCR) assays were done in a volume of 10 μL, containing 10 mM Tris-HCl (pH 8.4), 50 mM (NH4)2SO4, 1.5 mM MgCl2, 0.2 mM dNTPs, 0.2 μM of each primer, 50 ng of genomic DNA, and 1 unit Taq polymerase (Fermentas, Lithuania) for both nSSR and cpSSR. The amplification protocol was essentially similar to that described by Qiu for SSR and by Vendramin for cpSSR. The amplified products were separated in a 6% denaturing polyacrylamide gel and revealed using silver staining. A 25 bp DNA ladder (Promega, Madison, WI, USA) was used to estimate allele size.

Data analysis

Population genetic variation

For nuclear genomic SSR markers, a set of standard measures of population genetic diversity were calculated for all sampled populations using the program Microsatellite Analyzer (MSA) 4.05 (Dieringer and Schlotterer, 2003), including the mean number of alleles per locus (A), allele richness (AR) standardized to the smallest sample size, observed (HO) and expected heterozygosity (HE) for each locus and population, and inbreeding coefficient (FIS) for each population. The population genetic differentiation for each locus across all populations was estimated by Wright’s F-statistics (FST) according to Weir and Cockerham (1984) using the program GenALEx6 (Peakall and Smouse, 2006). Departure from Hardy-Weinberg expectations and linkage disequilibria between microsatellites were tested by Fisher’s exact tests implemented in GENEPOP, version 3.4 (Raymond and Rousset, 1995). Sequential Bonferroni adjustments were used to determine the significance levels of all statistics, with an initial ± level of 0.05.

In the three populations (FJM, FJT and HNZ) with large population size, each age class was considered as an individual temporal subpopulation and the genetic diversity measures (AR, HO and HE) were calculated separately for each subpopulation. Furthermore, the temporal subpopulations were redefined based on the cumulative gene pools of the different age classes (Table 2), and a hierarchical analysis was used to test the difference of genetic diversity (AR, HO and HE) between subpopulations under the assumption of random mating. Differences in the genetic diversity parameters among different age class subpopulations and among hierarchical age-class subpopulations were analyzed by ANOVA, respectively. FST values between each age class subpopulation and among hierarchical age-class subpopulations (as show in Table 2) were also examined to give a picture of how genetic structure changed over time.

cpDNA haplotypes were identified based on cpSSR variation. Length variants at each cpSSR locus were combined into one haplotype because of the non-recombinant nature of the chloroplast genome. The program ARLEQUIN (Excoffier ) was used to calculate haplotype diversity (H) (Nei, 1987) and the number of private haplotypes for each population and group (groups defined in the analysis of population genetic structure-see below). The global haplotype differentiation value GST was calculated using the program PERMUT (Pons and Petit, 1996; Burban ).

Migration-drift equilibrium analyses

The relative likelihood of a migration-drift equilibrium (gene flow model vs. drift model) between populations in different geographical localities was investigated using the program 2MOD based on nuclear SSR data sets (Ciofi ). 2MOD is based on the comparison of likelihoods for the two models, given the observed microsatellite frequency counts using coalescent theory and Markov Chain Monte Carlo simulation. The gene flow model assumes that gene frequencies within populations are determined by a balance between genetic drift and immigration, while the drift model assumes that an ancestral panmictic population became separated into several independent units which then started to diverge purely by genetic drift. A simulation with 500 000 iterations was carried out, and the first 10% of the output was discarded in order to avoid a bias resulting from the starting values.

Population genetic structure

The program STRUCTURE version 2.1 (Pritchard ) was used to infer the potential number of source population for all individuals sampled based on the nSSR data set. We examined the probabilities for a range of K values starting from 1 to 17. A burn-in of 106 and a Markov chain Monte Carlo run length of 500,000 iterations were used for each run of a K value to reach an approximate stabilization of the summary statistics based on the admixture ancestry model with correlated allele frequencies. Runs for each K value were independently replicated 20 times (Gilbert ). The outputs of Ln P(D) were then standardized using the method described by Evanno . An admixture analysis was next carried out for the K-inferred populations regardless of sampling populations in different locality.

Analysis of molecular variation (AMOVA) was used to determine the partitioning of molecular variance between and within populations for both the nSSR and cpSSR data sets using the program ARLEQUIN 3.11. Furthermore, population pairwise comparisons of FST were conducted using a genetic distance approach based on nSSR data under the infinite alleles model (IAM). This model assumes that genetic drift is the main driving force underlying population divergence (Weir and Cockerham, 1984; Michalakis and Excoffier, 1996). Significance of the variance components was tested with 10,000 permutations.

Results

Informative genetic markers and population genetic diversity

A total of 143 alleles from six nuclear microsatellite loci were found and the mean number of alleles per locus was 23.8, ranging from 6 (locus NT01) to 38 (loci NT03 and NT06). Moreover, mean gene diversity per locus (HE) varied from 0.093 (locus NT07) to 0.777 (locus NT03), with an average value of 0.605 (data not shown). The general genetic diversity measures of each population are summarized in Table 1. Intra-population microsatellite variation revealed the mean number of alleles per locus A = 7.6 (range from 1.7 to 11.3), allelic richness AR = 7.3 (range from 1.7 to 10.6), HO = 0.658 (range from 0.310 to 0.846), and HE = 0.704 (range from 0.291 to 0.851). High genetic diversity together with a high number of private alleles (rare alleles found only in one or two samples) were found in four populations of DY, NL, GZ F and GXM (Table 1). The level of inbreeding for each population also varied from −0.119 to 0.274 (Table 1). Significant inbreeding coefficients (FIS) were detected in two populations, FJT and GZF, and they were under Hardy-Weinberg disequilibrium (p < 0.05) across all loci analyzed. A significant linkage disequilibrium (p < 0.05) between SSR loci within each population was detected for 20 out of 150 comparisons. Since multiple tests were involved, a sequential Bonferroni correction was conducted by dividing the significance level by the number of items tested. As a result, none of these linkage disequilibria were significant after Bonferroni correction.

For chloroplast SSR markers, a total of 11 alleles were identified from the three cpDNA microsatellites, and genetic diversity across cpSSR loci varied greatly among groups (Table S1). The highest haplotype diversity He was found in the NL group (0.348), followed by the DY (0.271), GXM (0.174) and GZF (0.157) groups, with a mean He of 0.136 (Table S1). Thirteen different haplotypes were defined based on the combination of the 11 alleles from the three loci (Table S2). Haplotype frequencies varied from 0.3% (one individual for h13) to 62.1% (195 individuals for h11), with an average of 7.7%. The NL group had the highest number of haplotypes (8), followed by the DY (6), GZF (6) and GXM (4) groups. Interestingly, one haplotype (h11) occurred in every group (Figure 1). Moreover, three haplotypes (h6, h7and h8) uniquely occurred in the NL group, and three (h4, h10 and h13) were in the DY, GZF and GXM regions, respectively (Figure 1, Table S1).

Migration-drift equilibrium

The result of likelihood tests for the gene flow vs. drift models in the 2MOD analysis was strongly in favor of the gene drift mode (p (gene drift) = 1, Bayes factor = ∞), suggesting that the extant populations of N. longibracteata are genetically isolated across a large-scale spatial region and that there is only a small amount of migration between populations.

Population spatial genetic variation and structure

The assignment analysis implemented in the STRUCTURE program indicated that there are 11 fragmented populations in our sampling scheme (Figure 2), suggesting that there could be as many as 11 locations of populations during the Last Glacial Maximum. All pairwise comparisons of populations genetic differentiation (FST) were significant (p < 0.001). Moreover, of the total genetic variation partitioned in the AMOVA analysis based on the nSSR data set, the largest variance component (81.39%) was found among individuals within populations, while 11.52% and 7.1% variance were found between regions and among populations within the same region, respectively (Table 3), suggesting that the populations shared similar gene pools (common ancestry) until very recently. Similar results were also obtained from the cpSSR data, in which the genetic variation within populations, between regions, and among populations within region were 71.75%, 17.96% and 10.29%, respectively (Table 3).

Figure 2

Genetic structure of N. longibracteata populations as defined by STRUCTURE (K = 11) using prior information of population origin.

Table 3

Results from analysis of molecular variance (AMOVA) from nSSR and cpSSR markers for Nothotsuga longibracteata populations.

Analysis	Population	marker	Source of variation	d.f.	Percentage of variation	p value
Spatial	All	nSSR	Among populations	9	18.25	p < 0.01
			Within populations	614	81.75	p < 0.001
	All	cpSSR	Among populations	9	23.57	p < 0.001
			within populations	304	76.43	p < 0.001
Temporal	FJM	nSSR	Among age classes	5	1.89	p < 0.001
			Within age classes	292	98.11
	FJT	nSSR	Among age classes	3	1.8	p = 0.295
			Within age classes	331	98.2
	HNZ	nSSR	Among age classes	1	0.79	p < 0.05
			Within age classes	210	99.21

df = degrees of freedom.

In the entire geographical range of N. longibracteata sampled, no cpSSR genetic variation was found in the three populations JXS, HNQ and HNH (each population was fixed with a single haplotype) (Table S1). The global haplotype differentiation value of GST was 0.275. In the seven defined different geographical groups, the NL group showed the highest level of cpSSR diversity (with 8 haplotypes). Among the eight haplotypes in the NL group, h1 and h3 were shared with group DY in the east and h5 was shared with groups GXM and GZF in the west, while h2 was shared with group GZF (Figure 1).

Genetic diversity and differentiation of age-class subpopulations over time

Based on the ANOVA analysis result, the levels of nSSR genetic diversity (HE, HO or AR) were not significantly different (p = 0.372–0.972, data not shown) between the age-class subpopulations in the three large populations FJM, FJT and HNZ (Table 2). Furthermore, hierarchical age-class analysis also showed no significant differences in genetic diversity (p = 0.724 – 1.000, data not shown) for the cumulative gene pools of age-class subpopulations (Table 2).

Differentiation of FST between hierarchical age-classes was not significant among all the three populations (FJM, FJT and HNZ, data not shown). The AMOVA analysis showed that > 98% of the total variation occurred within the age class subpopulations (Table 3). This result was further confirmed when computing pairwise FST between age classes. In these subpopulations with juveniles (FJM, FJT and HNZ) and the ones with old trees (FJM and HNZ) were significant different in genetic differentiation (Table S3).

Discussion

Genetic diversity in N. longibracteata

Compared with the genetic diversity previously reported in other conifer species investigated by means of nuclear microsatellites markers, a slightly higher level of expected heterozygosity (HE = 0.291–0.851, mean = 0.704) was detected in N. longibracteata compared to Pinus strobus (Marquardt and Epperson, 2004), Picea abies (Meloni ), Taxus wallichiana var. mairei (Zhang and Zhou, 2013), Araucaria cunninghamii, Wollemia nobilis and Agathis robusta (Peakall ). But it was lower than in P. abies (Achere ), Tsuga mertensiana (Ally and Ritland, 2007) and P. pinaster (Mariette ). Notwithstanding, the allelic richness of this species (AR = 1.7–11.3, mean = 7.3) was much lower than in the other conifer species, including P. strobus (Marquardt and Epperson, 2004), P. abies Karst. (Meloni ), P. abies (Achere ) and P. pinaster (Mariette ), except for A. cunninghamii (Peakall ). Both theoretical and empirical evidence indicate that AR is more sensitive than HE to reflect the effects of short but severe genetic bottlenecks, and that it can more accurately represent current levels of genetic diversity within a fragmented population (Leberg, 2002). Accordingly, the relatively low value of AR but high HE suggest that recent bottlenecks have once occurred in the natural populations of N. longibracteata, probably caused by the great decline in the number of populations and population size, as evidenced by the fossil records (Wang, 2000). Fossil evidence indicates that N. longibracteata was widespread in Russia and Japan during the Pliocene (Karavaev, 1958; Miki, 1954) and Germany in the Miocene (Kunzmann and Mai, 2005), although its current geographical distribution is restricted to a few highly fragmented forest areas in subtropical China.

Chloroplast markers are good indicators of genetic bottlenecks due to their smaller (half) effective population size that is more sensitive to genetic drift and loss of alleles than the nuclear markers (Provan ; Zhao and Gong, 2012; Leigh ). In our study, the genetic diversity of cpDNA detected in N. longibracteata populations was lower (mean HE = 0.136) than in other conifer species (HE = 0.37, based on the average estimates of cpDNA genetic diversity measured from cpDNA microsatellite markers in eight coniferous species, Petit ), suggesting a possibly historical bottleneck in the species. Loss of chloroplast variation due to historical bottlenecks was also reported in other conifers (Provan ; Walter and Epperson, 2001; Vendramin ) and it is most likely the case for N. longibracteata.

Population genetic structure

Spatial

Conifer species generally display low population genetic differentiation due to the widespread occurrence of wind pollination and outcrossing (Hamrick ; Hamrick, 2004; Parchman ). For example, a mean GST value of 0.165 has been reported for 37 conifer species based on paternally inherited markers, and one of 0.116 for 33 conifer species based on biparentally inherited markers (Petit ). In the present study, a slightly higher level of the population genetic differentiation (GST value of 0.275 for cpSSR and 0.193 for nSSR) was found in N. longibracteata than those previously reported in other conifer species. Moreover, our AMOVA analysis indicated that a significant amount of genetic diversity was attributed to variation among populations (18.61% and 28.25% of the global genetic diversity between populations for nSSR and cpSSR, respectively, Table 3). Similar patterns of genetic structure have also been documented in Cathaya argyrophylla, an endangerd conifer species in subtropical China (22% of the total diversity at nuclear loci is attributed to variation among populations, Wang and Ge, 2006). This value is significantly higher than the average (GST = 0.116) summarized by Petit . The seeds of N. longibracteata are winged, and both pollen and seeds are wind-dispersed, suggesting a possibly high level of pollen and/or seed flow within and between populations. Thus, instead of the limitations of recent gene flow, the observed population genetic differentiation is more probably a result of the long-term evolutionary divergence of the populations in response to local environmental conditions and/or geological and climatic changes in the recent past. Furthermore, considering the low evolutionary rate of chloroplast DNA, the relatively high level of cpSSR genetic differentiation found between populations and regions agrees with the scenario of ancient genetic isolation. Fragmentation of populations, together with complex topographic features of mountains in the region, might constitute significant barriers to historical gene flow, and in turn might have promoted deep divergence of geographically isolated populations of N. longibracteata.

Temporal

Populations of long-lived organisms, composed of cohorts established at different times and occupying relatively large geographical areas, can be genetically differentiated both spatially and temporally (Linhart ). In our study, however, no significant difference was found among the estimates of genetic diversity between age classes in the three large populations FJM, FJT and HNZ (Table 2). Furthermore, the hierarchical age class analysis showed that the gene pools of the different age-class subpopulations did not change over time (Table 2). This was different from certain previous studies. For example, Ortego found that fragmentation contributed to reduce genetic variability and increase genetic differentiation in Quercus ilex saplings, indicating that the younger cohorts are suffering some negative genetic consequences of long-term population fragmentation. It was also found that elevated inbreeding occurs within seedling cohorts and elevated divergence occurs initially between seedling cohorts and older individuals, but that these differences decline in older cohorts due to selecting out of inbred individuals (Kitamura ; Chung ; Kelly ). When spatial genetic structures of species are profoundly shaped by global climate oscillations in a long time scale (Hewitt, 2000; Godbout ; Espíndola ), the temporal genetic variation generally occurred in a short time period due to local changes of land-use or natural fluctuations in the abundance and distribution of species (Parmesan and Yohe, 2003). Obviously, conifers combine various life history traits that promote constant genetic diversity under local changes, including large populations with wide geographical distribution, long-lived perennials, sexual reproduction, high fecundity, outcrossing, and relatively long distance pollen and seed movement (Hamrick ; Ledig, 1998; Parchman ). The lack of a temporal pattern indicates that N. longibracteata maintained a constant level of genetic diversity and structure over time, implying no disastrous environmental changes that selectively removed plants of certain genotypes in a short term.

Notwithstanding, we did found weak but significantly genetic differentiation among the young and most of the old cohorts of temporal groups in two of the three populations (FST of 1.9% for FJM, p = 0.000, and of 0.8% for HNZ, p = 0.028) (Tables 3 and S1). Assuming neutrality of microsatellite markers, the temporal changes in genetic composition of N. longibracteata could result from various forms of nonrandom mating, such as differences in reproductive success of a plant over time, or a sampling effect, or differences in relative recruitment of immigrant vs. resident individuals (Oddou-Muratorio ). In particular, the significant genetic differentiation in the old cohorts could be the result of extinction of their parental plants in the populations (Kitamura ), which led to the loss of alleles in the older hierarchical groups. The differences seen in young cohorts could be a result of a natural thinning process.

It is worthy of note that unlike the old growth populations FJM and FJT, which have not been disturbed, heavy logging occurred in the Nanling Mountain area (NL) where thousands of N. longibracteata trees were harvested during the period from 1960 to 1980. Nonetheless, the genetic diversity of the population HNZ in this area did not change over time (Table 2). This could be due to a high gene flow in this population, as suggested by STRUCTURE analysis (Figure 2). Furthermore, historical bottlenecks may have a long-term effect on genetic variation, to such an extent that even a recent population decline may leave no genetic signature (Taylor and Jamieson, 2008).

Potential refugia and conservation implications

The STRUCTURE results suggest that each of the sampled populations is largely genetically distinct, with very little signal of admixture, inferring that these populations have been isolated for some time, possibly in 11 separate glacial refugia. Although subtropical China has never been covered by ice sheets during the Quaternary, it has been suggested that the development of cooler and drier climates may have influenced the distribution and evolution of many plants in China and further leaded to the extinction or disappearance of some important species (Lu ; Shen ). The mountain regions including the Nanling and the adjacent Wuyi and Daiyun mountains in southern China have served as the main refugia regions for plant surviving and evolving during climatic oscillations (Ying, 2001; Wang and Ge, 2006). For example, during glacial expansion many previously dominant conifers in the northern part of East Asia were forced to migrate southwards into scattered refugia existing in small patches across mountain regions of southern China (Lu ). Ying (2001) localized three regions (the Hengduan range, the Central China and the Lingnan region) with high levels of plant diversity and endemism in China. As for N. longibracteata, the eleven possible Quaternary refugia suggested by our analysis are located in central China and the Lingnan region. Multiple refugia were also documented in other conifer and tree species in subtropical China resulting from climate changes (Lu ; Shen ; Wang and Ge, 2006; Wang ; Zou ). Our study thus provides further evidence of complex refugia in mountain regions of subtropical China during climatic oscillations.

Genetic variation is the foundation for adaptive evolution of a species, especially in forest trees that are undergoing fundamental population changes as a result of natural or anthropogenic disturbances (Jump ; Doi ; Alsos ). Consequently, for long fragmented populations, Hampe and Petit (2005) recommend preserving a maximum number of independent sites, as they are often important for the survival and evolution of biota. Considering the in situ conservation evolutionary potential of plant species, all the ten N. longibracteata populations deserve particular attention in genetic conservation programs for this species.