Making sense of multivariate community responses in global change experiments

Abstract

Ecological communities are being impacted by global change worldwide. Experiments are a powerful tool to understand how global change will impact communities by comparing control and treatment replicates. Communities consist of multiple species, and their associated abundances make multivariate methods an effective approach to study community compositional differences between control and treated replicates. Dissimilarity metrics are a commonly employed multivariate measure of compositional differences; however, while highly informative, dissimilarity metrics do not elucidate the specific ways in which communities differ. Integrating two multivariate methods, dissimilarity metrics and rank abundance curves (RACs), have the potential to detect complex differences based on dissimilarity metrics and detail the how these differences came about through differences in richness, evenness, species ranks, or species identity. Here we use a database of 106 global change experiments located in herbaceous ecosystems and explore how patterns of ordinations based on dissimilarity metrics relate to RAC-based differences. We find that combining dissimilarity metrics alongside RAC-based measures clarifies how global change treatments are altering communities. We find that when there is no difference in community composition (no distance between centroids of control and treated replicates), there are rarely differences in species ranks or species identities and more often differences in richness or evenness alone. In contrast, when there are differences between centroids of control and treated replicates, this is most often associated with differences in ranks either alone or co-occurring with differences in richness, evenness, or species identities. We suggest that integrating these two multivariate measures of community composition results in a deeper understanding of how global change impacts communities.

INTRODUCTION

Given the pressing concern of global change impacts on ecological communities (Franklin et al., 2016), there is an explosive growth of experiments studying global change. For example, globally distributed experiments alone (e.g., Borer et al., 2014; Henry & Molau, 1997; Yahdjian et al., 2021), such as the Nutrient Network, Drought-NET, and the International Tundra Experiment, result in 100s of experiments worldwide. A key aspect of global change experiments—to understand the impacts of global change on biodiversity—presents a challenge because studying biodiversity is complex, and there is no agreed upon best approach (Avolio et al., 2015; Hillebrand et al., 2018; Magurran, 2016; McGill et al., 2015). Two key aspects of biological community data are the species present and their associated abundances (Avolio et al., 2019; Foster & Dunstan, 2010). Univariate measures, such as species richness, boil all this complexity down to simply how many species are present, while most multivariate measures consider species identities and their abundances (Foster & Dunstan, 2010). Not surprisingly, multivariate measures have been shown to be more sensitive than richness to detecting experimental impacts on plant communities (Komatsu et al., 2019).

One common multivariate method is to use dissimilarity metrics, such as Bray–Curtis dissimilarity, to compare biological communities (Buckley et al., 2021; Legendre & Anderson, 1999; Legendre & Legendre, 2012). A benefit of this approach is that it summarizes complex community data into interpretable numbers—creating a scale with which to compare the similarity of biological communities. Another benefit is that these dissimilarity measures can be visualized in ordinations, typically called multivariate community space (Figure 1a) (Legendre & Legendre, 2012; McCune & Grace, 2002). When comparing two treatments, such as control and treated replicates, there are two aspects that can be reported: the distance between the centroids (representative of overall differences in community composition) (Anderson, 2001) and the dispersion of points (replicates) around the treatment centroid (representative of community variability among replicates) (Anderson, 2006) (Figure 1a).

FIGURE 1. Different methodological approaches to study community compositional differences. (a) Dissimilarity-based measures assess differences between centroids and the dispersion of points around the centroid. (b) Rank-abundance curve-based measures compare the richness, evenness, ranking of species, and species composition differences. Each colored point represents a species and S denotes the species richness. Figure modified from Avolio et al. (2019).

Dissimilarity-based measures are useful for visualizing and detecting community changes over time or difference over space. However, ordinations do not clearly detail how community composition is different between treatments, and RACs have been suggested to be a link between dissimilarity measures and specific community composition differences (Avolio et al., 2015, 2019). RACs are a type of species abundance distribution (McGill et al., 2007) that plots a species' abundance versus its abundance rank in a community (Figure 1b). Most work on RACs has been on quantifying the shape of the curve (e.g., Ulrich et al., 2010); however, important information is gained when attention is also paid to species identity (Mac Nally, 2007). When accounting for species identity, RACs can be used to compare two communities and detail all possible compositional differences (Figure 1b) (Avolio et al., 2019). Two communities can differ in richness (only measures the number of species, and ignores species identity and abundances), evenness (only measures distribution of abundances among species, ignores species identity), species identity (only measures species shared, and ignores abundances), and ranks (considers both species identities and their abundances) (Figure 1b). Thus, integrating RAC-based differences alongside dissimilarity-based measures of difference might be a powerful way to understand the impacts of experimental treatments on ecological community composition.

Avolio et al. (2015) presented evidence of six patterns, hereafter scenarios, of dissimilarity-based multivariate differences (Figure 2) that are found in global change experiments. These six scenarios consider both differences in centroid means of replicates in a treatment and dispersion of replicates around the centroid in a treatment (Figure 1a). Next, Avolio et al. (2015) suggested using RACs to understand what about community composition was different (Figure 1b), and outlined hypotheses for each of the six scenarios on how the communities differ between control and treatments (Table 1; Figure 2). Avolio et al. (2015) hypothesize that Scenario 1 (no difference in composition or dispersion) corresponds with no difference in richness, evenness, ranks, or species identities. They hypothesized Scenario 2 (no difference in composition but an increase in dispersion) and Scenario 3 (no difference in composition and a decrease in dispersion) were being driven by changes in rare species only, because compositional differences must be subtle for there to be no differences in centroids. They hypothesized that Scenario 4 (a difference in composition with no change in dispersion) could not be the result of differences in richness or evenness alone, because differences in centroid means would necessitate differences in ranks or species identities. Finally, they hypothesized Scenario 5 (a difference in centroids and an increase dispersion) and Scenario 6 (a difference in centroid means and a decrease in dispersions) required differences in ranks or species identities but that an increase in dispersion would occur when different species become highly abundant in replicates, and a decrease in dispersion would occur when the same species became highly abundant in replicates (Avolio et al., 2015).

FIGURE 2. Hypotheses of what mechanisms of compositional difference underlie patterns of dispersion and centroid differences. This figure is modified from Avolio et al. (2015); each scenario is numbered in red and matches the text accordingly.

TABLE 1 Hypothesized relationships between RACs and dissimilarity-based community patterns between control and treated replicates (differences in centroid and dispersion).

Scenario	Multivariate community pattern	RAC differences
		Richness	Evenness	Ranks	Species identity
Scenario 1	= composition, = dispersion	No	No	No	No
Scenario 2	= composition, ↑ dispersion	Yes	Yes	Yesa	Yes
Scenario 3	= composition, ↓ dispersion	Yes	Yes	Yesa	Yes
Scenario 4	Δ composition, = dispersion	No	No	Yesb	Yes
Scenario 5	Δ composition, ↑ dispersion	No	No	Yesc	Yes
Scenario 6	Δ composition, ↓ dispersion	No	No	Yesd	Yes

Note: For each dissimilarity-based scenario, Avolio et al. (2015) hypothesized whether differences in one aspect of the RAC alone would likely give rise to that pattern, and our unit of consideration is each replicate. Please note that these differences are not mutually exclusive, meaning multiple aspects can be different simultaneously. Abbreviation: RAC, rank abundance curve. a Will occur only if the response is limited to rare species. b Will occur only if the same species respond similarly in all replicates. c Will only occur if different rare species become dominant across replicates. d Will only occur if the same species become dominant across replicates. Source: Reprinted with modifications from Avolio et al. (2015).

Our goal is to explore how dissimilarity-based scenarios of community responses to experimental treatments are related to differences in RACs by testing the hypotheses presented in Avolio et al. (2015). For these analyses, we used the CoRRE database (Community Responses to Resource Experiments; www.corredata.weebly.com), which enabled us to compare control versus treatment community compositional differences across 106 experiments worldwide.

METHODS Database

The CoRRE database is a collection of global change experiments that manipulate at least one plant resource: CO₂, water, and nutrients (see Appendix S1: Table S1 for a list of experiments). Not all treatments in each experiment are required to be resource manipulations, and thus, there are non-resource treatments such as warming and herbivory. Additionally, each experiment in the database has species abundance data for every species recorded in each replicate, at least 3 years of data, and a minimum of four replicates per treatment; however, data collection methods vary by experiment. All replicates are assigned as either control or a treatment. The experiments included in this analysis ranged from 183 to 1568 mm mean annual precipitation and from −12 to 22°C mean annual temperature. For these analyses, we made control–treatment comparisons for each treatment within an experiment for each year of data, resulting in 2831 control–treatment comparisons. Komatsu et al. (2019) and Avolio et al. (2021) reported how plant communities responded to global change drivers; here we focus on how RAC-based measures are related to dissimilarity-based measures of community differences.

RAC-based measures of community differences

