Introduction

Habitat connectivity, which is defined as the degree to which habitat promotes or hinders species dispersal (Taylor et al. 1993), is a key factor for population persistence (Morelli et al. 2017). In the short terms, low habitat connectivity impedes ecological processes such as daily foraging, seasonal migration, successful reproduction, and juvenile dispersal (Rayfield et al. 2011; Braaker et al. 2014). In the long term, it increases extinction risk and decreases genetic diversity, because of population isolation, inbreeding depression, and demography stochasticity (Koen et al. 2014; Correa Ayram et al. 2016). Maintaining and restoring habitat connectivity is thus essential for biological conservation.

Habitat connectivity is a function of the area, quality, and configuration of habitat, as well as of the dispersal capability of species (Hodgson et al. 2009), and, hence, is sensitive to spatial variations in habitat, such as habitat loss, degradation, and fragmentation (Clauzel et al. 2015; Inoue and Berg 2017; Lechner et al. 2017). For example, the loss of habitat patches that serve as connectors for other patches can greatly reduce habitat connectivity and has an adverse impact on the abundance, diversity and stability of populations (Thompson et al. 2017; Peeters et al. 2020). Therefore, understanding the impact of habitat variations on habitat connectivity is imperative for effective connectivity conservation.

Climate change has been identified as the major driver of spatial variations in habitat (Mawdsley et al. 2009; Melles et al. 2011). Given that species’ distribution reflects their climatic niche, species will be forced to adjust their distribution in response to climate change, to match the new suitable areas (Inoue and Berg 2017). Such an adjustment may alter not only the range but also the configuration of their habitat. For instance, Dilts et al. (2016) have predicted that the habitat of the Mohave ground squirrel (Xerospermophilus mohavensis) in California will be greatly reduced and fragmented into three distinct clusters under future climate conditions. These alterations in habitat might have a profound impact on species dispersal patterns, leading to the loss of key components for connected habitats, the extension of dispersal routes, and the consequent decrease in habitat connectivity (Leblond et al. 2016; Morrison et al. 2016). For example, Peeters et al. (2020) have demonstrated that the continued loss of sea ice, which serves as connectors between island systems, would decrease habitat connectivity for Arctic ungulates. However, while changes in habitat distribution with respect to climate change have been examined for a wide range of taxa (e.g., Sharma et al. 2009; Attorre et al. 2011; Bambach et al. 2013; Fan et al. 2014; Luo et al. 2015; Qin et al. 2019; Yu et al. 2019; Singh et al. 2020), the subsequent changes in habitat connectivity are understudied.

In addition, although connectivity conservation has long been recommended as a critical adaptation strategy against climate change (Heller and Zavaleta 2009; Hannah 2011), it has rarely been integrated into conservation planning such as protected area (PA) design (Magris et al. 2014). Moreover, because the contribution of an area to habitat connectivity may change as species shift their range, PA design should also consider the uncertainties of climate-driven connectivity changes (Albert et al. 2017). However, little work has been done to examine the robustness of PAs against climate change in terms of connectivity conservation (but see Coleman et al. 2017).

The Tibetan Plateau, known as the third pole of the world and the water tower of Asia, plays a considerable role in maintaining the ecological security of the Northern Hemisphere (Xu et al. 2009). With its vast expanse and heterogeneous geography, the plateau sustains a distinct biome with numerous endemic species, especially endemic ungulate species such as the wild yak (Bos mutus), kiang (Equus kiang), Tibetan antelope (Pantholops hodgsonii), and Tibetan gazelle (Procapra picticaudata). However, along with other sensitive and fragile ecosystems, the Tibetan Plateau has been undergoing unprecedented climate change (Kang et al. 2010). According to the IPCC’s Fifth Assessment Report, the average annual mean temperature will become 0.9–4.9 °C warmer and the average annual precipitation will change by − 1 to + 32% on the Tibetan Plateau over the next 100 years (IPCC 2013). Species native to the Tibetan Plateau, physiologically and phenologically adapted to the local climate, will be forced to adjust their distribution to keep pace with these environmental changes (Sutherland et al. 2018; Wang et al. 2020). For instance, ungulates, which are keystone species in the Tibetan Plateau ecosystem, have been predicted to lose 30–50% of their distribution area and experience an average poleward range shift of 300 km because of climate change (Luo et al. 2015). Such processes would lead to dramatic changes in habitat configuration and, thus, habitat connectivity, which may have a profound impact on the long-term persistence of these species. However, to our knowledge, no study has assessed the consequent changes in the habitat connectivity of ungulates under climate change on the Tibetan Plateau. In addition, although the existing PAs cover a sizable portion of the Tibetan Plateau, they have rarely accounted for the need for connectivity conservation, nor have they considered the uncertainties in habitat connectivity associated with climate change. Thus, there is an urgent need to examine the effectiveness of PAs on the Tibetan Plateau in terms of connectivity conservation for ungulates and respond appropriately to the challenges imposed by climate change.

