Environmental Gradients Shape Mammal and Galliform Bird Communities in a Mountain Reserve Through Species Turnover and Niche Differentiation

1. Introduction

Since the Anthropocene [1], human impacts on ecosystems have intensified, while global biodiversity has undergone persistent decline [2]. Many researchers have suggested that the Earth has already entered a sixth mass extinction [3,4]. Against this background, maintaining biodiversity and designing effective conservation strategies have become increasingly complex [5]. Although protected areas are widely regarded as the core spatial units for biodiversity conservation [6,7], they are not internally homogeneous “safe habitats.” Instead, they usually encompass complex gradients shaped by environmental factors such as elevation, vegetation, and climate, while also being continuously influenced by multiple forms of human activity, including roads, settlements, grazing, poaching, collecting, and tourism [8,9,10]. Therefore, treating protected areas simply as relatively uniform conservation units is no longer sufficient to explain the spatial variation of biodiversity within them. In contrast, greater attention to the internal organization of communities within protected areas and their relationships with environmental gradients, as well as to how natural environmental heterogeneity and human disturbance jointly drive such spatial differentiation, is essential for understanding the mechanisms that maintain multispecies coexistence and for improving conservation management [11,12,13].

Community composition and its spatial variation along environmental gradients are closely linked to how biodiversity is maintained within protected areas and to the priorities of conservation management [14,15,16]. A key question is whether community composition and species associations are organized as clearly bounded discrete modules or instead vary continuously along environmental gradients. If communities are structured into relatively independent modules, conservation management may focus more on identifying and maintaining key habitat units. In such cases, differences among communities are more likely to reflect nestedness, with species-poor communities representing subsets of species-rich communities, making large habitat patches particularly important for conservation in fragmented systems [17,18]. By contrast, if communities primarily vary continuously along environmental gradients, differences among communities are more likely to arise from species turnover. In this case, conservation should place greater emphasis on maintaining environmental heterogeneity, habitat continuity, and connectivity, for example, through multi-scale strategies or complementary protected-area networks across environmental gradients [19,20,21].

Animal communities are often shaped jointly by environmental heterogeneity and human disturbance, a process that is particularly pronounced in mountain protected areas [22]. As elevation and topography change, environmental factors such as vegetation productivity, habitat structure, and hydrothermal conditions vary systematically, potentially leading different species to exhibit differentiated distributions along environmental gradients [23]. At the same time, human activities, including roads, built-up areas, and livestock grazing, continue to influence landscapes both within and around protected areas [24], thereby further altering species distributions and community structure [25]. Therefore, analyses based solely on overall community composition are insufficient to fully reveal the mechanisms underlying multispecies coexistence. It is also necessary to consider species-specific responses to environmental factors and human disturbance, and to examine whether species exhibit differentiated habitat selection and environmental preferences within shared spatial contexts [26,27]. By linking community-level patterns with species-level environmental responses, it becomes possible to more effectively determine whether communities vary continuously along environmental gradients or instead contain groups of species sharing similar ecological preferences [28], thereby providing insights into the mechanisms that maintain multispecies coexistence within protected areas [29].

Liziping Nature Reserve in Sichuan, China, provides an ideal setting for investigating the spatial organization of communities and their underlying drivers in mountain protected areas. The reserve harbors a high diversity of endangered, threatened, and endemic species [30], but conservation efforts have long been centered primarily on the giant panda and its habitat [31,32]. At the same time, the reserve contains pronounced environmental gradients, while surrounding agriculture, grazing, and tourism continue to influence its ecological processes [30]. Under these conditions, it remains unclear whether communities within the reserve are organized as discrete modules or instead vary continuously along environmental gradients. This distinction has direct implications for conservation management, determining whether priority should be given to key habitat patches or to the maintenance of environmental heterogeneity and habitat connectivity. Therefore, examining community organization and its drivers at the scale of Liziping Nature Reserve will not only improve understanding of the mechanisms underlying species coexistence in mountain ecosystems, but also provide a scientific basis for more refined management of mountain protected areas.

