Next Article in Journal
Efficacy of Sterculia diversifolia Leaf Extracts: Volatile Compounds, Antioxidant and Anti-Inflammatory Activity, and Green Synthesis of Potential Antibacterial Silver Nanoparticles
Next Article in Special Issue
Population Genetic Structure and Diversity of Cryptic Species of the Plant Genus Macrocarpaea (Gentianaceae) from the Tropical Andes
Previous Article in Journal
Research Progress on the Effects of Selenium on the Growth and Quality of Tea Plants
Previous Article in Special Issue
Genetic Structure of Native Blue Honeysuckle Populations in the Western and Eastern Eurasian Ranges
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Plants
Volume 11
Issue 19
10.3390/plants11192490
plants-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Pattern of Adaptive Divergence in Zingiber kawagoii Hayata (Zingiberaceae) along a Narrow Latitudinal Range

by
Yi-Shao Li
1,
Pei-Chun Liao
1,
Chung-Te Chang
2 and
Shih-Ying Hwang
1,*
1
School of Life Science, National Taiwan Normal University, 88 Tingchow Road, Section 4, Taipei 11677, Taiwan
2
Department of Life Science, Tunghai University, 1727 Taiwan Boulevard, Section 4, Taichung 40704, Taiwan
*
Author to whom correspondence should be addressed.
Plants 2022, 11(19), 2490; https://doi.org/10.3390/plants11192490
Submission received: 5 August 2022 / Revised: 16 September 2022 / Accepted: 20 September 2022 / Published: 23 September 2022
(This article belongs to the Special Issue Plant Molecular Evolution and Population Ecology)
Browse Figures
Versions Notes

Abstract

:
Ecological and evolutionary processes linking adaptation to environment are related to species’ range shifts. In this study, we employed amplified-fragment-length-polymorphism-based genome scan methods to identify candidate loci among Zingiber kawagoii populations inhabiting varying environments distributed at low to middle elevations (143–1488 m) in a narrow latitudinal range (between 21.90 and 25.30° N). Here, we show evidence of selection driving the divergence of Z. kawagoii. Twenty-six FST outliers were detected, which were significantly correlated with various environmental variables. The allele frequencies of nine FST outliers were either positively or negatively correlated with the population mean FST. Using several independent approaches, we found environmental variables act in a combinatorial fashion, best explaining outlier genetic variation. Nonetheless, we found that adaptive divergence was affected mostly by annual temperature range, and it is significantly positively correlated with latitude and significantly negatively correlated with the population mean FST. This study addresses a latitudinal pattern of changes in annual temperature range (which ranged from 13.8 °C in the Lanyu population to 18.5 °C in the Wulai population) and emphasizes the pattern of latitudinal population divergence closely linked to the allele frequencies of adaptive loci, acting in a narrow latitudinal range. Our results also indicate environmentally dependent local adaptation for both leading- and trailing-edge populations.

1. Introduction

The altitude-for-latitude model projected the scales of range retractions in altitude and latitude due to warming [1]. Locally adaptive alleles associated with environmental conditions in range shift margin populations are essential for species’ future resilience to climate change. Temperature and precipitation are the two most important climatic factors influencing the distribution, differentiation, and diversity of species [2]. Genetically based ecotypes may evolve corresponding to environmental changes and play a role in minimizing the extinction risk [3]. Genetic variation may vary along latitudinal and altitudinal clines [4]. Latitudinal environmental patterns are found to be correlated with intraspecific and interspecific diversity distributed in a large geographic scale [5]. Taiwan only covers a narrow latitudinal range of 385 km, lying between the north latitude 21.90 and 25.30, and the identification of latitudinal patterns in plant diversity and differentiation is important given the presence of ecological heterogeneities due to rugged topography and steep elevation in Taiwan.
Zingiber kawagoii is a species of herbaceous perennial plant in the Zingiberaceae family found in Taiwan and a small offshore island southeast of Taiwan. This species is endemic to Taiwan and is widely distributed along the west sides of the Hsuehshan Mountain Range and the Central Mountain Range in Taiwan, but it is sparsely distributed east of these mountain ranges, at low to middle elevations (140–1500 m), and from the southern to northern tip of Taiwan (Figure 1). Zingiber kawagoii was also found on a small island, Lanyu (orchid island), 90 km off the southeast coast of Taiwan. Populations of Z. kawagoii are found in different habitats including forest understories, forest edges, slopes, and valleys. Extremely low intraspecific variation and high population differentiation were found based on chloroplast DNA (cpDNA) variation [6]. A mean inbreeding coefficient of 0.656 was found, suggestive of a deficiency of heterozygotes that may have resulted from bottlenecks and/or inbreeding [7]. The latitudinal northerly expansion of Z. kawagoii after the last glacial maximum (LGM) has been inferred using ecological niche modeling [8]. This climate-induced northerly expansion may reduce genetic diversity and increase genetic as well as migration load [9,10], limiting the ability for adaptation and persistence in novel environments. However, locally adapted alleles may have been evoked during expansion, encountering novel selective regimes [4,11]. Although drift and migration both decrease local adaptation, smaller range-front populations may develop local adaptive divergence when selection is strong [12]. The leading-edge populations of northerly expanded Z. kawagoii may evolve locally adapted alleles in association with a leading-front environment. Nonetheless, trailing-edge populations are also important to the future survival of species [4,11,13].
Lower-latitude populations within species may show greater genetic divergence and evolutionary independence, contributing to reproductive isolation and speciation [5]. Moreover, adaptive genetic variation is widespread in herbaceous species with low levels of gene flow under strong selection pressures [14]. Because of its widespread distribution from the south to the north of Taiwan and the biogeographic history of latitudinal northerly expansion [8], the examination of the latitudinal pattern of Z. kawagoii population divergence related to environmental gradients will advance our understanding of species’ response to global change, particularly in a geographic area only spanning a narrow latitudinal range. Spatial heterogeneity in populations of Z. kawagoii may have played roles in driving local adaptation associated with environmental factors. Although genetic diversity in the latitudinal leading- and trailing-margin populations may be reduced, adaptive loci driven by natural selection may accumulate their frequencies at both range margins of distribution correlated strongly with environment [12,15,16]. No study has been conducted to test for local adaptation in Z. kawagoii. In this study, we surveyed the genetic variation in 212 individuals from 17 populations of Z. kawagoii using amplified fragment length polymorphisms (AFLPs) [17] and collected information regarding the environmental variables of sampling sites. The general aim of this study was to investigate the pattern of population divergence in Z. kawagoii related to adaptive evolution along a narrow latitudinal range. Specifically, we asked (1) how do population divergence and environmental variables vary along a narrow latitudinal range? (2) Which environmental variable(s) may have played the main role, and do environmental variables act in a combinatorial fashion in driving adaptive genetic divergence? (3) How do the allele frequencies of adaptive loci vary corresponding to the levels of population divergence?

2. Materials and Methods

2.1. Sampling, DNA Extraction, and AFLP Genotyping

Zingiber kawagoii individuals were collected from 17 populations (n = 212) which spanned a latitudinal range of 21.90–5.30° N and an altitudinal range of 143–1488 m (Figure 1). Within each population, samples were collected with a space at least 10 m apart. Because of latitudinal northerly expansion [7], populations distributed in higher latitudes can be recognized as leading-edge populations, and those distributed in the lower latitudes can be recognized as trailing-edge populations. Total genomic DNAs were extracted based on a cetyltrimethyl ammonium bromide (CTAB) procedure [18]. Total genomic DNA was ethanol precipitated and dissolved in 200 µL of TE buffer (pH 8.0). We quantified total genomic DNA using a NanoDrop spectrophotometer (NanoDrop Technology, Wilmington, DE, USA). In a total 10 µL reaction volume, 200 ng of total genomic DNA was mixed with 1 U EcoRI and 1 U MseI restriction enzymes incubated in 10X CutSmart buffer (New England Biolabs, Ipswich, MA, USA) at 37 °C for 1.5 h for restriction digestion. The reaction was deactivated at 65 °C for 15 min. The digested DNA products were ligated to AFLP adaptors (5 µM EcoRI and 50 µM MseI) with 5 U T4 DNA ligase (Thermo Scientific, Vilnius, Lithuania) and 5X ligation buffer (Thermo Scientific) at 22 °C for 1 h in a 10 µL ligation reaction mixture.
AFLP pre-selective amplification was performed using polymerase chain reaction (PCR). Adaptor ligated product (1:9 dilution with ddH2O) was mixed with 12 µM EcoRI (E00: 5′-GACTGCGTACCAATTC-3′) primer, 12 µM MseI (M00: 5′-GATGAGTCCTGAGTAA-3′) primer, 2.5 mM dNTPs, 1 U Taq DNA polymerase (Zymeset Biotech, Taipei, Taiwan), and 10X PCR buffer (Zymeset) in a 20 µL total volume. The PCR for pre-selective amplification was an initial holding at 72 °C for 2 min and pre-denaturation at 94 °C for 3 min, followed by 25 cycles of 30 s at 94 °C, 30 s at 56 °C, and 1 min at 72 °C, with a final 5 min holding at 72 °C. Either three or five additional bases were added at the ends of E00 and M00, and 90 primer combinations were screened initially for proper amplification, the number of fragments were amplified, and the genotyping error rate did not exceed 5%. Finally, 11 EcoRI-MseI (E00 and M00) primer combinations were used in selective amplification (Table S1). A labeled EcoRI selective primer (6-carboxyfluorescein or hexachloro-fluorescein dye labeled) was used in selective amplification. Selective amplification was performed in a 20 µL total volume containing 10 µM EcoRI and 10 µM MseI primers, 2.5 mM dNTPs, 2 U Taq DNA polymerase (Zymeset), 10X PCR buffer (Zymeset), and diluted pre-selective amplified products (1:19 dilution with ddH2O). Selective PCR was an initial holding at 94 °C for 3 min, followed by 13 cycles of 30 s at 94 °C, 30 s at 65–56 °C (decreasing the temperature by 0.7 °C each cycle), 1 min at 72 °C, then 23 cycles of 30 s at 94 °C, 30 s at 56 °C, and 1 min at 72 °C, with a final 5 min holding at 72 °C. Selective amplification products were electrophoresed on an ABI 3730XL DNA analyzer. We used Peak Scanner v.1.0 (Applied Biosystem, Foster City, CA, USA) to score amplification fragments. A fluorescent threshold set at 150 units was used to avoid background noise when scoring AFLP fragments in the range of 100–500 bp. We removed low peaks and fragments within one nucleotide in a ±0.8 base pair window which were recognized as the same fragment. Additionally, amplified fragments that scored higher than 99% or less than 1% of individuals were also removed. Three randomly chosen samples in each population were used to calculate the genotyping error rate per locus. Loci with an error rate per locus greater than 5% were removed [19]. The mean error rate per locus was 4.11% (Table S1).

