Isolation-by-environment as a driver of genetic differentiation among populations of the only broad-leaved evergreen shrub Ammopiptanthus mongolicus in Asian temperate deserts

Whether the effect of migration-selection-drift equilibrium on population structure is governed by spatial or environmental differences is usually elucidated by isolation-by-distance (IBD), isolation-by-environment (IBE), and isolation-by-resistance (IBR) tests. The population structure of Ammopiptanthus mongolicus, a broad-leaved evergreen psammophyte in eastern Central Asia, was previously thought to follow an isolation by distance pattern. However, recent studies have emphasized the effects of environmental factors on its growth and distribution, suggesting an important influence of local adaptation on the genetic structure of the species. Using inter-simple sequence repeat (ISSR) markers, we verified the previously inferred low intra-population variation and high inter-population differentiation. However, in contrast to previous studies, the results of partial Mantel tests and a maximum likelihood population effects mixed model (MLPE) suggested that local climate differences, rather than geographic distances or resistance distances, are the main factor affecting population differentiation. Further analysis with removal of multicollinear climatic variables and univariate MLPE found that summer and winter precipitation were crucial for shaping the current population genetic structure. Since local precipitation is related to the regeneration, colonization, and overwintering survival of A. mongolicus, its influence on demographic change may explain its effect on the population genetic structure. In addition, precipitation is related to terrain despite westward decreases, which explains the independence of genetic difference and geographic distance. The identified role of IBE suggests that collecting germplasm resources from genetically differentiated populations could be a more effective strategy to preserve the overall genetic diversity of the species than the establishment of corridors to enhance gene flow among populations.