In this study, we used camera-trap data to examine the spatial organization of mammal and galliform bird communities in the Liziping Nature Reserve at two complementary levels. First, using analyses of beta-diversity and community composition, we asked whether community structure across sites is better characterized as discrete modules or as a continuous gradient along environmental variation. Second, using species co-occurrence patterns and occupancy models, we asked whether species differ in their responses to environmental gradients and human disturbance, and whether similarity in species co-occurrence and activity structure is associated with similarity in environmental responses. By linking community-level structure with species-level environmental responses, we aimed to clarify the mechanisms underlying multispecies coexistence in a mountain protected area and to provide a basis for more community-inclusive conservation planning.

2. Materials and Methods

2.1. Study Area

Liziping Nature Reserve is located in the Xiaoxiangling Mountains on the eastern edge of the Himalaya–Hengduan Mountains region in Sichuan, China. The reserve is dominated by mid- to high-elevation mountainous terrain, with deeply dissected topography, pronounced vertical relief, and an elevational range of approximately 1,330–4,551 m [33]. Along this elevational gradient, vegetation transitions from montane evergreen broad-leaved forest and mixed conifer–broadleaf forest to subalpine coniferous forest, subalpine shrubland, and alpine shrub meadow, forming a distinct pattern of vertical vegetation zonation [30]. Several villages are distributed around the reserve, where local residents have long engaged in agriculture, livestock grazing, and tourism. As a result, human activities such as roads, built-up land, and grazing have exerted persistent influences on the landscape structure of the reserve.

2.2. Infrared Camera Deployment

From July 2018 to June 2019, infrared camera traps were deployed throughout the reserve, with a minimum spacing of 800 m between cameras. Cameras were preferentially placed at sites with frequent animal activity, such as mountain trails, areas near water sources, and open sites under forest cover, to maximize detection probability. Each camera was mounted on a tree trunk or other stable substrate at approximately 50 cm above the ground and oriented toward likely animal travel routes. Cameras operated continuously for 24 h per day and, when triggered, recorded three consecutive photographs and one 10 s video, with a 5 s interval between triggers. Consecutive records of the same species at the same camera within 30 min were treated as a single independent event. Camera identity, species identity, and detection time were extracted from each independent record for subsequent analyses.

To quantify variation in species activity across camera-trap sites, we calculated a relative abundance index (RAI) for each species at each camera location. RAI was defined as the number of independent detection events per 100 effective camera-trap days. Independent events were defined using a 30 min interval criterion between consecutive records of the same species at the same camera, and effective camera-trap days were calculated as the total number of days during which each camera was operational. The resulting RAI matrix was used to calculate Bray–Curtis dissimilarities among sites for site-based analyses and, after matrix transposition, among species for species-level analyses. For subsequent community and occupancy analyses, species with fewer than 20 independent detections and recorded at fewer than 10 camera sites were excluded to reduce instability associated with sparse data.

2.3. Environmental Data

Land-cover data were obtained from the ESA WorldCover 10 m 2021 v200 dataset [34], which provides global information on the distribution of major land-cover types at a spatial resolution of 10 m. Within the study area, relevant land-cover classes were extracted and the shortest distance from each camera-trap site to patches of different land-cover types was calculated to characterize habitat structure and human-modified environments. Based on the actual distribution of land-cover types in the study area and their ecological relevance, representative and widely distributed classes were selected, including forest (tree), shrubland (shrub), grassland (grass), cropland (crop), built-up areas (built), bare land (bare), and water bodies (water). Distance variables were calculated for each class, resulting in the predictors dist_tree, dist_shrub, dist_grass, dist_crop, dist_built, dist_bare, and dist_water, which were used to represent the potential influence of habitat types and human land use on community distribution patterns. To improve the characterization of aquatic environments, we also incorporated spatial data on rivers and streams provided by the Forestry Department of Sichuan Province [35]. The final water-distance variable was defined as the minimum distance from each site to either ESA WorldCover water bodies or mapped rivers and streams, thereby improving the completeness and accuracy of water-related spatial information.

