Predicting range shifts of pikas (Mammalia, Ochotonidae) in China under scenarios incorporating land use change, climate change and dispersal limitations

Two of the most important forces affecting biodiversity are land use change (LUC) and global climate change (GCC). Previous studies have modelled their impacts on species separately, and together, but few have done so for multiple species with dispersal limitations incorporated into the models.


| INTRODUC TI ON
Two of the most important contemporary forces affecting biodiversity are land use change (LUC) and global climate change (GCC) (Newbold, 2018;Powers & Jetz, 2019). Each, acting on its own, can harm or benefit individual species. For example, LUC can destroy or create habitats, depending on the ecology of the species in question (Gossner et al., 2016;Lawler et al., 2014). GCC can do the same by forcing species to retreat from current occupied habitats or, alternatively, by creating new areas of climate suitability that species can colonize (Román-Palacios & Wiens, 2020;Schwalm et al., 2016).
Moreover, the combination of LUC and GCC can dampen or amplify overall impacts on biodiversity (Northrup et al., 2019;Radinger et al., 2016), which makes predicting the responses of species to these forces very challenging.
An additional key factor in predicting the responses of species to these changes is species' dispersal ability because both habitat fragmentation and reduced climate connectivity can impede the ability of species to respond to environmental change (Gouveia et al., 2016;Senior et al., 2019). A species' dispersal range and dispersal frequency, as well as the presence of natural barriers (e.g., lakes or rivers for small mammals), all affect whether individuals of that species can reach suitable yet unoccupied areas and how long it will take them to do so (Franklin, 2010). Simply stated, new expanses of habitat potentially created by GCC or LUC will have little benefit for a species if the species cannot colonize them (Holloway et al., 2016). Previous studies have modelled the likely impacts of LUC and GCC on species separately (Gossner et al., 2016;Román-Palacios & Wiens, 2020) and collectively (Northrup et al., 2019;Radinger et al., 2016), but very few have attempted to predict the impacts of both LUC and GCC on multiple species with dispersal limitations incorporated into the models (but see Troia et al., 2019).
For a given species, GCC is likely to degrade or destroy habitat suitability at lower elevations but to create new suitable habitats at higher elevations, as GCC shifts suitable thermal conditions upslope (Freeman et al., 2018). Thus, high-elevation species may be especially vulnerable to GCC because it can cause the loss of a large proportion of their suitable habitat (Colwell et al., 2008). On the other hand, people typically settle and modify lower-elevation areas first and then expand upslope (Elsen et al., 2020;Nogués-Bravo et al., 2008).
Thus, low-elevation species are more likely to be affected by LUC than are their high-elevation counterparts (Freeman et al., 2018;Galbreath et al., 2009). For the same reason, one would predict that LUC is most likely to degrade or destroy a species' habitat at the lower elevations within its distributional range (although local conditions may upend this prediction) (Guo et al., 2018). Finally, the topography of mountain ranges will constrain how species respond to both LUC and GCC because the land area present at different elevations varies dramatically from one range to another (Elsen & Tingley, 2015). This constraint is likely to be especially important for species restricted to intact landscapes (Elsen et al., 2020).
Pikas (Mammalia, Ochotonidae, genus Ochotona) living in the Qinghai-Tibet plateau (QTP) and adjacent areas of China offer an especially rich system for exploring the interactions of LUC, GCC and dispersal along elevational gradients for several reasons: (a) Previous research has established that pikas are quite sensitive to GCC (Calkins et al., 2012;Wang et al., 2006); (b) the dispersal behaviour of several species has been reasonably well studied (Peacock, 1997;Smith, 1974) and, given the overall similarity of all members of the genus to each other, it is not unrealistic to assume that the other species behave similarly; (c) Pikas are important ecologically both as prey species for numerous predators and as burrowers that create habitats for other species (Smith et al., 2019). In some parts of the QTP, they are also important economically as abundant herbivores that may compete with livestock for forage (Retzer & Reudenbach, 2005;Smith et al., 2006); (d) the QTP lies at the junction of four biodiversity hotspots (Myers et al., 2000; see also an updated map developed by Conservation International [https://www. conse rvati on.org/]) and may therefore become an upslope refugium for range-restricted species; (e) despite the high elevation of the QTP, the human population there has increased from 2.3 million in the 1950s to 9.5 million now (Qi et al., 2020), with LUC an ongoing factor affecting the pikas (Chen et al., 2014); moreover, (f) the QTP is expected to experience significant changes in temperature as a result of GCC (Chen et al., 2014).
Here, we integrate species distribution models (SDMs) and a dispersal model to predict the future impacts of LUC and GCC on the ranges of five species of pikas in the QTP region ( Figure 1) (Chen et al., 2014;Li et al., 2017). The annual mean temperature of the QTP rose by 0.16°C per decade from 1955 to 1996, a rate that exceeds the warming rate found elsewhere at similar latitudes (Liu & Chen, 2000). We examine the interactive effects of LUC and GCC on these five pika species, consider how species' range shifts and overall vulnerability vary as a function of the elevations at which they occur, and explore how dispersal behaviour and natural barriers affect their ability to access suitable habitats in the future. In doing so, we provide insights into the conservation of montane species in the face of pervasive environmental change.  Karger et al., 2017) representing both current  and future (2021-2040, 2041-2060, 2061-2080) conditions. For land use variables, we summed five agricultural variables (showing percentages of the five types of agricultural land in each grid) into one "agricultural land" layer and removed the variables for urban area, primary forest and secondary forest, habitats that are not used by the five species and are not predicted to increase greatly in area through 2080. To keep the resolution the same across environmental variables, we resampled land use variable layers using a climate variable layer (bio01) as a reference. To reduce multicollinearity, we used the package "usdm" (Naimi et al., 2014) to select a subset of variables (seven land use variable and five bioclimate variables; Table S1) according to a collinearity threshold measured by the variance-inflation factor (VIF < 2.5). For climate change, we considered three representative concentration pathways (RCP4.5, RCP6.0, RCP8.5) representing different concentration trajectories from intermediate (RCP4.5) to worst case (RCP8.5). Correspondingly, for land use change, we considered three shared socioeconomic pathways (SSP2-4.5, SSP4-6.0, SSP5-8.5). For bioclimate variables representing future conditions, we considered five general circulation models (GCMs): BCC-CSM1.1, CSIRO-Mk3-6-0, FIO-ESM, GFDL-ESM2M, MRI-CGCM3, because these GCMs have proved to perform well in predicting past climates in the QTP region (Jia, Ruan, Yang, & You, 2019;Jia, Ruan, Yang, & Zhang, 2019;Su et al., 2013;Wei et al., 2019). To account for uncertainty in the climate change projections, we present our results based on the averaged GCM outputs. For land use variables representing future conditions, we used data from the recommended marker integrated assessment models (IAMs: MESSAGE for SSP2, GCAM for SSP4, MAGPIE for SSP5), which were included in the LUH2 data set (Hurtt et al., 2016). List and in LUH2 were reclassified as forest versus non-forest for matching purposes; Table S2). To avoid potential spatial autocorrelation among occurrence records or uneven sampling, we thinned the records using the package spThin, with a thinning distance of 10 km (Aiello-Lammens et al., 2015). See Table S3 for number of occurrence records for each species used in our analyses and Appendix S1: Supplementary Methods for detailed descriptions of environmental and species data.

| Species distribution models
We used the ensemble of small models (ESMs) approach to calibrate SDMs for the five species (Breiner et al., 2015(Breiner et al., , 2018. We used this approach to avoid overfitting in each small model, while keeping all important predictors within the final ESMs. For each species, 66 bivariate models (all possible bivariate combinations among 12 variables) were calibrated, evaluated and averaged to an ensemble weighted by Somer's D. We used the R package "ecospat" (Broennimann et al., 2021) to implement this method and to assemble results from multiple algorithms.
The algorithms we used included GLM (generalized linear model), ANN F I G U R E 1 The workflow of this study. GCC, global climate change; LUC, land use change; SDMs, species distribution models (artificial neural networks) and MAXENT.Phillips (maximum entropy), all of which have been recommended to be used in ESMs (Breiner et al., 2015(Breiner et al., , 2018. Individual algorithms were tuned for each species with the "BIOMOD_tuning" function in the R package "biomod2" (Thuiller et al., 2019); we did so by training (with occurrence records and current environment) and comparing models with all possible combinations of hyperparameters. We employed fivefold cross-validation (10 repeats) of our models by randomly splitting our occurrence records into two subsets: 80% of the data were used to calibrate the models, while the remaining 20% of the data were used for testing. We assembled bivariate models and ESMs (using different algorithms) by estimating the weighted mean of predicted distribution probabilities using Somers' D (rescaled version of AUC) as the weight (Breiner et al., 2015). We evaluated the performance of our ESMs using the metrics of "AUC," "TSS" and "Boyce" (Breiner et al., 2015). To avoid selecting a high proportion of less informative pseudo-absences (Barbet-Massin et al., 2012), we restricted the modelling to a bounding box that covered the current ranges of all five pika species plus an extended area outside that perimeter. This bounding box (61-130E, 17-62N) covered all grid cells that the five pika species could possibly appear in by 2080 (see Appendix S1: Supplementary Methods for more details).
We projected the probabilities of presence of the pika species using the ESMs described above for grids in the bounding box under current and future (2021-2040, 2041-2060, 2061-2080) environmental conditions. Our goal was to tease apart the impacts of LUC and GCC acting individually and together on pika distributions. Thus, we combined three sets of variables representing future conditions to create three scenarios: future land use + current climate (a scenario in which just land use is changing: "LUC" scenario); current land use + future climate (a scenario in which just climate is changing: "GCC" scenario); future land use + future climate (both land use and climate are changing: "BOTH" scenario). We also modelled a "Dispersal Only" scenario, in which the suitable habitats for the future periods (2021-2040, 2041-2060, 2061-2080) are the same as those for the current period, and environmental variables are assumed to be constant through time. This "Dispersal Only" scenario models the continued range expansion of each pika species over time (see "Dispersal model" below) into suitable but unoccupied habitats in a world that is not experiencing LUC or GCC (on the assumption that other, unknown factors do not prevent the pikas from colonizing that suitable habitat). We used the results of the "Dispersal Only" model as a baseline for quantifying predicted impacts of LUC and GCC on the pikas. See the Appendix S1: Supplementary Methods and Appendix S2: ODMAP for a detailed protocol for the SDMs.

| Dispersal model
A pika species' dispersal ability, coupled with the presence (or absence) of natural topographic barriers, will determine how much of the suitable habitat that appears in the future due to LUC and GCC can actually be colonized by that species. Therefore, we needed to create a dispersal model to simulate realistic movements from current occupied habitats into unoccupied but suitable habitats that appear in the future. Because we assumed that even today each species of pika occupies only a portion of the suitable habitat available to it (i.e., suitable habitat within the IUCN polygons), in our modelling we allowed each species to disperse to and colonize some suitable but currently unoccupied habitats in the future under the "Dispersal Only" scenario. Our determination of the distances pikas disperse (see Figure S6 for the dispersal kernel used) and how frequently they disperse, as well as what barriers they can or cannot cross is based on studies of better known pika species and described in greater detail in the Appendix S1: Supplementary Methods. For this step, we used the dispersal model MIGCLIM (R package; Engler et al., 2012).

| Abandoned habitats, colonized habitats and net change in occupied habitats
For each pika species, we identified grid cells that are currently occupied (according to our models) but predicted to become unsuitable in 2080, and classified them as "abandoned habitats." We also identified (a) grid cells that are currently suitable but unoccupied and are predicted to be colonized by 2080 and (b) grid cells that are currently unsuitable and are predicted to become suitable and be colonized by 2080, and we classified both types as "colonized habitats." We considered the interaction between LUC and GCC as synergistic or additive when the combined effect on the amount of abandoned or colonized habitat was larger than or equal to the sum of their separate effects ("BOTH" ≥ "LUC" + "GCC").
Otherwise, we considered the roles of LUC and GCC in creating new habitat or reducing current habitat to show at least partly cancelling effects when the combined effect was smaller than the sum of their separate effects ("BOTH" < "LUC" + "GCC") (Radinger et al., 2016). We also calculated the net change in the amount of occupied habitat for each species by 2080 by subtracting the amount of abandoned habitat from the amount of colonized habitat.

| Elevation of abandoned and colonized habitat
We extracted the elevations of abandoned and colonized habitats from a global digital elevation model product (GMTED2010; downloaded from EarthEnv; resolution = 30 arc seconds; Amatulli et al., 2018). As the assumption of variance homogeneity was not met, we ran Kruskal-Wallis rank-sum tests to compare the elevations of abandoned and colonized habitats under the different environmental change scenarios for each species. When the tests were significant, we performed multiple pairwise comparisons of the means of the scenarios using the Dunn test (Ogle et al., 2020). We calculated Cohen's d as a measure of effect size for each pairwise comparison and assessed the corresponding magnitude using the following thresholds: |d| < 0.2 "negligible," 0.2 < |d| < 0.5 "small," 0.5 < |d| < 0.8 "medium," otherwise "large" (package "effsize"; Makowski et al., 2019). We checked the homogeneity of variances and the normality of residuals with Levene's test (package "car"; Fox & Weisberg, 2019) and the Shapiro-Wilk test, respectively.

| Overall elevational range shifts
We estimated quantiles (10th, 20th, … 90th) of the elevations of occupied habitats under current conditions and in the future under our environmental change scenarios. We then calculated the shifts in elevation of each quantile to determine in which direction, and by how much, the elevational distributions in the different time periods differed. We estimated the quantiles using the Harrell-Davis quantile estimator (Harrell & Davis, 1982), which we implemented using the WRS package (Wilcox & Schönbrodt, 2019). We conducted all data analyses in R (R Core Team, 2019).

| Performance of SDMs and importance of variables
The ESMs performed well for all five pika species (all AUC > 0.943, all TSS > 0.784; all Boyce > 0.642; see Table S4 for performances of models using specific algorithms). Climate variables were more important than land use variables in explaining the distributions of O. cansus (p = .007), while no significant differences between the two types of variables were detected for the other four species (all p > .05; Figure S7; see Figures S1-S5 for response curves for each species).

| Sensitivity to the choice of representative concentration pathways and shared socioeconomic pathways
We present results using an intermediate representative concentration pathway (RCP6.0 for climate change) and the corresponding shared socioeconomic pathway (SSP4-6.0 for land use change) in the main text. Results using two other combinations of RCPs and corresponding SSPs (called "pathway sets" for short) are presented in the supplementary materials (Tables S7-S10; Figures S11-S20); these include RCP4.5+SSP2-4.5 and RCP8.5+SSP5-8.5. When presenting the results below, we have indicated if a given pattern (not exact values) did not hold across all pathway sets.

| Abandoned habitats
The relative amount of abandoned habitat under the "LUC" scenario versus the "GCC" scenario varied between species and pathway sets.

| Colonized habitats
The relative amount of colonized habitat under the "LUC" scenario versus the "GCC" scenario also varied between species and pathway sets. Under the [RCP6.0 + SSP4-6.0] pathway set, the "LUC" scenario resulted in more colonized habitat than the "GCC" scenario for O. dauurica, while the reverse was true for the other four species pathway set). In addition, the amount of habitat colonized under the "BOTH" scenario was consistently less than the sum under the "LUC" plus the "GCC" scenarios for all species and all pathway sets.
Compared to the "GCC" scenario, most of the colonized habitats under the "LUC" scenario occurred at lower elevations (all p < .001; Tables S6, S8, S10; Figure 3a

| Variables' contributions to the loss and gain of suitable habitats
Under the "LUC" scenario, the main land use changes that contributed to the loss or gain of suitable habitats were as follows: changes in secondary mean age (i.e., time period since the land was last logged or used to grow crops), non-forested secondary land and rangeland (Figures S8, S15, and S20). Under the "GCC" scenario, the main changes in climate variables that contributed to losses or gains of suitable habitat for pikas were related to increases in temperatures during specific seasons (bio08 and bio09).

| Net changes in amount of occupied habitat
We calculated the net change in the amount of occupied habitat However, recall that in our models, pikas were allowed to colonize some currently suitable but unoccupied habitats by 2080 in the absence of any LUC or GCC (our "Dispersal Only" scenario) because we did not wish to assume that each species fully occupies all of its currently suitable habitat. Therefore, in determining whether LUC or GCC helps or hurts a species, we also calculated the net changes in the amount of occupied habitat incorporating the "Dispersal Only" scenario as the baseline (dashed lines in Figure 4, Figures S13 and S18), and we compared the results to those under the "LUC," "GCC" and "BOTH" scenarios. Using the "Dispersal Only" scenario as the baseline, we found more species under more scenarios suffer a net which we calculated the net changes in the amount of occupied habitat assuming the pikas are able to colonize all suitable habitats in the future (grey bars in Figure 4, Figures S13 and S18).
With unlimited dispersal, our models predicted larger net gains or smaller net losses of occupied habitats for all species, scenarios and pathway sets.

F I G U R E 3
The amount and elevations of colonized habitat of the five pika species under three scenarios: only land use change ("LUC" scenario), only climate change ("GCC" scenario) and both land use change and climate change ("BOTH" scenario). (a) The amount of colonized habitat, expressed as a percentage in relation to the amount of currently occupied habitat. (b) The elevations of colonized habitats. In (b), the upper and lower dashed lines in each facet indicate the elevational range limits (0.5%-99.5%) of currently occupied habitats; the upper and lower hinges of the boxplots correspond to the first and third quartiles of colonized habitat; the upper and lower whiskers extend from the hinges to the highest and lowest values no farther than 1.5 × IQR from the hinge (where IQR is the distance between the first and third quartiles) F I G U R E 4 Net change in amount of occupied habitat for each of the five pika species under three scenarios: only land use change ("LUC" scenario), only climate change ("GCC" scenario) and both land use change and climate change ("BOTH" scenario). Two dispersal scenarios are considered: one incorporating realistic dispersal limitations (black) and one that assumes unlimited dispersal abilities (grey). The net change is expressed as a percentage in relation to the amount of currently occupied habitat. The dashed lines in each facet indicate the net change in amount of occupied habitat in the absence of any LUC or GCC (with realistic dispersal)

| Elevational shifts of occupied habitat
We estimated the overall elevational range shifts for each species by comparing the elevational quantiles of habitat in the present and future. Under the "Dispersal Only" scenario, which assumes no LUC or GCC over time, we found downslope shifts occurring in four out of the five species due to gradual colonization of suitable but un-

| D ISCUSS I ON
We combined species distribution models with a dispersal model to explore the separate and combined impacts of land use change (LUC) and global climate change (GCC) on the distributions of five species of pikas on the Qinghai-Tibet plateau (QTP) of China. By doing so, we were able to explore how two of the most pervasive threats to global biodiversity are likely to affect the distribution of species, as well as the role that dispersal limitations may play in exacerbating negative impacts of GCC and LUC. We found that the relative magnitude of the impacts of land use change and climate change on habitat abandonment and colonization varied among pika species. We further found that the amount of habitat abandoned or colonized under a "BOTH" scenario that incorporates the effects of LUC and GCC was consistently less than the sum of the "LUC" plus the "GCC" scenarios, indicating that the impacts of LUC and GCC are not necessarily additive or exacerbating. When we calculated net changes in the amount of occupied habitat predicted by the "BOTH" scenario in 2080 relative to the current distributions of the five species, we found that three of five species experience a net gain of habitat under all pathway sets, which would suggest that the combination of LUC and GCC causes no net harm to them. However, using as our baseline a "Dispersal Only" scenario (i.e., one that does not assume that the five species have fully occupied all of their currently suitable habitats and will therefore continue to expand their ranges in the absence of any anthropogenic disturbances), we found that all five pika species experience a net loss of occupied habitat under some pathway sets and that O. curzoniae experiences a net loss under all pathway sets. Thus, the choice of a baseline largely determines how one views the overall impacts of LUC and GCC to these species. By comparing elevational shifts relative to the "Dispersal Only" scenario, we also showed that GCC will likely push most species to higher elevations. By running the models with and without realistic dispersal limitations, we found that incorporating dispersal limitations resulted in smaller net gains and larger net losses of occupied habitats for all five species.

| Impacts of LUC and GCC on habitat abandonment and colonization
Although LUC and GCC have been implicated in range shifts in many species (Gossner et al., 2016;Lawler et al., 2014;Román-Palacios & Wiens, 2020;Schwalm et al., 2016), there is no general consensus as to the relative importance of each in places where both are operating (which, it should be noted, includes much of the earth's land surface; see Jetz et al., 2007). Previous studies have argued that LUC tends to affect species at lower elevations because humans tend to transform the landscapes at lower elevations first and then move upslope (Elsen et al., 2020;Nogués-Bravo et al., 2008). Our modelling predicted that the "LUC" scenario will result in more abandoned habitat than the "GCC" scenario for all five species under some pathway sets. Moreover, the "BOTH" scenario, which incorporates both LUC and GCC, yielded results for the amount of abandoned and colonized habitat that are closer to the "LUC" scenario than to the "GCC" scenario for most species and pathway sets. This may be because a large human population has been resident on the QTP led for hundreds, even thousands, of years, resulting in widespread and ongoing F I G U R E 5 Elevational range shifts (m) under the "Dispersal Only," "LUC," "GCC" and "BOTH" scenarios. See text for details. Data are divided into facets for each of the five pika species. In each facet, x values represent nine quantiles (10th, 20th, … 90th) of the elevations of currently occupied habitats, and y values indicate the shifts in elevation of quantiles by 2080 LUC (Qi et al., 2020). That being said, our modelling also predicted that the "GCC" scenario will result in more abandoned habitat than the "LUC" scenario for four of the five species under some pathway sets. According to empirical studies, pikas maintain high body temperatures (e.g., 40.1°C for O. princeps) and metabolic rates (Smith & Weston, 1990;Wang et al., 2006), which may render them notably sensitive to climate warming, especially to extreme temperatures, contingent on any behavioural adaptations to buffer these effects (Mathewson et al., 2016).

The interacting effects of LUC and GCC on species distributions
have been demonstrated in global and regional studies of many species (Northrup et al., 2019;Peters et al., 2019;Radinger et al., 2016).
Here, we found that the amount of habitat abandoned or colonized under the "BOTH" scenario was consistently less than the sum under the "LUC" plus the "GCC" scenarios. Similar interactions have been found in other modelling studies of LUC and GCC involving multiple taxonomic groups (Kuemmerlen et al., 2015;Northrup et al., 2019;Radinger et al., 2016). Previous studies have suggested that synergistic interactions between LUC and GCC are more likely to occur in warm regions and in aquatic ecosystems (Mantyka-Pringle et al., 2012. The intrinsic interactions between land use and climate may also lead to predominantly overlapping impacts of LUC and GCC (Oliver & Morecroft, 2014), as we found in our study.

| Elevational patterns of abandoned and colonized habitats
Relative to the baseline "Dispersal Only" scenario, the "GCC" and "BOTH" scenarios pushed the elevational ranges of all five pika species upslope with few exceptions. As many species of animals are known to have shifted their distributions towards historically cooler regions in response to climate warming (Chen et al., 2011;Freeman et al., 2018;Galbreath et al., 2009), one might expect that GCC would reduce habitat suitability at lower elevations and generate increased habitat suitability at higher elevations, which is in line with our findings for pikas. Interestingly, while studies of other species in other parts of the world have argued that LUC tends to induce habitat loss at lower elevations, because people alter the landscapes at lower elevations first and then move upslope (Elsen et al., 2020;Nogués-Bravo et al., 2008), we did not see this with the pikas (Figure 2b).
That we did not see the expected pattern with the QTP pikas may reflect the fact that the QTP is a high-elevation plateau that has been subject to heavy human usage.

| Net changes in occupied habitats
As noted above, when net gains in occupied habitats are calculated relative to a "Dispersal Only" scenario, we found more species under more scenarios suffer a net loss of occupied habitat by 2080 than we do when using the current distributions of the species as our baseline (which assumes they have a static distribution). This highlights the value of considering potential range shifts relative to a baseline of no environmental change that nonetheless allows continued dispersal by the focal species, in order to avoid overestimating colonized habitat under environmental change scenarios. However, a "Dispersal Only" scenario should not be employed under all conditions. In this study, we assumed that even today each species of pika occupies only a portion of the suitable habitat available to it, defined as the suitable habitat within its IUCN polygon. This is realistic only when the IUCN polygons (or any other range data used to delineate current occupation) are relatively up to date (as is true for the pikas in this study) and reliable. In addition, a species' dispersal ability (as measured in field studies) may not translate directly into its ability to move in a given circumstance. For example, biological interactions such as the presence of competitors or predators could impede a species from moving into a particular area of seemingly suitable habitat, a possibility we did not consider in our dispersal model (Godsoe & Harmon, 2012;Trainor et al., 2014).
Nonetheless, poor dispersal abilities coupled with natural barriers are known to restrict species' range shifts driven by environmental changes (Lawler et al., 2013;Senior et al., 2019). Uncolonized suitable habitats exist under all our model scenarios for all five species (Figure 6, Figures S9 and S10). When pikas were assumed to be fully capable of colonizing all suitable habitats in the future (i.e., unlimited dispersal abilities, which is unrealistic), we observed larger net gains or smaller net losses of occupied habitats for all F I G U R E 6 Shifts in habitat occupancy for pika species under the combined effects of land use change and climate change ("BOTH" scenario). The geographical extent of each map is 61-130E, 17-62N five species (Figure 4, Figures S13 and S18: grey bars). This underscores the importance of considering species' dispersal limitations when predicting range shifts in the face of environmental change (Engler et al., 2012;Sales et al., 2019). By 2080, uncolonized, suitable habitats are predicted to appear at the edges of the QTP (Figure 6, Figures S9 and S10), where topographic barriers will block the movement of pikas; new areas of suitable habitat will also appear in places that are quite distant from the pikas' current distributional ranges and will also be likely inaccessible to them. However, these areas could serve as potential targets for conservation translocation, should it ever be necessary (Berger-Tal et al., 2020).

| C AV E AT S
Although we considered multiple IAMs and GCMs for LUC and GCC data, respectively, uncertainty in both land use and climate models remains a concern (Hewitt et al., 2016). Future studies could incorporate more IAMs and GCMs to improve presumed accuracy in predicting species range shifts under LUC and GCC. Of the approximately 28 species of pikas found in our bounding box study area, only five were included in our analysis because there were too few reliable records of the others for use in our models. Unfortunately, species with fewer occurrence records are likely to be rarer or more localized and therefore more vulnerable to extinction, reflecting what has been called a "rare species modeling paradox" (Lomba et al., 2010). We also may over-estimate future occupied areas by not considering biological factors beyond habitat type in our SDMs; these other factors, including the presence of competitors, predators, food plants and pathogens, can limit the occupation of abiotically suitable habitats (Peterson, 2011).

| CON CLUS I ON S AND MANAG EMENT IMPLI C ATI ON
Our results reveal the complex potential negative and beneficial effects of LUC and GCC to pikas on the QTP. The combined effects of LUC and GCC are not simply additive or synergistic, leading to complicated patterns of habitat gain and loss for each of the five species we modelled. Moreover, an overall assessment as to whether a given species will be "worse off" or "better off" in terms of the amount of occupied habitat in the future in the face of LUC and GCC is dependent on the choice of baseline (specifically, whether one believes the current ranges of species reflect full occupancy of reachable, suitable habitats or whether the species are still colonizing new areas).
In terms of management practice, our work highlights the value of an integrative consideration of LUC and GCC when predicting potential range shifts of species. In addition, the three pathway sets we considered in this study resulted in different range-shift predictions for the five pika species. We therefore advocate consideration of multiple socioeconomic pathways and emission scenarios when predicting the combined effects of land use change and climate change on species distributions so as to generate a range of possible outcomes. We also advocate consideration of a "Dispersal Only" scenario rather than reflexively assuming that all currently suitable habitats are occupied by the species in question. Finally, our results make clear the importance of including species' dispersal limitations when modelling range shifts due to environmental change.

ACK N OWLED G EM ENTS
We thank Paul R. Elsen for his comments on the manuscript. We

CO N FLI C T O F I NTE R E S T
The authors report no conflict of interest.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13408.

DATA AVA I L A B I L I T Y S TAT E M E N T
Species occurrence records and R code used in this research are archived in Dryad (Li et al., 2021)