Global patterns and drivers of spatial autocorrelation in plant communities in protected areas

Stochastic and deterministic (biotic factors such as intraspecies competition or abiotic factors) processes affect the spatial patterns of ecological communities. These processes can be quantitatively assessed based on spatial autocorrelation (SAC) parameters. However, the global patterns of SAC and their differences across vegetation types remain unknown. We aimed to (1) quantitatively assess the SAC global pattern in plant communities and their variation from 1990 to 2020 and (2) identify the key drivers of SAC in plant communities in protected areas (PAs).


| INTRODUC TI ON
Species components and their proportions in plant communities are controlled by deterministic and stochastic factors (Brownstein et al., 2012;Hubbell, 2005).Deterministic processes, which emphasize the roles of biotic and environmental factors, have been extensively documented.By contrast, stochastic processes, which make the species composition of a community unpredictable due to biotic and environmental factors, have been rarely studied.However, the stochastic components of community structure remain difficult to assess despite being extensively studied, particularly in terms of neutral theory.Such stochasticity is caused by several factors, such as unpredictable dispersal by air or water, unpredictable dispersal by animals, unpredictable response to environmental variations, occasional disturbance, establishment and death rate of individual plants, and inter-or intraspecies competition (Brownstein et al., 2012).
Spatially autocorrelated vegetation would have an associated dissimogram that rises quickly to an asymptote that represents the highest average dissimilarity at the longest distances within the extent of the sample site (Lawrence Lodge et al., 2007).In a fitted curve of pairwise dissimilarity against the distance measuring the spatial autocorrelation (SAC) of plant communities, the y-intercept of the curve, i.e., the Nugget, where the distance is zero, all dispersal limitations and spatially correlated environmental differences are excluded, and only stochasticity remains; it represents the contribution of stochastic processes in controlling species composition (Brownstein et al., 2012;Scheller & Mladenoff, 2002).In this curve, the Partial Sill (PSill) represents the contribution of deterministic processes in plant communities.A small PSill value indicates relatively high similarity across plant communities along the environmental gradient, suggesting strong deterministic control.The Sill, which is the sum of the Nugget and PSill, represents the total variation in plant composition determined by stochastic and deterministic processes and reflects the spatial heterogeneity across plant communities.On a local scale, the Sill is affected by several non-mutually exclusive processes, including dispersal limitation, local stochastic effects, and species interactions (Morlon et al., 2008;Wang et al., 2011Wang et al., , 2015)).Consequently, SAC parameters, which indicate the contributions of stochastic and deterministic factors in plant assembly, influence plant community assembly mechanisms and biodiversity stability (Condit et al., 2002;Morlon et al., 2008;Wang et al., 2011).These parameters can also be applied to manage ecosystems (Biswas et al., 2017).SAC quantification is crucial to elucidate the mechanisms by which plant diversity patterns form in plant communities (Biswas et al., 2017).However, most studies exploring the drivers of plant community assembly have excluded the contribution of SAC.
The global patterns of SAC and their differences across vegetation types remain unknown.
Global warming has continued for nearly a century, and the rate of temperature rise in the 21st century is about double that in the 20th century (Karl et al., 2015), especially in high-latitude and high-altitude areas (Peñuelas et al., 2013).Many studies have analysed the effects of climate change on vegetation composition in a given region.A previous study demonstrated that species composition shifts with increasing temperature and cumulative N deposition (Boutin et al., 2017).Generally, the relationship between wood species and climate is highly consistent across spatial scales and organizational levels (Harrison et al., 2020;Kreft & Jetz, 2007).
Using a very large dataset of 6 million trees in more than 180,000 field plots in central Africa, a previous study demonstrated that dry forests are highly sensitive to climate change (Réjou-Méchain et al., 2021).Similarly, recent studies have consistently shown that compared with tropical wet forests, tropical dry forests exhibit larger functional variations in response to long-term drought in West Africa (Aguirre-Gutiérrez et al., 2019) and are likely to be more sensitive to global warming (Sullivan et al., 2020).A recent study in Amazonia has also found a peak in phylogenetic diversity at an intermediate level of precipitation, at which dry-and wetadapted lineages are mixed (Neves et al., 2020).By contrast, forests dominated by widespread tree taxa adapted to anthropogenic pressure show a relatively low sensitivity to climate change (Réjou-Méchain et al., 2021).Model prediction indicated that undisturbed semi-deciduous and transitional forests are phylogenetically more diverse and thus may be less sensitive to climate change than evergreen forests (Réjou-Méchain et al., 2021).However, these intensive studies have mainly focused on the effects of climate change on plant communities in a certain region, and studies on the global patterns of SAC across vegetation types are lacking.In addition, the key climatic drivers of SAC variations in plant communities remain unknown.
Field measurements of community composition and its SAC on a global scale are expensive and labour-intensive.For example, the similarity distance in plant communities between pairwise sites may exceed 10 km (Ma et al., 2021).Remote sensing techniques allow the simultaneous assessment of SAC in plant communities in large areas.Spectral variation hypothesis states that different plant species reflect different intensities of electromagnetic radiation; thus, spectral metrics can be used to assess the composition of and identify plant species (Cavender-Bares et al., 2017;Jetz et al., 2016;Oldeland et al., 2010).Many spectral indices indicate spectral species, but the most widely used is the normalized difference vegetation index (NDVI) (Gholizadeh et al., 2019;Hauser et al., 2021;Peng et al., 2019;Perrone et al., 2023;Pouteau et al., 2018;Somers et al., 2015;Tan et al., 2023).Thus, in the present study, we used this index to simulate the SAC patterns of plant communities.The rationale behind the use of this index is that in a control.An increase in the SAC in pairwise NDVI would mean that the number of species or evenness of species abundance distributions is similar among spatially adjacent locations, and vice versa (Brownstein et al., 2012).
In addition, they are ideal sites for evaluating the contributions of stochastic and deterministic processes in plant communities under the impact of climate change.In PAs, human disturbance is strictly restricted, thereby allowing the patterns and performance of stochastic and deterministic processes and their corresponding climatic drivers to be quantitatively explored.However, these dimensions are poorly understood in plant communities in PAs.Thus, using the same scale and spatial grain, we applied the SAC method to estimate the contributions of stochastic and deterministic processes in controlling plant assembly across vegetation types, including forests (such as deciduous broadleaf forests and deciduous needleleaved forests), shrubs, and grasslands varying in latitude, soil type, slope, aspect, and climate in PAs.This method was used to render the study results comparable across sites on Earth.We hypothesized that (1) the amount (3) larger than 20 square kilometres in the size, so we have enough amount in SAC analysis.These methods can guarantee the quality of Landsat images, sufficient space for plot sampling and SAC analysis, and the reliable comparisons across different vegetation types.We used ArcMap v.10.5 (Environmental Systems Research Institute, 2019) to select PAs, create buffer zones and core areas, and exclude non-vegetated PAs.From PAs that met the above criteria, 49 protected areas with a median size of 20 square kilometres were selected (Figure 1), representing a geographically stratified and broad selection of deciduous broadleaved forest (DBF), deciduous needleleaved forest (DNF), evergreen broadleaved forest (EBF), evergreen needleleaved forest (ENF), shrubs, and grasslands (see details in Appendix S1, Table S1).