Topographic variables were derived from a 30 m resolution digital elevation model obtained from the Geospatial Data Cloud of the Chinese Academy of Sciences (http://www.gscloud.cn). In addition to elevation (ele), terrain variables including slope, northness, terrain ruggedness index (TRI), and topographic wetness index (TWI) were calculated from the elevation raster. All topographic variables were retained at a spatial resolution of 30 m for subsequent analyses.

Vegetation productivity was characterized using the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) derived from Landsat 8 surface reflectance imagery (Collection 2 Level-2) available in Google Earth Engine. To represent peak growing-season conditions while reducing interannual variability, we used images acquired during June–September from 2016 to 2020. Cloud, cloud-shadow, and snow pixels were masked using the QA_PIXEL quality band, and surface reflectance bands were rescaled with standard Landsat scaling factors prior to index calculation. NDVI and EVI were calculated as follows:

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{e}\mathrm{d}}}$$

(1)

$$\mathrm{EVI}=2.5\times {\displaystyle \frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+6 \mathrm{Red}-7.5 \mathrm{Blue}+1}}$$

(2)

where NIR, Red, and Blue represent the surface reflectance in the near-infrared, red, and blue spectral bands, respectively. For each index, a median composite of all valid observations during the study period was generated to obtain stable estimates of vegetation productivity.

Recent anthropogenic disturbance (hereafter disturbance) was quantified at each site using a camera-trap-based relative abundance index (RAI) derived from independent detections of humans and livestock.

All continuous environmental predictors were standardized (mean = 0, SD = 1) prior to model fitting. The quadratic elevation term (ele2) was calculated as the square of the standardized elevation variable to reduce collinearity between linear and quadratic terms.

To reduce multicollinearity among environmental predictors, variance inflation factors (VIFs) were calculated for all candidate variables. Variables with VIF values greater than 5 were considered highly collinear and excluded from subsequent analyses. TRI and NDVI were initially included as candidate variables but were removed after collinearity screening. The remaining variables were retained for subsequent analyses (Figure S1, Table S1).

2.4. Analytical Framework and Data Structure

To distinguish community-level variation among sites from similarity among species, we used two complementary data structures in subsequent analyses. For community-level analyses, we used a site-by-species matrix based on species presence–absence data. For species-level analyses, we used the transposed species-by-site matrix, from which Jaccard dissimilarity was calculated from occurrence data and Bray–Curtis dissimilarity from species-specific RAI across sites.

2.5. Beta-Diversity and Null Model Analyses

Community dissimilarity among sites was quantified using Jaccard beta diversity based on species presence–absence data and partitioned into turnover and nestedness components following Baselga [36] using betapart [37]. To test whether the observed level of community differentiation deviated from random expectations, we randomized the community matrix while preserving site richness and species occurrence frequencies using the swap algorithm in vegan [38]. We performed 9,999 randomizations and compared the observed mean pairwise Jaccard dissimilarity with the null distribution.

2.6. Clustering Analysis

To examine whether species-level co-occurrence and relative activity patterns formed discrete modules, hierarchical clustering analyses were performed on the species-level Jaccard and Bray–Curtis distance matrices. Candidate solutions with k = 2–14 groups were evaluated using mean silhouette width to assess cluster separation. To determine whether the observed species arrangement departed from random expectations, we calculated a Δ statistic, defined as the difference between mean between-group and within-group dissimilarities for a given clustering solution. Larger Δ values indicate stronger differentiation among groups relative to within-group similarity, but do not by themselves imply the presence of well-separated discrete modules. Observed Δ values were compared with a null distribution generated by randomly assigning species to groups while preserving group sizes. Thus, Δ analysis was used to test whether the observed structure was more organized than expected by chance, whereas silhouette width was used to evaluate the strength of cluster separation. Significance was assessed using 9,999 randomizations.

2.7. NMDS Ordination and Environmental Fitting

NMDS ordinations were performed on the species-level Jaccard and Bray–Curtis distance matrices to examine patterns of species co-occurrence and relative activity structure. Analyses were conducted using the metaMDS function in vegan [38], and environmental vectors were fitted with envfit. Each ordination point represented a species rather than a camera-trap site. Correlation strength (r2) and significance were assessed using 9,999 permutations. Environmental predictors included the variables retained after VIF-based collinearity screening, representing topography, vegetation productivity, habitat structure, hydrological conditions, and human disturbance.

