Population structure and genetic diversity of Tamarix chinensis as revealed with microsatellite markers in two estuarine flats

Sample collection and DNA extraction

Samples were taken from nine natural populations of T. chinensis distributed from Beijing, Shandong Province and Zhejiang Province. Sampling mainly focused on two estuarine flats of different climate zones: Yellow River Delta of the warm-temperate zone, and Hangzhou Bay of the subtropical zone (Fig. 1, Table 1). These two estuarine flats have many T. chinensis populations of different sizes and ages and populations of Yellow River Delta usually have larger size than those of Hangzhou Bay. Four populations (YHK, YDG, YXX and YHD) from Yellow River Delta and three populations (HHJ, HLS and HCX) from Hangzhou Bay were sampled. One population (CY) of the Bohai Sea coast from Changyi City and one inland population (FS) from Beijing were also included. A total of 19–29 individuals per population were collected according to the population size and a total of 204 individuals from nine populations was sampled. We collected young leaves from individual trees spaced at least 50 m apart, in order to avoid sampling related individuals. Fresh leaves were desiccated insilica gel. Geographical coordinates and altitude were acquired with GPS (Table 1). A modified cetyltrimethyl ammonium bromide (CTAB) method (Su et al., 1998) was used to extract total genomic DNA. Quality and quantity of the extracted DNA were measured with 1.0% (w/v) agarose gel electrophoresis and a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), respectively. DNA was diluted to 20 ng/µl and stored at −20 °C for later use.

We used a set of 12 highly polymorphic SSR markers (Table 2) for population genetic analysis, six of which were EST-SSRs (Zhang et al., 2011) and the remaining six of which were genomic SSRs (Zhang, Wen & Xu, 2019). The polymerase chain reaction (PCR) was carried out according to (Zhang, Wen & Xu, 2019). Forward primers were labeled with four color fluorescent dyes (6-FAM, HEX, TAMRA and ROX) (Invitrogen, Shanghai, China).The PCR conditions of both kinds of SSR were as follows: 4 min denaturation at 94 °C; followed by 30 cycles of 94 °C denaturation for 30s, varied annealing temperatures ranging from 49.6–58.3 ° C for 30s (Table S1), and an elongation at 72 °C for 30s; and a final extension step at 72 °C for 5 min. PCR products were run on an ABI3730 XL Analyzer (Applied Biosystems, Foster City, CA, USA) and allele sizes were determined with GeneMapper 4.0 software (Applied Biosystems).

Data analysis

The occurrence of null alleles has a biased effect on the estimation of population genetic parameters, so we used MICROCHECKER v2.2.3 (Van Oosterhout et al., 2004) to check for the presence of null alleles for all loci. Hardy-Weinberg equilibrium exact tests and multilocus linkage disequilibrium tests were performed with GENEPOP (version 1.2) (Raymond & Rousset, 1995), and sequential Bonferroni correction tests were applied to adjust the significance levels of genotypic linkage disequilibrium. GenAlEx 6.502 (Peakall & Smouse, 2005) was used to calculate genetic diversity indices including observed heterozygosity (Ho), expected heterozygosity (He), number of allele (Na), number of effective alleles (Ne), inbreeding coefficient (F_IS), and Shannon diversity index (I). The genetic differentiation coefficient (Fst) across all populations at each locus and over all loci, gene flow (Nm = (1 − F_ST)/4 F_ST), and allelic richness (Ar) were calculated using FSTAT version2.9.3 software (Goudet, 2001). We used CERVUS version 3.0.3 (Kalinowski, Taper & Marshall, 2007) to calculate the polymorphism information content (PIC) value of each locus according to: $PIC=1-{\sum }_{i=1}^n{f}_i^2$ where f_i is the frequency of the ith allele and n is the allele number.

