Measuring change in biological communities: multivariate analysis approaches for temporal datasets with low sample size

Count the number of taxa, the number of samples and the number of empty samples and missing data points

As some multivariate methods are more susceptible than others to the influence of rare taxa, the consideration of the ratio of the number of taxa to the number of samples in a dataset is an important initial exploratory step (Lepš & Šmilauer, 2003). It is common for long-term datasets to contain missing values; these values should be denoted as missing (i.e. ‘NA’ in R) and not given zero or other numerical values. The importance of missing values and how to deal with them during data analysis depend on the number of missing data points and how they are distributed across the dataset. Multivariate analyses based on certain compositional dissimilarity measures cannot be applied to datasets containing any samples where no taxa are present (McCune & Grace, 2002; Legendre & Legendre, 2012). Many distance measures cannot be computed if there are empty samples or missing values. Deleting entire rows or estimating missing data so that the dataset can be analysed may be an option in some cases (McCune & Grace, 2002). However, other measures, such as the Gower coefficient, can be calculated while ignoring comparisons among samples with missing data (Legendre & Legendre, 2012; coefficient S15: 278). Comparing results based on presence-absence data with those generated using relative abundance data allows us to evaluate the role of taxon abundances in generating temporal patterns (Angeler & Johnson, 2012). In contrast, empty samples may be important for analyses where we are interested in understanding responses to disturbance and local extinctions.

Ensure that relative abundances are comparable

Estimating abundances of each population within the community has its own important considerations, such as variation in taxon detection probabilities, which have been addressed elsewhere (e.g. animals: Gaston & McArdle, 1994; amplicon data: Amend, Seifert & Bruns, 2010). However, once the taxon-sample matrix has been collated, abundance data from different taxa can be compared. If abundances were measured on different scales, the data should be standardised or transformed prior to analysis, for example in the case of aquatic communities, macroinvertebrate abundances are likely to have been estimated using different volumes of water to macroinvertebrates (Legendre & Legendre, 2012). In such cases, taxon abundance values within each sample can be divided by the total abundance within the sample, or can be standardised across all samples so that each taxon has a mean of zero and standard deviation of one (Borcard, Gillet & Legendre, 2011). Alternatively, a distance measure that incorporates a standardisation can be used, for example, Hellinger distance (Legendre & Gallagher, 2001; Anderson et al., 2011; Legendre & Legendre, 2012).

Examine taxon frequency distributions and response curves

Taxon frequency distributions can be visualised by plotting a histogram for each taxon of the number of individuals or their relative abundance across samples. This is important as some multivariate analyses rely on the assumption that these frequency distributions take particular shapes such as normal distribution or a linear response across an environmental gradient that samples were collected from; for example, canonical correspondence analysis (CCA) assumes that taxa have optima along environmental gradients (‘unimodal’ responses). Examination of these distributions can reveal whether taxa are responding strongly to environmental gradients in the dataset and whether the dataset conforms to analysis assumptions. Data transformations, such as log or square-root, can be applied to normalise abundance data, if required.

Determine the ratio of temporal to spatial variation in composition

The amount of spatial and temporal replication is important in determining the method of analysis because few samples with many measurements over time would be analysed differently from many samples measured only a few times. The ratio of the number of independent replicates to the number of explanatory variables is also important for temporal analyses because it affects the degrees of freedom available for statistical testing (Johnstone & Titterington, 2009). If the data have a nested structure this needs to be accounted for in the analysis, particularly if one wants to report significance values. For example if plant species were recorded in quadrats that were nested within transects that were nested within sites, then quadrats are not fully independent replicates. Similarly, information from technical replicates (such as PCR products) should be kept separate at this stage because they provide useful information about variability at different scales in the dataset. Thus, it is inadvisable to average or pool data until the significant sources of variation are understood. For instance, if there is a high degree of spatial variation among samples in the dataset, this can alter analysis choices or how results are presented and interpreted (Thrush et al., 2008). The relative amount of change in taxonomic composition can be estimated from the differences in sample scores on detrended correspondence analysis (DCA) ordination axes, for example, the ‘gradient length’ statistic quantifies the turnover in taxonomic composition (Lepš & Šmilauer, 2003).

Method 1: pairwise dissimilarity-based approaches