2.2. Environmental Variables

Three categories of environmental variables were used in the study. Nineteen bioclimatic variables for sample sites at 30 s spatial resolution (~1 km) were downloaded from the WorldClim 1.4 for information related to temperature and precipitation [20]. Topographic variables, including aspect, elevation, and slope at 20 m resolution were obtained from Taiwan Geospatial One Stop, Ministry of the Interior (https://data.gov.tw/dataset/138563, accessed on 12 June 2021). Apart from the bioclimatic variables downloaded from WorldClim and three topographic variables, twelve other environmental variables not directly defined as bioclimatic and topographic variables were grouped as ecological variables. The twelve ecological variables were the normalized difference vegetation index (NDVI), the enhanced vegetation index (EVI), the leaf area index (LAI), the fraction of absorbed photosynthetically active radiation (fPAR), the relative humidity (RH), cloud cover (CLO), sunshine hours (SunH), the number of rainfall days per year (RainD), the mean wind speed (WSmean), the soil pH, the annual total potential evapotranspiration (PET), and the annual moisture index (MI).
Data from a moderate resolution imaging spectroradiometer (MODIS) recorded during 2001–2020 in the Land Process Distributed Active Archive Center (http://lpdaac.usgs.gov, accessed on 12 June 2021) were obtained for NDVI and EVI (dataset MOD13A2, 1 km resolution), LAI and fPAR (MOD15A2 dataset, 500 m resolution), and PET (MOD16A3 dataset, 500 m resolution). The monthly mean values of these variables were computed using a maximum-value composite procedure. The monthly mean values of RH, CLO, SunH, RainD, and WSmean at 1 km resolution were calculated using a universal spherical model of the Kriging method in ArcGIS with data obtained from the Data Bank for Atmospheric & Hydrologic Research (https://dbahr.pccu.edu.tw/, recorded in 1991–2020; accessed on 12 June 2021). Sampling site soil pH values were obtained from the soil data of an island-wide, 1,150 site investigation conducted in 1969–1986 by the Agriculture and Food Agency of Taiwan. Annual potential evapotranspiration values derived from the annual mean temperature and annual precipitation were used to calculate the annual MI.
Environmental variables of the three categories were separately used in the calculation of variance inflation factor (VIF) with a correlation threshold of |0.7| using the vifcor function of the R package usdm [21] in the R environment [22]. The bioclimatic variables BIO7 (annual temperature range), BIO9 (mean temperature of the driest quarter), BIO12 (annual precipitation), and BIO19 (precipitation of the coldest quarter); topographic variables aspect, elevation, and slope; and the ecological variables CLO, EVI, LAI, NDVI, MI, PET, RH, soil pH, and WSmean were retained based on VIF values smaller than 10. Moreover, the forward selection procedure was used for these variables in each environmental category separately, and the final set of environmental variables retained were aspect, BIO7, BIO12, NDVI, PET, RH, and WSmean (Table 1) if more than 5% of outlier genetic variation (adjusted R2 ≥ 0.05) was explained by the individual variable (see forward selection below, Table S2). Pearson’s correlation coefficients of pairwise comparisons and VIFs (all < 5) are reported in Table S3.

2.3. Genetic Diversity

AFLP-SURV v.1.0 [23] software was used to estimate the population unbiased expected heterozygosity (uHE) [24] and the proportion of polymorphic loci (%P, 95% criterion) based on allele frequencies using the settings of the Hardy–Weinberg equilibrium and non-uniform prior distribution. The per locus uHE was estimated using ARLEQUIN v.6.0 [25]. The index of association IA [26] and the modified index of association (rD) [27] are measures of multilocus linkage disequilibrium. These two measures were calculated using the ia function of the R package poppr [28]. A linear mixed effect model (LMM) was used to estimate the difference of the mean uHE per locus among and between populations. In LMMs, population and locus were used as a fixed factor and a random factor, respectively, and they were analyzed using the lmer function of the R package lme4 [29] based on the reduced maximum likelihood method. Significance tested using the Anova function of the R package car was based on type II Wald χ2 statistics [30]. Pairwise population comparisons of the mean uHE per locus with Tukey’s post hoc test were assessed using the lsmeans function of the R package emmeans [31].

2.4. Genetic Differentiation, Clustering, and Relationships

The analysis of molecular variance (AMOVA) was used to estimate the level of genetic differentiation between populations (ΦST) using the poppr.amova function of the R package poppr. Significance was tested using the randtest function of the R package ade4 [32] with 9999 permutations. The pairwise population FST was computed using ARLEQUIN, and the significance was tested with 10,000 permutations. Additionally, the level of divergence for each population from the remaining populations was calculated as the mean value of the pairwise FST for each population against the rest of the populations (denoted as the population mean FST). The population mean FST can be used as a proxy of the level of one population diverging genetically from the remaining populations.
Genetic homogeneous groups of individuals were assessed using the sNMF algorithm of landscape and ecological association (LEA) [33] and discriminant analysis of principal components (DAPC) [34]. A clustering scenario of K = 1–18 based on least-squares optimization was estimated using the snmf function of the R package LEA [33]. The parameters including regularization, iterations, and repetitions in snmf were set to 100, 200, and 10, respectively, and other arguments were set to defaults. The find.clusters and dapc functions of the R package adegenet [35] were used in DAPC analysis setting K = 1–10. The mean minimal cross-entropy (CE) in LEA and the Bayesian information criterion (BIC) in DAPC were estimated to determine the optimal number of clusters.
A neighbor-joining (NJ) tree was used to assess the genetic relationships among individuals. The nei.dist function of the R package poppr was used to calculate pairwise Nei’s genetic distances [36]. Nei’s genetic distance matrix was used to generate an unrooted NJ tree using the nj function of the R package ape [37], and bootstrap support values (BSP) were calculated based on 1000 replicates using the aboot function of the R package poppr.

2.5. Test for FST Outliers

To detect the signature of selection on AFLP loci, FST outliers were identified using DFDIST and BAYSESCAN. DFDIST is a modification for dominant markers of the software developed by Beaumont and Nichols [38]. DFDIST estimates a distribution of observed FST versus uHE, and loci under selection were identified by comparing them to a simulated neutral distribution. DFDIST parameters include critical frequency = 0.99; Zhivotovsky parameter = 0.25; trimmed mean FST = 0.3 (excluding 30% of highest and 30% of lowest FST values); smoothing proportion = 0.04; 500,000 resamplings; and critical p = 0.05. AFLP Loci with observed FST against uHE falling above the 95% confidence level of simulated distribution were recognized as FST outliers under directional selection. BAYESCAN v.2.1 [39] uses a reversible-jump Markov chain Monte Carlo algorithm to estimate the ratio of posterior probabilities of selection over neutrality (the posterior odds (POs)). Two hundred pilot runs of 50,000 iterations followed by a sample size of 50,000 with a thinning interval of 20 among 106 iterations were performed in BAYESCAN. Selection is detected when locus-specific component (α) is significantly different from zero. A positive α suggests divergent selection, while negative values suggest balancing or purifying selection. We used a criterion of a logarithmic scale of log10 (PO) > 2 as decisive evidence, corresponding to posterior probabilities between 0.99 and 1 [40], for selection over neutrality for a locus under directional selection (α > 0).

2.6. Test for Associations of AFLP Loci with Environmental Variables

To assess the associations of all genetic loci with environmental variables, the latent factor mixed model (LFMM) [41] and Samβada [42] were employed in testing for significant correlations of genetic variation in each AFLP locus with environmental variables. A latent random factor was incorporated in the hierarchical Bayesian mixed effect model implemented in the LFMM. Considering the background level of the population structure due to the demographic history and isolation-by-distance pattern, a number of latent factors of 3 was used according to the DAPC result (see Results), and a matrix of genetic variation was used as a fixed factor. For each environmental predictor, ten LFMM runs with 10,000 iterations of the Gibbs sampling algorithm and a burn-in period of 5000 cycles were performed. We obtained Z-scores for each environmental predictor by combining the results of ten independent LFMM runs and p values adjusted using the genomic inflation factor (λ) [41]. A false discovery rate (FDR)-adjusted p value of 1% was further applied using the qvalue function of the R package qvalue [43]. Samβada was used to assess the correlations of allele frequencies of AFLP loci with values of environmental variables based on the multiple univariate logistic regression approach. A 1% FDR for p value adjustment for both Wald and G scores was used to assess the fit of the model with environmental variables against the null model without environmental variables.
A Bayesian logistic regression analysis implemented in the stan_glm function of the R package rstanarm [44] was employed to further justify the associations of the potential FST outliers, identified using both BAYESCAN and DFDIST, with environmental variables. In stan_glm, the weakly informative priors following Student’s t distribution with a mean of zero and seven degrees of freedom were used, and the scale of the prior distribution was 10 for the intercept and 2.5 for the predictors. All stan_glm models were run with four chains for 2000 warm-up and 2000 sampling steps. The posterior_interval function of the R package rstanarm was employed to estimate 95% credible intervals for the determination of significant correlations of potential FST outliers with environmental variables. In stan_glm analysis, the effective sample size values representing overall sampling efficiencies for each predictor estimated were between 1632.8 and 7199.4, and the convergence diagnostic statistics were all close to 1. These values indicate good priors applied and stable estimates obtained.

2.7. Relative Contribution of Environmental Variables Explaining Variation in Potential FST Outliers

To test for the most important environmental variables explaining outlier genetic variation, we used the forward.sel function of the R package adespatial [45]. Forward selection was stopped if either the conventional level of significance (p < 0.05) or the global adjusted R2 was exceeded to prevent the overestimation of the explained variance, and significance was determined based on 999 permutations. Environmental variables in the three environmental categories were analyzed separately, and variables explaining more than 5% of outlier genetic variation (adjusted R2 > 0.05) were retained as the final set of environmental variables in the study (Table S2). Variables most importantly influencing outlier genetic variation were assessed using functions within the R package MuMIn [46]. Generalized linear models (GLMs) with a logit link function and a binomial residual distribution were used to assess the relationships between the outlier variation and the final set of retained environmental variables, and McFadden’s pseudo-R2 for the fixed effect of the best predicting model explaining outlier variation was calculated using the pR2 function of the R package pscl [47]. GLMs that fit all possible models for each outlier (response variable) were used in the dredge function and the subsequent model averaging analyses based on the Akaike information criterion with a correction for small sample sizes (AICc) (ΔAICc ≤ 2, the model.avg function). The AICc was used to rank the models and to calculate the Akaike sum of weights (SW) for each model [48]. The SW index, calculated using the importance function, was used to assess the relative importance of environmental variables contributing to explaining variation in the outlier loci. However, the SW index is arguably not an appropriate measure, representing the importance of model selection [49]. Therefore, a 95% confidence interval (CI) for each environmental variable included in the best predicting model for the variation in each outlier was estimated. The allele frequencies of FST outliers were used to test for correlation with the population mean FST.

2.8. Mantel Test and Variation Partitioning

The mantel test was used to analyze the correlations of the outlier AFLP Euclidean distance matrix with the Euclidean distance matrix of environments using the mantel function of the R package vegan [50] and the Euclidean distance matrix of environments controlling for latitudinal difference using the mantel.partial function. Environmental variables and outlier variation were used in a redundancy analysis (RDA). An RDA estimates the relative contribution of environmental variables explaining the outlier variation using the varpart function of the R package vegan. The total outlier variation was partitioned into four fractions attributable to (a) a pure environmental effect, (b) a geographically structured environmental effect, (c) a pure geographic effect, and (d) a residual effect [51]. We tested the significance of these fractions using the anova.cca function of the R package vegan with 999 permutations. Sample site geographic coordinates were used as geographic effects in variation partitioning.

3. Results

3.1. Genetic Diversity and Structure

A total of 621 AFLP (mean ± SD: 60.09 ± 13.96) loci was obtained using 11 selective amplification primer combinations (Table S1). The percentage of polymorphism varied from 20.7 (population EFS) to 48.7 (population HDD), with a mean of 37.0 (Table 2). The average level of uHE was 0.123, ranging from 0.102 in population TRK to 0.151 in population BTWS. A strong departure from random association between AFLP loci based on the measures of multilocus LD, IA, and rD was found for all populations examined (Table 2). The linear mixed effect model (LMM) analysis showed significant differences in the mean uHE per locus among populations (χ2 = 116.2, p < 0.001) and in many between-population comparisons (Table S4). Using the total data, population differentiation was high based on the AMOVA (ΦST = 0.308, p < 0.001) (Table S5) and pairwise FST (average FST = 0.298, p < 0.0001) (Table S6).

3.2. Genetic Clustering and Relationships

The mean minimal CE in LEA (Figure S1a) and the BIC in DAPC (Figure S1b) were minimized at K = 18 and K = 8, respectively. However, we observed changes in both the mean minimal CE and BIC elbowed at K = 3, which is consistent with three genetically homogeneous groups observed in LEA and DAPC (Figure 2). The three genetic clusters revealed by DAPC were: cluster A, containing populations LY and SL; cluster B, containing populations JSY, BTWS, and WLS; and cluster C, containing populations AT, EFS, HDD, JS, KTS, NZ, RF, SBS, SML, THS, TRK, and WL (Figure 2b). Genetic homogeneous grouping is concordant when comparing LEA to DAPC, despite the gene flow between populations observed in the LEA result (Figure 2a). Individuals of DAPC clusters A and B were grouped together in the NJ tree (Figure 3). Individuals of different populations of DAPC clusters A and B were clearly separated into different clades in the NJ tree, albeit with low BSPs. However, individuals of populations grouped in DAPC cluster C showed intermingled relationships in the NJ tree.

3.3. Latitudinal Trend of Annual Temperature Range and Population Mean FST

Pearson’s correlation test found a significant negative relationship between the population mean FST and latitude (Table S7, Figure 4a) and a significant positive relationship between the annual temperature range and latitude (Table S7, Figure 4b). Therefore, a moderate negative correlation between the annual temperature range and population mean FST was found (Figure 4c). Significant relationships between the population mean FST, latitude, and annual temperature range were also observed when the LY population was excluded from the analysis (population mean FST vs. latitude: Pearson’s r = −0.5842, p = 0.01748; population mean FST vs. annual temperature range: Pearson’s r = −0.5786, p = 0.0189; latitude vs. annual temperature range: Pearson’s r = 0.9341, p < 0.0001).

3.4. FST Outliers and Relative Importance of Environmental Variables Explaining Outlier Variation

BAYESCAN and DFDIST identified 26 loci (4.18%) as potential FST outliers (Table 3). All 26 FST outliers identified by FST-based methods were strongly correlated with environmental variables assessed using Samβada, the LFMM, and rstanarm (Table 3). We found very high population genetic differentiation based on the outlier data using the AMOVA (ΦST = 0.628, p < 0.001; Table S5). Additionally, the annual temperature range and annual precipitation were the two most important environmental variables influencing outlier genetic variation (adjusted R2 = 0.192 and adjusted R2 = 0.098, respectively; Table 4).
The relative importance of environmental variables estimated using model averaging of the most parsimonious models (ΔAICc ≤ 2) and the 95% CIs for coefficients of environmental covariates in the best predicting models revealed that environmental variables acted in a combinatorial fashion, influencing the genetic variations in the 26 FST outliers (Table 5). Although the annual temperature range was not necessarily included in the best predicting models explaining outlier variation (Table 5), it was the only environmental variable significantly correlated with latitude (Table S7, Figure 4b).
The 95% CIs indicated that no environmental variable significantly explained genetic variations in three outlier AFLP loci (P12_1612, P19_2619, and P19_2812) (Table 5), whereas the other 23 FST outliers were significantly explained by a combination of environmental variables. Moreover, environmental variables with high SW values may show no significant effects (95% CIs bracket zeros) on outlier variation. Nine outlier loci were found to have either significant positive (P05_2291, P08_2919, P21_3013, and P35_1635) or negative (P01_1888, P08_2566, P13_2177, P21_1772, and P21_1955) relationships of allele frequency with the population mean FST (Figure 5). However, we found no allele frequency correlation with latitude (Table S8).
The total explainable outlier genetic variation by the seven retained environmental variables was 47.4% based on the 26 outlier loci, of which 24.0% and 18.9% were, respectively, attributed to pure environmental and geographically structured environmental effects. However, only 4.5% was attributed purely to geographic effects (Table 6).

4. Discussion

4.1. Pattern of Adaptive Divergence along a Narrow Latitudinal Range

Biological species richness and speciation rate are higher toward the equator, which may have been related to the higher genetic divergence and greater evolutionary independence of populations within species [5]. Clinal variation may arise from neutral drift processes or adaptation linked to local environmentally associated genotypes [52]. This study found no correlation between genetic diversity (uHE) and latitude (Pearson’s r = 0.021, p = 0.936; Table S7). However, a strong negative correlation between the latitude and population mean FST (Pearson’s r = −0.677, p = 0.003) indicates higher levels of population genetic divergence in the low-latitude populations of Z. kawagoii (Figure 4a, Table S7). Additionally, a strong positive correlation of the annual temperature range with latitude (Pearson’s r = 0.946, p < 0.001; Figure 4b, Table S7) is consistent with an environmental gradient along latitude. Correlations between the latitude and population mean FST and between the latitude and annual temperature range are evidence for local adaptation.
The LY population, located on a small island situated off the southeast coast of Taiwan, had the lowest genetic distance to the population JSY in southern Taiwan (FST = 0.300; Table S6) but clustered with the geographically closest SL population in DPAC (Figure 1 and Figure 2b). Additionally, the LY population was closely related to other geographic proximity populations in southern Taiwan, including JSY, BTWS, and WLS based on the NJ tree (Figure 1 and Figure 3). These results suggest genetic connectivity between populations across the sea barrier, which was also found in Pemphis acidula [53] and Setaria viridis [54]. Additionally, the direction and significant relationships between the population mean FST, latitude, and annual temperature range held when the LY population was excluded based on Pearson’s correlation test (see results). Thus, we included the LY population as one of the populations in the Z. kawagoii latitudinal distribution.
Multilocus LD (IA and rD, Table 1) measures the non-random association of alleles [26,27]. A significant departure from zero of these measures may result from recent bottlenecks because of mating among genetically close individuals within populations [55], but it can also be influenced by mutation, recombination, natural selection, genetic drift, gene flow, and population size [56]. Our findings of environmentally dependent genetic variation (Table 3 and Table 5) and changes in the allele frequencies of adaptive loci strongly correlated with population divergence (Figure 5) suggest that significant IA and rD detected in all populations examined (Table 2) could be owing in part to natural selection [14,57,58]. In this study, the exceptionally high level of population differentiation analyzed using the 26 adaptive loci (ΦST = 0.628, p < 0.001; Table S5) suggests that environmentally based divergent selection may have played important roles in generating population adaptive divergence (the mantel test of outlier genetic distance matrix against environmental distance matrix: rM = 0.505, p = 0.001). We identified environmental factors (the seven retained environmental variables, Table 1) significantly correlated with outlier genetic variation controlling for the latitudinal effect (partial mantel test: rM = 0.464, p = 0.001), suggesting that multiple environmental factors impose as selective drivers for local adaptation (Table 3 and Table 5). The strong adaptive differentiation can be attributed to the complexity of environmental factors, causing micro-evolutionary differentiation between Z. kawagoii populations [12,14,52].

4.2. Latitudinal Cline of Annual Temperature Range Is the Major Selective Driver for Local Adaptation

Temperature and precipitation are the two most important selective drivers for local adaptation in plants commonly found to influence fitness-related traits and survival [2,59]. In this study, the annual mean temperature was excluded from the final set of environmental variables, and it was not significantly correlated with the population mean FST (Pearson’s r = 0.116, p = 0.659) or with the latitude (Pearson’s r = 0.133, p = 0.611). The annual temperature range with a higher adjusted R2 than other environmental variables (Table 4) may have played the main role in driving outlier genetic variation, resulting in a significant latitudinal pattern (Figure 4b). The differential combinatorial effects of environmental factors [60] in the best predicting models also played crucial roles in influencing genetic variations in the 26 FST outliers, with high pseudo R2 values (Table 5). The co-optimizing environmental variables may invoke locally adaptive genetic variation, particularly in rugged topographic landscapes such as Taiwan.
A strong linear fit of annual temperature range to a trend of reduction in the population mean FST along latitude (Figure 4) is consistent with Martin and McKay [5], who found that lower-latitude populations displayed greater genetic divergence. Local genotypes may be better adapted to local conditions, and adapted gene frequencies increase as natural selection persists overtime, which may generate clines in allele frequencies along the environmental gradient [10,61]. The nine potential FST outliers displayed strong correlations of allele frequencies with the population mean FST (Figure 5), suggesting the strong impact of selective pressures on population differentiation and the maintenance of adaptive integrity against the opposing forces of maladapted gene flow [10,62].
Adaptive divergence associated with thermal plasticity has been observed in natural populations of Plantago lanceolata [63] and Cynodon dactylon [64] along large-scale latitudinal gradients. However, a study finding annual thermal plasticity associated with population adaptive divergence along small-scale latitudinal gradients, to our knowledge, has not been documented. Thermal plasticity was thought to be more adaptive at higher latitudes due to the greater thermal variation in higher latitudes for species distributed in large geographic scales [65]. However, our results suggest adaptive evolution is evoked in response to different thermal ranges at higher- and lower-latitude populations (Figure 5).

4.3. Not Only Leading- but Also Trailing-Edge Populations Are Important for Zingiber kawagoii Conservation

In the current context of climate change, ecological, evolutionary, and conservation studies have demonstrated that populations at both limits of species distribution range evolved distinct genetic and phenotypic features [4,12]. Selection along thermal gradients of the environment can lead to the local adaptation and acclimatization of thermal-tolerance limits among populations [66,67]. The result of adaptive population differentiation might have related to range expansion toward higher latitudes (Figure 4 and Figure 5), and environmental boundaries between populations sharply shaped latitudinal cline in genetic divergence [10].
An initiation of shifting poleward in latitude and upward in elevation after the LGM and under the current global warming is expected [1]. The degree of temperature variability can affect the thermal-tolerance margins of organisms and is crucial to locally adapted responses to warming [66,67]. The current level of genetic diversity is a key determinant of a population adapting to changing environments [2,4,10,14]. High-diversity populations have broader stress-mitigation responses than low-diversity populations [68]. However, the genetic diversity of Z. kawagoii estimated using AFLP was lower (average uHE = 0.123) than that of the Brazilian Z. officinale (average uHE = 0.312) [69], three Indian Zingiber species (Z. neesanum: 0.240; Z. nimmonii: 0.164, and Z. zerumbet: 0.367) [70], and the diversity of thirteen plant species (average uHE = 0.230) [71]. The relatively low level of genetic diversity in Z. kawagoii was also reflected in the low average percentage of polymorphism (Table 2) and low cpDNA variation [6]. The low population genetic diversity in Z. kawagoii (Table 2) is probably due to factors such as the nature of inbreeding [7] and the isolation of populations [2,3], which may reduce the potential of evolving local adaptation [57,59].
Although low-latitude rear edge populations are expected to be small in size and are hence characterized by low genetic diversity [4], our data do not meet this rear edge hypothesis (Pearson’s correlation test between uHE and latitude: r = 0.021, p = 0.936). It is likely that geographic variation in thermal tolerance limits, consisting of both spatial temperature gradient and warming, can influence the rate of range shifts [4]. While cool margins are expanding, warm margin populations may persist locally due in part to local topographic and ecological conditions [1,4]. Moreover, two contrasting patterns of allele frequency change which correlated strongly with population divergence (Figure 5) indicate an increase in the probability of the presence of adaptive loci associated with the leading- and trailing-edge environments (Table 3 and Table 5). Apart from the nine loci (Figure 5), other outlier AFLP loci may persist at intermediate frequencies (Figure S2) for long periods due to heterogeneous selective pressures. Spatial range expansions can generate allele frequency gradients attributed to distinct selective processes [72]. This study suggests that locally adapted trailing- and leading-edge populations (Figure 5) are important for the future survival of species such as Z. kawagoii.

5. Conclusions

We studied the population mean FST along a latitudinal gradient that extended from the south of the Z. kawagoii distribution to its northern distributional margin. The annual temperature range related to the thermal tolerance of local populations could be the major environmental factor influencing outlier genetic variation. Additionally, the combination of various environmental variables may also be important to the local adaptation of Z. kawagoii. This study identified the presence of natural selection acting on adaptive loci at a small latitudinal scale, contributing to understanding how herbaceous species respond to environmental changes. The results broaden the generality of the latitudinal population divergence which is closely linked to environmental gradients in both the latitudinal leading- and trailing-edge populations of Z. kawagoii. Ecological speciation may occur in the low-latitude populations of Z. kawagoii, particularly because of high genetic divergence against other populations with narrower thermal tolerance limits.

Supplementary Materials

The following are available online at: https://www.mdpi.com/article/10.3390/plants11192490/s1, Figure S1: Minimum cross-entropy (a) and Bayesian information criteria (b) for evaluation of clustering scenarios, respectively, analyzed using LEA and DAPC. Figure S2: Heatmap of allele frequencies of the 26 outlier loci identified. The sequence of populations was arranged according to degree of latitude (°N). Table S1: Primer combinations, number of markers, and error rate per locus in AFLP for investigation in Zingiber kawagoii. Table S2: Relative contribution (adjusted R2) and F test of environmental variables explaining outlier genetic variation in Zingiber kawagoii using a forward selection procedure. Table S3: Variance inflation factor (VIF) of the seven environmental variables and Pearson’s correlation coefficients between these variables. Table S4: Summary of Tukey’s post hoc pairwise population comparisons of the mean unbiased expected heterozygosity (uHE) per locus using a linear mixed effect model. In linear mixed effect model, population was treated as a fixed factor and locus as a random factor based on the total AFLP variation of Zingiber kawagoii populations. Table S5: Genetic differentiation between populations within species of the 17 populations of Zingiber kawagoii based on the total and outlier AFLP variation using analysis of molecular variance (AMOVA). Table S6: Pairwise FST between populations of Zingiber kawagoii based on the total AFLP data using ARLEQUIN with 10,000 permutations. All pairwise comparisons were found to be significant (p < 0.0001). See Table 1 for population code. Table S7: Summary of the results of Pearson’s correlation test of population mean FST and seven environmental variables against population latitude. Table S8: Summary of the results of Pearson’s correlation test of allele frequencies of the 26 outlier loci against population mean FST and against population latitude.

Author Contributions

Conceptualization, S.-Y.H.; methodology, S.-Y.H. and Y.-S.L.; filed work, P.-C.L.; formal analysis, S.-Y.H., Y.-S.L. and C.-T.C.; writing—original draft, S.-Y.H. and Y.-S.L.; writing—review and editing, S.-Y.H., Y.-S.L., C.-T.C. and P.-C.L.; acquisition of financing, S.-Y.H. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by the National Science and Technology Council grant project MOST109-2313-B-003-001-MY3 to S.-Y.H.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

The authors are grateful to Bing-Hong Huang and Yi-Ting Tseng for their assistance in plant material collection and identification.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Jump, A.S.; Mátyás, C.; Peñuelas, J. The altitude-for-latitude disparity in range retractions of woody species. Trends Ecol. Evol. 2009, 24, 694–701. [Google Scholar] [CrossRef]
  2. Franks, S.J.; Hoffmann, A.A. Genetics of climate change adaptation. Annu. Rev. Genet. 2012, 46, 185–208. [Google Scholar] [CrossRef]
  3. Booy, G.; Hendriks, R.J.J.; Smulders, M.J.M.; van Groenendael, J.M.; Vosman, B. Genetic diversity and the survival of populations. Plant Biol. 2000, 2, 379–395. [Google Scholar] [CrossRef]
  4. Hampe, A.; Petit, R.J. Conserving biodiversity under climate change: The rear edge matters. Ecol. Lett. 2005, 8, 461–467. [Google Scholar] [CrossRef]
  5. Martin, P.R.; McKay, J.K. Latitudinal variation in genetic divergence of populations and the potential for future speciation. Evolution 2004, 58, 938–945. [Google Scholar] [CrossRef]
  6. Huang, B.-H.; Lin, Y.-C.; Huang, C.-W.; Lu, H.-P.; Luo, M.-X.; Liao, P.-C. Differential genetic responses to the stress revealed the mutation-order adaptive divergence between two sympatric ginger species. BMC Genom. 2018, 19, 692. [Google Scholar] [CrossRef]
  7. Tseng, Y.-T.; Luo, M.-X.; Huang, B.-H.; Liao, P.-C. Development of transferable expressed sequence tag-simple sequence repeat (EST-SSR) markers for delimitating two recently diverged endemic to Taiwan. Taiwania 2019, 64, 209–216. [Google Scholar] [CrossRef]
  8. Lu, H.-P. Contrasting Pattern of Population Genetic Structure and Local Adaptation of two Taiwan Endemic Ginger Species. Master’s Thesis, National Taiwan Normal University, Taipei, Taiwan, 2020. [Google Scholar]
  9. Kirkpatrick, M.; Barton, N.H. Evolution of a species’ range. Am. Nat. 1997, 150, 1–23. [Google Scholar] [CrossRef]
  10. Polechová, J.; Barton, N.H. Limits to adaptation along environmental gradients. Proc. Natl. Acad. Sci. USA 2015, 112, 6401–6406. [Google Scholar] [CrossRef]
  11. Cao, J.-J.; Li, Y.-S.; Chang, C.-T.; Chung, J.-D.; Hwang, S.-Y. Adaptive divergence without distinct species relationships indicate early stage ecological speciation in species of the Rhododendron pseudochrysanthum complex endemic to Taiwan. Plants 2022, 11, 1226. [Google Scholar] [CrossRef]
  12. Hämälä, T.; Mattila, T.M.; Savolainen, O. Local adaptation and ecological differentiation under selection, migration, and drift in Arabidopsis lyrata. Evolution 2018, 72, 1373–1386. [Google Scholar] [CrossRef]
  13. Bennett, S.; Duarte, C.M.; Marbà, N.; Wernberg, T. Integrating within-species variation in thermal physiology into climate change ecology. Phil. Trans. R. Soc. B 2019, 374, 20180550. [Google Scholar] [CrossRef]
  14. Briggs, D.; Walters, S.M. Plant Variation and Evolution, 4th ed.; Cambridge University Press: Cambridge, UK, 2016; ISBN 9781139060196. [Google Scholar] [CrossRef]
  15. Chien, W.-M.; Chang, C.-T.; Chiang, Y.-C.; Hwang, S.-Y. Ecological factors generally not altitude related played main roles in driving potential adaptive evolution at elevational range margin populations of Taiwan incense cedar (Calocedrus formosana). Front. Genet. 2020, 11, 580630. [Google Scholar] [CrossRef]
  16. Wessely, J.; Gattringer, A.; Guillaume, F.; Hulber, K.; Klonner, G.; Moser, D.; Dullinger, S. Climate warming may increase the frequency of cold-adapted haplotypes in alpine plants. Nat. Clim. Change 2022, 12, 77–82. [Google Scholar] [CrossRef]
  17. Vos, P.; Hogers, R.; Bleeker, M.; Reijans, M.; van der Lee, T.; Hornes, M.; Friters, A.; Pot, J.; Paleman, J.; Kuiper, M.; et al. AFLP: A new technique for DNA fingerprinting. Nucl. Acids Res. 1995, 23, 4407–4414. [Google Scholar] [CrossRef]
  18. Doyle, J.J.; Doyle, J.L. A rapid DNA isolation procedure from small quantities of fresh leaf material. Phytochem. Bull. 1987, 19, 11–15. [Google Scholar]
  19. Bonin, A.; Bellemain, E.; Bronken, E.P.; Pompanon, F.; Brochmann, C.; Taberlet, P. How to track and assess genotyping errors in population genetics studies. Mol. Ecol. 2004, 13, 3261–3273. [Google Scholar] [CrossRef]
  20. Hijmans, R.J.; Cameron, S.E.; Parra, J.L.; Jones, P.G.; Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 2005, 25, 1965–1978. [Google Scholar] [CrossRef]
  21. Naimi, B.; Hamm, N.A.S.; Groen, T.A.; Skidmore, A.K.; Toxopeus, A.G. Where is positional uncertainty a problem for species distribution modelling? Ecography 2014, 37, 191–203. [Google Scholar] [CrossRef]
  22. R Development Core Team. R: A Language and Environment for Statistical Computing, Version 3.6.3; R Foundation for Statistical Computing: Vienna, Austria, 2020; Available online: https://www.gbif.org/tool/81287/r-a-language-and-environment-for-statisticalcomputing (accessed on 6 January 2022).
  23. Vekemans, X.; Beauwens, T.; Lemaire, M.; Roldán-Ruiz, I. 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. 2002, 11, 139–151. [Google Scholar] [CrossRef]
  24. Nei, M. Molecular Evolutionary Genetics; Columbia University Press: New York, NY, USA, 1987. [Google Scholar] [CrossRef]
  25. Excoffier, L.; Lischer, H.E. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 2010, 10, 564–567. [Google Scholar] [CrossRef] [PubMed]
  26. Brown, A.H.; Feldman, M.W.; Nevo, E. Multilocus structure of natural populations of Hordeum spontaneum. Genetics 1980, 96, 523–536. [Google Scholar] [CrossRef]
  27. Agapow, P.M.; Burt, A. Indices of multilocus linkage disequilibrium. Mol. Ecol. Notes 2001, 1, 101–102. [Google Scholar] [CrossRef]
  28. Kamvar, Z.N.; Brooks, J.C.; Grünwald, N.J. Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality. Front. Genet. 2015, 6, 208. [Google Scholar] [CrossRef] [PubMed]
  29. Bates, D.; Maechler, M.; Bolker, B.; Walker, S. Fitting linear mixed effects models using lme4. J. Stat. Soft. 2015, 67, 1–48. [Google Scholar] [CrossRef]
  30. Fox, J.; Weisberg, S. An R Companion to Applied Regression, 2nd ed.; Sage: Thousand Oaks, CA, USA, 2011; ISBN 9781412975188. Available online: http://socserv.socsci.mcmaster.ca/jfox/Books/Companion (accessed on 11 March 2019).
  31. Lenth, R. Emmeans: Estimated Marginal Means, aka Least-Squares Means. R Package Version 1.3.1. Available online: https://CRAN.R-project.org/package=emmeans (accessed on 7 April 2018).
  32. Dray, S.; Dufour, A.-B. The ade4 package: Implementing the duality diagram for ecologists. J. Stat. Soft. 2007, 22, 1–20. [Google Scholar] [CrossRef]
  33. Frichot, E.; Francois, O. LEA: An R package for landscape and ecological association studies. Methods Ecol. Evol. 2015, 6, 925–929. [Google Scholar] [CrossRef]
  34. Jombart, T.; Devillard, S.; Balloux, F. Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. BMC Genet. 2010, 11, 94. [Google Scholar] [CrossRef]
  35. Jombart, T.; Ahmed, I. Adegenet 1.3-1: New tools for the analysis of genome-wide SNP data. Bioinformatics 2011, 27, 3070–3071. [Google Scholar] [CrossRef] [Green Version]
  36. Nei, M. Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics 1978, 89, 583–590. [Google Scholar] [CrossRef]
  37. Paradis, E.; Schliep, K. ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 2019, 35, 526–528. [Google Scholar] [CrossRef]
  38. Beaumont, M.A.; Nichols, R.A. Evaluating loci for use in the genetic analysis of population structure. Proc. R. Soc. B Biol. Sci. 1996, 263, 1619–1626. [Google Scholar] [CrossRef]
  39. Foll, M.; Gaggiotti, O. A genome scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective. Genetics 2008, 180, 977–993. [Google Scholar] [CrossRef]
  40. Foll, M. Bayescan 2.1 User Manual. Available online: http://cmpg.unibe.ch/software/BayeScan/files/BayeScan2.1_manual.pdf (accessed on 23 June 2012).
  41. Frichot, E.; Schoville, S.D.; Bouchard, G.; François, O. Testing for associations between loci and environmental gradients using latent factor mixed models. Mol. Biol. Evol. 2013, 30, 1687–1699. [Google Scholar] [CrossRef]
  42. Stucki, S.; Orozco-terWengel, P.; Forester, B.R.; Duruz, S.; Colli, L.; Masembe, C.; Negrini, R.; Landguth, E.; Jones, M.R.; The NEXTGEN Consortium; et al. High performance computation of landscape genomic models integrating local indices of spatial association. Mol. Ecol. Resour. 2017, 17, 1072–1089. [Google Scholar] [CrossRef]
  43. Storey, J.D.; Bass, A.J.; Dabney, A.; Robinson, D. Qvalue: Q-Value Estimation for False Discovery Rate Control. R Package Version 2.14.1. Available online: http://github.com/StoreyLab/qvalue (accessed on 8 January 2019).
  44. Goodrich, B.; Gabry, J.; Ali, I.; Brilleman, S. Rstanarm: Bayesian Applied Regression Modeling via Stan. R Package Version 2.17.4. Available online: http://mc-stan.org/ (accessed on 3 April 2018).
  45. Dray, S.; Blanchet, G.; Borcard, D.; Guenard, G.; Jombart, T.; Larocque, G.; Legendre, P.; Madi, N.; Wagner, H.H.; Dray, M.S. Adespatial: Multivariate Multiscale Spatial Analysis. R Package Version 0.3–7. Available online: https://cran.r-project.org/web/packages/adespatial/index.html (accessed on 21 February 2020).
  46. Barton, K. MuMIn: Multi-Model Inference. R Package Version 1.40.4. Available online: https://CRAN.R-project.org/package=MuMIn (accessed on 15 January 2018).
  47. Jackman, S. pscl: Classes and Methods for R Developed in the Political Science Computational Laboratory. R Package Version 1.5.5. Available online: https://github.com/atahk/pscl/ (accessed on 11 January 2021).
  48. Burnham, K.P.; Anderson, D.R. Model Selection and Multimodel Inference, 2nd ed.; Springer-Verlag: New York, NY, USA, 2002; ISBN 9781441929730. [Google Scholar] [CrossRef]
  49. Galipaud, M.; Gillingham, M.A.F.; Dechaume-Moncharmont, F.-X. A farewell to the sum of weights: The benefits of alternative metrics for variable importance estimations in model selection. Methods Ecol. Evol. 2017, 8, 1668–1678. [Google Scholar] [CrossRef]
  50. Oksanen, J.; Blanchet, F.G.; Friendly, M.; Kindt, R.; Legendre, P.; McGlinn, D.; Minchin, P.R.; O’Hara, R.B.; Simpson, G.L.; Solymos, P.; et al. Vegan: Community Ecology Package. R Package Version 2.4–2. Available online: https://CRAN.R-project.org/package=vegan (accessed on 15 January 2018).
  51. Borcard, D.; Legendre, P. All-scale spatial analysis of ecological data by means of principal coordinates of neighbor matrices. Ecol. Model. 2002, 153, 51–68. [Google Scholar] [CrossRef]
  52. Räsänen, K.; Hendry, A.P. Disentangling interactions between adaptive divergence and gene flow when ecology drives diversification. Ecol. Lett. 2008, 11, 624–636. [Google Scholar] [CrossRef]
  53. Deng, S.-L.; Lu, F.-Y.; Sen, Y.-C.; Cheng, C.-Y.; Chang, K.-C. Population genetic variation of Pemphis acidula among Taiwan and its nearby islands. Quar. J. Chin. Forest. 2008, 41, 149–164. [Google Scholar] [CrossRef]
  54. Hsieh, W.-H.; Chen, Y.-C.; Liao, H.-C.; Lin, Y.-R.; Chen, C.-H. High differentiation among populations of green foxtail, Setaria viridis, in Taiwan and adjacent islands revealed by microsatellite markers. Diversity 2021, 13, 159. [Google Scholar] [CrossRef]
  55. Wall, J.D.; Andolfatto, P.; Przeworski, M. Testing models of selection and demography in Drosophila simulans. Genetics 2002, 162, 203–216. [Google Scholar] [CrossRef] [PubMed]
  56. Slatkin, M. Linkage disequilibrium—Understanding the evolutionary past and mapping the medical future. Nat. Rev. Genet. 2008, 9, 477–485. [Google Scholar] [CrossRef] [PubMed]
  57. Barrett, R.D.H.; Schluter, D. Adaptation from standing genetic variation. Trends Ecol. Evol. 2008, 23, 38–44. [Google Scholar] [CrossRef] [PubMed]
  58. Eckert, A.J.; Bower, A.D.; González-Martínez, S.C.; Wegrzyn, J.L.; Coop, G.; Neale, D.B. Back to nature: Ecological genomics of loblolly pine (Pinus taeda, Pinaceae). Mol. Ecol. 2010, 19, 3789–3805. [Google Scholar] [CrossRef]
  59. Hoffmann, A.A.; Sgrò, C.M. Climate change and evolutionary adaptation. Nature 2011, 470, 479–485. [Google Scholar] [CrossRef]
  60. Linhart, Y.B.; Grant, M.C. Evolutionary significance of local genetic differentiation in plants. Ann. Rev. Ecol. Syst. 1996, 27, 237–277. [Google Scholar] [CrossRef]
  61. Barton, N.H. Clines in polygenic traits. Genet. Res. 1999, 74, 223–236. [Google Scholar] [CrossRef]
  62. Holliday, J.A.; Suren, H.; Aitken, S.N. Divergent selection and heterogeneous migration rates across the range of Sitka spruce (Picea sitchensis). Proc. R. Soc. B Biol. Sci. 2012, 279, 1675–1683. [Google Scholar] [CrossRef]
  63. Marshall, M.M.; Batten, L.C.; Remington, D.L.; Lacey, E.P. Natural selection contributes to geographic patterns of thermal plasticity in Plantago lanceolate. Mol. Ecol. 2019, 9, 2945–2963. [Google Scholar] [CrossRef]
  64. Zhang, J.-X.; Chen, M.-H.; Gan, L.; Zhang, C.-J.; Shen, Y.; Qian, J.; Han, M.-L.; Guo, Y.-Z.; Yan, X.-B. Diversity patterns of Bermuda grass along latitudinal gradients at different temperatures in Southeastern China. Plants 2020, 9, 1778. [Google Scholar] [CrossRef]
  65. Ghalambor, C.K.; Huey, R.B.; Martin, P.R.; Tewksbury, J.J.; Wang, G. Are mountain passess higher in the tropics? Janzen’s hypothesis revisited. Integr. Comp. Biol. 2006, 46, 5–17. [Google Scholar] [CrossRef]
  66. Richardson, J.L.; Urban, M.C.; Bolnick, D.I.; Skelly, D.K. Microgeographic adaptation and the spatial scale of evolution. Trends Ecol. Evol. 2014, 29, 165–176. [Google Scholar] [CrossRef]
  67. Sunday, J.M.; Bates, A.E.; Dulvy, N.K. Global analysis of thermal tolerance and latitude in ectotherms. Proc. R. Soc. B Biol. Sci. 2011, 278, 1823–1830. [Google Scholar] [CrossRef]
  68. Reusch, T.B.H.; Ehlers, A.; Hammerli, A.; Worm, B. Ecosystem recovery after climatic extremes enhanced by genotypic diversity. Proc. Natl. Acad. Sci. USA 2005, 102, 2826–2831. [Google Scholar] [CrossRef]
  69. Blanco, E.Z.; Bajak, M.M.; Siqueira, M.V.B.M.; Zucchi, M.I.; Pinheiro, J.B. Genetic diversity and structure of Brazilian ginger germplasm (Zingiber officinale) revealed by AFLP markers. Genetica 2016, 144, 627–638. [Google Scholar] [CrossRef]
  70. Thomas, G.E.; Geetha, K.A.; Augustine, L.; Mamiyil, S.; Thomas, G. Analyses between reproductive behavior, genetic diversity and Pythium responsiveness in Zingiber spp. reveal an adaptive significance for hemiclonality. Front. Plant Sci. 2016, 7, 1913. [Google Scholar] [CrossRef]
  71. Nybom, H. Comparison of different nuclear DNA markers for estimating intraspecific genetic diversity in plants. Mol. Ecol. 2004, 13, 1143–1155. [Google Scholar] [CrossRef]
  72. Excoffier, L.; Foll, M.; Petit, R.J. Genetic consequences of range expansions. Annu. Rev. Ecol. Syst. 2009, 40, 481–501. [Google Scholar] [CrossRef]
Plants 11 02490 g001 550
Figure 1. Sampling localities of the 17 Zingiber kawagoii populations. The coordinates of sampling sites were used to plot population locations using Tools in ArcGIS v.10.8.1. Map was derived from the default map database in ArcGIS, and the 20 m digital elevation model was used in the generation of elevation gradients. See Table 1 for abbreviations of the population names.
Figure 1. Sampling localities of the 17 Zingiber kawagoii populations. The coordinates of sampling sites were used to plot population locations using Tools in ArcGIS v.10.8.1. Map was derived from the default map database in ArcGIS, and the 20 m digital elevation model was used in the generation of elevation gradients. See Table 1 for abbreviations of the population names.
Plants 11 02490 g001
Plants 11 02490 g002 550
Figure 2. Analysis of genetic homogeneous groups of 212 individuals of Zingiber kawagoii based on the total AFLP variation using LEA (a) and DAPC (b). The clustering scenarios for K = 2–4 were displayed in LEA. The two linear discriminants LD1 and LD2 of DAPC described 52.72% and 14.84% of the total AFLP variation, respectively.
Figure 2. Analysis of genetic homogeneous groups of 212 individuals of Zingiber kawagoii based on the total AFLP variation using LEA (a) and DAPC (b). The clustering scenarios for K = 2–4 were displayed in LEA. The two linear discriminants LD1 and LD2 of DAPC described 52.72% and 14.84% of the total AFLP variation, respectively.
Plants 11 02490 g002
Plants 11 02490 g003 550
Figure 3. Neighbor-joining tree of 212 individuals of Zingiber kawagoii based on Nei’s genetic distances calculated using the total AFLP variation. Branch tip labels for individuals of different populations are colored differently. For each node, bootstrap support values greater than 70%, between 50% and 70%, and smaller than 50% are coded with green, red, and blue, respectively.
Figure 3. Neighbor-joining tree of 212 individuals of Zingiber kawagoii based on Nei’s genetic distances calculated using the total AFLP variation. Branch tip labels for individuals of different populations are colored differently. For each node, bootstrap support values greater than 70%, between 50% and 70%, and smaller than 50% are coded with green, red, and blue, respectively.
Plants 11 02490 g003
Plants 11 02490 g004 550
Figure 4. Regression plots showing the relationships between population mean FST and latitude (a), between annual temperature range and latitude (b), and between annual mean temperature and population mean FST (c). Pairwise population FST was estimated using the total AFLP variation and used in calculation of population mean FST (Table S6) BIO7, annual temperature range.
Figure 4. Regression plots showing the relationships between population mean FST and latitude (a), between annual temperature range and latitude (b), and between annual mean temperature and population mean FST (c). Pairwise population FST was estimated using the total AFLP variation and used in calculation of population mean FST (Table S6) BIO7, annual temperature range.
Plants 11 02490 g004
Plants 11 02490 g005 550
Figure 5. Linear regression plots of nine FST outliers showing significant correlation relationships of allele frequency with population mean FST. Pearson’s correlation test results are reported in Table S8.
Figure 5. Linear regression plots of nine FST outliers showing significant correlation relationships of allele frequency with population mean FST. Pearson’s correlation test results are reported in Table S8.
Plants 11 02490 g005
Table
Table 1. The seven retained site environmental variables of the 17 populations of Zingiber kawagoii.
Table 1. The seven retained site environmental variables of the 17 populations of Zingiber kawagoii.
PopulationAspectBIO7BIO12 NDVIPETRHWSmean
Antong (AT)288.0515.619330.841380.8678.213.15
Beitawushan (BTWS)120.1114.946160.771497.8376.682.69
Erfenshan (EFS)314.4118.124940.851634.2278.232.73
Huangdidian (HDD)299.1618.434810.841419.2278.882.57
Jianshi (JS)164.9517.325390.861436.6478.912.51
Jinshuiying (JSY)80.9714.347490.801449.9675.922.65
Kantoushan (KTS)290.6917.031200.781622.8678.692.69
Lanyu (LY)235.5813.827600.771379.8787.587.45
Nanzhuang (NZ)269.4118.125640.841489.5178.672.60
Ruifang (RF)250.4318.432820.771398.8478.292.83
Shibishan (SBS)81.6616.027260.771855.0881.432.24
Shuangliu (SL)63.6814.431000.861830.6876.152.96
Sunmoonlake (SML)256.6516.522620.801757.0081.121.28
Tahsueshan (THS)323.0517.125690.811632.2078.312.54
Taroko (TRK)294.9316.522920.841453.2978.782.84
Wulai (WL)105.6718.532310.781477.3078.822.43
Weiliaoshan (WLS)358.4016.130930.841769.8877.342.78
Aspect (0–360°). BIO7, annual temperature range (°C); BIO12, annual precipitation (mm); NDVI, normalized difference vegetation index (unitless); PET, annual total potential evapotranspiration (kg/m2/year); RH, relative humidity (%); WSmean, mean wind speed (m/s).
Table
Table 2. Site properties and genetic parameters of the 17 sampled populations of Zingiber kawagoii estimated based on the total AFLP variation.
Table 2. Site properties and genetic parameters of the 17 sampled populations of Zingiber kawagoii estimated based on the total AFLP variation.
PopulationLatitude
Longitude		Altitude (m)N%PuHE
(SE)		IA
(p)		rD
(p)
Antong (AT)23.2847
121.3721		6101435.20.113
(0.006)		3.614
(0.001)		0.016
(0.001)
Beitawushan (BTWS)22.6148
120.7022		11921240.70.151
(0.007)		2.040
(0.001)		0.008
(0.001)
Erfenshan (EFS)24.3919
120.8240		7691220.70.115
(0.007)		1.753
(0.001)		0.009
(0.001)
Huangdidian (HDD)24.9894
121.6799		4321048.70.143
(0.007)		2.022
(0.001)		0.009
(0.001)
Jianshi (JS)24.7307
121.2895		8501332.20.105
(0.006)		3.914
(0.001)		0.023
(0.001)
Jinshuiying (JSY)22.4075
120.7564		14881436.60.126
(0.006)		1.137
(0.001)		0.005
(0.001)
Kantoushan (KTS)23.2671
120.5010		5831435.40.115
(0.006)		4.421
(0.001)		0.022
(0.001)
Lanyu (LY)22.0496
121.5257		3021333.90.123
(0.007)		6.455
(0.001)		0.031
(0.001)
Nanzhuang (NZ)24.5742
121.0436		4671145.50.126
(0.006)		6.256
(0.001)		0.029
(0.001)
Ruifang (RF)25.0861
121.8385		3491143.60.131
(0.007)		9.336
(0.001)		0.045
(0.001)
Shibishan (SBS)23.6077
120.7045		13471337.50.125
(0.006)		2.517
(0.001)		0.012
(0.001)
Shuangliu (SL)22.2140
120.7961		2551330.00.103
(0.006)		2.489
(0.001)		0.014
(0.001)
Sunmoonlake (SML)23.8519
120.8982		8161333.10.115 (0.006)7.748
(0.001)		0.036
(0.001)
Tahsueshan (THS)24.2326
120.9003		9371433.40.103
(0.006)		4.739
(0.001)		0.027
(0.001)
Taroko (TRK)24.1880
121.6382		9291531.00.102
(0.006)		7.054
(0.001)		0.034
(0.001)
Wulai (WL)24.8663
121.5498		1431046.70.145 (0.007)4.942
(0.001)		0.022
(0.001)
Weiliaoshan (WLS)22.8695
120.6571		6941044.30.144
(0.007)		2.460
(0.001)		0.011
(0.001)
Average 12.537.00.123
(0.006)
N, number of samples used; %P, the percentage of polymorphic loci; uHE, unbiased expected heterozygosity; IA, index of association; rD, modified index of association.
Table
Table 3. FST outliers identified via BAYESCAN and DFDIST and strongly associated with environmental variables. Codes below the environmental columns (aspect, BIO7, BIO12, NDVI, PET, RH, and WSmean) represent strong correlations between FST outliers and environmental variables identified using LFMM (L), Samβada (S), and rstanarm (R).
Table 3. FST outliers identified via BAYESCAN and DFDIST and strongly associated with environmental variables. Codes below the environmental columns (aspect, BIO7, BIO12, NDVI, PET, RH, and WSmean) represent strong correlations between FST outliers and environmental variables identified using LFMM (L), Samβada (S), and rstanarm (R).
LocusDFDIST
FST		BAYESCAN
log10 (PO)		AspectBIO7BIO12NDVIPETRHWSmean
P01_16120.3291000LSRLSS LSRRL
P01_18880.3721000SRSRS RR
P01_22130.3971000SRLSRSRRSR
P03_17600.3401000 SRLLSR R
P03_18900.4061000SRLSR R SR
P03_22000.3451000 LSR RR
P03_34750.2911000 LSRSRL LSR
P05_22910.3461000RLR RR
P08_25660.4931000 LSRS RRSR
P08_29190.4001000SRLSRSRSLRLR
P12_16120.2592.164RR
P12_19560.3231000 LSR R
P12_25910.3441000 RLSRS L
P13_18550.4071000R
P13_21770.4521000 LSR RSR
P13_22340.4341000SRSRLSR R R
P13_29910.3391000 LSRRLS
P19_21110.4871000 RSRSR RSR
P19_26190.3841000LSRLSRSR SR
P19_28120.2392.657S LSRS S
P21_17720.3841000RRRLSR RSR
P21_18650.4131000 RLSRSR LRR
P21_19550.4071000SRLSRSR RR
P21_30130.3661000SRSRSRRLRLR
P35_16350.3611000SLSRSRSRRR
P35_20140.3821000 R RRRR
Aspect (0–360°). BIO7, annual temperature range (°C); BIO12, annual precipitation (mm); NDVI, normalized difference vegetation index (unitless); PET, annual total potential evapotranspiration (kg/m2/year); RH, relative humidity (%); WSmean, mean wind speed (m/s).
Table
Table 4. Relative contribution (adjusted R2) and F test of environmental variables explaining outlier genetic variation in Zingiber kawagoii using a forward selection procedure.
Table 4. Relative contribution (adjusted R2) and F test of environmental variables explaining outlier genetic variation in Zingiber kawagoii using a forward selection procedure.
Environmental VariableAdjusted R2Cumulative Adjusted R2F Value (p)
BIO70.19160.191651.00 (0.001)
BIO120.09840.290030.11 (0.001)
NDVI0.03740.372413.68 (0.001)
RH0.03150.358911.09 (0.001)
WSmean0.02980.388710.16 (0.001)
PET0.02870.417411.63 (0.001)
Aspect0.01180.42925.27 (0.001)
Aspect (0–360°). BIO7, annual temperature range (°C); BIO12, annual precipitation (mm); NDVI, normalized difference vegetation index (unitless); PET, annual total potential evapotranspiration (kg/m2/year); RH, relative humidity (%); WSmean, mean wind speed (m/s).
Table
Table 5. Relative importance and significance of environmental variables explaining variations in the 26 FST outliers based on model averaging using MuMIn. Numbers in parentheses are the Akaike sum of weights (SW) of each environmental variable across all parsimonious predicting models (ΔAICc ≤ 2). In bold, variables receiving strong support (i.e., the 95% confidence interval did not overlap with zero). McFadden’s pseudo R2 was calculated with the variables (predictors) selected as the best model with the lowest AICc used in the generalized linear model. For variables that are part of the best model with the lowest AICc, the sign of regression coefficient is shown: +, positive; −, negative.
Table 5. Relative importance and significance of environmental variables explaining variations in the 26 FST outliers based on model averaging using MuMIn. Numbers in parentheses are the Akaike sum of weights (SW) of each environmental variable across all parsimonious predicting models (ΔAICc ≤ 2). In bold, variables receiving strong support (i.e., the 95% confidence interval did not overlap with zero). McFadden’s pseudo R2 was calculated with the variables (predictors) selected as the best model with the lowest AICc used in the generalized linear model. For variables that are part of the best model with the lowest AICc, the sign of regression coefficient is shown: +, positive; −, negative.
LocusPseudo R2AspectBIO7BIO12NDVIPETRHWSmean
P01_16120.4750.55 (8)0.5 (7)0.93 (14)+0.29 (5)1 (15)+0.32 (5)0.81 (11)+
P01_18880.3871 (4)+0.18 (2)0.21 (2)0.78 (3) −1 (4) −1 (4)+1 (4) −
P01_22130.5601 (3)+1 (3) −0.21 (1) 0.25 (1)1 (3) −1 (3) −
P03_17600.2481 (3) −0.18 (1)1 (3) − 1 (3) −1 (3) −0.41 (1) −
P03_18900.3641 (5)+1 (5)+0.51 (2) −0.86 (4) −0.15 (1)0.13 (1)1 (5) −
P03_22000.1971 (3) − 1 (3) −0.65 (2)1 (3) −0.66 (2)+1 (3) −
P03_34750.7261 (2)+1 (2) − 1 (2) − 0.38 (1)0.62 (1) −
P05_22910.6431 (5)0.11 (1)1 (5)+0.26 (1)0.21 (1)1 (5) −1 (5)+
P08_25660.6460.83 (2)+0.17 (1)0.53 (2)0.47 (1)+1 (3) −0.47 (1)+1 (3) −
P08_29190.749 0.26 (10)1 (3)+0.21 (1)0.74 (2)+ 1 (3)+
P12_16121.0000.2 (1)1 (5)+1 (5) −0.2 (1)0.2 (1) −0.2 (1)0.2 (1)
P12_19560.189 1 (3)+1 (3)+0.66 (2) − 0.18 (1)
P12_25910.3621 (3)+0.24 (1)1 (3)+1 (3) −1 (3) −0.25 (1)1 (3) −
P13_18550.5821 (3) −1 (3)+0.24 (1)1 (3)+1 (3)+0.23 (1)1 (3)+
P13_21770.4941 (2)+1 (2)+1 (2)+1 (2)0.25 (1)1 (2)+1 (2) −
P13_22340.2720.38 (2)0.34 (2)1 (5) −1 (5) −1 (5) −0.12 (1)1 (5) −
P13_29910.614 1 (3)+1 (3) −1 (3) −1 (3) −0.23 (1)0.24 (1)
P19_21110.5870.32 (2)0.47 (2)1 (4) −1 (4)+1 (4) −1 (4)+1 (4) −
P19_26190.8290.08 (1)0.07 (1)0.26 (3)0.24 (3)0.57 (7) −0.92 (10)+0.17 (2)
P19_28120.8460.2 (1)0.2 (1)0.8 (4)+0.4 (2) −0.2 (1)0.2 (1)
P21_17720.2950.34 (2)0.3 (2) 1 (6)+0.11 (1)0.13 (1)1 (6) −
P21_18650.3260.23 (1)1 (3)+1 (3)+1 (3) − 0.21 (1)1 (3) −
P21_19550.6520.63 (5)+0.59 (5)1 (8) −0.6 (5)0.59 (5) −0.49 (4)+0.5 (4) −
P21_30130.3740.78 (2) −1 (3) −1 (3)+1 (3)+ 0.21 (1)1 (3)+
P35_16350.5340.08 (1)0.56 (5)1 (8)+0.18 (2)0.92 (7)+0.4 (3)1 (8)+
P35_20140.2110.85 (4) −0.29 (2)1 (5) −1 (5) −1 (5) −1 (5) −0.4 (2)
Aspect (0–360°). BIO7, annual temperature range (°C); BIO12, annual precipitation (mm); NDVI, normalized difference vegetation index (unitless); PET, annual total potential evapotranspiration (kg/m2/year); RH, relative humidity (%); WSmean, mean wind speed (m/s).
Table
Table 6. The percentage of variation (adjusted R2)-explained outlier variation and variation 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 variations in 26 FST outliers. Fraction [a+b+c] represents total explainable variation.
Table 6. The percentage of variation (adjusted R2)-explained outlier variation and variation 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 variations in 26 FST outliers. Fraction [a+b+c] represents total explainable variation.
Adjusted R2
(Percentage of Total Explainable Variation)		Fp
Environment [a]0.240 (50.6%)14.660.001
Environment + Geography [b]0.189 (39.9%)
Geography [c]0.045 (9.5%)9.760.001
[a+b+c]0.47422.160.001
Residual [d]0.526
Environmental variables used in [a] were aspect; Geographic variable for [c] was calculated using geographical coordinates of sample sites.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, Y.-S.; Liao, P.-C.; Chang, C.-T.; Hwang, S.-Y. Pattern of Adaptive Divergence in Zingiber kawagoii Hayata (Zingiberaceae) along a Narrow Latitudinal Range. Plants 2022, 11, 2490. https://doi.org/10.3390/plants11192490

AMA Style

Li Y-S, Liao P-C, Chang C-T, Hwang S-Y. Pattern of Adaptive Divergence in Zingiber kawagoii Hayata (Zingiberaceae) along a Narrow Latitudinal Range. Plants. 2022; 11(19):2490. https://doi.org/10.3390/plants11192490

Chicago/Turabian Style

Li, Yi-Shao, Pei-Chun Liao, Chung-Te Chang, and Shih-Ying Hwang. 2022. "Pattern of Adaptive Divergence in Zingiber kawagoii Hayata (Zingiberaceae) along a Narrow Latitudinal Range" Plants 11, no. 19: 2490. https://doi.org/10.3390/plants11192490

APA Style

Li, Y. -S., Liao, P. -C., Chang, C. -T., & Hwang, S. -Y. (2022). Pattern of Adaptive Divergence in Zingiber kawagoii Hayata (Zingiberaceae) along a Narrow Latitudinal Range. Plants, 11(19), 2490. https://doi.org/10.3390/plants11192490

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop