Optimizing conservation of Ulmus wallichiana, threatened Himalayan species by combining population genetics and species distribution modeling

doi:10.21203/rs.3.rs-5003442/v1

Download PDF

Research Article

Optimizing conservation of Ulmus wallichiana, threatened Himalayan species by combining population genetics and species distribution modeling

https://doi.org/10.21203/rs.3.rs-5003442/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Global climate change and human action are dismantling many ecosystems at an alarming rate, causing an unprecedented erosion of biodiversity. In this context, investigation of the threats and evaluation of the remedies to protect biological resources are necessary. This study aims to define the appropriate conservation strategy for Ulmus wallichiana, a vulnerable Himalayan Elm species. The structure of genetic diversity was investigated within and among its populations using nuclear microsatellites (SSR). Further, ecological niche modelling was carried out to ascertain the climatic suitability in the Indian Western Himalayas. Our integrative approach identified genetic diversity hotspots and long-term suitable areas across the Western Himalayas for U. wallichiana, two key aspects required for its conservation in the future. Moreover, a negligible level of inbreeding and signs of a genetic bottleneck in several populations were observed. Germplasms from the genetic diversity hotspots for propagation in the climatic hotspots may ensure better genetic diversity under a safer climate in the future. We believe that this integrative approach will guide the formulation of cost, time, and labor-intensive effective strategies in restoring this Himalayan vulnerable species.

Population structure

Bioclimatic variables

Ecological niche modelling

Restoration

SSR markers

The Himalayas represent an important global biodiversity hotspot and an ideal location for studying the impact of climate change. The region has experienced warmer winters in recent years, with an increase in the mean temperature of the warmest month by 0.76 ± 0.25°C and an increase in the mean temperature of the coldest month by 3.65 ± 2°C. These temperature changes have led to geographical range shifts in 87% of the 124 endemic plant species from the Lachen and Lhonakh valleys of Sikkim (Telwala et al., 2013), but there is still a lack of spatio-temporal research on risk perception and management of climate change in the Himalayas (Schneiderbauer et al., 2021). In addition, human activities also impact the range shifts of species (Chen et al., 2011; Peterson et al., 2002; Zimmerer et al., 2019). While the current plant community structure has changed significantly from the past, it may potentially lead to the extinction of endemic species (Braje & Erlandson, 2013; Manne & Pimm, 2001). The changes in the current plant community structure in the Himalayas are likely influenced by several factors, including climate change, human activities, and habitat fragmentation. Climate change, as mentioned earlier, has led to warmer winters and shifts in temperature, which can affect the distribution and behaviour of plant species. Human activities, such as deforestation, overuse of plant resources (including medicinal purposes, timber, firewood, and fodder), and habitat destruction, can, directly and indirectly, affect the plant community structure. Habitat decline and loss of connectivity between populations may indeed lead to genetic erosion possibly affecting the long-term persistence of species, especially due to inbreeding depression (Young AG, 1996). For the scientific community, practitioners, and policymakers, a better understanding of processes at play is thus necessary to optimize conservation actions (Howden et al., 2007).

The Himalayan Elm (Ulmus wallichiana Planch) is a temperate deciduous tree that is classified as vulnerable in the IUCN Red Data List (IUCN, 1998). It grows in Afghanistan, Pakistan, India, and Nepal at an elevation from 800 to 3000 m (Batool et al., 2021). Overuse of U. wallichiana for medicinal purposes, timber, firewood, fodder, and thatching has caused a significant decline in its populations (Batool et al., 2021), especially in the context of habitat fragmentation in the Himalayas. Immediate attention is required to rescue this species from extinction in the wild. Population genetics of U. wallichiana may help to identify conservation units (i.e. differentiated populations, hotspots of diversity) and assess the impact of habitat loss and connectivity on their diversity, as recently shown in different European Ulmus species (Buiteveld et al., 2016; Goodall-Copestake et al., 2005; Nielsen & Kjær, 2010; Vakkari et al., 2009; Venturas et al., 2013). Such studies have revealed contrasting patterns of genetic diversity among taxa, which was high in Ulmus lamelosa C. Wang & S.L. Chang (Liu et al., 2016), and moderate in Ulmus crassifolia Nutt. and Ulmus pumila L. (Zalapa et al., 2008), and low in Ulmus glabra Huds. (Vakkari et al., 2009).

Integrating molecular data to evaluate population genetics with ecological niche modelling to ascertain future climatic suitability habitat would also provide better insights into developing a successful conservation program (Majeed et al., 2021). Molecular tools, such as simple sequence repeat (SSR, or microsatellite) markers can be used to study the effect of landscape structure on the geographic pattern of genetic variation (Litt & Luty, 1989; Storfer et al., 2007), and to identify the zones harbouring greater genetic diversity (genetic diversity hotspots) for U. wallichiana. Niche modelling can reveal the most suitable climatic zones (hereafter referred to as climatic hotspots) for U. wallichiana in the future. Additionally, GIS and remote sensing technology can also be used to analyze molecular data and understand the influence of the environment on genetic variation (Holderegger & Wagner, 2008). This study thus aims to identify and prioritize genetic diversity hot spots and climatic hotspots for effective and successful conservation of U. wallichiana.

2.1. Plant Material

Leaf samples of U. wallichiana were collected during the period 2016–2018 from the Indian Western Himalayas, covering the union territory of Jammu and Kashmir (J&K) and the state of Himachal Pradesh (HP; Fig. 1). The details of the sampled populations are presented in a supplementary file (S1). DNA of 179 samples from 14 populations (Supplementary file S1) was isolated with a CTAB method (Doyle, 1990) and stored at -20˚C for subsequent use.

2.2. Genetic diversity

Previously developed in-house 28 SSR loci (Singh et al., 2021) were utilized for genotyping of U. wallichiana populations using polyacrylamide gel electrophoresis followed by silver staining (Huang et al., 2018). The frequency-based genetic parameters like the number of alleles (Na), the effective number of alleles (Ne), observed heterozygosity (H_O), expected heterozygosity (H_e) and Fixation index (F_ST) were determined population-wise by GenAlEx v6.5 (Peakall & Smouse, 2006). Genetic diversity hotspots were identified through interpolation of He using the Inverse Distance Weighted (IDW) tool implemented in ArcGIS v10.5.

