Genetic Diversity and Structure of Natural Quercus variabilis Population in China as Revealed by Microsatellites Markers

Abstract

Quercus variabilis is a tree species of ecological and economic value that is widely distributed in China. To effectively evaluate, use, and conserve resources, we applied 25 pairs of simple sequence repeat (SSR) primers to study its genetic diversity and genetic structure in 19 natural forest or natural secondary forest populations of Q. variabilis (a total of 879 samples). A total of 277 alleles were detected. Overall, the average expected heterozygosity (H_e) was 0.707 and average allelic richness (AR) was 7.79. Q. variabilis manifested a loss of heterozygosity, and the mean of inbreeding coefficient (F_IS) was 0.044. Less differentiation among populations was observed, and the genetic differentiation coefficient (F_ST) was 0.063. Bayesian clustering analysis indicated that the 19 studied populations could be divided into three groups based on their genetic makeup, namely, the Southwest group, Central group, and Northeastern group. The Central group, compared to the populations of the Southwest and Northeast group, showed higher genetic diversities and lower genetic differentiations. As a widely distributed species, the historical migration of Q. variabilis contributed to its genetic differentiation.

1. Introduction

Fagaceae is one of the most important components of northern sub-tropical forests and temperate forests. Quercus spp. consists of about 450 different species [1] and is the largest Fagaceae genus. Moreover, these have important ecological and economic values. Also known as Chinese cork oak, Quercus variabilis Bl. is a species of oak in the section Quercus Sect. Cerris [2] and native to a wide area of China [3], as well as the Korean peninsula and the southwestern Japanese archipelago. In China, the Qinling and Dabie mountain land areas are considered the geographical distribution centers of Q. variabilis [4]. Having one of the widest amplitudes of distribution and most complex climate types in worldwide geographic distribution [5], Q. variabilis forests have formed many stable forest ecosystems in warm temperate and northern subtropical regions in China and are playing a significant role in the conservation and improvement of water and soil [6]. It is a precious timber and economic tree species, and also considered as an important resource for natural raw materials in cork production due to its most unique feature, i.e., its bark [7,8]. In recent years, with the development of cork, tannin, and other industries in China, the phenomenon of over harvesting of Q. variabilis has become a prominent agricultural problem that has resulted in extensive destruction and waste of natural resources [9]. Currently, the main distribution areas of Q. variabilis in China are mostly natural secondary forests [4]. The original forest has become extremely rare; scattered forests are found in Guizhou, Yunnan, and Fujian (on-site investigations). The destruction of the natural populations has been strongly associated with climate change, soil erosion, and a series of ecological problems [10], hence a shortage in natural resources for Q. variabilis has caused a sharp decline in raw material production.

Genetic diversity in forests is determined by gene flow, genetic drift, selection, mutation, and other factors [11,12,13], and the levels of genetic variation and structures are the combined effects of evolutionary history, distribution area, life-styles, and breeding forms [14]. For a tree with a relatively long life cycle, genetic diversity determines its ability to adapt to changing environments, which in turn serves as the basis for maintaining long-term stability of forest ecosystems. The evaluation of genetic variability, especially for forests with wide distribution ranges and high ecological and economic values, is of great significance for forest ecosystem conservation and management of genetic resources [15]. Therefore, investigating the evolutionary history and distribution, even including the development of protection strategies and measures based on forest diversity levels and genetic structures, is imperative.

Studies on Q. variabilis have mainly concentrated on population structure and classification, biodiversity characteristics, population dynamics, and ecological function [16,17,18]. In addition, various features of reproduction, seedling cultivation, and planting techniques have also been studied [19,20,21]. Various reports on the population genetics of Quercus such as Q. acutissima [22], Q. rubra and Q. ellipsoidalis [23], Q. mongolica [24], and Q. infectoria [25] have been published in recent years. There are relatively few genetic investigations on Q. variabilis, but Xu [6] was the first to conduct a preliminary study on its genetic diversity of Chinese five natural populations by using simple sequence repeat (SSR) markers. Recent studies have focused on Q. variabilis pedigree geography such as Chen [26], who used chloroplast simple sequence repeat (cpSSR) to relatively and comprehensively analyze the historical evolution and its geographical distribution. However, this study was limited by the number of samples, and may therefore not fully represent the genetic variation of existing populations. Therefore, based on the natural resource distribution of all existing populations, the establishment of a suitable means of detection, a sufficient number of markers, and more comprehensive sample size is necessary to study the genetic structure and mating system.

China is considered as the global distribution center for Q. variabilis distribution. Furthermore, research studies on its genetic structure based on a broader sampling strategy involving its main distribution area, are of high significance to better understand the distribution pattern and formation of existing populations. Furthermore, such studies may provide important molecular bases for the protection and exploitation of the genetic resources.

2. Materials and Methods

2.1. Population Sample Information

The tender leaves of mature trees were collected from 19 different sub-locations in nine provinces (Figure 1), which is representative of the entire distribution area covered from the Henan township Neixiang, Henan (33°50′ N 111°18′ E) to Leye, Guangxi (24°47′ N 106°57′ E). Annual temperature varies from 8.3 °C to 18.4 °C, spanning two climatic zones: subtropical monsoon and Temperate monsoon. More details on the number of samples and geographical features of each group are listed in

Table 1. Geographic locations and mean values of climate parameter in 19 populations of Q. variabilis.

No.	Population	Code	Latitude (°N)	Longitude (°E)	N	Annual Temperature (°C)	Weather Patterns	Annual Precipitation (mm)
1	Ankang, Shaanxi	AK	32°40′	109°08′	45	15.3	S	803
2	Chuzhou, Anhui	CZ	32°17′	118°17′	39	15.3	S	1009
3	Macheng, Hubei	MC	31°17′	115°00′	45	16.4	S	1217
4	Xingshan, Hubei	XS	31°02′	110°07′	45	16.5	S	1498
5	Zhumadian, Henan	ZM	32°08′	114°01′	45	15.4	S	1098
6	Neixiang, Henan	NX	33°50′	111°18′	48	8.3	T	871
7	Leye, Guangxi	LY	24°47′	106°57′	45	18.4	S	1314
8	Nanjing, Jiangsu	NJ	32°03′	118°52′	45	15.6	S	1017
9	Xiaoxian, Anhui	XX	34°12′	116°56′	45	9.9	T	756
10	Pengze, Jiangxi	PZ	29°54′	116°34′	48	17.4	S	1460
11	Chengkou, Chongqing	CK	31°59′	108°40′	48	13.1	S	1165
12	Fengjie, Chongqing	FJ	31°04′	109°31′	48	13.6	S	1071
13	Pengshui, Chongqing	PS	29°12′	108°12′	48	15.6	S	1296
14	Wuxi, Chongqing	WX	31°25′	109°36′	48	12.9	S	1131
15	Jingshan, Hubei	JS	31°02′	113°07′	45	16.3	S	1035
16	Yingshan, Hubei	YS	30°45′	115°34′	48	17.1	S	1459
17	Yuanshan, Hube	YA	31°00′	111°36′	48	16.6	S	1018
18	Suizhou, Hubei	SZ	31°36′	113°18′	48	16.1	S	1007
19	Longquan, Zhejiang	LQ	28°02′	119°05′	48	17.6	S	1793
ALL					879

Population: sub-location, province; N: sample size; S: Subtropical monsoon; T: Temperate monsoon.

. The studied forests are typical natural secondary forests, and the plants of each population were separated by at least 50 m. The collected samples were brought back to the lab, dried by silica, and stored at room temperature until analysis.

2.2. Experimental Methods

Total DNA was extracted from the samples by using a modified cetyl trimethylammonium bromide (CTAB) lysis-silica beads adsorption method [27]. The SSR primers used in population identification are shown in

Table 2. Source, repeat motif, and polymorphism information for the 25 microsatellite loci analyzed in the 879 Q. variabilis trees.

