Impact Factor 3.517 | CiteScore 3.60
More on impact ›

Original Research ARTICLE

Front. Genet., 09 August 2019 | https://doi.org/10.3389/fgene.2019.00742

Testing the Effect of Mountain Ranges as a Physical Barrier to Current Gene Flow and Environmentally Dependent Adaptive Divergence in Cunninghamia konishii (Cupressaceae)

Yi-Shao Li1, Kai-Ming Shih1†, Chung-Te Chang2†, Jeng-Der Chung3 and Shih-Ying Hwang1*
  • 1School of Life Science, National Taiwan Normal University, Taipei, Taiwan
  • 2Department of Life Science, Tunghai University, Taichung, Taiwan
  • 3Division of Silviculture, Taiwan Forestry Research Institute, Taipei, Taiwan

Populations can be genetically isolated by differences in their ecology or environment that hampered efficient migration, or they may be isolated solely by geographic distance. Moreover, mountain ranges across a species’ distribution area might have acted as barriers to gene flow. Genetic variation was quantified using amplified fragment length polymorphism (AFLP) and 13 selective amplification primer combinations used generated a total of 482 fragments. Here, we tested the barrier effects of mountains on gene flow and environmentally dependent local adaptation of Cunninghamia konishii occur in Taiwan. A pattern of genetic isolation by distance was not found and variation partitioning revealed that environment explained a relatively larger proportion of genetic variation than geography. The effect of mountains as barriers to genetic exchange, despite low population differentiation indicating a high rate of gene flow, was found within the distribution range of C. konishii. Twelve AFLP loci were identified as potential selective outliers using genome-scan methods (BAYESCAN and DFDIST) and strongly associated with environmental variables using regression approaches (LFMM, Samβada, and rstanarm) demonstrating adaptive divergence underlying local adaptation. Annual mean temperature, annual precipitation, and slope could be the most important environmental factors causally associated with adaptive genetic variation in C. konishii. The study revealed the existence of physical barriers to current gene flow and environmentally dependent adaptive divergence, and a significant proportion of the rate of gene flow may represent a reflection of demographic history.

Introduction

Understanding connectivity among populations and pattern of gene flow has been important for deciphering how species adapt in response to environments and has clear practical implications for conservation of forest genetic resources (Allendorf et al., 2010; Aitken and Bemmels, 2016). Low levels of population genetic differentiation are commonly found in coniferous species due to efficient gene flow based on paternally inherited plastid DNA (ptDNA) variation (Neale et al., 1991; Hipkins et al., 1994; Hwang et al., 2003) and biparentally inherited molecular markers, such as allozyme (Hamrick et al., 1992; Lin et al., 1993; Hamrick and Godt, 1996; Lin et al., 1998). Although gene flow maintains genetic diversity and is critical to population resilience and persistence (Allendorf et al., 2013), high rate of gene flow among closely related sexual populations can preclude local adaptation by eroding divergence driven by natural selection due to the homogenizing effect of dispersing individuals mating among populations (García-Ramos and Kirkpatrick, 1997), resulting in no population adaptive divergence (Li et al., 2018). Nevertheless, selection invoked by environmental heterogeneity can counteract the effect of gene flow (Lenormand, 2002). The failure of detecting variation under selection might have also related to past demographic histories and stochastic mechanisms (Lande, 1988; Wang et al., 2017). It is likely that the balance between migration and selection is species dependent and is related to the realized ecological niche where species can tolerate in response to ecological factors (Pulliam, 2000; Bruno et al., 2003).

The accumulation of genetic variation can provide as raw materials for evolutionary potential under natural selection (Petit and Hampe, 2006; Barrett and Schluter, 2008). The investigation involved orthologous coding sequences in gymnosperms and angiosperms showed that the average synonymous mutation rate was found to be higher in angiosperms in contrast to gymnosperms (Buschiazzo et al., 2012). However, the rate of nonsynonymous substitutions in protein-coding genes was found to be higher in conifers compared with angiosperms suggesting that conifers harbored higher number of fixed adaptive mutations than angiosperms (Buschiazzo et al., 2012), and it is common to observe population local adaptation of conifers in response to environmental heterogeneity (e.g., Mimura and Aitken, 2010; Grivet et al., 2011; Chen et al., 2012; Fang et al., 2013; Shih et al., 2018). However, no adaptive genetic divergence between populations of Taiwania cryptomerioides occurring in Taiwan was detected probably due to random neutral drift because of the high rate of gene flow among populations that eroded the adaptive local mutations (Li et al., 2018).

Natural selection driven by ecological factors will result in the development of ecological adaptation and divergence, and selection can act on genetic variation (e.g., Holderegger et al., 2008; Chen et al., 2012; Fang et al., 2013). Genetic variation in natural populations of a species can be quantified using methods involving next-generation sequencing (NGS), such as restriction site-associated sequencing (Shih et al., 2018). Although less powerful than NGS technologies, amplified fragment length polymorphism (AFLP) is an efficient approach allows generation of hundreds of molecular markers from genome sequences of nonmodel organisms to identify candidate genetic variation involved in adaptive evolution that derived from DNA sequence variation (Holderegger et al., 2008; Chen et al., 2012; Fang et al., 2013).

Cunninghamia konishii Hayata (Cupressaceae) is a coniferous species disjunctly distributed in northern and central Taiwan and along part of the border between Vietnam and Lao People’s Democratic Republic (Nga et al., 2016; http://threatenedconifers.rbge.org.uk/taxa/details/cunninghamia-konishii). C. konishii is closely related to Cunninghamia lanceolata (Lamb). Hook. distributed in mainland China (Lin et al., 1998; Chung et al., 2004) and it may represent a montane form rather than a distinct species (http://threatenedconifers.rbge.org.uk/taxa/details/cunninghamia-konishii). In Taiwan, C. konishii is a dominant species among other conifers associated mainly with Chamaecyparis, Pinus spp., and Pseudotsuga wilsoniana at elevations of 1,300–2,800 m (Liu, 1966). Human interferences have greatly impacted this species, e.g., a pure stand of C. konishii(ca. 60 hectares) in Shiangshanshan, Hsinchu County, northwest Taiwan, was entirely vanished by logging. Low level of population differentiation was revealed for C. konishii based on allozyme (Lin et al., 1998) and ptDNA (Hwang et al., 2003). However, Chung et al. (2004) found high level of genetic differentiation among C. konishii populations based on 357 AFLP markers derived from three AFLP selective amplification primer pairs. Populations of C. konishii are separated not only by major mountain ranges including the Hsuehshan Mountain Range (HMR) and the Central Mountain Range (CMR) trending from north to south, but also by mountains in between of the HMR and the CMR (Figure 1). Mountains were revealed to be the third important factor, after tectonics and climate, as a physical barrier to dispersal in the biogeographic differentiation of a species (Antonelli, 2017). It is likely that mountains may have had an active role in population dispersal of C. konishii. In addition, the relationship between genetic variation in C. konishii and environmental variables implying local adaptation has not been tested. Investigation assesses the association of genetic variation with environment is important for the conservation program because identification of environmentally dependent ecotypes can be crucial in the assisted migration program of C. konishii, particularly in the face of global climate change (Aitken et al., 2008; Aitken and Whitlock, 2013; Aitken and Bemmels, 2016).

FIGURE 1
www.frontiersin.org

Figure 1 Geographic distribution of the 11 populations of C. konishii occur in Taiwan. See Table 1 for abbreviations of the 11 populations of C. konishii.

Environment may play an important role in driving population divergence (Wang et al., 2013; Sexton et al., 2014). However, population differentiation may simply reflect a correlation between genetic and geographic distance independent of environmental conditions (Hutchison and Templeton, 1999). C. konishii populations occupy a wide elevational band (∼1,500 m), though within a small latitudinal range, with wide variety of habitats because of varied geographical topologies and climates (Su, 1984; Li et al., 2013). It is likely that mountains and environment both have been critically involved in shaping the population structure of C. konishii. Therefore, genetic variation of 115 individuals from 11 populations was surveyed using 13 AFLP selective primer pairs. Genetic variation and sample site coordinates were used in testing the hypothesis of mountain ranges as a physical dispersal barrier to genetic exchange and genetic variation together with information of environmental variables were used to test environmentally dependent adaptive divergence underlying local adaptation. The specific aims of the present study were to 1) test the barrier effects of mountains to current gene flow and 2) assess the relationships of genetic variation with environmental variables underlying local adaptation, in the light of usefulness in the future conservation of C. konishii.

Materials and Methods

Sampling and Genotyping

