Community stability is related to animal diversity change

Abstract

Understanding the drivers of community stability in times of increasing anthropogenic pressure is an urgent issue. Biodiversity is known to promote community stability, but studies of the biodiversity–stability relationship rarely consider the full complexity of biodiversity change. Furthermore, finding generalities that hold across taxonomic groups and spatial and temporal scales remains challenging because most investigations have narrow taxonomic, spatial, and temporal scopes. We used organismal data collected through the National Ecological Observatory Network (NEON) at sites across the contiguous United States to evaluate linkages between community stability and biodiversity change for four taxonomic groups: small mammals, ground beetles, fish, and freshwater macroinvertebrates. We defined community stability as constancy of aggregate species' abundance. We quantified change in biodiversity as (1) dissimilarity in community taxonomic and functional composition and species replacement and richness change components of that dissimilarity and (2) change in species' abundance distributions as captured by change in species rank, richness, and evenness. We found that community stability increased with species replacement and with contribution of species replacement to overall dissimilarity for all taxonomic groups, but declined with increasing change in species richness and evenness. This is consistent with the notion that temporal fluctuations in species abundance can help stabilize community properties. We also found that community stability was highest when change in community functional composition was either lower or higher than expected given reshuffling of each community's taxonomic composition. This suggests that long‐term community stability can result from fluctuations of functionally similar species in assemblages with high taxonomic reshuffling. On the contrary, the functional uniqueness of fluctuating species compensates for lower taxonomic reshuffling to drive stabilization of community properties. Our study provides an initial assessment of the relationship between community stability and biodiversity change and illustrates the utility of fine temporal resolution data collected across ecosystems and biomes to understand the general mechanisms underlying biodiversity–stability relationships.

INTRODUCTION

In times of increasing anthropogenic pressure on many ecosystems (Ma et al., 2017) and the resulting global biodiversity loss on par with Earth's past mass extinctions (Barnosky et al., 2011), understanding the drivers of ecological stability is a pressing issue (Hautier et al., 2015). Biodiversity has long been credited with promoting stability (Pennekamp et al., 2018; but see May, 1972 and subsequent work by, e.g., Allesina & Tang, 2012 and Jacquet et al., 2016), and its loss has been shown to reduce temporal stability of communities (Allan et al., 2011; Schnabel et al., 2019) and ecosystem functions (Cardinale et al., 2012). Biodiversity change, however, is a complex phenomenon that goes beyond net gains or losses in species richness and abundance (i.e., alpha diversity). For example, alpha diversity often remains constant even when community composition undergoes dramatic shifts (Dornelas et al., 2014; Jarzyna & Jetz, 2017) such as those resulting from replacement of rare or specialist species with common or generalist ones (Buckland et al., 2017). As a result, the conclusions reached about the maintenance of ecological stability using solely alpha-diversity metrics might contrast with those that account for species identities (Magurran & Henderson, 2010). Investigations of the biodiversity–stability relationship, however, rarely consider this full complexity of biodiversity change. While ecological stability is a multidimensional concept that describes different aspects of system dynamics and response to perturbations (Donohue et al., 2016; Hillebrand et al., 2018), here we focus on community (i.e., compositional) stability defined as temporal constancy of aggregate species' abundance (Hillebrand et al., 2018; Micheli et al., 1999).

The ubiquity of the biodiversity crisis calls for finding generalities that hold across taxonomic groups and spatial and temporal scales. Challenges in accomplishing this persist because most experimental and empirical work on the relationship between community stability and biodiversity has poor replication across ecosystems or geographies. Grassland and aquatic systems remain disproportionately represented (Donohue et al., 2016), constraining the taxonomic scope largely to plants (e.g., García-Palacios et al., 2018; Sasaki & Lauenroth, 2011; Usinowicz et al., 2017), phytoplankton (Rocha et al., 2011; Vallina et al., 2017), and zooplankton (Howeth & Leibold, 2010; Pennekamp et al., 2018), with comparatively fewer studies focused on invertebrate and vertebrate fauna (e.g., fish—Shimadzu et al., 2013; small mammals—Kalies & Covington, 2012). Most studies additionally rely on experimental designs with limited or no spatial replication (Hector et al., 2010; Tilman et al., 2006; but see García-Palacios et al., 2018, Oehri et al., 2017, Wilcox et al., 2017, for examples for plants), or use controlled mesocosm environments that may not resemble the environmental complexity inherent to natural or seminatural ecosystems (Cottingham et al., 2001; Fried-Petersen et al., 2020). To uncover generalities in the relationships between stability and different dimensions of biodiversity change, studies should encompass broader spatial extents and longer temporal scales, and a diverse range of ecosystems and taxonomic groups, and be designed with standardized protocols that allow comparative investigations (Donohue et al., 2016; Fried-Petersen et al., 2020). Ultimately, these studies can improve our ability to understand the consequences of biodiversity change for community stability at larger spatial scales that are more relevant for ecosystem management.

The National Ecological Observatory Network (NEON) collects organismal data for several sentinel taxa across a broad (i.e., continental, spanning >1000 km) geographic sampling extent that spans varied ecological and climatic domains (Goodman et al., 2015; Kao et al., 2012) and a planned long temporal duration (i.e., >30 years), enabling a rich exploration of the relationship between multidimensional biodiversity change and community stability across macrosystems. NEON is a network of observational sites comprising 47 terrestrial (20 core and 27 gradient) and 34 freshwater aquatic (24 wadeable streams, seven lakes, and three rivers) locations distributed across 20 ecoclimatic domains representative of distinct landforms, vegetation types, climatic zones, and ecosystem dynamics across the contiguous United States, Alaska, Hawaii, and Puerto Rico (Keller et al., 2008; Schimel et al., 2007). Organismal information collected by NEON encompasses a wide range of taxa from both terrestrial and freshwater ecosystems, while standardized data collection for each taxonomic group enables spatial and temporal replication. Here, we use NEON organismal data to explore links between community stability, defined as temporal constancy of aggregate species' abundance (Micheli et al., 1999), and different dimensions of biodiversity change. While NEON can still be considered at the beginning of its life span by the time of this publication, sufficiently long (i.e., 4–7 years) time series are already available for a number of taxonomic groups, allowing investigation of biodiversity change through time, and establishing and testing the protocol to be re-employed when more data become available through continued sampling.

The full complexity of biodiversity change can be captured via a suite of complementary metrics. Temporal dissimilarity (also termed temporal beta diversity; Baselga, 2010; Carvalho et al., 2012) detects shifts in community composition that might occur even when species richness undergoes little to no change. Temporal dissimilarity can be a result of two different phenomena: species replacement (i.e., species replacing one another from one time period to another) and richness change (i.e., species composition of one time period is a strict subset of another). By decomposing dissimilarity into its two underlying components, the role of turnover in community composition change can be disentangled from that of species losses or gains (Baselga, 2010; Carvalho et al., 2012; Villéger et al., 2013). Likewise, more subtle shifts in the distribution of relative species abundance, such as changes in the dominance of species within communities or patterns of evenness in species abundance, rarely strongly affect alpha diversity but can potentially underlie community stability through its relationship to dominant species (Sasaki & Lauenroth, 2011).

Finally, comprehensive characterizations of biodiversity change must include species' intrinsic functional characteristics. Measures that leverage information on species traits, such as functional diversity, provide the opportunity for better insights into community stability by directly linking to species' niches and their functional roles in ecosystems (McGill et al., 2006). For instance, life-history attributes such as reproductive strategies and morphological traits such as body dimensions are better proxies of environmental disturbance regimes than species identities alone (Blanck et al., 2007). In forests (Morin et al., 2014) and grasslands (Fischer et al., 2016), functional diversity contributes to stability via species complementarity in resource use. While others have investigated the relationship between functional diversity and stability in plants (Craven et al., 2018), freshwater macroinvertebrates (Bonada et al., 2007; Péru & Dolédec, 2010; Schriever et al., 2015), or plankton (Korhonen et al., 2010), attempts to incorporate the full complexity of functional diversity change (e.g., functional dissimilarity) remain rare.

Here, we link community stability to change in biodiversity for four animal taxonomic groups: small mammals, ground beetles, fish, and freshwater macroinvertebrates, across a continental scale. We quantify change in biodiversity in two ways: as (1) dissimilarity in community taxonomic and functional composition and species replacement and richness change components of that dissimilarity and (2) change in species' abundance distributions as captured by change in species rank, richness, and evenness. We predict that stable communities will be characterized by strong compositional shifts (i.e., high taxonomic dissimilarity due to species replacement) but without net changes in species richness (i.e., low taxonomic dissimilarity due to species richness change). This is because community properties—such as species richness—are often stabilized by temporal fluctuations in species abundance (Micheli et al., 1999; Schmid et al., 2008; Yachi & Loreau, 1999). We further predict that stable communities will show low functional dissimilarity due to species replacement despite their strong shifts in taxonomic composition. This is because species whose populations fluctuate through time are thought to fill similar functional roles, ensuring stability (Loreau, 2010). Lastly, we expect community stability to be related to shifts in the distribution of relative species abundance. Specifically, we predict that stable communities will show strong changes in species abundance rank because reshuffling of species rank is an inevitable consequence of fluctuating species, but that evenness will remain invariant for most stable communities. Our work illustrates how fine temporal resolution data systematically collected from diverse ecosystems illuminate the general consequences of biodiversity dynamics to stability of communities across biomes.

MATERIALS AND METHODS Data NEON organismal data

Organismal information collected by NEON encompasses a wide range of taxa from terrestrial and freshwater ecosystems. Here, we focus on four taxonomic groups: small mammals, ground beetles (family Carabidae), fish, and freshwater macroinvertebrates. We chose these taxonomic groups because they represent the animal taxa with the longest time series in the nascent observatory.

