Global impoverishment of natural vegetation revealed by dark diversity

Abstract

Anthropogenic biodiversity decline threatens the functioning of ecosystems and the many benefits they provide to humanity¹. As well as causing species losses in directly affected locations, human influence might also reduce biodiversity in relatively unmodified vegetation if far-reaching anthropogenic effects trigger local extinctions and hinder recolonization. Here we show that local plant diversity is globally negatively related to the level of anthropogenic activity in the surrounding region. Impoverishment of natural vegetation was evident only when we considered community completeness: the proportion of all suitable species in the region that are present at a site. To estimate community completeness, we compared the number of recorded species with the dark diversity-ecologically suitable species that are absent from a site but present in the surrounding region². In the sampled regions with a minimal human footprint index, an average of 35% of suitable plant species were present locally, compared with less than 20% in highly affected regions. Besides having the potential to uncover overlooked threats to biodiversity, dark diversity also provides guidance for nature conservation. Species in the dark diversity remain regionally present, and their local populations might be restored through measures that improve connectivity between natural vegetation fragments and reduce threats to population persistence.

Fig. 1. Estimating dark diversity and related biodiversity metrics in ecological communities.

a, Data included a local study site where certain species were present, but many species sampled elsewhere in the region were absent. To estimate the probability that a species that is absent from the site but present in the region belongs to the dark diversity of the site, we used information about species co-occurrences at other sites in the region. b, We calculated an indicator matrix in which each present species indicated the ecological suitability of each absent species for the study site. We compared the observed number of co-occurrences with the number of co-occurrences expected at random (according to the hypergeometric distribution) and standardized the difference using the standard deviation from the hypergeometric distribution. c, By averaging across all observed species, each absent species was assigned a probability of belonging to the dark diversity for the study site. Consequently, the dark diversity was a fuzzy set to which species belonged to varying degrees. d, Several biodiversity metrics were characterized for each site in the region. Alpha diversity was the number of species recorded at the site, and gamma diversity was the total number of species recorded in a region. The size of dark diversity was estimated as the sum of the probabilities of absent species belonging to the dark diversity of the study site. Alpha and dark diversity together formed the site-specific species pool, and gamma diversity not falling into this category was considered the unsuitable part of gamma diversity; that is, belonging to the species pools of other sites. We investigated the percentage of the species pool that was present among the alpha diversity (community completeness) and the turnover of species pools in the region, expressed as the percentage of gamma diversity that was unsuitable for the study site (beta diversity).

Fig. 2. Plant diversity in natural vegetation in relation to human effects in the surrounding regions.

a, Relationship between community completeness in natural vegetation and the human footprint index in the surrounding area, defined by a radius of 300 km. The prediction line from a multiple linear regression model is shown with the 95% confidence intervals. Note that community completeness values on the y axis are back-transformed from the logit scale. The symbol tones indicate forest cover (0–100%). R² value of the model and two-tailed P value of the relationship are shown; n = 116 regions. The distribution of community completeness is shown in the histogram on the right (median, 25%). b, Left, model summaries linking community completeness to the human footprint index and its components across spatial scales. Human influence was averaged over various spatial scales around the study regions (radii 10 km, 50 km, 100 km, 200 km, 300 km and 400 km), and the respective models were compared using the Akaike information criterion (AIC). Filled symbols indicate significant relationships (P < 0.05), and the large symbol indicates the set of best significant models (ΔAIC < 2). Right, from the best model (the smallest scale at which ΔAIC < 2), the effect of the human footprint index or one of its components is shown as a standardized coefficient (dot) with a 95% confidence interval (CI; line); n = 116 regions. Filled symbols and bold confidence interval lines indicate significant effects. c, Map of sampling regions, with community completeness indicated by symbol size and the underlying map showing the global variation in the human footprint index (the highest value within each grid cell of around 0.25° × 0.25°). The inset shows part of Europe containing a large number of study regions. Triangles indicate regions in which only woody species were sampled. Symbol tones indicate the percentage of forests in regions.

Extended Data Fig. 1. Distribution of the 119 DarkDivNet study regions in relation to mean annual temperature and annual precipitation.

Lines indicate ranges within a radius of 100 km. Approximate broad biomes are shown. Triangles indicate regions in which only woody species were sampled.

Extended Data Fig. 2. Using independent data to test the dark diversity method that relies on species co-occurrences to estimate the ecological suitability of absent species.