We collected 115 individuals of C. konishii from 11 populations (Table 1, Figure 1). Leaf samples dried with silica gel were used for DNA extraction (Dehestani and Kazemi Tabar, 2007). We quantified genetic variation using AFLP (Vos et al., 1995). In selective amplification, 13 EcoRI-MseI primer combinations (E00: 5’-GACTGCGTACCAATTC-3’ and M00: 5’-GATGAGTCCTGAGTAA-3’) with additional five and three bases were added, respectively, at the ends, were used (Supplementary Table 1). Fragments amplified by polymerase chain reaction were electrophoresed on an ABI 3730XL DNA analyzer (Applied Biosystem, Foster City, CA, USA) and scored with Peak Scanner v.1.0 (Applied Biosystem). The presence and absence of amplified fragments were scored with the fluorescent threshold set at 150 units. Amplified fragment peaks in the range of 150–500 base pair (bp) separated by less than one nucleotide in a ± 0.8 bp window were scored as the same fragment. Markers scored higher than 99% or less than 1% of individuals were removed. The AFLP dataset was deposited in Supplementary Data Sheet 1.

TABLE 1
www.frontiersin.org

Table 1 Site properties and population genetic parameters based on 482 amplified fragment length polymorphism (AFLP) loci of the sampled C. konishii populations.

Genotyping error rate of each primer combination was estimated based on the ratio of mismatches in three amplification replicates of three samples in each population. Loci with error rates greater than 5% were removed (Bonin et al., 2004) and the mean error rate was 1.91% (Supplementary Table 1).

Genetic Diversity

The percentage of polymorphic loci (%P) at the 95% level with rarefaction to the smallest sample size (n = 3) was calculated using AFLPDIV software v.1.1 (Coart et al., 2005). The private band richness (PBr) was estimated separately for each population using ADZE software v.1.0 (Szpiech et al., 2008) with rarefaction (Kalinowski, 2004). We calculated index of association (IA) (Brown et al., 1980) and modified index of association (rD) (Agapow and Burt, 2001) using the ia function of R poppr package (Kamvar et al., 2014; Kamvar et al., 2015) in the R environment (R Development Core Team, 2018) to assess the level of multilocus linkage disequilibrium (LD), and significant departure from zero was tested with 999 permutations. Unbiased expected heterozygosity (uHE) (Nei, 1987) within a population were computed using AFLP-SURV v.1.0 (Vekemans et al., 2002). In AFLP-SURV, allele frequencies were estimated assuming Hardy–Weinberg equilibrium with nonuniform prior distribution (Zhivotovsky, 1999) for uHE computation. Per locus uHE was calculated using ARLEQUIN v.6.0 (Excoffier and Lischer, 2010). Assessment of significant difference of mean uHE per locus between populations was performed using a linear mixed effect model (LMM) with population and locus treated as fixed and random effect, respectively, based on reduced maximum likelihood estimation using the lmer function of R lme4 package (Bates et al., 2015). Significance was assessed based on type II Wald χ2 test using the Anova function of R car package (Fox and Weisberg, 2011), and P values were adjusted with Tukey’s post hoc method using the lsmeans function of R emmeans package (Lenth, 2018).

Genetic Differentiation

The level of genetic differentiation among populations was analyzed via analysis of molecular variance (AMOVA) using the poppr.amova function of R package poppr, and significance tested using the randtest function of R package ade4 with 9,999 permutations (Dray and Dufour, 2007). Across population FST was computed using AFLP-SURV with 9,999 permutations. ARLEQUIN was used to compute pairwise FST with 10,000 permutations. Moreover, HICKORY v.1.1 (Holsinger and Lewis, 2003) was used to estimate θII (an analog to FST), considering the uncertainty associate with the inbreeding coefficient (f) for dominant markers. To check the convergence of parameters, two HICKORY runs were performed with sampling and chain length parameters include burnin = 5,000, samples = 100,000, and thinning = 20. Genetic data was evaluated fitting to the full model, f = 0 model, θII = 0 model, and f-free model, and the best fitting model in estimating θII was determined by the deviance information criterion (DIC). The f-free model was used to estimate f within populations because other models produced unreliable f estimates. The f-free model selects f values from a uniform prior distribution without generating a posterior distribution of f independent of assumptions.

Genetic Clustering and Testing the Barrier Effects of Mountains on Gene Flow

Genetic homogeneous groups of C. konishii individuals were assessed based on sparse nonnegative factorization (sNMF) algorithm (Frichot and Francois, 2015) and discriminant analysis of principal components (DAPC) (Jombart et al., 2010). Individual assignments based on the least-squares optimization using the snmf function of R LEA package (Frichot and Francois, 2015) was performed for K = 1–11. The regularization parameter, iterations, and repetitions in snmf were set to 100, 200, and 10, respectively, with other arguments set to defaults. The best K in LEA was evaluated with the means of minimal cross-entropy (CE). The find.clusters and dapc functions of R adegenet package (Jombart and Ahmed, 2011) were used in DAPC analysis, in which a principal component analysis (PCA) was first performed and followed by a discriminant analysis that maximize variance between groups. The best K in DAPC was indicated by an elbow in the curve of Bayesian information criterion (BIC), as suggested by the authors of the R adegenet package, estimated using the find.clusters function.

Given the heterogeneous nature of C. konishii populations distributed across mountains, Monmonier’s maximum difference algorithm (Monmonier, 1973; Manni et al., 2004) was used to assess sharp genetic discontinuities in the geographic distribution area of C. konishii, based on the Euclidean distance of AFLP and sample coordinates, using the optimize.monmonier function of R adegenet package. In optimize.monmonier, 10 different starting points were used to find the largest sum of local distances that explains genetic distances among populations.

Environmental Variables and Heterogeneity