2.8. Occupancy Analysis

We used camera-trap data to model species occupancy across the study area. Detection histories were constructed by dividing the monitoring period into 14-day sampling occasions, with each species coded as detected or not detected at each site during each occasion. To examine spatial variation in occupancy, the model included 13 environmental predictors: elevation (ele), elevation squared (ele2), slope (slope), northness (northness), topographic wetness index (twi), vegetation productivity (evi), distance to forest (dist_tree), shrubland (dist_shrub), grassland (dist_grass), bare land (dist_bare), cropland (dist_crop), built-up areas (dist_built), and recent anthropogenic disturbance (disturb). Since elevation was standardized prior to analysis and the quadratic term (ele2) was calculated as the square of standardized elevation, the signs of the linear and quadratic elevation coefficients can be used to classify species-specific elevational response curve types (Table S2). Occupancy probability was modeled as:

$$\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}\left({\mathsf{\psi }}_{\mathrm{s}\mathrm{i}}\right)={\mathsf{\alpha }}_{\mathrm{s}}+\sum _{\mathrm{v}=1}^{\mathrm{V}}{\mathsf{\beta }}_{\mathrm{s}\mathrm{v}}{\mathrm{Z}}_{\mathrm{i}\mathrm{v}}$$

(3)

where ψsi represents the occupancy probability of species s at site i, αs is the species-specific intercept, βsv denotes the regression coefficient describing the effect of environmental predictor v on species s, and Ziv represents the value of environmental variable v at site i. Detection probability (pi) was modeled as a function of camera-trap operational days, distance to water (dist_water), and camera type:

$$\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}\left({\mathrm{p}}_{\mathrm{s}\mathrm{i}\mathrm{j}}\right)={\mathsf{\gamma }}_{0\mathrm{s}}+{\mathsf{\gamma }}_{1\mathrm{s}}\hspace{0.17em}\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{k}\mathrm{d}\mathrm{a}{\mathrm{y}}_{\mathrm{i}\mathrm{j}}+{\mathsf{\gamma }}_{2\mathrm{s}}\hspace{0.17em}\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}\_\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}{\mathrm{r}}_{\mathrm{i}}+{\mathsf{\gamma }}_{3\mathrm{s}}\hspace{0.17em}\mathrm{c}\mathrm{a}\mathrm{m}\_\mathrm{t}\mathrm{y}\mathrm{p}{\mathrm{e}}_{\mathrm{i}}$$

(4)

where psij represents the probability of detecting species s at site i during sampling occasion j, γ0s is the intercept, and γ1s–γ3s are regression coefficients describing the effects of detection covariates on detection probability. Here, workdayij denotes the number of operational camera-trap days at site i during occasion j.

2.9. Clustering of Species Environmental Response Coefficients

To examine whether species could be grouped according to their environmental response patterns, we performed hierarchical clustering using species-specific environmental regression coefficients derived from the occupancy models. To reduce the influence of coefficient uncertainty on clustering results, these coefficients were subjected to uncertainty-weighted shrinkage prior to clustering. Specifically, each coefficient was multiplied by a shrinkage factor:

$$\mathrm{E}({\mathsf{\beta }}_{\mathrm{i}\mathrm{v}}\mid {\widehat{\mathsf{\beta }}}_{\mathrm{i}\mathrm{v}})={\widehat{\mathsf{\beta }}}_{\mathrm{i}\mathrm{v}}\cdot {\displaystyle \frac{{\mathsf{\tau }}^2}{{\mathsf{\tau }}^2+{\mathsf{\sigma }}_{\mathrm{i}\mathrm{v}}^2}}$$

(5)

