Study area
The Himalaya a global biodiversity hotspot represents the world's premier mountain system on earth and is thereby known to possess varied ecological, biogeographical, and evolutionary attributes giving rise to tremendous biodiversity of endemic origin (Joshi and Pant 2012; White et al., 2019; Rana et al., 2021). The Himalayan biodiversity hotspot falls within the geo-coordinates (34–36°N, 27°E) to (27–28° N, 90°E) involving north-western Pakistan and Northern Myanmar respectively, and is snuggled between the south Indian subcontinent and the northern highland of Tibetan (https://www.cepf.net) (Fig. 1). It has been divided into two categories biogeographically such as Western Himalaya including regions of Jammu and Kashmir, Uttarakhand and Himachal Pradesh while Eastern Himalaya comprises of north-eastern Indian states (Sikkim, Assam, West Bengal, and Arunachal Pradesh), southeast Tibet, northern Myanmar, Nepal, and Bhutan (Hamid et al., 2019; Rather et al., 2022; Wani et al., 2022). It is located at a biogeographic crossroads amidst the Palearctic and Indo-Malayan regions. Furthermore, the Himalayas' geological, climatic, and altitudinal differences, besides complexity in topography (Troll 1972; Zurick and Pacheco 2006; Miehe et al. 2015), bestow Himalaya with abundant biodiversity.
Currently, Himalayan biodiversity hotspot is facing increased anthropogenic influences such as land use change, habitat loss, spreading of invasive alien species, and non-judicious ways of gathering natural resources which ultimately led to biodiversity loss, deprivation of fragile ecosystems as well as compromising ecosystem services and products (Ren et al., 2007; Khuroo et al., 2012; Roy and Srivastava, 2012; Palmer et al., 2014; Ahmad et al., 2019b; Uprety et al., 2021)
Target species
Aquilegia fragrans Benth, an endangered perennial medicinally important herb belonging to the family Ranunculaceae, are distributed predominantly in the Western Himalaya (Hussain et al., 2011; Baig et al., 2014). Aquilegia fragrans, prevalently known as “sweet scented columbine” and in Kashmir locally possess different names “Janglikuth”, “Kalumb” and “Zaoneel”. The plant grows naturally in India (Kashmir valley) at elevated altitudes (Mathew and Sinnott, 2003). Its stem is 30–80 cm in length, branched, densely hairy with glands below the flower. Leaflets are wedge-shaped to obovate and glaucous. Flowers are several, somewhat horizontal to slightly nodding. Sepal’s 20–30×9–12 mm, whitish to purple, slightly velvety (Fig. 2). This plant finds its traditional use in treating diverse inflammatory ailments and wounds (Baig et al., 2014). Besides, the rhizome is found to be curative in kidney stones, eczema, psoriasis, blood sugar, cystitis, gout, and bovine mastitis (Ganie et al., 2019, Rani et al., 2013). For headache alleviation, dried flower paste is utilized, and jaundice and pneumonia are treated traditionally with seeds and the whole plant respectively (Khan et al., 2018). Recently its root has been found effective ethanomedicinally as an analgesic and antipyretic (Majid et al., 2020). For the treatment of boils plant flower paste in refined butter is made and subsequently massaged for one to two weeks (Lal and Singh, 2008). Overexploitation and habitat fragmentation due to the immense medicinal properties of A. fragrans has put enormous strain on the population of this endangered species in the Himalaya.
Species distributional records
A Comprehensive field survey was done from (2019–2021) in the entire Kashmir Himalaya to collect occurrence records of Aquilegia fragrans, which are supplemented with geographical records from herbarium reconstruction (KASH) and some of the well-known databases like BIEN, IBP, GBIF, and BBP (Botanical information Facility, Indian Biodiversity Portal, Global Biodiversity Information Facility and Bhutan Biodiversity Portal respectively). All the above-mentioned databases were accessed on 30 February 2022.
A total of 75 georeferenced records for its distribution were gathered from the aforementioned sources (S1) and subsequently reduced to 32 after clipping for the research area (Himalaya). The species distributional records are frequently oversampled from those areas which are physically suited and easily accessible (Hijmans et al., 2005). This led to geographic sample bias (Ahmad et al. 2019a; Hamid et al., 2019; Rather et al., 2022). We employed spatial thinning to reduce geographical autocorrelation and sample bias by dividing the research region into grid cells of 1×1 km and arbitrarily picking a lone point among every cell having single or multiple points of occurrence (Ahmad et al. 2019a; Hamid et al., 2019; Rather et al., 2022). The ultimate dataset is taken for modelling of A. fragrans comprised of 30 geo-referenced records
Environmental data
To forecast the present distribution pattern of plant species under study across the research area, WorldClim database version 1.4 was utilized for the extraction of climatic variables at a geographical resolution of 30 arc seconds (Hijmans et al., 2005). Additionally, greater collinearity across bioclimatic variables may reduce model performance and can produce unreliable interpretations (Dormann et al., 2013). To solve this issue Pearson's correlation analysis was taken into account wherein we utilized only one variable from each pair of a highly interrelated variable having a correlation coefficient of (r > 0.75). Following correlation analysis, five variables were selected for the distribution modelling of A. fragrans under current climatic conditions (Table 1). For the prediction of distribution in future with respect to selected species, Hadley Global Environment Model 2-Earth System (HADGEM2-ES) simulations for two representative concentration pathways (RCP4.5 and RCP8.5) for time periods 2050 and 2070 were obtained from the Intergovernmental Panel on Climate Change (IPCC) fifth assessment report (AR5) (Moss et al., 2010). Future distributions of the target species were predicted using the same set of environmental variables that were used to estimate their current distribution.
Table 1
Bioclimatic variable for modelling the distribution of A. fragrans in the present study. Selected variables are in bold.
BIO-1 (Annual Mean Temperature) |
BIO-2 (Mean Diurnal Range) |
BIO-3 (Isothermality) |
BIO-4 (Temperature Seasonality) |
BIO-5 (Max Temperature of Warmest Month) |
BIO-6 (Min Temperature of Coldest Month) |
BIO-7 (Temperature Annual Range) |
BIO-8 (Mean Temperature of Wettest Quarter) |
BIO-9 (Mean Temperature of Driest Quarter) |
BIO-10 (Mean Temperature of Warmest Quarter) |
BIO-11 (Mean Temperature of Coldest Quarter) |
BIO-12 (Annual Precipitation) |
BIO-13 (Precipitation of Wettest Month) |
BIO-14 (Precipitation of Driest Month) |
BIO-15 (Precipitation Seasonality) |
BIO-16 (Precipitation of Wettest Quarter) |
BIO-17 (Precipitation of Driest Quarter) |
BIO-18 (Precipitation of Warmest Quarter) |
BIO-19 (Precipitation of Coldest Quarter) |
Modelling technique
During the present study species distribution modelling and ensemble forecasting was executed by utilising “biomod2” package (Thuiller et al., 2009) within R software version 4.0.3 (R Development Core Team, 2021). A total of nine algorithms were utilised (Table 2) in the current study. Presence/absence data is required for modelling but acquiring accurate absence data is challenging. Subsequently, in accordance with Barbet-Massin et al. (2012), Guisan et al. (2017), and Rather et al., (2022) we generated five hundred pseudo-absences at random within the study area. Since our pseudo-absence selection procedure uses a stochastic approach that selects pseudo-absences at random, we carried out the operation 3 times for removing bias in the sample created during the production of pseudo-absence (Guisan et al., 2017; Rather et al., 2022). We used 80% of the data (training set) to develop the models and evaluated on the remaining 20% (validation set). The whole approach was repeated thrice. As a consequence, we got a total of 81 models for each climate situation and time period combination (3 replicate pseudo-absence datasets ×9 algorithms × 3 replicates) (Rather et al., 2022). Two assessment metrics, such as area under the curve (AUC) of receiver operating characteristics and true skills statistics (TSS), were utilised for examining performance of model using a repeated data-splitting technique called cross-validation (Allouche et al., 2006; Marmion et al., 2009; Beaumont et al., 2016; Rather et al., 2022).
Using committee averaging as well as weighted-mean approach, the final ensemble model was constructed for each climate scenario along with time period combination based on various modelling results (Marmion et al., 2009). Those models having TSS scores of ≥ 0.8 were retained for the creation of final ensemble models. Consequently, five ensemble projections were created, one reciprocates with current climatic suitability and the remaining ones for anticipation of future prediction on habitat suitability for the time periods of 2050 and 2070 periods), representing 2 typical concentration trajectories RCP4.5 and 8.5 (Kumari et al., 2022, Rather et al., 2022).
Table 2
List of algorithms used in the present study
S. No | Model | Type | Reference |
1 | ANN Classification Tree Analysis (CTA) | Machine learning | (Breiman et al., 1984), |
2 | Generalised Linear Model (GLM) | Machine learning | (McCullagh and Nelder, 1989) |
3 | Surface Response Envelope (SRE) | Profile | (Busby, 1991) |
4 | Flexible Discriminant Analysis (FDA) | Statistical regression | (Hastie et al., 1994), |
5 | Artificial Neural Networks (ANN) | Machine learning | (Ripley, 1996) |
6 | Generalised Boosted Models (GBM) | Machine learning | (Ridgeway, 1999) |
7 | Random Forest (RF) | Machine learning | (Breiman, 2001) |
8 | Maximum Entropy (MAXENT) | Machine learning | (Phillips et al.,2006) |
9 | Generalised Additive Models (GAM) | Statistical regression | (Hastie and Tibshirani, 1990) |
Variable importance
We employed a permutation method to determine the relative impact of each climatic condition in affecting the distribution of A. fragrans (Elith et al., 2005, Rather et al., 2022). Predictions were generated during this process by using a predetermined algorithm after changing just a solitary target variable whilst holding other variables constant. The significance of a particular variable is estimated by subtracting the score of correlation (1-correlation score) from the difference between the original forecast and the permuted variable prediction (Rather et al., 2022; Ahmad et al., 2019a). Correspondingly, a zero value indicates, the predictor variable holds no connotation toward the model, whereas higher values indicate a predictor variable of more importance.
Species Range Change
In the biomod2 R package RangeSize () function was employed to measure and determine the range change shift for A. fragrans under projected future climate change (Thuiller et al., 2009). We get two outputs from the RangeSize() function: a table with summarized information on species range shifts and a geographical map showing whether species may gain or lose appropriate habitats (Guisan et al., 2017; Rather et al. 2022; Kumari et al., 2022).
Taking these output types into consideration one will come to know about absolute metrics such as "loss", "absent" and "stable". So far "loss" is concerned it denotes pixel in number forfeited by species under study owing to the fact of climate alterations; "absent" represents pixel resolution in number presently not engrossed through species under study besides having no prediction for suitability under a peculiar environmental scenario; and “stable," represents a pixel in number engrossed via species under research and additionally predicted to be appropriate beneath a particular climatic scenario(Guisan et al., 2017; Rather et al. 2022; Kumari et al., 2022).
We also get three additional relative metrics like: "percentage loss" indicating a loss in the percentage of sites inhabited presently and measured via (loss/ (loss + stable); "percentage gain" which is figured out through (gain/(loss + stable); simply can be put as a percentage of new sites taking into account the present distribution size of species and "range change" expresses comprehensive projection output equivalent to percentage gain– percentage loss (Guisan et al., 2017; Rather et al. 2022).