Locus	Repeat Motif	Na	Ho	He	PIC	Source
2p24	(CA)14	8	0.429	0.717	0.794	Alexis, R.S. et al. [29]
E71-72	(GA)46	8	0.212	0.771	0.803	Qin, Y.Y. et al. [30]
PIE040	(TTC)8	10	0.157	0.724	0.747	Alexis, R.S. et al. [29]
GOT040	(GA)11	5	0.389	0.468	0.782	Durand, J. et al. [31]
G0T009	(TC)7	6	0.208	0.590	0.673	Durand, J. et al. [31]
FIR053	(GTG)7	8	0.294	0.751	0.786	Durand, J. et al. [31]
FIR039	(CT)7	9	0.459	0.754	0.782	Durand, J. et al. [31]
FIR004	(CT)18	8	0.472	0.760	0.778	Alexis, R.S. et al. [29]
G11	(TC)22	4	0.324	0.551	0.615	Xu, X.L. et al. [6]
PL111-112	(TC)9	6	0.534	0.694	0.720	Qin, Y.Y. et al. [30]
PL229-230	(AG)15	9	0.387	0.669	0.689	Qin, Y.Y. et al. [30]
VIT107	(TA)13	5	0.306	0.452	0.524	Durand, J. et al. [31]
DN949726	(GAT)6	15	0.384	0.861	0.878	Saneyoshi, U. et al. [24]
E11-12	(GA)32	14	0.578	0.851	0.887	Qin, Y.Y. et al. [30]
E79-80	(TC)18	24	0.626	0.841	0.889	Qin, Y.Y. et al. [30]
EE812	(AG)7	20	0.393	0.804	0.856	Zhang, Y.Y. et al. [22]
G7	(TC)17	21	0.447	0.784	0.883	Xu, X.L. et al. [6]
G16	(AG)21	20	0.606	0.829	0.860	Xu, X.L. et al. [6]
PL127-128	(AG)12	18	0.408	0.846	0.874	Qin, Y.Y. et al. [30]
Q16	(GA)18	26	0.704	0.875	0.911	Xu, X.L. et al. [6]
EE802	(CT)8	7	0.427	0.748	0.790	Zhang, Y.Y. et al. [22]
EE856	(GGT)6	4	0.308	0.418	0.423	Zhang, Y.Y. et al. [22]
FIR048	(CT)9	8	0.438	0.758	0.784	Durand, J. et al. [31]
FIR110	(TC)20	6	0.186	0.552	0.602	Alexis, R.S. et al. [29]
PIE125	(GGAAGC)3	8	0.353	0.619	0.664	Durand, J. et al. [31]
mean		11	0.401	0.707	0.760
min		4	0.157	0.418	0.423
max		26	0.704	0.875	0.911

N_a: number of alleles, H_o: observed, H_e: expected heterozygosities, PIC: polymorphism information content.

. The total volum e of each single locus polymerase chain reaction (PCR) was 20 μL, which consisted of 20 ng of DNA, 50 mM KCl, 20 mM Tris-HCl (pH 8.3), 2 mM MgCl₂, 0.2 mM dNTP, 0.5 μM of each primer (forward and reverse), and 1 U of Taq enzyme (Qiagen Inc., Hilden, Germany) PCR thermal cycling was performed on a GeneAmp PCR System 9700 (Applied Biosystems Inc., Carlsbad, CA, USA) using the following conditions: an initial denaturation step of 94 °C for 4 min; followed by 30 cycles of denaturation at 94 °C for l min, annealing at N °C for 2 min (N is the specific annealing temperature for each pair of primers), and extension at 72 °C for 1 min; and a final extension at 72 °C for 8 min and final storage at 4 °C. The PCR products were detected by using capillary electrophoresis (Bioptic100-automated nucleic acid analyzer, BiOptic Inc., Tucson, AZ, USA). Electrophoresis results were interpreted using Q-editor software [28].

2.3. Data Analysis