where ${\widehat{\mathsf{\beta }}}_{\mathrm{i}\mathrm{v}}$ is the posterior mean of the regression coefficient for species i and environmental predictor v, ${\mathsf{\sigma }}_{\mathrm{i}\mathrm{v}}^2$ is the posterior variance of that coefficient, and ${\mathsf{\tau }}^2$ is the overall variance of regression coefficients for predictor v across species. This procedure shrinks highly uncertain coefficients toward zero, thereby reducing the influence of statistical noise on multivariate clustering structure and allowing groups to better reflect stable ecological response patterns. Clustering was based only on occupancy-level environmental regression coefficients and did not include detection parameters. The shrinkage-adjusted beta coefficient matrix was first standardized, and pairwise distances among species were then defined as one minus the correlation coefficient. Hierarchical clustering was subsequently performed using the Ward.D2 method. Candidate clustering solutions were evaluated using mean silhouette width, and the optimal number of clusters was selected as the k value with the highest silhouette score. PERMANOVA was then performed on the corresponding distance matrix, and 9,999 permutations were used to assess whether species groups occupied distinct regions of environmental response space.

2.10. Mantel Analysis

To compare species-level patterns derived from different data types, we calculated pairwise Jaccard distances among species based on site-level occurrence data, Bray–Curtis distances based on site-level RAI values, and distances based on environmental response coefficients (beta). Mantel tests with 9,999 permutations were used to assess concordance among these species-level distance matrices.

All statistical analyses were conducted in R 4.5.2 [39], using the packages vegan [38], betapart [37], cluster [40], and rstan [41].

3. Results

A total of 148 camera-trap sites were deployed in this study, yielding 24,678 camera-trap days and 4,065 independent detections. After excluding species with fewer than 20 independent detections and records from fewer than 10 camera-trap sites, 15 species were retained for analysis (Table S3). Species richness at individual camera-trap sites ranged from 1 to 10, with a mean of 3.9348 (Figure 1). Species richness varied markedly among sites and showed an overall right-skewed, long-tailed distribution (Figure S2).

3.1. Community Beta Diversity

Beta-diversity decomposition showed that differences among communities were driven primarily by species turnover (mean turnover = 0.7233), whereas the nestedness component was relatively low (mean nestedness = 0.1791). Turnover accounted for 77.09% of total beta diversity, indicating that variation among sites was dominated by species replacement. To evaluate whether the observed level of community differentiation deviated from random expectations, we compared the observed mean pairwise Jaccard dissimilarity with a null distribution generated by the swap algorithm while preserving both site richness and species occurrence frequencies. The observed value did not differ significantly from the null expectation (Pr = 0.7623), suggesting that the overall magnitude of community differentiation was consistent with random expectations given site richness and species occurrence frequencies.

3.2. Species Similarity and Clustering

We further examined whether species co-occurrence patterns formed discrete modules. Hierarchical clustering analyses were conducted using Jaccard dissimilarity based on species occurrence across camera-trap sites (Figure 2a) and Bray–Curtis dissimilarity based on species-specific RAI values across sites (Figure S3a). Across the range of k = 2–14, mean silhouette values were consistently below 0.1, indicating weak separation among groups and providing little support for well-defined discrete modules. At the same time, Δ analysis (Figures S4 and S5) showed that the observed grouping structure differed significantly from random expectations at all k values (p < 0.001), indicating that similarity structure among species was non-random. Taken together with the low silhouette values, these results suggest that the observed structure is better interpreted as weak differentiation or continuous variation rather than strongly separated modular groups.

3.3. NMDS Ordination and Environmental Gradients

NMDS ordination further supported this interpretation, indicating that species co-occurrence and relative activity patterns were better characterized by continuous variation associated with environmental gradients than by clearly separated discrete modules. Two-dimensional ordinations based on Bray–Curtis and Jaccard dissimilarities both showed species distributed continuously in ordination space, with stress values of 0.2004 and 0.1969, respectively, suggesting gradual variation in co-occurrence relationships and relative activity structure among species (Figure 2b and Figure S3b). Environmental fitting analysis (Table S4) showed that environmental variables explained both ordination patterns consistently and relatively strongly. Elevation was the strongest correlate of species ordination structure, with r2 values of 0.3135 and 0.3343 for the Jaccard- and Bray–Curtis-based ordinations, respectively (both p < 0.001). Vegetation productivity (EVI) was the second strongest correlate, with r2 values of approximately 0.21 in both ordinations (both p < 0.001). Other variables significantly associated with species ordination structure included habitat variables, such as distance to forest, bare land, grassland, and shrubland, as well as human-related variables, including distance to built-up areas and recent anthropogenic disturbance, with r2 values ranging from 0.06 to 0.12 (all p < 0.001). In contrast, TWI, slope, northness, and distance to water showed relatively weak explanatory power (all r2 < 0.03) and were mostly non-significant (p > 0.05).