We used two tests. In the Vicinity test, we examined whether species absent from the site (100 m²) but present in the immediate vicinity (2500 m²) have higher estimated suitabilities than absent species found further away. The sample area and vicinity area are assumed to share relatively similar ecological conditions. In the Expert test, we compared whether species absent from the site but assessed by expert opinion to be ecologically suitable (i.e., belong to the site-specific species pool) have higher calculated suitabilities than those absent species that were evaluated as unsuitable. In both cases, we calculated the log response ratio of the mean suitability of species in the respective groups. Positive log response ratios indicate agreement between assessments of suitability calculated from co-occurrences and from the independent information considered in the tests. The length of the lines (vertical for the Vicinity test and horizontal for the Expert test) shown at study region locations indicates the magnitude of the log response ratio; negative values are in red and positive values are in blue. Both tests comprised data from a subset of study regions. The box plot on the left (centre line, median; box limits, upper and lower quartiles; whiskers, the range, excluding outlying points that exceed the quartiles by more than 1.5× the interquartile range) shows the results of single-sample two-sided t-tests (difference from zero), with log response ratios significantly larger than zero in both cases, n = 115 regions for the Vicinity test and n = 76 regions for the Expert test. Hexagons on the map (made with Natural Earth; free vector and raster map data; https://www.naturalearthdata.com/) delimit the spatial blocks used in cross-validation.

Extended Data Fig. 3. Bivariate plots of biodiversity metrics derived using alternative sampling methods.

a, Gamma diversity from 30 sites compared with extrapolation up to complete sampling coverage. b–f, Sites described in a 2500 m² area in addition to the DarkDivNet standard of 100 m². g–k, Biodiversity metrics when 60 sites were used to estimate co-occurrences in addition to the DarkDivNet standard of 30. The scatter plots show mean values for regions where the respective sampling scheme was applied (n = 119 regions for a, n = 116 regions for b–f and n = 27 for g–k). The 1:1 lines are shown as diagonals. Estimates of Spearman correlation for each comparison are shown above the panels. Comparisons where the alternative sampling method did not influence the metric (i.e. gamma diversity when using a larger sample area or alpha diversity when using more sites) are not shown.

Extended Data Fig. 4. Relationships between biodiversity metrics and the human footprint index in the surrounding regions.

a, Alpha diversity. b, Beta diversity. c, Gamma diversity. d, Dark diversity. e, Species pool size. A similar graph for community completeness is shown in Fig. 2a. For each metric, the relationships from the spatial scale producing the best multiple linear regression model is shown (n = 116 regions, see Extended Data Table 1 for details of the models). The prediction lines are shown with 95% confidence intervals. The solid line indicates a significant relationship (two-tailed p < 0.05); the dashed lines indicate non-significant trends. Note that the range of the human footprint index varies when averaged at different spatial scales. Diversity values on y-axes are back-transformed from log or logit scales.

Extended Data Fig. 5. Correlations between four principal components and explanatory variables.

Correlations with raw environmental variables (used in the PCA; top) and the human footprint index or its components (not used in the PCA; bottom).

Extended Data Fig. 6. Improvement in multiple linear regression models describing community completeness provided by including quantiles of the human footprint index values found in regions at various spatial scales.

Filled symbols indicate significant relationships (n = 116 regions, two-tailed p < 0.05), and the large symbols indicate models where the Akaike information criterion (AIC) is lower than the minimal value among the models using the mean human footprint index value (Extended Data Table 1). The asterisk indicates a combination of quantile and spatial scale that yielded a considerably better model (AIC value lower by more than 2 units).

Extended Data Fig. 7. Bivariate relationships between diversity metrics.

a-e, Alpha diversity. a, f-i, Dark diversity. b, f, j-l, Species pool size. c, g, j, m-o, Gamma diversity. d, h, k, m, o, Community completeness. e, i, l, n, o, Beta diversity. Lines indicate variation within regions (99% quantiles, i.e. omitting outliers), crossing at median values within regions. Several metrics are inherently related (see Fig. 1), and strong relationships are expected. However, variation within regions can be large, indicating the importance of site-specific metrics. The colours of lines reflect the different biodiversity metrics (see Fig. 1) to facilitate comparisons across panels.

Extended Data Fig. 8. Normalized root mean square error values from linear and nonlinear models predicting various biodiversity metrics.

a, Community completeness. b, Alpha diversity. c, Dark diversity. d, Species pool size. e, Gamma diversity. f, Beta diversity. We used fivefold spatial cross-validation (see Extended Data Fig. 2) with bootstrapping to estimate the variation. The box plots (centre line = median; box limits = upper and lower quartiles; whiskers = the range, excluding outlying points that exceed the quartiles by more than 1.5× the interquartile range) show that while the nonlinear models had lower errors for the training set, the test data were predicted with lower error by the linear models.