Pairwise dissimilarity methods are the most common way biologists quantify the difference in the identities of taxa occurring in two samples (Buckley et al., 2021). The result is a single number that measures the distance in compositional space between the two samples. Dissimilarity can also be represented as similarity, that is, 1—dissimilarity (Whittaker, 1960; Anderson et al., 2011; Magurran, 2013). When applied to multiple samples, these methods result in a matrix of all possible comparisons (a triangular dissimilarity matrix).

A community dataset may be an abundance matrix (number of individuals), relative abundance matrix (e.g. percent cover, sequence abundances), rank-abundance matrix (ordinal abundance ranks), or a taxon-occurrence (presence-absence) matrix. These different response variables can influence the choice of analysis methods, including which dissimilarity measure to use, so it is important to explore different metrics to determine which method gives the best representation of differences in composition. For instance, the commonly-used Euclidean distance measure is not usually appropriate for taxon abundance data, which are often characterised by non-linear relationships among taxa, abundance values measured on different scales, or by the ‘double-zero problem’, which is where samples that are missing the same taxa have lower dissimilarity (Borcard, Gillet & Legendre, 2011).

Different authors often use different names for the same dissimilarity measure, so we recommend Legendre & Legendre (2012 Table 7.4: 324) for a comprehensive and consistent guide. However, other definitive sources for overviewing and comparing dissimilarity indices are provided by Anderson et al. (2011), Legendre & De Cáceres (2013), and Magurran (2013). The development of dissimilarity indices is an active area of research, including those specifically applicable to temporal data (Podani & Schmera, 2011; Baselga, 2012; Shimadzu, Dornelas & Magurran, 2015) and some contributions aid in improving dissimilarity measures for use in certain circumstances, such as imperfect detection of taxa (Chao et al., 2006) or where phylogenetic information is available from sequencing that has been used to determine microbial taxa in samples (unifrac; Lozupone & Knight, 2005; Lozupone et al., 2007).

Method 2: rank abundance

Rank abundance distributions (RADs, a.k.a. ‘rank abundance curves’), can provide relatively sophisticated visualisations of community dynamics, and rank abundances can be used to quantitatively measure the shift in taxon ranks within samples taken at different time points (Table 2; Code S4; Avolio et al., 2015, 2019; Hallett et al., 2016). RADs are plots of the relative abundances (on the y-axis) of all species in a community ranked (on the x-axis) from highest to lowest. They are closely related to k-dominance curves, in which species ranks are shown on the x-axis, but cumulative proportional abundance (dominance) is plotted on the y-axis (Lambshead, Platt & Shaw, 1983). If samples taken at different times are plotted separately, the set of RADs can illustrate changes in community structure over time. RADs also can be used to explore community dynamics if a dissimilarity measure, such as Bray–Curtis, is calculated for the species’ ranks at each time point. The resulting temporal dissimilarity of the RADs can be explored by plotting rank dissimilarity against time lag or by ordinating rank dissimilarities in the same way that spatial compositional dissimilarities are examined. For example, temporal changes in RAD dissimilarity can illustrate changes in species evenness over time as might be expected if compositional homogenisation and a decline in evenness occurs because of invasion by an increasingly dominant species (Ruhl, 2008). Rank clocks can be used to visualise changes in composition over longer time series than Venn Diagrams (Collins et al., 2008; Hallett et al., 2016). Here, the rank order of abundance of each species is plotted at each time point around a central point to form a circle; no change in rank orders over time would be indicated by concentric circles (Collins et al., 2008). This method is available in the ‘codyn’ package in R (Hallett et al., 2016; Avolio et al., 2019). RADs can also be combined with measures of dissimilarities to understand which taxa are responsible for changes in composition; this and the comprehensive use of RADs for testing hypotheses of community change are available in Avolio et al. (2015).

Method 3: time cores and venn diagrams

Some methods are more descriptive than statistical, but are appropriate for testing for specific types of temporal patterns. For example, ‘time cores’ is a method that identifies ‘core’ groups of taxa that remain present over time (Vallès et al., 2014; D’Souza & Hebert, 2018). Others, such as Venn diagrams (Code S2) and measures of compositional overlap, are useful only for very small numbers of temporal comparisons. These methods are useful for assessing hypotheses of taxon replacement while specifically identifying the taxa driving compositional shifts (Code S2; Table 2).

Method 4: zeta diversity