Genetic parameters, calculated by the Fstat software, were: average number of alleles: N_a; observed heterozygosity: H_o; expected heterozygosity or gene diversity, H_e [32]; allelic richness, AR [33]; the fixation index (F_IS: the inbreeding coefficient) of every point and group level; the differentiation index between pairwise populations (F_ST, G_ST, R_ST) and the matrix of F_ST between various groups on which F-statistics were based [34,35]; Hardy-Weinberg equilibrium testing and estimation of allele frequency were calculated [36]. The polymorphic information index (PIC) and Linkage disequilibrium were calculated by using the PowerMarker v3.23 software (Bioinformatics Research Center; North Carolina State University, Raleigh, NC, USA) [37]. P values were adjusted by using a Bonferroni correction. Private alleles were calculated by using the GenAlEx software [38,39]. Gene flow (N_m) was calculated using the formula: N_m = (1 − F_ST)/4F_ST [40]. To investigate genetic diversity and geography-space-climate relationships [41], the average AR value was calculated using a total of 25 loci, whereas the relationships between H_e value and each geographic population (latitude and longitude), annual average temperature, and annual precipitation in climate parameters (obtained through DIVA-GIS software [42]) were calculated by using the Kendall rank-correlation test as provided in the SAS v9.1 software (SAS Institute 2001, Wallisellen, Switzerland). Detection of bottleneck effects was calculated by using BOTTLENECK version 1.2 software [43]. A two-phase mutation model (TPM) and a stepwise mutation model (SMM) were used for Wilcoxon signed-rank tests. The parameter settings included a 90% SMM and 10% TPM with a variance of 12%, and 1000 repeats [44].

For genetic relationships at the population level, a factorial correspondence analysis (FCA) was performed by using the Genetix 4.05 software (CNRSUMR 5000; Universite Montpellier II, Montpellier, France) [45]. Population genetic structure was analyzed based on Bayesian clustering using the STRUCTURE software [46]. For clustering from K = 1 to K = 20 (populations + 1), an admixture ancestry model and correlated allele frequency model were used to perform a Markov chain Monte Carlo simulation algorithm (MCMC) [47]. The length of the burn-in period at start time was set as 100,000; MCMC after the length of the burn-in period was set as 100,000, and for each of K value, simulation calculation was repeated 10 times. The method from Evanno was used to determine the optimal K value [48]. The relative proportions of the geographical distribution diagram were plotted by using the Arc-GIS software (Environmental Systems Research Institute, Inc., Redlands, CA, USA). Molecular variance analysis (AMOVA) for population genetic structures [49] and a Mantel test to detect the presence of geographic segregation between populations by means of logarithmic normalization of geospatial distance and genetic distance [50] were calculated by using the Arlequin software (CMPG, Institute of Ecology and Evolution; University of Berne, Berne, Switzerland).

3. Results

3.1. Genetic Diversity

A total of 277 alleles were amplified from 25 SSR primers. On average, each microsatellite loci could be amplified 11 times, and the number of alleles ranged from a minimum of 4 (G11, EE856) to a maximum of 26 (Q16). The H_o of each microsatellite loci ranged from 0.157 (PIE040) to 0.704 (Q16), and the average observed heterozygosity was 0.401. For each microsatellite loci, H_e was higher than H_o, and ranged from 0.418 (EE856) to 0.875 (Q16), and the average H_e was 0.707. AR ranged from 3.50 (EE856) to 21.12 (Q16), and the average AR was 9.51. The PIC of SSR microsatellite loci ranged from 0.423 to 0.911, and the average PIC was 0.760 (Table 2). After detecting significance (P < 0.05), we did not find the existence of linkage disequilibrium between any two points.