Results
Low intra-population genetic variation and high inter-population differentiation. From a total of 200 samples collected from 10 populations ( Fig. 1 and Table 1), 105 sharp and clear bands (loci) of the ISSR marker were recorded, of which 71 loci were polymorphic. The genetic diversity estimated from these 71 among-population polymorphic loci revealed that the percentage of within-population polymorphic loci (%P) ranged from 9.86 (EJNYG) to 29.58% (ALSY and WH). The overall expected heterozygosity (H E ) was below 0.12 and the Shannon index (I) was below 0.16 in all studied populations, with two distant populations (WLTH and ALSY) exhibiting the largest values (Table 2). Although the genetic diversity was slightly higher in the whole species than within populations, it was still low, especially H E (I = 0.452 ± 0.024, H E = 0.296 ± 0.019, total populations). The high genetic diversity of the total population relative to each single population suggests high differentiation among the populations. AMOVA was used to assess the population genetic structure and revealed that 76.58% of the genetic variation was partitioned among populations, while the remaining 23.42% was attributed to differences between individuals within populations ( Table 3). The inference of high genetic differentiation was also confirmed by both the neighbor-joining (NJ) tree (Fig. 2) and discriminant analysis of principal components (DAPC, Fig. 3). Our results supported the existence of ten genetic clusters, indicating that each population had its own genetic signature. In DAPC, seven principal components (PCs) were retained according to the 1000-run K-means algorithm assessed by the Bayesian information criterion (BIC), and the optimal number of clusters was 10, which corresponds to the number of sampled populations. Figure 3 shows both the component and scatter plots; all samples were clearly assigned to their own populations, except one sample in population ALSYSHT with roughly half of its genetic component from HLS and a few samples of HLS with small proportions of genetic admixture with ALSYSHT (the component plot of Fig. 3).
STRUCTURE analysis showed that the optimal grouping number (K) of genetic components was two based on the logarithmic probability change rate of successive K-value data. When K = 2, the populations ALSZNRG, EJNYG, ALYSHT, and HLS clustered together, while the remaining populations formed another group (Fig. 4). (a) The map shows the relative locations of the distribution of A. mongolicus in the inland of temperate Asia; (b) the detailed sampling sites in this study and the topographic variation of the distribution. The current altitude layer is publicly available from WorldClim version 2.0 67 (www.worldclim.org), and the map was generated with the package raster 72 (http://www.rspatial.org/) in R 58 . (2019) 9:12008 | https://doi.org/10.1038/s41598-019-48472-y www.nature.com/scientificreports www.nature.com/scientificreports/ Except for two ALYSHT samples, only weak genetic admixture was detected between these two groups ( Table 3). The second best K was three, with the WH and WLT populations forming the third group (Fig. 8). The grouping pattern of STRUCTURE was consistent with that of the NJ tree (Fig. 2). When K = 10, almost every population had its own unique genetic components, except for a composite component in ALZBURG similar to a part of EJNYG (Fig. 4). Although there were slight differences, the STRUCTURE, DAPC and NJ tree analyses yielded congruent inferences of obvious genetic differentiation.
IBE explains the population structure of A. mongolicus. Models of IBD, IBE, and IBR were tested to explain the genetic differentiation patterns among populations of A. mogolicus. The results of the partial Mantel test suggested that the population genetic structure could be explained by IBE (r = 0.609, p = 0.001) instead of IBD (r = 0.137, p = 0.212), IBR clim (r = −0.031, p = 0.254), or IBR alt (r = 0.179, p = 0.250) ( Table 4 and Fig. 5). However, there was a marginally significant correlation between geographic distance and environmental difference (Mantel test, r = 0.368, p = 0.073), implying that the farther the geographic distance, the greater the environmental difference. Model selections for the maximum likelihood population effect mixed effect (MLPE) revealed that the IBE was the first-ranked model explaining the population genetic structure according to the ranking    Tables 4 and 5). Both IBE and IBR clim attribute population genetic structure to the climatic effect; the former explains the impact of climate differences on the survival and reproduction of colonizers, while the latter emphasizes the facilitation or inhibition of the migration (gene flow) process of organisms by climate differences. However, despite small ΔAIC between IBE and IBR clim , the effect size of the fixed effect (climatic composite resistance distance) is small (fixed estimate = 0.020) in IBR clim , suggesting that the environmental resistance during migration contributes less than the selective pressure after colonization.
The Mantel test has been criticized for high Type 1 error due to multicollinearity 9,28 . Therefore, we removed the bioclimatic factors with multicollinearity and conducted the partial Mantel test again using each retained single factor (Fig. 6) to further explore which bioclimatic variable is the key factor affecting population genetic differentiation. Four bioclimatic factors, bio3, bio4, bio6, and bio18, were retained; only bio18 (precipitation of the warmest quarter) was positively correlated with the genetic distance (r = 0.553, p = 0.007, Fig. 7). Since some bioclimatic factors were removed due to collinearity with bio18, we tested environmental distances based on these individual factors (bio12, bio13, and bio16) for correlations with the genetic distance among populations, which confirmed their positive correlations with genetic distance (bio12: r = 0.563, p = 0.005; bio13: r = 0.554, p = 0.003; bio16: r = 0.567, p = 0.004, by partial Mantel test, conditioning on geographic distance). Bio12 (annual precipitation), bio13 (precipitation of the wettest month), bio16 (precipitation of the wettest quarter), and bio18 (precipitation of the warmest quarter) are all bioclimatic dimensions related to precipitation. According to the monthly precipitation records (Fig. 8a), the annual precipitation mostly accumulates from June to September. The regional precipitation not only decreases westward but is also obviously related to topography (Fig. 8b).
Univariate MLPE regression was also conducted to test the IBE model with each of the 19 bioclimate distances as the fixed variable and the population effect as the random variable. The ranked AIC revealed that both models with bio17 and bio19 as the fixed variable had the smallest AIC values and significantly better fits than the other models (ΔAIC > 5, Table 6). Bio17 and bio19 are the precipitation of the driest and coldest quarters, respectively. In our study area, the coldest and driest seasons are the same, resulting in the same estimates in bio17-and bio19-univariate MLPE. Although the most crucial climatic factor affecting the genetic distance differed between the partial Mantel test (summer precipitation) and MLPE (winter precipitation), both analyses suggest that the regional precipitation difference is the key factor affecting the genetic structure of A. mongolicus.   The test of isolation-by-distance (IBD); (b) the test of isolation-by-environment (IBE); (b) the test of isolation-by-resistance in climate (IBR clim ); (d) the test of isolation-by-resistance in altitude (IBR alt ); (e) the test of correlation between geographic and environmental distance. Among these linear relationships, only the climatic distance was significantly correlated with genetic distance (i.e. IBE), as supported by the Mantel test, partial Mantel test, and MLPE (Table 4). (2019) 9:12008 | https://doi.org/10.1038/s41598-019-48472-y www.nature.com/scientificreports www.nature.com/scientificreports/ Discussion IBE is the best model on population structure. The population genetic structure of A. mongolicus was previously suggested to fit the IBD model 20 , which implies an inverse proportion of effective dispersal to geographical distance [29][30][31] . Over the past decade, accumulating studies have indicated that geographic distance or geographic barriers may not be the only factor affecting gene flow. Environmental differences may be the key factor underlying effective migration 9,10,31 . In this study, we suggest that the adaptability of A. mongolicus to local climate affects its seed germination and colonization. The effect of selection pressure on population differentiation is usually faster than that of drift and could occur at a small geographic scale [32][33][34] . Notably, environmental differences were marginally correlated with geographic distance. We therefore suggest that the previous inference of IBD 20 could be due to the intercorrelation between geographic and environmental differences. The increasing number of open databases is now helping to clarify ecological and evolutionary phenomena. A meta-analysis showed that 74.3% of phylogeographic studies (52 of 70 studies) revealed significant IBE patterns, including 37.1% (27 studies) revealing spatial autocorrelation (i.e. covariates with IBD) 9 . Similarly, from 106 IBE studies, Shafer and Wolf 10 reported effect sizes of 0.34 (95% CI 0.24-0.42) and 0.26 (95% CI 0.13-0.37) for a mixed-effect model with and without controlling spatial autocorrelation, respectively, suggesting that spatial autocorrelation reduces IBE correlations for environmental variables. These studies indicated the relevance of environmental autocorrelation for the spatial effect (i.e. IBD). That is, the previous inference that the population differentiation of A. mongolicus aligns with geographic distance 20 probably reflects differential adaptation to the local climate. Differential adaptability to heterogeneous environments provides a better explanation than IBD in A. mongolicus, i.e. divergent selection is more important than neutral processes. Genetic draft explains low genetic diversity. The low estimates of genetic diversity are consistent with the previous estimation by Ge et al. 20 , which included populations located farther south but no populations in Alashan (ALSY, ALSZNRG, ALSYSHT, and ALSZCHE). The genetic diversity of A. mongolicus was also lower than that of other desert species estimated by ISSR, e.g. Achillea fragrantissima in Egypt 35 Table 5. The order of the models corresponds to their ranking from best (smallest AIC and BIC) to worst. The first-ranked model is marked in bold.  Table 5. Model selection for 15 MLPE models. The order of the models corresponds to their ranking from best (smallest AIC and BIC) to worst. The first-ranked model is marked in bold. D gen , D geo , D env , R clim , and R alt denote the genetic distance (i.e. F ST /(1 − F ST )), geographic distance, environmental difference, and resistance distances estimated from climatic composite resistance surface and from altitudinal resistance surface among populations, respectively.  Environmental heterogeneity would reduce the chance of dispersal 38,39 , and the constraint of the range of distribution may lead to deleterious erosion of genetic diversity due to increased inbreeding and genetic draft 39 . Precipitation is an important limiting factor for the reproductive success of desert plants. The growth pattern of A. mongolicus is similar to that of desert deciduous plants or summer annuals, with blossoming and germination during the high-rainfall season 40 . Rapid blooming is advantageous for plant reproductive success in the desert 40 . However, pollinators tend to visit flowers of the same or adjacent plants instead of distant flowers in the short blooming season, which may reduce the outcrossing rate of A. mongolicus (inbreeding coefficient F IS > 0 in all loci 21 ).