Environmental variables used were 19 bioclimate, 2 topological, and 12 ecological variables. Nineteen bioclimate variables for sample sites at 30-s spatial resolution (∼1 km) were downloaded from the WorldClim v.1.4 (Hijmans et al., 2005). Topographic variables including aspect and slope at 30-m resolution were obtained from Aster Global Digital Elevation Map (https://asterweb.jpl.nasa.gov/gdem.asp). We obtained ecological factors including normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) derived from moderate resolution imaging spectroradiometer (MODIS) dataset MOD13A2 (1-km resolution), and leaf area index (LAI) and fraction of absorbed photosynthetically active radiation (fPAR) derived from MOD15A2 dataset (500-m resolution). The annual total potential evapotranspiration (PET) was calculated based on MOD16A3 dataset (500-m resolution). All the MODIS datasets were acquired from Land Process Distributed Active Archive Center (http://lpdaac.usgs.gov) during 2001‒2013, and monthly mean values of NDVI, EVI, LAI, and fPAR were computed using a maximum-value composite procedure (Huete et al., 2002). Monthly mean values of the other five ecological factors including relative humidity (RH), cloud cover (CLO), time of sunshine (SunH), number of rainfall days per year (RainD), and mean wind speed (WSmean) were also calculated for data obtained from the Data Bank for Atmospheric & Hydrologic Research (https://dbahr.pccu.edu.tw/, recorded in 1990–2013) at spatial resolution of 1 km using a universal spherical model of the Kriging method in ArcGIS (Chang et al., 2014). Additionally, we acquired soil pH values of sample sites based on an island-wide soil investigation (n = 1,150) conducted in 1969‒1986 from the Agriculture and Food Agency of Taiwan (Chang et al., 2009). Annual precipitation and annual PET (derived from annual mean temperature) were used to calculate annual moisture index (Thornthwaite, 1948).

The cor function of R was used to calculate correlation coefficients between environmental variables. Variance inflation factor (VIF) was calculated using the vif function of R package usdm (Naimi et al., 2014). We retained eight environmental variables for further use based on environmental variables with VIF < 10 and which not strongly correlated with other variables (|r| < 0.8). The eight retained environmental variables were aspect, BIO1 (annual mean temperature), BIO7 (annual temperature range), BIO12 (annual precipitation), NDVI, PET, RainD, and slope (Supplementary Table 2). Environmental Euclidean distance matrix was used in a permutational multivariate analysis of variance (PERMANOVA) to assess environmental heterogeneity among sample sites and among population genetic clusters (revealed by LEA and DAPC analyses, see Results) using the adonis function of R package vegan (Oksanen et al., 2017). Pairwise population and cluster comparisons were also performed using the pairwise.perm.manova function of R package RVAideMemoire (Herve, 2018). Significant pairwise comparisons were tested with 999 permutations and a false discovery rate (FDR) of 5%.

Mantel test was used to assess genetic isolation by distance (IBD) by analyzing the correlation of the population Euclidean distance matrix of AFLP with the population Euclidean distance matrix of geography (latitude and longitude) using the mantel function of R vegan package.

Partitioning of Genetic Variation Explained by Environmental Variables

The eight retained environmental variables were used in a redundancy analysis (RDA) to assess the relative contribution of environmental variables explaining the total AFLP variation using the varpart function of R package vegan, and significance tested using the anova.cca function with 999 permutations (Oksanen et al., 2017). The total variation was partitioned into four fractions: pure environmental variables (fraction [a]), geographically structured environmental variables (fraction [b]), pure geographic variables (fraction [c]), and residual effects (fraction [d]) (Borcard et al., 1992; Borcard and Legendre, 2002). Adjusted R2 value was used to represent the amount of variation explained in each fraction (Peres-Neto et al., 2006). Longitude and latitude of sample sites were used as geographic variables in the analysis.

Test for FST Outliers

Two genome scan methods, BAYESCAN and DFDIST, were used to identify FST outliers indicating signature of selection across populations. Hierarchical Bayesian method implemented in BAYESCAN v.2.1 (Foll and Gaggiotti, 2008) uses a reversible-jump Markov chain Monte Carlo algorithm to estimate the ratio of posterior probabilities of selection over neutrality [the posterior odds (PO)]. We used 100 pilot runs of 50,000 iterations followed by a sample size of 50,000 with thinning interval of 20 among 106 iterations in BAYESCAN analysis. A logarithmic scale for model choice of selection over neutrality can be defined as: substantial (log10PO > 0.5); strong (log10PO > 1.0); very strong (log10PO > 1.5); and decisive (log10PO > 2) (Jeffreys, 1961; Foll, 2012). We considered a locus with log10PO > 0.5 as a potential selective outlier under selection.

DFDIST incorporates the Beaumont and Nichols model modified for AFLP (Beaumont and Nichols, 1996). The probability of a locus that may be under selection by observed FST (Weir and Cockerham, 1984) and uHE (Zhivotovsky, 1999) compared to simulated neutral distributions was estimated. Parameters for running DFDIST were: critical frequency = 0.99; Zhivotovsky parameters = 0.25; trimmed mean FST = 0.3 (excluding 30% of highest and 30% of lowest FST values); smoothing proportion = 0.06; 500,000 resamplings; critical P = 0.05; and level of differentiation (target average θ) = 0.063254. Empirical loci considered to be the outliers potentially under directional selection were those with FST values significantly greater (P < 0.01) than the simulated distribution.

Test for Associations of Genetic Variation With Environmental Variables

Associations of genetic variation with environmental variables were assessed using Samβada v.0.8.1 (Stucki et al., 2017) and latent factor mixed model (LFMM) (Frichot et al., 2013). A multiple univariate logistic regression approach was employed in Samβada to test for significant correlations of allele frequencies with the values of environmental variables. Models with and without environmental variables were compared, and significant fit was determined based on both Wald and G scores with an FDR cutoff of 0.05. LFMM is a method uses a hierarchical Bayesian mixed model considering background levels of population structure as random effects due to demographic history and IBD patterns, and is introduced via latent factors, K. In LFMM, genetic data matrix was used as fixed effect in a testing procedure based on Z-scores. The number of latent factors, K, was set to 3 (based on the LEA and DAPC clustering results). We ran LFMM five times for each value of K with 10,000 iterations in Gibbs sampling algorithm and a burn-in period of 5,000 cycles. Z-scores from five independent replicate runs were combined using Fisher–Stouffer method (Brown, 1975), and the resulting P values were adjusted using the genomic inflation factor (λ). P values were further adjusted based on an FDR correction of 1% using the R qvalue package (Storey et al., 2019).

For those potential selective outliers identified by either genome scan method (BAYESCAN or DFDIST), a Bayesian logistic regression approach implemented in the stan_glm function of R rstanarm package (Goodrich et al., 2018) was also used to test for the associations of the potential selective FST outliers with environmental variables. In stan_glm analysis, Student’s t distribution with mean zero and 7 degrees of freedom was used as the weakly informative prior. The scale of the prior distribution was 10 for the intercept and 2.5 for the predictors. All models were run with four chains for 1,000 warm-up and 1,000 sampling steps. Credible intervals (95% and 99%) of estimate were calculated using the posterior_interval function of R rstanarm package.

Results

Diversity and Differentiation

We obtained 482 AFLP loci (mean ± SD: 37.08 ± 7.87) for C. konishii (Supplementary Table 1). The average value of PBr ± SD was 0.017 ± 0.007 (ranged from 0.0098 for the population YH to 0.0296 for the population DY) (Table 1). The average value of %P ± SD was 53.9 ± 15.1 (ranged from 32.2 for the population DY to 74.3 for the population KW). The average level of uHE ± SD was 0.252 ± 0.016 and ranged between 0.227 (population DT) and 0.284 (population CT). Within C. konishii, the level of uHE per locus was significantly different among populations using LMM (χ2 = 41.38, P < 0.0001). However, significant difference in the level of uHE per locus was not found for all pairwise comparisons (Supplementary Table 3). No positive correlation of sample size with population uHE (S = 305.78, ρ = −0.390, P = 0.2358) and PBr (S = 251.29, ρ = −0.142, P = 0.6766), respectively, was found based on Spearman’s rank correlation test. However, significant positive correlation was found between sample size and the percentage of polymorphism (S = 14.12, ρ = 0.936, P < 0.001). Using HICKORY, an estimate of f within population was 0.4969 (95% CI: 0.022–0.976) using the f-free model representing the contemporary reproductive mode of C. konishii (Supplementary Table 4). The measures of multilocus LD, IA and rD, showed significant departure from random association between AFLP loci in 7 (populations AL, DT, DY, KW, SK, TJ, and TS) of the 11 populations examined (Table 1).

The full model was the best model fitting to the genetic data based on DIC in HICKORY analysis (Supplementary Table 4). Low, albeit significant, level of genetic differentiation of C. konishii populations was found based on AMOVA and θIIST = 0.0827, θII = 0.0848; Table 2, Supplementary Table 4). Across population FST estimated using AFLP-SURV was 0.0656 (P < 0.001). Significant population pairwise FST estimated using ARLEQUIN was commonly found (Supplementary Table 5).

TABLE 2
www.frontiersin.org

Table 2 Summary of the analysis of molecular variance (AMOVA).

Genetic Clustering and Test for IBD

In the LEA analysis, the means of minimal CE was minimized at K = 3 (Supplementary Figure 1A). LEA clustering results were depicted for K = 2–3 (Figure 2). The lowest BIC was found at K = 5, but elbowed at K = 3, in DAPC analysis (Supplementary Figure 1B). The eigenvalues for the first two retained PCs were 1.04 and 0.64 (Supplementary Figure 1C) and the eigenvalues for the first two DAPC linear discriminants were 28.77 and 14.63 (Supplementary Figure 1D). Three population clusters can be distinguished based on DAPC: cluster 1 contains population AL; cluster 2 contains populations DT, TS, and YH; and cluster 3 contains populations AM, CT, DY, KW, SK, SL, and TJ. However, a few individuals from one population can be grouped with individuals of other populations in a separate cluster (Figure 3). Figure 4 displayed the spatial boundaries of genetic exchange among C. konishii populations using Monmonier’s algorithm. Mantel test revealed no significant genetic isolation by geographic distance (rM = 0.251, P = 0.107).

FIGURE 2
www.frontiersin.org

Figure 2 Individual assignments of 115 individuals from 11 populations of C. konishii analyzed using LEA. The clustering scenarios for K = 2–3 were displayed.

FIGURE 3
www.frontiersin.org

Figure 3 Clustering results analyzed using discriminant analysis of principal components (DAPC).

FIGURE 4
www.frontiersin.org

Figure 4 Barriers to genetic exchange identified using Monmnier’s algorithm. The levels of physical barrier effect of gene flow were represented by the thickness of blue lines. See Table 1 for abbreviations of the 11 populations of C. konishii.

Variation Partitioning of Genetic Variation Explained by Environment and Geography

PERMANOVA revealed significant environmental difference among (P = 0.001) and between the three genetic clusters (P = 0.001). Significant pairwise differences between populations were found (Supplementary Table 6), albeit PERMANOVA revealed no environmental difference across populations (P = 1). Total explainable genetic variation was 7.38% (Table 3). This value was relatively small as compared to the amounts of unexplained variation (92.62%). Nevertheless, significant amount of genetic variation was explained by both pure environment (5.18%, fraction [a]) and pure geography (1.84%, fraction [c]), albeit small, based on F tests of RDA model. Moreover, environment explained a relatively larger proportion of genetic variation than geography.

TABLE 3
www.frontiersin.org

Table 3 The percentage of variation explained in genetic loci accounted for by non-geographically structured environmental variables [a], shared (geographically structured) environmental variables [b], pure geographic factors [c], and undetermined component [d] analyzed based on the eight retained environmental variables.

Potential Outliers Strongly Correlated With Environmental Variables

