Skip to main content

Combining the responses of habitat suitability and connectivity to climate change for an East Asian endemic frog



Understanding the impacts of past and contemporary climate change on biodiversity is critical for effective conservation. Amphibians have weak dispersal abilities, putting them at risk of habitat fragmentation and loss. Both climate change and anthropogenic disturbances exacerbate these risks, increasing the likelihood of additional amphibian extinctions in the near future. The giant spiny frog (Quasipaa spinosa), an endemic species to East Asia, has faced a dramatic population decline over the last few decades. Using the giant spiny frog as an indicator to explore how past and future climate changes affect landscape connectivity, we characterized the shifts in the suitable habitat and habitat connectivity of the frog.


We found a clear northward shift and a reduction in the extent of suitable habitat during the Last Glacial Maximum for giant spiny frogs; since that time, there has been an expansion of the available habitat. Our modelling showed that “overwarm” climatic conditions would most likely cause a decrease in the available habitat and an increase in the magnitude of population fragmentation in the future. We found that the habitat connectivity of the studied frogs will decrease by 50–75% under future climate change. Our results strengthen the notion that the mountains in southern China and the Sino-Vietnamese transboundary regions can act as critical refugia and priority areas of conservation planning going forward.


Given that amphibians are highly sensitive to environmental changes, our findings highlight that the responses of habitat suitability and connectivity to climate change can be critical considerations in future conservation measures for species with weak dispersal abilities and should not be neglected, as they all too often are.


Changing climatic conditions have begun altering or removing the conditions that are necessary for many species to exist on the Earth’s surface, leaving a clear fingerprint on global biodiversity [1,2,3]. Climate change induces shifts in the geographic ranges, abundance, phenology, individual behaviour or physiology, and genetic diversity of natural populations [4,5,6,7,8]. Furthermore, the rapid warming rates projected for the planet [1, 9] are expected to doom many species to extinction because climatically suitable habitats may disappear or become inaccessible due to geographic barriers or conditions that make it unsuitable for species to disperse [2, 3]. Regardless of other factors, the impacts of climate change on species are becoming increasingly apparent [5, 10]. It is clear that immense challenges to conservation and environmental management in the face of climate change exist [7, 11, 12].

Habitat degradation, fragmentation and loss are the primary causes of biodiversity loss in most of the world’s ecosystems [13]. Habitat fragmentation/isolation impedes dispersal, increasing inbreeding, genetic drift, and, hence, the likelihood of local extirpation due to demographic stochasticity [14, 15]. These factors might limit species’ abilities to adapt to environmental changes and can thus cause extinctions [13]. Habitat connectivity, the degree to which a landscape facilitates or impedes dispersal among habitat patches, confers ecosystems with greater resilience towards disturbances and thereby facilitates population viability [16]. As a consequence of the effects of climate change and human activities, habitat connectivity is often degraded, causing negative impacts to global biodiversity [12, 17].

Amphibians are highly sensitive to environmental changes due to their low mobility and strict physiological constraints [2, 18]. They are at disproportionately high risk of extinction, with many species’ populations rapidly declining. Amphibians can thus act as useful indicators of environmental changes across space and time [2, 18, 19]. One amphibian species that is facing a dramatic population decline is the giant spiny frog (Quasipaa spinosa), an endemic species to East Asia [20, 21]. Over-harvesting, observed distribution shrinkage, and ongoing habitat destruction and degradation have caused wild populations of Q. spinosa to significantly decrease, with an estimated decline of more than 30% over the last 15 years (roughly three Q. spinosa generations) [20, 22]. Their main range falls within a densely populated region of China, where, despite a long history of human occupation, recent developments have put new pressure on the species. There are currently no protective regulations for Q. spinosa except for collection prohibitions within nature reserves since the frog is not yet listed as a National Protected Wild Animal Species under Chinese law. All these factors create an urgent need to better understand the interactive threats of climate change and fragmentation to this species.

Here, we inferred the impacts of historical and future climate changes on the landscape connectivity of habitats of Q. spinosa. We simulated the environmental suitability for this species over time and measured its potential distribution as a function of the varying climatic conditions from the last interglacial period (LIG, 120–140 kyr BP) to the projected future climate, i.e., for the years 2050 and 2070. Then, we characterized landscape connectivity among habitat patches and determined how climate change could affect the dispersal of individuals. Overall, these efforts can provide novel insight into how habitat connectivity contributes to conserving amphibian species in a changing world.


Study species

Quasipaa spinosa is a large-sized anuran. Adult males are characterized by keratinized spines on their chests that appear during breeding seasons [22]. The distribution of Q. spinosa covers 12 provincial-level regions across central, southern and southwestern China (including Yunnan, Guizhou, Guangxi, Guangdong, Hong Kong, Fujian, Jiangxi, Hunan, Hubei, Anhui, Jiangsu and Zhejiang) and northern Vietnam [20, 22]. It inhabits rocky streams in evergreen forests and open countrysides/fields on hills and mountains. Its elevation range is 200–1500 m a.s.l. [20, 22, 23]. This species is considered an important food and medicinal resource according to traditional views [22]. It has been collected for consumption throughout its range for many decades, a practice that persists. To make matters worse, habitats of Q. spinosa have been degraded and/or destroyed by agricultural pollution and the construction of dams for hydropower projects [20, 22, 24]. This species is currently classified as “vulnerable (A2abc)” on the IUCN Red List [20]. Despite the initiation of captive breeding programmes in the 1980s, low fertilization and hatching rates, diseases, and high overwinter mortality mean that the cultivation of Q. spinosa strongly depends on the supply of tadpoles or adults from wild populations [24]. Thus, other conservation measures to arrest the population decline of this species are urgently required.

Data collection and processing