Models
In addition, rainfall restrictions in deserts may also result in strong selection on A. mongolicus. With a selective sweep, genetic variation of adjacent genes decreases along with adaptive loci, which will even expands to most genome regions. Compared with other plants that may also be affected by genetic draft, such as Dactylis glomerata L. in the plateau of Central Asia and Western China 41 , A. mongolicus exhibits extremely low intra-population genetic variation, suggesting that regional environmental pressures (especially precipitation) in the desert have a more severe impact on this broad-leaved green plant. Rainfall-induced declines in outcrossing opportunities and strong selective sweeps could explain the low genetic variation of A. mongolicus, which may also be resistant to the rescue effect of gene flow among populations. Differential local precipitation is the key to population differentiation. Summer rainfall almost completely determines the annual precipitation in the distribution of A. mongolicus. In general, the annual precipitation tends to decline in a southeast-to-northwest direction across the Asian continent 42 , but fluctuations in terrain (e.g. the Hetao Plain, Helan Mountains, and Mongolian Highlands) make the local climate more complicated. Such local differential precipitation may have long been the selective pressure not only for the breeding and dispersion of A. mongolicus but also for the water supply in the dry season.
Several studies have indicated that water is the key factor affecting the seed germination 43 and seedling growth 44 of A. mongolicus. In summer (July and August), the legume of A. mongolicus is ripe and dehiscent, and seeds fall off, quickly absorb water and germinate 43 . In a manipulation experiment, 85% of seedlings wilted in a 5-day drought treatment 44 , indicating that the demand for water is a limiting factor for the regeneration of A. mongolicus. Due to the lack of defoliation in winter, supplementing evapotranspiration with some precipitation may also affect A. mongolicus survival in winter. Although the local precipitation is small and varies little in winter (the cumulative precipitation ranges from 0 to 6.06 mm in Dec~Feb), such differences may cause local adaptation. www.nature.com/scientificreports www.nature.com/scientificreports/ Differential adaptability to precipitation among populations might accelerate the process of population differentiation by stalling maladaptive immigration.
As described above, differential local precipitation might also affect pollinators' species and visiting frequency 45 . Differences in precipitation could vary the ratio of bee pollinators to fly pollinators; the former require dry soil for nesting, whereas the latter require moist environments for larval growth and metamorphosis 46 . Changes in the abundance of pollinators could affect the success of pollination and seed yield, even though the connectivity of the plant and pollinator relationship may not be disturbed by precipitation 47 . In addition, precipitation can also affect rhizobium symbionts 48 and pathogen infectivity 49 , thereby affecting plant health and population regeneration. The presence of some endosymbiotic fungi (dark septate endophytes) can facilitate the growth of A. mongolicus under drought conditions 50 . We have not explored differences in soil microorganisms and endophytes among different populations, but distributional differences of these endophytes may also cause differential adaptability to precipitation among populations. The local adaptation caused by differential precipitation may lead to divergent directions of genetic draft, which may explain apparent population differentiation in A. mongolicus.