To enable a robust characterization of community stability and biodiversity change, we included NEON data starting in 2013 and available through 2019 (Appendix S1: Table S1). Sites included in this analysis were located in the contiguous United States, had data for at least 4 years, and were surveyed for at least two (freshwater macroinvertebrates, fish), three (beetles), and four (small mammals) survey bouts within each of those years. We define a “bout” as a single observation event per taxon and site. Differences in the minimum number of bouts required for inclusion in the analysis reflect differences in sampling protocols for each taxon. For fish and freshwater macroinvertebrates, we retained only wadeable stream sites because species composition in lakes is substantially different from those in wadeable streams, and thus is not directly comparable (Irz et al., 2006; Jackson et al., 2001). Twenty-six NEON sites met these criteria for small mammals, 30 for ground beetles, 14 for freshwater macroinvertebrates, and nine for fish. We ultimately analyzed 676 unique bout–site combinations for the 2013–2019 time period for small mammals, 1720 (2013–2019) for ground beetles, 71 (2016–2019) for fish, and 160 (2014–2019) for freshwater macroinvertebrates. Data were downloaded in November 2020 and are available through the Environmental Data Initiative (https://doi.org/10.6073/pasta/7f0e0598132e3fea1bfd36a4257af643; Appendix S1: Table S1). Data from 2020 or 2021 were not included in the analysis due to incomplete sampling resulting from the COVID-19 pandemic.

For small mammals, we included only target species (i.e., those that the NEON protocol was designed to capture) and excluded all recaptures, all individuals not identified to a species level, and those escaped or found dead, for a total of 77 species. A small mammal sampling bout comprises three consecutive (or nearly consecutive) nights of trapping and is based on the lunar calendar, with timing of sampling constrained to occur within 10 days before or after the new moon. The number of bouts per year is determined by site type, and most sites contain four to six bouts per year (minimum = 1, maximum = 8). We used raw abundance and pooled all trapping grids from each NEON site into a single community at bout level because most small mammal trapping grids are separated by only 1–2 km, with the median distance between small mammal grids of ~1200 m (Read et al., 2018). Therefore, small mammal assemblages within a NEON site can thus be considered as a single interacting community.

For ground beetles, we included only target species for a total of 509 species. Sampling for ground beetles occurs biweekly throughout the sampling season (i.e., when the 10-day average low temperature at the site is >4°C), with no single trap being sampled more frequently than every 12 days, thus setting the duration of a bout to 12–14 days. The exact duration of a bout differs slightly among NEON sites, resulting in an effort bias. We accounted for this potential bias by calculating catch per unit effort (CPUE) for each species at a given bout and NEON site as observed abundance normalized by the total number of trap days for all plots in a site sampled during a bout.

Fish were sampled via backpack electrofishing. We excluded all individuals not identified to the species level, which resulted in a total of 79 species included in the analysis. Because the precise stream-reach lengths varied among NEON sites, we standardized abundance by calculating CPUE per sampling reach length at each site. CPUE was calculated as the relative abundance normalized by average electrofishing time. We excluded all bouts with incomplete sampling.

For freshwater macroinvertebrates, specimens not identified to the genus level were excluded, resulting in a total of 568 taxonomically processed morphospecies included in the analysis. Multiple benthic areas were sampled exhaustively at a site with equipment suitable to the habitat. Sampled habitat types included riffles, runs, pools, and step pools; a dominant habitat and a second dominant habitat at each site were sampled. Because freshwater benthos are heterogeneous and different methods were used to sample, we calculated each species' CPUE as sample density (individuals per square meter of stream bottom sampled) at a given bout for each NEON site. Bouts were defined by month of sampling. We refer the reader to the NEON website (https://www.neonscience.org/data-collection/observational-sampling) for details of sampling methodology for each of the taxonomic groups.

Trait data

We explored functional diversity for small mammals and fish only. For small mammals, trait data were obtained from EltonTraits 1.0 (Wilman et al., 2014), which contains data for four traits: body mass, diet, foraging (substrate use) niche, and activity time. The diet trait includes 10 categories: proportions of invertebrates, mammals and birds, reptiles and amphibians, fish, vertebrates unknown, carrion, fresh fruits, nectar and pollen, seeds, and other plant materials in a species' diet. The foraging stratum trait includes five categories: proportion of use of marine, terrestrial or aquatic ground level, scansorial, arboreal, and aerial. Activity time (i.e., nocturnal, crepuscular, and diurnal) was excluded from the analysis because NEON small mammal sampling targets nocturnal species, and we only used nocturnal small mammals.

For fish, trait data were obtained by combining two databases: FishTraits database (http://www.fishtraits.info/; Angermeier & Frimpong, 2011, Xie et al., 2013) and FishBase (https://www.fishbase.de/; Froese & Pauly, 2020). We included the following four functional and morphometric traits: feeding guild, life-history attributes, habitat association, and migratory behavior. The five feeding guild trait categories (dietary) included algae and phytoplankton; macrovascular plants; detritus (partly decomposed organic matter); invertebrates and fish larvae (excluding large crustaceans); and fish, other vertebrates, and large crustaceans. We treated all five categories as binomial (true/false) traits and allowed a single species to be included in more than one feeding guild. Life-history attributes included four numeric traits—maximum possible total length, age at first reproduction, maximum life span, and the mean number of eggs per female—and reflected information on growth, sexual maturity, longevity, and reproductive biology. To account for habitat associations, we selected 10 different types of aquatic substrate categories (muck, clay or silt, sand, gravel, cobble or pebble, boulders, bedrock, aquatic vegetation, organic debris or detrital substrate, and large woody debris) and treated those as binomial traits while permitting a single species to fall into more than one substrate category. Last, to reflect their migratory patterns, we classified each fish species as migratory (diadromous or potamodromous) or residential. The selected traits have a long history of being used in fisheries and fish biology (Blanck et al., 2007; Mims et al., 2010; Villéger, Brosse, et al., 2017). We did not incorporate intraspecific trait variation into measures of functional diversity for either small mammals or fish because many of the traits used in this analysis lack information beyond species mean trait values.

Community stability

We defined community stability as the mean abundance of a given taxonomic group divided by the standard deviation in species' abundance (μ/σ)—that is, the inverse of coefficient of variation (Craven et al., 2018; Hillebrand et al., 2018; Micheli et al., 1999). To evaluate the influence of the temporal duration on the relationship between stability and biodiversity change, for each NEON site, community stability was quantified for an increasing number of bouts (i.e., increasing temporal durations). That is, at each site, we quantified community stability across (2, 3, …, N) consecutive bouts, where N is the number of bouts within a site for a given taxonomic group (e.g., Tilman et al., 2006). We deemed two consecutive bouts to not be a long enough duration to allow sound quantification of community stability and thus retained only those stability measurements that spanned at least three bouts. Ultimately, we ended up with 8359 community stability measurements for small mammals, 48,682 for ground beetles, 182 for fish, and 753 for freshwater macroinvertebrates.

Temporal change in community taxonomic composition

For each taxonomic group, we quantified change in community composition for each pairwise combination of all bouts within each site (following literature on pairwise beta diversity; Baselga, 2013; Marion et al., 2017). First, we quantified the Jaccard-based temporal dissimilarity (beta diversity) in community taxonomic composition (TD_β) using function beta in the R package BAT (Cardoso et al., 2015; Appendix S2: Table S1). We decomposed the Jaccard dissimilarity into its replacement (TD_βRep) and richness change (TD_βRic) components (Cardoso et al., 2014), and also quantified their relative contributions to dissimilarity (cTD_βRep, cTD_βRic). The Jaccard dissimilarity and its replacement and richness change components vary between 0 and 1, with higher values of each measure indicating larger changes across time and larger relative contributions. Mathematically, species replacement and richness change components sum to dissimilarity. We note that the BAT package extends the original framework by Cardoso et al.  (2014, 2015) to include abundance by directly replacing incidence counts with summed abundance (personal correspondence with P. Cardoso), and the resulting dissimilarity is not equivalent to the Jaccard dissimilarity as it would be for incidence data. While the implications of this extension are not obvious and potentially nontrivial, further analysis of the behavior of these measures is outside the scope of this paper. While a similar framework, developed by Baselga (2010), has been extensively used in biodiversity research (Jarzyna & Jetz, 2017; Li et al., 2020), it does not currently allow quantifying dissimilarity based on abundance for both taxonomic and functional diversity (Baselga & Orme, 2012), and was thus not used here. We only retained those pairwise comparisons that corresponded to measurements of community stability (8359 unique bout–bout comparisons for small mammals, 48,682 for ground beetles, 182 for fish, and 753 for freshwater macroinvertebrates)—that is, we removed pairwise comparisons for two consecutive bouts. Having such high replication in pairwise bout comparisons minimizes the potential risk of mischaracterization of compositional change in cases when the two bouts might not be fully representative of the state of biodiversity in time periods surrounding it. We refer the reader to Appendix S3: Figures S1–S3 for an example illustration of the connection between abundance data and metrics of community stability and biodiversity change.

We additionally quantified change in species abundance distributions by calculating change in evenness (TD_Ev), rank (TD_Rank), and richness (TD_Ric) for the same pairwise combinations of bouts as above, using function RAC_change in the R package codyn (Avolio et al., 2019). TD_Ric and TD_Ev vary between −1 and 1, and their positive and negative values of TD_Ev indicate increases and declines, respectively, in richness and evenness over time. Change in species rank was calculated as the absolute value of the mean change in species ranks between two bouts divided by the total number of unique species in both bouts. Species absent from one bout were given the S + 1 rank in that sample, where S is the number of species. TD_Rank varies between 0 and 0.5, and larger values indicate larger change in average species rank across time. Change in richness was calculated as change in richness between bouts standardized by the total number of unique species in both bouts. By standardizing by the number of unique species, TD_Ric provides information complementary to that of TD_βRic, which focuses solely on the contribution of richness change to dissimilarity. Positive and negative values of TD_Ric indicate increases and declines, respectively, in species richness over time. Together, temporal dissimilarity measures (TD_β, TD_βRep, TD_βRic, cTD_βRep, and cTD_βRic) and measures of change in species abundance distributions (TD_Ev, TD_Rank, and TD_Ric) reflect the complete spectrum of the potential differences among communities, in terms of their taxonomic composition, within a site.

Temporal change in community functional composition

For small mammals and fish, we also quantified temporal dissimilarity in community functional composition using the R package BAT (Cardoso et al., 2014, 2015), which assesses functional diversity with ultrametric trees depicting species' functional relationships (i.e., a functional dendrogram). We used BAT to quantify change in community functional composition in the same way as we did for taxonomic composition. That is, for the same pairwise combinations of bouts, we quantified the Jaccard-based temporal dissimilarity (FD_β), decomposed it into its replacement (FD_βRep) and richness change (FD_βRic) components, and quantified their relative contributions to dissimilarity (cFD_βRep and cFD_βRic). To obtain a functional dendrogram, we first calculated multivariate trait dissimilarity using Gower's distance for each pairwise combination of all species. Gower's distance handles quantitative, semiquantitative, and qualitative variables, as well as missing values, and assigns different weights to individual traits (Pavoine et al., 2009). Equal weights were given to each of the traits and to each category within the trait. For small mammals, each category of the diet trait was given a weight of 1/10 (for a total weight of 1 for diet); each category of the foraging niche trait was given a weight of 1/5 (for a total weight of 1 for foraging niche); and body mass was given a weight of 1. For fish, each category of the diet trait was given a weight of 1/5 (for a total weight of 1 for diet); each category of the habitat association trait was given a weight of 1/10 (for a total weight of 1 for habitat association); and total length, age at first reproduction, maximum life span, the mean number of eggs per female, and migratory behavior were all given a weight of 1. Functional dendrograms for all species were then built for small mammals and fish using the unweighted pair group method with arithmetic mean (UPGMA) clustering. UPGMA clustering has the highest cophenetic correlation coefficient (i.e., a measure of how faithfully a dendrogram preserves the pairwise distances between the original unmodeled data points) among the most popular clustering methods (i.e., ward, single, complete, weighted pair group method with arithmetic mean, weighted pair group method with centroid averaging, and unweighted pair group method with centroid averaging clustering methods) and the lowest two-norm index (Mérigot et al., 2010), ensuring the most faithful preservation of the original distances in the dissimilarity matrix. Because using a dendrogram-based approach might in some cases artificially increase the functional distance between species that have similar trait values, we quantified the mean squared deviation (mSD) to assess the congruence between the functional distance as given by Gower's metric and the cophenetic distance on the functional dendrogram (Maire et al., 2015; Villéger, Maire, & Leprieur, 2017). mSD was 0.002 and 0.003 for small mammals and fish, respectively, suggesting the average absolute deviations between Gower's and the dendrogram-based distances were approximately 2% and 3%. The correlation between Gower's distance matrix and the distance matrix based on the functional dendrogram equaled 0.93 and 0.85 for small mammals and fish, respectively, further indicating minimal loss of information.

Community functional composition is typically strongly correlated with community taxonomic composition. To statistically control for this association, we built null models to evaluate whether the different measures of change in community functional composition deviated from those expected given change in taxonomic composition. Null values (i.e., values expected given change in taxonomic composition) were obtained by randomly selecting species from a regional species pool and randomizing abundance for small mammals and CPUE values for fish, while keeping species richness and species abundance distribution at the assemblage level constant (Swenson, 2011). Regional species pools included all species found within a given NEON domain. Randomizations were performed 1000 times. To estimate the direction and magnitude of the deviation of the observed from the expected (randomized) values for each of the measures (FD_β, FD_βRep, FD_βRic, cFD_βRep, and cFD_βRic), we calculated the standardized effect sizes (SES) as (observed value – mean[null distribution])/SD(null distribution) (Swenson, 2011) and denoted those as sFD_β, sFD_βRep, sFD_βRic, scFD_βRep, and scFD_βRic. We refer to those indices corrected for community taxonomic composition as SES-based indices. SES > 0 and SES < 0 indicated change in community functional composition larger and smaller, respectively, than expected by chance given change in community taxonomic composition.

Relationships between community stability and change in community composition

In total, we derived 18 measures of change in community taxonomic and functional composition (Appendix S2: Table S1). We explored the relationships of community stability (the dependent variable) with each measure of change (the independent variable) using a model that followed linear (model type 1; 18 models) and quadratic (model type 2; 18 models) relationships (Appendix S4: Table S1). To account for site-level effects, a random effect on a site was included in all models. The Watanabe-Akaike information criterion (WAIC) was used to select the best-fitting model (Model 1 or 2), with a ΔWAIC >7 used to determine the best model (Burnham & Anderson, 2004). WAIC is preferable to similar alternatives because it averages over the posterior distribution rather than using only point estimates (Gelman et al., 2014). We did not use the coefficient of determination, R², to evaluate goodness-of-fit of each model because R² does not integrate across the full posterior probability and thus ignores uncertainty in parameter estimates (Gelman et al., 2019). Each model was fitted in a Bayesian framework using INLA (Rue et al., 2009) through the R-INLA package for R (https://www.r-project.org/). Note that the problem of multiple comparisons (type 1 error rates) disappears when viewed from a hierarchical Bayesian perspective, and thus, the Bonferroni correction was not necessary (Gelman et al., 2012). All data were processed and analyzed in R (R Core Team, https://www.r-project.org/). All codes are available from Zenodo (https://doi.org/10.5281/zenodo.5541938).

Influences of temporal duration on relationships between community stability and change in community composition

Because temporal duration of the experiment strongly affects the relationship between biodiversity and ecosystem and community stability (Cardinale et al., 2012), we additionally fitted models containing an interaction term between temporal duration and linear (model type 3; 18 models) and quadratic (model type 4; 18 models) effects of biodiversity change on community stability (e.g., Tilman et al., 2006). A random effect on a site was included in Models 3 and 4. We compared Models 3 and 4 using WAIC.

RESULTS Relationships between community stability and change in community composition

Broadly, community stability decreased linearly with dissimilarity for all taxonomic groups (TD_β; Figure 1a, Appendix S4: Table S1), dissimilarity due to richness change (TD_βRic; Figure 1c, Appendix S4: Table S1), and relative contribution of richness change to dissimilarity (cTD_βRic; Figure 1e, Appendix S4: Table S1). Community stability generally increased with increasing dissimilarity due to species replacement (TD_βRep; Figure 1b), though small mammals and freshwater macroinvertebrates showed a quadratic relationship, indicating that community stability is maximized at intermediate levels of species replacement. For all taxonomic groups, community stability increased linearly with contribution of species replacement to dissimilarity (cTD_βRep; Figure 1d, Appendix S4: Table S1). For small mammals and fish, dissimilarity in community taxonomic composition was mostly due to changes in species richness rather than species replacement (Figure 1d,e). For ground beetles and freshwater macroinvertebrates, species replacement contributed more strongly to dissimilarity in community taxonomic composition than richness change (Figure 1d,e).

FIGURE 1. Relationships between community stability and change in community taxonomic composition as given by temporal dissimilarity (a), dissimilarity due to species replacement (b), dissimilarity due to species richness (c), relative contributions of species replacement to dissimilarity (d), and relative contributions of richness change to dissimilarity (e) for four taxonomic groups: small mammals, ground beetles, fish, and freshwater macroinvertebrates. Density plots in (d) and (e) show the frequency of relative contributions of species replacement and richness change, respectively, to dissimilarity in community taxonomic composition. Only relationships whose credible intervals did not overlap zero are shown

We found fewer significant relationships between community stability and changes in species abundance distributions (Figure 2, Appendix S4: Table S1). Nonetheless, for small mammals, community stability was maximized when standardized change in species richness was the lowest (Figure 2a, Appendix S4: Table S1). Small mammals and ground beetles also showed highest community stability at lowest levels of evenness and change in species rank (Figure 2b,c, Appendix S4: Table S1).

FIGURE 2. Relationships between community stability and change in the distribution of species abundance as given by standardized change in species richness (a), change in species abundance rank (b), and change in evenness of species abundance (c) for small mammals and ground beetles. Only relationships whose credible intervals did not overlap zero are shown

For small mammals, community stability increased with increasing SES-based dissimilarity (sFD_β; Figure 3a, Appendix S4: Table S1) and dissimilarity due to species replacement (sFD_βRep; Figure 3b, Appendix S4: Table S1). That is, communities that showed both lower and higher changes in community functional composition than expected given change in community taxonomic composition were more stable. Contributions of species replacement and richness change to functional dissimilarity on average did not deviate significantly from the expectation for either small mammals or fish (Figure 3c,d), and neither showed significant relationships with stability (Appendix S4: Table S1). Expectedly, absolute changes in community functional composition (FD_β, FD_βRep, FD_βRic, cFD_βRep, and cFD_βRic) followed closely those of taxonomic composition (Appendix S3: Figure S4, Appendix S4: Table S1).

FIGURE 3. Relationships between community stability and change in community functional composition as given by the standardized effect sizes (SES) for temporal dissimilarity (a) and temporal dissimilarity due to species replacement (b). Only relationships whose credible intervals did not overlap zero are shown. Density plots in (c) and (d) show the frequency of relative contributions of species replacement and richness change, respectively, to functional dissimilarity in community taxonomic composition

Influences of temporal duration on relationships between community stability and change in community composition

Temporal duration over which community stability and biodiversity change were measured had a strong influence on the relationship between both variables. Community stability decreased with increasing dissimilarity (for all taxa but fish), dissimilarity due to richness change (all taxa), and contributions of richness change to dissimilarity (all taxa) when measured over shorter temporal durations, but this relationship reversed and became positive when temporal duration increased (Figure 4a,c, Appendix S3: Figures S5–S7, Appendix S4: Table S1). When measured over shorter temporal durations, community stability was either maximized at intermediate levels of species replacement (small mammals; Figure 4b) or increased with it and with contributions of species replacement to dissimilarity (ground beetles, fish, and freshwater macroinvertebrates; Appendix S3: Figures S5–S7). Conversely, when measured over longer time spans, community stability was highest at lowest (ground beetles, fish, and freshwater macroinvertebrates; Appendix S3: Figures S5–S7) or minimized at intermediate (small mammals; Figure 4b) levels of species replacement. Likewise, contributions of species replacement to dissimilarity measured over longer temporal duration consistently showed negative relationship with community stability (Figure 4e, Appendix S3: Figures S5–S7).

FIGURE 4. Influence of the temporal duration on the relationship between community stability and change in community taxonomic composition as given by temporal dissimilarity (a), dissimilarity due to species replacement (b), dissimilarity due to species richness (c), standardized change in species richness (d), change in species abundance rank (e), and change in evenness of species abundance (f) for the example taxonomic group, small mammals

When measured across shorter time periods, community stability was highest at lowest levels of richness (small mammals, ground beetles, and freshwater macroinvertebrates), species' rank (small mammals and ground beetles), or evenness (small mammals and ground beetles; Figure 4d–f, Appendix S3: Figures S5–S7) change. These relationships weakened or reversed when stability and change in species abundance distributions were measured over longer time frames (Figure 4d–f, Appendix S3: Figures S5–S7).

Over short temporal durations, small mammal communities that showed both lower and higher functional dissimilarity (sFD_β), functional dissimilarity due to species replacement (sFD_βRep), and higher functional dissimilarity due to richness change (sFD_βRic) than expected given change in community taxonomic composition were more stable (Figure 5). Over longer time frames, stability was highest for communities whose dissimilarity and replacement in functional composition did not deviate from the expectation (Figure 5a,b) and whose functional dissimilarity due to richness change was lower than expected (Figure 5c). Fish did now show significant relationships between stability and change in community functional composition (Appendix S4: Table S1).

FIGURE 5. Influence of the temporal duration on the relationship between community stability and change in community functional composition as given by the standardized effect sizes (SES) for temporal dissimilarity (a), temporal dissimilarity due to species replacement (b), and temporal dissimilarity due to richness change (c) for the example taxonomic group, small mammals

DISCUSSION

We take advantage of the continent-wide organismal data collected by NEON and show that community stability is closely linked to multiple dimensions of change in animal diversity. Consistently across all taxonomic groups analyzed here, community stability generally increased with species replacement and its contribution to overall dissimilarity, supporting our predictions. This corroborates findings from studies suggesting that community-level characteristics are maintained invariant through high replacement of species identities (Allan et al., 2011; Gonzalez & Loreau, 2009; Micheli et al., 1999; Yachi & Loreau, 1999). Despite this general increase, however, we also found that for some taxa community stability declined slightly when highest levels of species replacement were reached. Likewise, stable communities tended to show lower levels of changes in species abundance rank. Together, our work suggests that some degree of compositional reshuffling indeed promotes stable communities, but a complete turnover in community composition might impede stability.

Stable communities also displayed the lowest levels of change in species richness or its contribution to overall dissimilarity, suggesting that sustaining a constant level of diversity is crucial to stability. Likewise, stability tended to be higher for communities characterized by lower change in evenness, further implying that most stable communities maintain not only their richness but also their species abundance distributions despite the relatively high levels of species turnover. Together, our results point to the importance of the maintenance of community-level attributes such as richness or evenness. Further, this constancy can be underpinned by considerable temporal variation in species identities without a detriment to long-term community stability. While we did not explicitly test for asynchronous population dynamics (Loreau & de Mazancourt, 2008), our results are consistent with the idea that temporal fluctuations in species abundance can help stabilize community properties across taxonomic groups and geographies, corroborating theoretical, observational, and empirical studies carried out at more local scales and taxonomic scopes (e.g., Allan et al., 2011; Isbell et al., 2009; Jiang & Pu, 2009; Magurran & Henderson, 2010; Pennekamp et al., 2018; Tilman & Downing, 1994; but see Cardinale et al., 2013; Grman et al., 2010).

Our findings supported our predictions that stability will be highest for communities whose functional dissimilarity and replacement are lower than expected given reshuffling of their taxonomic composition. This is consistent with the notion that long-term community stability relies on fluctuations of species that provide similar functional roles within the ecosystem (Loreau, 2010; Srivastava & Vellend, 2005). Contrary to our predictions, however, we also found communities to be stable at the other end of the functional change spectrum—that is, when functional dissimilarity and replacement were higher than expected given reshuffling in species identities. This agrees with some literature on stability in plant biomass. For example, Allan et al. (2011) found that strong functional turnover promotes stability because species replacing one another provide high biomass production despite belonging to different functional groups. In our system, we reconcile these seemingly conflicting results by noting that communities with highest functional dissimilarity due to replacement are also ones with low reshuffling in species identities (Appendix S3: Figure S8). Perhaps high functional uniqueness of fluctuating species compensates for lower taxonomic reshuffling to drive stabilization of community properties.

Potential implications of our methodological choices warrant discussion. First, we disentangled the effects of functional diversity from those of richness by obtaining measurements independent of change in taxonomic composition. Without such statistical correction, functional turnover might be overestimated when fluctuating species have strong functional redundancy and underestimated when they have strong functional uniqueness. Results from previous studies that did not correct for taxonomic composition (e.g., Allan et al., 2011) should be reconsidered in light of this potential biasing factor. Secondly, because the dimensionality of potential ecological functions is larger in animals than in plants (Pigot et al., 2020), it is not yet fully understood which animal functional traits are strongly mechanistically linked to abundance, and thus community stability (Violle et al., 2007). Though the selected traits have a long history of being used in both fish (Blanck et al., 2007; Mims et al., 2010; Villéger, Brosse, et al., 2017) and small mammal (Kohli et al., 2021) studies and have been vetted for their importance in providing vital functions in ecosystems, a selection of alternative traits might have led to a different inference.

Ecological communities are perturbed by environmental fluctuations over a range of temporal scales (Gonzalez et al., 2020). We found significant effects of temporal duration on community stability itself and the relationship between community stability and biodiversity change, suggesting that the processes maintaining community stability are temporally scale-dependent. Notably, the broad relationships uncovered here held mostly at shorter temporal durations and weakened or reversed when stability and change were measured over longer time spans. At longer time spans, stable communities were characterized by higher taxonomic dissimilarity, stronger changes in richness, and more pronounced species replacement and change in species abundance rank. Given that the range of environmental conditions increases with increasing time span (Gonzalez et al., 2020), one would indeed expect species fluctuations to become more pronounced with increasing temporal extent, potentially leading to higher community stability. In fact, our initial prediction on high stability being linked to reshuffling in species abundance rank—though not supported for the broad dataset or at shorter time spans—gained support at longer temporal durations, suggesting that the role of fluctuations in species identities and abundance increases with time span. Our findings are consistent with a growing body of research suggesting that the effects of biodiversity on stability grow stronger with increasing temporal duration (Cardinale et al., 2007; Stachowicz et al., 2008; Tilman et al., 2001). We caution, however, that community stability and biodiversity change for taxonomic groups that have longer generation times might not be currently fully captured by NEON data and encourage that long time series data are used to further explore these relationships. Moving forward, NEON data will provide distributed monitoring of the relative abundance over decadal timescales, which will better capture fluctuations and long-term trends in populations of all taxonomic groups to better manage threats to these species (Degrassi et al., 2019).

Significant temporal turnover in community composition has been reported across terrestrial and marine ecosystems and annual to decadal time spans (Dornelas et al., 2014; Jarzyna & Jetz, 2017; Korhonen et al., 2010; Li et al., 2020). We provide initial insights into how to place findings from these studies in the context of community stability. High reshuffling of taxonomic composition reported by others (Dornelas et al., 2014) may not necessarily imply loss of community stability through time, though the implications might differ depending on change in community functional composition and the temporal scale of the investigation.

We defined community stability as stability of aggregate species' abundance (e.g., Micheli et al., 1999). However, other definitions of stability such as resistance (i.e., capacity to withstand change) or resilience (i.e., rate at which a system returns to its equilibrium; Donohue et al., 2016) can be applied to ecosystem-level attributes (e.g., biomass; Hillebrand et al., 2018; Hoover et al., 2014; Pennekamp et al., 2018) rather than—as we did here—community-level properties. Temporal community- and ecosystem-level stability, resistance, and resilience might show conflicting relationships with one another and with their putative drivers (Hillebrand et al., 2018), and these relationships can be significantly modified or decoupled by environmental perturbations (Donohue et al., 2013). For example, community-level resistance and resilience in freshwater mesocosm experiments were not related to one another, but their ecosystem-level equivalents showed a negative association (Hillebrand et al., 2018). Likewise, high compositional resistance promoted temporal stability of communities, but ecosystem resistance was not associated with temporal stability ecosystem-level properties (Hillebrand et al., 2018). The current temporal extent of the NEON data products and limited access to biomass data limits the investigation of community- or ecosystem-level resistance, resilience, or reliability, or how these alternative measures of stability are related to temporal change in community composition across our study systems. This study nonetheless provides an initial observatory- and continent-wide assessment of links between community stability and change in biodiversity, providing an important baseline for continued data collection across all NEON sites, which will enable pursuing these and other questions related to stability of communities and ecosystems to environmental disturbances as the observatory matures.

ACKNOWLEDGMENTS

The authors thank two anonymous reviewers for comments on the previous versions of this manuscript. The National Ecological Observatory Network is a program sponsored by the National Science Foundation and operated under cooperative agreement by Battelle Memorial Institute. This material is based in part on work supported by the National Science Foundation through the NEON Program. The authors acknowledge support from the NSF Award #1906144 to attend the 2019 NEON Science Summit. Additionally, the authors acknowledge support from the NSF DEB 1926598 to Marta A. Jarzyna, NSF DEB 1926341 to Jalene M. LaMontagne, NSF DEB 1926568 to Sydne Record, NSF DEB 1926567 to Phoebe L. Zarnetske, and the US Department of Energy Computational Sciences Graduate Fellowship DE-FG02-97ER25308 to Kari E. A. Norman.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

Code and data (Norman et al., 2021) are available from Zenodo: https://doi.org/10.5281/zenodo.5541938 and are additionally archived with the Environmental Data Initiative: https://doi.org/10.6073/pasta/7f0e0598132e3fea1bfd36a4257af643