In this study, we focused on four endemic and endangered ungulate species of the Tibetan Plateau, namely, the wild yak, kiang, Tibetan antelope, and Tibetan gazelle. These ungulates are the most representative species on the Tibetan Plateau (Jiang et al. 2018), and, hence, are extremely important for understanding the adaptation of the Tibetan Plateau ecosystem to climate change. Specifically, we first developed ecological niche models (ENMs) to assess the spatial variations in habitat under climate change for each focal species. Then, we conducted a network-theoretical connectivity analysis to estimate consequent changes in habitat connectivity. Furthermore, we characterized the dispersal pattern of each species in each climate scenario based on the circuit theory. Finally, we used gap analyses to estimate the contribution of the current PAs on the Tibetan Plateau to the conservation of the four species across different climate scenarios.

Our objectives were as follows: (1) to assess how habitat connectivity change following climate-driven spatial variations in habitat for four ungulate species endemic to the Tibetan Plateau; (2) to identify key areas for maintaining habitat connectivity for the four ungulate species; and (3) to evaluate the robustness of PAs on the Tibetan Plateau for the conservation of the four species under climate change.

Methods

Study area and species occurrence records

Our study area covered the entire Tibetan Plateau (25.99°–39.83°N, 73.45°–104.67°E), extending from the Kunlun and Qilian mountains in the north to the greater Himalayas in the south, and from the Pamir and Hindu Kush Himalayas in the west to the Hengduan Mountains in the east (Fig. 1). The study area was spread across four provinces and autonomous regions in southwestern China, namely, the Qinghai Province, Sichuan Province, Tibet Autonomous Region, and Xinjiang Autonomous Region, with a total area of approximately 2.6 × 106 km2. We converted the entire study area into 2,589,784 grid cells, each with a 1 × 1 km resolution, and all the analyses were conducted based on these grid cells.

Fig. 1
figure 1

Location of the study area and the occurrence records of the four ungulate species

Full size image

The occurrence records of the four ungulate species, namely, the wild yak, kiang, Tibetan antelope, and Tibetan gazelle, were mostly derived from the Second National Survey on Terrestrial Wildlife Resources in China, which used transect lines for investigation. We also collected records from available literature by searching the scientific names of these four species in Google Scholar (Su et al. 2014; Wu et al. 2014; Dong et al. 2015; Luo et al. 2015; Liang et al. 2016). Records that lacked precise coordinates or were repeated in the same grid cell were removed to ensure modeling accuracy and eliminate potential bias of spatial autocorrelation (Shcheglovitova and Anderson 2013; Luo et al. 2015). In total, 2551 occurrence records, including 2386 from surveys and 165 from the literature, were used in this study. These records divided into 390 for the wild yak, 640 for the kiang, 868 for the Tibetan antelope, and 653 for the Tibetan gazelle (Fig. 1).

Environmental variables

We used 16 environmental variables falling into four groups, namely, climate, habitat, topography, and anthropogenic influence, to develop ENMs for the four species (Table 1).

Table 1 Environmental variables used in the ecological niche models for the four ungulate species
Full size table

The climate variables included the annual mean temperature, temperature seasonality, maximum temperature of the warmest month, minimum temperature of the coldest month, temperature annual range, annual precipitation, precipitation of the wettest month, precipitation of the driest month, and precipitation seasonality. All the climate variables were derived from the WorldClim database v 2.0 (www.worldclim.org).

