The study was carried out in the Qiyunshan National Nature Reserve (Qiyunshan) located in Jiangxi Province, China (113°54'-114°07' E, 25°41'-25°54' N) (Fig. 1). Officially established as a protected area in 1997, cover a total area of 17,105 ha with a remarkable forest coverage rate of 97.6% (Liu et al. 2010). The Reserve features a humid subtropical monsoon, characterised by a warm climate and abundant rainfall. The average annual precipitation is 1,750 mm and the average temperature is 17℃ (Liu et al. 2010, Qian et al. 2022). The minimum altitude of the Reserve is 300 m and the maximum is 2,061.3 m. Vegetation types include evergreen broadleaved forest, mixed coniferous broadleaved forest and meadow (Huang et al. 2020, Qian et al. 2022). Over the past decade, there have been no notable changes in land use within Qiyunshan, due to the stringent protection measures implemented by the Chinese government. Despite its ecological importance, biodiversity studies on Qiyunshan remain limited. Most of the studies have focused on migratory birds (Huang et al. 2008, Liu et al. 2010), with relatively few studies focusing on mammals (Qian et al. 2022) and butterfly diversity (Huang et al. 2020).

Figure 1.  

Location of the study area.

Bird monitoring in Qiyunshan was conducted annually from 2014 to 2022 using the line transect method (Bibby et al. 1992). Ten transects, each measuring 2-4 km in length, were established, spaced 200-500 m apart from each other. The number of transects, along with their starting and ending points, remained constant for all survey years (Fig. 1). The bird observations were carried out by trained ornithologists, with two survey periods annually, the first from 1 April to 1 May and the second observation period from 1 June to 30 June. The bird surveys were conducted in the morning from 07:00 h to 10:00 h and from 17:00 h to 19:00 h, aligning with bird activities using binoculars. Bird species and the number of individuals were recorded within 50 m on both sides of transects. Surveys were not conducted during unfavourable weather conditions, such as rain or heavy fog. The identification of bird species was based on Zheng (2023). We referred to the List of State Key Protected Wild Animals in China (http://www.forestry.gov.cn) and the International Union for Conservation of Nature (IUCN) Red List of Threatened Species (http://www.iucnredlist.org) to identify the protection levels and threat levels of birds. Only three observers conducted all sampling events to ensure standardisation of accuracy in species detection.

To evaluate the trend of taxonomic diversity, we assessed the adequacy of bird sampling using the 'iNEXT' package in R. This allowed us to generate rarefaction and extrapolation curves, based on sample size and coverage (Chao et al. 2014, Li et al. 2024). A sample coverage exceeding 0.90 indicates sufficient sampling (Li et al. 2024) (Fig. 2). We assessed bird diversity through species richness, Shannon-Wiener Diversity Index and Simpson Index following Hill numbers (Li et al. 2024), where "q = 0" represents species richness, "q = 1" represents the Shannon-Wiener Index reflecting common species and "q = 2" represents the Simpson Index indicating dominant species (Chao et al. 2014, Li et al. 2024).

Figure 2.  

The rarefaction and extrapolation of bird species diversity in Qiyunshan were measured using the Hill number in various years (2014-2022). The solid lines represent rarefaction, while the dashed lines represent extrapolation; shaded regions, 95% confidence intervals; " q = 0" represents species richness; "q = 1" represents the Shannon-Wiener Index reflecting common species; "q = 2" represents the Simpson Index indicating dominant species; SC represents sample coverage. a 2014; b 2015; c 2016; d 2017; e 2018; f 2019; g 2020; h 2021; i 2022.

To quantify functional diversity, we obtained data on six functional traits from Wang et al. (2021) that characterise the niche of the species, offering insights into how species utilise and compete for resources within their habitats (Luck et al. 2012, Li et al. 2024). We used the following traits: body mass, beak length, body length, wing length, feeding guild (insectivores, carnivores, nectarivorous, granivores and omnivores) and nest site (ground, water, scrub, canopy, rock wall). Amongst these traits, body mass, beak length, body length and wing length were continuous variables (Li et al. 2024). The feeding guild and the nest site were categorical variables. Due to computational limitations, we converted categorical variables into ordinal variables. Four functional diversity metrics were calculated: functional richness (FRic), evenness (FEve), divergence (Fdiv) and dispersion (FDis) (Oliveira and dos Anjos 2022). Functional diversity is defined as follows (Mason et al. 2005, Villéger et al. 2008, Mouchet et al. 2010, Laliberté and Legendre 2010, Li et al. 2024): functional richness measures the extent of functional space utilised by a group of species. Functional evenness indicates the uniformity of species abundances across functional space (Li et al. 2024). Functional divergence quantifies the distance between high species abundances and the centre of functional space (Li et al. 2024). The three indices complement each other. Additionally, functional divergence, functional evenness and functional dispersion are not influenced by species richness, making it possible to compare communities with varying species richness without bias (Villéger et al. 2008, Li et al. 2024). Functional dispersion measures the average distance of each species to the centre of all species, taking into account their respective weights (Li et al. 2024). The analysis of functional diversity was conducted in the R "fundiversity" package (R Core Team 2024).

To perform phylogenetic analyses, we used a vector of scientific species names to generate a phylogeny from megatrees (see details in Li (2023) and Li et al. (2024)). Based on the bird species checklist from our survey, we constructed a phylogenetic tree by classifying species into orders, families and genera using supertrees (Li et al. 2024) (Suppl. material 1). We constructed a phylogenetic tree using the BirdTree database (http://birdtree.org) (Jetz et al. 2012, Li 2023, Li et al. 2024). This phylogenetic tree enables the calculation of Faith’s index of phylogenetic diversity (Faith 1992), mean pairwise distance (MPD) and mean nearest taxon distance (MNTD) for analysis of phylogenetic diversity (Webb 2000). To quantify phylogenetic patterns in community structure, we used a null model with an independent swapping algorithm to stochastically generate species richness and occurrence frequencies (Li et al. 2024). Subsequently, mean FD, MPD and MNTD values for the null model are calculated and compared with the observed values (Li et al. 2024). The standard effect size (SES) is calculated using the following formula:

SES=（M_obs - M_null）/ SD_null

where M_obs is the observed value of FD/MPD/MNTD, M_null is the average of the 999 null model FD/MPD/MNTD values generated randomly and SD_null is the standard deviation of the 999 random values (Li et al. 2024). A negative SES.MPD/MNTD indicates a clustered community phylogenetic structure, whereas a positive SES suggests an overdispersion structure (Webb 2000, Li et al. 2024). SES.MPD or MNTD greater than 1.96 (p < 0.05) signifies a significant overdispersion community structure, potentially attributed to competitive exclusion (Li et al. 2024). Conversely, an SES less than -1.96 (p < 0.05) indicates a significantly clustered community structure, likely influenced by environmental filtering. An SES falling within the range of -1.96 to 1.96 (p > 0.05) implies that the community assembly follows a random process (Webb 2000, Li et al. 2024). The analysis was conducted using the R package "picante" (R Core Team 2024).

Beta diversity reflects two distinct ecological processes, spatial species turnover and nestedness of assemblages (Baselga 2010), which result from two antithetic processes, namely species replacement and species loss, respectively (Baselga 2010). This method also allows for the computation of dissimilarity between locations over time 1 and time 2 (Baselga et al. 2015). We computed the dissimilarity of bird composition for Qiyunshan between 2014 and 2022, considering the turnover and nestedness components of temporal change and the sum of both values (overall change) (Baselga 2010). Following this method, we utilised the Sørensen pairwise dissimilarity index for evaluating the overall beta diversity (β_sor) (Zhao et al. 2021). Then, we employed the Simpson dissimilarity index to quantify the turnover (β_sim) and nestedness (β_sne) components (Baselga 2010, Baselga 2012). Pairwise and multiple-site dissimilarity partitioning frameworks were performed in R using the "betapart" package (Baselga and Orme 2012). To strengthen our results more robustly, we used temporal beta diversity (TBI) to measure species composition differentiation through time (Legendre 2019). The temporal beta diversity was divided into losses and gains of species independently for each site to ascertain the direction of change (Legendre 2019, Jiménez et al. 2024). TBI analysis is especially interesting in species-rich communities where we cannot examine the changes in every species individually (Legendre and Condit 2019). We used the "%difference" method to compute the percentage difference index, this index also being known as the Bray-Curtis dissimilarity (Legendre 2019). Temporal beta diversity was carried out in the "adespatial" package in R (R Core Team 2024).

The non-parametric Mann-Kendall test was used to assess the monotonic trend in the time series of bird diversity for alpha diversity metrics. The non-parametric Mann-Kendall test is commonly employed to detect monotonic trends in a series of environmental data. To explore the potential non-linear trends in bird diversity over the past nine years in Qiyunshan, Local Polynomial Regression was fitted to the alpha diversity metrics. The Mann-Kendall test and Local Polynomial Regression were performed in R using the "trend" package (R Core Team 2024). Generalised Dissimilarity Modelling (GDM) was used to estimate beta diversity across multiple sites (Ferrier et al. 2007). GDM is a powerful and unique method for characterising and predicting beta diversity (Mokany et al. 2022), the change in biodiversity over space, time and environmental gradients (Mokany et al. 2022). GDM uses I-spline basis functions to transform each of the predictor variables (Mokany et al. 2022), so that summed absolute difference in the transformed predictor values (Ferrier et al. 2007, Mokany et al. 2022), combined with the model intercept, provide a predicted ecological distance. The negative exponential link function ensures predicted dissimilarities increase and saturate with increases in the predicted ecological distance (Ferrier et al. 2007, Mokany et al. 2022). Since nestedness is not a central issue in our study, we considered the overall beta diversity, turnover components and temporal beta diversity. GDM was carried out in the "gdm" package in R (R Core Team 2024).

Consistent with studies in the tropical forest nature reserve (Stouffer et al. 2020, Oliveira and dos Anjos 2022), our study showed taxonomic changes in bird communities in Qiyunshan over time. We identified taxonomic turnover as the primary mechanism driving annual variations in species composition. However, despite variations in species richness and composition, both FD and PD remained stable throughout the study period. This non-linear pattern highlights the dynamic nature of bird communities, where temporal shifts in species co-existence, resulting in fluctuations in species richness and composition, even in undisturbed nature reserves (Oliveira and dos Anjos 2022). Our findings suggest that protected areas like Qiyunshan have the potential to attract novel bird species that share similar functional and genetic relationships with existing communities. This dynamic contributes to functional redundancy amongst species, ensuring that FD and PD remain largely unaffected despite taxonomic changes (Villéger et al. 2008, Mouchet et al. 2010, Oliveira and dos Anjos 2022). For instance, we documented Asian Barred Owlet (Glaucidium cuculoides) in 2014, 2018 and 2019 and we recorded Collared Owlet (Glaucidium brodiei) from 2018 to 2022. Although these two carnivorous birds were recorded in different years, their similar functional traits and phylogenetic relationship prevented any niche gap. Furthermore, no significant trends were observed in bird taxonomic composition (Suppl. material 2), indicating that the overall taxonomic structure remained stable.

Furthermore, we observed no distinct boundary between the forested areas within and outside Qiyunshan. This lack of separation, combined with birds' strong dispersal ability enables them to effectively colonise and move across reserves, facilitating their widespread distribution (Tanentzap and Lloyd 2017). Surprisingly, we observed that the species loss showed a sharp upward trend. This phenomenon may be attributed to spillover effects, which amplify the conservation benefits by facilitating species dispersion beyond the boundaries of the protected area. For example, the spillover effects of fenced reserves in New Zealand, where avian species breed within the reserve and disperse into surrounding areas demonstrate this process (Tanentzap and Lloyd 2017). Since we do not have data on birds around Qiyunshan, our conclusion should be interpreted cautiously. Policy-makers and conservation authorities should consider these limitations when applying our findings to broader conservation strategies.

Our findings indicate that the patterns observed in all measures of biodiversity are not always consistent. Numerous studies have demonstrated the asynchrony of TD, FD and PD in birds (Villéger et al. 2008, Devictor et al. 2010, Morelli et al. 2017). For example, breeding birds in France from 1989 to 2012 showed noticeable variability both within and amongst the temporal trends of each aspect of bird diversity, revealing a notable increase in TD and PD (Monnet et al. 2014), while FD showed no discernible pattern (Monnet et al. 2014). Similarly, in China, water-bird species richness and SES.MNTD declined from the 1950s to the 2010s without significant trends (Che et al. 2021), whereas SES.MPD and FD experienced a significant decline during the same period (Che et al. 2021). Previous studies have shown that PD is a better measure of biodiversity than taxonomic counts (Miller et al. 2018). PD is particularly valuable when many traits exhibit a phylogenetic signal, indicating that PD can assess the functional trait space of a community (Srivastava et al. 2012). However, the effectiveness of PD relies on the robustness of phylogenetic signals to both traits and interactions (Srivastava et al. 2012). If future studies can integrate these intricacies, there is potential for conceptual unification of evolutionary biology and ecosystem ecology through exploring the relationships between PD and ecosystem function (Srivastava et al. 2012).

Our findings indicate that the composition of bird communities is shaped by both environmental filtering and neutral processes. For community assembly rules, phylogenetic clustering may result from the local extinction of species that are phylogenetically dissimilar to those present or the colonisation of species that are phylogenetically and functionally similar to the residents (Li et al. 2015). However, according to the neutral theory of biodiversity, ecological communities are dynamic and open systems where species co-occurrence is unstructured (Chave 2004, Hubbell 2011, Oliveira and dos Anjos 2022). Communities are randomly structured and independent with an assembly process driven by neutral dispersal, local stochastic extinction and ecological drift, rather than species-specific characteristics (Hubbell 2011, Oliveira and dos Anjos 2022).

Previous studies suggest that PD and MNTD are effective metrics for detecting habitat filtering, while MPD is more suited to identifying competitive exclusion (Miller et al. 2016). For community assembly, our findings are consistent with previous studies in the Ailao Mountains in southwest China (He et al. 2018), Nyungwe National Park in southwest Rwanda (Rurangwa et al. 2021) and Mt. Kenya of Kenya (Gao et al. 2024). Environmental filtering suggests that only species with similar traits can survive and reproduce in a given abiotic environment (Miller et al. 2016). Thus, if these traits are evolutionarily conserved, habitat selection will lead local communities that are more phylogenetically similar than expected by chance, based on phylogenetic clustering (Miller et al. 2016). However, long-term patterns of plant phylogenetic diversity have indicated a decrease in phylogenetic clustering as succession progresses, with colonisation playing a primary role in driving this transition (Li et al. 2015). Early colonists were closely related and functionally similar to resident species, whereas later colonists showed decreasing similarity to the existing species (Li et al. 2015). In this study, the phylogenetic clustering pattern may be influenced by the resource availability of Qiyunshan. The lack of clear patterns of the phylogenetic assemblage structures suggests complex interactions between ecological and evolutionary processes (Gao et al. 2024).

Understanding the temporal dynamics of biodiversity is essential for developing effective environmental management policies in PAs. Our findings indicate a significant increase in temporal beta diversity, with the turnover component being larger than the nestedness component, suggesting potential changes in bird community composition in PAs. Our findings are consistent with previous studies in southwest France (Baselga et al. 2015) and the Atlantic Forest protected area in southern Brazil (Oliveira and dos Anjos 2022). Generally, temporal changes in species composition may result from deterministic processes driven by environmental changes, random processes of colonisation and local extinction (Baselga et al. 2015). No significant land-use changes were observed in Qiyunshan during the study period. Even in anthropogenic landscapes, prior studies have demonstrated that alterations in land cover within agricultural environments have had minimal effects on the temporal beta diversity of bird assemblages (Baselga et al. 2015). More broadly, previous studies suggest that stochastic processes may play a role in explaining temporal changes in biodiversity (Oliveira and dos Anjos 2022) and that bird species appeared and disappeared from specific localities in a random-like way (Baselga et al. 2015). In other words, this phenomenon likely reflects random variation driven by the observed random losses and gains of species (Legendre 2019, Legendre and Condit 2019).