Concluding remarks.
In conclusion, the selective pressure of the environmental gradient (differences in precipitation) is strong for A. mongolicus and likely explains the low genetic variation within populations and high population differentiation. Most individuals carry not only locally adapted genes but also homogenized genomic variation, which decreases successful emigration to populations with different environments, i.e. selection against maladapted dispersers 12,31 . A. mongolicus is the only evergreen shrub in the desert of Northwest China and is an important wintering place for several small animals, i.e. an umbrella species. Given the low genetic variation within populations and maladapted gene flow among populations, every population is a unique evolutionarily significant unit and should be considered as a unique management unit for conservation. The high dependence of adaptability on precipitation is not propitious for effective gene flow among populations. Therefore, the establishment of ecological corridors 51-53 may not be an appropriate strategy for conservation. Germplasms from different populations should be actively preserved to maintain the complete gene pool and increase the evolutionary resilience 9 of A. mongolicus in the face of increasingly severe climate change.

Materials and Methods
Species studied and sampling. The genus Ammopiptanthus was suggested to have originated from the broad-leaved evergreen Tethyan flora 54 , as supported by molecular dating indicating that the genus Ammopiptanthus split from its sister taxa in the early Miocene (chloroplast DNA matK sequences:19.6 Mya; nuclear ITS sequences: 21.8 Mya) 55 . Ammopiptanthus mongolicus is discontinuous distributed in western Inner Mongolia, northern Ningxia and Northern Gansu in China, ranging from 36°27′N-42°01′N, 102°36′E-108° 49′E 26 . The sampling area of this study was the core distribution of A. mongolicus in Inner Mongolia, China. We chose 10 populations covering the main distribution range of A. mongolicus (Table 1). Fresh leaves were sampled from 20 individuals per population, and each sampled plant was distant from other plants by at least 20 meters. A total of 200 individuals were sampled. The sampled leaves were placed immediately in a liquid nitrogen tank and stored in a −20 °C refrigerator after carrying to the laboratory.  Table 6. Summary results of the MLPE and the model selection. The order of the models corresponds to their ranking from best (smallest AIC and BIC) to worst. The first-ranked models are marked in bold.
(TIANGEN Biotech Co., Ltd., Beijing, China). DNA quality was checked by agarose gel electrophoresis and by the DNA absorbance ratio (OD 260 /OD 280 : 1.7~1.9) in a WD-9403C UV Viewing Cabinet (BEIJING LIUYI Biotechnology Co., Ltd., Beijing, China). ISSR amplification was performed using fifteen primers (Table 7) with the following PCR procedure: denaturation at 95 °C for 5 min, followed by 35 cycles of denaturation at 95 °C for 1 min, annealing at the proper temperature for 1 min (Table 7), and extension at 72 °C for 1 min, with a final 10-min extension at 72 °C. PCR was conducted in an MJMini personal thermal cycler (Bio-Rad, Hercules, USA) and T100 ™ thermal cycler (Bio-Rad, Hercules, USA). All PCR products were checked by 1.5% agarose gel electrophoresis, and the appearance of bands was read. Ghost bands were excluded by comparison with a negative control in which water was used as the template with the same ISSR protocol. The ISSR experiments were repeated twice to ensure that the peak signals affirming the bands (loci) were not PCR errors. Only loci that were consistently present or absent in all preliminary tests were read in the formal experiment.
Genetic diversity and population genetic structure. The genetic diversity was estimated by the indices of percentage of polymorphic loci (%P), average number of different alleles per locus (N A ), effective number of alleles per locus (N E ), Shannon's information index (I), expected heterozygosity (H E ), and unbiased heterozygosity (UH E ) using GenAlEx v. 6.5 56 . The contributions of genetic variation between and within populations were assessed by analysis of molecular variance (AMOVA). The significance of genetic differentiation between populations was estimated by Φ ST under 999 permutations. We also conducted DAPC to determine if the spatially structured population was also genetically structured using the package adegenet 57 in R 58 . The best clustering number (k) was inferred by the k-means algorithm with 10 6 simulations evaluated by BIC (the elbow in the BIC curve and the smallest BIC). The optimal number of PCs retained for DAPC was evaluated by a-score optimization. A component plot and scatter plot were drawn to illustrate the population clustering pattern of A. mongolicus. Patterns of genetic admixture were assessed by Bayesian clustering analysis, a population model-based approach based on Hardy-Weinberg and linkage equilibria 59 , with the assistance of STRUCTURE 2.3.4 60 . We estimated the posterior probability of the grouping number (K = 1-20) by 10 independent runs using 10 6 steps of Markov chain Monte Carlo (MCMC) replicates after 10% burn-in for each run to evaluate consistency. The best grouping number was evaluated by ΔK 61 in STRUCTURE HARVESTER ver. 0.6.94 62 . In addition, to understand the relationships between each of the lineages, we transformed the number of differences in ISSR loci between individuals into a triangular matrix and then constructed an NJ tree using MEGA6 63 .
Testing IBD, IBE, and IBR. To test the effects of geographic distance and environmental differences on genetic structure, the partial Mantel test was conducted using the R package vegan 64 . We calculated the pairwise genetic distances among populations using F ST /(1 − F ST ) 65 . Euclidean distances of geographic distance were calculated using the R package fossil 66 . We also collected 19 standard bioclimatic variables of 10 sampling sites as environmental data from WorldClim version 2.0 67 . We considered the 19 bioclimatic variables as different environmental space vectors and used the Canberra distance to calculate the distance between populations in this vector space. To test IBD, the genetic distance was used as the response, the geographic distance as the predictor, and the environmental distance as the condition factor. To test IBE, the roles of environmental distance and geographic distance were interchanged. In addition, to test whether these environmental factors impeded gene flow, the climatic composite resistance surface was transformed from raster layers of the bioclimatic variables using Circuitscape 4.0 68 . We also transformed the altitudinal layer into an altitudinal resistance surface. These two  www.nature.com/scientificreports www.nature.com/scientificreports/ resistance surfaces were used to test IBR, namely IBR clim and IBR alt , respectively. The Mantel statistic was based on Spearman's rank correlation with 9999 permutations.
Mantel and partial Mantel tests have been strongly criticized for inflated type-I error, potential collinearity between environmental variables when building an environmental dissimilarity matrix, low power, etc. 9,28,69 . Therefore, we fit linear mixed-effects models using the MLPE parameterization, which has been found to perform better than other regression-based statistical approaches 70 , to account for the non-independence of values within pairwise distance matrices and to distinguish the effects of multiple independent variables. Mixed-effects models were fit by maximum likelihood to test the effects of fixed factors (geographic, bioclimatic, and two resistance distances) with the random effect of populations. AIC and BIC were used as the objective criteria to evaluate model fit from four models of single fixed factor, six combinations of double fixed factors, four combinations of triple fixed factors, and the full model (the combinations of all fixed factors, Table 5).
Since IBE was suggested as the first-ranked model by both the partial Mantel test and model selection for MLPE (see Results), we further identified the most crucial environmental factors affecting genetic distance using two strategies. First, we re-executed the partial Mantel test by calculating the distance of each environmental factor between populations. To avoid unnecessary weighting due to intercorrelations among bioclimatic variables, we used variance inflation factor (VIF) analysis to reduce multicollinearity 71 . We discarded variables with high VIF values (VIF > 10) and then calculated the distances of the retained bioclimatic variables among populations to test the IBE hypothesis (Fig. 6). This remaining factor is the most likely environmental factor affecting the population genetic structure of A. mongolicus. Testing one variable by a Mantel (or partial Mantel) test has been suggested to be more credible than testing multiple variables 28 . Second, we performed model selection to evaluate 19 models with every single bioclimatic distance as the fixed factor in MLPE. These single-bioclimate-distance IBE models were ranked by AIC, and the model with the lowest AIC was suggested as the best one for the prediction of population genetic structure.

Data Aavailability
All genetic and environmental data used in this study are available in the Supplementary Data.