Analysis of molecular variance (AMOVA) was performed using ARLEQUIN 3.5 (Excoffier & Lischer, 2010) to partition the genetic variability into within- and among-population components. The significance was tested based on 1,000 permutations for the AMOVA results (Excoffier, Laval & Schneider, 2005). Bottleneck events of the studied populations were tested using BOTTLENECK 1.2.02 (Piry, Luikart & Cornuet, 1999) with stepwise mutation model (SMM) and two-phased mutation model (TPM). Recently bottlenecked populations have higher observed heterozygosity (unbiased gene diversity) than expected from the observed number of alleles across microsatellite loci. Variance for TPM was set to 30 and the proportion of SMM in TPM was set to 80%. The “mode-shift” of allele frequency distribution was also used to identify bottlenecked populations (Luikart et al., 1998).

Bayesian clustering analysis was performed using STRUCTURE v2.3.4 software (Pritchard, Stephens & Donnelly, 2000) to determine the number of genetically homogeneous groups of the collected individuals and to evaluate the amount of admixture between individuals with the admixture model and allele frequencies correlated. We run the program with 100,000 burn-in steps followed by 100,000 Markov chain Monte Carlo (MCMC) steps, and 10 independent runs. The optimal K was determined through the online program STRUCTURE HARVESTER (Earl & Von Holdt, 2012) for each K from 1 to 10. Graphics were displayed with the DISTRUCT program (Rosenberg, 2004).We also implemented the Discriminant Analysis of Principal Components (DAPC) (Jombart, Devillard & Balloux, 2010) using the “adegenet” package (Jombart, 2008) for software R version 3.4.2 (R Core Team, 2017) to evaluate population structure. DAPC is a multivariate approach that is free of assumptions of Hardy–Weinberg equilibrium (HWE) and maximizes the genetic differentiation between groups with unknown prior clusters, thus improving the discrimination of populations. A dendrogram was also constructed using the unweighted pair-group method with arithmetic means (UPGMA) algorithm based on Nei’s genetic distance among populations with the NTSYS version 2.1 software (Rohlf, 2000).