Based on statistics of population genetic diversities of 25 microsatellite loci from 19 geographic populations (

Table 3. Mean values of genetic diversity statistics for 25 microsatellite loci in 19 Q. variabilis populations.

Code	Na	AR	He	FIS	TPM	SMM	FST	GST	RST
AK	8.20	8.03	0.710	0.038	ns	ns
CZ	7.12	7.12	0.683	0.054	ns	0.037 *
MC	8.32	8.16	0.720	0.033	ns	ns
XS	8.68	8.56	0.723	0.045	ns	ns
ZM	7.84	7.65	0.711	0.043	ns	ns
NX	7.92	7.75	0.709	0.042	ns	ns
LY	8.44	8.23	0.707	0.046	ns	ns
NJ	7.96	7.81	0.725	0.049	ns	ns
XX	7.04	6.90	0.690	0.041	ns	ns
PZ	8.76	8.54	0.725	0.042	ns	ns
CK	8.48	7.54	0.745	0.044	ns	ns
FJ	7.12	6.95	0.623	0.045	0.045 *	0.002 *
PS	8.68	7.73	0.699	0.050	ns	0.019 *
WX	7.44	7.29	0.701	0.055	ns	ns
JS	7.92	8.45	0.705	0.039	ns	ns
YS	8.24	8.00	0.694	0.049	ns	ns
YA	8.04	7.86	0.726	0.041	ns	ns
SZ	8.16	7.93	0.728	0.043	ns	ns
LQ	7.76	7.54	0.716	0.039	ns	ns
mean	8.01	7.79	0.707	0.044			0.063	0.060	0.073
min	7.04	6.90	0.623	0.033
max	8.76	8.56	0.745	0.055

N_a: number of alleles; AR: allelic richness; H_e: expected heterozygosity; F_IS: fixation index. TPM: a two-phase mutation model, SMM: a stepwise mutation model; * Significant deviation from Wilcoxon signed-rank tests (P < 0.05). F_ST, G_ST differentiation among populations according to Weir and Cockerham [34]; R_ST, Slatkin [51].

), the average number of alleles in each population was 8.01; the population with the highest number of alleles was PZ (8.76), and the lowest was XX (7.04). The average number of effective alleles was 4.52; the highest was LY, and the lowest was FJ. The overall average AR was 7.79; the highest one (8.56) was distributed in XS, whereas the lowest was in XX (6.90). Q. variabilis manifested a significant loss of heterozygosity, and the mean of inbreeding coefficient was 0.044. The average H_e was 0.707; the highest one was 0.745 with a distribution in CK, whereas the lowest one was 0.623 in FJ. BOTTLENECK analysis showed that in the 19 analyzed populations, three groups could be identified: CZ (SMM, P < 0.05), PS (SMM, P < 0.05), and FJ (TPM, SMM, P < 0.05), which may have recently undergone a severe decline in population size, i.e., the bottleneck effect had occurred.

Kendall rank-correlation analysis of genetic diversity of geographic populations with latitude/longitude and climatic parameters showed that in general, the AR of Q. variabilis geographic populations has a significant positive correlation with higher average annual temperatures (r = 0.507, P < 0.01), whereas no correlations among H_e, geographic latitude, and climate were observed.

3.2. Genetic Differentiation and Genetic Structure

The overall population differentiation degree among geographical populations of Q. variabilis: F_ST = 0.063, G_ST = 0.060 (P < 0.001) have been reported in Table 3. Gene flow: N_m = 3.648. AMOVA analysis [52] showed that 6.3% (P < 0.001) of the genetic variations was among populations. The greatest percentage of variation was contained within populations. The F_ST matrix between every two populations is shown in

Table 4. The pairwise F_ST for all populations of Q. variabilis.

Code	AK	CZ	MC	XS	ZM	NX	LY	NJ	XX	PZ	CK	FJ	PS	WX	JS	YS	YA	SZ	LQ
AK	0.000
CZ	0.051	0.000
MC	0.067	0.079	0.000
XS	0.040	0.042	0.049	0.000
ZM	0.059	0.044	0.062	0.048	0.000
NX	0.062	0.033	0.078	0.064	0.044	0.000
LY	0.055	0.040	0.073	0.040	0.058	0.061	0.000
NJ	0.063	0.048	0.051	0.035	0.050	0.053	0.050	0.000
XX	0.076	0.068	0.111	0.072	0.077	0.061	0.082	0.075	0.000
PZ	0.059	0.055	0.047	0.047	0.051	0.076	0.050	0.048	0.099	0.000
CK	0.027	0.050	0.049	0.040	0.050	0.057	0.047	0.041	0.077	0.039	0.000
FJ	0.077	0.081	0.105	0.078	0.108	0.112	0.059	0.091	0.135	0.082	0.079	0.000
PS	0.039	0.056	0.069	0.040	0.064	0.074	0.047	0.057	0.081	0.057	0.050	0.068	0.000
WX	0.078	0.082	0.101	0.066	0.091	0.096	0.053	0.078	0.101	0.076	0.079	0.082	0.053	0.000
JS	0.058	0.076	0.067	0.071	0.086	0.093	0.084	0.065	0.113	0.049	0.045	0.093	0.079	0.111	0.000
YS	0.040	0.040	0.049	0.039	0.041	0.058	0.043	0.035	0.088	0.029	0.033	0.065	0.036	0.070	0.060	0.000
YA	0.048	0.047	0.068	0.054	0.063	0.061	0.054	0.050	0.082	0.057	0.036	0.082	0.066	0.099	0.064	0.054	0.000
SZ	0.052	0.079	0.056	0.052	0.077	0.089	0.062	0.051	0.105	0.046	0.028	0.085	0.062	0.083	0.054	0.053	0.049	0.000
LQ	0.052	0.050	0.050	0.031	0.052	0.069	0.049	0.043	0.082	0.039	0.037	0.085	0.054	0.084	0.060	0.038	0.045	0.053	0.000