There are four RAC-based measures of community compositional differences (Figure 1b)—differences in richness, evenness, rank, and species identity—which are detailed in Avolio et al. (2019). Briefly, richness difference is a proportional measure that ranges from −1 to 1 and is the difference in the number of species between the control and treated plots divided by the total number of species in the control and treated plots. Evenness difference ranges from −1 to 1 and is the difference in evenness measured using Evar (Smith & Wilson, 1996) between control and treated plots. “Rank difference” ranges from 0 to 0.5 and is the average difference in rank for each species between control and treated plots divided by the total number of species in treated and control plots. Species identity difference ranges from 0 to 1 and is the difference in species identity (or turnover) between two samples only, removing richness differences caused by species nestedness (Baselga, 2010; Carvalho et al., 2012).

Statistical analysis

All analyses were performed in R (R Core Team, 2002) using an α of 0.05, using data archived on Enivronmental Data Initiative (Avolio, 2022a) and code archived on Zenodo (Avolio, 2022b). We chose to conduct analyses for each year separately because we are interested in the relationships between measures of community composition rather than global change effects per se, and this way we have the greatest possible control–treatment comparisons to detect patterns. To investigate patterns of community differences using dissimilarity metrics, we calculated the difference in community composition and dispersion for each treatment in an experiment between control and treatment plots for all years of an experiment using Bray–Curtis dissimilarity in the multivariate_difference() function in library(codyn) (Avolio et al., 2019; Hallett et al., 2020). We also determined whether there were significant differences between control and treatment centroids using the adonis() function and differences in dispersion around the centroids using the betadisper() function in library(vegan) (Oksanen et al., 2019), again based on Bray–Curtis dissimilarity. We corrected for multiple comparisons within each dataset, adjusting the p-value for the total number comparisons including four RAC difference measures, number of years, and number of treatments using Benjamini–Hochberg correction (Benjamini & Hochberg, 1996) with p.adjust() in R. For example, if an experiment had four treatments plus a control and five years of data, we corrected for 80 comparisons (4 RAC distance measures × 4 treatment–control comparisons × 5 years). When there was a significant difference in dispersion between controls and treated plots, we noted whether the treatment had either greater or less dispersion than control plots, using the output from multivariate_difference(). Depending on whether there was a significant difference in centroids (p value from adonis <0.05) or dispersion (p value from betadisper <0.05) (Figure 2), and whether dispersion increased or decreased, we classified patterns of differences in centroids and dispersion according to the six scenarios outlined in Avolio et al. (2015). Thus, the different scenarios only reflect patterns based on dissimilarity-based multivariate measures.

To link these dissimilarity-based metrics of differences between control and treatment plots with differences in RACs, we calculated the RAC differences (richness, evenness, rank, and species identity) between control and treated plots for each year of an experiment using the RAC_difference() function in library(codyn). Specifically, we were interested in whether differences between control and treated plots were significantly different from what would be expected from just natural heterogeneity, which we calculated as differences among control plots only. To determine the difference in control and treated plots, each control plot was compared with each treatment plot for all treatments in the experiment; because most experiments were not blocked, we did not account for blocks in this analysis. We then took the average of the control–treatment differences for each control plot in each treatment. Thus, if there were six control plots and six treatment plots, for each control plot the average difference between it and each treatment plot was calculated, resulting in n = 6. We repeated this analysis, but for the control plots only, as a comparison for natural heterogeneity among control plots for each year of an experiment. For each control plot, we calculated the average of how different it was from all other control plots. Thus, similar to the control–treatment comparisons, if there were six control plots, for each control plot the average difference between it and each control plot was calculated, resulting in n = 6. We then ran t tests between control–treatment difference versus control–control difference for each treatment for each year of the experiment. We again corrected for multiple comparisons for each dataset using the Benjamini–Hochberg method. We binned these responses as significant (p < 0.05) or not, with significant effects demonstrating that the control–treatment differences were more different than the control–control differences. We first looked at whether each RAC-based measure alone was different for each scenario and then combined these measures to encompass all the possible combinations. Because Avolio et al. (2015) hypothesized that differences in richness and evenness alone could not result in a difference between centroids, we first combined these two measures only, and then combined them with differences in ranks and differences in species identities. Finally, we combined differences in ranks and species identities with differences in richness and/or evenness.

We hypothesized that rare species drive changes in dispersion when there is no concurrent change in community centroids (Scenarios 2 and 3). To test this hypothesis, we re-ran the adnois(), betadisper(), and multivariate_difference() functions excluding rare species. Rare species were defined as those species that have less than 5% relative cover, averaged over all control plots for all years of the experiment. If the hypothesis is correct, the percent of studies that fit those scenarios would decrease with rare species removed. Finally, to test whether the same or different species were becoming more abundant in Scenarios 5 and 6, we used the abundance_difference() function in codyn to compare the treatment with the control plots. For each species, we calculated the difference in relative abundance between treatment and control plots. We then calculated the proportion of treatment plots where a species increased by at least 20% in their relative abundance relative to the control plots (similar results were obtained with a 10% increase). If the same species was increasing in many/most plots, similar to the hypothesis for Scenario 6, this would be a high proportion. In contrast, if different species were increasing in different replicates, similar to the hypothesis for Scenario 5, this would be a lower proportion. We performed a t test to study the difference in proportion of treatment plots a species became dominant in Scenario 5 versus Scenario 6, expecting it to be higher in 6 than 5.