The habitat variables included land cover, the normalized difference vegetation index (NDVI), and the distance to the nearest surface water. Land cover was determined using the GlobeLand30 database (www.globallandcover.com). The NDVI was obtained from the Ministry of Environment Protection of the People’s Republic of China (www.zhb.gov.cn), and the average of 3 years (i.e., 2012–2014), which corresponded to the period of the field surveys, was used. The distance to the nearest surface water was calculated as the linear distance between the center of any given pixel and the nearest surface water (i.e., lakes or rivers).

The topography variables included elevation and slope, which were calculated from GDEM 30-m digital elevation data that were derived from the International Scientific and Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences (www.gscloud.cn).

The anthropogenic influence variables included the human footprint index and distance to anthropogenic disturbances. The human footprint index was obtained from the Global Footprint Network (www.footprintnetwork.org). The distance to anthropogenic disturbances was calculated as the linear distance between the center of any given pixel and the nearest potential anthropogenic disturbance source (i.e., cities, villages, railways, and roads).

Future climate scenarios

We derived future climate scenarios using four global climate models (GCMs), namely, HadGEM2-ES, IPSL-CM5A-LR, MIROC-ESM-CHEM, and NorESM1-M, in combination with four representative concentration pathways (RCPs) emission scenarios, namely, RCP2.6, RCP4.5, RCP6, and RCP8. These four GCMs are widely used in ENMs because of their capacity to project a wide range of temperatures and precipitation levels (Caminadea et al. 2014), and have a good predictive performance for the Tibetan Plateau (Su et al. 2013). The corresponding data were obtained from the WorldClim dataset (www.worldclim.org) for the 2050s (i.e., the average for the years 2041–2060) and the 2070s (i.e., the average for the years 2061–2080).

For each climate variable, we averaged the sixteen projections (i.e., four GCMs × four RCPs) for both periods (i.e., the 2050s and 2070s) in a GIS environment. Notably, there is no dataset available for the future habitat and anthropogenic influence variables of the Tibetan Plateau, and these variables are contingent on a range of socioeconomic drivers; hence, any simple estimations, such as extrapolations from past trends, are likely to be misleading (Hu et al. 2010). Therefore, following Luo et al. (2015), we assumed that habitat and anthropogenic influence variables were stable and used the current data to represent the future situation.

Ecological niche models

As strict systematic sampling is difficult to conduct on most of the Tibetan Plateau because of harsh environmental conditions, absence data were unavailable for the four species. Therefore, we used the maximum entropy model (Maxent) with presence-only data to develop our ENMs. We built ENMs for each species based on their occurrence records and the aforementioned environmental variables. Then, we projected these ENMs under future climate scenarios to predict the future habitats where these ungulates might persist.

To reduce the influence of collinearity between variables in the models, we filtered some correlated environmental variables for each species. To this end, we first calculated the Pearson correlation coefficient between all pairs of environmental variables. Then, we produced an initial model to examine the percentage contribution of every variable for each species, and, finally, removed the variable with a lower percentage contribution in pairs with a correlation coefficient greater than 0.8 (Luo et al. 2015; Table 1, Online Appendix Table 1).

ENMs were produced with default settings (i.e., feature selection automatic; regularization multiplier: 1; maximum iterations: 500; and convergence threshold: 10–5) and logistic output format (range from 0 to 1, representing the habitat suitability) in the Maxent software (v 3.3.3). The model accuracy was evaluated with a 25% random subset of the occurrence records, using the area under the curve (AUC) of the receiver operating characteristic curve, which ranged from 0.5 (random accuracy) to 1.0 (perfect discrimination; Phillips et al. 2006). For each species and period (i.e., the current period, 2050s, and 2070s), we ran ten cross-validation replicates and averaged the outcomes to obtain a more robust outcome of the model performance (Bambach et al. 2013).

Spatial variation assessment

In order to acquire core habitat patches with high suitability and spatial integrity from the results generated by the ENMs, we first conducted a hotspot analysis with the Getis-Ord \(G_{i}^{*}\) statistic to identify significant (p < 0.05) suitability clusters (Bagstad et al. 2016; Schank et al. 2017). The Getis-Ord \(G_{i}^{*}\) statistic was computed in a GIS environment, as follows (Valck et al. 2016):