See Table 1 for population information.

; the differentiation coefficient between populations CK and SZ was the smallest, (F_ST = 0.028); that between XX and FJ was the largest (F_ST = 0.135).

A FCA on genetic structure differences of pattern detection at the population level showed that the populations could be divided into three groups based on geographic location: (a) the Southwest populations: LY, PS, FJ, and WX; (b) the Central populations (the most resource-rich center of Q. variabilis distribution in China): AK, MC, XS, PZ, CK, JS, YS, YA, SZ, and LQ; and (c) the Northeast populations: NJ, CZ, ZM, NX, and XX. The Central populations were adjacent to the Southwest populations and Northeast populations in the FCA clustering groups, whereas the Southwest populations were distantly located from the Northeast populations in the FCA clustering groups (see Figure 2a,c). Small genetic differences between the Central populations and the other two populations were observed, whereas larger genetic differences were detected between the Southwest populations and the Northeast populations.

For a more detailed understanding of the genetic structures, using a ΔK value that Evanno proposed to determine a reasonable K, the ΔK value reached a maximum when K = 3 (Figure 3). The 19 geographic populations of Q. variabilis could then be divided into three genetically distinct groups (Figure 4). In addition, we used the relative proportions of the geographical distribution diagram to visually compare the distribution ratio of each group within the population (Figure 5). The 19 populations as a whole were divided into three groups based on geographic distribution (Figure 4 and Figure 5); the first group was the Southwest group (pop1–4), the second group was the Central group (pop5–14), and the third group was the Northeast group (pop15–19). Average distribution ratios [53] were 76%, 72% and 73%, respectively. The Central group had the lowest distribution ratio, i.e., the highest mixing degree among groups. XX (Northeast group, ratio: 87.9%), JS (Central group, ratio: 83.4%), and WX (Southwest group, ratio: 88.9%) were the source populations in each group, which had the highest distribution ratio.

The results of a FCA and STRUCTUE analysis indicated that the population as a whole could be divided into three genetic groups. Furthermore, AMOVA analysis showed 1.75% variation between groups (P < 0.001), 5.59% variation within groups (P < 0.001), and 92.66% variation within populations (P < 0.001). All genetic groups were ordered according to the calculated F_ST value: Northeast group (F_ST = 0.070) > Southwest group (F_ST = 0.039) > Central Group (F_ST = 0.027). The standard genetic distance and geographic distance of F_ST/(1 − F_ST) were calculated using the Mantel test, which revealed that the changes in the F_ST values did not follow the geographical distance separating mode (R = 0.0013, P = 0.371). The geographical distances between populations AK, CK and FJ, WX were very small, but FCA and Bayesian analyses did not cluster those into the same groups, and AMOVA analysis showed that most of the variations that existed within groups resulted in significant genetic structural differences in the populations AK, CK and FJ, WX.

4. Discussion

4.1. Genetic Diversity

