Multidimensional β-diversity responses to global change: a meta-analysis highlighting divergent effects on plant communities

Literature search

We conducted a literature search in Web of Science, Google Scholar, and the China Knowledge Resource Integrated (CNKI) between 29 March and 31 November 2023. The search strategies for each database were as follows:

• Web of Science: We used the following search query in the Topic field (TS = ): TS= (beta diversity OR β diversity) AND (biotic homogen* OR biotic differentiat* OR biotic heterogen* OR biological homogen* OR biological differentiat* OR biological heterogen*).

• Google Scholar: The same keyword combination was entered directly into the general search bar: (beta diversity OR β diversity) AND (biotic homogen* OR biotic differentiat* OR biotic heterogen* OR biological homogen* OR biological differentiat* OR biological heterogen*).

• CNKI: We searched in the Topic field using the Chinese equivalents of the above terms, including the phrases corresponding to beta diversity OR β diversity OR biotic homogenization OR biotic heterogenization.

Across the three databases, our search initially retrieved 35,206 records after deduplication within each database (see Supplementary Fig. 1). Since some studies appeared in multiple databases, we subsequently merged the datasets and conducted cross-database deduplication to eliminate overlapping records. In addition, we manually reviewed the reference lists of key review papers and relevant studies that cited them^{16,17,19,35,36,37}, which yielded 229 additional studies. After combining the database records with these citation-tracked studies, we performed a final round of deduplication across all sources. This process resulted in a total of 14,977 unique studies for screening. Only published, peer-reviewed articles were retained for the final analysis.

Literature screening and exclusion criteria

First, all studies were independently screened by two researchers (HLL and YH) based on titles, keywords, and abstracts to identify studies that might focus on the effects of human-induced global change on plant β-diversity, excluding 4766 papers. Then, two researchers (HLL and MN) independently conducted a full-text review of the remaining 1779 articles. Unsuitable publications were excluded based on the following criteria: (1) no consideration of specified GCFs, (2) review papers and studies lacking quantitative analysis, (3) articles predicting changes related to homogenization or differentiation in plant communities composition based on changes in plant α-diversity but not directly analyzing the impact of GCFs on plant β-diversity, (4) articles lacking statistical comparisons between the effects of GCFs and control groups, and (5) articles focusing on intraspecific variation caused by GCFs. (6) To reduce the potential biases in the results, the observations with sample sizes smaller than 3 were also excluded. Any uncertainties or disagreements were resolved through discussion with a third reviewer (PL). The PRISMA flowchart (Supplementary Fig. 1) illustrates the procedures we used to select the studies suitable for analysis.

Finally, our search yielded 1604 test results from 256 studies (Supplementary Data 1) across all continents except Antarctica (Fig. 4A, D). This dataset included information pertaining to TB (1,305 tests from 246 studies), FB (136 tests from 28 studies), and PB (163 tests from 20 studies). The tests conducted in these studies spanned terrestrial and freshwater biomes across six vegetated continents, covering elevations from less than 5 to 4715 m (Fig. 4C) and wide temperature and precipitation gradients (Fig. 4B).

Critical appraisal

To assess the quality of studies included in this review, we applied an adapted version of the framework proposed by Mupepele et al.³⁸, which was specifically designed to evaluate the strength of evidence in studies related to ecosystem services and biodiversity conservation. This framework considers aspects such as study design, data collection, analytical methods, and reporting quality (Supplementary Table 1). In selecting and adapting the specific criteria for our study, we also referred to the approaches used by Huynh et al.³⁹ and Su et al.⁴⁰.

Each included study was critically appraised using this framework, and assigned a quality score. Based on the total score, studies were categorized into four levels of evidence: very strong (>75%), strong (50–74%), moderate (25–49%), and weak (<24%). All studies included in the meta-analysis were rated as ‘very strong’ evidence (quality score >75%) based on our appraisal framework (see Supplementary Data 2 for details).

Data extraction

