Identifying core habitats and corridors for giant pandas by combining multiscale random forest and connectivity analysis

Abstract Habitat loss and fragmentation are widely acknowledged as the main driver of the decline of giant panda populations. The Chinese government has made great efforts to protect this charming species and has made remarkable achievements, such as population growth and habitat expansion. However, habitat fragmentation has not been reversed. Protecting giant pandas in a large spatial extent needs to identify core habitat patches and corridors connecting them. This study used an equal‐sampling multiscale random forest habitat model to predict a habitat suitability map for the giant panda. Then, we applied the resistant kernel method and factorial least‐cost path analysis to identify core habitats connected by panda dispersal and corridors among panda occurrences, respectively. Finally, we evaluated the effectiveness of current protected areas in representing core habitats and corridors. Our results showed high scale dependence of giant panda habitat selection. Giant pandas strongly respond to bamboo percentage and elevation at a relatively fine scale (1 km), whereas they respond to anthropogenic factors at a coarse scale (≥2 km). Dispersal ability has significant effects on core habitats extent and population fragmentation evaluation. Under medium and high dispersal ability scenarios (12,000 and 20,000 cost units), most giant panda habitats in the Qionglai mountain are predicted to be well connected by dispersal. The proportion of core habitats covered by protected areas varied between 38% and 43% under different dispersal ability scenarios, highlighting significant gaps in the protected area network. Similarly, only 43% of corridors that connect giant panda occurrences were protected. Our results can provide crucial information for conservation managers to develop wise strategies to safeguard the long‐term viability of the giant panda population.

random forest habitat model to predict a habitat suitability map for the giant panda.
Then, we applied the resistant kernel method and factorial least-cost path analysis to identify core habitats connected by panda dispersal and corridors among panda occurrences, respectively. Finally, we evaluated the effectiveness of current protected areas in representing core habitats and corridors. Our results showed high scale dependence of giant panda habitat selection. Giant pandas strongly respond to bamboo percentage and elevation at a relatively fine scale (1 km), whereas they respond to anthropogenic factors at a coarse scale (≥2 km). Dispersal ability has significant effects on core habitats extent and population fragmentation evaluation. Under medium and high dispersal ability scenarios (12,000 and 20,000 cost units), most giant panda habitats in the Qionglai mountain are predicted to be well connected by dispersal. The proportion of core habitats covered by protected areas varied between 38% and 43% under different dispersal ability scenarios, highlighting significant gaps in the protected area network. Similarly, only 43% of corridors that connect giant panda occurrences were protected. Our results can provide crucial information for conservation managers to develop wise strategies to safeguard the long-term viability of the giant panda population.

K E Y W O R D S
factorial least-cost path, multiscale habitat selection, Qionglai mountain, random forest, resistant kernel, UNICOR

T A X O N O M Y C L A S S I F I C A T I O N
Conservation ecology; Landscape ecology; Landscape planning; Spatial ecology