We used the mantel function in the “vegan” R package (Dixon, 2003) to test the correlation between genetic distance (pairwise F_ST) and Euclidean geographical or environmental distance. The climatic layers of 19 bioclimatic variables under current periods were downloaded from the WorldClim database (http://www.worldclim.org/) with a resolution of 2.5 arc-minutes to evaluate the effects of environment on genetic variation. The Vifstep function of the R package usdm (Naimi et al., 2014) was used to select the geographical and environmental variables with variance inflation factor (VIF) below 10 for Euclidean distance calculation. VIF measures the collinearity of the variables and the larger the VIF, the stronger the linear relationship between two or more predictor variables. After forward selection, five bioclimate variables (Bio2, Bio5, Bio8, Bio15 and Bio18) were chosen to constitute the environmental variable matrix. Euclidean distances of the five chosen environmental variables from different populations constituted the environmental variable matrix for the Mantel test. The significance of the correlation tests was evaluated with 9,999 permutations. To determine whether geographical or environmental variables may have affected genetic distance, we also conducted a partial Mantel test using the mantel.partial function in the R package vegan with 9,999 permutations.

Genetic diversity

The micro-checker analysis revealed one locus in population FS and population HHJ, three loci in population HCX, four loci in population YHD, five loci in population YHK with null alleles, while no null alleles were detected in four populations (CY, YDG, YXX, and HLS) (Table S2). Specifically, loci 2, 4, 5, 9, 11 revealed null alleles in two populations; loci 3, 8 and 10 in one population, while loci 1, 6, 7, 12 showed no null alleles. Hardy-Weinberg exact tests showed two loci in one populations (YXX), three loci in four populations (CY, FS, YHD, and HCX), four loci in two populations (YDG and HHJ), five loci in one population (YHK), and seven loci in one population (HLS) that deviated from equilibrium at p < 0.05 or at p < 0.01 (Table S3A). At population level, only one population (CY) was in HWE, and other populations deviated from equilibrium at p < 0.05 or p < 0.01 (Table S3B). Except for populations YDG, YXX and YHD, genotypic linkage disequilibrium (LD) was detected with a range of one locus pair (for population HYD) to 16 locus pairs (for population HHJ) after Bonferroni correction (Table S4A). Across all populations, 12 out of 132 locus pairs showed linkage disequilibrium, with loci 2 and 4 having the lowest number of pairs (1 pair) in LD (Table S4B).

We used 12 pairs of SSR primers to generate 162 alleles from 204 individuals of nine T. chinensis populations (Table 2). The results showed an average value of 13.50 for all alleles ranging from three (Gssr4) to 33 (AB298390) per locus. The effective number of alleles (Ne) ranged from 2.312 (Gssr4) to 10.636 (AB298390), with an average of 4.607. The observed heterozygosity (Ho) for each locus varied from 0.3939 (Gssr4) to 0.8469, with an average of 0.6533 and the expected heterozygosity (He) varied from 0.4073 (Gssr4) to 0.9083 (AB298390), with an average of 0.7111. Shannon’s information index (I) ranged from 0.5401 (Gssr4) to 2.7643 (AB298390), with an average of 1.569. The genetic differentiation coefficient (Fst) (0.036 (Gssr5)–0.110 (Gssr3)) and the polymorphic information content (PIC) values (0.406 (Gssr4) −0.899 (AB298390)) were detected at all loci. Gene flow (Nm) ranged from 2.023 (Gssr3) to 6.694 (Gssr5), with an average of 4.245.

The genetic diversity analysis of nine T. chinensis population of was shown in Table 3. Allele richness ranged from 8.56 (HHJ) to 10.38 (YXX) with an average of 9.37. The number of alleles (Na) varied from 7.25 for the HHJ population to 9.5 for the YXX population with an average of 8.042 and the effective number of alleles (Ne) varied from 3.978 for the HHJ population to 5.915 for the YXX population with an average of 4.811. Shannon’s information index (I) varied from 1.547 (HHJ) to 1.854 (YXX), with a mean of 1.659. Mean value of expected heterozygosity (He = 0.696) was slightly higher than that of observed heterozygosity (Ho =0.674). The fixation index (Fis) ranged from −0.099 for the HHJ population to 0.177 for the HHK population, with an average of 0.025. Out of the nine studied populations, population YXX from Yellow River Delta showed the highest values of genetic diversity parameters except for the number of private allele (Np =7), which was only lower than that of population HCX (Np =16). Altogether, 54 private alleles were identified and populations YHK and CY had the lowest number (Np =2) while population HCX displayed most private alleles (Np =16).

Averagely, higher values of genetic diversity indices (Na =8.5625, Ne =5.294, I = 1.73125, Ho =0.680, He =0.710) for populations from Yellow River Delta were obtained than those of populations from Hang Zhou Bay (Na =7.7916, Ne =4.528, I = 1.631, Ho =0.666, He =0.695) except for the number of private allele (Np =4.75 for populations from Yellow River Delta and Np =9 for populations from Hangzhou Bay) (Table S8). Inbreeding coefficients of populations from Yellow River Delta (−0.072−0.017, average =0.028) are lower than those of populations from Hangzhou Bay (−0.099−0.148, average = 0.0367).

Under SMM and TPM, our Wilcoxon tests detected no significant heterozygote excess for the nine populations (Table S7). Likewise, mode-shift tests revealed an L-shaped distribution of alleles, suggesting the absence of a recent bottleneck.

Population differentiation and genetic structure

Low population genetic differentiation ( F_ST = 0.0532, p < 0.001) (Table 4) was detected across T. chinensis populations. F_ST values of the 12 loci ranged from 0.036(for Gssr5) to 0.110(for Gssr3) (Table 2). AMOVA test revealed significant differentiation detected among regions (Fct =0.0252, p < 0.05), among populations within regions (Fsc =0.028, p < 0.001), and 94.68% of variation was within-populations. The pairwise F_ST values between populations are shown in Table S5 . The genetic differentiation between the CY and HCX populations was highest at 0.110.

Bayesian-based structure analysis supported the populations studied were classified into two clusters almost without admixture with a maximum delta value of 72.40 at K = 2 (Table S6). Populations from Yellow River Delta, CY population, FS population from Beijing, and one population from Hangzhou Bay (HHJ) formed the first cluster, and the two other populations, Hangzhou Bay HCX and HLS, formed the second cluster (Fig. 2A). At K = 3, HHJ population of the first cluster was assigned to the third cluster only with 2 admixed individuals (Fig. 2B), which is consistent with the DAPC results. At K = 4, HHJ population and two other Hangzhou Bay populations (HLS and HCX) formed two distinct clusters almost without admixture, while the individuals from the remaining populations were mostly admixed with two other clusters with approximate proportion (Fig. 2C).

The nine T. chinensis populations were grouped into two groups based on UPGMA (Fig. 3), two populations from Zhejiang Province, and the remaining populations. HHJ population was grouped into the second group with populations from Shandong Province and Beijing, but with the shortest branch length to populations of the first group. DAPC identified three distinct groups (Fig. 4), which is consistent with structure result at K = 3.

Mantel test

Mantel tests showed the studied populations conformed to the patterns of isolation by distance (IBD) (R² = 0.5784, p = 0.012) (Fig. 5A) and isolation by environment (IBE) (R² = 0.5923, p = 0.015) (Fig. 5B) at significant levels. Given the autocorrelation of geographical distance and environmental distance (R² = 0.9357, p = 0.002) (Fig. 5C), partial mantel tests were performed, which detected highly significant IBD ( r = 0.7324, p = 0.002), but no IBE patterns in the populations studied.

High genetic diversity of T. chinensis

With highly polymorphic microsatellite markers, we detected high genetic diversity of T. chinensis (12 nSSRs, He =0.7111) (Table 2) comparable to the widely distributed long-lived tree species including Quercus acutissima (12 nSSRs, He = 0.760) (Zhang et al., 2013), Q. liaotungensis (19 nSSRs, He =0.801) (Guo et al., 2021), Q. rubra (15 SSRs, He =0.78) (Lind & Gailing, 2014), and Q. ellipsoidalis (15 SSRs, He =0.76) (Lind & Gailing, 2014). Similarly, population genetic analysis through dominant markers including RAPD and ISSR also detected high genetic diversity in this species (Nei index =0.4061 for RAPD analysis and Nei index =0.276 for ISSR analysis) (Jiang, Chen & Bao, 2012; Zhao et al., 2008). The high level of genetic diversity of T.chinensis may be attributable to its high outcrossing rate, wide distribution, and long life span. Generally, the widely distributed long-lived plant species often maintain relatively high genetic diversity compared with narrow geographical ranges or endangered species (Wang, Bernhardsson & Ingvarsson, 2020; Hamrick, Godt & Shermann-Broyles, 1992). T. chinensis is an insect-pollinated plant with high outcrossing rate, although sometimes its pollen is dispersed via wind to the stigma of another plant. Our results showed there was no clonal reproduction in this species (data not shown) despite that it can be reproduced by cuttings and root suckers. It has a short generation time of 1-2y and its raceme flowers produce a large number tiny hairy seeds which can easily be transported by wind and/or water. Its seed can tolerate high salt concentration (0.15 g/l) (Zhang & Zhang, 2019) for a long time.

Our measurement of population genetic diversity showed significantly higher level than that of Liang et al. (2018a) where the genetic diversity ranged from 0.366 to 0.740 for populations along the Yellow River at six SSR markers and that of Zhu et al. (2016) where an average population heterozygosity of 0.493 was detected with 7 SSR markers in 26 populations from Yellow River Delta, and also that of its congener T. taklamakanensis (7SSRs, He = 0.432) (Su et al., 2017). The relatively low genetic diversity detected by those T. chinensis studies may be due to their usage of the microsatellite markers sets which developed through biotinylated-oligonucleotide capture method from T. ramosissima Ledeb. (Gaskin, Pepper & Manhart, 2006). These 10 microsatellite markers actually showed relatively low polymorphism, altogether two of which amplified length-variant alleles when PCR products were separated on PAGE gels (Zhang, 2011).

Population genetic differentiation and structure

T. chinensis populations showed a low level of population differentiation (F_ST = 0.0532) (Table 2), in comparison with the result from its congener T. taklamakanensis (Fst =0.21) (Su et al., 2017). The level of population differentiation was consistent with findings from two population genetic studies of T. chinensis with SSR markers: (Zhu et al., 2016) (F_ST = 0.0533) and Liang et al. (2018a) (F_ST = 0.058). The low level of population differentiation has also been revealed by using RAPD markers (Zhao et al., 2008) and ISSR markers (Jiang, Chen & Bao, 2012). Overall, high outcrossing rate and gene flow facilitate the formation of continuously distributed homogenous populations across its natural range of warm-temperate and subtropical zones, thus resulting in comparatively low population genetic differentiation. With highly polymorphic SSR markers, we robustly identified population genetic structure corresponding to populations’ geographical locations. Our results detected three genetic clusters: the northern cluster including the four populations of Yellow River Delta, the fifth Bohai costal population CY, and Beijing population, the second genetic cluster including two populations of Huang Zhou Bay, and the third genetic cluster containing the remaining Hang Zhou Bay population. The northern (warm-temperate zone) populations were clustered together as one genetic group and populations of Hang Zhou Bay fell into two other genetic clusters, which may imply the different origins of the studied populations of the two different climate zones. Previous studies revealed more nucleus ribosome DNA (nrDNA) haplotypes in the populations of the lower reaches of the Yellow River, implying this species may have dispersed eastward along the direction of the flowing water (Sun et al., 2019). The continuously distributed populations in the Bohai Sea coastal flats are possibly the results of the northward and southward dispersal of seeds from Yellow River Delta by sea currents, which has yet to be tested by extensive sampling along coastal flats of the Bohai Sea and new polymorphic markers (e.g., chloroplast microsatellite markers) for population genetic analysis. The inland population FS from Beijing was clustered together with populations of Bohai Bay, which is consistent with a study where a Beijing population is genetically closer to a population of the Bohai Sea coastal flat from Korea (Lee, Gaskin & Young, 2018b). The subtropical (the East Sea) populations were genetically different from populations from the Bohai Sea coastal flats, which may be the result of seed dispersal via the Yellow River at different times when 150 years ago the Yellow River gorged into the Huanghai Sea instead of the Bohai Sea. Yet we did not include populations from the old estuary flat of the Yellow River in the Huanghai Sea coast in this study.

Isolation pattern of T. chinensis

Using Mantel tests, we detected a highly significant correlations between genetic distance and geographical and environmental variables. Because the autocorrelation between geographical and environmental variables is also high (R² = 0.9357), we performed partial Mantel tests to disentangle the role of geographic distance and environmental variables in shaping the genetic variation detected by microsatellite markers. Our high correlation (R² = 0.5224 for IBD and R² = 0.5225 for IBE) between genetic distance and geographical and environmental variables contrasts with no or low significant correlation (R² = 0.1773, p < 0.01, (Liang et al., 2018a); p = 0.39, (Zhu et al., 2016)) detected before, which may partly be due to their relatively small sampling areas (Yellow River Delta for Zhu et al., 2016, and the Yellow River reaches for Liang et al., 2018a) and the usage of less variable microsatellite loci.