We independently extracted measures of central tendency (mean and median), dispersion (standard deviations, standard error, confidence intervals, and amplitude), and sample sizes for both the control (sites not exposed to GCFs) and statistical treatment (sites affected by GCFs) groups in the selected studies by two groups (HLL and WWX combined as group one, and YC and YH combined as group two). In cases where studies did not report outcomes correctly or lacked sufficient original data, they were excluded from our analysis. In the case of studies that did not provide means and standard deviations, we estimated these values from medians and interquartile range values and amplitudes based on Hozo et al.⁴¹ and Wan et al.⁴² or from standard errors and confidence intervals as described in Lajeunesse⁴³. We used WebPlotDigitizer (v.4.7) (https://apps.automeris.io/wpd/index.zh_CN.html) to extract values from the graphs in the selected studies. If a study involved more than two levels of GCFs, we extracted data for as many levels as possible.

We also collected data from each study for the following seven categories: (1) GCF types, which included biological invasion (BI), climate change (CC), land use change (LUC), land management (LM), and multiple factors (MF); (2) biodiversity dimension, encompassing TB, FB, and PB; (3) biogeographic context, including latitude, longitude, region, realm, elevation, and whether the study site was within a biodiversity hotspot; (4) ecological and climatic context, covering climatic zone, ecosystem, flora realm, mean annual precipitation (MAP, mm), and mean annual temperature (MAT, °C); (5) study design and sampling framework, specifying whether the study focused on temporal or spatial design, sample size, control sample type, temporal extent, and spatial extent; (6) data structure, including data sources and whether the data were abundance-weighted; and (7) focal taxa and functional groups, which included the examined plant taxa, plant life forms, data sources, and whether the data were abundance-weighted. If the referenced study did not report the latitude and longitude for the study area, the approximate coordinates were derived by geocoding the site name in Google Earth 7.0 (free version). If the referenced study did not report MAT, MAP, or elevation, the values were derived from WorldClim using the site’s geographic location (i.e., latitude and longitude). Because some studies covered wide areas (e.g., entire countries or even the entire globe) and did not include specific latitude and longitude coordinates, they were not assigned environmental metrics using the above-mentioned methods, accounting for 6% of the total.

In our classification, “CC” encompassed events such as wildfires, temperature increases, and broader climatic shifts; “LUC” referred to the effects of fragmentation, urbanization, habitat loss, road construction, the erection of wind towers, and mining activities; “LM” pertained to a range of human actions including controlled burning, ecological restoration, logging operations, grazing practices, plantations, comprehensive forest management, and agricultural activities; and finally, “MF” denoted instances where more than two human activities were simultaneously present. In this study, all these factors were referred to as GCFs.

Effect size calculation

Beta diversity is broadly characterized as the ratio (or difference) between gamma diversity and alpha diversity, or the variation (or dissimilarity) in composition among assemblages within a defined spatial area¹⁹. Correspondingly, the methods used for the calculation of β-diversity can be categorized into three main groups: multiplication (β = γ/\(\bar{a}\)), addition (β = γ − \(\bar{a}\)), and dissimilarity measures (such as the Jaccard or Sørensen index for presence/absence data, or percentage difference, such as the Bray–Curtis index of dissimilarity for variations in the relative abundances of species).

We estimated the log-response ratio (ln RR) as a measure of effect size. This parameter is the natural logarithm of the ratio of mean β-diversity measured during statistical treatment (numerator) to that of the control (denominator)⁴⁴. To minimize bias in effect size estimation and variance calculation, we adopted the multivariate Delta method proposed by Lajeunesse⁴⁵, which addresses the small-sample bias associated with conventional log response ratio (ln RR) estimators. Specifically, we used a second-order Taylor series expansion to adjust both the mean (Eq. 3) and sampling variance (Eq. 4) of ln RR, enhancing the robustness and accuracy of effect size estimates, particularly under the small-sample conditions common in ecological research. Negative and positive effect sizes indicated that GCFs decreased and increased β-diversity, respectively (denoting biotic homogenization and biotic differentiation, respectively). Equations are defined as follows:

where Xt and Xc represent the mean plant diversity under GCFs and in the control, respectively.

Furthermore, ordination analysis, which includes principal component analysis, redundancy analysis, correspondence analysis, principal coordinate analysis, and nonmetric multidimensional scaling, serves as a pivotal method for analyzing ecological community data⁴⁶. In our dataset, there are 47 studies that used ordination plots to visually display the effects of human activities on β-diversity. To include the data from these 47 studies in our meta-analysis, we followed the method described by Zhou et al.⁴⁷ to convert the community data from two-dimensional ordination plots to one-dimensional data in order to extract the raw data from these articles. Initially, we extracted the sample positions on the first two ordination axes. Then, we computed the Euclidean distances among different samples using the packages vegan⁴⁸ v.2.6.8 in R, obtaining the distances within the control (Dc) and treatment (Dt) groups, and calculated the means, standard deviations, and sample sizes of both groups. Finally, we determined the RR of β-diversity (RR_Beta) using Eq. (5) as follow:

where \(\overline{{Dc}}\) and \(\overline{{Dt}}\) are the means for the Dc and Dt groups, respectively. Their corresponding variances were calculated as Eq. (5). Overall, ln RR_Beta values < 0 indicated that GCFs decreased β-diversity, while values > 0 indicated that they increased this parameter. Its variances (\(v\)) were calculated as:

where n_t and n_c are the sample sizes of the variable in the treatment and control, respectively; s_t and s_c are the standard deviations of the variable in the treatment and control, respectively. The percentage change in plant β-diversity induced by GCFs was calculated using the equation (e ^{ln RR}−1) × 100%.

According to the test results, the response of plant β-diversity to GCFs did not differ significantly between studies that used different computational methods (Supplementary Fig. 2), suggesting that studies that adopt such methodological approach can be used to assess the response of plant β-diversity to GCFs.

Statistical analysis

Testing the overall impact of GCFs on plant β-diversity

To evaluate whether GCFs generally drive biotic homogenization in plant communities, we tested the overall direction of plant β-diversity responses to GCFs using a mixed-effects model with restricted maximum likelihood (REML) estimation to compute the weighted mean of ln RR, implemented via the rma.mv() function in the metafor⁴⁹ package (v.4.6.0). Since some studies contributed more than one effect size and there were cases where multiple datasets shared a single control group, we conducted a hierarchical meta-analysis to account for the non-independence of our data. Specifically, we treated Control ID (nested within Study ID) as a random effect to control for the potential correlation of effect sizes within the same study or control group. The random effect of Study ID accounted for potential autocorrelation arising from regional environmental conditions, study-specific methodologies, and biases. The model is as follows:

Where Effect Size_ij represents the effect size for the i-th observation within the j-th study. β₀ is the fixed effect (the total effect of GCF on plant beta diversity). (1∣Study ID) accounts for variation between studies. (1∣Study ID: Control ID) accounts for variation due to control groups nested within each study. \({\epsilon }_{i,j}\) is the residual error term.

We computed the overall mean effect size across all studies included in the meta-analysis. To assess heterogeneity among effect sizes, we reported both the total Q statistic (Q_t) and the I² statistic. The Q_t statistic evaluates whether the variability in observed effect sizes exceeds that expected by chance alone, with a significant result indicating the presence of true heterogeneity. The I² statistic quantifies the proportion of total variation attributable to true heterogeneity rather than sampling error, with values of 25%, 50%, and 75% commonly interpreted as low, moderate, and high heterogeneity, respectively⁵⁰. The overall model showed a high total heterogeneity (Q_t = 308760, p < 0.0001, I² = 99.93%), indicating that nearly all differences in the observed effect were attributable to between-study differences in intervention outcomes.

Consistency analysis of GCFs effects on TB, FB, and PB

To assess the effects of specific GCFs on the three dimensions of plant β-diversity, we performed a two-level subgroup analysis. At the first level, we categorized the data by diversity dimension (taxonomic, functional, and phylogenetic) and GCF type (biological invasion, climate change, land-use change, land management, and multiple factors), resulting in eight subgroups. The second level further combined these categories, yielding 15 unique subgroups to evaluate the interactive effects of GCF type and diversity dimension.

To further explore the consistency of responses in TB, FB, and PB to GCFs, we conducted multilevel bivariate meta-analyses. These analyses focused on studies that included two or more biodiversity dimensions simultaneously. The model accounted for the nested structure of the data by incorporating random effects for observation identity (~1|Obser_ID) and study-level variation across dimensions (~Dimension|Study_ID), using an unstructured variance-covariance matrix.

Identifying drivers of variation in β-diversity responses

To evaluate the relative contribution of potential moderators—including GCF type, biodiversity dimension, geographic location, environmental conditions, data characteristics, research subject, and study design—we employed a random forest approach. Prior to model fitting, we assessed the extent of missing data across moderator variables. The variable “temporal extent” was excluded due to >50% missingness. The remaining variables were imputed using the missForest⁵¹ package (v1.5), a non-parametric random forest-based imputation method. Imputation accuracy was high (NRMSE = 0.1192; PFC = 0), and similarly acceptable across all GCF-specific subsets (BI: 0.1557; CC: 0.0013; LUC: 0.1369; LM: 0.1112; MF: 0.1187; all PFC = 0), satisfying the recommended thresholds (NRMSE < 0.2, PFC < 0.1^51,52).

Random forest models were fitted using the randomForest() function from the randomForest⁵³ package (v4.7-1.2) in R. To identify the most influential moderators and assess statistical significance, we applied two-sided permutation-based significance testing using the rfPermute⁵⁴ package (v2.5.2) with 1,000 iterations. This method evaluates the importance of each predictor by assessing its contribution to model accuracy through random permutations, with false discovery rate (FDR) correction applied for multiple comparisons. Model performance was evaluated on the independent test set using R² (calculated as 1 – SS_residual/SS_total) and RMSE.

To further investigate the influence of continuous moderators on effect sizes, we conducted additional meta-regression analyses using the rma.mv() function from the metafor⁴⁹ package (v4.6.0). The models were specified with a nested random-effects structure (Study_ID/Control_ID) and fitted using restricted maximum likelihood (REML). To account for potential violations of model assumptions, we applied cluster-robust standard errors with CR2 adjustment—an approach recommended for enhancing the reliability of meta-regression results under conditions where standard assumptions may not be fully met^55,56.

Publication bias

To assess potential publication bias, we conducted Egger’s regression test⁵⁷ to detect funnel plot asymmetry and applied the trim-and-fill method⁵⁸ to estimate and adjust for potentially missing studies. Although small publication bias was observed for some variables, the corrected effect sizes were largely consistent with the original values and retained the same direction of effect (Supplementary Tables 2 and 3). However, there were exceptions in the global-scale effects of LUC and LM, as well as the response of TB to GCFs. This, along with the high heterogeneity observed in the initial subgroup analysis, prompted us to conduct a two-level subgroup analysis, which showed that our results were not influenced by publication bias (Supplementary Tables 2 and 3).

All statistical analyses were conducted in R v.4.4.1⁵⁹.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.