Six and nine nonoverlapped loci, respectively, were identified as FST outliers by BAYESCAN and DFDIST (Table 4). Allele frequencies of these outlier loci across populations arranged either latitudinally or longitudinally were depicted in Figure 5. Samβada and LFMM found significant associations of 18 and 83 loci, among the 482 loci scored, respectively, with at least one environmental variable. The Wald and G scores in the Samβada analysis were reported in Supplementary Table 7. The corresponding Z-score, −log10(P value), and adjusted P value of candidate outlier loci identified using LFMM were summarized in Supplementary Table 8. For the 15 loci identified either by BAYESCAN or DFDIST, Samβada found five (P1_1715, P5_2456, P9_1014, P11_1715, and P15_1918) and LFMM found seven (P1_1409, P1_1715, P5_2456, P11_1715, P12_2853, P13_1547, and P15_1446) loci, respectively, strongly correlated with environmental variables (Table 4). Of the 15 loci identified either by BAYESCAN or DFDIST, nine (P1_1715, P4_1326, P5_1061, P5_2456, P9_1014, P11_1715, P13_1547, P15_1918, and P18_1421) were found to be strongly associated with environmental variables using Bayesian logistic regression (the stan_glm function of R package rstanarm). In summary, 12 AFLP loci (2.49%) were identified either by BAYESCAN or DFDIST and associated strongly with environmental variables using Samβada, LFMM, and rstanarm. Results showed that annual mean temperature, annual precipitation, and slope were closely associated with a large proportion of the 12 FST outliers identified by BAYESCAN and DFDIST; and aspect, annual temperature range, NDVI, PET, and RainD associated strongly only with a minor proportion of these FST outliers. (Table 4).

TABLE 4
www.frontiersin.org

Table 4 Potential outliers associated with environmental variables.

FIGURE 5
www.frontiersin.org

Figure 5 Heatmap of allele frequencies of the 15 outlier loci identified by either BAYESCAN or DFDIST. The sequence of populations was arranged according to (A) latitude or (B) longitude.

Discussion

Genetic Diversity

The average uHE ± SD (= 0.252 ± 0.016) in the present study was found to be higher compared to that of the previous study examined for C. konishii (uHE = 0.184; Chung et al., 2004). The higher number of amplification primer combinations used in the present study could be the cause of the discrepancy compared with the previous AFLP study, which may result in higher interlocus variance (Nei, 1987). Although there is no direct evidence for the relationship between the number of amplification primer pairs and the level of genetic diversity, higher interlocus variance was found to contribute more to the total variance associated with the level of genetic diversity compared to intralocus variance in Pinus pinaster (Ribeiro et al., 2002). Moreover, more amplified AFLP loci, though with higher error rate, could be positively correlated with higher level of genetic diversity (Zhang and Hare, 2012). Additionally, higher number of amplification primer pairs used could have amplified fragments spanning more chromosomal loci and contributed to the higher level of genetic diversity.

Based on the result of the present study, the average level of C. konishii genetic diversity appeared to be comparable to other conifers occurring in Taiwan that have been examined based on AFLP such as Keteleeria davidiana var. formosana (average uHE = 0.233, Fang et al., 2013) and T. cryptomerioides (average uHE = 0.236, Li et al., 2018). The average level of C. konishii genetic diversity was also found to be comparable to the average level of AFLP diversity (= 0.230) for 13 plant species in general summarized in Nybom (2004). Despite C. konishii covers only a small latitudinal range, it distributed across complex topography and habitat heterogeneity in mountain ranges (Su, 1984; Li et al., 2013) and hence could have contributed to the accumulation of genetic diversity (Hewitt, 1996; Petit et al., 2003; Hewitt, 2004).

Mountain Ranges as a Physical Barrier to Genetic Exchange

The spatial distribution of genetic variability within natural plant populations can be shaped by population demographic history and geography through the combined effects of genetic drift, gene flow, and selection (Hoffmann and Sgrò, 2011). The discrepancy of the level of genetic differentiation based on AFLP was also observed as that of the level of genetic diversity between the present (Table 3) and the previous study (Chung et al., 2004). The ΦST, FST, and θII estimates are consistent with the results of clustering analyses using LEA (Figure 2) and DAPC (Figure 3). In the present study, the level of population genetic differentiation was much smaller compared with the previous AFLP study based on AMOVA analysis (ΦST = 0.0827 vs. ΦST = 0.2460, respectively). However, the level of genetic differentiation of the present study is consistent with the levels of genetic differentiation found based on allozyme (FST = 0.029) (Lin et al., 1998) and ptDNA (GST = 0.073) (Hwang et al., 2003). The higher level of genetic differentiation among populations found in the previous study (Chung et al., 2004) could have arisen in part due to the smaller number of genomic regions amplified, compared with the present study, and harvested a relatively smaller fraction of ancestral polymorphisms contributing to the higher level of population differentiation. Nonetheless, the levels of population genetic structure estimated (ΦST, FST, and θII) in the present study were consistent with the results found in other conifers occurring in Taiwan including K. davidiana var. formosana (Fang et al., 2013) and T. cryptomerioides (Li et al., 2018) based on AFLP. In addition, these results are in accordance with the general realization that conifers typically have high rates of effective pollen dispersal resulting in low population differentiation because of long distance wind pollination based on allozyme data (Hamrick et al., 1992; Hamrick and Godt, 1996).

Although low levels of across population and population pairwise FST were found (Table 2, Supplementary Table 5), genetic boundaries have been established by mountain ranges (Figure 4) and current gene flow among populations might have also been restricted by spatial environmental heterogeneity as suggested by PERMANOVA analysis (Supplementary Table 6). Additionally, Mantel test revealed no significant correlation of genetic distance matrix with geographic distance matrix suggesting no IBD pattern. These results suggest that the pattern of gene flow in C. konishii may not follow a stepping stone migration model (Kimura, 1953) due to the presence of physical dispersal barriers and sharp discontinuities in gene frequencies can be observed in populations in close proximity geographically (Figures 1 and 5). Our results suggest that the presence of physical barriers to gene flow and migration–drift equilibrium may not be reached, and a significant proportion of the observed low level of genetic differentiation reflects historical population expansion (Hwang et al., 2003) rather than current level of gene flow (Latta and Mitton, 1999; Turgeon and Bernatchez, 2001; Ringbauer et al., 2018). In the present study, the barrier effects of the HMR and the CMR and mountains in between of the HMR and the CMR could have played roles in shaping population differentiation of C. konishii. Moreover, mountains surrounding the AL population could have played a damping effect on dispersal between the AL and all other C. konishii populations, particularly between the AL and its nearby CT population, genetically (Figure 4). Therefore, natal dispersal within different geographic areas could have played an important role in shaping the genetic structure of C. konishii. However, environment could also be an important factor contributed to C. konishii population divergence because genetic variation is structured more by pure environment than by pure geography as revealed by the analysis of variation partitioning (Table 3). Our results suggest that environmental factors accompanied with physical separations have shaped the divergence between C. konishii populations. Ecotypes are likely to have arisen at local scales and the emergence of clusters of locally beneficial mutations forming genomic islands of divergence (Via, 2009; Wolf et al., 2010; Strasburg et al., 2012).

The Most Important Environmental Variables for Adaptive Genetic Divergence in C. konishii

Current mating of C. konishii individuals among populations could be limited as suggested by significant population-level IA and rD in 7 of the 11 populations examined, indicating nonrandom association of loci within these populations (Table 1). Local reduction in gene flow mediated by divergent selection and/or immigrant inviability that reduced survival and reproduction of immigrants might have also contributed to the detection of multilocus LD (Kawecki and Ebert, 2004; Nosil et al., 2005; Jump and Peñuelas, 2006). However, detection of LD does not ensure a lack of linkage equilibrium (Slatkin, 2008) and LD within populations may be related to the demographic history of species (Wall et al., 2002).

Our results based on analyses using BAYESCAN, DFDIST, Samβada, LFMM, and rstanarm found adaptive divergence of genetic variation (Table 4). Although BAYESCAN and DFDIST found no common outlier loci, regression approaches (Samβada, LFMM, and rstanarm) can be effective in testing their associations with environmental gradients (Stucki et al., 2017). Results showed that outlier AFLP loci were mainly driven by environmental variables including annual temperature range, annual precipitation, and slope (Table 4), and these environmental factors could have played major roles in driving outlier AFLP variation for local adaptation in C. konishii. Environmental variables including aspect, annual temperature range, NDVI, PET, and RainD could have played only minor roles as selective drivers in shaping adaptive divergence in C. konishii.

Strong correlations between environmental variables and population genetic variability can be provided as evidence for local adaptation and selection (Linhart and Grant, 1996). Temperature and precipitation have been found to be the two most important selective drivers for local adaptation in conifers (e.g., González-Martínez et al., 2006; Fang et al., 2013; Shih et al., 2018). Temporal and spatial variation in temperature and precipitation can influence fitness-related traits (Parmesan and Yohe, 2003; Hoffmann and Sgrò, 2011; Franks and Hoffmann, 2012) and consequently influence the survival of conifers (Brodribb et al., 2014). Topographic gradients over short distances in the rugged geographic landscape can play significant roles in shaping species composition of forest communities (Kitagawa et al., 2015). Topographic factors such as slope are important predictors of forest attributed to differences in radiation exposure and have a strong influence on the microclimate (Rosenberg et al., 1983; Bennie et al., 2008; Brousseau et al., 2015). Aspect and slope were found to be closely correlated with genetic variation within and between species (Nakazato et al., 2010; Monahan et al., 2012; Bothwell et al., 2013; Brousseau et al., 2015; Huang et al., 2015a; Huang et al., 2015b; Li et al., 2018). In the present study, slope could be the major topographic factor that played an important role in driving outlier AFLP variation for local adaptation (Table 4).