| INTRODUC TI ON
The giant panda (Ailuropoda melanoleuca) is a rare and protected wildlife endemic to China and a flagship species of biodiversity conservation globally (Swaisgood et al., 2010). The giant panda once roamed throughout most of the lowlands of eastern and southern China, northern Vietnam, and northern Myanmar (Pan et al., 2014).
Due to human activity, climate change, and natural disasters, the habitat of giant pandas has been continuously lost and fragmented.
Currently, giant pandas are only distributed in part of mountainous areas in Sichuan, Shaanxi, and Gansu (State Forestry Administration, 2015). The Chinese government has made great efforts to protect this charming species by establishing 67 reserves, Grain to Green Project, and Natural Forest Protection Program (Wei et al., 2015).
These conservation efforts have led to some achievements, including population growth and habitat expansion (State Forestry Administration, 2015). Results of the Fourth National Giant Panda Survey (hereinafter the fourth survey) revealed that there is an estimated population size of 1,864 individuals in the wild and showed a 16.8% population increase compared to the third survey, which was conducted from 1998 to (State Forestry Administration, 2015. Based on the observed population increase, the International Union for Conservation of Nature (IUCN) changed the status of the giant panda from "endangered" to "vulnerable" (Swaisgood et al., 2016).
However, the panda habitat is becoming increasingly fragmented (Xu et al., 2017), and small populations will face high extinction risks (Kong et al., 2021). According to the fourth survey, the panda's range is estimated to be subdivided into about 33 subpopulations separated by mountain ranges, rivers, roads, forest clearings, and human settlements (State Forestry Administration, 2015). Despite the enormous efforts that have been put in panda conservation, approximately 46% of the habitat (33% of the panda population) remains unprotected (State Forestry Administration, 2015). It is urgently needed for a knowledge-based metapopulation management strategy for the long-term viability of giant panda subpopulations . Establishing and protecting core habitat patches and the connectivity networks that connect them is one of the ways to ensure the long-term survival of large terrestrial mammals at a regional scale (Kaszta et al., 2020). Additionally, it is unlikely to protect all the landscape with limited financial resources, large human populations, and complicated land ownership, making it crucial to identify core habitats and corridors. Several previous studies have identified habitat connectivity for giant pandas in different mountain ranges using least-cost path analysis or circuit theory Qi et al., 2012;Wang et al., 2021). Researchers usually first mapped habitat patches based on habitat suitability and then simulated corridors among patches in these studies. There were some shortcomings in them. For example, it is a simplification using habitat patches as source points rather than species occurrences, besides few studies took giant panda dispersal ability into account, though these two aspects are foundations for reliable prediction of corridor networks . Furthermore, it was hard to prioritize corridors in these assessments given corridor construction and restoration are projects that consume both huge manpower and money (Kang & Li, 2016). Therefore, there is a need to apply more comprehensive approaches to map corridors.
Landscape resistance is a crucial component in connectivity modeling. It is challenging to quantify resistance to movement in a large extent because movement data are usually unavailable. Given the lack of movement or genetic data, habitat suitability is frequently used as a proxy to reflect landscape resistance (Zeller et al., 2012).
Therefore, habitat suitability models may have essential impacts on resistance estimation. Research on the habitat of giant pandas facilitates our understanding of the resource needs and ongoing threats and is also a necessary basis for conservation decision-making (Hull et al., 2014). The habitat selection of animals is multidimensional, and the response of animals to different environmental factors often occurs at multiple hierarchical levels and over a range of spatial and temporal scales Wiens, 1989). When describing the relationship between species and habitat, it is necessary to determine the suite of the covariates relevant to habitat selection by the species and determine the scale of interaction between species and habitat (Graf et al., 2005). Incorrect insight into the nature and significance of relationships between species responses and environmental variables may result from ignoring scale in habitat modeling . As an iconic species of global biodiversity conservation, the research on giant pandas and their habitats has received extensive attention (Bai et al., 2020;Hull et al., 2014Hull et al., , 2016Wei et al., 2015). Surprisingly, almost all of these studies were conducted using a single-scale model framework that all covariates are measured at the same spatial scale. These spatial scales are frequently determined arbitrarily by researchers or justified based on expert biological knowledge of the giant panda, such as 20 * 20 m plot size or 250 * 250 m raster cell size (e.g., Feng et al., 2009;Wang et al., 2010). Modeling habitat suitability under a multiscale modeling framework is superior to the single-scale models in terms of model predictive ability and proportion of deviance explained for some species (Bellamy et al., 2013;Timm et al., 2016).
In addition, the multiscale analysis provided new insight on the relationship between species response and habitat covariates that single-scale model did not detect (Mateo-Sánchez et al., 2014;Timm et al., 2016;Wasserman et al., 2012). Therefore, studying habitat selection using a multiscale framework is of significance for deepening our understanding of the relationship between giant pandas and their habitat and formulating more targeted conservation and management strategies.
For more than a decade, traditional statistical methods such as logistic regression have been the dominant method in multiscale habitat modeling . However, in ecological modeling, well-defined issues, such as complicated nonlinear interactions, spatial autocorrelation, high-dimensionality, nonstationarity, and scale, make it difficult for the collected ecological data to meet the assumptions of traditional statistical models (e.g., independence, homogeneity of variance, and multivariate normality), thereby reducing the robustness of the model results (Olden et al., 2008). In the face of these problems, the advantages of machine learning methods gradually emerge and are increasingly widely used in habitat modeling (Evans et al., 2011). Among them, the random forest model stands out because it can deal with a large number of predictors and find sound signals from data with noise. Random forest is an algorithm that developed out of classification and regression tree (CART) and bagging approaches (Breiman, 2001). It builds a classification and regression tree through repeated resampling to form a weak classifier and ensemble many weak classifiers to develop into a strong classifier. Many studies have shown that random forest outperforms traditional statistical methods in terms of model predictive ability (Cushman et al., 2017;Cushman & Wasserman, 2018;Mi et al., 2017), whereas there are few applications of random forest in giant panda habitat modeling, in which the most commonly used are maximum entropy model and habitat suitability index model (e.g., Songer et al., 2012;Xu et al., 2006).
This study tries to identify core habitats and corridors for giant pandas using the improved method by combining multiscale random forest habitat modeling and connectivity analysis. There are three specific objectives of this study: (1) we combine an extensive giant panda occurrence dataset and the multiscale random forest habitat modeling framework to delineate a habitat suitability map for giant pandas in the Qionglai mountain; (2) we used resistant kernel approaches and factorial least-cost path analysis to identify core habitats and corridors; and (3) we assessed the representation of the predicted core habitats and corridors in the protected area network.
Although using giant panda in the Qionglai mountain as a case study, the approaches are expected to provide crucial information for conservation managers to develop more effective conservation strategies for other wildlife species.