$$G_{i}^{*} = \frac{{\mathop \sum \nolimits_{j = 1}^{n} w_{ij} x_{j} - \overline{X}\mathop \sum \nolimits_{j}^{n} w_{ij} }}{{S\sqrt {\frac{{n\mathop \sum \nolimits_{j = 1}^{n} w_{ij}^{2} - \left( {\mathop \sum \nolimits_{j = 1}^{n} w_{ij} } \right)^{2} }}{n - 1}} }}$$

where \(x_{j}\) represents the grid cell j, \(w_{ij}\) is the spatial weight between the grid cells i and j, n is the total number of grid cells, \(\overline{X} = \frac{{\mathop \sum \nolimits_{j = 1}^{n} x_{j} }}{n}\), and \(S = \sqrt {\frac{{\mathop \sum \nolimits_{j = 1}^{n} x_{j}^{2} }}{n} - \left( {\overline{X}} \right)^{2} }\). We then used a 0.5 probability threshold, following Buuveibaatar et al. (2016), to exclude low suitability areas from the suitability clusters. Patches with an area of less than 100 km2 were removed because they might not be large enough to sustain ecological processes of ungulate herd in the Tibetan Plateau (Hu et al. 2018). Besides, the removal of these patches would lead to a substantial decrease in computational time (Dickson et al. 2014).

For each species, we calculated several metrics to characterize the spatial variations in habitat in response to climate change. To assess climate-driven range shift, we calculated the mean latitude and mean longitude by averaging the latitudes and longitudes of the weighted centroids of each core habitat patch. To estimate changes in habitat amount, we calculated the number of patches and the total weighted area of all patches. Notably, we used weighted area, that is, the sum of the grid values within core habitat patches, rather than area, to account for the uneven distribution of habitat suitability (Dilts et al. 2016). To explore the potential effects of climate change on habitat configuration, we calculated the area weighted mean patch size by summing the weighted area across all patches multiplied by the proportional abundance of the corresponding patch (i.e., the weighted area of the patch divided by the total weighted area), as well as the proportion of the largest patch by dividing the weighted area of the largest patch by the total weighted area.

Habitat connectivity estimation

To quantify and interpret changes in habitat connectivity following climate-driven spatial variations, we used the equivalent connected area (ECA) index, a network theoretical metric that quantifies connectivity by measuring the size of a single habitat patch that would provide the same probability of connectivity (i.e., the probability that two randomly placed points in the landscape fall into habitat areas that are reachable from each other) than the actual habitat configuration in the studied landscape (Saura et al. 2011). This index was calculated as follows:

$$ECA = \sqrt {\mathop \sum \limits_{i = 1}^{n} \mathop \sum \limits_{j = 1}^{n} a_{i} a_{j} p_{ij}^{*} }$$

where \(a_{i}\) and \(a_{j}\) correspond to the weighted areas of the core habitat patches i and j respectively, \(p_{ij}\) represents the probability of direct dispersal between the patches i and j, and \(p_{ij}^{*}\) is the maximum product probability of all possible dispersal between the patches i and j, including both direct and indirect connections (Saura et al. 2011). In our study, \(p_{ij}\) was translated from the length of the least cost path (LCP), using a negative exponential decay function with 0.5 for a corresponding distance range from 25 to 75% of the distance between the most distant patches (Bunn et al. 2000; Hanski and Ovaskainen 2000; Blazquez-Cabrera et al. 2014). The LCP was calculated in Linkage Mapper (McRae and Kavanagh 2011) with resistance surfaces created from the inverse of our ENM outcomes (Online Appendix Fig. 2), and the ECA calculations were performed using the Conefor 2.6 software (Saura and Torné 2009).

Because the ECA index has the same unit as core habitat patches, the relative change in ECA (dECA) can be directly compared with the relative variation in the total weighted areas of the patches (dWA) under the same climate change scenarios. Following Dilts et al. (2016), we defined this comparison as relativized ECA (rECA), which was calculated by dividing the dECA by the dWA. The rECA can provide straightforward insights into how connectivity changes occur in relation to changes in habitat area. An rECA value greater than 1 indicates that habitat changes result in a disproportionately large change in habitat connectivity, while a value lower than 1 indicates connectivity changes due to random habitat changes (Saura et al. 2011; Dilts et al. 2016).