At the species level, NMDS ordinations based on Jaccard and Bray–Curtis dissimilarities showed broadly concordant patterns (Figure 2b and Figure S3b). Most species occupied similar relative positions in the two ordination spaces, indicating that the main axes of ecological differentiation were largely consistent across distance metrics. For example, blood pheasant (Ithaginis cruentus), Chines red panda (Ailurus styani), and forest musk deer (Moschus berezovskii), and Chinese serow (Capricornis milneedwardsii) were positioned primarily along the high-elevation and forest-associated gradient, whereas Tibetan macaque (Macaca thibetana), Asiatic black bear (Ursus thibetanus), and Lady Amherst's pheasant (Chrysolophus amherstiae) were more closely associated with the vegetation productivity (EVI) gradient. Although a few species showed minor shifts between the two ordinations, the overall pattern was highly similar, indicating that species co-occurrence and relative activity structure exhibited a stable pattern of differentiation along environmental gradients.

3.4. Occupancy Analysis

The occupancy model converged well, with all environmental variable parameters having Rhat values below 1.01, indicating stable estimation. Of the 15 species, all but the red panda, forest musk deer, and leopard cat showed significant effects for at least one environmental variable (Figure 3). Although some species did not meet traditional significance thresholds, posterior probabilities of the direction of effects indicated that all species, except the forest musk deer, showed clear directional responses to at least one environmental variable (posterior probability > 0.9 or < 0.1), suggesting distinct environmental response differentiation.

Elevation-related variables had the strongest influence on species differentiation. Standardized elevation (ele) and its squared term (ele2) revealed different response patterns, such as low elevation preference, high elevation preference, and hump-shaped responses at mid-elevation. Habitat structure and human activity variables also showed species-specific effects, while slope, aspect, EVI, and TWI had significant effects only for some species. These findings suggest ecological niche differentiation rather than a single, consistent environmental filtering pattern.

Specifically, in the study region, Lady Amherst’s pheasant (Chrysolophus amherstiae) and Tibetan macaque (Macaca thibetana) showed hump-shaped responses with optima at relatively low elevations, whereas Chinese serow (Capricornis milneedwardsii) exhibited a hump-shaped response with an optimum at higher elevations. In contrast, blood pheasant (Ithaginis cruentus) showed a monotonic increase in occupancy with elevation, while Asiatic black bear (Ursus thibetanus) and Temminck’s tragopan (Tragopan temminckii) displayed hump-shaped responses with optima at mid-elevations. For habitat structure variables, Temminck’s tragopan was associated with sites closer to bare land and farther from forest, whereas wild boar (Sus scrofa) was associated with areas closer to grassland. Tufted deer (Elaphodus cephalophus) showed opposite responses to distances from shrubland and grassland. Regarding human-related variables, distance to built-up areas had a negative effect on Asiatic black bear but positive effects on tufted deer, wild boar, and Malayan porcupine (Hystrix brachyura). Recent disturbances such as livestock grazing had a positive effect on giant panda (Ailuropoda melanoleuca), but negative effects on tufted deer and Malayan porcupine.

Hierarchical clustering based on shrinkage-adjusted species-specific environmental regression coefficients showed that mean silhouette width was highest at k = 2, with a value of 0.3311, indicating that a two-cluster solution provided the best separation among the candidate groupings (Figure 4). PERMANOVA further suggested that these two clusters occupied distinct regions of multivariate environmental response space (R² = 0.4355, F = 10.0279, p < 0.001), indicating that the identified groups corresponded to different environmental response patterns. The first cluster included Lady Amherst’s pheasant, leopard cat, Tibetan macaque, red fox, Temminck’s tragopan, yellow-throated marten, tufted deer, wild boar, and Malayan porcupine, whereas the second cluster included giant panda, Asiatic black bear, forest musk deer, Chinese red panda, blood pheasant, and Chinese serow.