Environmental variation has been correlated with forest ecosystem properties in various geographic regions, suggesting the sensitivity of species composition responding to environmental conditions (Vandewalle et al., 2010; Roscher et al., 2012; Zhang et al., 2014). NDVI is a proxy to photosynthetic activity representing the level of vegetation greenness and has been shown to be correlated with intraspecific and interspecific adaptive divergence (Nakazato et al., 2010; Huang et al., 2015a; Huang et al., 2015b; Chen et al., 2017; Li et al., 2018). RainD represents the number of rainfall days per year, which was found to be associated with adaptive genetic variation of an angiosperm species, Rhododendron oldhamii, occur in Taiwan (Hsieh et al., 2013). PET was identified to be closely correlated with adaptive genetic variation between genetic lineages of a coniferous species, T. cryptomerioides, occurring in Taiwan, mainland China, and Vietnam (Li et al., 2018).

Our results suggest that ecologically relevant selective drivers involved in population adaptive divergence of C. konishii. Spatial environmental heterogeneity might have invoked genetic divergence among populations and in consequence the formation of local adaptation.

Conclusions

The maintenance of genetic ecotypes adapting to varying environmental conditions can be critical for the conservation of species in the face of global climate change. Environmental variables might have exerted effects in shaping adaptive evolution of genetic variation. In the present study, we found annual mean temperature, annual precipitation, and slope could be the most important environmental factors that played crucial roles in shaping adaptive genetic divergence. Nonetheless, aspect, annual temperature range, NDVI, PET, and RainD, though might have played only minor roles, could still be important selective drivers of local adaptation in C. konishii populations. In addition, geographic barriers to gene flow exerted by mountain ranges could have played an important role in restricting adaptive genetic divergence at local scales. The three genetic clusters revealed by LEA and DAPC can be employed as evolutionarily significant units for the future conservation program of this species.

Author Contributions

S-YH proposed, funded, and designed the research; Y-SL, K-MS, J-DC, and S-YH collected samples; Y-SL, K-MS, and C-TC performed research; S-YH, Y-SL, K-MS, and C-TC analyzed data; S-YH and Y-SL wrote the paper. All authors have read and approved the final manuscript.

Funding

This work was funded by the Taiwan Ministry of Science and Technology under grant number of MOST 106-2313-B-003-001-MY3.

All plant materials collected conformed to the regulations of the Taiwan Ministry of Science and Technology.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be constructed as a potential conflict of interest.

Supplementary Materials

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00742/full#supplementary-material

Data Sheet 1 | Comma-separated value file for scored AFLP markers of the 11 populations (n = 115) of Cunninghamia konishii. The first two columns recorded the information of sample ID and population code. The rest of columns represent presence/absence of markers.

References

Agapow, P. M., Burt, A. (2001). Indices of multilocus linkage disequilibrium. Mol. Ecol. Notes 1, 101–102. doi: 10.1046/j.1471-8278.2000.00014.x

CrossRef Full Text | Google Scholar

Aitken, S. N., Bemmels, J. B. (2016). Time to get moving: assisted gene flow of forest trees. Evol. Appl. 9, 271–290. doi: 10.1111/eva.12293

PubMed Abstract | CrossRef Full Text | Google Scholar

Aitken, S. N., Whitlock, M. C. (2013). Assisted gene flow to facilitate local adaptation to climate change. Annu. Rev. Ecol. Evol. Syst. 44, 367–388. doi: 10.1146/annurev-ecolsys-110512-135747

CrossRef Full Text | Google Scholar

Aitken, S. N., Yeaman, S., Holliday, J. A., Wang, T., Curtis-McLane, S. (2008). Adaptation, migration or extirpation: climate change outcomes for tree populations. Evol. Appl. 1, 95–111. doi: 10.1111/j.1752-4571.2007.00013.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Allendorf, F. W., Hohenlohe, P. A., Luikart, G. (2010). Genomics and the future of conservation genetics. Nat. Rev. Genet. 11, 697–709. doi: 10.1038/nrg2844

PubMed Abstract | CrossRef Full Text | Google Scholar

Allendorf, F. W., Luikart, G., Aitken, S., (2013). Conservation and the genetics of populations. 2nd Edn. Chichester: Wiley-Blackwell.

Google Scholar

Antonelli, A. (2017). Biogeography: drivers of bioregionalization. Nat. Ecol. Evol. 1, 0114. doi: 10.1038/s41559-017-0114

CrossRef Full Text | Google Scholar

Barrett, R. D. H., Schluter, D. (2008). Adaptation from standing genetic variation. Trends Ecol. Evol. 23, 38–44. doi: 10.1016/j.tree.2007.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Bates, D., Maechler, M., Bolker, B., Walker, S. (2015). Fitting linear mixed-effects models using lme4. J. Stat. Soft. 67, 1–48. doi: 10.18637/jss.v067.i01

CrossRef Full Text | Google Scholar

Beaumont, M. A., Nichols, R. A. (1996). Evaluating loci for use in the genetic analysis of population structure. Proc. Roy. Soc. B-Biol. Sci. 263, 1619–1626. doi: 10.1098/rspb.1996.0237

CrossRef Full Text | Google Scholar

Bennie, J., Huntley, B., Wiltshire, A., Hill, M. O., Baxter, R. (2008). Slope, aspect and climate: spatially explicit and implicit models of topographic microclimate in chalk grassland. Ecol. Model. 216, 47–59. doi: 10.1016/j.ecolmodel.2008.04.010

CrossRef Full Text | Google Scholar

Bonin, A., Bellemain, E., Bronken, E. P., Pompanon, F., Brochmann, C., Taberlet, P. (2004). How to track and assess genotyping errors in population genetics studies. Mol. Ecol. 13, 3261–3273. doi: 10.1111/j.1365-294X.2004.02346.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Borcard, D., Legendre, P. (2002). All-scale spatial analysis of ecological data by means of principal coordinates of neighbor matrices. Ecol. Model. 153, 51–68. doi: 10.1016/S0304-3800(01)00501-4

CrossRef Full Text | Google Scholar

Borcard, D., Legendre, P., Drapeau, P. (1992). Partialling out the spatial component of ecological variation. Ecology 73, 1045–1055. doi: 10.2307/1940179

CrossRef Full Text | Google Scholar

Bothwell, H., Bisbing, S., Therkildsen, N. O., Crawford, L., Alvarez, N., Holderegger, R., et al. (2013). Identifying genetic signatures of selection in a non-model species, alpine gentian (Gentiana nivalis L.), using a landscape genetic approach. Conser. Genet. 14, 467–481. doi: 10.1007/s10592-012-0411-5

CrossRef Full Text | Google Scholar

Brodribb, T. J., McAdam, S. A. M., Jordan, G. J., Martins, S. C. V. (2014). Conifer species adapt to low-rainfall climates by following one of two divergent pathways. Proc. Natl. Acad. Sci. U.S.A. 111, 14489–14493. doi: 10.1073/pnas.1407930111

PubMed Abstract | CrossRef Full Text | Google Scholar

Brousseau, L., Foll, M., Scotti-Saintagne, C., Scotti, I. (2015). Neutral and adaptive drivers of microgeographic genetic divergence within continuous populations: the case of the neotropical tree Eperua falcata (Aubl). PLoS ONE 10, e0121394. doi: 10.1371/journal.pone.0121394

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, M. B. (1975). A method for combining non-independent, one-sided tests of significance. Biometrics 31, 987–992. doi: 10.2307/2529826

CrossRef Full Text | Google Scholar

Brown, A. H., Feldman, M. W., Nevo, E. (1980). Multilocus structure of natural populations of Hordeum spontaneum. Genetics 96, 523–536. doi: 10.2307/2529826

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruno, J. F., Stachowicz, J. J., Bertness, M. D. (2003). Inclusion of facilitation into ecological theory. Trends Ecol. Evol. 18, 119–125. doi: 10.1016/S0169-5347(02)00045-9

CrossRef Full Text | Google Scholar

Buschiazzo, E., Ritland, C., Bohlmann, J., Ritland, K. (2012). Slow but not low: genomic comparisons reveal slower evolutionary rate and higher dN/dS in conifers compared to angiosperms. BMC Evol. Biol. 12, 8. doi: 10.1186/1471-2148-12-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, C.-T., Lin, T.-C., Lin, N.-H. (2009). Estimating the critical load and the environmental and economic impact of acid deposition in Taiwan. J. Geogr. Sci. 56, 39–58.

Google Scholar