To gain a more explicit understanding of how climate change could affect the dispersal patterns of the four species, and to identify key areas for their habitat connectivity under different climate scenarios, we used concepts from the circuit theory, a connectivity model integrating network theory and random walks (McRae et al. 2008). In the circuit theory, the entire landscape is treated as a conductive surface (i.e., a raster grid consisting of nodes and resistors), through which current (analogous to species movement) will flow when it is injected into a source node and a target node is tied to the ground. The current density of each intervening grid cell can be interpreted as the probability that a random walker would pass through the cell on its way to the target node. The more alternative pathways that exist between two nodes, the lower current densities of the intervening grid cells, and vice versa. High current density areas, the so-called “pinch points”, denote areas where species have a high likelihood or necessity of passing during their dispersal as alternative pathways are not available, and, thus, represent key areas for connectivity (McRae et al. 2008; Rayfield et al. 2011; Braaker et al. 2014).

We calculated the spatial distribution of the current density for each species under each climate scenario using the Circuitscape 4.0 software (McRae et al. 2013) with conductive surfaces derived from the ENM outcomes. To account for the heterogeneity in habitat suitability and estimate the current density within the core habitat patches, we used weighted centroids as nodes, rather than the entire patch (Dickson et al. 2017). We used Circuitscape’s pairwise mode iterated across all centroid pairs, injecting a 1 A current into one centroid and setting the other to the ground. Current densities were added up across all iterations to produce a cumulative current density map, which reflected the contribution of each grid cell for maintaining the connectivity among all the possible pathways across the entire landscape.

Gap analysis

For each climate scenario, we averaged the suitability maps (Guerin and Lowe 2013) and summed the cumulative current density maps (Dickson et al. 2017) of all the ungulate species, to identify key areas for multi-species suitability and connectivity respectively. Moreover, we overlaid the multi-species suitability and connectivity maps on the PAs of the Tibetan Plateau, and summed the suitability value and current density within each PA to assess the contributions of these PAs for ungulate conservation and identify potential conservation gaps. The extent of the PAs on the Tibetan Plateau was extracted from the World Database on Protected Areas (WDPA; www.protectedplanet.net).

Results

Spatial variations in habitat

All ENMs performed well in predicting the habitats of the four ungulate species, with AUC values of more than 0.8. Most predicted habitats were distributed in the western part of the Tibetan Plateau, where alpine desert-steppe vegetation predominates (Fig. 2, Online Appendix Fig. 1). For all species, the core habitat patches were predicted to shift northwards, with average increases in the mean latitude of 0.75° and 1.08° by the 2050s and the 2070s, respectively (Figs. 2, 3). The core habitat patch tended to decrease for all species, with average decreases in the patch number of 20.82% by the 2050s and 23.44% by the 2070s, and average decreases in the weighted area of 28.22% by the 2050s and 33.25% by the 2070s. As for variations in habitat configuration, they were not uniform across species. For the wild yak and the kiang, both the proportion of the largest patch and the area weighted mean patch size decreased in response to climate change. Conversely, the two metrics increased with climate change for the Tibetan antelope and the Tibetan gazelle.

Fig. 2
figure 2

Core habitat patches of the wild yak (a, e, i), Tibetan antelope (b, f, j), kiang (c, g, k), and Tibetan gazelle (d, h, l) under current (a–d) and future climate scenarios for the 2050s (e–h) and 2070s (i–l)

Full size image
Fig. 3
figure 3

Habitat characteristics of the four ungulate species under three climate scenarios. The mean latitude, mean longitude, patch number, weighted area, area weighted mean patch size, and proportion of the largest patch are displayed for each species and under each climate scenario

Full size image

Changes in habitat connectivity

The ECA index increased with dispersal distance, indicating a positive correlation between the habitat connectivity and dispersal capability of the four species. For all four species, the ECA indices of the future habitats were significantly lower than that of the current habitat at all dispersal distances (Fig. 4). The average dECA was predicted to range from − 0.25 to − 0.47 by the 2050s and from − 0.25 to − 0.62 by the 2070s (Table 2). When relativizing the relative change in habitat connectivity to the change in weighted area, all four species consistently displayed connectivity loss that was higher than habitat loss, with rECA values greater than 1 under all climate change scenarios (Table 2). The largest difference between connectivity loss and habitat loss was detected in the Tibetan gazelle, with an rECA value of 1.65 in the 2050s and 1.92 in the 2070s, which indicates that the predicted connectivity loss would exceed the habitat loss by more than 50%.

