Global change on the roof of the world: Vulnerability of Himalayan otter species to land use and climate alterations

Climate Change Vulnerability Assessment (CCVA) prescribes the quantification of species vulnerability based on three components: sensitivity, adaptive capacity and exposure. Such assessments should be performed through combined approaches that integrate trait‐based elements (e.g., measures of species sensitivity such as niche width) with correlative tools quantifying exposure (magnitude of changes in climate within species habitat). Furthermore, as land use alterations may increase climate impacts on biodiversity, CCVAs should focus on both climate and land use change effects. Unfortunately, most of such assessments have so far focused exclusively on exposure to climate change.


| INTRODUC TI ON
Climate change has recently emerged as a key threat to biodiversity and ecosystems (Chen et al., 2011;Pacifici et al., 2015). In fact, it is expected to surpass habitat loss as the prime threat to global biodiversity in the future (Leadley, 2010). Besides climate alterations, land use change is among the most critical drivers of global change and biodiversity loss (Lambin, 1999) as it influences many aspects of geo-environmental systems (Guzha et al., 2018;Kiruki et al., 2017;Rogora et al., 2018), hampering the capacity of species to move to suitable habitats (Oliver et al., 2015;Pereira et al., 2012). Despite that, the majority of studies exploring species' responses to global change have largely disregarded the role of land use characteristics, assuming that species are limited only by the shift of climatic conditions (Pacifici et al., 2015;Sirami et al., 2017;but see Di Febbraro et al., 2019). Overall, our current knowledge of the implications of climate and land use change interactions on ecological systems is limited and mainly founded on broad assumptions rather than on observed interactions (Brook et al., 2008;Oliver & Morecroft, 2014).
As land use change is expected to increase climate change impacts on biodiversity, it is critical to integrate these two components to understand likely consequences of such pressures on biodiversity (Sirami et al., 2017;Titeux et al., 2017).
From this perspective, an effective assessment of global change threats to biodiversity relies on identifying highly vulnerable species, together with the drivers of such vulnerability. One of the most comprehensive frameworks developed to assess species vulnerability to climate change is the "Climate Change Vulnerability Assessment" (CCVA; Foden et al., 2019). This analytical framework uses an array of different approaches (i.e., biological traits, correlative and mechanistic models) to quantify species vulnerability based on three main components: (a) sensitivity (i.e., degree to which a species is affected by climate change as a consequence of intrinsic attributes such as niche specialization, environmental tolerance, etc.); (b) adaptive capacity (i.e., ability to adjust to climate change by virtue of high levels of phenotypic plasticity, dispersal ability or genetic diversity) and; (c) exposure (i.e., rate and magnitude of climate changes within species habitat; Foden et al., 2019). Despite CCVAs having been widely used in recent years, most of these assessments have focused solely on the exposure component (Thompson et al., 2015), usually quantified through correlative approaches (i.e., Species Distribution Models, "SDMs" hereafter; Pacifici et al., 2015). On the other hand, assessments involving a strictly trait-based approach (i.e., using species' biological characteristics, such as niche width or dispersal capabilities, as predictors of climate-driven extinction risk), mostly quantified species vulnerability according to direct observations, published literature or expert judgment (Borges et al., 2019), though precise vulnerability thresholds associated with each trait are often unknown (Pacifici et al., 2015). Several authors underlined the advantages of using combined approaches able to assess vulnerability by integrating for example, correlative and trait-based elements (Foden et al., 2019;Pacifici et al., 2015). Yet, such compound frameworks have thus far only been used in a limited number of studies (but see e.g., Garcia et al., 2014;Warren et al., 2018), especially if compared to those based solely on correlative models (Pacifici et al., 2015).
Along with this general recommendation towards the use of combined approaches, a deeper understanding of how climate change and other pressures, such as land use alterations, could interact to affect the future survival of species was highlighted as another major challenge, as well as a better integration of this interplay in CCVA techniques (Foden et al., 2019).
Vulnerability assessments are particularly crucial in mountain environments, which are amongst the most sensitive ecosystems to global change (Beniston, 2003). Temperatures in the Himalaya, for instance, are increasing at twice the average warming rate of the northern hemisphere (Chen et al., 2009). Increasing climate change impacts, followed by rising resource demands of growing populations, cause serious concerns about the future of Himalayan biodiversity (Chakraborty et al., 2018). Under the current combination of accelerating global warming and habitat degradation and loss, many Himalayan species may not be able to adapt to such overwhelming pressures (Xu et al., 2009).
In mountain ecosystems, global climate change can exert its strongest effects on freshwater habitats (e.g., increased extreme flooding and erosion; Cianfrani et al., 2018). These systems are already heavily affected by anthropogenic factors, such as pollution, dam construction and habitat disturbance. Otters are highly vulnerable species, sensitive to habitat loss and water pollution (Kruuk, 2006;Ruiz-Olmo et al., 1998). As keystone aquatic predators, otters are considered as indicator species of healthy freshwater ecosystems (Kruuk, 2006;Ruiz-Olmo et al., 1998). The Himalaya encompasses the altitudinal range edge of three otter species ( Figure S1): the Small-clawed otter (Aonyx cinereus), the Eurasian otter (Lutra lutra) and the Smoothcoated otter (Lutrogale perspicillata; Duplaix & Savage, 2018). The IUCN Red List classifies L. lutra as Near Threatened (NT), and A. cinereus and L. perspicillata as Vulnerable (VU; Duplaix & Savage, 2018; see also Appendix S1). Effects of climate change on the three species have been modelled at large scale by Cianfrani et al. (2018), combining species sensitivity, habitat exposure and other metrics to calculate a global vulnerability index. Also, land use and anthropogenic disturbance were considered as static factors by Cianfrani et al. (2011) to evaluate the effects of climate change on the future distribution of the Eurasian otter in Europe. Investigating how these species could respond to future climate and land use change in a highly vulnerable area such as the Himalaya represents a crucial task for delineating effective conservation strategies against global change impacts. In K E Y W O R D S climate change, climate change vulnerability assessment, climate niche factor analysis, habitat exposure, land use change, niche width, SDMs, species distribution models addition, this particular context represents a unique opportunity to test innovative, combined approaches to assess species vulnerability to multiple global change drivers (i.e., climate and land use), in keeping with the prevailing future perspectives delineated for CCVAs (Foden et al., 2019).
In this scenario, the present study aims to explore determinants of otter species vulnerability to future climate and land use change in the Himalaya in terms of species sensitivity (i.e., a measure of niche specialization) and habitat exposure (i.e., dissimilarity between current and future environmental conditions within species habitat).
Specifically, we formulated two hypotheses: (i) including land use change in vulnerability assessments yields higher vulnerability values than assessments relying only on climate change; (ii) species sensitivity exerts a significantly higher contribution than habitat exposure in determining vulnerability values. To test these hypotheses, we have adopted the recently proposed "CNFA" (Climate Niche Factor Analysis) framework (Rinnan & Lawler, 2019) combined with stateof-the-art SDMs. In particular, we explicitly integrated SDMs into the "CNFA" framework by calculating species sensitivity and habitat exposure within the presence/absence maps predicted by SDMs.
Specific objectives of this study were: (a) to model the current potential distribution of the three otter species in the Himalaya, considering different combinations of topographic, climatic and land use variables, (b) to predict their potential distribution in 2050 under climate and land use change scenarios, (c) to assess the percentage of gain/loss in range area and the degree of shift experienced by the three species under such scenarios, (d) to assess species vulnerability to these global change scenarios by coupling SDMs and the "CNFA" framework and (e) to quantify the relative contribution of species sensitivity and habitat exposure in shaping vulnerability.

| Analytical framework
To evaluate the effect of climate and land use change on the three species of interest, we modelled their distribution considering present-day environmental conditions under two alternative setups: (i) including only topographic and climate predictors (climate SDMs), and (ii) adding also land use variables (full SDMs). Subsequently, we projected climate SDMs under two 2050 climate change scenarios (IPCC, 2013), and full SDMs by including two additional 2050 land use change scenarios (Li et al., 2017). Predictions from both climate and full SDMs were then used alternatively to quantify range net change and shift, as well as included into the "CNFA" framework to calculate species sensitivity, habitat exposure and vulnerability.

| Study area
We set our study in the Himalayan region encompassing the following Asian countries: Bhutan, Nepal and some parts of India and Myanmar.
The Himalaya is a vast mountain region located between the Indian Subcontinent in the south and the Tibetan Highland in the north. The region extends from Afghanistan in the northwest (36°N and 70°E) to China in the southeast (26°N and 100°E), covering an area of more than 1,000,000 km 2 with a length of 3,000 km and a maximum width of 400 km. The climate ranges from a tropical zone in the lowlands of India and Nepal, to permanent ice and snow at the highest elevations, and from a more continental regime in the northwest to a more oceanic regime in the southeast. Its abruptly rising mountains result in a high diversity of ecosystems and species, characterizing this area as a unique biodiversity hotspot (Xu et al., 2019).

| Species occurrence data
Occurrence records were gathered from the literature and from online databases, as well as published technical reports (full data source list is provided in Appendix S2). We cross-checked all the occurrences excluding those with a minimum geographical precision coarser than 1 km and collected before 1970, as well as those referring to captive specimens (Cianfrani et al., 2018). Moreover, we drop any occurrence falling farther than 1 km from the hydrographical network (as provided by Lehner et al., 2008). Although opportunistic records may yield accurate predictions of species distribution (Tiago et al., 2017), they are often spatially autocorrelated and/or discontinuous due to a usually unknown and unbalanced sampling effort that can vary noticeably across space (van Strien et al., 2013).
Hence, the occurrence dataset for SDMs calibration, which originally included 229 records among the three species, underwent a further filtering procedure to remove clustered data (Appendix S3). After this step, we obtained a final dataset of 12 records for A. cinereus, 59 records for L. lutra and 43 records for L. perspicillata ( Figure 1).

| Environmental predictors
As environmental variables, we considered the 19 Worldclim predictors, along with two topographical and seven land use variables, all rasterized at a resolution of about 1 km. We derived land use variables from the GeoSOS global database (Li et al., 2017; http://geosi mulat ion.cn/Globa lLUCC Produ ct.html), which describes remotely sensed land use categories as detected by MODIS in 2010 (product MCD12Q1; https://lpdaac.usgs.gov/). In particular, for each pixel we calculated the Euclidean distance from each of these land-use categories. In addition, we calculated the alpha diversity (i.e., count) of land use categories by applying a moving window approach. Specifically, the count of each category was calculated within a 10-km radius circular window moved throughout each pixel of the land use map.
Such radius length was coherent with the average home range size of the three species (Kruuk, 2006). Topographical variables included elevation (gathered from Jarvis et al., 2008) and slope (derived from the elevation map). Predictors were further subselected by checking for multicollinearity (VIF ≤5), retaining the following 13 covariates: elevation, mean diurnal temperature range, temperature seasonality, precipitation of the driest month, precipitation of the warmest quarter, precipitation of the coldest quarter, diversity of land use categories and Euclidean distance from farmlands, forests, grasslands, urban areas, barren lands and water bodies (Table S1). As specified above, elevation and climate variables were used to calibrate climate SDMs, while the entire predictors set was used for full SDMs.

| Species distribution models
To calibrate climate and full SDMs, we applied the so-called "ensemble of small models" approach to avoid model overfitting (Breiner et al., 2015), an issue arising when few records are available with respect to the number of predictors (Guisan & Zimmermann, 2000).
Accordingly, we trained a set of bivariate models for each species, that is, considering all possible combinations of the environmental covariates taken two at a time (10 combinations per species in climate SDMs and 78 combinations per species in full SDMs), and then averaging each model's results. Bivariate models were calibrated by using two algorithms: generalized linear models (GLMs) and maximum entropy models (Maxent; Phillips et al., 2006). We chose these two algorithms as they reported the highest predictive accuracies within an ensemble of small models context (Breiner et al., 2015(Breiner et al., , 2018. For each species, we randomly placed a set of 10,000 background points over a region identified by all the WWF Terrestrial Ecoregions (Olson et al., 2001) intersecting the species records gathered in the Himalayan study area ( Figure 1; Barve et al., 2011). To evaluate predictive performance of the models, we performed a "block" cross-validation (Muscarella et al., 2014) by splitting data into four geographically non-overlapping bins of equal numbers of occurrences, corresponding to each corner of the entire geographical space. Occurrences falling into three out of four bins were used in turn for model calibration, while the data from the held-out bin were used for evaluation purposes. This method proved useful in assessing model transferability, that is, the ability to extrapolate predictions into new areas/times (Roberts et al., 2017), as well as to penalize models based on biologically meaningless predictors (Fourcade et al., 2018). The predictive performance of each model was assessed by measuring the area under the receiver operating characteristic curve (AUC; Hanley & McNeil, 1982) and the true skill statistic (TSS; Allouche et al., 2006). Models were optimally tuned through an intensive procedure as implemented by Breiner et al. (2018). Specifically, for each bivariate model, we varied parameters and complexity as generated by alternative settings and then chose the configuration that yielded the highest AUC under the block cross-validation scheme (Breiner et al., 2018). For GLMs, we tested the shape of the relationship (i.e., linear, quadratic or cubic) and the interaction level (absent or present). For Maxent, we evaluated regularization values between 0.5 F I G U R E 1 Study area, species occurrence data and ecoregions considered for SDMs calibration (green: A. cinereus; blue: L. lutra; red: L. perpsicillata; purple polygons indicate overlapping ecoregions by A. cinereus and L. lutra). The Digital Elevation Model (grey shadowed) and the hydrographic network (yellow) are reported as background layers and 4, with 0.5 steps, and tested the following feature classes: linear, linear + quadratic, hinge, linear + quadratic + hinge, linear + quadratic + hinge + product and linear + quadratic + hinge + product + threshold (Fourcade et al., 2018). The settings configurations that emerged as the optimal ones for most of the tested bivariate models were then used to calibrate final GLMs and Maxent bivariate models (Breiner et al., 2018). Among these, we dropped poorly calibrated models (i.e., achieving AUC <0.7; Di Febbraro et al., 2019) from the subsequent analyses. Ensemble models were obtained by averaging individual GLMs and Maxent models projections weighted by their respective AUC scores (Marmion et al., 2009). For each species, we calculated the relative importance of predictors from the ensemble model using the functionality provided in the "biomod2" package . Lastly, we calculated the spatial autocorrelation in ensemble models residuals through Moran's I correlograms (Appendix S4).
Differences in predictive performance between climate and full SDMs were assessed through t-test.

| Future scenarios
Climate SDMs were projected on two climate change scenarios (i.e., SDMs projections were binarized to obtain range maps according to four threshold approaches (i.e., "equalize sensitivity and specificity," "maximize TSS," "mean occurrence probability" and "minimize receiver operating characteristic plot distance"; Liu et al., 2013), as to account for the effect of using different binarization schemes (Di

| Quantification of the effects of climate and land use change on otter distribution and vulnerability
For each species, model and scenario, we calculated two metrics of climate and land use change effect on binary range maps: net change (in terms of gain/loss percentage) and geographical shift (expressed as 100-percentage overlap between current and future range maps; Franklin et al., 2013). We also assessed the direction of geographical shift by measuring changes in the position of range centroids (Aguirre-Gutiérrez et al., 2016) and margins (Taheri et al., 2021) along longitude, latitude and elevation between current and 2050 range maps (further details are provided in Appendix S6). Species vulnerability to climate and land use change was calculated by following the above mentioned CNFA approach, which was recently implemented by Rinnan and Lawler (2019). This approach relies on ecological niche factor analysis (ENFA, a technique that calculates the environmental niche of a species using presence-only data; Hirzel et al., 2002), to derive two metrics: (i) species sensitivity and (ii) habitat exposure to global change. In particular, sensitivity is calculated as the array of environmental conditions a species tolerates (the so-called specialization vector within the ENFA framework) with respect to a given environmental variable. The rationale behind this metric is that the smaller the range of environmental conditions a species tolerates for a given variable, the smaller its niche width and, consequently, the more sensitive the species will be to strong deviations from these optimal conditions (e. where s and e, respectively, indicate species sensitivity and habitat exposure to the j-th environmental variable. According to the CNFA framework, both sensitivity and exposure were measured within the current distribution of the species (Rinnan & Lawler, 2019). Hence, we alternatively used binary predictions by climate and full SDMs to delineate spatially such current habitat. Specifically, we calculated a different vulnerability value for each binary map deriving from each scenario, global circulation model and binarization threshold. All the vulnerability analyses were carried out using the "CENFA" R package (Rinnan & Lawler, 2019). Differences in vulnerability values between each species pair were assessed through analysis of variance and Tukey post-hoc test. To evaluate the relative contribution of species sensitivity and habitat exposure in determining species vulnerability to climate and land use future scenarios, we regressed vulnerability scores against sensitivity and exposure, respectively, considering each environmental variable. For this purpose, we used linear mixed models (LMM), setting sensitivity and exposure in turn as fixed terms and global circulation model, global change scenario and binarization threshold as random intercepts. LMMs were fitted using the "lme4" R package (Bates et al., 2015).
Specifically, while full SDMs achieved good predictive performances sensu Swets (1988) and Landis and Koch (1977), with AUC values >0.80 and TSS values >0.40 for all species (Table S3), climate SDMs were less accurate, with just one species (i.e., L. perpsicillata) reporting an AUC >0.80 and two (i.e., A. cinereus and L. perpsicillata) achieving TSS values >0.40. In light of that, we hereafter showed and discussed response curves and spatially explicit predictions only for full SDMs ( Figure S3-S10). According to these models, both climate and land use variables were listed amongst the most important predictors shaping the distribution of the three species ( Figure   S4). For example, distance from forests was among the three most contributing predictors for all the species (along with distance from grasslands in L. lutra), showing an inverse relationship with their suitability. Moreover, precipitation during the warmest-driest periods was maximally important in L. perspicillata, and within the four most important predictors for A. cinereus and L. lutra, overall indicating moderate precipitation conditions as the most preferred ( Figure S4).
Elevation was highly important for L. lutra (i.e., highest suitability at 1,000-3,500 m asl) and, to a lesser extent, for A. cinereus (i.e., suitability peak around 1,000 m asl), while predictors describing temperature fluctuations (i.e., mean diurnal range and temperature seasonality) were among the five most contributing variables for all the three species, showing the otters to prefer a rather high thermal stability ( Figure S4). Full SDMs for all the three species exhibited a negligible spatial autocorrelation in their residuals (Appendix S4), as well as a minimal extrapolation effect in models predictions (see Appendix S5 and Figures S7-S9).

| Quantification of the effects of climate and land use change on otter distribution and vulnerability
Climate SDMs, by accounting for climate change scenarios only, predicted far less severe alterations of otter distribution than full SDMs, which also involve land use change scenarios (Figure 2  showed a moderate sensitivity to temperature seasonality and distance from forests ( Figure 5). In contrast to sensitivity, the three species exhibited habitat exposure scores of comparable magnitude among each other. In fact, their current habitat is similarly exposed to alterations in mean diurnal range, precipitation of the warmest quarter and distance from grasslands ( Figure 5). Taking together the differential contribution of sensitivity and exposure among global circulation models, global change scenarios and binarization thresholds, LMMs outcomes highlighted sensitivity to exert a higher influence than habitat exposure on vulnerability scores. Although statistically significant for all the species, the preponderant role of sensitivity emerged most clearly for L. perspicillata ( Figure 6). On the other hand, the statistical relationship between habitat exposure and vulnerability is overall weaker, and not even significant for L. perspicillata ( Figure 6).   Table S1 Granado -Lorencio, 1996;Prenda et al., 2001). Moreover, Cianfrani et al. (2011, 2018 carried out modelling studies reporting precipitation during the driest month among the most important drivers affecting L. lutra geographical distribution. The relevance of this predictor is likely linked to its effects on hydrological regimes and river characteristics, which in turn affect L. lutra habitat availability, especially during the dry season (Kruuk, 2006;Scorpio et al., 2016).

| D ISCUSS I ON
Overall, elevation also contributed to shape otter distribution, especially for L. lutra. While this species exhibited moderate to high suitability values along a wide altitudinal gradient, confirming its environmental flexibility (Kruuk, 2006;Remonti et al., 2009), A. cinereus suitability sharply decreased toward high elevations (Duplaix & Savage, 2018), although we cannot entirely exclude such pattern to be altered by a lower sampling effort at higher elevation, or by competitive exclusion possibly exerted by L. lutra. Raha and Hussain (2016), exploring the habitat use of the three species in Western Ghats (i.e., a tropical plain and mountain area in southwestern India), found L. perspicillata and L. lutra to prefer dense forested areas at low to intermediate elevations, while A. cinereus tends to occur in waterbodies at higher elevations surrounded by dense grassland. In fact, such evidence is only partly coherent with our findings, though it is important to remark that our study context (i.e., the Himalayan chain) is far different from the area investigated by Raha and Hussain (2016). In addition, such partial inconsistencies may depend upon the great variability of habitat types the three species are known to use throughout their range (Duplaix & Savage, 2018).
Interestingly, our predictions showed land use change to remarkably override the future role played by climate change. In fact, A2 land use change scenario was systematically associated with the severest range contraction and shift in all the three species, whether in combination with low-or high-intensity climate change scenarios.
Under A2 scenarios, Himalayan landscape will increase farmlands and more urbanized areas in lieu of forests and grasslands ( Figure   S2), therefore impacting two of the most important land use categories driving otter distribution ( Figure S4). Despite a noticeable scarcity of studies designed to forecast both climate and land use change impacts on species distribution (Titeux et al., 2017), a pattern where land use change effects overtake those by climate change has been observed several times (Sirami et al., 2017). For instance, Eglington The importance of considering both climate and land use in predictive studies on biodiversity has been repeatedly underlined (Sirami et al., 2017;Titeux et al., 2017), though it started to be adequately accounted for only in rather recent times (Newbold, 2018;Peters et al., 2019 Cianfrani et al. (2018) showed an overall increase in suitable area under climate change for L. lutra (+10% to +11%), whereas we predicted a systematical loss for this species, especially under the high-intensity land use change scenarios. That said, we cannot exclude that such inconsistencies may partly depend on the different geographical extent adopted in Cianfrani et al. (2018;i.e., global scale) with respect to our study setup. It is also noteworthy that we focused on the Himalayan region, which encompasses the margins of the altitudinal ranges of the three species. Accordingly, resident otters might present local F I G U R E 6 Estimates and statistical significance of the coefficients in the linear mixed models fitted between otter vulnerability and sensitivity or exposure, in turn (p value legend: *, <0.05; **, <0.01; ***, <0.001) adaptations to mountain environments, thus likely being more susceptible to climate and land use alterations if compared to lowland populations. This pattern seems particularly plausible for L. lutra, whose several subspecies are distributed in discrete spatial enclaves, especially along the Himalaya (Hung & Law, 2016). In fact, a solid body of literature has shown that there can be important geographical intraspecific variation in the sensitivity and response to changing environmental conditions (Nice et al., 2019;Peterson et al., 2019).
Such intraspecific variation in environmental requirements substantially derives from local adaptation or adaptive plasticity along the gradient of environmental conditions (Macdonald et al., 2018).
Consequently, predictions of species response to global change that neglect to account for such intraspecific differences may strongly diverge from those generated from local partitions of species distributions (Martin et al., 2020). This evidence could further explain the difference between our predictions of otter future distribution and those provided by Cianfrani et al. (2018).  (Foden et al., 2019), is not necessarily consistent with niche measures of species intrinsic sensitivity (Dickinson et al., 2014). This tendency was also highlighted for otters within the global scale vulnerability assessment performed by Cianfrani et al. (2018). As for the three species we focused on, such assessment reported L. lutra and L. perspicillata as the least and the most vulnerable species, respectively. Interestingly, although Cianfrani et al. (2018) showed a higher exposure for A. cinereus than for L. perspicillata (in contrast to our results), the latter was consistently more vulnerable by virtue of a far higher sensitivity. This body of evidence substantially emphasizes the importance of assessing vulnerability by including trait-based elements, for example, measures of niche width, into "classical" correlative approaches usually adopted to predict range dynamics (Foden et al., 2019;Pacifici et al., 2015).
As for the main limitations of this study, we did not consider the third component usually included in the CCVA theoretical framework, that is, the adaptive capacity (Foden et al., 2019). Such ability to adjust to environmental change should be inferred by traits as for example, phenotypic plasticity or dispersal abilities, which are often unavailable for many species. This lack of information makes adaptive capacity the most challenging component to be quantified in CCVAs, which often simply neglect it (Cianfrani et al., 2018;Foden et al., 2019). It is worth mentioning that species sensitivity and habitat exposure can be given different weights as to tune their relative contribution to the final vulnerability score (Rinnan & Lawler, 2019).
Although the way such weights are defined should bias, at least in theory, the vulnerability score toward sensitivity or exposure, Cianfrani et al. (2018) designed a specific analysis showing that different weighting schemes on the single vulnerability components altered the final score only marginally (especially for the three species analysed here). Furthermore, we were interested in quantifying the relative contribution of sensitivity and exposure on vulnerability accounting for several variability sources (e.g., global circulation models, global change scenarios and binarization thresholds). Accordingly, we found the most obvious choice was to assign the same weight (i.e., 1) to both vulnerability components, instead of artificially introducing any weighting scheme that lacked of a strong ecological rationale.
Lastly, we did not include any projection of climatic extreme events (e.g., drought or floods) in our vulnerability assessment. Despite a well-established evidence that such events will likely have stronger effects on ecological systems than simple shifts in mean values (Thompson et al., 2013), modelling their frequency and occurrence is particularly difficult due to their rarity and stochastic nature. As a result, such phenomena have usually been neglected in conservation studies until very recently (Palmer et al., 2017).

| CON CLUS IONS
Notwithstanding these drawbacks, our study underlines how coupling climate and land use change in CCVAs can generate profoundly diverging results compared to approaches relying on solely climate change. Such ever-mounting evidence appears particularly relevant both to reassess more cautionarily the species ranked as insensitive or favoured by climate change (e.g., L. lutra), and, even more, to increase attention on those species (e.g., African otters A. congicus and A. capensis) that were already classified as highly vulnerable under climate change scenarios only (Cianfrani et al., 2018). The outcomes of this study also highlight the potentialities of performing CCVAs through a "combined" framework, thus emphasizing the importance of including trait-based elements (e.g., species sensitivity) into purely correlative approaches (Foden et al., 2019;Pacifici et al., 2015). In fact, such intrinsic elements proved significantly more important in determining species vulnerability than extrinsic metrics such as habitat exposure. Our analyses show that such a pattern still holds when accounting for several variability sources usually linked to correlative models (e.g., alternative global circulation models and binarization schemes). Future research efforts should focus on developing tools to quantify adaptive capacity in a more reliable and generalizable way than currently available approaches (Williams et al., 2008).
Including this "great absent" in CCVAs, along with a deeper focus on modelling climate variability and extreme events, will surely bring this conservation tool towards a more advanced level.

ACK N OWLED G EM ENTS
We are indebted to the Department of Biosciences and Territory at the University of Molise for their support during the study phase. We are grateful for support from the IUCN/SSC Himalayan Otter Network. We are also grateful to the Jammu and Kashmir Department of Wildlife Protection for granting permission and providing the required logistical help and cooperation for this research in the Ladakh area. Lastly, we would like to express our gratitude to two anonymous reviewers for their time and expertise in improving the quality of this manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad at https://doi.org/10.5061/dryad.zw3r2 287v.