3.5. Mantel Test

Mantel tests showed no significant correlation between Jaccard distances based on presence–absence data and Bray–Curtis distances based on species-specific RAI values (r = 0.0389, p = 0.3932). In contrast, Jaccard distances were positively correlated with distances derived from shrinkage-adjusted environmental response coefficients (r = 0.4085, p < 0.001), indicating concordance between species co-occurrence structure and similarity in environmental responses. Bray–Curtis distances were not significantly correlated with beta distances (r = −0.0771, p = 0.7593), suggesting that relative activity structure did not vary in parallel with either co-occurrence structure or environmental response differentiation.

4. Discussion

Understanding the spatial organization of communities within protected areas and its environmental drivers is important for elucidating the mechanisms that maintain biodiversity and for improving conservation management [16,42]. This study shows that the mammal and galliform bird communities in Liziping Nature Reserve do not exhibit distinct boundary-separated structures, but instead resemble a continuous gradient pattern shaped by environmental factors. Beta diversity was primarily driven by species turnover, indicating that the community differences between sites were largely due to species replacement rather than simple nestedness. Clustering analyses based on species co-occurrence and relative activity intensity did not identify well-separated species modules. Furthermore, NMDS ordination revealed that both species co-occurrence and relative activity structures changed continuously along environmental gradients. Together, these results suggest that the communities in the study area are not composed of clearly distinct, discrete units but are more likely the result of continuous environmental selection processes [43,44]. Occupancy models further indicated that species responded differently to environmental gradients, suggesting that coexistence in this area is more reliant on niche differentiation and environmental heterogeneity than on simple discrete partitioning. Taken together, the findings suggest that the community in Liziping Nature Reserve is better described as a continuous gradient structure rather than a discrete modular one. Environmental heterogeneity and the niche differentiation it promotes are likely key mechanisms maintaining multispecies coexistence in this reserve.

4.1. Environmental Gradients and Community Turnover

The turnover-dominated pattern observed in Liziping Nature Reserve indicates substantial environmental heterogeneity within the reserve. When considered together with the ordination and environmental fitting results, elevation emerged as the primary factor structuring species association patterns in the study area, consistent with a broad body of research on mountain communities [45,46]. Vegetation productivity, habitat structure, and human-related variables also contributed to continuous environmental variation within the reserve. Because species differ in their environmental tolerances and resource-use patterns, different sites are more likely to support distinct species assemblages, leading to continuous species replacement along environmental gradients rather than simple gains or losses of species. Such turnover may help maintain community diversity within protected areas [47,48].

The null model results further suggest that this turnover pattern does not imply strong discrete community partitioning within the reserve, but is more likely the outcome of the combined effects of species occurrence frequencies, site-level richness, and environmental gradients. This finding also implies that conservation management should not focus solely on local richness hotspots, but should also prioritize the environmental heterogeneity maintained by elevational, habitat, and disturbance gradients. Avoiding landscape homogenization or overly uniform management may therefore be important for preserving the diverse habitat conditions required by different species [49,50].

4.2. Species-Specific Environmental Responses and Niche Differentiation

The occupancy model results indicate that species in the study area do not follow a single, consistent pattern of environmental filtering; instead, they exhibit pronounced species-specific environmental responses. Species responded differently to elevation, vegetation productivity, habitat structure, and human disturbance, suggesting that coexistence is more likely maintained through differentiated environmental preferences and resource-use strategies within a shared landscape [51,52]. For example, species exhibited contrasting responses along the elevational gradient, including preferences for low elevations, high elevations, or hump-shaped responses at mid-elevations. Responses to habitat structure and human disturbance variables were likewise inconsistent or even opposite among species. These patterns suggest that multispecies coexistence in the reserve is more likely supported by niche differentiation under environmentally heterogeneous conditions rather than by a uniform environmental filtering process.