Zeta diversity, a measure of similarity, analyses change in community composition among samples in the dataset (Codes S3 and S4). It is calculated by taking the average number of taxa that are shared between samples that are grouped in pairs, triples, quadruples, and so on (Table 2; Hui & McGeoch, 2014; McGeoch et al., 2019). Zeta diversity estimates compositional turnover from multiple (n ≥ 2) samples, rather than the pairwise comparisons used to estimate beta diversity. What is obtained by the method is a series of zeta values (similarity) for different orders of zeta: zeta₁ to zeta_n, where n is the total number of time points. Each of these values represents the average of the number of taxa shared by each set of times in the dataset. For example, zeta₂ (zeta order 2) represents the average number of taxa shared by all possible pairs of samples in the dataset (equivalent to 1—beta diversity or mean similarity), zeta₃ represents the average number of taxa shared by all possible subsets of three samples, and so on. Zeta₁ (zeta order 1) represents the average alpha diversity (taxon richness) of all the samples in the dataset.

If there is complete turnover in community composition in time, at some point, zeta_i (where i ≠ 1) will equal zero when there are no taxa shared by a given sized subset of samples. In contrast, if community composition never changes through time, zeta_n = zeta₁. Therefore, the rate of decline in zeta with increasing ‘zeta order’ tells us about the rate of change in community composition in time (‘zeta decline’). ‘Zeta decay’ can be considered analogous to time-lag graphs and time-lag regression analyses, which we described earlier for pairwise dissimilarity measures (Table 2). This means it experiences the same issue of few time points at higher lags, and thus higher uncertainty. The temporal scale(s) of community dynamics can be illustrated by plotting the change in sample zeta diversity of a particular zeta order against the temporal lag of samples, meaning that this method is useful for detecting both directional changes and cyclical patterns. For instance, if a multi-year time series in composition is cyclical (e.g. shows seasonality), zeta would seasonally decrease and then increase. These analyses can be done for different sets of samples to make spatial or temporal comparisons and, if necessary, zeta values can be normalised using the gamma diversity (total number of taxa in the dataset) to account for differences in taxon richness across samples (McGeoch et al., 2019). Further, different methods of creating the subsets for calculating zeta can be used, such as nearest neighbours or relative to a time-zero sample (McGeoch et al., 2019). This enables the user to assess taxonomic turnover at different orders of zeta and understand drivers of turnover in rare and common taxa (Latombe, Hui & McGeoch, 2017). This method could be extended to assess temporal turnover at different orders of zeta by adding time as an ‘environmental’ variable. Code and worked examples for the zeta diversity method is provided in the Codes S3 and S4 supplements and comes with its own R package, ‘zetadiv’ (Latombe et al., 2017).

Method 5: synchrony

Synchronicity refers to similarities and differences in the temporal dynamics among groups of taxa. A useful measure of this, called ‘synchrony’, was developed by Loreau & De Mazancourt (2008). This standardised value measures how closely changes in taxon abundances within a community track one another in time. It ranges from zero (complete asynchrony of taxon abundances in the community) to one (perfect synchrony; all taxon abundances changing in the same way) and, by exploring different sized ‘time windows’, can be used to detect scales of temporal variation. Overall significance tests of synchronicity in community change can be obtained by performing null model analysis (Pedersen et al., 2017). We provide code and worked example of this method in the Code S4 supplement, which comes with its own R package, ‘synchrony’ (Gouhier & Guichard, 2014).

A similar method (Houlahan et al., 2007) uses the sum of the temporal variances and pairwise covariances of all taxa in the community to distinguish if the community shows ‘compensatory dynamics’ over time. This method results in either a positive, negative or zero community-level covariance value, indicating that taxa are increasing or decreasing in synchrony (positive community covariance), some taxa are increasing, but this is compensated for by decreases in other taxa (negative community covariance), or that no change is evidence (zero community covariance) (Houlahan et al., 2007; Hallett et al., 2016).

Method 6: turnover rates

The degree of temporal variation in community composition can be assessed by calculating the turnover rate of the community (a.k.a. ‘temporal turnover’ or ‘species turnover’), which is a measure of the rate of change in taxonomic composition for a sample over time (Code S3; Aguirre et al., 2003). Turnover rates measure the magnitude of compositional change, but do not tell us about the direction of change or the taxa involved in compositional shifts. This method can detect when, and how rapidly, temporal shifts occur in a temporal sequence to detect complex temporal patterns such as perturbations, tipping points and pivot points (Table 2). Code for a worked example of this method can be found in Code S3.