The study area constituted central and southern China and adjacent countries in Southeast Asia, including Vietnam, Laos, Thailand, Cambodia, the Philippines, and Myanmar (10–35°N, 100–125°E; Fig. 1). We collected a total of 1714 occurrence records of Q. spinosa from field expeditions conducted by Chengdu Institute of Biology (CIB) personnel (2010–2015), supplemented with location data from the literature [24,25,26], the Global Biodiversity Information Facility (GBIF;, and georeferenced specimen records in the Herpetological Museum of the Chengdu Institute of Biology (CIB). To avoid georeferencing errors and over-fitting in the following modelling, we checked all location data in ArcGIS 9.2 (ESRI, Redland, USA) and removed duplicate occurrences at a spatial resolution of 1 × 1-km so that each grid cell had only a single record (in each grid cell, we randomly selected one occurrence and deleted all the other occurrences) [27]. We obtained a total of 136 occurrences composing the final dataset (Fig. 1).

Fig. 1

Study area and occurrence records of Quasipaa spinosa

Based on the habitats of Q. spinosa [20, 22] and our field experiences, we considered five types (i.e., climate, habitat, biogeography, topography, and human impact) of environmental variables, including 22 parameters (Table S1 in Additional File 1 [28,29,30,31,32,33];) that may impact the habitat use and distribution of Q. spinosa in our modelling processes. For several variables of each type, we can use dimension-reduction techniques (e.g., correlation analysis, clustering algorithms) by pruning the variables with possible redundancy [7, 34]. This procedure was to minimize the risk of overfitting species-environment relationships. Specifically, we retained those variables that gave a higher value in the jackknife analysis [35], which consequently yielded the 14 climatic variables and eight other types of variables (Table S1 in Additional File 1). All variables were transformed into 1 × 1-km equal-area grid rasters in ArcGIS 9.2 and projected onto the UTM WGS 1984 projection.

To capture the uncertainties in future climate projections [9] and represent different temperature sensitivities and greenhouse gas emissions, we considered six global climate models (GCMs: NCC-NorEsm1-M, MRI-CGCM3, BCC-CSM1.1, CSIRO_MK_3_6_0, MIROC_ESM_CHEM, and IPSL_CM5A_LR) for four representative greenhouse gas concentration pathways (RCPs: RCP2.6, RCP4.5, RCP6.0, and RCP8.5) from the IPCC 5th Assessment report [1]. We obtained the 14 climatic variables for 2050 and 2070 from the WorldClim 1.4 dataset [36], that is, we collected 24 projections (6 GCMs × 4 RCPs) for each variable at each time slice to illustrate the temperature and precipitation conditions at that time. For the past climate, variable layers were obtained for the time frames of the Mid-Holocene (MidH, ~ 6 kyr BP) [36], the Last Glacial Maximum (LGM, ~ 22 kyr BP) [37], and the LIG [38]. Three GCMs (CCSM4, Mirco-ESM, and MPI-ESM-P) were used for both the MidH and LGM, but only one GCM (CCSM) was available for the LIG.

For each climate change scenario, climatic variables were used to make projections into the past and the future. The non-climatic variables (i.e., AI, aridity index; PET, annual potential evapo-transpiration; Landcover, land cover type; NPP, net primary productivity; Biome, terrestrial ecoregion; HFP, human footprint index) were assumed to be constant over time, as past and future estimates of these factors were, unfortunately, not available. In addition, these variables are likely to be determined not only by climatic drivers but also a wide range of socioeconomic drivers, and any simple estimates or extrapolations could be misleading [cf. 39]. Thus, their constancy is an assumption that presents a conservative prognosis and limits additional uncertainties [cf. 39]. To determine the influence of the non-climatic data, we predicted and analysed the environmental suitability for Q. spinosa by comparing two competing models [39,40,41]: the full model (using all 22 variables) and a climate-only model (using only the 14 bioclimatic variables).

Ensemble distribution modelling and climate change impacts

We predicted the potential distributions of Q. spinosa under current climatic conditions using a maximum entropy approach in MaxEnt v3.3.3 k [42, 43]. MaxEnt has been shown to have good performance and consistently outperforms many other methods, especially with small sample sizes, noisy input data, and correlated parameters [43, 44]. As an inappropriate model complexity (e.g., inappropriate combination of feature classes, level of regularization, or variable selection) or data organization (e.g., occurrence and background localities partition) can typically reduce the accuracy of inferred habitat quality, we selected model settings and conducted model choices using the ENMeval 0.1.0 package in R 4.0.3 software [45,46,47,48]. Following the method outlined in Muscarella et al. [47], five feature class values (FC: linear, quadratic, product, threshold, hinge), eight regularization multiplier values (RM: 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4), five locality partition methods (jackknife, randomkfold, block, checkerboard1, checkerboard2) were calculated; the mean values of the area under the receiver operating characteristics curve using the testing data across all data bins (mean AUCtest) and the mean values of the difference between the training and testing AUCs across all data bins (mean AUCdiff) were computed as the evaluation metrics for the goodness-of-fit and the degree of overfitting of the model. For model selection, we used the Akaike Information Criteria corrected for a small sample size (AICc) to compare different variable combinations [46, 49]. The existence of sampling bias in background points could decrease modelling effectiveness; thus, we used a kernel density estimator (KDE) surface to create 10,000 random background points in the Software for Automated Habitat Modelling (SAHM) [50]. As the models incorporating all environmental factors (i.e., 22 variables for the full model, 14 variables for the climate-only model) and with the conditions FC = linear, RM = 1, and the randomkfold (k = 10) partition method showed the best performances (see Results), 80% of the localities were randomly selected to generate a training set, and the remaining 20% of the localities were withheld to use as an evaluation set. Other model settings were adopted, including the maximum number of iterations (500) and the convergence threshold (10− 5) [43]. We selected the logistic output format with suitability values ranging from 0 to 1 and conducted jackknife procedures to evaluate the relative contribution of each variable [43].

We then projected the contemporary species-environment relationships into different time slices (the past: LIG, LGM and MidH; the future: 2050 and 2070). For each time slice and each environmental variable group, we ran ten cross-validation replications of the model and weighted them by their AUCs to obtain an ensemble distribution prediction [51, 52]. Next, we averaged the outputs for each past (3GCMs each for LGM and MidH) and future time slice (6 GCMs × 4 RCPs) to achieve a suitability map for each timeframe.

To convert continuous outputs into presence-absence maps, we extracted the suitability values of the occurrences of Q. spinosa from the raw modelling results and calculated the mean values as the thresholds (full model: 0.61, climate-only model: 0.58) [53]. That is, areas with predicted suitability values above and below the threshold were considered “present” (suitable) and “absent” (unsuitable), respectively. Considering the home ranges, movements and dispersal abilities of anurans [54], we excluded habitat patches < 10 km2 in size from the subsequent analyses. To quantify the impacts of climate change, we first calculated the total area of suitable habitats, the number of suitable habitat patches, and the area of the largest and smallest habitat patches for each time slice. Then, the latitudes, longitudes, and elevations of the “present” grid cells were extracted, and one-way ANOVAs were conducted to detect their variations among time slices. In addition, to show the effectiveness of the existing protected area network on protecting Q. spinosa habitats, we overlaid the distributions of protected areas ([55]; only national and international protected areas were included as they have relatively reliable management efforts [56, 57]) and suitable habitat patches at present and in the future (for the current time, 2050, and 2070; for the full model and the climate-only model). We calculated the percentages of the area and the percentages of the number of patches of suitable habitats covered by the protected areas for each time slice.

Landscape connectivity assessment

Based on the response curves, which show the relationships among environmental variable gradients and the predicted suitability of species, we ranked each variable into rasters ranging from 1 to 10, with higher values indicating higher habitat quality (Table S2 in Additional File 1). Low-quality grid cells have high movement costs for organisms, whereas good habitats are associated with lower resistance to movement [14, 58]. We then reversed the integer categories of each variable for habitat quality to reclassify the movement cost score of Q. spinosa by grid cell (Table S2 in Additional File 1). To create a movement cost surface for each time period, we weighted the reclassified rasters of environmental variables with training gains with only the variables obtained from the jackknife procedure (Table S2 in Additional File 1), and then averaged the movement cost surfaces from all the GCMs and RCPs (3 GCMs for LGM and MidH, 6 GCMs × 4 RCPs for 2050 and 2070).

Landscape connectivity was assessed using a pairwise cost distance method, the Cost Distance Tool in ArcGIS 9.2 [59]. This method evaluates the best potential dispersal linkage by summing the movement cost of each grid cell between two core habitat patches. The best potential dispersal linkage is the one that has the least resistance (the least cost patch-to-patch path, LCP). Based on the movement cost surfaces and suitable habitats, we created pairwise LCPs among habitat patches for each time slice with the Corridor Design Tool [59]. We accounted for the potential limits to dispersal distances by using four LCP length caps: 20 km, 50 km, 100 km, and no limit. To investigate the dynamics of landscape connectivity as consequences of climate change over time, we compared the numbers, mean lengths, and mean movement costs of the LCPs among the time slices.


The full model with 22 environmental variables (mean ± SE: AUCtest = 0.948 ± 0.033, AUCdiff = 0.023 ± 0.002, AICc = 712.328, AICw = 0.412) and the climate-only model with 14 variables (AUCtest = 0.922 ± 0.044, AUCdiff = 0.029 ± 0.003, AICc = 3005.074, AICw = 0.407) both had great performances (FC = linear, RM = 1, randomkfold partition method: k = 10). Precipitation in the dry season (Bio14 & Bio17), annual precipitation (Bio12) and precipitation seasonality (Bio15) were the most important contributing parameters for model development, and the aridity index (AI) and ecoregion type (Biome) had obvious impacts on predicted suitability in the full model (Fig. S1 in Additional File 2).

Over time, the suitable habitats of Q. spinosa moved noticeably as a consequence of climatic variations among the different time slices. Relatively wider ranges of suitable habitats were shown in the MidH and current time periods than in the other time slices, with regions of higher suitability concentrated in the central-south and south-eastern parts of China and in the eastern part of the Indo-Chinese Peninsula. More restricted highly suitable ranges were revealed in the Quaternary LIG and LGM and in the future (Fig. 2). The total area, number, and mean size of habitat patches decreased from the LIG to the LGM, followed by obvious increases in the MidH and the present before declining under projected future conditions (i.e., in 2050 and 2070; Table 1). Compared with the full model, the outputs projected by the climate-only model were similar but showed larger areas of high suitability. Overall, we found a significant latitudinal increase in the area of suitable habitat over time (ANOVA: full model, P = 0.01; climate-only model, P = 0.02; Fig. 3), indicating northward range expansion across the climate change scenarios. No significant shift along either the elevational or the longitudinal gradient was detected (all P > 0.05; except for a slightly upward trend for the full model, P = 0.04; Fig. 3). Our results showed that, based on the full model (the climate-only model), the existing protected area network covered only 12.1% (10.3%) of the total Q. spinosa habitat area, and only 12.2% (14.3%) of the suitable habitat patches overlapped or partly overlapped with protected areas at the present (Fig. 4). For 2050, the habitat area and patch number percentages were 29.2 and 12.5% based on the full model and 22.7 and 21.6% based on the climate-only model, respectively (Fig. 4). For 2070, the habitat area and patch number percentages were 31.8 and 15.6% based on the full model and 30.1 and 30.0% based on the climate-only model, respectively (Fig. 4).

Fig. 2

Habitat suitability predictions for Quasipaa spinosa. Top two rows: outputs are from the full model (all 22 environmental variables are included). Bottom two rows: outputs are based on the climate-only model (only the 14 bioclimatic variables are included)

Table 1 Total suitable habitat areas, numbers of habitat patches, and mean areas of habitat patches for Quasipaa spinosa over time
Fig. 3

Distribution shifts of Quasipaa spinosa at the elevational, latitudinal, and longitudinal dimensions under climate change (over the time periods of the LIG, the LGM, the MidH, Current, 2050, and 2070). abc: full model (all 22 environmental variables are included), df = 5, 218,228, α = 0.05; def: climate-only model (only the 14 bioclimatic variables are included), df = 5, 589,914, α = 0.05. Black squares show the means, solid lines show the medians, the edges of boxes are the quartiles, the crosses are the lower and upper adjacent limits, and the whiskers are the standard deviations

Fig. 4

Overlaps between the existing protected area network and the patches of suitable habitats for Quasipaa spinosa

Our results showed high landscape connectivity in the MidH and current time frames compared to those of the past (the LIG and LGM) and future (2050 and 2070; Fig. 5). The LCP values among habitat patches generally decreased across the ice ages (from the LIG to LGM) and increased after glaciations (from the LGM to the MidH and the current period) and then decreased again by more than half by 2050 and by nearly 75% by 2070 (Table 2). Accordingly, the mean LCP length and movement cost among patches revealed an increase-decrease-increase pattern over time (Table 2). This finding was robust to assumptions about maximum LCP lengths.

Fig. 5

Suitable habitat patches and pairwise patch-to-patch least cost paths (LCPs) of Quasipaa spinosa. Top two rows: outputs are from the full model (all 22 environmental variables are included). Bottom two rows: outputs are based on the climate-only model (only the 14 bioclimatic variables are included)

Table 2 The numbers, mean lengths, and mean movement costs of patch-to-patch LCPs (least cost paths) for Quasipaa spinosa over time under different LCP length limits based on the full and climate-only models


Our results predicted how the locations and extents of suitable habitats for Q. spinosa varied and will vary as functions of past and future climate changes. The most striking finding is that the inferred future loss of habitat connectivity could lead to a decrease in individual dispersal, increasing the extinction risk of this threatened frog species. In summary, this study presents robust, spatially explicit predictions of historical and future suitable habitats for Q. spinosa.

Global paleoclimate fluctuations constitute dramatic cycles of environmental change on Earth. These fluctuations have profoundly shifted the distributions of many plant and animal species [60, 61]. Characterized by at least 30 intermittent glacial-interglacial cycles, the Pleistocene climatic oscillations are believed to have caused many iterations of massive species extirpations over large portions of their ranges, dispersal events to new locations, the trapping of species in refugia, and postglacial expansions [60, 61]. In East Asia, the climate is characterized as a unique monsoon system with dry winters and abundant precipitation in summers thanks to the uplift of the Qinghai-Tibetan Plateau. From the LIG to the MidH, East Asia experienced expanded ice caps and extensive glacier valley systems rather than extensive ice sheets, especially during the LGM [60, 62,63,64]. Such relatively moderate conditions can be conducive for species to retain hospitable habitats or to make relatively easy migrations to newly suitable areas [65, 66]. Previous studies suggested that during cold periods, some amphibian species in southern and eastern China might disperse through streams that provide connectivity to refugia, including some mountainous regions, which can be topographically heterogeneous with abundant suitable habitats with relatively stable microclimates [67, 68]. We found that there was a north-easterly shift of suitable Q. spinosa habitats from the LIG to the LGM as a result of warming conditions and the angle of the coastline with the South China Sea. This shift may have moved frog populations closer to mountains that could function as refugia over geologic time, with the frogs using the abundant streams and river branches within the water systems (including the Yangtze River network and the Pearl River network) in South China as their dispersal passages [69]. The suitable habitats generally moved from the regions around the southern coasts of China and the eastern countries of south-eastern Asia to some mountainous areas in eastern China, including Huangshan, Dabieshan-Wuyishan-Luoxiaoshan, Nanling, Wushan-Xuefengshan, and the Qinling-Dabashan Mountains (Fig. 2). As the temperatures rose after glaciation (in the MidH and the current period), suitable habitats were greatly enlarged and expanded into surrounding regions. In the coming decades, the continuously changing climate will contract and restrict the potential distribution of the frogs around the hilly areas of eastern China (Fig. 2).

We found that precipitation in the dry season, together with temperature variables and rainfall abundance in eastern Asia, acted as the main drivers of the dynamics of suitable habitats (Fig. S1 in Additional File 2). We found that not only did the cold and dry climatic conditions that occurred during the ice ages pose challenges for the survival of Q. spinosa, but also that the “overwarm” climate in the future could be a threat to these environmentally sensitive and non-heat-tolerant frogs. This phenomenon has also been found in other amphibians [61, 67, 69]. The dual climate challenges mean that the present conditions are optimal for Q. spinosa when compared with all the time slices considered. Our findings suggest that hilly areas around the lower reaches of the Yangtze River, eastern parts of the Yunnan-Guizhou Plateau (also named the Yunnan-Kweichow Plateau), the south-eastern mountains in China, and the Sino-Vietnamese transboundary areas in northern Vietnam (Fig. 2) are high-priority areas for protection. Therefore, on-going and future protection actions should emphasize landscape connectivity concurrently with habitat area and the number of habitat patches.

Our results showed that landscape connectivity was low during the cold periods in the geologic past (from the LIG to the LGM), followed by dramatic improvements post-glaciation. However, we also found that connectivity is likely to decline again in the face of additional warming. Considering the dispersal capacity of anurans [54], our projections of connectivity based on the assumption of a 20-km path-length limit could be the closest to reality and could provide valuable references to conservation planning for both the present and the future [70]. Regardless of the exact maximum length of dispersal, frogs may not use these paths due to the limited mobility and strict physiological constraints of anurans. As a result, additional research is needed to more fully inform the population management of Q. spinosa [15, 71]. Even assuming we correctly identified inter-patch passageways with the highest potential and permeability for individuals to disperse successfully, constraints to animal movements are complex, with a number of stochastic ecological and landscape factors that may have important effects [15, 19]. Hence, other possible paths that could promote individual exchanges among populations or patches should also be seriously considered in protection programmes. Small isolated fragments of suitable habitats, which serve as “stepping stones” between large patches and dispersal stopovers for frogs, cannot be ignored [cf. 71]. In addition, to reduce possible edge effects and anthropogenic influences, buffering the LCPs into width-sufficient patch-to-patch corridors based on the species’ home range requirements could be a critical step, and the functions or effectiveness of these buffers should be tested with population monitoring data and empirical gene flow studies among patches [14, 15].

Given the threats that already exist [20, 22, 24] and those that will most likely emerge under climate change and due to other anthropogenic impacts in the future, we suggest that appropriate protection strategies and strengthened conservation efforts should be developed immediately and be modified as habitat degradation and fragmentation continue. Specifically, based on the results of this study, we propose the following conservation recommendations for Q. spinosa. As the first step towards efficient protection and population sustainability actions [72], several representative areas and critical refugia (i.e., habitats, patches, or local populations) identified here should be focused on as key suitable habitats, sources of genetic diversity, and conservation priorities. These areas include the hilly and mountainous regions in southern, south-eastern, and south-western China and the transboundary region between China and Vietnam, which requires transboundary cooperation for species protection. Second, the LCPs in our map of the present and future distributions of Q. spinosa habitats should be emphasized in conservation planning. Actions that enhance inter-patch dispersal and the exchange of frog individuals between patches should be considered. For example, preserving LCPs as sufficiently wide ecological connections and preserving habitat stability within corridors should be pursued as conservation priorities. Because of future reductions in landscape connectivity, long-term resources should be allocated to monitor and analyse habitat quality, individual dispersal, and gene flow among patches to assess conservation effectiveness. Third, the networks of protected areas within the extant range of Q. spinosa should be reassessed in light of our findings. As the coverage ratios of habitat areas and patch numbers in the protected area network are relatively low (< 30% or approximately 30%), selecting new protected areas and adjusting existing protected areas may be urgent in order to ensure that sufficient habitat areas are protected for the species, especially in the eastern part of its distribution region. Finally, to help counteract population decline, environmental education in local communities, the refinement of protection laws and management systems, and the restriction of human activities may also be necessary.

We recognize that some caveats exist in this study. Our results only predict and quantify the dynamics of suitable habitats and landscape connectivity at the landscape scale but do not consider the potential impacts of climate change at a finer resolution. Second, more efforts, such as analysing gene flow, monitoring individuals’ dispersals among populations, and comparing different landscape models, are needed to empirically validate some of the assumptions made in our model. The third caveat is that the influences of some other ecological factors, such as predation risk and local conservation efforts, should be considered to explain local population dynamics [14, 73]. Furthermore, if we took into account the human impacts on and microhabitat requirements of Q. spinosa when modelling habitat suitability, the potentially suitable habitats would be more fragmented (Fig. 2; Table 1). An additional problem is the assumption of the constancy of the non-climatic factors (i.e., AI, PET, Landcover, NPP, Biome, HFP) over time; that is, we use current data as conservative predictions of these variables for the future and past. This could lead to a caveat in our modelling results, especially those from the full models (some of the non-climatic variables showed relatively high contributions to habitat suitability). However, we considered the method used in this study to be a reasonable way to restrict additional uncertainty, as these variables could be influenced by a variety of factors (not only climatic but also anthropogenic-mediated factors, etc.) and are thus difficult to estimate, and timely analyses need to be conducted when updated data are available [37]. Finally, to better capture the consequences of future climate change, i.e., range shifts, it would be necessary to quantify genetic diversity and associated ecological drivers of geographic populations [74]. Despite these caveats, by exploring shifts in the suitable habitat and landscape connectivity of an anuran species that occur as a consequence of climate change, our findings can serve as a template for conservation planning for a variety of amphibians that face the interlocking threats of habitat fragmentation and climate change.


Our results reveal that suitable habitats of Q. spinosa obviously declined and shifted northwardly during the Last Glacial Maximum and expanded after the ice age. In the future, suitable habitats are projected to decrease again due to “overwarm” climatic conditions. The dynamics of suitable habitats also lead to habitat connectivity changes in this species during post-glaciation periods. In the face of future climate change, an increase in habitat fragmentation and a decrease in habitat connectivity are predicted. We propose that the mountainous areas in the southern China and Sino-Vietnamese transboundary regions can be considered conservation priorities for this frog species. Given that amphibians are highly sensitive to environmental changes, this study highlights that the dynamics of environmental suitability and habitat connectivity should be considered to guide large-scale conservation management measures for endangered amphibian species under climate change.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.


  1. 1.

    IPCC. Climate change 2014: impacts, adaptation, and vulnerability. Part a: global and sectoral aspects. Contribution of working group II to the fifth assessment report of the intergovernmental panel on climate change. Cambridge: Cambridge University Press; 2014.

    Google Scholar 

  2. 2.

    Hof C, Araujo MB, Jetz W, Rahbek C. Additive threats from pathogens, climate and land-use change for global amphibian diversity. Nature. 2011;480(7378):516–9.

    CAS  Article  PubMed  Google Scholar 

  3. 3.

    Thomas CD, Cameron A, Green RE, Bakkenes M, Beaumont LJ, Collingham YC, Erasmus BFN, de Siqueira MF, Grainger A, Hannah L, Hughes L, Huntley B, van Jaarsveld AS, Midgley GF, Miles L, Ortega-Huerta MA, Townsend Peterson A, Phillips OL, Williams SE. Extinction risk from climate change. Nature. 2004;427(6970):145–8.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Hua F, Hu J, Liu Y, Giam X, Lee TM, Luo H, Wu J, Liang Q, Zhao J, Long X, Pang H, Wang B, Liang W, Zhang Z, Gao X, Zhu J. Community-wide changes in intertaxonomic temporal co-occurrence resulting from phenological shifts. Glob Chang Biol. 2016;22(5):1746–54.

    Article  PubMed  Google Scholar 

  5. 5.

    Garcia RA, Cabeza M, Rahbek C, Araújo MB. Multiple dimensions of climate change and their implications for biodiversity. Science. 2014;344(6183):1247579.

    CAS  Article  PubMed  Google Scholar 

  6. 6.

    Scheffers BR, De Meester L, Bridge TCL, Hoffmann AA, Pandolfi JM, Corlett RT, et al. The broad footprint of climate change from genes to biomes to people. Science. 2016;354:aaf7671.

    Article  Google Scholar 

  7. 7.

    Hu J, Huang Y, Jiang J, Guisan A. Genetic diversity in frogs linked to past and future climate changes on the roof of the world. J Anim Ecol. 2019;88(6):953–63.

    Article  PubMed  Google Scholar 

  8. 8.

    Yang C, Tang S, Luo Z. Distribution changes of Chinese skink (Eumeces chinensis) in China: the impacts of global climate change. Asian Herpetol Res. 2020;11:132–8.

    Google Scholar 

  9. 9.

    Collins M, Chandler RE, Cox PM, Huthnance JM, Rougier J, Stephenson DB. Quantifying future climate change. Nat Clim Chang. 2012;2(6):403–9.

    Article  Google Scholar 

  10. 10.

    Urban MC. Accelerating extinction risk from climate change. Science. 2015;348(6234):571–3.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Piaggio AJ, Segelbacher G, Seddon PJ, Alphey L, Bennett EL, Carlson RH, et al. Is it time for synthetic biodiversity conservation? Trends Ecol Evol. 2016;32:97–107.

    Article  Google Scholar 

  12. 12.

    Dawson TP, Jackson ST, House JI, Prentice IC, Mace GM. Beyond predictions: biodiversity conservation in a changing climate. Science. 2011;332(6025):53–8.

    CAS  Article  PubMed  Google Scholar 

  13. 13.

    Fahrig L. Effects of habitat fragmentation on biodiversity. Annu Rev Ecol Evol S. 2003;34(1):487–515.

    Article  Google Scholar 

  14. 14.

    Sawyer SC, Epps CW, Brashares JS. Placing linkages among fragmented habitats: do least-cost models reflect how animals use landscapes? J Appl Ecol. 2011;48(3):668–78.

    Article  Google Scholar 

  15. 15.

    Crooks KR, Sanjayan M. Connectivity conservation. Cambridge: Cambridge University Press; 2006.

    Book  Google Scholar 

  16. 16.

    Olds AD, Connolly RM, Pitt KA, Maxwell PS. Habitat connectivity improves reserve performance. Conserv Lett. 2011;5:56–63.

    Article  Google Scholar 

  17. 17.

    McGuire JL, Lawler JJ, McRae BH, Nuñez TA, Theobald DM. Achieving climate connectivity in a fragmented landscape. P Natl Acad Sci USA. 2016;113(26):7195–200.

    CAS  Article  Google Scholar 

  18. 18.

    Wake DB, Vredenburg VT. Are we in the midst of the sixth mass extinction? A view from the world of amphibians. P Natl Acad Sci USA. 2008;105(Supplement 1):11466–73.

    Article  Google Scholar 

  19. 19.

    Li Y, Cohen JM, Rohr JR. Review and synthesis of the effects of climate change on amphibians. Integr Zool. 2013;8(2):145–61.

    Article  PubMed  Google Scholar 

  20. 20.

    Lau MWN, Geng B, Gu H, Dijk PPv, Bain R. Quasipaa spinosa. The IUCN Red List of Threatened Species 2004: e.T58439A11781309. Accessed 25 May 2020.

  21. 21.

    Hu J, Li C, Xie F, Jiang J. Endemic amphibians and their distribution in China. Asian Herpetol Res. 2012;3:163–71.

    Article  Google Scholar 

  22. 22.

    Fei L, Hu S, Ye C, Huang Y. Fauna Sinica. Amphibia. Vol. 3. Anura Ranidae. Beijing: Science Press; 2009.

    Google Scholar 

  23. 23.

    Hu J, Xie F, Li C, Jiang J. Elevational patterns of species richness, range and body size for spiny frogs. PLoS One. 2011;6(5):e19817.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  24. 24.

    Chan HK, Shoemaker KT, Karraker NE. Demography of Quasipaa frogs in China reveals high vulnerability to widespread harvest pressure. Biol Conserv. 2014;170:3–9.

    Article  Google Scholar 

  25. 25.

    Yu D, Zheng R, Lu Q, Yang G, Yao F, Zhang Y. Genetic diversity and population structure for the conservation of giant spiny frog (Quasipaa spinosa) using microsatellite loci and mitochondrial DNA. Asian Herpetol Res. 2016;7:75–86.

    Google Scholar 

  26. 26.

    Van Ngo B, Lee, Ngo CD. Variation in dietary composition of granular spiny frogs (Quasipaa verrucospinosa) in Central Vietnam. Herpetol J. 2014;24:245–53.

    Google Scholar 

  27. 27.

    Hu J, Jiang J. Inferring ecological explanations for biogeographic boundaries of parapatric Asian mountain frogs. BMC Ecol. 2018;18(1):3.

    Article  PubMed  PubMed Central  Google Scholar 

  28. 28.

    The CGIAR consortium for spatial information (CGIAR-CSI). CGIAR International Research Centers. 1999. Accessed 25 May 2020.

  29. 29.

    Olson DM, Dinerstein E, Wikramanayake ED, Burgess ND, Powell GVN, Underwood EC, D'amico JA, Itoua I, Strand HE, Morrison JC, Loucks CJ, Allnutt TF, Ricketts TH, Kura Y, Lamoreux JF, Wettengel WW, Hedao P, Kassem KR. Terrestrial ecoregions of the world: a new map of life on earth: a new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity. BioScience. 2001;51(11):933–8.[0933:TEOTWA]2.0.CO;2.

    Article  Google Scholar 

  30. 30.

    Thematic Database of Human-earth System. Institute of Geographic Sciences Natural Resources Research, Chinese Academy of Sciences. 2002. Accessed 25 May 2020.

  31. 31.

    Global and Cover 2000 Database. European Commission, Joint Research Centre, Institute for Environment and Sustainability 2003. Accessed 25 May 2020.

  32. 32.

    Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37(12):4302–15.

    Article  Google Scholar 

  33. 33.

    Gridded Population of the World, Version 4 (GPWv4): Land and Water Area. New York: NASA Socioeconomic Data and Applications Center (SEDAC) Palisades. Center for International Earth Science Information Network-CIESIN-Columbia University. 2016. Accessed 25 May 2020.

  34. 34.

    Synes NW, Osborne PE. Choice of predictor variables as a source of uncertainty in continental-scale species distribution modelling under climate change. Glob Ecol Biogeogr. 2011;20(6):904–14.

    Article  Google Scholar 

  35. 35.

    Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Model. 2006;190(3-4):231–59.

    Article  Google Scholar 

  36. 36.

    Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25(15):1965–78.

    Article  Google Scholar 

  37. 37.

    Collins WD, Bitz CM, Blackmon ML, Bonan GB, Bretherton CS, Carton JA, Chang P, Doney SC, Hack JJ, Henderson TB, Kiehl JT, Large WG, McKenna DS, Santer BD, Smith RD. The community climate system model version 3 (CCSM3). J Clim. 2006;19(11):2122–43.

    Article  Google Scholar 

  38. 38.

    Otto-Bliesner BL, Marshall SJ, Overpeck JT, Miller GH, Hu A. CAPE last interglacial project members. Simulating Arctic climate warmth and icefield retreat in the last interglaciation. Science. 2006;311(5768):1751–3.

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    Thuiller W, Broennimann O, Hughes G, Alkemade JRM, Midgley GF, Corsi F. Vulnerability of African mammals to anthropogenic climate change under conservative land transformation assumptions. Glob Chang Biol. 2006;12(3):424–40.

    Article  Google Scholar 

  40. 40.

    Hu J, Jiang Z. Predicting the potential distribution of the endangered Przewalski's gazelle. J Zool. 2010;282(1):54–63.

    Article  Google Scholar 

  41. 41.

    Luo Z, Zhou S, Yu W, Yu H, Yang J, Tian Y, Zhao M, Wu H. Impacts of climate change on the distribution of Sichuan snub-nosed monkeys (Rhinopithecus roxellana) in Shennongjia area. China Am J Primatol. 2015;77(2):135–51.

    Article  PubMed  Google Scholar 

  42. 42.

    Phillips SJ, Dudík M, Schapire RE. Maxent software for modeling species niches and distributions (Version 3.4.1). 2017.

    Google Scholar 

  43. 43.

    Phillips SJ, Dudík M. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography. 2008;31(2):161–75.

    Article  Google Scholar 

  44. 44.

    Elith J, Graham CH, Anderson RP, Dudík M, Ferrier S, Guisan A, et al. Novel methods improve prediction of species' distributions from occurrence data. Ecography. 2006;29(2):129–51.

    Article  Google Scholar 

  45. 45.

    Peterson AT, Papeş M, Soberón J. Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecol Model. 2008;213(1):63–72.

    Article  Google Scholar 

  46. 46.

    Warren DL, Seifert SN. Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecol Appl. 2011;21(2):335–42.

    Article  PubMed  Google Scholar 

  47. 47.

    Muscarella R, Galante PJ, Soley-Guardia M, Boria RA, Kass JM, Uriarte M, Anderson RP. ENMeval: an R package for conducting spatially independent evaluations and estimating optimal model complexity for MAXENT ecological niche models. Methods Ecol Evol. 2014;5(11):1198–205.

    Article  Google Scholar 

  48. 48.

    The R Core Team. R: A language and environment for statistical computing (Vienna, Austria). 2020.

    Google Scholar 

  49. 49.

    Burnham KP, Anderson DR. Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods Res. 2016;33:261–304.

    Article  Google Scholar 

  50. 50.

    Phillips SJ, Dudík M, Elith J, Graham CH, Lehmann A, Leathwick J, Ferrier S. Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecol Appl. 2009;19(1):181–97.

    Article  PubMed  Google Scholar 

  51. 51.

    Marmion M, Parviainen M, Luoto M, Heikkinen RK, Thuiller W. Evaluation of consensus methods in predictive species distribution modelling. Divers Distrib. 2009;15(1):59–69.

    Article  Google Scholar 

  52. 52.

    Araújo MB, New M. Ensemble forecasting of species distributions. Trends Ecol Evol. 2007;22(1):42–7.

    Article  PubMed  Google Scholar 

  53. 53.

    Hu J, Liu Y. Unveiling the conservation biogeography of a data-deficient endangered bird species under climate change. PLoS One. 2014;9(1):e84529.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Wells KD. The ecology and behavior of amphibians. Chicago: The University of Chicago Press; 2007.

    Book  Google Scholar 

  55. 55.

    UNEP-WCMC (United Nations Environment Programme-World Monitoring Center). World Database on Protected Areas (WDPA). UNEP-WCMS, New York, USA. 2011. Accessed 07 January 2021.

  56. 56.

    Luo Z, Wei S, Zhang W, Zhao M, Wu H. Amphibian biodiversity congruence and conservation priorities in China: integrating species richness, endemism, and threat patterns. Biol Conserv. 2015;191:650–8.

    Article  Google Scholar 

  57. 57.

    Maiorano L, Falcucci A, Boitani L. Size-dependent resistance of protected areas to land-use change. P Roy Soc B-Biol Sci. 2008;275:1297–304.

    Google Scholar 

  58. 58.

    Saura S, Pascual-Hortal L. A new habitat availability index to integrate connectivity in landscape conservation planning: comparison with existing indices and application to a case study. Landscape Urban Plan. 2007;83(2-3):91–103.

    Article  Google Scholar 

  59. 59.

    Majka D, Jenness J, Beier P. Corridor Designer: ArcGIS Tools for Designing and Evaluating Corridors. In: Northern Arizona University: Environmental Research, Development and Education for the New Economy (ERDENE); 2007.

    Google Scholar 

  60. 60.

    Hewitt G. The genetic legacy of the quaternary ice ages. Nature. 2000;405(6789):907–13.

    CAS  Article  PubMed  Google Scholar 

  61. 61.

    Qiu Y, Fu C, Comes HP. Plant molecular phylogeography in China and adjacent regions: tracing the genetic imprints of quaternary climate and environmental change in the world’s most diverse temperate flora. Mol Phylogenet Evol. 2011;59(1):225–44.

    Article  PubMed  Google Scholar 

  62. 62.

    An Z, Zhang P, Wang E, Wang S, Qiang X, Li L, et al. Changes of the monsoon-arid environment in China and growth of the Tibetan plateau since the Miocene. Quat Sci. 2006;26:678–93.

    Google Scholar 

  63. 63.

    Shi Y, Ren B, Wang J, Derbyshire E. Quaternary glaciation in China. Quat Sci Rev. 1986;5:503–7.

    Article  Google Scholar 

  64. 64.

    Owen LA, Caffee MW, Finkel RC, Seong YB. Quaternary glaciation of the Himalayan-Tibetan orogen. J Quater Sci. 2008;23(6-7):513–31.

    Article  Google Scholar 

  65. 65.

    Song G, Qu Y, Yin Z, Li S, Liu N, Lei F. Phylogeography of the Alcippe morrisonia (Aves: Timaliidae): long population history beyond late Pleistocene glaciations. BMC Evol Biol. 2009;9(1):143.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  66. 66.

    Harrison SP, Yu G, Takahara H, Prentice IC. Diversity of temperate plants in East Asia. Nature. 2001;413(6852):129–30.

    CAS  Article  PubMed  Google Scholar 

  67. 67.

    Li J, Wei S, Hu M, Luo Z, Zhao M, Wu H. Reflection of paleoclimate oscillations and tectonic events in the phylogeography of moustache toads in southern China. J Zool. 2018;305(1):17–26.

    Article  Google Scholar 

  68. 68.

    Zhang H, Yan J, Zhang G, Zhou K. Phylogeography and demographic history of Chinese black-spotted frog populations (Pelophylax nigromaculata): evidence for independent refugia expansion and secondary contact. BMC Evol Biol. 2008;8(1):21.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  69. 69.

    Wu Y, Wang Y, Jiang K, Hanken J. Significance of pre-quaternary climate change for montane species diversity: insights from Asian salamanders (Salamandridae: Pachytriton). Mol Phylogenet Evol. 2013;66(1):380–90.

    Article  PubMed  Google Scholar 

  70. 70.

    Binzenhöfer B, Schröder B, Strauss B, Biedermann R, Settele J. Habitat models and habitat connectivity analysis for butterflies and burnet moths-the example of Zygaena carniolica and Coenonympha arcania. Biol Conserv. 2005;126(2):247–59.

    Article  Google Scholar 

  71. 71.

    Luo Z, Yu H, Pu Y, Yang J, Mei H, Wang D, Zhu Z, Zhao M, Wu H. Assessment of habitat fragmentation and corridors for an isolated subspecies of the Sichuan golden snub-nosed monkey, Rhinopithecus roxellana hubeiensis. Inter J Primatol. 2016;37(3):438–59.

    CAS  Article  Google Scholar 

  72. 72.

    Araújo MB, Williams PH. Selecting areas for species persistence using occurrence data. Biol Conserv. 2000;96(3):331–45.

    Article  Google Scholar 

  73. 73.

    Beier P, Majka DR, Spencer WD. Forks in the road: choices in procedures for designing wildland linkages. Conserv Biol. 2008;22(4):836–51.

    Article  PubMed  Google Scholar 

  74. 74.

    Liu Y, Webber S, Bowgen K, Schmaltz L, Bradley K, Halvarsson P, Abdelgadir M, Griesser M. Environmental factors influence both abundance and genetic diversity in a widespread bird species. Ecol Evol. 2013;3(14):4683–95.

    Article  PubMed  PubMed Central  Google Scholar 

Download references


We thank Ka Wah Leung from Chengdu Institute of Biology, CAS, allowed us to use his photograph of Quasipaa spinosa in Fig. 1.


This study was supported by the National Natural Science Foundation of China (No. 32071544, 31770427, 32071497, 31770568), Technology Basic Work Programme (2013FY111500), and the ‘Light of West China’ Programme of the Chinese Academy of Sciences. The funding bodies had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information




Z.L. and J.H. conceived and designed the study; Z.L. and J.H. prepared the data with assistance from X.W., S.Y., X. C and Y.L; Z.L. and J.H. analyzed the data; Z.L. and J.H. wrote the article with assistance from X.W., S.Y., X. C and Y.L; and all authors read and approved the final manuscript.

Corresponding author

Correspondence to Junhua Hu.

Ethics declarations

Ethics approval and consent to participate

No animal was killed for the purpose of this study. All specimens collecting and processing were in accordance with the Law of the People’s Republic of China on the Protection of Wildlife and approved by the Animal Care Committee of CIB, CAS.

Consent for publication

Not Applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Environmental variables used in potential distribution modelling for Quasipaa spinosa. Table S2. Habitat quality score and movement cost score of each environmental variable and their weights in movement cost surface for Quasipaa spinosa.

Additional file 2: Figure S1.

Jackknife analyses on the contributions of environmental variables when modelling potential distributions for Quasipaa spinosa: (a) full model (all 22 environmental variables are included), (b) climatic-only model (only 14 bioclimatic variables are included).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Luo, Z., Wang, X., Yang, S. et al. Combining the responses of habitat suitability and connectivity to climate change for an East Asian endemic frog. Front Zool 18, 14 (2021).

Download citation


  • Climatic condition
  • Giant spiny frog
  • Habitat connectivity
  • Dispersal ability
  • Suitable habitat