2.3. Population connectivity

The complex topography of the Himalayas may influence the gene flow, thus, BARRIER v2.2 (Manni et al., 2004) was used to assess the genetic barriers across the landscape. The tool uses genetic and geographical distance matrix as input, and Monmonier's maximum difference algorithm to identify the genetic barriers. The output from BARRIER was georeferenced and digitized using ArcGIS for visualization of the genetic barriers. This approach provides the directions and locations of barriers to show where the geographical patterns of the variables differ.

2.4. Population differentiation and structure

The hierarchical partitioning of genetic variance within and among the populations was examined using an AMOVA test implemented in GenAlEx (Peakall & Smouse, 2006). In this analysis, the Euclidean genetic distance matrix was used as input, Codom–Allelic mode as a method of calculation, and the number of permutations used was 999. We also used STRUCTURE v2.3.4 (Pritchard et al., 2000) to reveal the actual number of genetic clusters in our sampling with an admixture model and allele frequencies correlated. We used a Burnin period of 100,000 and the number of MCMC repeats after Burnin of 100,000 for K = 1 to 14, with 10 replicates each. The output from the STRUCTURE was evaluated through STRUCTURE HARVESTER (Earl & VonHoldt, 2012), which determines the optimum number of genetic stocks (separate populations), based on the highest value of ΔK (Evanno et al., 2005). We also used CLUMPAK (Kopelman et al., 2015) to visualize the results of STRUCTURE and check the consistency of the outcome of each replicate for K = 1 to 14. Moreover, we also used a Principal Coordinate Analysis (PCoA) for population clustering. The Euclidean genetic distance matrix was used as an input for clustering.

2.5. Assessment of isolation by distance (IBD) and isolation by environment (IBE)

To evaluate the relationship between geographical and genetic distance, we used Mantel's test (Mantel, 1967) implemented in GenAlEx with the expectation that populations are expected to become more distant with increasing geographical distance (IBD). To evaluate the relationship between environmental variables and population differentiation, we used a correlation analysis between F_ST and Bioclimatic variables, which may assess the IBE effect. The bioclimatic variables were retrieved from the WorldClim v2 at a spatial resolution of 30 seconds (Fick & Hijmans, 2017). To reduce collinearity among the variables and improve the statistical power of the analysis, we have performed a correlation analysis on the bioclimatic variables. Collinearity occurs when two or more variables are highly correlated, making it challenging to distinguish their individual effects on the response variable, we have set a threshold for collinearity (e.g., -0.7 < r < 0.7) and retained only those variables that were less correlated with each other. By removing highly correlated variables, they could focus on a subset of environmental factors that provide unique information for the analysis. Distance matrices of bioclimatic variables and F_ST matrix were used as input for regression. The distance matrix of bioclimatic variables was created using the dist function in R, whereas the distance matrix of F_ST was created using GenAlEx. Moreover, we used a regression analysis to evaluate whether population differentiation in U. wallichiana is also governed by elevation or not. The altitude data was recorded during field visits and Digital elevation model (DEM) data were downloaded from the digital elevation data set of the United State Geological Survey (USGS). A correlation between altitude and genetic diversity of U. wallichiana populations was carried out to assess the altitudinal influence on the diversity of populations.

2.6. Bottleneck analysis

Genetic bottleneck reflects a strong reduction of genetic diversity in a transitory small population (Chakraborty & Nei, 1977; Mayr, 1963), which could have serious consequences on its survival and evolution. So, the quantification of this diversity loss is essential prior to planning any conservation measures. We used BOTTLENECK v1.2.02 (Piry et al., 1999) to test for any recent reduction in genetic diversity for each population. BOTTLENECK takes heterozygosity excess into account. This heterozygosity excess should not be confused with the excess of heterozygotes resulting from deviations from Hardy–Weinberg equilibrium expectations (Sastre et al., 2011). This heterozygosity excess test asserts that, in a contracted population, allelic diversity reduces faster than observed heterozygosity; the latter is thus expected to be transitory larger than expected at mutation-drift equilibrium (Cornuet & Luikart, 1996). The presence of a recent bottleneck can be tested with three methods: the Two-Phase Model (TPM), which allows multiple-step mutation (recommended for SSR markers), the Infinite Allele Model (IAM), and the Stepwise Mutation Model (Di Rienzo et al., 1994; Haanes et al., 2011; Kimura & Ohta, 1978; Luikart et al., 1998). For heterozygosity tests, the TPM model was run for each population with the proportion of SMM in TPM = 0.003 and geometric distribution variance for TPM = 0.36 for 1000 iterations. Under the Two-Phase Model (TPM), sign tests, standardized differences tests, and Wilcoxon sign rank tests were conducted. Visual representation of the bottleneck was performed through a mode-shift indicator test using BOTTLENECK.

2.7. Ecological niche modelling