Fig. 4
figure 4

Equivalent connected area (ECA) index values for the habitats of the four ungulate species under three climate scenarios. The ECA was calculated based on dispersal distances ranging from 25 to 75% of the distance between the most distant patches

Full size image
Table 2 Average relative variation of the equivalent connected area (ECA) index (dECA), relative variation of the total weighted areas (dWA), and relativized ECA index (rECA), under two climate change scenarios, for the habitats of the four ungulate species
Full size table

Changes in key connectivity areas

The circuit model showed that the current flows of the four ungulate species were mostly concentrated in the northwestern part of the Tibetan Plateau, including northern Tibet, western Qinghai, and southern Xinjiang (Fig. 5a–d). Within this range, the highest current densities of the four species were observed along the west–east orientated valleys of the Qiangtang Basin, which is located between the Kunlun Mountains and the Gangdise Mountains. The current flows of the wild yak, kiang, and Tibetan antelope were predicted to shift northwestwards in response to climate change, with the highest current densities emerging in northwestern Tibet and southern Xinjiang (Fig. 5e–g, i–k). In contrast, the current flows of the Tibetan gazelle were predicted to be separated into two parts, because of the absent of patches in the Selincuo Lake Basin, under the future scenarios (Fig. 5h, j).

Fig. 5
figure 5

Current density of the wild yak (a, e, i), Tibetan antelope (b, f, j), kiang (c, g, k), and Tibetan gazelle (d, h, l) under current (a–d) and future climate scenarios for the 2050s (e–h) and 2070s (i–l)

Full size image

Dynamics of conservational contribution

The multi-species suitability and connectivity maps exhibited a generally similar pattern, with key suitability and connectivity areas concentrated in the central part of Qiangtang Basin (Figs. 6, 7). Both key suitability and connectivity areas were predicted to shift to the western part of the Qiangtang Basin and the Altun Mountains under future climate scenarios.

Fig. 6
figure 6

Habitat suitability for all four ungulate species under current and future climate scenarios for the 2050s and 2070s. The boundaries of the protected areas on the Tibetan Plateau are shown. The letters indicate the five protected areas contributing the most to the connectivity conservation of the habitat of the ungulate species. A the Qiangtang Nature Reserve; B the Arjinshan Nature Reserve; C the Sanjiangyuan National Park; D the Selincuoheijinghe Nature Reserve; E the Qomolangma Nature Reserve

Full size image
Fig. 7
figure 7

Habitat connectivity for all four ungulate species under current and future climate scenarios for the 2050s and 2070s. The boundaries of the protected areas on the Tibetan Plateau are shown. The letters indicate the five protected areas contributing the most to the connectivity conservation of the habitat of the ungulate species. A the Qiangtang Nature Reserve; B the Arjinshan Nature Reserve; C the Sanjiangyuan National Park; D the Selincuoheijinghe Nature Reserve; E the Qomolangma Nature Reserve

Full size image

There are 109 PAs on the Tibetan Plateau that could contribute to the conservation of the four ungulate species in terms of habitat suitability and connectivity. Among them, we found that the Qiangtang Nature Reserve accounted for the highest contribution to suitability and connectivity (27.49 and 34.23%, respectively), followed by the Arjinshan Nature Reserve (5.67 and 7.99%, respectively) and the Sanjiangyuan National Park (8.70 and 7.60%, respectively; Figs. 6, 7, Table 3). Notably, 53.39 and 46.64% of the land that could contribute to the conservation of the four ungulate species in terms of habitat suitability and connectivity, respectively, is not under protection. These unprotected areas would increase to 55.30and 52.56% by the 2050s and 55.78 and 53.27% by the 2070s, as the four species would shift their range northwards under the future climate scenarios (Table 3).

Table 3 Top five protected areas with the highest proportion of accumulative suitability and accumulative current density on the Tibetan Plateau
Full size table

Discussion

Alteration of the habitat of ungulates by climate change on the Tibetan Plateau