The level of species genetic diversity is often associated with specific characteristics such as the length of life cycle, mating system and reproduction, the size of the geographic range, and genetic exchange [11,54,55,56]. The H_e (0.707) of Q. variabilis in our study was slightly lower than that of Xu’s previous study (0.806) [6]. The main reason for this difference is that only five Chinese populations of Q. variabilis were studied in Xu’s study, whereas we examined representative populations within the distribution range from the South to the North. Therefore, the population genetic variation of the species reflected in the present study was more comprehensive. Compared to the overall genetic variation level of Q. variabilis with other Fagaceae Quercus genera species such as the summer oak (Q. robur) (H_e = 0.764) [29] and Liaodong oak (Q. liaotungensis) (H_e = 0.754) [30], the levels of diversity were about the same yet slightly higher than that of the Sawtooth oak (Q. acutissima) (H_e = 0.660) [22] and Mongolian oak (Q. mongolica) (H_e = 0.630) [24]. Compared to the other endangered species or narrow domains that were left over from the Quaternary Ice Age in China such as Ginkgo biloba (H_e = 0.241) [57], the overall genetic diversity level of Q. variabilis was high, which was mainly due to the extensive distribution of the species distributed from 19° N to 42° N across the temperate forests. Because of huge differences in climate and habitat conditions among populations in the distribution area, Q. variabilis may have differentiated adapted ecological and genetic types, thereby resulting in a wide range of genetic variations. In addition, affected by the fact that the populations were sampled from natural secondary forest, heterozygote deficits were shown (F_IS = 0.044).

H_e and AR of each population by Kendall rank-correlation analysis revealed that neither followed the geographical space separation mode. However, changes in AR clearly showed significant correlations with changes in the average annual temperature, indicating that temperature has a significant impact on Q. variabilis. Among the 19 populations studied, the genetic diversity of the central populations (AR = 8.06, H_e = 0.718) was higher than that of the Southwest populations (AR = 7.55, H_e = 0.683) and the Northeast populations (AR = 7.45, H_e = 0.704). The two populations with the highest genetic diversity—XS (AR = 8.56) and CK (H_e = 0.745)—both belonged to the Central group; these results supported the conclusions of Chen [26] in the populations of Q. variabilis in the Central region, which experiences an average annual temperature of 15 °C or above and has a higher level of genetic diversities than other areas. Bottleneck detection revealed that the Southwest populations FJ and PS, and Northeast population CZ experienced a genetic bottleneck, which may be the main reason for the low level of genetic diversity in all three populations [14,43,58].

Mayr proposed a “core–periphery” hypothesis that states that compared to the core groups in general [59], the genetic diversity of peripheral groups decreased and genetic differentiation increased under the pressures of the bottleneck effect, genetic drift, and selection pressures [60]. Based on the genetic diversity of five Q. variabilis populations, Xu concluded that the genetic diversity of the central groups was higher than that of the peripheral groups [6], but because of the limitations of the sample population scales, these findings could only be used as a reference. In our study, the samples of 19 populations were collected from a representative distribution in the main distribution areas. For the populations located in the geographical centers of existing distribution areas (the Central group), overall genetic diversities were higher than those in the peripheral regions (the Southwest group and the Northeast group). In addition, private alleles were only found in the southwestern population, such as LY; this indicated that there was genetic variation in the peripheral population in adaptation to the climate and environment. Accordingly, the diversity distribution of Q. variabilis was in line with the “core–periphery” hypothesis.

4.2. Genetic Differentiation

The genetic differentiation level of Q. variabilis was moderately low (F_ST = 0.063), which was in agreement with the results of AMOVA analysis (among populations variation = 6.3%). The genetic differentiation level of Q. variabilis is similar to that of the summer oak (Q. robur) (F_ST = 0.080) [29] and Mongolian oak (Q. mongolica) (F_ST = 0.077) [24]. Its genetic differentiation level is close to that reported by Hamrick [61], who investigated diversity in widespread species, including long-lived woody perennials and wind-pollinated outcrossing species. As a perennial tree, Q. variabilis is a typical wind-pollinated and outcrossing species, which undoubtedly increases gene flow between populations and reduces differentiation among populations.

