Data reporting
The data underpinning this study is a compilation of existing datasets and therefore, no statistical methods were used to predetermine sample size, the experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment. The measurements were taken from distinct samples, repeated measurements from the same sites were averaged in the main analysis.
Inclusion & ethics
Data were primarily collected from individual archives of contributing co-authors. The data collection initiative was openly announced via the mailing list of the 10th International Seminar on Apterygota and via social media (Twitter, Researchgate). In addition, colleagues from less explored regions (Africa, South America) were contacted via personal networks of the initial authors group and literature search. All direct data providers who collected and standardised the data were invited as co-authors with defined minimum role (data provision and cleaning, manuscript editing and approval). For unpublished data, people who were directly involved in sorting and identification of springtails, including all local researchers, were invited as co-authors. Principal investigators were normally not included as co-authors, unless they contributed to conceptualisation and writing of the manuscript. All co-authors were informed and invited to contribute throughout the research process—from the study design and analysis to writing and editing. The study provided an inclusive platform for researchers around the globe to network, share and test their research ideas.
Data acquisition
Both published and unpublished data were collected, using raw data whenever possible entered into a common template. In addition, data available from Edaphobase47 was included. The following minimum set of variables was collected: collectors, collection method (including sampling area and depth), extraction method, identification precision and resources, collection date, latitude and longitude, vegetation type (generalized as grassland, scrub, woodland, agriculture and other for the analysis), and abundances of springtail taxa found in each soil sample (or sampling site). Underrepresented geographical areas (Africa, South America, Australia and Southeast Asia) were specifically targeted by a literature search in the Web of Science database using the keywords ‘springtail’ or ‘Collembola’, ‘density’ or ‘abundance’ or ‘diversity’, and the region of interest; data were acquired from all found papers if the minimum information listed above was provided. All collected datasets were cleaned using OpenRefine v3.3 (https://openrefine.org) to remove inconsistencies and typos. Geographical coordinates were checked by comparing the dataset descriptions with the geographical coordinates. In total, 363 datasets comprising 2783 sites were collected and collated into a single dataset (Supplementary Fig. 1).
Calculation of community parameters
Community parameters were calculated at the site level. Here, we defined a site as a locality that hosts a defined springtail community, is covered by a certain vegetation type, with a certain management, and is usually represented by a sampling area of up to a hundred metres in diameter, making species co-occurrence and interactions plausible. To calculate density, numerical abundance in all samples was averaged and recalculated per square metre using the sampling area. Springtail communities were assessed predominantly during active vegetation periods (i.e., spring, summer and autumn in temperate and boreal biomes, and summer in polar biomes). Our estimations of community parameters therefore refer to the most favourable conditions (peak yearly densities). This seasonal sampling bias is likely to have little effect on our conclusions, since most springtails survive during cold periods38,48. Finally, we used mean annual soil temperatures49 to estimate the seasonal mean community metabolism (described below) and tested for the seasonal bias in additional analysis (see Linear mixed-effects models).
All data analyses were conducted in R v. 4.0.250 with RStudio interface v. 1.4.1103 (RStudio, PBC). Data was transformed and visualised with tidyverse packages51,52, unless otherwise mentioned. Background for the global maps was acquired via the maps package53,54. To calculate local species richness, we used data identified to species or morphospecies level (validated by the expert team). Since the sampling effort varied among studies, we extrapolated species richness using rarefaction curves based on individual samples with the Chao estimator51,52 in the vegan package53. For some sites, sample-level data were not available in the original publications, but site-level averages were provided, and an extensive sampling effort was made. In such cases, we predicted extrapolated species richness based on the completeness (ratio of observed to extrapolated richness) recorded at sites where sample-level data were available (only sites with 5 or more samples were used for the prediction). We built a binomial model to predict completeness in sites where no sample-level data were available using latitude and the number of samples taken at a site as predictors: glm(Completeness~N_samples*Latitude). We found a positive effect of the number of samples (Chisq = 1.97, p = 0.0492) and latitude (Chisq = 2.07, p = 0.0391) on the completeness (Supplementary Figs. 17–19). We further used this model to predict extrapolated species richness on the sites with pooled data (435 sites in Europe, 15 in Australia, 6 in South America, 4 in Asia, and 3 in Africa).
To calculate biomass, we first cross-checked all taxonomic names with the collembola.org checklist55 using fuzzy matching algorithms (fuzzyjoin R package56) to align taxonomic names and correct typos. Then we merged taxonomic names with a dataset on body lengths compiled from the BETSI database57, a personal database of Matty P. Berg, and additional expert contributions. We used average body lengths for the genus level (body size data on 432 genera) since data at the species level were not available for many morphospecies (especially in tropical regions), and species within most springtail genera had similar body size ranges. Data with no genus-level identifications were excluded from the analysis. Dry and fresh body masses were calculated from body length using a set of group-specific length-mass regressions (Supplementary Table 1)58,59 and the results of different regressions applied to the same morphogroup were averaged. Dry mass was recalculated to fresh mass using corresponding group-specific coefficients58. We used fresh mass to calculate individual metabolic rates60 and account for the mean annual topsoil (0–5 cm) temperature at a given site61. Group-specific metabolic coefficients for insects (including springtails) were used for the calculation: normalization factor (i0) ln(21.972) [J h−1], allometric exponent (a) 0.759, and activation energy (E) 0.657 [eV]60. Community-weighted (specimen-based) mean individual dry masses and metabolic rates were calculated for each sample and then averaged by site after excluding 10% of maximum and 10% of minimum values to reduce impact of outliers. To calculate site-level biomass and community metabolism, we summed masses or metabolic rates of individuals, averaged them across samples, and recalculated them per unit area (m2).
Parameter uncertainties
Our biomass and community metabolism approximations contain several assumptions. To account for the uncertainty in the length-mass and mass-metabolism regression coefficients, in addition to the average coefficients, we also used maximum (average + standard error) and minimum coefficients (average—standard error; Supplementary Table 1) in all equations to calculate maximum and minimum estimations of biomass and community metabolism reported in the main text. Further, we ignored latitudinal variation in body sizes within taxonomic groups62. Nevertheless, latitudinal differences in springtail density (30-fold), environmental temperature (from −16.0 to +27.6 °C in the air and from −10.2 to +30.4 °C in the soil), and genus-level community compositions (there are only few common genera among polar regions and the tropics)55 are higher than the uncertainties introduced by indirect parameter estimations, which allowed us to detect global trends. Although most springtails are concentrated in the litter and uppermost soil layers20, their vertical distribution depends on the particular ecosystem63. Since sampling methods are usually ecosystem-specific (i.e. sampling is done deeper in soils with developed organic layers), we treated the methods used by the original data collectors as representative of a given ecosystem. Under this assumption, we might have underestimated the number of springtails in soils with deep organic horizons, so our global estimates are conservative and we would expect true global density and biomass to be slightly higher. To minimize these effects, we excluded sites where the estimations were likely to be unreliable (see data selection below).
Data selection
Only data collection methods allowing for area-based recalculation (e.g. Tullgren or Berlese funnels) were used for analysis. Data from artificial habitats, coastal ecosystems, caves, canopies, snow surfaces, and strong experimental manipulations beyond the bounds of naturally occurring conditions were excluded (Supplementary Fig. 1). To ensure data quality, we performed a two-step quality check: technical selection and expert evaluation. Collected data varied according to collection protocols, such as sampling depth and the microhabitats (layers) considered. To technically exclude unreliable density estimations, we explored data with a number of diagnostic graphs (Supplementary Table 2; Supplementary Figs. 12–20) and filtered it, excluding the following: (1) All woodlands where only soil or only litter was considered; (2) All scrub ecosystems where only ground cover (litter or mosses) was considered; (3) Agricultural sites in temperate zones where only soil with sampling depth <10 cm was considered. Additionally, 10% of the lowest values were individually checked and excluded if density was unrealistically low for the given ecosystem (outliers with density over three times lower than 1% percentile within each ecosystem type). In total, 237 sites were excluded from density, and 394 sites from biomass, and community metabolism analyses based on these criteria (Supplementary Figs. 15 and 16). For the local species richness estimates, we removed all extrapolations based on sites with fewer than three samples and no (morpho)species identifications (647 sites; Supplementary Figs. 1 and 20).
Data expert evaluation
We performed manual expert evaluation of every contributed dataset. Evaluation was done by an expert board of springtail specialists, each with extensive research experience in a certain geographic area: Anatoly Babenko—high latitude regions in both north and south hemispheres; Bruno Bellini—Central and South America; Jean-François Ponge—Central and Western Europe; Louis Deharveng—Africa and Asia; Lubomir Kovac—Southern Europe; Mikhail Potapov and Natalia Kuznetsova—Eastern and Northern Europe. Each dataset was scored by the experts separately for density and species richness estimation as either trustworthy, acceptable, or unreliable. Density estimation quality was assessed using information about the sampling and extraction method and the density estimation itself. Species richness estimation quality was assessed using information about the identification key, experience of the person who identified the material, species (taxa) list, and the species richness estimation itself. Based on the expert opinions, unreliable estimates of density (together with biomass and community metabolism) and species richness were excluded (Supplementary Fig. 1). The resulting final dataset included 2470 sites and 43,601 samples64 with a median of six samples collected at each site. The dataset comprised 2210 sites with density estimation (69–2,181,600 individuals m−2), 2,053 sites with mean fresh body mass (1.8–3110 µg), mean metabolic rate (0.028–2.4 mJ h−1), dry biomass (0.5–93,000 mg m−2), fresh biomass (1.6–277,000 mg m−2) and community metabolism estimations (0.03–1000 J h−1), and 1735 sites with local species richness estimation (1–136.7 species; Supplementary Figs. 1 and 2). The dataset covered all major biomes (Supplementary Fig. 3), years 1970–2019, and all months: 8% of the samples were taken between December and February, 14% between March and May, 55% between June and August, and 23% between September and November (see Data availability).
Data transformation
All parameters except for extrapolated local species richness were highly skewed (e.g., density had a global median of 21,016 individuals m−2 and a mean of 60,454 individuals m−2) and we applied log10-transformation prior to analysis. This greatly improved the fit of all statistical analyses.
Latitudinal and ecosystem trends
To explore changes in springtail communities with latitude, we sliced the global latitudinal gradient into 5-degree bins and calculated average parameters across sites in each bin after trimming to ensure the same statistical weight for each latitudinal bin while plotting the gradient. The latitudinal gradient was plotted with ggplot265, and quadratic smoothers were used to illustrate trends. Mean parameters of springtail communities were compared across ecosystem types using a linear model and multiple comparisons with the Tukey HSD test using HSD.test in the agricolae package66. Habitats were classified according to the vegetation types. Climates were classified as polar (beyond the polar circles, i.e., more than 66.5 and less than −66.5 degrees), temperate (from the polar circles to the tropics of Capricorn/Cancer, i.e. to 23.5 and −23.5 degrees) and tropical (in between 23.5 and −23.5 degrees). Habitats and climates were combined to produce ecosystem types. For the analysis, only well-represented ecosystem types were retained: polar scrub (n = 253), polar grassland (n = 39), polar woodland (n = 28), temperate woodland (n = 907), temperate scrub (n = 104), temperate grassland (n = 445), temperate agriculture (n = 374), tropical agriculture (n = 68) and tropical forest (n = 141; Supplementary Fig. 3).
Selection of environmental predictors
To assess the predictors of global distributions of springtail community metrics, we pre-selected variables with a known ecological effect on springtail communities (based on expert opinions) and constructed a hypothetical relationship diagram (Supplementary Fig. 9a). Environmental data were very heterogeneous across the springtail studies, so we used globally available climatic and other environmental layers. Overall, we included global layers bearing the following information: climate (mean annual air temperature, air temperature seasonality, air temperature annual range, mean annual precipitation, precipitation seasonality, precipitation of the driest quarter67, inversed aridity index68), topography (elevation, roughness69), vegetation and land cover (aboveground biomass70, tree cover71, Net Primary Production, Normalized Difference Vegetation Index [NDVI]72), topsoil physicochemical properties (0–15 cm depth C to N ratio, pH, clay, sand, coarse fragments, organic carbon, bulk density73) and human population density74. Some of environmental layers could not be included due to the lack of appropriate data. For example, soil phosphorus and nitrogen concentrations had to be omitted. While the global distribution of soil nitrogen concentration is available73, it is a modelled product, which strongly correlates with soil carbon concentration, and thus cannot be used as an independent predictor.
Geospatial global projections
To create global spatial predictions of springtail density, species richness, biomass, and community metabolism, we followed the approach previously used for nematodes10,75 that is based on spatial associations of community parameters with global environmental information. The analysis for geospatial modelling was done in Python version 3.6.5 (Python software foundation). A Random Forest algorithm was applied to identify the spatial associations and extrapolate local observations to the global scale at the 30 arcsec resolution (approximately 1 km2 pixels)18,75. After retrieving the environmental variable values for each location, we trained 18 model versions, each with different hyperparameter settings, i.e., variables per split (range: 2–7); minimum leaf population (range: 3–5). To minimize the potential bias of a single model, we used an ensemble of the top 10 best-performing models, selected based on the coefficient of determination (R2), to create global predictions of each of the community parameters.
Model performance was assessed by 10-fold cross validation, with folds assigned randomly. The R2 values for each of the five response variables were in the range of 0.30–0.57 (density: 0.567, dry biomass: 0.463, community metabolism: 0.359, extrapolated species richness: 0.302). For some of the modelled variables we observed positive spatial autocorrelation: at ranges below 150 km for density, below 100 km for community metabolism and below 150 km for extrapolated species richness (Supplementary Note). Yet, the Moran’s I values were very close to zero (the highest value was 0.07), indicating that the effect of spatial autocorrelation was very weak. These results were obtained by performing Moran’s I tests using the spatialRF package in R76. To investigate the effect of spatial autocorrelation on model performance we performed a buffered leave-one-out cross-validation tests (described in detail as an alternative performance statistic for models with potential spatial autocorrelation77,78). As expected, the predictive power declined with increasing buffer sizes. At the scales at which we observe positive autocorrelation, i.e., where we have significant Moran’s I values, coefficient of determination remained positive.
To reduce potential artifacts produced by extrapolation, geographical regions with climatic conditions poorly represented by our sites and without NPP data were excluded from the extrapolation (e.g., Sahara, Arabian desert, Himalayas). We evaluated our extrapolation quality based on spatial approximations of interpolation versus extrapolation75. In this approach, we first determined the range of environmental conditions represented by the observations. Next, we classified all pixels to fall within or outside the training space, in univariate and multivariate space. For the latter, we first transformed the data into principal component space, and selected the first 11 PC axes, collectively explaining 90% of the variation. Finally, we classified pixels to fall within or outside the convex hulls drawn around each possible bivariate combination of these 11 PC axes; pixels that fell outside the convex hulls in >90% of cases were masked on the main maps; for the map with density-species richness visualisation, two corresponding masks were applied (Fig. 2).
To estimate spatial variability of our predictions while accounting for the spatial sampling bias in our data (Fig. 1a) we performed a spatially stratified bootstrapping procedure. We used the relative area of each IPBES79 region (i.e., Europe and Central Asia, Asia and the Pacific, Africa, and the Americas) to resample the original dataset, creating 100 bootstrap resamples. Each of these resamples was used to create a global map, which was then reduced to create mean, standard deviation, 95% confidence interval, and coefficient of variation maps (Supplementary Figs. 4–7).
Global biomass, abundance, and community metabolism of springtails were estimated by summing predicted values for each 30 arcsec pixel10. Global community metabolism was recalculated from joule to mass carbon by assuming 1 kg fresh mass = 7 × 106 J80, an average water proportion in springtails of 70%58, and an average carbon concentration of 45% (calculated from 225 measurements across temperate forest ecosystems)81. We repeated the procedure of global extrapolation and prediction for biomass and community metabolism using minimum and maximum estimates of these parameters from regression coefficient uncertainties (see Parameter uncertainties).
Path analysis
To reveal the predictors of springtail communities at the global scale, we performed a path analysis. After filtering the selected environmental variables (see above) according to their global availability and collinearity, 13 variables were used (Supplementary Fig. 9b): mean annual air temperature, mean annual precipitation (CHELSA database67), aridity (CGIAR database68), soil pH, sand and clay contents combined (sand and clay contents were co-linear in our dataset), soil organic carbon content (SoilGrids database73), NDVI (MODIS database72), human population density (GPWv4 database74), latitude, elevation69, and vegetation cover reported by the data providers following the habitat classification of European Environment Agency (woodland, scrub, agriculture, and grasslands; the latter were coded as the combination of woodland, scrub, and agriculture absent). Before running the analysis, we performed the Rosner’s generalized extreme Studentized deviate test in the EnvStats package82 to exclude extreme outliers and we z-standardized all variables (Supplementary R Code).
Separate structural equation models were run to predict density, dry biomass, community metabolism, and local species richness in the lavaan package83. To account for the spatial clustering of our data in Europe, instead of running a model for the entire dataset, we divided the data by the IPBES79 geographical regions and selected a random subset of sites for Eurasia, such that only twice the number of sites were included in the model as the second-most represented region. We ran the path analysis 99 times for each community parameter with different Eurasian subsets (density had n = 723 per iteration, local species richness had n = 352, dry biomass had n = 568, and community metabolism had n = 533). We decided to keep the share of the Eurasian dataset larger than other regions to increase the number of sites per iteration and validity of the models. The Eurasian dataset also had the best data quality among all regions and a substantial reduction in datasets from Eurasia would result in a low weight for high-quality data. We additionally ran a set of models in which the Eurasian dataset was represented by the same number of sites as the second-most represented region, which yielded similar effect directions for all factors, but slightly higher variations and fewer consistently significant effects. In the paper, only the first version of analysis is presented. To illustrate the results, we averaged effect sizes for the paths across all iterations and presented the distribution of these effect sizes using mirrored Kernel density estimation (violin) plots. We marked and discussed effects that were significant at p < 0.05 in more than a given number of iterations (arbitrary thresholds were set to 25%, 50%, 75% and 95% of iterations; Fig. 3).
Linear mixed-effects modelling
To test if our results are biased by seasonal effects, sampling methodology, and/or species richness extrapolation, we selected a subset of sampling events with known sampling year and month (2997 sampling events representing 1703 sites) and ran linear mixed-effects models for springtail density, species richness (both raw and extrapolated), biomass, and community metabolism. The models were run using the lme4 package84. Data were transformed as described above and analysed using Gaussian distributions except for raw species richness, which was analysed using generalised models with Poisson distribution. Sampling site was included as random effect to account for the dependence of the sampling events coming from the same sites. We included mean monthly air temperatures (offset from the annual mean) and the sum of monthly precipitation at the sampling month as additional climatic predictors. We also included total collection area and the presence of litter (or any other soil cover such as mosses) and soil in the sample to account for methodological biases. All models were run using the full dataset (n = 2884 sampling events for density, n = 2540 for raw species richness; n = 1708 for extrapolated species richness; n = 2462 for dry biomass; n = 2289 for community metabolism). To test if the effect of temperature on species richness is non-linear, we additionally ran the same model including quadratic function poly(MAT, 2).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Source: Ecology - nature.com