It is also notable that the proportion of environmental effects reaching conventional significance thresholds in the occupancy models was relatively limited, whereas directional differentiation among species was widespread. Clustering based on shrinkage-adjusted β coefficients further showed that species could be divided into two groups with contrasting environmental response patterns, indicating substantial heterogeneity in ecological responses within the community. Together, these results support the interpretation that multispecies coexistence in the study area is structured by differentiated environmental responses among species.

4.3. Multiple Ecological Dimensions Revealed by Different Data Types

Mantel tests further indicated that species relationships inferred from different data types were not entirely consistent. Jaccard distances were significantly correlated with distances derived from environmental response coefficients, whereas Bray–Curtis distances were not correlated with either of them. This suggests that the species co-occurrence structure characterized by presence–absence data is broadly consistent with patterns of environmental response differentiation, whereas similarity based on relative activity intensity reflects a partially independent ecological dimension.

In other words, the environmental response patterns that determine whether a species can occur at a given site do not necessarily determine how intensively it uses that site once present. Because both the occupancy model and the Jaccard analysis are based on presence–absence data, occurrence information can effectively reveal environmental suitability and broad distribution patterns [53], but may not fully capture the degree to which species actually use available habitats [54]. Consequently, when evaluating habitat quality and conservation priorities, relying solely on occurrence data may underestimate the complexity of community organization, and complementary information on habitat use or activity intensity should also be considered.

4.4. Human Disturbance

Although natural environmental gradients remained the primary drivers of community structure, the influence of human disturbance should not be overlooked. In this study, distance to built-up areas and recent disturbance intensity represented two different types of anthropogenic pressure. The former mainly reflects long-term landscape modification associated with roads, settlements, and agricultural land, which occurs primarily near the reserve's boundaries. In contrast, recent disturbances represent short-term or recurrent human activities, such as livestock grazing and other human presence. These disturbances can penetrate more diffusely into the interior of the reserve and influence local habitat conditions and species activity at finer spatial scales. The relatively weak correlation between these two variables suggests that different types of human disturbance may operate through largely independent spatial distributions and ecological mechanisms.

Combining the co-occurrence ordination and occupancy results shows that species responded differently to these two types of disturbance. For example, Asiatic black bears exhibited a significant negative response to proximity to built-up areas, whereas tufted deer and Malayan porcupines showed significant negative responses to recent disturbance. These patterns indicate that human activities still exert clear suppressive effects on some species. Previous studies have reported that livestock grazing can pose substantial threats to wildlife such as the giant pandas in some reserves [55,56], although other studies suggest that livestock presence does not necessarily have significant effects on species such as giant pandas, Chinese serows, or forest musk deer [56]. In the present study, giant pandas did not show a clear negative response to recent disturbance. One possible explanation is that human activities within panda habitats in Liziping Nature Reserve have been relatively well managed during the study period, so disturbance intensity may not have reached the threshold required to trigger strong avoidance behavior. Under the current conditions of the study area, the effects of human disturbance are therefore more likely to depend on disturbance type, intensity, and duration rather than manifest uniformly as a negative effect across species.

5. Conclusions

Overall, mammal and galliform bird communities in Liziping Nature Reserve did not form clearly bounded discrete modules, but were instead primarily organized as continuous variation along environmental gradients, with community differences driven mainly by species turnover rather than nestedness. The pronounced differentiation in species-specific environmental responses suggests that multispecies coexistence in the study area is more likely maintained by niche differentiation under environmentally heterogeneous conditions than by fixed community modules or spatial segregation. Accordingly, conservation management should not focus solely on a limited number of richness hotspots or flagship-species core habitats, but should place greater emphasis on maintaining the integrity of environmental gradients, habitat connectivity [57], and microhabitat diversity in order to better encompass broader community patterns and multispecies conservation needs within the reserve. Although human activities did not produce a uniformly strong negative effect at the community level in the present study, they had clear effects on some species, indicating that continued regulation of livestock grazing [58] and control of the cumulative impacts of human disturbance remain necessary.