Our study predicted significant spatial changes in the habitats of ungulate species on the Tibetan Plateau because of climate change. In general, our four focal species were predicted to shift their range northwards and lose a considerable part of their core habitats under future climate conditions. Such predictions are consistent with the findings of Luo et al. (2015). However, the four species showed different responses in terms of their habitat configuration. The habitats of the wild yak and kiang were predicted to fragment with climate change, which is similar to what has been anticipated for a variety of taxonomic groups (e.g., Wasserman et al. 2012; Jackson et al. 2016; Kang et al. 2016; Inoue and Berg 2017). Conversely, the habitats of the Tibetan antelope and Tibetan gazelle were predicted to be more cohesive under future climate scenarios. These differences are probably attributed to the degrees of fragmentation of the habitats of the four species (Thompson et al. 2017). Compared to the current habitats of the wild yak and kiang, those of the Tibetan antelope and Tibetan gazelle are more fragmented and consist of many small patches, which are considered more vulnerable to environmental changes (Tscharntke et al. 2002; Archibald et al. 2011). Therefore, rather than the fragmentation of large habitats, the impact of climate change on the habitats of the Tibetan antelope and Tibetan gazelle might rather consist in the loss of small patches, which increases the mean patch size and largest patch proportion.

Connectivity changes following climate-driven habitat variation

Our study predicted a remarkable decrease in the habitat connectivity for all four ungulate species following changes in habitat distribution and configurations. This decrease may be related to the habitat loss of these species because such loss could lead to fewer habitat patches available for connection (Saura et al. 2011). However, the connectivity decreases outpaced the habitat loss for all species, especially for the Tibetan gazelle, whose connectivity reduction was projected to be nearly twice as much as the habitat loss per unit area by the 2070s. Notably, the habitat of the Tibetan gazelle is more fragmented, with numerous small patches, which may disappear under the future climate conditions. Our results suggest that some of the habitat patches of the Tibetan gazelle habitat that were predicted to disappear are located in a crucial topological position. That is, their absence may make the originally connected habitats become unreachable, and, thus, lead to an impact on connectivity disproportionately larger than that expected solely based on the area they provide. Such an interpretation aligns with our circuit model, which predicted that current flows between the eastern and western habitats of the Tibetan gazelle would be inhibited due to the absence of patches in the Selincuo Lake Basin under the future scenarios. Moreover, as indicated by Saura et al. (2014), such an impact could hardly be compensated through other measures, such as increasing the habitat amount or quality. Therefore, conservation efforts in relation to the impact reduction of climate change should focus not only on the fragmentation and shrinkage of large habitat patches, but also on the disappearance of the small patches that are essential to connect the entire habitat. The restoration of such essential patches can greatly improve the connectivity performance of habitat, which would be a great means of adaptation to climate change. For example, restoring the degenerate habitats of the Selincuo Lake Basin could help enhance the habitat connectivity of the Tibetan gazelle, and facilitate the dispersal of this species between the western and eastern Tibetan Plateau.

Key areas for the habitat connectivity of ungulates on the Tibetan Plateau

The circuit theory is a useful tool to model species dispersal patterns. When interpreting the maps resulting from circuit models, the most important concern is related to pinch points, that is, areas with high current density due to the lack of alternative pathways nearby (McRae et al. 2008). Our study found that pinch points for all four ungulate species were mostly presented in the central part of the Qiangtang Basin. Our results indicate that this area serves as a “traffic hubs” where all four species need to pass by during their dispersal on the Tibetan Plateau. For example, it has been reported that Tibetan antelopes regularly travel from northwestern Tibet, through the central part of the Qiangtang Basin, to the Kekexili and Arjinshan regions to lamb every year (Zhuge et al. 2015). Therefore, the center of the Qiangtang Basin is crucial for maintaining connectivity for all four ungulate species, and deserve to be granted conservation priority, because once it is blocked, the dispersal and gene flow of all four species will be inhibited. This inhibition would decrease population viability, with serious consequences for the biodiversity of the Tibetan Plateau. In response to climate change, the pinch points were predicted to shift northwestwards to the Altun Mountains and the western part of the Qiangtang Basin. This prediction is consistent with the results of Luo et al. (2015), who found an increase in the species richness of ungulates in the northwestern part of the Tibetan Plateau, which indicates that the importance of this area for the conservation of ungulates would increase as more ungulates immigrate to it.