Turnover rates are calculated using combinations of colonisation, immigration, extinction, mortality, recruitment, and survival. One of the commonly used methods is that of Diamond (1969) where the percentage turnover rate, T, is calculated as T = 100 × (E + C)/(S1 + S2), where E is the number of taxa that went extinct between two time points, C is the number of taxa that colonised between the same two time points, S1 is the number of taxa present at time 1 and S2 is the number of taxa present at time 2. T gives the percent change in taxon identities over the time period and ranges from 0, representing no turnover, to 1, indicating complete turnover. This measure can be calculated only for a single location at two time points; pairwise turnovers must be calculated separately for datasets with multiple time periods or multiple samples. T values can be compared through time, at different time intervals, or at different spatial scales (Boulinier et al., 2001; Tworek, 2004). They also can be used in a regression-style analysis to assess the effects of a treatment variable or other explanatory factor on compositional turnover (Aguirre et al., 2003; Tworek, 2004). Using T values in regression in this way is similar to the measure of ‘species turnover’ described above, based on pairwise dissimilarities (Baselga, 2010), because it represents the rate of appearance of new species over time, relative to the total number of species.

The ability to use turnover methods also is tied to the temporal extent of the study for the given system and the average lifespan of taxa. For example, complete turnover in prokaryotes can be observed over a period of days (Kara & Shade, 2009), whilst many decades may be required to detect any species turnover in some plant communities (Kardol et al., 2010). However, because it can be standardised by the number of taxa in the community, T is a potentially comparable measure of compositional change among different datasets, if they are on similar timescales and encompass communities with a similar distribution of life-histories.

Method 7: joint species distribution models

There are a range of regression approaches that simultaneously model individual species’ community dynamics, such as multispecies N mixture models and joint species distribution models (Warton et al., 2015). These methods can be used to model explicit hypotheses of community change, with a focus on interactions among taxa (Dorazio, Connor & Askins, 2015; Ovaskainen et al., 2017a, 2017b) and can work for datasets with reasonably high numbers of species (Thorson et al., 2016), but not yet for the very high numbers of taxa often resulting from molecular analyses. Some of these methods allow for useful additions, such as accounting for uncertainty in the detection of taxa, variation in traits and phylogenetic relatedness. In Hierarchical Modelling of Species Communities (HMSC), species-to-species association matrices are able to be estimated at multiple spatial or temporal scales. This enables a researcher to assess if temporal variation in community composition is due to intrinsic or external factors. Modelling is conducted in a Bayesian framework. We do not illustrate this method here because comprehensive R code and explanations are presented in the R package ‘HMSC’ (Tikhonov et al., 2020) and demonstrative examples are given in the literature cited above.

Case study 1—ants: a single sample measured multiple times (Code S1)

Record et al. (2018) measured ant community composition over time in experimentally treated, 90 × 90-m forest plots. We selected one plot, measured annually for 13 years (2003–2015), to demonstrate methods appropriate for detecting compositional change for one sample over time. The focal plot was part of the Harvard Forest Hemlock Removal Experiment (Ellison et al., 2010) whereby Tsuga canadensis (eastern hemlock) was removed by logging in 2005 after 2 years (2003–2004) of baseline data collection. An ordination trajectory plot (PCoA based on dissimilarity measure D10, Legendre & Legendre, 2012) showed there was little overall change in composition over the entire sample period (Fig. 1). Turnover analysis showed that compositional changes were driven initially by species losses, followed by species gains (Fig. 2) and the analysis of compositional overlap (Table 1) showed that at the end of the sampling period, the plot contained all the same species as in baseline sampling, along with an additional three species. Time lag analysis showed that the largest compositional differences appeared across approximately 4 years; that is, species composition changed the most at approximately 4-year intervals (Fig. 3).

Case study 2—plants: controlled experiment with treatments (Code S2)

Porensky et al. (2016a, 2016b) measured plant species composition annually from 2007 to 2014 in quadrats on 24 transects to which six grazing treatments were applied since 1982 (four transects in each treatment): none (exclosure), light, moderate, heavy, heavy-to-light (lightly grazed since 2007; ‘new light’) and heavy-to-none (new exclosure since 2007). They calculated the mean cover across 25 quadrats for each species on each of the 24 transects at each time (24 samples at each time).