A total of 179 GPS coordinates of U. wallichiana were recorded through field visits from Indian West Himalayan from 2016 to 2018 (Supplementary file S1). Additionally, three occurrence data points for U. wallichiana were obtained from the Global Biodiversity Information Facility (GBIF) database (Global Mountain Biodiversity Assessment (GMBA). The study area for the niche modelling consisted of the Indian Western Himalayas (J&K and HP). The environmental variables were retrieved from the WorldClim database for the current and future periods (2070) at a spatial resolution of 30 seconds. The accuracy and certainty of the spatial distribution modelling are compromised by the collinearity among the variables. Generally, collinearity increases uncertainty and decreases statistical power (De Marco & Nóbrega, 2018). We used a correlation test on the 19 bioclimatic variables, followed by the selection of only less collinear variables (-0.7 < r < 0.7) to reduce the uncertainty and to ensure good statistical power (Mateo et al., 2013; Meier et al., 2010; Syfert et al., 2013; Zhu & Peterson, 2017). Then, the Maxent model was used for predicting the suitability of the environment for U. wallichiana (Phillips et al., 2006). Maxent quantifies statistical relationships between predictor variables at the observed and background locations of the area of interest. However, recent studies revealed that the default Maxent settings are not always appropriate and sometimes can result in non-optimal models, especially under the scenario of a small number of occurrence points. Evaluation of the best potential combination of parameters in terms of feature classes (FC) and regularization multiplier (RM) is, therefore, recommended to select the most appropriate model (Morales et al., 2017). We used ENMeval (Muscarella et al., 2014) to determine the best-fitting FC and RM. Further, processing of the environmental predictors was performed in R. Twenty percent of the occurrence points were withheld to validate the accuracy of prediction and only 80% were used as the training data to run the model. A receiver operator curve (AUC) was generated, whose value ranged from 0 (no reliability) to 1 (perfectly reliable), to evaluate the prediction accuracy by assessing true positives as a function of false positives. The current and future predicted environment suitability file was georeferenced on the boundary map of J&K and HP using ArcGIS v10.5. The shapefile was created and edited with an editor toolbar. The shapefile of genetic hotspots was used to measure the hotspot area. Current and future predicted suitability shapefiles were imported into Earth Explorer to identify the exact locations of suitable environments in current and future climate conditions.

2.8. Modelling the temperature determinants of genetic diversity in the study area

To predict the spatial relationship between the land surface temperature (LST) and genetic diversity, geographically weighted regression (GWR) and Ordinary Least Squares (OLS) methods were used. These methods assume a constant relationship between the variables. We used GWR and OLS implemented in ArcGIS v5.0 for performing spatial modelling to ascertain the relationship of LST and emissivity with genetic diversity across the Indian Western Himalayas. The MODIS (Moderate Resolution Imaging Spectroradiometer) LST and emissivity data were retrieved from NASA LPDAAC (Land Process Distributed Active Archive Centre) collections (Centre) from January to June 2008 to 2019. The genetic diversity (H_O) represented a dependent variable, whereas the LST and emissivity of January to June represented the six independent variables. The Akaike Information Criterion (AIC) was used to choose the best model for our data. If the difference between the AIC values were > 4, the model with the lowest AIC was chosen as the best, whereas if the difference between AIC values was < 4, both have equal potential to explain the results (Charlton et al., 2009).

3.1. Genetic diversity, barriers, and recent bottleneck

The details of the population-wise primary genetic parameters are presented in Table 1. The number of alleles (Na) per locus ranged from 2 to 3.27, with an average value of 2.75, whereas the effective number of alleles (Ne) ranged from 1.76 to 2.63, with an average of 2.15. The observed heterozygosity per population (H_O) ranged from 0.43 (Pop1) to 0.83 (Pop14), with an average value of 0.63, and expected heterozygosity (H_e) ranged from 0.313 to 0.581 with an average of 0.470, indicating substantial genetic diversity in U. wallichiana. A total of 78 alleles were amplified on 18 SSR loci. The allelic pattern across the populations showed a total of 12 private alleles (7 in pop3, 1 in pop7, and 4 in pop9). The Fixation index (F_ST) ranged from − 0.47 (pop14) to -0.16 (pop3), with an average value of -0.33 per population.

Table 1

Population wise genetic indices and Population bottleneck signs
Pop Id	Na	Ne	Ho	He	Fis	Sign test	SD test	Wilcoxon Test	Mode shift curve
Pop1*	2	1.765	0.43	0.313	-0.375	0.01	0.001	0.00073	Shifted mode
Pop2*	2.556	2.098	0.699	0.490	-0.406	0.007	0.00083	0.00168	Shifted mode
Pop3	2.667	1.938	0.496	0.426	-0.162	0.401	0.07	0.2216	normal L-shaped
Pop4*	2.833	2.268	0.605	0.500	-0.226	0.00426	0.0009	0.00001	Shifted mode
Pop5*	3	2.63	0.772	0.581	-0.325	0.00002	0	0	Shifted mode
Pop6*	2.833	2.219	0.708	0.492	-0.437	0.0076	0.00384	0.00067	Shifted mode
Pop7	2.889	2.217	0.725	0.493	-0.428	0.0802	0.06178	0.05419	normal L-shaped
Pop8*	2.333	1.919	0.551	0.394	-0.408	0.03653	0.01152	0.00116	Shifted mode
Pop9	2.667	1.953	0.502	0.415	-0.204	0.43905	0.19271	0.16447	normal L-shaped
Pop10	3	2.245	0.657	0.498	-0.289	0.13301	0.07524	0.06071	normal L-shaped
Pop11*	2.778	2.143	0.689	0.473	-0.448	0.41552	0.01945	0.0145	Shifted mode
Pop12*	2.611	2.125	0.529	0.435	-0.226	0.02214	0.04148	0.01767	Shifted mode
Pop13*	3.111	2.24	0.68	0.509	-0.334	0.01815	0.00513	0.02527	Shifted mode
Pop14*	3.278	2.413	0.83	0.565	-0.454	0.00529	0.01046	0.00032	Shifted mode
Mean	2.754	2.155	0.634	0.470	-0.337
Na = number of alleles, Ne = effective number of alleles, Ho = observed heterozygosity, He = Expected Heterozygosity, Fis = inbreeding coefficient, SD = Standardized Difference test, * = the populations undergone genetic bottlenecks, the values under Sign test, SD test and Wilcoxon test represent the probability values for heterozygosity excess.

The areas harboring greater genetic diversity (genetic diversity hotspots) are represented in Fig. 2. Among populations, the genetic diversity of pop2, pop4, pop5, pop10, pop11, pop13, and pop14 is higher compared to pop1, pop3, pop8, pop9 and pop12. We identified one genetic barrier hindering the gene flow across the populations on its opposite sides (Fig. 2), which lies on the ranges of Ghariankalan, Marothi, Licha Top, Banjala, Tipri, Udrana, Jalga, Durga Nagar, Amritgarh, Lass, Pora, Keeru, and Brammah. The observed barrier between populations using Toposheets and Earth Explorers are represented in supplementary files S2 and S3, respectively. Bottleneck analysis showed evidence of a recent genetic bottleneck as the function of heterozygosity excess in ten populations of U. wallichiana. With p-values < 0.05 in two of the sign, standardized difference, and Wilcoxon tests, we found a significant excess of heterozygosity in these populations (Table 1). Furthermore, the distribution of alleles in these populations showed a shifted mode from the normal L-shaped curve, further supporting the existence of a recent bottleneck.

3.2. Population differentiation and structure

According to the AMOVA test, 19% of the total genetic variation in the populations of U. wallichiana was found among populations, while 81% was found within populations. The mean population differentiation showed moderate genetic differentiation (F_ST = 0.18). The pairwise population genetic differentiation showed that pop4 and pop6 were the least differentiated (F_ST = 0.068), while pop1 and pop8 were the most differentiated (F_ST = 0.321; Supplementary file S4). The model-based clustering approach through STRUCTURE, revealed two optimal genetic clusters (Fig. 3), which was also supported by the PCoA analysis. The PCoA, which was based on genetic distance, grouped the populations into two clusters and explained 24.05% and 18.41% of the variation in the first two axes, respectively (Fig. 4A). The individual replicates were also analyzed through CLUMPAK, which showed consistency across individual runs of each replicate. Consistency across individual runs of each replicate through CLUMPAK.

3.3. Isolation by distance (IBD), isolation by environment (IBE), and altitudinal impact on population differentiation

The geographical distance between populations was found to be correlated with the genetic distance, as expected by the IBD model. Both the PCoA and STRUCTURE analyses showed that the clustering of the populations corresponded with their geographical locations, and this was statistically confirmed by Mantel's test, which showed a significant positive correlation between the genetic distance and the geographic distance (r = 0.126, p = 0.04). The bioclimatic variables bio1 (annual mean temperature, r = − 0.21, p = 0.04), bio7 (temperature annual range, r = 0.27, p = 0.009), bio17 (precipitation of driest quarter, r = -0.29, p = 0.004), bio18 (precipitation of warmest quarter, r = 0.33, p = 0.001), and bio19 (precipitation of coldest quarter, r = -0.26, p = 0.01) were found to have a significant correlation with population differentiation (Fig. 4B to 4F). Among these variables, only bio7 and bio18 showed population clustering similar to the genetic distance according to the PCoA (Fig. 4C, 4E). The PCoA analysis using the distance matrices of these bioclimatic variables showed that they explained a significant amount of the variation in the populations (37.27% and 27.25% for bio1, 55.91% and 19.11% for bio7, 45.13% and 26.25% for bio17, 73.30% and 14.73% for bio18, and 73.30% and 14.73% for bio19 in the first two axes). These variables may possibly be associated with isolation by environment (IBE). The analysis also showed a weak but significant relationship between H_o and altitude (R² = 0.26, p-value = 0.03) indicating that altitude explains approximately 26% of the variation in observed heterozygosity. This relationship is still relatively weak, as indicated by the low R-squared value. This observed relationship could be due to factors, such as genetic drift, natural selection, or gene flow.

3.4. Ecological niche modelling

In this study, bioclimatic variables including bio1, bio2, bio3, bio14, bio18, and bio19 were used to predict the current and future suitability of the environment for U. wallichiana. After using ENMeavel to determine the best model parameters, it was found that the Linear, Quadratic, Hinge, Product (LQHP) was the most suitable FC and that an RM value of 1 was the most suitable for our dataset, as measured by the lowest Akaike Information Criterion (AIC). The relative contribution of the variables in predicting the suitability of the environment is presented in Table 2. After conducting a collinearity test and filtering highly correlated variables (Supplementary file S5), a common model was created and found to be reliable, with an area under the curve (AUC) close to 1 for both the test and training datasets (Supplementary file S6). The predicted suitability of the environment for U. wallichiana in the Western Himalayas is shown in Figs. 5A and 5B for the current and future datasets, respectively. These figures reveal that pop1, pop4, and pop5 are expected to receive a more favorable environment, whereas pop3, pop13, and pop14 are expected to receive a less favorable environment in the future. The comparison of the predicted suitability of the current and future models (Supplementary file S7) also showed a 5323 km² current environmentally suitable area for U. wallichiana and a 5223 km² suitable area in the future, with an estimated overall reduction of 100 km² in environmental suitability. Details of the analysis can be found in supplementary files S8 and S9.

Table 2

Relative contribution of bioclimatic variables
Sr.No.	Bioclimatic Variables	Contribution in current	Contribution in future
Sr.No.	Bioclimatic Variables	prediction (%)	prediction (%)
1	Bio1(Annual Mean Temperature)	0.04	10.94
2	Bio2(Mean Diurnal Range (Mean of monthly (max temp- min temp))	0.05	16.95
3	Bio3(Isothermality)(Bio2/Bio7) (x100)	11.66	0.45
4	Bio14(Precipitation of Driest Month)	79.69	43.21
5	Bio18(Precipitation of Warmest Quarter)	0.07	0.03
6	Bio19(Precipitation of Coldest Quarter)	8.64	28.39

3.5. Modelling the temperature determinants of genetic diversity in the study area

The average land surface temperature (LST) and visualization of LST on the landscape of the 14 sampled populations were analyzed over 12 years from 2008 to 2019 (Supplementary files 11, 12 respectively), and the results are depicted in Figs. 6A and 6B. The highest average LST of 34.9°C and the lowest average LST of -0.88°C were observed in May and January, respectively. Both the geographically weighted regression (GWR) model and the ordinary least squares (OLS) model were used to examine the relationship between LST and genetic diversity across the study area. It was found that both models had similar explanatory power, with the highest proportion of variation in genetic diversity (18% and 15% for GWR and OLS, respectively) being explained by June LST. In order to ensure the reliability of the OLS model, the spatial autocorrelation of the residuals (the difference between the observed and fitted values) was analyzed using Moran's I statistic. The results showed no significant spatial autocorrelation, with a value of I = -0.10 (p-value = 0.79) for the OLS model, and I = -0.07 (p-value = 0.85) for the GWR model. This indicates that the model is accurately predicting the relationship between LST and genetic diversity. The mapping of standard residuals (Fig. 6C) does not show unusually high or low residuals, further supporting the reliability of the model.

4.1 Genetic diversity and recent bottleneck

The genetic diversity of Ulmus species has been studied by various researchers. In a microsatellite loci-based study Zalapa et al. (2008) found that the average genetic diversity of U. pumila from 53 locations in China was relatively low (H_e < 0.5) (Zalapa et al., 2008). Similarly, the haplotype diversity of U. lamellosa in northern China was found to be reduced (Hou et al., 2020). Low genetic diversity was also observed in populations of U. parvifolia in eastern China and U. glabra in the Iberian Peninsula using Single Nucleotide Polymorphism (SNPs) and microsatellite markers, respectively (Lyu et al., 2020; Martín del Puerto et al., 2017). In our results H_o is found consistently lesser than H_e and hence negative F_st values. The negative Fixation index values, indicating minimal differentiation have been consistently observed across U. wallichiana populations. Similar results have been observed in some of the populations of other Himalayan species (Majeed et al., 2021; Radosavljević et al., 2015; Wu et al., 2023). This pattern challenges traditional models of genetic differentiation. Continuous habitats, ecological continuity and adaptive uniformity contribute to this pattern. Negative Fst values indicates continuity in population and high gene flow among them. This suggests that any habitat separation that has been noticed during field trips is probably quite recent, and it is expected to result in more population differentiation soon.

In contrast, this study found substantial genetic diversity (mean H_o = 0.634) and moderate differentiation between populations of U. wallichiana (mean F_ST = 0.18). The presence of substantial genetic diversity and the absence of inbreeding suggest that the evolutionary potential of this species is retained, despite being threatened. Therefore, present scenario offers a great chance to protect this species' genetic diversity. Certainly, genetically diverse population (pop2, pop5, pop10, pop11, and pop14) represents valuable genetic resources with greater environment resilience, that can be utilized in species expansion programs into future climatic hotspots. By incorporating germplasm from these diverse populations into conservation initiatives, we can enhance the genetic variability of U. wallichiana populations in regions expected to experience favorable climates in the future. Moreover, it is critical to acknowledge the significance of maintaining climate-specific adaptations that have developed among various populations. The presence of private alleles suggests localized genetic variation, emphasizing the importance of population-specific conservation strategies. Region specific adaptations that enable populations to flourish in various climatic regimes can be preserved by protecting populations in their natural habitats and maintaining the integrity of their genetic make-up. The presence of recent genetic bottleneck signs in this study is a cause for concern. The most frequent explanation for the high frequency of genetic bottleneck is the dispersed populations of this species in the Western Himalayan mountainous terrain. Like this research, it has been observed that a genetic bottleneck has affected most populations of various species in the Himalayan regions (Fan et al., 2018; Majeed et al., 2021; Nag et al., 2015). Due to their remote locations, geographical isolation, and the negative effects of climate change, mountainous species primarily face genetic bottlenecks. Genetic bottlenecks can reduce genetic diversity due to inbreeding (Young AG, 1996). However, the observation of substantial genetic diversity and no signs of inbreeding in U. wallichiana suggests that its population has faced recently contraction to a significant level which causes a pronounced decrease in homozygosity and an increase in heterozygosity. However, a decrease in genetic diversity and an increase in inbreeding may become apparent in natural population in the future (Majeed et al., 2021).

4.2 Gene flow barriers and population structure

Due to the complex topography of the Himalayas, geographic barriers can affect gene flow and may have important implications for conservation strategies. The results of this study revealed the presence of a prominent genetic barrier separating the populations of U. wallichiana. This barrier restricts gene flow between the populations on either side, leading to fragmentation of the gene pool. To mitigate the effects of this barrier, it is suggested that germplasm be propagated around its vicinity, allowing gene flow to continue on either side. Our analysis using STRUCTURE and principal coordinate analysis (PCoA) revealed two genetically distinct populations, possibly determined by the observed genetic barrier. As the geographical distance between populations increases, the genetic distance also increases. A weak relationship was observed between physical distance, altitude, and genetic differentiation, similar to what has been observed in Betula pubescens (Truong et al., 2007). Bioclimatic variables bio7 and bio18 were found to have a significant positive correlation with population differentiation and showed a similarity in population clustering based on these variables to the genetic distance. This suggests that these variables possibly may be involved in isolation by environment (IBE), a phenomenon that occurs when populations are separated due to differences in environmental conditions. Microclimatic variations may cause population differentiation through asynchronous blooming of geographically separated populations (Hirao & Kudo, 2004; Majeed et al., 2021; Perry & Starrett, 1980), which can be reflected through SSR markers.

4.3 Geographical isolation and environmental influence

The ecological niche model for the current and future datasets was found to be particularly sensitive to bio1 (annual mean temperature), bio2 (mean diurnal range), bio3 (isothermality), bio14 (precipitation of driest month), bio18 (precipitation of warmest quarter), and bio19 (precipitation of coldest quarter). Niche modelling results revealed an overall reduction of climatic suitability in the future for U. wallichiana, however, certain zones are expected to receive favorable environments in the future, therefore, can represent potential zones for safe conservation programs. Germplasm from genetically more diverse populations (pop2, pop5, and pop14) can be planted in areas with favorable climates. This approach could enhance the performance and diversity of U. wallichiana in these areas. Further, our analysis showed that the climate in certain areas (pop3, pop14, and pop13) is expected to become less suitable for this species in the future, posing a greater possible threat to this already vulnerable species, which might result in the loss of diversity or even the loss of accessions from such areas. In this study, the genetic diversity was shown to have a modest correlation with altitude, suggesting that there may be adaptive responses to altitudinal gradients. However, further studies regarding altitudinal relationship with genetic diversity are recommended. In view of preservation, it is recommended that the germplasm of U. wallichiana from these climatically degrading areas be propagated in climatically suitable regions. Moreover, ex-situ conservation also plays a crucial role in the conservation of vulnerable species, therefore maintenance of ex-situ conservation banks of these top-performing germplasm is suggested for future use. In conclusion, considering future climate change in the design and implementation of conservation programs, the strategy of germplasm propagation from genetically diverse populations in suitable environments may be effective in maintaining the genetic diversity of this species. The conservation of U. wallichiana requires diligent and dedicated efforts to ensure its continued survival and flourishing. The species, with its unique and irreplaceable presence in nature, must be protected and preserved for the future. Overall, the conservation of U. wallichiana is important for well-being, biodiversity, and environmental benefits, as well as for the cultural and societal benefits it provides.

4.4 Temperature dynamics in shaping genetic diversity

Ulmus wallichiana typically begins blooming in the first week of April, with the pollen season and duration varying across the landscape. Some regions experience earlier pollination and longer pollen seasons than others, which can have important consequences for the genetic diversity of the species. The most important environmental factor influencing the blooming timing and duration of a species is temperature (Fitter & Fitter, 2002). Our analysis of the temporal relationship between land surface temperature (LST) and genetic diversity in U. wallichiana showed that temperature variation in June had the greatest influence on population genetic diversity. The average temperature for blooming is likely optimal from March to June, with temperature variations in June being the most important predictor variable, followed by May and April (Urhausen et al., 2011). However, it should be noted that this is an indirect estimate of the temporal relationship between LST and genetic diversity and further research incorporating blooming and pollination data may be necessary to obtain more robust and accurate results.

Emphasizing the importance and benefits of conserving U. wallichiana as an alternative mitigation to climate change and highlighting future perspectives to avoid species extinction would further strengthen the study's implications and potential contributions to conservation efforts. Himalayan Elm, being a temperate deciduous tree species, plays a crucial role in mitigating climate change and its impacts. Additionally, mature trees of Himalayan Elm can provide important ecosystem services, such as regulating microclimates, reducing soil erosion, and maintaining water quality, all of which are vital in the face of climate change impacts. Furthermore, protecting and restoring U. wallichiana habitats can contribute to preserving overall biodiversity, which is essential for maintaining resilient ecosystems that can adapt to changing climatic conditions.

To mitigate the risk of extinction for U. wallichiana, proactive conservation methods must be undertaken, combining insights from genomic analysis and ecological niche modelling. These strategies should encompass habitat restoration, the creation of protected areas, and targeted replanting projects. The predictive power of ecological niche modelling can help guide conservation efforts by identifying areas that are likely to become climatically suitable for the species in the future. Active engagement with local people is critical to the long-term viability of conservation efforts, involving collaboration across a wide range of stakeholders, including government agencies, non-governmental organizations, researchers, and community members. Regular monitoring of population dynamics and genetic diversity is paramount for assessing the success of conservation measures in response to changing environmental conditions. By incorporating these future perspectives, preventive steps can be implemented to save U. wallichiana from extinction and preserve its genetic variety for future generations. The outcomes of this study highlight the vital necessity of incorporating both genomic and climatic data into conservation planning, emphasizing the need for holistic approaches to ecosystem conservation. Protecting habitats with high genetic diversity and prioritizing conservation efforts in climate change-prone areas are critical steps towards ensuring the long-term survival of U. wallichiana in the face of environmental challenges.

AIC

Akaike Information Criterion

AMOVA

Analysis of Molecular Variance

AUC

Area under the receiver operating characteristics curve5

CTAB

Cetyltrimethylammonium bromide

GBIF

Global Biodiversity Information System Facility

GPS

Global Positioning System

GWR

Geographical Weighted Regression

Himachal Pradesh

IBD

Isolation by distance

IBE

Isolation by environment

IDW

Inverse distance weightage

IUCN

International Union for Conservation of Nature

J&K

Jammu and Kashmir

LST

Land Surface Temperature

MODIS

Moderate Resolution Imaging Spectroradiometer

MODIS

Moderate resolution Imaging Spectroradiometer

NASA LPDAAC

The Land Process Distributed Active Centre

OLS

Ordinary least squares

PCoA

Principal Coordinate Analysis

SSR

Simple Sequence repeats

UPGMA

Unweighted pair group method with arithmetic mean

USGS

United States Geographical Survey

UTM

Universal Transverse Mercator

WGS

World Geodetic System

Funding: This work was financially supported by the National Mission on Himalayan Studies, Ministry of Environment, Forest and Climate Change (NMHS MoEF& CC), India, under the grant NMHS/SG-2016/011.

Acknowledgments

The authors acknowledge the forest department of J&K and HP for permitting the collection of samples. AS acknowledges MoEF&CC for financial support during the tenure of this study.

Author contributions: PB conceptualized and designed the study; and arranged the funding. AS and AM performed sampling and isolated the DNA. AS performed genotyping and data analysis, and wrote the first draft of the manuscript with AM. PB edited and finalized the final draft of the manuscript. All authors read and approved the final version manuscript.

Batool, N., Bibi, Y., & Ilyas, N. (2021). 01. Current status of Ulmus wallichiana: Himalayan endangered Elm. Pure and Applied Biology (PAB), 3(2), 60–65.
Braje, T. J., & Erlandson, J. M. (2013). Human acceleration of animal and plant extinctions: A Late Pleistocene, Holocene, and Anthropocene continuum. Anthropocene, 4, 14–23.
Buiteveld, J., Vanden Broeck, A., Cox, K., & Collin, E. (2016). Human impact on the genetic diversity of Dutch field elm (Ulmus minor) populations in the Netherlands: implications for conservation. Plant Ecology and Evolution, 149(2), 165–176.
Centre, U. L. P. D. A. A. MODIS (Moderate Resolution Imaging Spectroradiometer) LST and emissivity data. https:/earthexplorer.usgs.gov/
Chakraborty, R., & Nei, M. (1977). Bottleneck effects on average heterozygosity and genetic distance with the stepwise mutation model. Evolution, 347–356.
Charlton, M., Fotheringham, S., & Brunsdon, C. (2009). Geographically weighted regression. White paper. National Centre for Geocomputation. National University of Ireland Maynooth, 2.
Chen, I.-C., Hill, J. K., Ohlemüller, R., Roy, D. B., & Thomas, C. D. (2011). Rapid range shifts of species associated with high levels of climate warming. Science, 333(6045), 1024–1026.
Cornuet, J. M., & Luikart, G. (1996). Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics, 144(4), 2001–2014.
De Marco, P., & Nóbrega, C. C. (2018). Evaluating collinearity effects on species distribution models: An approach based on virtual species simulation. PloS one, 13(9), e0202403.
Di Rienzo, A., Peterson, A., Garza, J., Valdes, A., Slatkin, M., & Freimer, N. (1994). Mutational processes of simple-sequence repeat loci in human populations. Proceedings of the National Academy of Sciences, 91(8), 3166–3170.
Doyle, J. J. (1990). Isolation of plant DNA from faesh tissue. Focus, 12, 13–15.
Earl, D. A., & VonHoldt, B. M. (2012). STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation genetics resources, 4, 359–361.
Evanno, G., Regnaut, S., & Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular ecology, 14(8), 2611–2620.
Fan, L., Zheng, H., Milne, R. I., Zhang, L., & Mao, K. (2018). Strong population bottleneck and repeated demographic expansions of Populus adenopoda (Salicaceae) in subtropical China. Annals of Botany, 121(4), 665–679.
Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International journal of climatology, 37(12), 4302–4315.
Fitter, A., & Fitter, R. (2002). Rapid changes in flowering time in British plants. Science, 296(5573), 1689–1691.
Goodall-Copestake, W., Hollingsworth, M. L., Hollingsworth, P. M., Jenkins, G., & Collin, E. (2005). Molecular markers and ex situ conservation of the European elms (Ulmus spp.). Biological Conservation, 122(4), 537–546.
Haanes, H., Røed, K. H., Perez-Espona, S., & Rosef, O. (2011). Low genetic variation support bottlenecks in Scandinavian red deer. European Journal of Wildlife Research, 57, 1137–1150.
Hirao, A., & Kudo, G. (2004). Landscape genetics of alpine-snowbed plants: comparisons along geographic and snowmelt gradients. Heredity, 93(3), 290–298.
Holderegger, R., & Wagner, H. H. (2008). Landscape genetics. Bioscience, 58(3), 199–207.
Hou, H., Ye, H., Wang, Z., Wu, J., Gao, Y., Han, W., Na, D., Sun, G., & Wang, Y. (2020). Demographic history and genetic differentiation of an endemic and endangered Ulmus lamellosa (Ulmus). BMC Plant Biology, 20(1), 1–14.
Howden, S. M., Soussana, J.-F., Tubiello, F. N., Chhetri, N., Dunlop, M., & Meinke, H. (2007). Adapting agriculture to climate change. Proceedings of the National Academy of Sciences, 104(50), 19691–19696.
Huang, L., Deng, X., Li, R., Xia, Y., Bai, G., Siddique, K. H., & Guo, P. (2018). A fast silver staining protocol enabling simple and efficient detection of SSR markers using a non-denaturing polyacrylamide gel. JoVE (Journal of Visualized Experiments)(134), e57192.
IUCN. (1998, 2021-22). The IUCN Red List Threatened species. World Conservation Monitoring Centre. Retrieved 13 October from https://www.iucnredlist.org/species/30410/9545816
Kimura, M., & Ohta, T. (1978). Stepwise mutation model and distribution of allelic frequencies in a finite population. Proceedings of the National Academy of Sciences, 75(6), 2868–2872.
Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., & Mayrose, I. (2015). Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Molecular ecology resources, 15(5), 1179–1191.
Litt, M., & Luty, J. A. (1989). A hypervariable microsatellite revealed by in vitro amplification of a dinucleotide repeat within the cardiac muscle actin gene. American journal of human genetics, 44(3), 397.
Liu, L., Chen, W., Zheng, X., Li, J., Yan, D.-T., Liu, L., Liu, X., & Wang, Y.-L. (2016). Genetic diversity of Ulmus lamellosa by morphological traits and sequence-related amplified polymorphism (SRAP) markers. Biochemical Systematics and Ecology, 66, 272–280.
Luikart, G., Sherwin, W., Steele, B., & Allendorf, F. (1998). Usefulness of molecular markers for detecting population bottlenecks via monitoring genetic change. Molecular ecology, 7(8), 963–974.
Lyu, Y.-z., Dong, X.-y., Huang, L.-b., Zheng, J.-w., He, X.-d., Sun, H.-n., & Jiang, Z.-p. (2020). SLAF-seq uncovers the genetic diversity and adaptation of Chinese Elm (Ulmus parvifolia) in Eastern China. Forests, 11(1), 80.
Majeed, A., Singh, A., & Bhardwaj, P. (2021). Transcribing molecular and climatic data into conservation management for the Himalayan endangered species, Taxus contorta (Griff.). Conservation Genetics, 22, 53–66.
Manne, L. L., & Pimm, S. L. (2001). Beyond eight forms of rarity: which species are threatened and which will be next? Animal conservation forum,
Manni, F., Guérard, E., & Heyer, E. (2004). Geographic patterns of (genetic, morphologic, linguistic) variation: how barriers can be detected by using Monmonier's algorithm. Human biology, 173–190.
Mantel, N. (1967). The detection of disease clustering and a generalized regression approach. Cancer research, 27(2_Part_1), 209–220.
Martín del Puerto, M., Martínez García, F., Mohanty, A., & Martín, J. P. (2017). Genetic diversity in relict and fragmented populations of Ulmus glabra Hudson in the central system of the Iberian Peninsula. Forests, 8(5), 143.
Mateo, R. G., Vanderpoorten, A., Muñoz, J., Laenen, B., & Désamoré, A. (2013). Modeling species distributions from heterogeneous data for the biogeographic regionalization of the European bryophyte flora. PloS one, 8(2), e55648.
Mayr, E. (1963). Animal species and evolution. Harvard University Press.
Meier, E. S., Kienast, F., Pearman, P. B., Svenning, J. C., Thuiller, W., Araújo, M. B., Guisan, A., & Zimmermann, N. E. (2010). Biotic and abiotic variables show little redundancy in explaining tree species distributions. Ecography, 33(6), 1038–1048.
Morales, N. S., Fernández, I. C., & Baca-González, V. (2017). MaxEnt’s parameter configuration and small samples: are we paying attention to recommendations? A systematic review. PeerJ, 5, e3093.
Muscarella, R., Galante, P. J., Soley-Guardia, M., Boria, R. A., Kass, J. M., Uriarte, M., & Anderson, R. P. (2014). ENM eval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods in ecology and evolution, 5(11), 1198–1205.
Nag, A., Ahuja, P. S., & Sharma, R. K. (2015). Genetic diversity of high-elevation populations of an endangered medicinal plant. AoB Plants, 7, plu076.
Nielsen, L. R., & Kjær, E. D. (2010). Fine-scale gene flow and genetic structure in a relic Ulmus laevis population at its northern range. Tree Genetics & Genomes, 6, 643–649.
Peakall, R., & Smouse, P. E. (2006). GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular ecology notes, 6(1), 288–295.
Perry, D. R., & Starrett, A. (1980). The pollination ecology and blooming strategy of a neotropical emergent tree, Dipteryx panamensis. Biotropica, 307–313.
Peterson, A. T., Ortega-Huerta, M. A., Bartley, J., Sánchez-Cordero, V., Soberón, J., Buddemeier, R. H., & Stockwell, D. R. (2002). Future projections for Mexican faunas under global climate change scenarios. Nature, 416(6881), 626–629.
Phillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecological modelling, 190(3–4), 231–259.
Piry, S., Luikart, G., & Cornuet, J.-M. (1999). BOTTLENECK: a program for detecting recent effective population size reductions from allele data frequencies. Journal of Heredity, 90(4), 502–503.
Pritchard, J. K., Stephens, M., & Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics, 155(2), 945–959.
Radosavljević, I., Satovic, Z., & Liber, Z. (2015). Causes and consequences of contrasting genetic structure in sympatrically growing and closely related species. AoB Plants, 7, plv106.
Sastre, N., Vila, C., Salinas, M., Bologov, V. V., Urios, V., Sánchez, A., Francino, O., & Ramírez, O. (2011). Signatures of demographic bottlenecks in European wolf populations. Conservation Genetics, 12, 701–712.
Schneiderbauer, S., Pisa, P. F., Delves, J. L., Pedoth, L., Rufat, S., Erschbamer, M., Thaler, T., Carnelli, F., & Granados-Chahin, S. (2021). Risk perception of climate change and natural hazards in global mountain regions: A critical review. Science of the total environment, 784, 146957.
Singh, A., Majeed, A., & Bhardwaj, P. (2021). Transcriptome characterization and generation of marker resource for Himalayan vulnerable species, Ulmus wallichiana. Molecular biology reports, 48(1), 721–729.
Storfer, A., Murphy, M., Evans, J., Goldberg, C., Robinson, S., Spear, S., Dezzani, R., Delmelle, E., Vierling, L., & Waits, L. (2007). Putting the ‘landscape’in landscape genetics. Heredity, 98(3), 128–142.
Syfert, M. M., Smith, M. J., & Coomes, D. A. (2013). The effects of sampling bias and model complexity on the predictive performance of MaxEnt species distribution models. PloS one, 8(2), e55158.
Telwala, Y., Brook, B. W., Manish, K., & Pandit, M. K. (2013). Climate-induced elevational range shifts and increase in plant species richness in a Himalayan biodiversity epicentre. PloS one, 8(2), e57103.
Truong, C., Palmé, A. E., & Felber, F. (2007). Recent invasion of the mountain birch Betula pubescens ssp. tortuosa above the treeline due to climate change: genetic and ecological study in northern Sweden. Journal of evolutionary biology, 20(1), 369–380.
Urhausen, S., Brienen, S., Kapala, A., & Simmer, C. (2011). Climatic conditions and their impact on viticulture in the Upper Moselle region. Climatic Change, 109(3–4), 349–373.
USGS. Digital elevation model (DEM) data. https://earthexplorer.usgs.gov/
Vakkari, P., Rusanen, M., & Kärkkäinen, K. (2009). High genetic differentiation in marginal populations of European white elm (Ulmus laevis).
Venturas, M., Fuentes-Utrilla, P., Ennos, R., Collada, C., & Gil, L. (2013). Human-induced changes on fine-scale genetic structure in Ulmus laevis Pallas wetland forests at its SW distribution limit. Plant Ecology, 214, 317–327.
Wu, Q., Dong, S., Zhao, Y., Yang, L., Qi, X., Ren, Z., Dong, S., & Cheng, J. (2023). Genetic diversity, population genetic structure and gene flow in the rare and endangered wild plant Cypripedium macranthos revealed by genotyping-by-sequencing. BMC Plant Biology, 23(1), 254.
Young AG, M. D., Robertson A. (1996). Genetic bottlenecks and population management. Nature, 373(6514), 215–220.
Zalapa, J. E., Brunet, J., & Guries, R. P. (2008). Genetic diversity and relationships among Dutch elm disease tolerant Ulmus pumila L. accessions from China. Genome, 51(7), 492–500.
Zhu, G.-P., & Peterson, A. T. (2017). Do consensus models outperform individual models? Transferability evaluations of diverse modeling approaches for an invasive moth. Biological Invasions, 19, 2519–2532.
Zimmerer, K. S., de Haan, S., Jones, A. D., Creed-Kanashiro, H., Tello, M., Carrasco, M., Meza, K., Amaya, F. P., Cruz-Garcia, G. S., & Tubbeh, R. (2019). The biodiversity of food and agriculture (Agrobiodiversity) in the anthropocene: Research advances and conceptual framework. Anthropocene, 25, 100192.

No competing interests reported.

Download PDF

Editorial decision: Revision requested
18 Nov, 2024
Reviews received at journal
15 Nov, 2024
Reviewers agreed at journal
01 Nov, 2024
Reviews received at journal
16 Oct, 2024
Reviewers agreed at journal
07 Oct, 2024
Reviewers invited by journal
07 Oct, 2024
Editor assigned by journal
30 Aug, 2024
Submission checks completed at journal
30 Aug, 2024
First submitted to journal
30 Aug, 2024

You are reading this latest preprint version

Optimizing conservation of Ulmus wallichiana, threatened Himalayan species by combining population genetics and species distribution modeling

Status:

Version 1

Abstract

Figures

1. Introduction

2. Materials and Methods

2.1. Plant Material

2.2. Genetic diversity

2.3. Population connectivity

2.4. Population differentiation and structure

2.5. Assessment of isolation by distance (IBD) and isolation by environment (IBE)

2.6. Bottleneck analysis

2.7. Ecological niche modelling

2.8. Modelling the temperature determinants of genetic diversity in the study area

3. Results

3.1. Genetic diversity, barriers, and recent bottleneck

3.2. Population differentiation and structure

3.3. Isolation by distance (IBD), isolation by environment (IBE), and altitudinal impact on population differentiation

3.4. Ecological niche modelling

3.5. Modelling the temperature determinants of genetic diversity in the study area

4. Discussion

4.1 Genetic diversity and recent bottleneck

4.2 Gene flow barriers and population structure

4.3 Geographical isolation and environmental influence

4.4 Temperature dynamics in shaping genetic diversity

5. Future perspectives and Conclusions

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1