Chang, C.-T., Wang, S.-F., Vadeboncoeur, M. A., Lin, T.-C. (2014). Relating vegetation dynamics to temperature and precipitation at monthly and annual timescales in Taiwan using MODIS vegetation indices. Int. J. Remote Sens. 35, 598–620. doi: 10.1080/01431161.2013.871593

CrossRef Full Text | Google Scholar

Chen, J., Källman, T., Ma, X., Gyllenstrand, N., Zaina, G., Morgante, M., et al. (2012). Disentangling the roles of history and local selection in shaping clinal variation of allele frequencies and gene expression in Norway spruce (Picea abies). Genetics 191, 865–881. doi: 10.1534/genetics.112.140749

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J.-H., Huang, C.-L., Lai, Y.-L., Chang, C.-T., Liao, P.-C., Hwang, S.-Y., et al. (2017). Postglacial range expansion and the role of ecological factors in driving adaptive evolution of Musa basjoo var. formosana. Sci. Rep. 7, 5341. doi: 10.1038/s41598-017-05256-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, J.-D., Lin, T.-P., Tan, Y.-C., Lin, M.-Y., Hwang, S.-Y. (2004). Genetic diversity and biogeography of Cunninghamia konishii (Cupressaceae), an island species in Taiwan: a comparison with Cunninghamia lanceolata, a mainland species in China. Mol. Phylogenet. Evol. 33, 791–801. doi: 10.1016/j.ympev.2004.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Coart, E., Glabeke, S. V., Petit, R. J., Bockstaele, E. V., Roldań-Ruiz, I. (2005). Range wide versus local patterns of genetic diversity in hornbeam (Carpinus betulus L). Conserv. Genet. 20056, 259–273. doi: 10.1007/s10592-004-7833-7

CrossRef Full Text | Google Scholar

Dehestani, A., Kazemi Tabar, S. K. (2007). A rapid efficient method for DNA isolation from plants with high levels of secondary metabolites. Asian J. Plant Sci. 6, 977–981. doi: 10.3923/ajps.2007.977.981

CrossRef Full Text | Google Scholar

Dray, S., Dufour, A.-B. (2007). The ade4 package: implementing the duality diagram for ecologists. J. Stat. Soft. 22, 1–20. doi: 10.18637/jss.v022.i04

CrossRef Full Text | Google Scholar