| NDVIs derived from Landsat images
Prior to extraction of spectral diversity indices, we pre-processed the Landsat images.Cloud-free Landsat satellite images with a spatial resolution of 30 m obtained at 1990, 2000, 2010, and 2020 for selected 49 PAs were downloaded from the Global Land-cover Facility site hosted by the University of Maryland (glcfa pp.umiacs.umd.edu).All Landsat images were radiometrically and atmospherically corrected using IDRISI GIS (Levin et al., 2007).We corrected the images for shading effects caused by topography using Aster GDEM (Global Digital Elevation F I G U R E 1 Geographical distribution of 147 samples in 49 protected areas in forests, shrubs, and grassland. Model, 30 m, http:// gdem.ersdac.jspac esyst ems.or.jp).All images are from the growing season, when it is easier to separate non-vegetation features from the PAs and minimize problems of image incompatibility due to seasonal and annual differences.
Radiance values were converted to surface reflectance in order to account for differences in exoatmospheric irradiance and solar zenith angles, which reduces spectral variability between images acquired at different dates (Duro et al., 2014;Rocchini, 2007).All image processing was carried out using ERDAS Imagine software.
After pre-processing, the value of NDVI was calculated, and the circle areas with 10-kilometre radius and NDVI > 0.2 were derived as the samples.There are three samples for every PA representing triplicate.There are totally 147 samples (circles) from 49 PAs were selected for further analysis.

| SAC analysis
Semivariograms were produced based on NDVI dataset to evaluate SAC of plant communities at each sample within every PA.
Semivariance is a measure of dissimilarity, which therefore increases with increasing distance; if the semivariance curve is flat rather than increasing trend, this indicates insignificant spatial autocorrelation or that the data are randomly distributed in space (Lee et al., 2011).Semivariance increases over distance because in nature, dissimilarity tends to increase as the distance between sampled objects increases (Lee et al., 2011).Beyond a certain distance called the range (m), the semivariance tends to flatten out, indicating that the variable is spatially independent beyond that distance range; at distances less than the range, the data are considered to be autocorrelated.Nugget, Partial Sill (PSill), and Sill in every semivariogram were derived, which indicates the amount of random (chance or stochasticity), deterministic process (mainly climatic-driven process in the current study), and the variability of spectral species in plant communities, respectively.For most areas, the exact spatial scale of SAC was uncertain (Fahrig, 2013).In a study focused on the determinants of vascular species richness in Finland meadow, 1-2 km radius buffer zones around each study site were selected (Raatikainen et al., 2018).Based on these tests and SAC in various plant communities in the literatures (Engen et al., 2011;Qian & Ricklefs, 2012;Wang et al., 2018), we chose 100 m as the length of step and 100-200 as the account of steps to test ideal SAC models, which represent the possible maximum of range of autocorrelation within each sample.We used these selected distances correspond to reasonable dispersal distances for plant species living on the corresponding habitats (Cousins & Lindborg, 2008;Helm et al., 2005).Therefore, in each semivariogram, there were at least 200 pairs of data points per distance interval.To determine the range of SAC, we created semivariograms of model residuals via random cross-validation and leave-one-out cross-validation tests.Semivariograms, ordinary block Kriging, and contour plots were generated with GS + VERSION 9.0 (Gamma Design Software, MI, USA).

| Importance of environmental drivers
Random forests were used to quantify the relative importance of each of the 132 environmental variables to explain the variation in SAC.Models were fitted using the 'randomForest' package in R (Liaw & Wiener, 2002;R Core Team, 2019).Models were fitted at a global scale for the SAC variables for the four forest types (DBF, DNF, EBF, and ENF), shrubs, grasslands, and the total vegetation (i.e., the seven vegetations combined) in relation to the suite of environmental covariates detailed above.Random forest models also employed in above seven vegetation types in the Northern Hemisphere, since monthly climate variables are comparable in a same hemisphere.
We initially fitted a random forest with n set to 200 and m set to 1.We used the coefficient of determination (R 2 , from random tenfold cross-validation) to assess model performance across multiple trees and repeated iterative amount.A further 200 trees were then added to the model, and the R 2 was recalculated on the evaluation data.If the additional trees improved model performance by more than 1%, they were accepted.This was repeated iteratively until model performance was not improved by additional trees.We then repeated this process with m set to 2 and 3.This process of fitting models with varying values of nt and m was repeated 10 times.The values of m and nt that maximized R 2 across the 10 iterations were then used to fit the final set of models.This resulted in seven random forest models of SAC variables for global vegetations and seven models for Northern Hemisphere.The first 12 environmental variables with the highest importance values in the models were regarded as the key drivers of SAC distribution.

| The relationship of environmental drivers with changes in SAC
To further identify the influence of individual environmental variables on SAC parameters, we used Pearson correlation coefficients.
Temporal changes of SAC and environmental variables were calculated as Δvariable over 5 years.The dependent variable ΔSAC was calculated from the difference in SAC variables at 5-year interval during 1990-2020.The average SAC changes (ΔSAC) were related to the change in key environmental variables (selected by above random forest models) for the same period across all vegetation types, and every vegetation type, by implementing linear regressions using the lm function in R (R Core Team, 2019).

| Spatial and temporal shift in SAC variables
There was considerable variation in SAC parameters from 1990 to 2020 within each vegetation type (Figure 2a-c), which can rival the amount of variation observed across types.Nugget has decreased averagely by 92.36% since 1990 across vegetation types, 89.68% for forest, 99.63% for shrub, and 93.36% for grass.Within forests, the highest variations are found in deciduous needleleaved forest regions (by 98.28%, Figure 2a).From low to high latitude, Nugget value deceased mostly around 39° across all vegetation types, more obvious for grassland (Figure S1), whilst those of deciduous broadleaved or deciduous needleleaved forest have decreased mostly around 3° and 42° (Figure S1).Decreased rate of Partial Sill is 58.36% across all vegetation types (Figure 1b), with the most in grasslands (79.42%) and the least in evergreen broadleaved forest (0.22%), which decreased around 31° and increased around 4.5° (Figure S2).Decreased rate of Sill is 60.61% across all vegetation types (Figure 1c), with the most in grasslands (80.42%) and the least in evergreen broadleaved forest (4.47%), which decreased around 66°and 34° and increased around 53° (Figure S2).The performance of PSill and Sill is similar.That decreased values exhibit the same trends across vegetation types is consistent with the hypothesis that plant communities are becoming similar under climate change along latitude or geographic gradients.
Across the globe, we found that vegetations exhibit significant and substantial shifts in SAC from 1990 to 2020 along latitude (Appendix S2, Figures S1-S3).These shifts represent changes in the shapes of SAC distributions that directly correspond to distinct aspects of plant composition and diversity.All SAC shifts from 1990 to 2020 exhibit obvious latitude-oriented trends (Figures S1-S3).
Along absolute latitude gradient, the decrease in Nugget of all plant communities demonstrates a single peak around the latitude from 12 to 50°.In addition, those of EBF and DBF demonstrated another peaks around 2-4°.Grasslands only have one peak around 35-43° (Figure S1).These trends suggest that the amount of random process in plant communities has strongly decreased in the most recent 40 years, with uneven distribution across vegetation types and disproportionately to the latitudinal gradient.
For PSill, which is believed reflect the change in the amount of climatic-driven process in plant communities, we demonstrated both increase and decrease shifts (Figure S2), with more increased shifts in evergreen forests (EBF and ENF).For deciduous forests, decreased trends dominant in all samples.However, those of grasslands were mainly decreased, indicating spatial similarity in plant communities became stronger from 1990 to 2020 (Figure S2).In total, the shift trends in Sill were same as those of PSill, demonstrating that distance between dissimilar communities is generally decreased.These latitude-oriented changes in plant communities suggest that evergreen forests are less changeable under climate change than deciduous forests, and grasslands demonstrate high sensitivity to climate change, and they all fluctuated along latitude gradient, consistent with the 'favourability hypothesis' (Swenson et al., 2011;Wieczynski et al., 2019), stating that environmental filtering becomes more pronounced as environmental conditions become less favourable.Taken together, these results demonstrate the shifts in random and climatic-driven process in plant communities are uneven across vegetation types, and unproportionable with latitude gradient.Between vegetation types, environmental drivers can explain 92.46% of the variation in PSill in grassland communities, followed by DNF (81.02%) and DBF (67.61%).The explain amount of the variation in Sill for grasslands is 90.68% and that of DNF and DBF is 81.90% and 67.98%, respectively.However, the explain amount of the variation in Nugget is only 56.19% for grasslands, 82.07%for DNF, and 69.64% for ENF (Figure 3).It is noticeable that annual temperature and precipitation are not the dominant drivers for all vegetation types.

| The environmental drivers of SAC parameters
Monthly climate variables may provide more details on environmental drivers for the variation in SAC parameters.Since the months between Northern and Southern Hemispheres are not matched, we used the samples in Northern Hemisphere to identify the contribution of monthly climate variables to the variation in SAC parameters.
The results indicate that monthly variables have obviously improved the explain ability for the variation in SAC parameters (Figure 4).
For Nugget, 82.98% in the variation can be explained by environmental monthly drivers, whilst it is only 48.04% using models based on yearly climatic variables.For PSill, it increased from 69.36% to 80.22%; for Sill, it increased from 72.68% to 82.51%.The explain amount of variation in SAC parameters for every vegetation type also has been considerably increased based on environmental monthly variables.Over vegetation types, those of grasslands can be explained mostly (>91%), and shrubs is the least (<50%).NDVI-CV and WET7-WET9 dominant the drivers of SAC distribution, as well as CLD7-9 and drt3-4 dominant for Nugget, and PRE6, PRE8, and DTR12 dominant for PSill and Sill (Figure 4).For Nugget of grasslands, the key drivers are WET3, VAP11, and DTR4, whilst for ENF, the key drivers are CLD12, CLD8, CLD9, and CV.For EBF, they are CLD7 and CLD9.It is noticeable to point out that annual temperature and precipitation are not the dominant drivers for all vegetation types, only monthly precipitation matters in controlling PSill and Sill across vegetation types.

| Relationship between environmental variables and SAC parameters
NDVI-CV demonstrates significant positive relationship (p < .05)with nearly all SAC parameters across vegetation types, except of Nugget in EBF (Figure 3).Several climate variables administrate consistent relationship with SAC parameters across vegetation types.VAP, TMX, TMP, and TMN are significantly negatively related to all SAC parameters in EBF, whilst FRS positively related to SAC in EBF.
Wet has a significant positive relationship with SAC of grasslands.
PSill and Sill demonstrate nearly the same relationship with environmental variables across vegetation types.
When monthly climatic variables are included, they exhibit stronger relationships with SAC than yearly ones.For example, although SAC parameters of EBF are related to yearly VAP, they are actually strongly related to VAP in July, August, and December when applying monthly climatic variables, not a yearly VAP as a whole (Figure 4).However, we also find interactions between these drivers can obscure the simple directional effects (Figure 5, Figure S4), may lead to discrepancy in interactions between temperature and precipitation.In areas where annual temperature ranged from −10 to −4°C, precipitation exhibited a positive effect on SAC parameters across all vegetation types (Figure 5, Figure S4), probably because in cold areas, the chance of random distribution of plant communities has been increased by increased water supply.In contrast, in areas with annual temperature ranges from 14 to 30°C, precipitation demonstrated a significant negative effect, suggesting increased water supply has enhanced the functional similarity across plant communities.Effects of temperature also have been interacted by precipitation, temperature positively related to SAC in the areas with annual precipitation ranges from 200 to 400 mm, although they negatively associated in other areas.It is noticeable that SAC sharply decreased with increased temperature in the areas where annual precipitation is more than 1600 mm.These strong interactions between temperature and precipitation probably result from geographic variation in the interplay among water, temperature, and local conditions.

| Climatic change and SAC shift
We also examine whether climate change (Figures S5 and S6) is driving changes in SAC over time.The SAC changes (ΔSAC) were related to the change in selected key climate variables (ΔWET, ΔVAP, and ΔPRE) for the same period across deciduous needleleaved forests, shrubs, and grasslands, respectively, by implementing linear regression model (Figure 6, Figure S5).These models demonstrated considerable different relationships within same climate variables at monthly scales.Although the change in annual WET (ΔWET) has a strong negative effect on DNF Nugget changes (Figure 6, Figure S5), its effect is not consistent and significant across all vegetation types (Figure 6, Figure S5) and is stronger positive in grasslands compared to the effect of monthly change.However, change in month variable VAP12, VAP in December has a strong negative effect on ΔPSill of shrubs, is stronger than that of yearly VAP.This effect is also not consistent across vegetation types.For ΔSill of DNF and grasslands,  (Dembicz et al., 2021).Using Landsat images from1982 to 2017, the evolution of month-and annual-climate variables and spectral plant diversity indices in four grasslands in China indicated that the beta diversity has sharply decreasing (Bai et al., 2021).At global scale, an obvious biotic homogenization in plant communities was also predicted in the future decades in various types of grasslands (Peng et al., 2022).
Other vegetation types showed low SAC, with no clear distinctions among different vegetation types.The lowest SAC was found in EBFs and DNFs, which are located around the equator and near the Arctic Circle, respectively.This SAC pattern confirms that plant communities under strong environmental control become similar as they apply the same adaptive mechanisms in response to environmental stress (Wieczynski et al., 2019).Similar biological homogenization was observed by other studies in the same vegetation types as ours.In Tornionjoen area, which near the Arctic Circle, many of the evergreen needleleaved trees likely possess high persistence potential to climate change, subsequently dominant in the changed plant community and contributed to biological homogenization (Bisbing et al., 2021).Under climate change, the plant communities in Taman Negara (around the equator) iteratively maximize homogeneity in species composition based on models considering all environmental variables simultaneously (Macdonald et al., 2019).

F I G U R E 5
Relationships between climate variables and SAC parameters: Nugget, PSill, and Sill.The effects of annual mean temperature (TMP, above) and annual precipitation (PRE, below) on SAC parameters along different temperature and precipitation gradients, respectively.For above three panels, SAC is shown in response to annual mean temperature along six mean annual precipitation gradients (90-200, 200-400, 400-800, 800-1600, and 1600-3000 mm); for below three panels, SAC is shown in response to annual precipitation along four annual mean temperature gradients (−10 to −4, −4 to 4, 4 to 14, and 14 to 30°C).Trend lines are interpolated from a simple linear regression model and are in many cases significant at p < .05(see Figure S4 for more details).
The results of the present study demonstrated that NDVI-CV, a proxy for plant diversity, contributed to the variation in SAC, with a significantly positive relationship across vegetation types.
These results are also consistent with those of previous studies.Brownstein et al. (2012) suggested that communities with low species richness have fewer stochastic factors affecting the Nugget than those with high species richness; consequently, Nugget values would be lower in communities with low species richness than in communities with high species richness.As the number of species interactions increases, the composition becomes more unpredictable and disorderly, and the component of chance increases.
The same researchers found that randomness is positively correlated with whole-community species richness but not with plot-level richness (Brownstein et al., 2012); that is, chance plays a more important role in communities with high species richness than in those with low species richness.Other studies also obtained a high Nugget value in communities with high local species diversity (Anderson et al., 2004;Lieberman & Lieberman, 2007).
Variations in species richness are likely to show strong and positive SAC (Bivand & Wong, 2018).Coniferous species in cold areas are more sensitive than others to global warming (Dyderski et al., 2018;Fisichelli et al., 2014).
With global warming, cold-adapted species are gradually being replaced by warm-favoured species (Elmendorf et al., 2012;Schwörer et al., 2016).Such biotic homogenization trends have been reported in areas similar to the south-central Alaskan and western Canadian boreal forests (Trugman et al., 2018), southern Canada F I G U R E 6 Relationships between changes in selected climate variables by random forest models and changes in SAC parameters: Nugget, PSill and Sill of DNF (deciduous needleleaved forest), shrubs, and grass (grasslands).The effects of changes in wet day frequency (WET, days), vapour pressure (VAP, hectopascals (hPa)), and precipitation (PRE, millimetres) on SAC parameters changes differs considerably among monthly variables.Trend lines are interpolated from a simple linear regression model and are in many cases significant at p < .05(See Figure S5 for more details).and northern America boreal forests (Fisichelli et al., 2014), Landes de Gascogne PA shrubs (Mora et al., 2014), and Olympic PA DNFs (Albright et al., 2013).
This biotic homogenization might be caused by the high mortality of sensitive species and an increase in the abundance of certain general species under the impact of climate change, which has also been reported in other temperate forests (Bose et al., 2017).Under the impact of climate change, climate-related stress selects against highly specialized species.Specialist species are currently declining worldwide at a higher rate than generalist species owing to environmental changes (Olden et al., 2004;Rooney et al., 2004).This trend has drastic effects on species diversity because the replacement of specialist species by generalist species leads to functional homogenization (Clavel et al., 2010;Olden et al., 2004;Olden & Rooney, 2006).Species pairs with high trait similarity tended to be spatially aggregate across spatial scales, suggesting that climate adaptation can help assemble functionally similar species in plant communities (He et al., 2020).This phenomenon can be easily confirmed using spatially aggregated NDVI data (He & Biswas, 2019).Habitat filtering, dispersal limitations, and interspecific interactions may also synergistically influence spatial community dissimilarities.For example, increased dispersal may homogenize community composition (Myers et al., 2013).The present results demonstrate that the contribution of the Nugget (e.g., random or plant-plant interaction processes) to plant community assembly is weaker than that of the PSill, indicating a climatic-driven process.Thus, patterns of spatial community dissimilarities may reflect the outcomes and interplay of multiple community assembly processes.

| Climatic drivers of SAC variations
Annual climatic factors explained a smaller proportion of SAC variations than monthly climatic variables, suggesting that unbalanced energy-related processes and other seasonal fluctuation factors strongly influence SAC patterns.For example, significant positive and negative relationships were observed between ΔPRE and ΔSill across different vegetation types.As a result, the total relationship was weak for all vegetation types.In addition, these effects were regulated by interactions between multiple climatic conditions.For example, ΔPRE and ΔSill showed a significant positive correlation in areas at 4-14°C, whilst a significant negative correlation in areas at 14-30°C.When many or approximately equal numbers of correlation coefficient values are strongly positive and strongly negative, the total coefficient can become null or weak because positive and negative values cancel each other out.Global annual climatic factors exert a weak influence on SAC for functional diversity across plant formations (Bruelheide et al., 2018).In addition, the effects of temperature also are mediated by other factors as water availability (Quan et al., 2019;Sun et al., 2019), which also has been confirmed by our study.However, studies confirming the importance of monthly climatic effects are still insufficient.Therefore, the present results support that monthly climatic variables can explain SAC variations better than annual ones in local communities in PAs.
With the increasing intensity and frequency of global warming, the proportion of warming-tolerant species also increases, with intolerant species relegated to interspaces between large communities.This process is expected to increase the size and number of homogeneous patches with different species combinations relegated to areas further apart and decelerate the increase in community dissimilarities with distance.Therefore, warming in temperate grasslands may be a deterministic process that overrides other stochastic processes contributing to communities (e.g., dispersal, other types of disturbance, and chance establishment and death) and even other deterministic processes (e.g., assembly rules and species interactions).

| Limitations and implications
Despite including several types of vegetation across all continents and latitudes, our dataset lacked sufficient information about some widely distributed vegetation, such as Australian grasslands and shrubs, owing to the unideal remote sensing data from these areas.
For example, China's nature reserves only capture 13.1% of the entire habitat area for threatened plants (Xu et al., 2017).Therefore, our results related to the SAC patterns of shrubs and grasslands could be further refined in the future by incorporating data from currently missing PAs.Nevertheless, our aim was not to present a complete census of global vegetation regions but to assess their SAC patterns and the corresponding drivers using a representative sample.In this respect, the collection of georeferenced vegetation samples presented in this study encompasses all continents and a wide range of latitudes, representing regions with markedly different biogeographical units and vegetation types.We also noted that using NDVI data extracted from satellite remote sensing datasets might create issues related to species discrimination.In fact, we do not need to discriminate among species when analysing SAC patterns across communities.In metacommunities, relations are determined by the degree to which a single community shares (similarity) or differs (distances) in the characteristics of their species or functional compositions (Büchi & Vuilleumier, 2014).
Dissimilarities in community composition are one of the most fundamental and conspicuous features by which different vegetation types may be distinguished.Traditional estimates of community dissimilarity are based on differences in species incidence or abundance (e.g., Jaccard, Sørensen, and Bray-Curtis dissimilarity indices).Consequently, dissimilarities among communities can be accurately detected through NDVI extraction.This similarity across neighbouring communities demonstrated that the dispersal ability of species was the same as their functional adaptation to environmental stress, considering that they perform similarly along the environmental gradient.Specifically, our results are mostly consistent with the results of previous studies at the same areas, indirectly indicating the remote sensed data are reliable at a considerable degree.
The large unexplained variance at the site level may be ascribed to site-related factors, such as local environmental factors, and vegetation history (Dembicz et al., 2021).It can also be attributed to the fact that not all environmental factors were measured or related to community composition on a nonlinear scale of a type not allowed for by the analysis (Brownstein et al., 2012;Yuan et al., 2019).
Intensive studies have suggested that many drivers of biodiversity cannot be generalized and depend only on the regional context.The dataset lacked information on local environmental properties, such as soil fertility.Hence, temporal variations in vegetation in PAs were used in the analyses.Environmental factors, such as soil, landform, slope, and elevation, were assumed constant, and only the impact of climate change was considered.Thus, this study quantified the different actions of plant community assembly mechanisms in a different manner from and complementary to manipulative experiments.

| CON CLUS ION
The SAC distribution of plant communities in PAs worldwide, explained by annual energy metrics, such as maximum, average, and

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this study.

PE E R R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/ api/ gatew ay/ wos/ peer-review/ 10. 1111/ ddi.13793 .
plant diversity, protected areas, spatial autocorrelation, vegetations | 121 PENG et al. fitted curve of pairwise NDVI dissimilarity against the distance, the y-intercept of the curve, that is, the Nugget, where the distance is zero, represents the contribution of stochastic processes to NDVI composition.In this curve, PSill indicates the amount of NDVI dissimilarities across communities controlled by environmental factors, such as climatic variables.A small PSill value indicates a strong environmental control, given that the NDVI dissimilarities across communities will decrease under a strong environmental of stochastic processes decreased and that of deterministic processes (mainly climate factors) increased in PAs worldwide under the impact of climate change from 1990 to 2020, (2) monthly climatic variables can explain the variations in SAC patterns better than annual climatic variables, and (3) SAC performance and climate-SAC relationships differ considerably across vegetation types. 2 | MATERIAL S AND ME THODS 2.1 | Study areas We selected PAs from the World Database on Protected Areas (WDPA) (https:// www.prote ctedp lanet.net) falling in various International Union for Conservation of Nature (IUCN) reserve classifications.The selection follows three criterions: (1) at least one PA for every vegetation type; (2) having Landsat images in the growing season within the study period (1990-2020) with cloud cover less than 10%; Twelve climatic variables describing average temperature, precipitation, cloud cover, and their monthly variation for 1990, 2000, 2050, and 2070 were extracted at each sample from climate maps.These variables provide an adequate description of local climate for all plots.The average NDVI and its CV of every sample were derived as environmental factors.In total, 132 environmental layers, including 130 climate and 2 NDVI variables, were used as covariates in our analyses.All layers were standardized to 10-kilometre radius circle to match the sample extent.Twelve climatic variables are cloud cover (CLD, percentage (%)), diurnal temperature range (DTR, degrees Celsius), frost day frequency (FRS, days), potential evapotranspiration (PET, millimetres per day), precipitation (PRE, millimetres), daily mean temperature (TMP, degrees Celsius), average daily minimum temperature (TMN, degrees Celsius), average daily maximum temperature (TMX, degrees Celsius), vapour pressure (VAP, hectopascals (hPa)), and wet day frequency (WET, days).Their yearly values are marked by y in the name, and monthly values are marked by the number of the month; for example, FRSy indicates yearly FRSy and FRS5 indicates fry in May.
To generate a mechanistic understanding of the environmental factors driving the global variation in SAC parameters in PAs, we explored the directional effects of environmental annual factors (Figure3).Our analyses show that, across all vegetation types, NDVI-CV can explain 33.67%, 44.22%, and 47.80% of the variation in Nugget, PSill, and Sill, respectively, higher than other environmental factors, with a positive relationship.Climate factors such as DTRy, PREy, CLDy, and FRSy are also the dominant drivers of Nugget across vegetation types, as well as NDVI, PREy, and WETy dominant in PSill, and NDVI, DTRy, and CLDy dominant in Sill variation.
Indeed, monthly climatic variables can explain more than yearly F I G U R E 2 (a-c) Boxplots of observed autocorrelation parameters per vegetation type during 1990-2020.Top graphs indicate the differences in autocorrelation parameters between 1990 and 2020 (a: Nugget, b: Partial Sill, c: Sill).For each variable in each plot, the horizontal line is the average value across the samples, and the box indicates the interquartile range (IQR), points highlight outliers.dbf, deciduous broadleaved forest; dnf, deciduous needle forest; ebf, evergreen broadleaved forest; enf, evergreen needleleaved forest; grass, grasslands; shrub, shrubs.Same as the followings.F I G U R E 3 Average values (top), variation explain amount (%) of environmental variables (middle), and variable importance scores in annual environmental variable from random forest models (below) of SAC parameters: Nugget, PSill, and Sill in different vegetation types.The black columns in below panels indicate importance scores in every environmental variable.Note: Different colours in the below panels indicate Pearson's correlation coefficients between every annual environmental variable and Nugget, PSill, and Sill, respectively.variables in the variation of SAC.More specifically, DTR4, VAP11, and WET6 are among the strongest climatic factors related to mean and variance of Nugget across vegetation types.PSill and Sill exhibit very similar relationship between climate factors, with both shrub and DNF are strongly positively related to VAP7, yet both DBF and EBF are negatively related to VAP7.Likewise, PRE8 is strongly positively related to DNF, EBF, shrubs, and grasslands, and negatively related to DBF.

F
Figure S5), whilst a negative effect in shrubs.In contrast to the positive effect of ΔPRE, precipitation change in June (ΔPRE6) demonstrated a significant negative effect on ΔSill of grasslands.These inconsistences in effects on SAC shifts between annual and monthly climate change probably result from great seasonal and yearly divergence in water supply, temperature, and local conditions across vegetation types.
Vegetation under the impact of climate change in 1990-2020, particularly DNFs, grasses, and shrubs, demonstrated biological homogeneity, as indicated by a sharp decrease in SAC.Similarly, Peng et al. (2022) observed a local homogenization of global plant assemblages due to the impact of climate change in grasslands, and Testolin et al. (2021) highlighted climate change as a driver of biotic homogenization of woody plants in the Atlantic forest.
minimum temperatures, differed across vegetation types.It reached the highest levels in grasslands and ENFs whilst the lowest levels in EBFs and DNFs around the equator and near the Arctic Circle, which are regions under extreme climate conditions.Monthly climatic variables can explain more SAC variations than annual ones, suggesting that the annual variables may neglect the unbalanced distribution within a year and the monthly variables can capture the key information in identifying the forces driving the variation in SAC patterns.From 1990 to 2020, the SAC parameters decreased sharply across vegetation types.On average, the Nugget decreased by 92.36%, the PSill by 58.36%, and the Sill by 60.61%, indicating obvious biological homogeneity.With this trend, the amount of stochastic processes in plant community assembly decreased and that of climatic control processes increased, demonstrating a geographical variation across vegetation types.These changes were strongly related to monthly climatic changes, in particular, precipitation.The change weakly related to annual climatic changes resulting from a total coexistence of significant positive and negative relationships by a single monthly climatic variable.These relationships also differed across vegetation types and SAC parameters.Our results support that the single and joint effects of annual and monthly climatic variations on SAC patterns are overarching in plant communities in PAs worldwide, but the specific mechanisms underlying such effects remain to be elucidated.Therefore, future research on research should consider local information about soil biotic and abiotic composition, topographical features, and microclimatic variation on a regional scale.Furthermore, efforts should be directed toward the extraction of plant community data from underrepresented regions.Indeed, this work could serve as a basis for defining global inhomogeneous trends of plant communities, even those in PAs, and further investigations should include patterns of endemism, functional variation, and phylogenetic diversity in the field.This information can improve our understanding of SAC patterns and serve as a guide for the effective management of PAs in response to climate change.ACK N O WLE D G E M ENTS This research was funded by the National Natural Science Foundation (32271555) and the National Key Research and Development Program of China (Project No. 2022YFF1303005).