| Study area
The Qionglai mountain is located in the west of the Sichuan Basin and is the geographical boundary between Sichuan Basin and Tibetan Plateau. The study area comprises eight counties with a total area around 15,712 km 2 (between 102.26° E and 103.82° E longitude and 29.82°N and 31.72°N latitude; Figure 1). There are six major vegetation/elevation zones in the Qionglai mountain range (Xu et al., 2006): (1) subtropical evergreen broad-leaf forests below 1,600 m; (2) mixed forests of evergreen and deciduous broad-leaf forests between 1,600 and 2,000 m; (3) coniferous and broad-leaf mixed forests between 2,000 and 2,600 m; (4) subalpine coniferous forests between 2,600 and 3,600 m; (5) scrub meadows and alpine talus vegetation between 3,600 and 4,400 m; (6) screes and permanent snow belts above 4,400 m. Many rare wild animals coexist with giant pandas in the Qionglai mountain range, such as golden snubnosed monkey (Phinopithecus roxellana), red panda (Ailurus fulgens), sambar (Cervus unicolor), takin (Budorcas taxicolor), among others .

F I G U R E 1
The geographic extent of the study area and distribution points of the giant panda in the Qionglai Mountain The Qionglai mountain is where the giant panda was scientifically discovered (Hu, 2001). It also is the second-largest tract of the habitat of the six mountain ranges occupied by giant pandas and is home to about 30% of the entire wild giant panda population (State Forestry Administration, 2015). The Qionglai panda population is estimated to be divided into five subpopulations, including Xiaojin, Wolong-Caopo, Xiling-Jiajin, Baishahe, and Sanhe, by major roads