Excoffier, L., Lischer, H. E. (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, J.-Y., Chung, J.-D., Chiang, Y.-C., Chang, C.-T., Chen, C.-Y., Hwang, S.-Y. (2013). Divergent selection and local adaptation in disjunct populations of an endangered conifer, Keteleeria davidiana var. formosana (Pinaceae). PLoS ONE 8, e70162. doi: 10.1371/journal.pone.0070162

PubMed Abstract | CrossRef Full Text | Google Scholar

Foll, M. (2012). Bayescan 2.1 user manual. Available online at: http://cmpg.unibe.ch/software/BayeScan/files/BayeScan2.1_manual.pdf.

Google Scholar

Foll, M., Gaggiotti, O. (2008). A genome scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics 180, 977–993. doi: 10.1534/genetics.108.092221

PubMed Abstract | CrossRef Full Text | Google Scholar

Fox, J., Weisberg, S. (2011). An R companion to applied regression, 2nd Edn. Thousand Oaks, CA: Sage. Available online at: http://socserv.socsci.mcmaster.ca/jfox/Books/Companion.

Google Scholar

Franks, S. J., Hoffmann, A. A. (2012). Genetics of climate change adaptation. Annu. Rev. Genet. 46, 185–208. doi: 10.1146/annurev-genet-110711-155511

PubMed Abstract | CrossRef Full Text | Google Scholar

Frichot, E., Francois, O. (2015). LEA: an R package for landscape and ecological association studies. Methods Ecol. Evol. 6, 925–929. doi: 10.1111/2041-210X.12382

CrossRef Full Text | Google Scholar

Frichot, E., Schoville, S. D., Bouchard, G., François, O. (2013). Testing for associations between loci and environmental gradients using latent factor mixed models. Mol. Biol. Evol. 30, 1687–1699. doi: 10.1093/molbev/mst063

PubMed Abstract | CrossRef Full Text | Google Scholar

García-Ramos, G., Kirkpatrick, M. (1997). Genetic models of adaptation and gene flow in peripheral populations. Evolution 51, 21–28. doi: 10.1111/j.1558-5646.1997.tb02384.x

PubMed Abstract | CrossRef Full Text | Google Scholar

González-Martínez, S. C., Krutovsky, K. V., Neale, D. B. (2006). Forest-tree population genomics and adaptive evolution. New Phytol. 170, 227–238. doi: 10.1111/j.1469-8137.2006.01686.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodrich, B., Gabry, J., Ali, I., Brilleman, S. (2018). rstanarm: Bayesian applied regression modeling via Stan. R package version 2.17.4. Available online at: http://mc-stan.org/.

Google Scholar

Grivet, D., Sebastiani, F., Alía, R., Bataillon, T., Torre, S., Zabal-Aguirre, M., et al. (2011). Molecular footprints of local adaptation in two Mediterranean conifers. Mol. Biol. Evol. 28, 101–116. doi: 10.1093/molbev/msq190

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamrick, J. L., Godt, M. J. W. (1996). Effects of life history traits on genetic diversity in plant species. Philos. Trans. R. Soc. Lond. B Biol. Sci. 351, 1291–1298. doi: 10.1098/rstb.1996.0112

CrossRef Full Text | Google Scholar

Hamrick, J. L., Godt, M. J. W., Sherman-Broyles, S. L. (1992). Factors influencing levels of genetic diversity in woody plant species. New For. 6, 95–124. doi: 10.1007/BF00120641

CrossRef Full Text | Google Scholar

Herve, M. (2018). RVAideMemoire: testing and plotting procedures for biostatistics. R Package Version 0.9-69. Available online at: https://CRAN.R-project.org/package=RVAideMemoire.

Google Scholar

Hewitt, G. M. (1996). Some genetic consequences of ice ages, and their role in divergence and speciation. Biol. J. Linnean Soc. 58, 247–276. doi: 10.1111/j.1095-8312.1996.tb01434.x

CrossRef Full Text | Google Scholar

Hewitt, G. M. (2004). Genetic consequences of climatic oscillations in the Quaternary. Phil. Trans. R. Soc. Lond. B 359, 183–195. doi: 10.1098/rstb.2003.1388

CrossRef Full Text | Google Scholar

Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. doi: 10.1002/joc.1276

CrossRef Full Text | Google Scholar

Hipkins, V. D., Krutovskii, K. V., Strauss, S. H. (1994). Organelle genomes in conifers: structure, evolution, and diversity. For. Genet. 1, 179–189.

Google Scholar

Hoffmann, A. A., Sgrò, C. M. (2011). Climate change and evolutionary adaptation. Nature 470, 479–485. doi: 10.1038/nature09670

PubMed Abstract | CrossRef Full Text | Google Scholar

Holderegger, R., Herrmann, D., Poncet, B., Gugerli, F., Thuiller, W., Taberlet, P., et al. (2008). Land ahead: using genome scans to identify molecular markers of adaptive relevance. Plant Ecol. Div. 1, 273–283. doi: 10.1080/17550870802338420

CrossRef Full Text | Google Scholar

Holsinger, K. E., Lewis, P. O. (2003). Hickory: a package for analysis of population genetic data v1.1, (Storrs, USA: Department of Ecology and Evolutionary Biology, University of Connecticut), Available online at: http://www.academia.edu/1839794/.

Google Scholar

Hsieh, Y.-C., Chung, J.-D., Wang, C.-N., Chang, C.-T., Chen, C.-Y., Hwang, S.-Y. (2013). Historical connectivity, contemporary isolation and local adaptation in a widespread but discontinuously distributed species endemic to Taiwan, Rhododendron oldhamii (Ericaceae). Heredity 111, 147–156. doi: 10.1038/hdy.2013.31

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C.-L., Chang, C.-T., Huang, B.-H., Chung, J.-D., Chen, J.-H., Chiang, Y.-C., et al. (2015a). Genetic relationships and ecological divergence in Salix species and populations in Taiwan. Tree Genet. Genom. 11, 39. doi: 10.1007/s11295-015-0862-1

CrossRef Full Text | Google Scholar

Huang, C.-L., Chen, J.-H., Tsang, M.-H., Chung, J.-D., Chang, C.-T., Hwang, S.-Y. (2015b). Influences of environmental and spatial factors on genetic and epigenetic variations in Rhododendron oldhamii (Ericaceae). Tree Genet. Genom. 11, 823. doi: 10.1007/s11295-014-0823-0

CrossRef Full Text | Google Scholar

Huete, A. R., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., Ferreira, L. G. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 83, 195–213. doi: 10.1016/S0034-4257(02)00096-2

CrossRef Full Text | Google Scholar

Hutchison, D. W., Templeton, A. R. (1999). Correlation of pairwise genetic and geographic distance measures: inferring the relative influences of gene flow and drift on the distribution of genetic variability. Evolution 53, 1898–1914. doi: 10.1111/j.1558-5646.1999.tb04571.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hwang, S.-Y., Lin, T.-P., Ma, C.-S., Lin, C.-L., Chung, J.-D., Yang, J.-C. (2003). Postglacial population growth of Cunninghamia konishii (Cupressaceae) inferred from phylogeographical and mismatch analysis of chloroplast DNA variation. Mol. Ecol. 12, 2689–2695. doi: 10.1046/j.1365-294X.2003.01935.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeffreys, H. (1961). Theory of probability. Oxford: Oxford University Press.

Google Scholar

Jombart, T., Ahmed, I. (2011). Adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinformatics 27, 3070–3071. doi: 10.1093/bioinformatics/btr521

PubMed Abstract | CrossRef Full Text | Google Scholar

Jombart, T., Devillard, S., Balloux, F. (2010). Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 11, 94. doi: 10.1186/1471-2156-11-94

PubMed Abstract | CrossRef Full Text | Google Scholar

Jump, A. S., Peñuelas, J. (2006). Genetic effects of chronic habitat fragmentation in a wind-pollinated tree. Proc. Natl. Acad. Sci. U.S.A. 103, 8096–8100. doi: 10.1073/pnas.0510127103

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalinowski, S. T. (2004). Counting alleles with rarefaction: private alleles and hierarchical sampling designs. Conserv. Genet. 5, 539–543. doi: 10.1023/B:COGE.0000041021.91777.1a

CrossRef Full Text | Google Scholar

Kamvar, Z. N., Tabima, J. F., Grünwald, N. J. (2014). Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. Peer J. 2, e281. doi: 10.7717/peerj.281

CrossRef Full Text | Google Scholar

Kamvar, Z. N., Brooks, J. C., Grünwald, N. J. (2015). Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality. Front. Genet. 6, 208. doi: 10.3389/fgene.2015.00208

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawecki, T. J., Ebert, D. (2004). Conceptual issues in local adaptation. Ecol. Lett. 7, 1225–1241. doi: 10.1111/j.1461-0248.2004.00684.x

CrossRef Full Text | Google Scholar

Kimura, M. (1953). Stepping stone model of population. Ann. Rep. Nat. Inst. Genet. Japan 3, 62–63.

Google Scholar

Kitagawa, R., Mimura, M., Mori, A. S., Sakai, A. (2015). Topographic patterns in the phylogenetic structure of temperate forests on steep mountainous terrain. AoB Plants 7, 1–134. doi: 10.1093/aobpla/plv134

CrossRef Full Text | Google Scholar

Lande, R. (1988). Genetics and demography in biological conservation. Science 241, 1455–1460. doi: 10.1126/science.3420403

PubMed Abstract | CrossRef Full Text | Google Scholar

Latta, R. G., Mitton, J. B. (1999). Historical separation and present gene flow through a zone of secondary contact in ponderosa pine. Evolution 53, 769–776. doi: 10.1111/j.1558-5646.1999.tb05371.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lenormand, T. (2002). Gene flow and the limits to natural selection. Trends Ecol. Evol. 17, 183–189. doi: 10.1016/S0169-5347(02)02497-7

CrossRef Full Text | Google Scholar

Lenth, R. (2018). Emmeans: estimated marginal means, aka least-squares means. R package version 1.3.1. Available online at: https://CRAN.R-project.org/package=emmeans.

Google Scholar

Li, C.-F., Chytrý, M., Zelený, D., Chen, M.-Y., Chen, T.-Y., Chiou, C.-R., et al. (2013). Classification of Taiwan forest vegetation. Appl. Veg. Sci. 16, 698–719. doi: 10.1111/avsc.12025

CrossRef Full Text | Google Scholar

Li, Y.-S., Chang, C.-T., Wang, C.-N., Thomas, P., Chung, J.-D., Hwang, S.-Y. (2018). The contribution of neutral and environmentally dependent processes in driving population and lineage divergence in Taiwania (Taiwania cryptomerioides). Front. Plant Sci. 9, 1148. doi: 10.3389/fpls.2018.01148

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, T.-P., Lu, C.-S., Chung, Y.-L., Yang, J.-C. (1993). Allozyme variation in four populations of Taiwania cryptomerioides in Taiwan. Silv. Genet. 42, 278–284.

Google Scholar

Lin, T.-P., Wang, C.-T., Yang, J.-C. (1998). Comparison of genetic diversity between Cunninghamia konishii and C. lanceolata. J. Hered. 89, 370–374. doi: 10.1093/jhered/89.4.370

CrossRef Full Text | Google Scholar

Linhart, Y. B., Grant, M. C. (1996). Evolutionary significance of local genetic differentiation in plants. Ann. Rev. Ecol. Syst. 27, 237–277. doi: 10.1146/annurev.ecolsys.27.1.237

CrossRef Full Text | Google Scholar

Liu, T. (1966). Study on the phytogeography of the conifers and taxads of Taiwan. Bull. Taiwan For. Res. Inst. 122, 1–33.

Google Scholar

Manni, F., Guerard, E., Heyer, E. (2004). Geographic patterns of (genetic, morphologic, linguistic) variation: how barriers can be detected by using Monmonier’s algorithm. Human Biol. 76, 173–190. doi: 10.1353/hub.2004.0034

CrossRef Full Text | Google Scholar

Mimura, M., Aitken, S. N. (2010). Local adaptation at the range peripheries of Sitka spruce. J. Evol. Biol. 23, 249–259. doi: 10.1111/j.1420-9101.2009.01910.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Monahan, W. B., Pereira, R. J., Wake, D. B. (2012). Ring distributions leading to species formation: a global topographic analysis of geographic barriers associated with ring species. BMC Biol. 10, 20. doi: 10.1186/1741-7007-10-20

PubMed Abstract | CrossRef Full Text | Google Scholar

Monmonier, M. S. (1973). Maximum-difference barriers: an alternative numerical regionalization method. Geograp. Anal. 3, 245–261. doi: 10.1111/j.1538-4632.1973.tb01011.x

CrossRef Full Text | Google Scholar

Naimi, B., Hamm, N. A. S., Groen, T. A., Skidmore, A. K., Toxopeus, A. G. (2014). Where is positional uncertainty a problem for species distribution modelling? Ecography 37, 191–203. doi: 10.1111/j.1600-0587.2013.00205.x

CrossRef Full Text | Google Scholar

Nakazato, T., Warren, D. L., Moyle, L. C. (2010). Ecological and geographic modes of species divergence in wild tomatoes. Am. J. Bot. 97, 680–693. doi: 10.3732/ajb.0900216

PubMed Abstract | CrossRef Full Text | Google Scholar

Neale, D. B., Marshall, K. A., Harry, D. E. (1991). Inheritance of chloroplast and mitochondrial DNA in incense-cedar (Calocedrus decurrens). Can. J. For. Res. 21, 717–720. doi: 10.1139/x91-100

CrossRef Full Text | Google Scholar

Nei, M. (1987). Molecular evolutionary genetics. New York, NY: Columbia University Press. doi: 10.7312/nei-92038

CrossRef Full Text | Google Scholar

Nga, N. T. T., Dung, N. A., Chung, N. T., Thai, T. H., Hung, N. D. (2016). The distribution and some ecological characteristics, and essential oil of Cunninghamia konishii Hayata in Pu Hoat nature reserve, Nghe An province, Vietnam. KKU Eng. J. 43, 121–124. doi: 10.14456/kkuenj.2016.37

CrossRef Full Text | Google Scholar

Nosil, P., Vines, T. H., Funk, D. J. (2005). Reproductive isolation caused by natural selection against immigrants from divergent habitats. Evolution 59, 705–719. doi: 10.1111/j.0014-3820.2005.tb01747.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Nybom, H. (2004). Comparison of different nuclear DNA markers for estimating intraspecific genetic diversity in plants. Mol. Ecol. 13, 1143–1155. doi: 10.1111/j.1365-294X.2004.02141.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2017). vegan: community ecology package. R Package Version 2.4-2. Available online at: https://CRAN.R-project.org/package=vegan.

Google Scholar

Parmesan, C., Yohe, G. (2003). A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42. doi: 10.1038/nature01286

PubMed Abstract | CrossRef Full Text | Google Scholar

Peres-Neto, P. R., Legendre, P., Dray, S., Borcard, D. (2006). Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology 87, 2614–2625. doi: 10.1890/0012-9658(2006)87[2614:VPOSDM]2.0.CO;2

PubMed Abstract | CrossRef Full Text | Google Scholar

Petit, R. J., Aguinagalde, I., de Beaulieu, J.-L., Bittkau, C., Brewer, S., Cheddadi, R., et al. (2003). Glacial refugia: hotspots but not melting pots of genetic diversity. Science 300, 1563–1565. doi: 10.1126/science.1083264