Connectivity priorities versus species richness priorities

Identifying priority areas is an important strategy for conservation planning, such as PA design. Species richness has long been used as an indicator to identify areas with a high conservation value (Luo et al. 2015; Liang et al. 2017; Quan et al. 2017). However, the priorities defined by species richness may diverge from those related to connectivity (Crouzeilles et al. 2013; Albert et al. 2017; Fung et al. 2017), because the areas crucial for species dispersal do not necessarily coincide with the areas suitable for species survival (Braaker et al. 2014). For example, our study indicates that the valley in Taxkorgan will become an important corridor to connect the suitable habitats of the kiang in northern Tibet and the Pamirs by the 2050s and the 2070s, even though the valley itself would not be a suitable habitat for the kiang. Notably, species richness reveals areas that are suitable for multispecies survival, whereas habitat connectivity highlights areas indispensable to the dispersal of these species assemblages. Consequently, conservation planning that merely represents species richness priorities may miss the information provided by connectivity priorities, and, thus, may not ensure long-term persistence and viability of the species assemblage (Magris et al. 2014).

Potential gaps in ungulate connectivity conservation

The Tibetan Plateau is the region with the largest PA coverage in China. In particular, in the hinterland of the plateau, three adjacent vast PAs, namely, the Qiangtang Nature Reserve, Arjinshan Nature Reserve, and Sanjiangyuan National Park, compose one of the largest PA groups in the world, making up over 20% of the total area of the plateau. However, our study found that more than 45% of the areas that could contribute to the conservation of the four ungulate species in terms of habitat suitability and connectivity are not under protection. Specifically, in western Gangdise, a large area that is key for habitat suitability and connectivity of ungulates, including several pinch points extending southward along the valley into the western Himalayas, was exposed outside the southern boundary of the Qiangtang Nature Reserve. These areas should be considered as conservation gaps that deserve particular attention in future conservation planning. Moreover, as the four species were predicted to shift their range northwards in response to climate change, more key suitability and connectivity areas will fall outside the boundaries of existing PAs, which may further expand the aforementioned conservation gaps.

Adjusting the boundary and strengthening the management of certain key PAs will be an effective way to address the conservation gaps. For example, expanding the southwestern border of the Qiangtang Nature Reserve will help incorporate the new key suitability and connectivity areas into conservation efforts. In addition, since the Qiangtang Nature Reserve, Arjinshan Nature Reserve and Sanjiangyuan National Park play an important role in ungulate conservation in terms of both habitat suitability and connectivity, more conservation efforts, such as controlling graziery expansion and human activities and establishing biological corridors, are urgently needed to secure the key suitability and connectivity areas in these PAs.

Conclusions

In the face of ongoing climate change, species are forced to alter their habitat to adapt to the new climate pattern, which may have a profound impact on their habitat connectivity. Our case study suggests that the habitat connectivity of the wild yak, kiang, Tibetan antelope, and Tibetan gazelle will significantly decrease following climate-driven habitat variations on the Tibetan Plateau. Moreover, our study revealed that even a small change in habitat could cause a disproportionate connectivity loss if this alteration affected a location that is essential for maintaining connectivity within the entire habitat. In this context, we highlighted the importance of the central Qiangtang Basin for the habitat connectivity of all four ungulate species, which needs particular attention in conservation planning. We also found notable conservation gaps between key suitability and connectivity areas and the existing PAs on the Tibetan Plateau. These gaps could further expand as the connectivity priorities shift northwestwards because of climate change. Adjusting the extent and enhancing the management of specific key PAs may help address these conservation gaps.

Given that the dynamics of habitat (e.g., land cover) and anthropogenic influence (e.g., human footprint) factors associated with climate change were not included in the present study, our results may be somewhat conservative. Notably, this could lead to inadequate decisions on conservation management or planning. To improve the success of conservation, we recommend: (a) exploring the combined effects of climate, land cover, and anthropogenic influences on species; (b) incorporating more detailed species information, such as dispersal limitation, biotic interactions, and population density, into connectivity analysis; and (c) conducting fine-scale monitoring with camera traps and GPS collaring to verify the proposed connectivity priority.