RESULTS Linking multivariate scenarios to RAC-based measures Scenario 1 (= composition, = dispersion)

This scenario was the most common scenario (Figure 3). We hypothesized that when no community differences were detected based on multivariate measures, there would also be no differences in richness, evenness, ranks, or species identity of control and treated replicates. For the most part, we find support for this hypothesis, as we detected the fewest control–treatment differences with measures based on RACs (Figure 4). The most common control–treatment differences were in richness or evenness only (Figure 4). However, there was some mismatch, as we detected significant control–treatment RAC differences in 20% of comparisons, suggesting that some control–treatment differences are occasionally not detected by multivariate measures.

FIGURE 3. For each scenario, the percentage of times it was observed out of all 2831 control–treatment comparisons.

FIGURE 4. For each RAC difference measure alone and in combination, the percent of studies that found differences between control and treated replicates based t tests. We binned community differences as differences in richness (R), evenness (E), rank (Ra), and species identity (S) alone, then differences in richness and evenness only (R + E), rank differences co-occurring with any combination of richness and evenness (Ra + R and/or E), species identity difference co-occurring with any combination of richness and evenness (S + R and/or E), and finally, ranks and species identity difference co-occurring with any combination of richness and evenness (Ra + S + R and/or E). Scenarios (1–6) are denoted in the upper left corner of each panel. RAC, rank abundance curve.

Scenarios 2 (= composition, ↑ dispersion) and 3 (= composition, ↓ dispersion)

Scenarios 2 and 3 rarely occurred (Figure 3). We hypothesized that Scenarios 2 and 3 were driven by rare species. When we removed rare species, there was no change in the percent of studies that fit Scenario 2 and an increase to 7% of the number of control–treatment comparisons that fit Scenario 3. Thus, we found no evidence that these scenarios are being driven by rare species; in fact, rare species appear to dilute the ability of betadisper() to pick up significant decreases in treatment plot dispersions (Scenario 3). For Scenario 2, most cases involved differences in species ranks either alone or co-occurring with other measures (Figure 4). For Scenario 3, most cases exhibited differences in richness with a few others showing differences in species identities co-occurring with other measures (Figure 4). Interestingly, while both Scenarios 2 and 3 reflect community differences between treated and control replicates, again, there was a mismatch between what was detected by RAC measures. When dispersion was significantly different, measures based on RACs detected significant differences only 49% of the time for Scenario 2 and only 27% of the time for Scenario 3.

Scenario 4 (Δ composition, = dispersion)

Scenario 4 was the second most common multivariate scenario of control–treatment differences detected (Figure 3). When there were only differences in centroid means between control and treated replicates, there were most often differences in ranks co-occurring with richness or evenness differences, and less commonly with any one measure alone (Figure 4). Generally, there was agreement between multivariate and RAC measures of community differences; RAC differences were detected in 85% of the studies with significant control–treatment centroid differences.

Scenarios 5 (Δ composition, ↑ dispersion) and 6 (Δ composition, ↓ dispersion)

Scenarios 5 and 6 were also uncommon (Figure 3). Scenario 5 was most associated with differences in ranks co-occurring with other RAC differences, while Scenario 6 was most associated with differences in ranks and species identity differences co-occurring with differences in richness and evenness (Figure 4). There was also the strongest agreement between detecting significant differences based on dissimilarity measures versus RAC-based measures; RAC-based measures detected at least one aspect of control–treatment differences: 95% of the time for Scenario 5 and 92% of the time for Scenario 6. Avolio et al. (2015) hypothesized that the difference between Scenarios 5 and 6 was dependent on the proportion of replicates in which the same species increased in abundance. For Scenario 5, when there is divergence, it was hypothesized that a species would become abundant in a small proportion of treatment replicates. In contrast, it was hypothesized that in Scenario 6, when there is convergence, the same species would become abundant in many of the treatment replicates. We found evidence for both hypotheses, where similar species became dominant in only 25.7% ± 0.3% of replicates in Scenario 5, whereas similar species became dominant in significantly more replicates (59% ± 2.2%) in Scenario 6 (t = −14.75; p < 0.001).

DISCUSSION

