Elevational Shifts of Freshwater Communities Cannot Catch up Climate Warming in the Himalaya

Climate warming threatens biodiversity at global, regional and local levels by causing irreversible changes to species populations and biological communities. The Himalayan region is highly vulnerable to climate warming. This calls for efficient environmental management strategies because biodiversity monitoring is costly, particularly for the developing countries of the Himalaya. Species distribution modeling (SDM) represents a tool that can be used to identify vulnerable areas where biodiversity monitoring and conservation are required most urgently and can be prioritized. Here, we investigated the potential present-day community compositions of river invertebrates in the central and eastern Himalayas and predicted changes in community compositions in future decades using SDMs. We then quantified the climate-induced range shifts of benthic invertebrates along the elevational gradient and tested whether the predicted community shift is fast enough to fully compensate for the projected climate warming. Our model predicts future increases in benthic invertebrate taxonomic richness. Further, projected community shifts are characterized by the movement of warm-dwellers to higher elevations and losses in cold-dwellers. The predicted model shows that benthic invertebrate communities would not be able to compensate climate warming through uphill migration and thus would accumulate climatic debts. Our findings suggest that the ongoing warming effect would cause continued elevational range shifts of mountain river communities.


Introduction
As a response to climate warming, many species have been observed or are expected to shift their distribution ranges to higher latitudes or elevations [1][2][3][4].However, there are considerable variations in how species respond to climate warming.Species adapted to or preferring relatively high temperatures (warm-dwellers) are often anticipated to expand their distribution ranges as warm habitats become increasingly available under a projected warming climate.At the same time, species Water 2016, 8, 327 2 of 12 adapted to or preferring relatively low temperatures (cold-dwellers) are expected to experience range losses.This dynamic could lead to new species interactions and novel community compositions [5][6][7].
There is a need to understand the interactions between the speed with which species and communities can shift ranges, the speed of climate warming in a community's habitat, and the related potential losses of biodiversity.For example, species and communities inhabiting areas that are expecting more rapid climate warming need either to shift their ranges more quickly or are more likely to be stressed or extirpated by failing to track their preferred climatic niche.The speed with which a species or community can shift its range to track a preferred climatic niche is thus decisive for the species' and ultimately a community's potential ability to respond to climate warming, yet it remains largely unknown.We applied the recently developed community temperature index (CTI) [8,9] to assess precisely this situation for benthic invertebrates in Himalayan rivers.CTI accounts for species-specific sensitivity to temperature changes.For example, while a changing community may support a constant number of species, its CTI may change due to a proportional shift of cold-and warm-dwellers.This climatic sensitivity makes CTI a particularly useful method to detect climate-related community changes [10].If a given species shift its range more slowly than the mean temperature (i.e., the isotherm) shifts, the term of "climatic debt" is used to represent how much community change lags behind climate warming [8].Due to the expected impacts of climate warming on community composition, the CTI and assessment of climatic debt represent biologically more meaningful measures of community change than diversity measures alone.
We conducted our study in the Himalayas because (a) compared to lowlands, mountains exhibit steep temperature gradients across a small geographic range requiring that species only need to migrate short distances to track their climatic preferences [11], and we can expect that at the scale of our study, distance alone will not limit the movement of benthic invertebrates; (b) this region was already identified as a suitable place for tracing climate-induced elevational shifts in several species [12,13]; (c) the Himalayas support a high proportion of global biodiversity, including many endemic and endangered species [14]; and (d) species in the Himalaya are known to be particularly vulnerable to warming effects [15,16] because the progressive warming rate at high elevations of the Himalaya is about three times faster than the global average [17].As the water tower of Asia, the Himalayas source ten large Asian rivers that sustain the livelihoods of about 1.3 billion people [17].Climate warming in this region is therefore thought to be a matter of global concern.In the Himalaya, climate-induced responses of terrestrial plants, insects and mammals have been individually surveyed [15,[17][18][19][20]].However, information on how climate warming will affect freshwater benthic invertebrate communities is largely lacking, despite of their being prevalent in these mountains' rivers and inherently sensitive to environmental changes around the globe [21][22][23].
Historical occurrence records have been widely used to investigate climate warming-induced species' range shifts, but such biological data is lacking in most parts of Himalaya.Species distribution models (SDM) can be used as an approximation of species potential habitat for such a lack of data and have been increasingly applied for projecting likely future changes in species' geographical distributions in various taxonomic groups [20, [23][24][25][26][27].We selected pristine or near-pristine catchments in our study because any physical constraints (e.g., man-made dams and habitat fragmentations) on elevational species' shifts are negligible.
In contrast to previous studies in the Himalaya where several individual species were modeled to gain a better understanding on habitat distribution under different climate warming scenarios [20, 28,29], our study estimates the potential distributions of benthic invertebrates with respect to future climate warming at the community level for the first time.Climate warming associated elevational shifts of freshwater species have not been measured or observed and rarely modelled in the Himalayas (but see [20]).Here we compare the projected shifts of both benthic invertebrate communities and isotherms in two distinct Himalayan regions (i.e., central and eastern Himalayas).The regions were selected as we expect the responses of communities to climate warming will vary among Himalayan regions in terms of pronounced differences in climate, topography and biotic composition, as well as political rule, culture, infrastructure development, deforestation, and agricultural activities [16].To address our questions, we performed species distribution modelling of benthic communities to project present-day and future distributions in central and eastern regions in the Himalayas.We investigate if and how taxonomic richness and community compositions are projected to change in response to future climate warming.We specifically assess if (1) benthic invertebrate ranges are projected to shift uphill in response to the expected climate warming; (2) the communities' uphill shifts coincide with a relative increase in warm-dwellers and a decrease in cold-dwellers; (3) communities shift their ranges uphill fast enough to track temperature preference under climate warming; and (4) the differences in landscape and biota between eastern and central Himalayan rivers affect the projected community responses.

Study Region and Data Collection
We studied two catchments in the central and eastern Himalayas, the Langtang and Gangqu catchments, respectively.The Langtang catchment is located in Langtang National Park in north-central Nepal, which drains the high, gentle Langtang valley [22].The Gangqu catchment is located in Shangrila Natural Reserve in Yunnan Province, Southwest China [30].Both catchments exhibit limited human disturbance, and all of the sampling sites are pristine or near-pristine [22,30].These two catchments are located within a similar latitudinal zone (central: 27.83-28.29˝N; eastern: 27.88-28.61˝N; Figure 1), and they both have alpine climatic characteristics.
Water 2016, 8, 327 3 of 12 projected to change in response to future climate warming.We specifically assess if (1) benthic invertebrate ranges are projected to shift uphill in response to the expected climate warming; (2) the communities' uphill shifts coincide with a relative increase in warm-dwellers and a decrease in colddwellers; (3) communities shift their ranges uphill fast enough to track temperature preference under climate warming; and (4) the differences in landscape and biota between eastern and central Himalayan rivers affect the projected community responses.

Study Region and Data Collection
We studied two catchments in the central and eastern Himalayas, the Langtang and Gangqu catchments, respectively.The Langtang catchment is located in Langtang National Park in northcentral Nepal, which drains the high, gentle Langtang valley [22].The Gangqu catchment is located in Shangrila Natural Reserve in Yunnan Province, Southwest China [30].Both catchments exhibit limited human disturbance, and all of the sampling sites are pristine or near-pristine [22,30].These two catchments are located within a similar latitudinal zone (central: 27.83-28.29°N; eastern: 27.88-28.61°N; Figure 1), and they both have alpine climatic characteristics.The Langtang catchment was surveyed between April and June in 2012 and 2013, while the Gangqu catchment was surveyed in May 2005.A multi-habitat sampling approach was used to collect benthic samples in 34 rivers in Langtang catchment.A multi-habitat sample comprises ten habitat specific kick samples over a 100-m river stretch collected with a kick net (0.25 × 0.25 m 2 frame size and a 500 μm mesh size), resulting in 0.63 m 2 of river bed being sampled.In Gangqu catchment, a kick net (1 × 1 m 2 frame size and a 500 μm mesh size) was employed to collect samples from 32 rivers.Here the three dominant microhabitats within a 100 m river stretch were sampled for five minutes each, resulting in 3 m 2 of river bed being sampled.The organisms were sorted and identified to the The Langtang catchment was surveyed between April and June in 2012 and 2013, while the Gangqu catchment was surveyed in May 2005.A multi-habitat sampling approach was used to collect benthic samples in 34 rivers in Langtang catchment.A multi-habitat sample comprises ten habitat specific kick samples over a 100-m river stretch collected with a kick net (0.25 ˆ0.25 m 2 frame size and a 500 µm mesh size), resulting in 0.63 m 2 of river bed being sampled.In Gangqu catchment, a kick net (1 ˆ1 m 2 frame size and a 500 µm mesh size) was employed to collect samples from 32 rivers.Here the three dominant microhabitats within a 100 m river stretch were sampled for five minutes each, resulting in 3 m 2 of river bed being sampled.The organisms were sorted and identified to the genus or species level, except Chironomidae and Oligochaeta, which were identified to the sub-family or family level.Prior to analysis, taxa occurring in less than 10% of the samples (ď3 occurrences) were eliminated to reduce stochastic occurrences and to maintain high performances of SDMs.This elimination resulted in 52 and 60 taxa being selected in Langtang and Gangqu catchments, respectively (Supplementary Table S1).
For each site, we extracted three topographical variables, namely elevation, slope, and aspect (i.e., the compass direction that a slope faces) from a digital elevation model (http://reverb.echo.nasa.gov;resolution: 30 m) and obtained 19 bioclimatic variables from the WorldClim dataset (http://www.worldclim.org;resolution: 1 km).The 19 bioclimatic variables for present-day (average of  and future decades (2030s, 2050s, and 2080s) were derived from the Canadian Centre for Climate Modelling and Analysis (CCCMA) global circulation model with the emission scenarios of the representative concentration pathway (RCP) 4.5 (the possible range of radiative forcing value in the year 2100 relative to the pre-industrial value is +4.5 W¨m ´2).We selected RCP 4.5 because it includes a balanced mix for all energy sources and closely matches the most likely realistic situation for national energy sources [31].Prior to SDM calculations, we excluded highly correlated variables (|r| > 0.75 in the Pearson correlation analysis) to reduce the redundancy of environmental variables.We retained seven abiotic variables (elevation, slope, aspect; bio06 = the minimum temperature of the coldest month; bio11 = the mean temperature of the coldest quarter; bio15 = precipitation seasonality; and bio19 = precipitation of coldest quarter) as the final abiotic variables for SDMs (Supplementary Table S2).

Species Distribution Model
We performed SDMs with the present-day geographical records of each taxon and the selected seven abiotic variables in four time periods (present, i.e., 2000s; 2030s; 2050s; and 2080s) to predict species' potential distributions in the present and future decades.We used the BIOMOD2 package [32] in R [33] with eight algorithms in the SDM: artificial neural network (ANN), classification tree analysis (CTA), flexible discriminant analysis (FDA), general boosting model (GBM), generalized linear model (GLM), maximum entropy (MAXENT), multiple adaptive regression splines (MARS), and random forest (RF).For each algorithm, three fold cross-validation and 1000 pseudo-absences were used, and the prevalence value was set as 0.5.
Prior to analysis, the occurrence of each taxon was randomly split into two sets, a training set (70%) and a testing set (30%).The training set was used to calibrate the model, and the testing set was used to validate the calibrated model.Models were evaluated with the true skill statistic (TSS), which has been proven to be superior to quantify the performance of SDM with presence prediction data [34].TSS is an integrity measure incorporating sensitivity (true positive prediction, namely presence) and specificity (true negative prediction, namely absence).TSS ranges from 0 to 1, where 0 represents no discrimination and 1 represents perfect discrimination.We retained only models with TSS ě 0.6 for the ensemble model to produce more robust predictions [35].The ensemble model was generated based on the weighted averages of the predictive performance of the algorithm.The relative importance of each algorithm was calculated by multiplying the averaged TSS value by a weight decay of 1.6 [36].Finally, the distribution probability maps were converted into binary maps by applying a cut-off value which minimizes the difference between sensitivity and specificity [36].The SDM projection was restricted within river networks since only water courses were considered as potential habitats for benthic invertebrates.To account for variability and uncertainty in predicted taxonomic richness, 500 random sites (elevation > 1000 m) along the river network were created for each catchment, and all analyses with predicted values were performed based on these 500 random sites.Random sites Water 2016, 8, 327 5 of 12 were selected using the package sp [37] in R [33].We assessed the temporal trends of overall taxonomic richness and richness in three elevational bands (low (ď3000 m), medium (3000-4000 m), and high elevations (>4000 m); Table 1) for 2000s-2080s to assess if benthic invertebrates will migrate uphill in response to future climate warming.We evaluated the differences in taxonomic richness among decades and elevational bands with analysis of variance (ANOVA), and Tukey's multiple comparison test was used to detect significant differences (α = 0.05).
Table 1.Summary of the physical characteristics, the regression of taxonomic richness against period, and the analysis of variance (ANOVA) of taxonomic richness among periods across the entire catchment and three elevational bands in the central and eastern Himalayan regions.The slope and intercept are obtained from linear regressions of taxonomic richness against periods (2000s, 2030s, 2050s, and 2080s are expressed as 0, 3, 5, and 8 in linear regressions, respectively).We further assessed if community shifts to higher elevation coincide with relative increases in warm-dwellers and decreases in cold-dwellers by examining the temporal trends in the relative richness of cold-and warm-dwellers during the 2000s and 2080s.Cold-and warm-dwellers were defined as taxa belonging to the 10th and 90th percentiles of species temperature preference index (STPI) of all the observed taxa.STPI of a given species was calculated as the average present-day air temperature conditions (bio01 in the WorldClim dataset) over the 500 random selected grids from the modelled species distribution.Air temperature was used as a surrogate for river temperature because:

Region
(1) species preference in river temperature is not available for the study areas; and (2) air and river water temperature tend to be linearly correlated in many study regions [38].STPI may not reflect the absolute temperature preference, but a relative temperature preference can be derived for each taxon using this approach if physiological or experimental evidence is lacking.
CTI is a measure of average thermal condition for a given ecological community, with a low value indicating a large proportion of cold-dwellers and a high value representing more warm-dwellers.We compared temporal trends in CTI from the 2000s to 2080s to assess if the communities can shift their range fast enough to keep up with climate warming or alternatively develop a climatic debt.The CTI of a given site was calculated as the average of STPI weighted by species composition (presence = 1 and absence = 0).The rate of change in the community's thermal preference (i.e., the observed rate of change) was defined as the slope of CTI over the four distinct time periods using a linear regression.The rate of climate warming was defined as the slope of annual average air temperature (AAAT; i.e., bio01 in the WorldClim dataset) over the four distinct time points using a linear regression, and it represents the required rate of change in the thermal preference that a community would need to undergo to avoid a climatic debt.We assessed the significant differences between the observed and required rates of change using a Wilcoxon test with 9999 random bootstraps in R [33].These analyses were repeated ten times, and the final p value was the median of these repeats.A climatic debt is predicted if the difference between the observed and required temporal speed is significant.
The results of all of the modeled taxa were thus used in the following analyses.
The predicted taxonomic richness increased over time in both central (ANOVA test: F 3, 1996 = 221, p < 0.0001; Figure 2a) and eastern regions (F 3, 1996 = 36, p < 0.0001; Figure 2b).However, the increase occurred quicker in the central (1.59 taxon¨decade ´1) than in the eastern region (0.91 taxon¨decade ´1).After splitting the entire 500 random sites into three elevational bands, this difference was more pronounced (Figure 2a,b).Interestingly, the speed with which taxonomic richness increased was much slower at low elevations (0.85 taxon¨decade ´1) than at high elevations (2.96 taxon¨decade ´1; Table 1) in the central region, whereas the trend was opposite in the eastern region (low elevation: 1.79 taxon¨decade ´1, high elevation: 0.63 taxon¨decade ´1; Table 1).
Water 2016, 8, 327 6 of 12 After splitting the entire 500 random sites into three elevational bands, this difference was more pronounced (Figure 2a,b).Interestingly, the speed with which taxonomic richness increased was much slower at low elevations (0.85 taxon•decade −1 ) than at high elevations (2.96 taxon•decade −1 ; Table 1) in the central region, whereas the trend was opposite in the eastern region (low elevation: 1.79 taxon•decade −1 , high elevation: 0.63 taxon•decade −1 ; Table 1).
(a) (b) In the central region, both cold-and warm-dwellers were projected to significantly increase over time (0.07 and 0.32 taxon•decade −1 , respectively), with warm-dwellers exhibiting a quicker increase (Figure 3a; Supplementary Table S3).In the eastern region, cold-dwellers were projected to significantly decrease (−0.05 taxon•decade −1 ) and warm-dwellers were projected to increase (0.02 taxon•decade −1 ) over time (Figure 3b, Supplementary Table S3).In the central region, both cold-and warm-dwellers were projected to significantly increase over time (0.07 and 0.32 taxon¨decade ´1, respectively), with warm-dwellers exhibiting a quicker increase (Figure 3a; Supplementary Table S3).In the eastern region, cold-dwellers were projected to significantly decrease (´0.05 taxon¨decade ´1) and warm-dwellers were projected to increase (0.02 taxon¨decade ´1) over time (Figure 3b, Supplementary Table S3).
Both CTI and AAAT were projected to increase over time, but the temporal speed of CTI was 0.14 ˝C¨decade ´1, and significantly lower (Wilcoxon test, p < 0.05) than the required temporal speed of isotherm (0.45 ˝C¨decade ´1) in central region, causing a climatic debt of 0.31 ˝C¨decade ´1 (Figure 4a, Supplementary Table S3).An even greater climatic debt (0.40 ˝C¨decade ´1) was projected in eastern region, where the temporal speed of CTI and AAAT were 0.02 and 0.42 ˝C¨decade ´1, respectively (Figure 4b, Supplementary Table S3).
In the central region, both cold-and warm-dwellers were projected to significantly increase over time (0.07 and 0.32 taxon•decade −1 , respectively), with warm-dwellers exhibiting a quicker increase (Figure 3a; Supplementary Table S3).In the eastern region, cold-dwellers were projected to significantly decrease (−0.05 taxon•decade −1 ) and warm-dwellers were projected to increase (0.02 taxon•decade −1 ) over time (Figure 3b, Supplementary Table S3).S1; statistical results of linear regressions are available in Supplementary Table S3.
Both CTI and AAAT were projected to increase over time, but the temporal speed of CTI was 0.14 °C•decade −1 , and significantly lower (Wilcoxon test, p < 0.05) than the required temporal speed of isotherm (0.45 °C•decade −1 ) in central region, causing a climatic debt of 0.31 °C•decade −1 (Figure 4a,  S1; statistical results of linear regressions are available in Supplementary Table S3.Supplementary Table S3).An even greater climatic debt (0.40 °C•decade −1 ) was projected in eastern region, where the temporal speed of CTI and AAAT were 0.02 and 0.42 °C•decade −1 , respectively (Figure 4b, Supplementary Table S3).

Discussion
Unlike many studies that have projected future shifts of a few surveyed species for the Himalayas, we quantitatively assessed range shifts of the whole community to draw a more complete picture of the impacts of climate warming on river biodiversity.Our study represents a first approximation about how changes in climatic condition alter river invertebrate communities on the top of the world.River invertebrate communities are projected to shift their ranges uphill in response to future climate warming.Our analyses further indicate that taxonomic richness will increase in both central and eastern Himalayan regions, and not decrease as often warned [17].In our study these changes in taxonomic richness are based on an increase in warm-dwellers and less severe and ambiguous changes in cold-dweller's richness.A clear climatic debt is then projected for communities in both regions.With a greater climatic debt, more or less stable richness in warm-dwellers and losses in cold-dwellers, we therefore expect that the eastern region will face higher climate-induced stress compared to the central regions of the Himalaya.This reflects a general vulnerability of mountain communities to future climate warming, but accentuates the importance of regional variation in these threats.
Our study reveals a clear uphill shift of benthic invertebrates in response to future climate warming.This is in line with many studies e.g., on plants in the Swiss Alps [39], butterflies in central  S3.

Discussion
Unlike many studies that have projected future shifts of a few surveyed species for the Himalayas, we quantitatively assessed range shifts of the whole community to draw a more complete picture of the impacts of climate warming on river biodiversity.Our study represents a first approximation about how changes in climatic condition alter river invertebrate communities on the top of the world.River invertebrate communities are projected to shift their ranges uphill in response to future climate warming.Our analyses further indicate that taxonomic richness will increase in both central and eastern Himalayan regions, and not decrease as often warned [17].In our study these changes in taxonomic richness are based on an increase in warm-dwellers and less severe and ambiguous changes in cold-dweller's richness.A clear climatic debt is then projected for communities in both regions.With a greater climatic debt, more or less stable richness in warm-dwellers and losses in cold-dwellers, we therefore expect that the eastern region will face higher climate-induced stress compared to the central regions of the Himalaya.This reflects a general vulnerability of mountain communities to future climate warming, but accentuates the importance of regional variation in these threats.
Our study reveals a clear uphill shift of benthic invertebrates in response to future climate warming.This is in line with many studies e.g., on plants in the Swiss Alps [39], butterflies in central Spain [2], insects in Bavarian Forest in Germany [40], and a combined taxonomic groups at a global scale [4], and river invertebrates in central Europe [23,41,42].While these former studies may be confounded by heavy anthropogenic impact on river ecosystems, the effect of human disturbance on the signal of an uphill shift in the Himalayan streams is negligible.The lack of longitudinal barriers in our study rivers is particularly beneficial to river species' migrating uphill through the widely connected river channels.River species may use this dendritic network structure as a highway to higher elevations in response to climate warming [43].However, our study also indicates that the speed of uphill range shifts may vary at different elevations and among regions within a mountain system.In case, a quicker uphill movement was projected at higher elevations in the central region, suggesting that these areas currently harboring a high taxonomic diversity [39,44] are likely to face a mountain summit trap effect in the future [20,41].Unexpectedly our results suggest a slower speed of increasing taxonomic richness at high elevations in eastern region.This may be the result of the summit trap effect, where species' ability to track climate conditions uphill is necessarily stopped at the summit of mountain, leading to abrupt losses in biodiversity [41].Nevertheless, the extinction of species and summit trap effect are interconnected, and the potential summit trap effect is inversely related to the height of mountain [6,45].The effect should be less evident in the higher central Himalayan mountains and more evident in the relatively lower elevation of Gangqu [16].Beside the effect of landscapes, the difference in biotic composition may be another important driver underlying the pattern of the slowly increasing taxonomic richness at the high elevations in the eastern region.Under an expected climate warming, many warm-dwellers may expand their range to higher elevations thereby compensating for the loss of taxonomic richness caused by the extirpation of cold-dwellers.This will likely not happen in areas with only few warm-dweller species, e.g., the eastern region, where only one Hirudinea and two Oligochaeta prefer warm conditions.
Devictor et al. [9] reported that the latitudinal shifts in French bird communities were driven by an increase in warm-dwellers from 1989 to 2006.In our case, the losses or very limited gains of cold-dwellers and largely gains of warm-dwellers in the eastern region are in line with previous studies [5,9].In the high mountains of the central Himalayas, the colonization was dominated by warm-dwellers from low elevations and cold-dwellers from medium elevations.Nevertheless, it is worthy to note that the continued loss of previous co-occurrences and the gain of novel co-occurrences may extensively alter community compositions.
A significant increase in CTI of benthic invertebrates was detected in both regions, which is consistent with similar studies in plants, birds and butterflies in Switzerland [10] as well as birds and butterflies at an European continental scale [8].Based on the rule of thumb that air temperature decreases by 0.6 ˝C at every 100 m of higher elevation, the speed of CTI changes are equivalent to 23.3 and 3.3 m¨decade ´1 in central and eastern regions, respectively.In many other systems the communities' range shifts are lagging behind isothermal shifts [4,14].This is also true in our study where the isotherm is projected to shift uphill by 75.0 and 70.0 m¨decade ´1 in the central and eastern region, consequently causing climatic debts of 51.7 and 66.7 m¨decade ´1 in the central and eastern region, respectively.These climatic debts are similar to those observed in butterflies (62.4 m¨decade ´1) and birds (56.7 m¨decade ´1) in the Alps [10], but smaller than those reported for birds (117.8 km¨decade ´1 along the latitudinal gradient, equivalent to 117.8 m¨decade ´1 along the elevational gradient) and butterflies (equivalent to 75.0 m¨decade ´1) by Devictor et al. [8].Actual climatic debts, i.e., not those inferred from SDMs, could arise when the low dispersal ability of species does not allow species to migrate with their suitable climate space [46], a change in species interactions [14] or a combination of both.Furthermore a perceived climatic debt might simply result from phenotypic plasticity that allows a species to persist under a different climate without a need to migrate [47].In addition, extreme climatic events (i.e., flash floods, drought, and extreme hot or cold conditions) may also affect aquatic species' climate-induced range shifts [48].
Similar to most studies forecasting climate-induced range shifts, some limitations are inherent in our approach.First, the sampling years were different between central (2012 and 2013) and eastern catchments (2005).However, the difference in sampling years is comparably smaller than the projection periods (i.e., 2030s, 2050s, and 2080s), and we are convinced that such a difference has limited effects on our results.Second, we conducted SDM based on existing field surveys, but the full distribution areas of species may not be captured, especially given the very low number of occurrences available for several species.Third, due to a lack of taxonomical studies in these two regions, most taxa were identified to the genus level, which may mask species-specific differences in responses to climate warming [22].Fourth, rare species were excluded from the analysis to maintain a high model performance, which is subject to the application restrictions of SDM [25].Yet these rare species may be ecologically particularly sensitive and could provide more detailed insights into species responses to climate warming, although their effect on the community are likely minimal due to their limited occurrence.Fifth, beyond climatic and topographical variables, other potentially important variables were not included in our study, such as substrate and water quality metrics.Unfortunately, these data do not exist at the scale of our study.Furthermore, biotic interactions and species traits were not modelled to avoid over-parameterization of the models.Sixth, the changes in hydrological conditions that species must deal with when moving uphill through a river network are not considered in our study.Finally, the spatial resolution (1 km 2 ) used in our study may fail to extract climatic variations at the fine spatial scale most relevant to the investigated species.It should be noted that no modelling studies can overcome such limits entirely, and that directly monitoring the impact of climate warming on ecological communities is costly and requires years of data.Yet greater efforts of field surveys would improve the strength of the model and inferences.Thus, as more community data become available, follow-up studies may consider incorporating more of the above factors at finer spatial resolutions.

Conclusions
Our findings have global conservation implications, as considerable ecological changes of river communities are expected by the 2080s in the high mountain regions of the Himalayas, where most of the recognized biodiversity is endemic.Importantly, our findings highlight that ongoing climate warming will likely lead to severe upward shifts in benthic communities.An uphill shift in warm-dwelling species may be superficially beneficial for mountain communities, but the summit trap effect is likely to threaten high mountain biodiversity in the long run.Indeed, the projected opposite trends in cold-and warm-dwellers as well as large climatic debts will likely form novel communities with taxa occurring outside of their present geographical ranges, and this could affect river ecosystem functioning and present a major challenge for conservation planning.We propose that extensive datasets are vital for interpreting biodiversity changes, and that establishing robust monitoring systems for various ecological communities in the Himalaya will improve our ability to understand and manage climate-induced community-level changes [22].

Figure 1 .
Figure 1.Geographical locations of the study catchments and sampling sites in central and eastern part of Hindu-Kush Himalayan (HKH) region.

Figure 1 .
Figure 1.Geographical locations of the study catchments and sampling sites in central and eastern part of Hindu-Kush Himalayan (HKH) region.

Figure 2 .
Figure 2. Predicted taxonomic richness of benthic invertebrates in the (a) central; and (b) eastern Himalayan regions.Dots represent mean, whiskers are standard error, and the letter indicates significant differences among periods within each elevational band (Tukey's multiple comparison test, p < 0.05).

Figure 3 .
Figure 3. Violin plots of temporal trends in cold-and warm-dwellers over time in the (a) central; and (b) eastern Himalayan regions.Dots are mean, thick black bars within the violin are interquartile

Figure 2 .
Figure 2. Predicted taxonomic richness of benthic invertebrates in the (a) central; and (b) eastern Himalayan regions.Dots represent mean, whiskers are standard error, and the letter indicates significant differences among periods within each elevational band (Tukey's multiple comparison test, p < 0.05).

Figure 3 .
Figure 3. Violin plots of temporal trends in cold-and warm-dwellers over time in the (a) central; and (b) eastern Himalayan regions.Dots are mean, thick black bars within the violin are interquartile range, and shaded areas denote 95% confidence intervals of linear regressions.The taxonomic lists of cold-/warm-dwellers are available in Supplementary TableS1; statistical results of linear regressions are available in Supplementary TableS3.

Figure 3 .
Figure 3. Violin plots of temporal trends in cold-and warm-dwellers over time in the (a) central; and (b) eastern Himalayan regions.Dots are mean, thick black bars within the violin are interquartile range, and shaded areas denote 95% confidence intervals of linear regressions.The taxonomic lists of cold-/warm-dwellers are available in Supplementary TableS1; statistical results of linear regressions are available in Supplementary TableS3.

Figure 4 .
Figure 4. Violin plots of temporal trends in community temperature index and annual average air temperature in the (a) central; and (b) eastern Himalayan regions.Dots are mean and shaded areas denote 95% confidence intervals.Statistical results of linear regressions are available in Supplementary TableS3.

Figure 4 .
Figure 4. Violin plots of temporal trends in community temperature index and annual average air temperature in the (a) central; and (b) eastern Himalayan regions.Dots are mean and shaded areas denote 95% confidence intervals.Statistical results of linear regressions are available in Supplementary TableS3.