| Presence and pseudo-absence points of giant pandas
The giant panda occurrence data used in this study were from the fourth survey conducted between 2011 and 2014 (Tang et al., 2015).
Based on the third panda survey (1999)(2000)(2001)(2002)(2003) results, the survey area in the fourth survey was classified into key survey areas (2 km 2 survey cell size) and general survey areas (6 km 2 survey cell size). In each survey cell, one line-transect generally greater than 0.75 km was placed to cover all types of panda habitat. Walking along line transects, investigators collected evidence of giant panda occurrences, including feces, foraging sites, dens, footprints, or entities, and used GPS (Global Position System) to record the coordinates of these occurrences (Tang et al., 2015). To reduce the effect of spatial autocorrelation between occurrence points on model performance, we used "SDMtoolbox" to spatially filter panda occurrences (Brown & Anderson, 2014); the filter radius was set to 1.2 km according to the average home range size of giant pandas estimated based on GPS collar study in Wolong reserve (4.4 ± 1.2 km 2 ; Hull et al., 2015).
This means the minimum distance among the filtered occurrences is 1.2 km. Finally, 403 out of 528 giant panda occurrences were retained for modeling ( Figure 1).
The survey results of giant pandas only include the occurrence points and do not include the absence data of giant pandas. In fact, due to the elusive behavior of giant pandas and dense vegetation, it is not practical to confirm the absence of giant pandas in a 2 km 2 or 6 km 2 survey cell . To apply the random forest model in the absence of reliable absence points, we generate a set of pseudo-absence points based on the giant panda occurrence locations (Wang et al., 2010). To fully unleash the power of random forest, we used a random selection of geographically stratified pseudo-absences approach recommended by (Barbet-Massin et al., 2012). The pseudo-absences should lie outside the 3-km radius buffer zone (based on the maximum territory size of the giant panda, which is around 30 km 2 ; Hu, 2001) of panda occurrences. The pseudo-absences are also limited in areas with elevation <4,000 m and slope <50°as giant pandas often avoid high elevation and steep areas (Wang et al., 2010). The minimum distance among pseudo-absences is also set to 1.2 km to alleviate spatial autocorrelation. We randomly select ten sets of the same number of pseudo-absences as available presences (403 in this study) because random forest model is sensitive to class imbalance (Barbet-Massin et al., 2012).