Gaining insights into how biodiversity is impacted by global change is a top research priority (O'Connor et al., 2021), and numerous studies have pointed to the necessity of including both species identity and their abundances in our attempts to understand community changes (Gotelli & Colwell, 2001; Heip et al., 1998; Hill, 1973; Wilsey et al., 2005). Experiments are an important tool for studying global change (Schlesinger, 2006), and here, we compare and integrate two different methods for studying multivariate community compositional differences between control and treated replicates—dissimilarity metrics (and their ordinations) (Anderson et al., 2006; Bray & Curtis, 1957) and RACs (Avolio et al., 2019; Foster & Dunstan, 2010). We assessed relationships between scenarios of community compositional differences based on dissimilarity metrics with compositional differences based on RACs. Overall, we find support for the hypotheses put forth by Avolio et al. (2015). For example, differences in richness and evenness alone rarely resulted in differences in centroids. Further, we found intuitive links between dissimilarity-based scenario of control–treatment differences and RAC-based measures of community differences.

When comparing control and treatment replicates, if there is no difference between their centroids (Anderson et al., 2008), one would conclude there is no compositional difference between control and treatment communities (Scenarios 1–3; Figure 2). Overall, when there was no difference between control–treatment centroids, we found many fewer control–treatment differences based on RACs. Differences in dispersion were associated with all RAC-based measures of community composition differences except species differences alone. Interestingly, we found much fewer significant control–treatment differences based on RAC-based measures, demonstrating that dispersion is a subtle aspect of community difference. Detecting differences in variance is important for ecological properties (Benedetti-Cecchi, 2003), and dissimilarity metrics appear to be most sensitive to detecting these differences.

When there are differences in control–treatment centroids (Scenarios 4–6; Figure 2), one can conclude that there are compositional differences between these communities. Avolio et al. (2015) hypothesized that rarely would richness or evenness alone result in differences in centroids, which we confirmed. Richness has been recognized as an insensitive measure of community compositional differences (Magurran, 2016), with evenness less investigated. We found differences in centroids were mostly associated with differences in ranks either alone or co-occurring with other measures of differences, such as richness, evenness, or species identity. Shifts in the rank order of species in a community is an underappreciated mechanism by which community composition can change (Jones et al., 2017), in comparison to species identity differences and turnover. In the context of global change, while species migrations will occur as species track their more ideal environmental conditions (Neilson et al., 2005), for terrestrial plant species, this will be a slow process (Aitken et al., 2008; Jump & Penuelas, 2005). Instead changes in local competitive conditions will result in different species hierarchies as species increase or decrease in dominance. Thus, paying attention to shifts in ranks may be a more important process on shorter time scales than shifts in species identities. By contrast, we found that differences between control and treatment communities were rarely the result of species identity differences. Dissimilarity metrics have been used as a proxy for compositional differences, or turnover (Magurran et al., 2010); however, our findings here suggest that compositional differences inferred from ordinations based on Bray–Curtis dissimilarity rarely reflect species identity differences and instead reflect differences in ranks (Avolio et al., 2019). Instead, different indices that focus on turnover may be necessary to directly study species identity differences (Vellend, 2001).

The number of ways to study community composition is large and increasing. Here, we demonstrate that ordination patterns based on dissimilarity metrics, a common approach to comparing two or more communities, intuitively reflect underlying community compositional differences. This link between RACs and dissimilarity measures will hopefully yield insights into mechanisms underlying compositional responses of communities to global change. For example, knowing that differences in ranks are an important process, studies can then investigate which species have different ranks and try to pinpoint the mechanisms associated with these differences. One method may be to use rank clocks to visualize which species are undergoing the most drastic shifts in abundance (Collins et al., 2008). Our work stresses the importance of utilizing multiple measures of community differences to holistically study biological communities from a variety of angles. This will provide more detailed management targets and predictions for cascading effects of altered community composition on ecosystem function.

ACKNOWLEDGMENTS

We thank the LTER Network for funding our synthesis working groups in 2012 (NSF EF 1545288 to Meghan L. Avolio and Kimberly J. Komatsu) and 2016 (NSF EF 0553768 to Kimberly J. Komatsu, Meghan L. Avolio, and Kevin R. Wilcox). Importantly, we are immensely grateful to the researchers who provided data for this manuscript and their research teams who assured that the long-term integrity of these global change experiments are/were maintained.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

Data (Avolio, 2022a) are available from the Environmental Data Initiative Data Portal: https://doi.org/10.6073/pasta/200b9c492419ec14c43a68d4411ad076. Code (Avolio, 2022b) is available from Zenodo: https://doi.org/10.5281/zenodo.6689891.