Gene flow determines the genetic structure and survival potential of future populations of a species [12]. The present study showed that Q. variabilis populations are characterized by a relatively large gene flow (N_m = 3.648) that may inhibit genetic drifting and prevent population genetic differentiation [11]. The large gene flow is mainly determined by its biological means of spreading via pollen and seeds. Xu [6] reported a higher gene flow (N_m = 5.239) than this study, which objectively reflected a decreasing gene flow for Q. variabilis in the last decade; the original natural forest of Q. variabilis has degenerated into a natural secondary forest, which may have caused a reduction in gene flow.

The clustering results of STRUCTURE Bayesian clustering showed that the 19 populations could be divided into three groups: the Southwest group, the Central group, and the Northeast group. These findings were in agreement with the results of a FCA, wherein the mix in the population genetic structure of the Southwest group and the Northeast group was relatively low, and there was a large difference in the genetic composition between the two groups. On the other hand, the Central populations were the largest; they showed a high mixing degree and the lowest differentiation. In general, the population that is capable of maintaining its level of genetic diversity is often the largest population group [62]. Therefore, as expected, the Central group showed the highest capacity to maintain its level of genetic diversity. The Mantel test showed that correlations of genetic and geographic distances between the groups were not significant (R = 0.0013, P = 0.371), which showed that geographical distance is not the main reason for genetic differentiation among populations.

As a widely distributed species and perennial tree, the genetic structure formation of Q. variabilis should not be analyzed from a single geographic distribution and distance. During the Quaternary Ice Age, which was affected by the global glacial climate, Q. variabilis in northern China was unable to adapt and eventually disappeared. Furthermore, the South of Qinling and Dabie Mountains had become barriers that blocked the Quaternary glaciers, with their complex landforms and diverse ecological environments rendering them sources of plant heterogeneous differentiations [63]. The area thus served as a sanctuary for Q. variabilis during Quaternary glaciation [64] and eventually became the geographical distribution center for Q. variabilis [4]. At the end of the Ice Age and with the recovery of temperature as well as other natural conditions, it is assumed that some populations in the distribution center differentiated into varying types which adapted to a cold, dry climate and began to migrate to high latitudes in the northeast, while some others adapted to warm and moist types and began to spread southwestward, and gradually formed the geographical distribution pattern of Q. variabilis (in the present study, there were significant correlations between AR and the average annual temperature). This evolutionary geographic history has been verified by simple haplotypes in peripheral populations, and populations in the central region represented almost all haplotypes of Q. variabilis [26]. It is in this complex historical dynamic migration process that the genetic structure of the existing distribution of Q. variabilis was formed.

In addition, human activities may also have an impact on the genetic structure of various species [65,66,67]. However, the results of the present study showed that among the 19 Q. variabilis populations, the population in the center of the distribution did not experience a bottleneck effect, whereas three populations at the periphery (FJ, PS, CZ) did, and three genetic lineages (cluster I,II and III) obtained from STRUCTURE showed distributions in every population (Figure 4 and Figure 5). Based on these results, we inferred that Q. variabilis has been a widely distributed species for several historical periods, and the history of its population dynamics and genetic variations stretches back further than the time of human activities [68], so human activity is not the major factor that has affected its current genetic structure.

5. Conclusions

This represents the first research carried out on the genetic diversity and structure of Q. variabilis based on a broader sampling strategy involving its main distribution area. Q. variabilis, in general, was in line with the “core–periphery” hypothesis. The overall genetic diversity level of Q. variabilis was relatively high; the genetic diversity level of populations in the geographical distribution center was higher than that of peripheral populations. The accessions from southwest, center and northeast areas clustered into three separate groups, and the genetic structure of Q. variabilis was mainly affected by the preferable adaptability to the climate and environment in the complex historical dynamic migration process. Due to the development of cork, tannin, and other industries, the distribution scale of Q. variabilis has gradually reduced; we should establish conservation measures to prevent population decline; peripheral populations FJ, PS and CZ, in particular, suffered the bottleneck effect. Additionally central populations with abundant genetic variation could be used as the preferred germplasm resource in plantations. Furthermore, private alleles were found only in peripheral population LY; they are of great research value to the maintenance and evolution of species in disadvantageous habitats.