| Environmental variables
According to previous studies, we selected a set of environmental variables that may affect giant panda habitat selection or distribution Wang et al., 2010;Xu et al., 2006). These variables can be summarized into four categories: topographic, land cover, vegetation, and anthropogenic disturbance (Table 1).
Topographic variables included elevation, degree of slope, aspect, and terrain ruggedness index. A 90-m resolution digital elevation model (DEM) product was downloaded from the Shuttle Radar Topography Mission (SRTM; http://srtm.csi.cgiar.org). We calculated these four topographic variables using the Gradient and Geomorphometric Modeling Toolbox in ArcGIS .
To avoid the circular issue of aspect, we transformed the aspect from the range 0-360 to the range 0-1 using the method developed by Roberts and Cooper (1989).
We then used FRAGSTATS v4. 2 (McGarigal, 2002) to calculate four landscape-level metrics (Aggregation Index, AI; Edge Density, ED; Patch Density, PD; Shannon Diversity Index, SHDI) and two class-level metrics for the two dominant forest types (i.e., CNF and CBF) (Largest Patch Index, LPI_; Percentage of Landscape, PLAND_; Table 1) to characterize landscape composition and configuration.
We used a remotely sensed measure of net primary productivity (NPP) obtained from the MODIS (Moderate Resolution Imaging Spectroradiometer) satellite image at 500-m resolution (Running et al., 2004). We calculated a 4-year average of NPP from 2011 to 2014 to be consistent with the time of the panda survey. We also included the bamboo distribution as a vegetation variable as the giant panda is primarily dependent on bamboo (Hu et al., 1985). Ground-based surveys are unavailable for detailed information on bamboo's spatial distribution across large extents. Then, we modeled bamboo distribution using the method developed by Tuanmu et al. (2010). This method extracts 11 phenology metrics from a time series of MODIS satellite images and combines these metrics with the maximum entropy modeling (MaxEnt; a machine learning-based species distribution model) to model the probability of bamboo presence. We used this method to model the probability of bamboo presence and converted it to a binary map of bamboo distribution (i.e., bamboo vs. no bamboo) using the threshold that maximizes the summation of model sensitivity and specificity (Liu et al., 2013). Finally, we calculated the proportion of bamboo coverage within a range of moving windows. See the bamboo modeling details in Appendix A and a binary distribution map for bamboo in Figure A1.
We calculated the Euclidean distance to the village, major road, and minor road for spatial measures of anthropogenic disturbance.
We also calculated the density of villages and roads (major and minor) across various spatial scales. Shapefile of villages and roads were ob- Categorical variables were resampled using the nearest neighborhood method, whereas continuous variables were resampled using the bilinear interpolation method.

| Multiple scale variables
Scale optimization plays a vital role in habitat modeling . We transformed all variables but the distance-based variables (i.e., distance to villages or roads) to multiple scale variables. We considered six spatial scales in the present study, including 1,000, 2,000, 3,000, 4,000, 5,000, 6,000 m; these scales correspond to a spatial extent of 3.14-113 km 2 , which include the average size of giant panda home range (Hull et al., 2015)

| Multiscale random forest habitat model
We used the random forest approach developed by Evans and Cushman (2009) (Svetnik et al., 2004). We compared all subset models and selected the lowest total OOB error as the final model. Before any random forest modeling, we evaluated the minimum number of trees needed by evaluating 2,000 bootstrap samples and observed when the OBB error rate stopped improving. The result showed that OOB error rate ceased to improve after 200 trees ( Figure A3), but we used 500 trees in all models to be conservative as Evans et al. (2011) recommended. Model building and selection were performed using the R package "rfUtilities" (Evans & Murphy, 2014) and "randomForest" (Liaw & Wiener, 2002). Model predictions for the random forest model were generated by creating a habitat suitability map using a ratio of the majority in the votes matrix. We repeated the above steps for the ten presence and pseudo-absence datasets. The final habitat suitability prediction is averaged over the ten models (Barbet-Massin et al., 2012;; this method is called equal-sampling random forest. We also plotted the partial plots for the selected variables by plotting the range of a variable against the estimated probability while keeping other variables at their mean.

| Model evaluation
We assessed the predictive performance of the equal-sampling multiscale random forest model using the area under the total operating characteristic curve (AUC). AUC is a threshold-independent evaluation metrics, it measures the ability of the model to discriminate presences from pseudo-absences (Pearce & Ferrier, 2000).
Presences and pseudo-absences of giant pandas were randomly divided into a training set (70%) and a validation set (30%). This procedure was repeated ten times, and we calculated the mean AUC.

| Landscape resistance layer
A study has shown that the relationship between the resistance species moving in the landscape and the habitat suitability is usually an exponential function rather than a linear function (Keeley et al., 2016). We converted the predicted habitat suitability map from multiscale random forest to the landscape resistance layer using an exponential function (Equation 1): where R is the resistance value, and HS is the predicted habitat suitability. We then rescaled the resistance values to the range between 1 and 100 using linear interpolation, such that the resistance values equal 1 when HS is 1 and 100 when HS is 0. Such transformation means that most pixels in the studied landscape receive low resistance values, and only areas with very low habitat suitability receive high resistance values (Keeley et al., 2016).

| Identification of core habitats and corridors
We used the resistant kernel method (Compton et al., 2007) and factorial least-cost path analysis    (Compton et al., 2007). This surface is a function of source points, landscape resistance, and dispersal ability (Cushman, McRae, et al., 2013). The factorial least-cost path analysis uses Dijkastra's algorithm to calculate the least-cost path from every species occurrence to every other occurrence location in the landscape . These simulated least-cost paths were then buffered based on kernel density estimation, and the Gaussian function was selected in our study. All buffered least-cost paths were summed to produce a map of corridor intensity. The value of a pixel in this map represents the frequency of least-cost paths passing through it.
We used the spatially filtered giant panda occurrences as source points and the transformed habitat suitability map as the resistance layer. A previous study demonstrated that dispersal ability has a significant effect on the estimation of population connectivity . However, we do not have a certain knowledge of giant pandas' dispersal ability. To account for uncertainties regarding giant panda dispersal ability, we used three distance thresholds in the resistant kernel analysis: 6,000, 12,000, 20,000 cost units, which indicate movement distance of 6, 12, 20 km, respectively, in ideal low resistance habitat. We selected 6 km because the biggest home range recorded was about 30 km 2 (Hu, 2001); if the home range was seen as a circle, its diameter is approaching 6 km. In addition, genetic (1) R = 1000 (−1 * HS) studies indicated that the spatial extent of the genetic structure of one population occurred within about 12 km (Hu et al., 2010;Zhan et al., 2007). Furthermore, Pan et al. (2014) and Zhan et al. (2007) reported several giant panda dispersal events with a distance exceeding 20 km. For factorial least-cost paths analysis, we set the upper limit of dispersal ability to 50,000 cost units to model long-distance connections, as the maximum dispersal distance recorded was 50 km (Swaisgood et al., 2010). We used the resistant kernel maps to identify core habitats for the giant panda, as in . We defined core habitats as contiguous areas with resistant kernel values greater than 5% of the highest value.

| Effectiveness of protected areas
To assess the effectiveness of the current protected area network on protecting core habitat patches and corridors for the giant pandas in the Qionglai mountain, we quantified the extent and proportion of predicted core habitats and corridors within protected areas.

| Scale optimization
The univariate random forest optimization analysis showed that optimized scales vary by variables (Table 1; Figure A2); bamboo cover, elevation, net primary productivity, patch density, and largest patch index of closed broad-leaf forest were all strongly related to giant panda occurrence at fine scale (1,000 m), while aspect, terrain ruggedness index, village density, percentage of closed needle-leaf forest were strongly related at much coarser scale (≥4,000 m). Four variables (PLAND_CBF, ED, PLAND_CNF, SLP) were screened due to their high correlation with other variables.

| Multivariate random forest model
We selected the most parsimonious model based on MIR for each presence and pseudo-absence dataset. The most significant variables were the percentage of bamboo cover, elevation, and net primary productivity, while other variables had a relatively low influence on giant panda occurrence ( Figure A4).

| Core habitats and corridor network
We presented resistant kernel simulation results for the giant panda in the Qionglai mountain at three dispersal scenarios (i.e., 6,000, 12,000, 20,000 cost units; Figure 4). Our connectivity simulation showed that high predicted rates of panda movement were mainly concentrated in the northern and central parts of the study area and a relatively small area in the southern region. Dispersal ability showed a significant effect on population connectivity simulation for the giant panda, with a broader range of connected area produced at high dispersal ability scenario (20,000 cost units; Figure 4c) than at low dispersal ability scenario (6,000 cost units; Figure 4a). Under the 6,000 cost units scenario, the giant panda population in the Qionglai mountain was predicted to be broken up into more than ten core patches (three large patches and several small patches; Figure 4a). Under the 12,000 cost units scenario, there were predicted to be two large patches and a few small patches (Figure 4b). Under the highest dispersal ability scenario (i.e., 20,000 cost units), most of the giant panda population was predicted to be connected within one dominant patch while a few small patches were isolated (Figure 4c). However, patches on each side of G318 could not be connected under all dispersal scenarios.
The extent of predicted core habitats varied between 3,451 km 2 and 5,450 km 2 along with dispersal ability, with more remarkable dispersal ability indicating a more significant predicted core habitats extent ( Table 2).
The extent and percentage of protected core habitats differed along with dispersal ability scenarios. The area of protected core habitats varied between 1,485 km 2 and 2,074 km 2 , with the protection rate ranging between 43% and 38% ( Table 2).
The factorial least-cost paths map ( Figure 5) showed that dominant pathway density lies in the mountain area's northern and central parts. The extent of simulated corridors was 3,234 km 2 , of which protected areas covered 1,394 km 2 (43%).

| Scale dependence of habitat selection
Scale is a vital component to consider in ecological research, and grain size is one of its key facets (Connor et al., 2018). Connor et al. (2018) showed that predictor grain size significantly impacts species distribution model accuracy and area of species presence prediction. But there are differences between their study and ours. They constructed species distribution models using variables that are all calculated within the same size moving window (i.e., same spatial scale) at one time (Connor et al., 2018). Strictly, their study should be considered as multiple single-scale models rather than multiscale models because in the latter different variables can be included in the final model at variable scales . The missing step between them is scale optimization.
Scale optimization is critical for robust habitat models, which is crucial in developing conservation and management strategies for endangered species

| The effect of predictors on the distribution of giant pandas
We used an equal-sampling multiscale random forest habitat modeling framework to delineate the relative habitat suitability map for the giant panda in the Qionglai mountain. Random forest is a tree-based method based on "bagging" and is demonstrated to outperform traditional statistical models in the field of species distribution models (Cushman et al., 2017;Evans et al., 2011). In addition, the habitat suitability map derived from random forests is more discriminatory, with higher spatial heterogeneity than predictions from traditional generalized linear models like logistic regression (Cushman & Wasserman, 2018), providing a more robust delineation of priority suitable areas.
Different from the general implementation of random forest in modeling species spatial distribution (Cushman & Wasserman, 2018;Dar et al., 2021;Rather et al., 2020) that only randomly select one set of pseudo-absence sample of size equal to the number of presences, in this study, we used the equal-sampling method. The idea of equal-sampling is to fit n different random forest models (where n is As expected, the percentage of bamboo cover is the most crucial predictor driving the distribution of giant pandas. The giant panda is F I G U R E 3 (a) The habitat suitability map shows giant panda's predicted occurrence based on equal-sampling multiscale random forest habitat modeling in the Qionglai mountain. (b) The landscape resistance map shows the movement resistance for the giant panda, which is transformed from habitat suitability using an exponential function a specialist species with bamboo comprising about 99% of its diet, and it may spend up to 14 h/day foraging bamboo (Schaller, 1985).
Therefore, including biotic interaction with bamboo will improve the performance of habitat models of the giant panda. Modeling understory bamboo distribution in dense forests to a large extent is challenging; however, phenological variables derived from time-series remote sensing images (e.g., MODIS) provide a way to address such an issue .
We found that road density and village density had more significant impacts on the giant panda occurrence probability than the distance to road or village. This finding indicated panda's relative habitat suitability is more related to human activities in the landscape than the proximity to linear roads. Human activity is the primary deterrent to giant pandas' road use; sometimes, low-use roads such as abandoned logging roads were positively related to panda's habitat selection (Qi et al., 2011).
The largest patch index of closed needle forest (LPI_CBF) was positively associated with panda occurrence probability, highlighting that pandas prefer large dense forest patches, a result similar to the result of Wang et al. (2010). Of the selected variables, the landscape composition variables (i.e., LPI_CBF_1000 and LPI_CNF_4000) were more important than the variables reflecting landscape configuration (i.e., AI_4000). This result is similar to other studies on Ursidae (e.g., brown bear; Mateo-Sánchez et al., 2014) and is in agreement with the general pattern that habitat extent is more important than habitat configuration (Cushman & McGarigal, 2002). The giant panda has a high dependence on forest cover and has poor movement ability, suggesting that habitat composition should dominate its habitat relationships.

| Giant panda population connectivity and corridors
Through many years of protection, the population number and habitats of wild giant pandas have increased (State Forestry Administration, 2015), significant species conservation results have been achieved in China (Kang & Li, 2016). However, habitat fragmentation has always been the key factor threatening their survival and is getting worse (Xu et al., 2017). A few studies assessed population connectivity and proposed corridors to link fragmented habitat patches Qi et al., 2012;Wang et al., 2021). The method usually used in these studies was the least-cost analysis Qi et al., 2012), which simulated narrow linear corridors and the structural connectivity among habitat patches. Giant pandas may not use those simulated corridors because there may be no individuals in predicted habitat patches or giant pandas are hard to traverse long corridors due to limited mobility. We integrated panda's dispersal ability into connectivity analysis, which previous studies usually ignore. One strength of the resistant kernel approach is its explicit and realistic incorporation of species dispersal ability . There are predicted to be >10 core habitat patches if the dispersal ability of giant panda is limited to 6,000 cost units, but with the dispersal of 12,000-20,000 cost units, it would result in seven to four patches ( Figure 4). This result highlights that the extent and fragmentation of connected habitats are highly dependent on the dispersal ability of the focal species (also see Ashrafzadeh   (Kang & Li, 2016). Furthermore, spatialexplicit resistant kernel estimation can provide more information in the zoning or the effectiveness evaluation of protected areas  than simple habitat suitability distribution map, which is often used in such assessments Wang et al., 2021). Under all dispersal ability scenarios, there is predicted to be a large proportion (57%-62%) of core habitats that are not protected by the current nature reserve network (Figure 4), highlighting a great potential to establish new protected areas.
The factorial least-cost path analysis identified optimal routes between giant panda occurrences to facilitate connectivity. General corridor simulation methods (e.g., least-cost path or least-cost corridor) take habitat patches as the "source" and then calculate the path with the least cumulative cost between source patches; its result only reveals the location of the corridor (Cushman, McRae, et al., 2013). However, the factorial least-cost path approach simulates corridor network based on species occurrence, and the simulation result provides the location and intensity of the corridor. The corridor intensity is a kernel density estimation based on the number of least-cost paths. Such information is important for corridor priority planning as corridor building and restoration usually need to invest a lot of money and manpower. We recommend paying more attention to the corridors linking small populations with much higher extinction risk to large populations. For example, subpopulations at the southernmost of the study area need more attention as they are predicted to be isolated from large core habitat patches under all dispersal ability scenarios. In previous evaluations, these subpopulations were also believed to be isolated from other subpopulations by national roads (G318) (State Forestry Administration, 2015;Xu et al., 2006). In addition, areas with high predicted least-cost paths frequency outside protected areas also need prior protection. We can identify barriers that may impede giant panda dispersal based on the corridor pathway, like major roads. Combined with corridor density, we can further locate and rank the intersection of corridors with roads, which can provide crucial information for conservation practitioners to implement road mitigation measures such as warning signs, reduced speed limits, fencing, and construction of crossing structures (Cushman, Lewis, et al., 2013;Zeller et al., 2020).

| Caveats and limitations
Our findings must be interpreted with regard to several major considerations. First, giant panda occurrences used in this study were from the fourth survey conducted between 2011 and 2014. Natural and socio-economical conditions have changed since then (Xu et al., 2017). For example, the reduction of total and agricultural population and the increasing of infrastructure development (e.g., hydropower stations and road construction should also be noted that habitat suitability is not a good proxy for landscape resistance, as habitat suitability reflects habitat selection in home range while species may use the landscape differently during dispersal movements (Keeley et al., 2017). It would be better to develop resistance models with movement (Zeller et al., 2018) or gene flow (Cushman et al., 2006) data. Lastly, when designing new protected areas or corridors, it is also necessary to consider other large carnivores, as long-term monitoring studies have revealed a wide distribution range retreat of large carnivore populations across the giant panda distribution range (Li et al., 2020).

ACK N OWLED G M ENTS
We gratefully thank all participants in the Fourth Giant Panda Surveys

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

AUTH O R CO NTR I B UTI O N S
Xue Sun: Conceptualization (lead); data curation (equal); formal analysis (equal); methodology (equal); writing -original draft (lead). Zexu Long: Conceptualization (equal); formal analysis (equal); methodology (lead); software (equal); visualization (equal); writing -review and editing (equal). Jingbo Jia: Data curation (equal); funding acquisition (equal); resources (lead); supervision (equal); writing -review and editing (equal).  & Jönsson, 2012;Tuanmu et al., 2010;Viña et al., 2010). The 3-year averages of these phenology metrics from 2011 to 2013 were used to map bamboo distribution. In addition, we included the elevation as a variable for modeling bamboo distribution because vegetation type is closely related to elevation zones in the Qionglai mountain range (Xu et al., 2006). We used the software Maxent (version 3.4.1, Phillips et al., 2006) to build a model for mapping the overall bamboo distribution as in Tuanmu et al. (2010). We did not have random bamboo presence localities across the study area, and we used two types of points as surrogate bamboo occurrence. One is the giant panda occurrence localities associated with bamboo from the Fourth Giant Panda Survey in Sichuan Province. The other is 500 points randomly selected from a coarse bamboo distribution map obtained from the 4th National Giant Panda Survey. We randomly selected 75% of occurrence data as the training dataset, and the remaining 25% were used as the validation dataset. The procedure was repeated ten times to get an average prediction of bamboo occurrence probability. We set other parameters in Maxent as default. To get a binary map of bamboo distribution, we used the threshold that maximizes the sum of sensitivity and specificity to convert the continuous occurrence probability map into a binary map (Liu et al., 2013).

F I G U R E A 1
Binary map for bamboo distribution in Qionglai mountain predicted by using MaxEnt and phenological metrics

F I G U R E A 2
The frequency scale was selected as the optimal scale in 10 univariate equal-sampling random forest models for each variable

F I G U R E A 3
The number of trees versus the error rate plot F I G U R E A 4 Bar plot of variable importance based on the mean Model Improvement Ratio (MIR) from random forests of 10 equal-sampling presencepseudo absence datasets