PubMed Abstract | CrossRef Full Text | Google Scholar

Petit, R. J., Hampe, A. (2006). Some evolutionary consequences of being a tree. Annu. Rev. Ecol. Evol. Syst. 37, 187–214. doi: 10.1146/annurev.ecolsys.37.091305.110215

CrossRef Full Text | Google Scholar

Pulliam, H. R. (2000). On the relationship between niche and distribution. Ecol. Lett. 3, 349–361. doi: 10.1046/j.1461-0248.2000.00143.x

CrossRef Full Text | Google Scholar

R Development Core Team. (2018). R: a language and environment for statistical computing. Version 3.5.2. Vienna: R Foundation for Statistical Computing. Available online at: https://www.r-project.org/.

Google Scholar

Ribeiro, M. M., Mariette, S., Vendramin, G. G., Szmidt, A. E., Plomion, C., Kremer, A. (2002). Comparison of genetic diversity estimates within and among populations of maritime pine using chloroplast simple-sequence repeat and amplified fragment length polymorphism data. Mol. Ecol. 11, 869–877. doi: 10.1046/j.1365-294X.2002.01490.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ringbauer, H., Kolesnikov, A., Field, D. L., Barton, N. H. (2018). Estimating barriers to gene flow from distorted isolation by distance patterns. Genetics 208, 1231–1245. doi: 10.1534/genetics.117.300638

PubMed Abstract | CrossRef Full Text | Google Scholar

Roscher, C., Schumacher, J., Gubsch, M., Lipowsky, A., Weigelt, A., Buchmann, N., et al. (2012). Using plant functional traits to explain diversity–productivity relationships. PLoS ONE 7, e36760. doi: 10.1371/journal.pone.0036760

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenberg, N. J., Blat, B. L., Verma, S. B., (1983). Microclimate: the biological environment. New York, NY: Wiley.

Google Scholar

Sexton, J. P., Hangartner, S. B., Hoffmann, A. A. (2014). Genetic isolation by environment or distance: which pattern of gene flow is most common? Evolution 68, 1–15. doi: 10.1111/evo.12258

PubMed Abstract | CrossRef Full Text | Google Scholar

Shih, K.-M., Chang, C.-T., Chung, J.-D., Chiang, Y.-C., Hwang, S.-Y. (2018). Adaptive genetic divergence despite significant isolation-by-distance in populations of Taiwan cow-tail fir (Keteleeria davidiana var. formosana). Front. Plant Sci. 9, 92. doi: 10.3389/fpls.2018.00092

PubMed Abstract | CrossRef Full Text | Google Scholar

Slatkin, M. (2008). Linkage disequilibrium—understanding the evolutionary past and mapping the medical future. Nat. Rev. Genet. 9, 477–485. doi: 10.1038/nrg2361

PubMed Abstract | CrossRef Full Text | Google Scholar

Szpiech, Z. A., Jakobsson, M., Rosenberg, N. A. (2008). ADZE: a rarefaction approach for counting alleles private to combinations of populations. Bioinformatics 24, 2498–2504. doi: 10.1093/bioinformatics/btn478

PubMed Abstract | CrossRef Full Text | Google Scholar

Storey, J. D., Bass, A. J., Dabney, A., Robinson, D. (2019). qvalue: Q-value estimation for false discovery rate control. R package version 2.14.1. Available online at: http://github.com/StoreyLab/qvalue.

Google Scholar

Strasburg, J. L., Sherman, N. A., Wright, K. M., Moyle, L. C., Willis, J. H., Rieseberg, L. H. (2012). What can patterns of differentiation across plant genomes tell us about adaptation and speciation? Philos. Trans. Biol. Sci. 367, 364–373. doi: 10.1098/rstb.2011.0199

CrossRef Full Text | Google Scholar

Stucki, S., Orozco-terWengel, P., Bruford, M. W., Colli, L., Masembe, C., Negrini, R., et al. (2017). High performance computation of landscape genomic models integrating local indices of spatial association. Mol. Ecol. Resour. 17, 1072–1089. doi: 10.1111/1755-0998.12629

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, H. J. (1984). Studies on the climate and vegetation types of the natural forest in Taiwan. (II). Altitudinal vegetation zones in relation to temperature gradient. Q. J. Chin. For. 17, 57–73.

Google Scholar

Thornthwaite, C. W. (1948). An approach toward a rational classification of climate. Geogr. Rev. 38, 55–94. doi: 10.2307/210739

CrossRef Full Text | Google Scholar

Turgeon, J., Bernatchez, L. (2001). Clinal variation at microsatellite loci reveals historical secondary intergradation between glacial races of Coregonus artedi (Teleostei: Coregoninae). Evolution 55, 2274–2286. doi: 10.1111/j.0014-3820.2001.tb00742.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Vandewalle, M., de Bello, F., Berg, M. P., Bolger, T., Dolédec, S., Dubs, F., et al. (2010). Functional traits as indicators of biodiversity response to land use changes across ecosystems and organisms. Biodivers. Conserv. 19, 2921–2947. doi: 10.1007/s10531-010-9798-9

CrossRef Full Text | Google Scholar

Vekemans, X., Beauwens, T., Lemaire, M., Roldán-Ruiz, I. (2002). Data from amplified fragment length polymorphism (AFLP) markers show indication of size homoplasy and of a relationship between degree of homoplasy and fragment size. Mol. Ecol. 11, 139–151. doi: 10.1046/j.0962-1083.2001.01415.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Via, S. (2009). Natural selection in action during speciation. Proc. Natl. Acad. Sci. U.S.A. 106, 9939–9946. doi: 10.1073/pnas.0901397106

PubMed Abstract | CrossRef Full Text | Google Scholar

Vos, P., Hogers, R., Bleeker, M., Reijans, M., van der Lee, T., Hornes, M., et al. (1995). AFLP: a new technique for DNA fingerprinting. Nucl. Acids Res. 23, 4407–4414. doi: 10.1093/nar/23.21.4407

CrossRef Full Text | Google Scholar

Wall, J. D., Andolfatto, P., Przeworski, M. (2002). Testing models of selection and demography in Drosophila simulans. Genetics 162, 203–216.

PubMed Abstract | Google Scholar

Wang, I. J., Glor, R. E., Losos, J. B. (2013). Quantifying the roles of ecology and geography in spatial genetic divergence. Ecol. Lett. 16, 175–182. doi: 10.1111/ele.12025

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Feng, C., Jiao, T., Von Wettberg, E. B., Kang, M. (2017). Genomic signatures of adaptive divergence despite strong nonadaptive forces on edaphic islands: a case study of Primulina juliae. Genom. Biol. Evol. 9, 3495–3508. doi: 10.1093/gbe/evx263

CrossRef Full Text | Google Scholar

Weir, B. S., Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population-structure. Evolution 38, 1358–1370. doi: 10.2307/2408641

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, J. B., Lindell, J., Backström, N. (2010). Speciation genetics: current status and evolving approaches. Philos. Trans. Biol. Sci. 365, 1717–1733. doi: 10.1098/rstb.2010.0023

CrossRef Full Text | Google Scholar

Zhang, J., Cheng, K., Zang, R., Ding, Y. (2014). Environmental filtering of species with different functional traits into plant assemblages across a tropical coniferous-broadleaved forest ecotone. Plant Soil 380, 361–374. doi: 10.1007/s11104-014-2088-7

CrossRef Full Text | Google Scholar

Zhang, H., Hare, M. P. (2012). Identifying and reducing AFLP genotyping error: an example of tradeoffs when comparing population structure in broadcast spawning versus brooding oysters. Heredity 108, 616– 625. doi: 10.1038/hdy.2011.132

CrossRef Full Text | Google Scholar

Zhivotovsky, L. A. (1999). Estimating population structure in diploids with multilocus dominant DNA markers. Mol. Ecol. 8, 907–913. doi: 10.1046/j.1365-294x.1999.00620.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: AFLP, barriers to gene flow, conservation, Cunninghamia konishii, mountain ranges

Citation: Li Y-S, Shih K-M, Chang C-T, Chung J-D and Hwang S-Y (2019) Testing the Effect of Mountain Ranges as a Physical Barrier to Current Gene Flow and Environmentally Dependent Adaptive Divergence in Cunninghamia konishii (Cupressaceae). Front. Genet. 10:742. doi: 10.3389/fgene.2019.00742

Received: 01 February 2019; Accepted: 15 July 2019;
Published: 09 August 2019.

Edited by:

Melanie April Murphy, University of Wyoming, United States

Reviewed by:

Roser Vilatersana, Spanish National Research Council (CSIC), Spain
Zhiyong Zhang, Jiangxi Agricultural University, China

Copyright © 2019 Li, Shih, Chang, Chung and Hwang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Shih-Ying Hwang, hsy9347@ntnu.edu.tw

These authors have contributed equally to this work.