The first two freshwater data sets originate from a state-wide lake survey program in Iowa, United States (Figure 2), with samplings taking place three times a year in early, mid and late summer (Arbuckle and Downing, 2001). Phytoplankton data were collected using integrated water column samples from the epilimnion (2 m maximum depth) and identified to genus level using light microscopy. Zooplankton samples were collected by vertically towing a Wisconsin net (63 μm mesh) from the thermocline to water surface and identified to the most practical taxonomic level. Biomass values are given in μgL⁻¹ (zooplankton) and mgL⁻¹ (phytoplankton). Filstrup et al. (2014a) and Arbuckle and Downing (2001) provide a comprehensive overview on further methodological details and the study area. The study region covers the entire state of Iowa (146,000 km²). The distance between lake pairs varies between 1 and 484 km. Since missing data points in the biomass time series of the different lakes would have introduced bias when calculating temporal stability of aggregate biomass values, we created a subset of lakes with complete sampling records only. Our subset contains data from 49 lakes (Figure 2) covering all three annual sampling rounds from 2001 to 2007 resulting in 21 data points per time series.

The third data set was collected as part of a monitoring program consisting of annual autumn samples on North Sea macrozoobenthos communities off the Northern German coast (Figure 2). Again, we reduced the data set to include only stations with complete sampling records, resulting in a subset of 72 stations over a sampling period of 7 years (2005–2011). Since the original data consisted of abundance records only, we converted it to biomass values (wet weight gm⁻²) using individual biomass estimates from the same species if possible or otherwise the same genus obtained from other North Sea macrozoobenthos sampling campaigns [source: data information systems MARLIN (Federal Maritime and Hydrographic Agency, https://lindevmarlin61.bsh.de/MARLINDMZ/publicSites/MainAppPublic.jsf), CRITTERBASE (Alfred Wegener Institute, https://critterbase.awi.de/)]. Across the 72 stations and 7 years, 195 species were identified. In this data set, sampling sites are arranged across a regular grid with a distance of roughly 10 km between two neighboring stations and a maximum distance of 115 km.

Patterns of local and regional biodiversity and community structure are distinct in each of the monitored species communities. Calculations of Pielou’s evenness index [ranging from complete dominance (0) to complete evenness (1)] revealed on average rather uneven community structures for both phytoplankton (mean = 0.348, std. = 0.053) and macrozoobenthos (mean = 0.399, std. = 0.043) at the local scale. In the phytoplankton data set, the uneven community structure is mainly driven by the strong dominance of the Cyanobacteria genus Microcystis. Similarly, in the North Sea macrozoobenthos data, two dominant species Echinocardium cordatum and Ensis leei in combination with a large number of species with very little and infrequent biomass contributions lead to the low evenness values (Supplementary Figure S1). Dominance patterns in local zooplankton communities were less pronounced (mean = 0.565, std. = 0.024). The average biomass-based Bray–Curtis dissimilarity between pairs of local communities also varied among the phytoplankton (mean = 0.573, std. = 0.064), zooplankton (mean = 0.454, std. = 0.035) and macrozoobenthos (mean = 0.556, std. = 0.041) communities.

We used two approaches to describe the dependence of temporal stability of aggregate community biomass measurements on the spatial extent of the studied area: (1) The invariability-area relationship (IAR, Wang et al., 2017), and (2) a stability-synchrony framework introduced by Wang et al. (2019).

The invariability-area relationship describes the relationship of the invariability I (i.e., stability) of a measured ecological variable and the spatial extent of the study area (A). It is formulated as the inverse of the squared coefficient of variation:

where 𝐴 is the area of aggregation and 𝐶𝑉 is the coefficient of variation of the ecological variable of interest.

In order to investigate how spatial context, i.e., spatial autocorrelation, might affect the slope of the IAR, we apply two different approaches for area increase: (i) the addition of the next closest sampling site for a nested increase of spatial scale and (ii) random addition of sampling sites irrespective of spatial context. We then calculated IARs for all sampling sites. In our study, spatial extent is represented by the number of sampling sites over which the community biomass is aggregated so the increase in spatial extent across the IARs of different stations varies slightly depending on which sampling site was used as starting point.

For the second analysis, we follow the framework introduced by Wang et al. (2019) that allows the distinction between the relative contributions of species vs. spatial insurance effects that increase stability with the level of spatial data aggregation. This partitioning approach results in four measures of variability, represented by the squared coefficient of variation (sensu Wang and Loreau, 2014) at two levels of ecological [species (S), community (C)] as well as spatial [local (L), regional (R) scale (Figure 3)]. These are (1) the variability of local populations (CV_S,L) which is calculated as the weighted average of local population biomass variability across species and sampling locations and stands for the lowest level of spatial as well as ecological aggregation, (2) the variability of local communities (CV_C,L) calculated as the weighted average of community biomass variability across sampling locations, which represents the next higher or intermediate level of ecological aggregation, (3) the meta-population variability (CV_S,R) calculated as the weighted average of the pooled meta-population biomass variability, which serves as a measure for the intermediate level of spatially aggregated variability, and (4) the variability at the meta-community level (CV_C,R) calculated as the weighted overall meta-community variability, i.e., pooled biomass values across all sites and species. For all measures, variability is weighted by the respective overall biomass of a species or community. The up-scaling of variability is possible either by first aggregating biomass values across space and then across levels of ecological organization, i.e., CV_S,L → CV_S,R → CV_C,R, or vice versa, i.e., CV_S,L → CV_C,L → CV_C,R. The variability at the next higher level of data aggregation can be calculated as the product of lower level variability and the synchrony between its components, i.e., at the species (ρ_S) or community level (ρ_C): CV_S,L * ρ_{S → C, L} = CV_{C, L}; CV_S,L * ρ_{S, L → R} = CV_S,R; CV_C,L * ρ_{C, L → R} = CV_C,R; and CV_S,R * ρ_{S → C, R} = CV_C,R.

We then generated null model distributions of uncorrelated species population dynamics using cyclic shift permutation (Hallett et al., 2016; Lamy et al., 2019). This approach creates a null model community matrix by assigning random starting points to each individual time series and therefore removing correlations between population dynamics across species and sites while preserving most of the temporal autocorrelation of each individual population. The algorithm has been criticized for not completely preserving autocorrelation structure of species-specific population dynamics (Kalyuzhny, 2020), but for the purpose of comparing between independent population dynamics and dynamics that are driven by ecological mechanisms affecting synchrony patterns across space, these shortcomings are of minor importance.

As sites and species will likely contribute differently to stabilizing or destabilizing mechanisms across space, we ran analyses on all data subsets each time excluding one species (or sampling site), and calculated log response ratios of meta-population and spatial community synchrony to quantify the contribution of each species (or sampling site) to synchrony [i.e., LRR = ln(rho’/rho), sensu Lamy et al., 2019]. We then related the sign and degree of synchrony to site and species-specific characteristics. For sites we used the total, mean and variability of temporal community biomass, effective number of species (ENS), as well as a sites’ compositional uniqueness, represented by its average Bray Curtis dissimilarity from all other sampling locations. As indicators of species characteristics, we calculated a species’ total biomass across sites and sampling occasions, as well as temporal biomass variability, and the heterogeneity of a species’ distribution across the study area, i.e., spatial variability of a species’ biomass.

All statistical analyses were performed using R statistical computing (R Core Team, 2021). The code for the analyses is accessible under https://github.com/DorotheeHodapp/Code-Spatial-Scaling-of-Ecology.

Overall, invariability was lower in the phytoplankton dataset than for the other two organism groups, which was consistent across spatial scales, i.e., local population biomass as well as total metacommunity biomass showed higher temporal fluctuations in the phytoplankton samples (Figure 4). Differences between temporal biomass stability was less pronounced between the species groups at the larger spatial scale indicating a lower degree of stabilization with increasing spatial scale for the zooplankton and macrozoobenthos meta-communities (Figure 4A). For all three species groups, the slope of the IAR was less than one indicating that the observed species dynamics were more synchronous across space than would be expected in case of completely independent biomass fluctuations. The 1:1 line represents the limiting case of completely uncorrelated local species communities whereas a horizontal line indicates complete synchrony across local patches and consequently no increase in invariability with increasing spatial scale. At the local scale, i.e., at a single site, invariability can vary depending on local biomass fluctuations. In our data, both zooplankton and macrozoobenthos communities showed considerably higher invariability than the highly variable phytoplankton communities. This can be caused by either higher synchrony of local species populations or higher overall variability in biomass measurements. The comparison between the two approaches of area increase by aggregating either the next closest sites (Figure 4A) or adding them randomly (Figure 4B) revealed differences in the synchrony decay with increasing level of site aggregation indicating spatially structuring elements in all three data sets.

The second analysis confirmed the IAR results in that temporal variability in local species populations was highest for phytoplankton, followed by zooplankton and macrozoobenthos, which was consistent across all levels of ecological and spatial scale (Figure 3). The reduction in variability when aggregating at the meta-community level was similar for all species groups. Biomass weighted CV² decreased from local species population variability (CV_S,L) to overall meta-community variability (CV_C,R) by a factor of just over three.

Local and regional scale species synchrony yielded synchrony coefficients that were on average twice as large as the spatial synchrony coefficients leading to stronger reductions in variability between spatial levels of aggregation (local ➔ regional) in comparison to the two levels of ecological aggregation (species ➔ community). While the phytoplankton communities showed the highest degrees of species synchrony, they exhibited lower levels in spatial synchrony than the other two species groups. Synchrony measures of zooplankton and macrozoobenthos were similar, although synchrony among zooplankton communities was only slightly but consistently smaller.

The comparison of the synchrony coefficients between the biomass measurements in the three data sets and completely uncorrelated simulated biomass time series yielded higher synchrony levels within and between the natural communities than for the randomly fluctuating population biomasses (Figure 5).

Overall, absolute values of single species and single site contributions to the stabilization of meta-community biomass follow a clear pattern (Figure 6). A species can have disproportionately high effects on the spatial scaling of stability if it has high total biomass values across space and time or high overall temporal or spatial variability in biomass compared to the other species (Figures 6A–C; Supplementary Table S7). Similarly for sites, high values of temporal mean and variability in biomass can result in greater effects on stabilization or destabilization of temporal biomass values across spatial scales (Figures 6D,E). Importantly, the effects can be positive or negative (Supplementary Figure S6). Sign and magnitude varied in our data depending on species and community characteristics. The effects of community characteristics, i.e., compositional uniqueness, temporal mean and variability of the effective number of species (ENS) on stability patterns across scales differed between the three species groups (Figures 6F–H), most of the relationships were non-significant (Supplementary Table S7). Phytoplankton and macrozoobenthos communities showed higher variability with increasing compositional uniqueness of the community and the mean and temporal variability of species diversity (ENS). For the zooplankton data, there was neither a clear trend in slope nor variation (Supplementary Figures S6F–H).

A comparison of IAR shapes across the three organism groups showed on average lower stability levels for phytoplankton communities, i.e., higher temporal variability in biomass than for the two other groups. This difference was consistent at local and regional scales. Potential reasons are the extreme dominance and destabilizing effect of the cyanobacteria Microcystis and the short doubling times of plankton species. In addition, annual bloom formation and succession would likely introduce some synchrony in biomass fluctuations across sampling sites, which could partly be captured by the Iowa lake sampling scheme with samples taken in early, mid and late summer. Zooplankton communities also follow seasonal patterns, but show less pronounced seasonal dynamics as a consequence of longer generation times in comparison to phytoplankton. Their seasonal dynamics are also lower in eutrophic lakes (Sommer et al., 1986) especially if the dominant phytoplankton species are of low edibility or produce toxins like many bloom-forming cyanobacteria species do (Paerl and Otten, 2013). Using the same dataset, Filstrup et al. (2014a) and Heathcote et al. (2016) showed less zooplankton biomass per unit phytoplankton biomass when communities were dominated by few taxa (largely Microcystis). In addition, the lower levels of synchrony in zooplankton and macrozoobenthos communities might be caused by their often patchy distributions in space (Gray, 2002; Armonies and Reise, 2003; Gutow et al., 2020). For the macrozoobenthos data, the yearly sampling scheme and the variable dominance patterns, with different species dominating standing stock biomass across years (Reiss and Kröncke, 2005; Dannheim and Rumohr, 2012; Hodapp et al., 2014), reduced synchrony between sampling locations. For instance, in our data set, the two dominating macrozoobenthos species Echinocardium cordatum and Ensis leei only had peak abundances in a small subset of years over the sampling period and the spatial location of maximum standing stock values also varied considerably between the years (Supplementary Figure S5).

The mean curves of all three species groups (Figure 4A) resemble IAR patterns caused by an exponential synchrony decay, where the correlation between sites abruptly decreases to zero exceeding a certain distance “characteristic correlation length” (Wang et al., 2017). In contrast, our second approach of randomly adding sampling sites irrespective of geographic proximity of single sampling locations yielded IAR curves following a linearly increasing pattern on log–log scale (Figure 4B). This difference concurs with hypotheses on effects of spatially structuring heterogeneity on inter-site synchrony decay across space (Liebhold et al., 2004; Walter et al., 2017). These effects are likely stronger than depicted in our analyses as the necessary averaging across all possible IARs for our set of lakes dilutes effects of individual IAR curves. On the other hand, mean and variability of pairwise correlations across sites are similar for any given distance in the data sets (see Supplementary Figure S3). Reasons for similarities independent of catchment or spatial proximity include the type of land use or location in a watershed. Indeed, there are broad dissimilarities among lakes and watersheds in intensity of the drivers of plankton dynamics (Downing et al., 2008; Filstrup et al., 2016). The macrozoobenthos communities on the other hand are usually associated with specific sediment types introducing local synchrony (Salzwedel et al., 1985; Fiorentino et al., 2017). The patch size of different sediment types can vary greatly thus leading to environmental heterogeneity at much smaller scales than represented by the sampling scheme. One environmental factor likely enhancing synchrony across the whole study area is winter temperature since many macrozoobenthos organisms show lower biomass following cold winters (Neumann et al., 2008). Other common drivers of synchrony across space, such as dispersal (Ranta et al., 2008), community similarity and species interactions, can enhance or even enable the synchronizing effects of the other factors (Vasseur and Fox, 2009). However, these might only play a minor role in our example data sets. Dispersal rates among local populations of macrozoobenthos larval stages are indeed high across space due to recurring exchange of water masses during tidal cycles. However, successful establishment after dispersal driven by differential mortality across species due to physical constraints and biotic interactions (Kröncke and Reiss, 2010; Dannheim and Rumohr, 2012) which combined can lead to low community similarity values.

All three data sets show phases of decreasing stability around intermediate levels of spatial scale. Especially, many of the individual zooplankton and macrozoobenthos IAR relationships show sudden declines in stability after adding specific (highly unstable) sites. It seems that for these two species groups, high variability in biomass more frequently coincides in sites with high biomass values on average and in total (see Supplementary Figure S4) than is the case for the phytoplankton data. This suggests that the correlation between variability and magnitude of biomass across sites (and also species), which can be substantial, likely influences the shape of the IAR curve.

It should be noted that the IARs are affected by a number of study specific factors. As the distance between sites varies, substituting area increase with the addition of a site introduces variability. In addition, as our samples are point estimates, i.e., only represent a fraction of the actual area covered by the respective monitoring scheme, the slope of the IAR might be overestimated as shown in Wang et al. (2017). This is supported by the weak relationship between pairwise temporal correlations of biomass at two sampling locations and their distance (Supplementary Figure S3). Despite these deviations from the model system applied in Wang et al. (2017), the general patterns are consistent across the three data sets and differ from the IARs generated by our second approach of randomly adding sampling sites.

For all three data sets, we found stronger stabilizing effects as a consequence of asynchronous biomass fluctuations across space than across species. As previously discussed, seasonal and annual environmental forcing and the sampling scheme likely enhance synchrony across species outweighing much of the temporal variability in biomass dynamics due to species-specific responses. Consequently, depending on the sampling scheme the results could vary as effects of spatial and temporal grain and extent of sampling greatly determine which ecological processes can be detected. In contrast to species-specific fluctuations, local environmental conditions seem to stabilize meta-community biomass by introducing variability in biomass dynamics across space. Drivers of local differences are likely caused by varying land use type (Arbuckle and Downing, 2001) and catchment size which influences nutrient delivery to the lakes (Downing et al., 2008). As a consequence, the lakes and their communities exhibit differing responses to limiting nutrient concentrations. Previous regional studies have also shown that regional characteristics, such as land use, can modify biomass response curves to limiting nutrients (Filstrup et al., 2014a; Filstrup et al., 2014b). In addition, bloom formation depends on several factors other than nutrient availability (Purz et al., 2021) and can therefore be restricted by lake or community specific aspects.

Other examples in the literature found stronger species insurance effects in a kelp forest meta-community, which they explain by niche differentiation among algae species leading to diverse population dynamics due to, e.g., shading tolerance, susceptibility to grazing and wave disturbance (Lamy et al., 2019). In addition, they suspected the high spatial synchrony to be caused by large scale variation in environmental forcing, e.g., sea surface temperature. Another study on spatial scaling of stability in benthic fish communities reported that spatial compensatory dynamics were prevalent in their data and were mainly related to changes in aggregate biomass variability of a certain region (Thorson et al., 2018).

A comparison between spatial and species synchrony in our data sets and simulated meta-communities with randomly fluctuating population biomasses revealed the same or stronger synchrony patterns in the natural communities than expected by chance (also shown in the lower levels of IAR compared to 1:1 line in Figure 4). Given that natural communities at the scale of our study are usually exposed to a set of common environmental drivers, it seems plausible that biomass fluctuations will not be entirely random.

Our analyses show that the magnitude and fluctuations of dominant species’ biomass influence which effects single species and sites have on stabilizing processes. Although an addition of sites and species is generally assumed to result in an increase in stability as suggested by the diversity-stability hypothesis, an addition of a highly fluctuating high biomass species can also negatively influence the stability of the aggregated biomass. This is in agreement with the mass-ratio hypothesis (Grime, 1998) stating that community characteristics are usually determined by the traits of the most dominant species. It is also partly inherent in the way stability is measured in both frameworks. In communities with differences in individual biomasses of several magnitudes and highly fluctuating abundances across time and space, species with high total biomass are likely more influential in terms of overall fluctuation of biomass values, i.e., stability. Similar findings have been reported for macroalgal species communities (Lamy et al., 2019) and subtropical forests (Yu et al., 2020), where the influence of the dominant species also outweighed effects of biodiversity on stability. Further, combined effects of single species contributions to stability patterns are not additive, i.e., they cannot be predicted from the single species effects on stability (White et al., 2020). Our research illustrates how dominance is not the sole driver of stability effects. Whether biomass fluctuations of a dominant species result in stabilization or destabilization of community biomass very much depends on the biomass distribution across species and their fluctuations. In our phytoplankton data set, one dominant species is the main driver of community stability across scales with extremely negative effects on community stability. In contrast, the two dominant species in the macrozoobenthos community both stabilized community biomass fluctuations through partly asynchronous dynamics. Consequently, species loss can likewise have stabilizing as well as destabilizing effects on biomass stability and its scaling relationship with space.

In summary, our analysis and comparison of the relationship between ecological stability and increasing spatial extent for three different organism types highlights a few aspects of the spatial scaling of stability in natural species communities. (1) Natural communities are usually more synchronous than expected from completely random dynamics. (2) Despite overall similarities in shape and magnitude of the spatial scaling of stability across data sets of different organisms, our analyses showed that species communities with highly dynamic dominance patterns across space and time foster strong impacts of individual species or locations on stability patterns across space. (3) In addition to single species dynamics, community structure and spatial aggregation of species influence the shape and amount of stability increases with spatial scale. Therefore, our small subset of ecosystem and organism types already highlights a number of variations and potentially associated drivers of scaling of stability across space.