The data exploration phase showed that the 24 transects were dominated across all measurements by ten common species out of a total of 92 species. Species turnover through time was high; mostly 30% or greater for time period-treatment combinations (Fig. 4). PERMANOVA showed that transects were significantly compositionally distinct and changed significantly over time in different ways (P < 0.05; Code S1: PERMANOVA outputs; based on dissimilarity measure D10, Legendre & Legendre, 2012). Partitioning spatial and temporal beta diversity showed that two of the control (exclosure) treatments differed significantly (Fig. 5). Principal Response Curves showed that the grazed treatments differed from the control and that the relative abundances of four species in particular were influential in driving these differences among treatments, but did not show increasing divergence over time (Fig. 6). This result is consistent with the analysis of temporal coherence, which showed that transects were not synchronous in their temporal changes (icc = −0.04). Thus, we conclude that treatment differences in turnover and shifts in composition were not consistent or directional over time relative to the control.

Case study 3—plants: multiple spatial and temporal replicates (Code S3)

The 50-ha forest census plot on Barro Colorado Island (BCI), Panama, has been censused eight times between 1982 and 2015 (Condit, 1998; Hubbell et al., 1999; Hubbell, Condit & Foster, 2005, http://ctfs.si.edu/webatlas/datasets/bci/, 2005). This plot was established in part to test for neutral community dynamics (Hubbell, 1979) by fully censusing all woody stems greater than 1 cm diameter at breast height every five years. The data exploration phase showed there were 324 species generating a highly right-skewed abundance distribution (most species were rare). A line graph of species’ abundance values over time showed very little change in total abundance over time. We found that change in species composition was lower than expected if individual trees were randomly distributed in time (keeping their spatial positions fixed): (1) Zeta diversity analysis showed that the presence of species across census times was more consistent than expected (Fig. 7) and (2) total turnover in species abundances at each of the seven time periods was lower than expected, particularly when the plot data were at larger spatial scales (based on dissimilarity measure D14, Legendre & Legendre, 2012). Thus, we conclude that the observed lack of pattern is consistent with neutral community dynamics; however, when researchers used basal area instead of total abundance in a trajectory analysis, some compositional change after drought was detected (De Cáceres et al., 2019).

Case study 4—bacteria: single site with seasonal sampling (Code S4)

Gilbert et al. (2012) measured bacterial taxonomic composition monthly over 6 years using pyrosequencing of the 16S rRNA gene for water samples from a single site off the coast of Plymouth, United Kingdom (source data can be obtained from the journal’s online Supplemental File). The data exploration phase showed that the dataset consisted of 74 operational taxonomic units (OTUs) across 62 temporal samples. Most taxa were common through time. Proportional turnover between temporal samples varied from 4% to 48%, which is relatively low, and was consistent with the detrended correspondence analysis (DCA) result showing a short gradient length of 1.8 for axis 1. Similarly, PCoA graphs showed little consistent change within years (Fig. 8). However, 13 of the OTUs were deemed to be highly variable in that they changed by more than 500 reads across the sampling period (mean number of reads = 296 ± 829 S.D.: row and column summaries, relative abundance distributions and abundance plots).

Time-lag analysis with Bray–Curtis dissimilarities (D14, Legendre & Legendre, 2012) showed that although there was a significant increase in dissimilarity within increasing time lag, variability was high, showing that temporal changes did not occur in a simple divergent or convergent manner (Fig. 9). Temporal synchrony analysis (sensu Loreau & De Mazancourt, 2008) showed that some subsets of OTUs in the community significantly changed in similar ways (Synchrony = 0.27; P < 0.05). When patterns in these variable taxa were examined through time using zeta diversity analysis (Hui & McGeoch, 2014), which is based on the presence of species in samples, winter samples shared significantly more species across sets of temporal samples than autumn, spring and summer (Fig. 10; ANOVA, P < 0.01). Further analysis of the most variable taxa using the raw dissimilarity values based on abundance data calculated for temporal samples within seasons showed clear cyclical patterns (Fig. 11). Thus, we conclude that the patterns are consistent with seasonal dynamics in community composition, largely driven by a small subset of OTUs that varied in their abundance, rather than presence, over time.