A new method for detecting and interpreting biodiversity and ecological community thresholds

Methods in Ecology and EvolutionVolume 1, Issue 1 p. 25-37 Free Access A new method for detecting and interpreting biodiversity and ecological community thresholds Matthew E. Baker, Corresponding Author Matthew E. Baker Department of Geography and Environmental Systems, University of Maryland-Baltimore County, Baltimore, MD 21205, USA Correspondence authors. E-mail: [email protected], [email protected] (Joint first authors)Search for more papers by this authorRyan S. King, Corresponding Author Ryan S. King Center for Reservoir and Aquatic Systems Research, Department of Biology, Baylor University, One Bear Place #97388, Waco, TX 76798, USA Correspondence authors. E-mail: [email protected], [email protected] (Joint first authors)Search for more papers by this author Matthew E. Baker, Corresponding Author Matthew E. Baker Department of Geography and Environmental Systems, University of Maryland-Baltimore County, Baltimore, MD 21205, USA Correspondence authors. E-mail: [email protected], [email protected] (Joint first authors)Search for more papers by this authorRyan S. King, Corresponding Author Ryan S. King Center for Reservoir and Aquatic Systems Research, Department of Biology, Baylor University, One Bear Place #97388, Waco, TX 76798, USA Correspondence authors. E-mail: [email protected], [email protected] (Joint first authors)Search for more papers by this author First published: 23 February 2010 https://doi.org/10.1111/j.2041-210X.2009.00007.xCitations: 322 Correspondence site: http://www.respond2articles.com/MEE/ AboutSectionsPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Summary 1. Existing methods for identifying ecological community thresholds are designed for univariate indicators or multivariate dimension-reduction of community structure. Most are insensitive to responses of individual taxa with low occurrence frequencies or highly variable abundances, properties of the vast majority of taxa in community data sets. We introduce Threshold Indicator Taxa ANalysis (TITAN) to detect changes in taxa distributions along an environmental gradient over space or time, and assess synchrony among taxa change points as evidence for community thresholds. 2. TITAN uses indicator species scores to integrate occurrence, abundance and directionality of taxa responses. It identifies the optimum value of a continuous variable, x, that partitions sample units while maximizing taxon-specific scores. Indicator z scores standardize original scores relative to the mean and SD of permuted samples along x, thereby emphasizing the relative magnitude of change and increasing the contributions of taxa with low occurrence frequencies but high sensitivity to the gradient. TITAN distinguishes negative (z−) and positive (z+) taxa responses and tracks cumulative responses of declining [sum(z−)] and increasing [sum(z+)] taxa in the community. Bootstrapping is used to estimate indicator reliability and purity as well as uncertainty around the location of individual taxa and community change points. 3. Using two simulated data sets, TITAN correctly identified taxon and community thresholds in more than 99% of 500 unique versions of each simulation. In contrast, multivariate change-point analysis did not distinguish directional taxa responses, resulting in much wider confidence intervals that in one instance failed to capture thresholds in 38% of the iterations. 4. Retrospective analysis of macroinvertebrate community response to a phosphorus gradient supported previous threshold estimates, although TITAN produced narrower confidence limits and revealed that several taxa declined at lower levels of phosphorus. Re-analysis of macroinvertebrate responses to an urbanization gradient illustrated disparate change points for declining (0·81–3·3% urban land) and increasing (6·8–26·6%) taxa, whereas the published threshold estimate (20–30%) missed the declining-taxa threshold because it could not distinguish their synchronous decline from the gradual increase in ubiquitous taxa. 5. Synthesis and applications. By deconstructing communities to assess synchrony of taxon-specific change points, TITAN provides a sensitive and precise alternative to existing methods for assessing community thresholds. TITAN has tremendous potential to inform conservation of rare or threatened species, develop species sensitivity models, identify reference conditions and to support development of numerical regulatory criteria. Introduction Ecologists have become increasingly interested in analytical methods for detecting and quantifying ecological thresholds (Brenden, Wang, & Su 2008; Andersen et al. 2009; Sonderegger et al. 2009). Ecological thresholds may be defined as transition points or zones of relatively rapid change between alternate ecosystem states or ecological condition, often in response to small, continuous changes in one or more causal variables (Toms & Lesperance 2003; Huggett 2005; Groffman et al. 2006). Ecological thresholds may be particularly relevant in the context of anthropogenic environmental gradients because such gradients can represent novel physical and chemical conditions falling outside those experienced by species in evolutionary time. Coevolved communities of interacting species possess unique morphological, behavioural and physiological adaptations, often corresponding to a narrow range of environmental conditions (sensuShelford 1913). Species distributions under otherwise unaltered conditions may simultaneously and abruptly change at a critical level of a novel environmental gradient (May 1977; King & Richardson 2003; Sasaki et al. 2008). Detection and description of such ecological community thresholds has important implications both for ecological theory and application (Townsend, Uhlmann, & Matthaei 2008; Martin & Kirkman 2009). Ecological community thresholds may be distinct from ecosystem-level, univariate thresholds described by current models (Groffman et al. 2006; Andersen et al. 2009; Suding & Hobbs 2009) because community responses are multivariate (i.e., one dimension for each taxon). Such thresholds are theoretically relevant to ecologists because of the evolutionary implications of a synchronous response of species to environmental pressures (Huggett 2005; Okland, Skarpaas, & Kausrud 2009). Potential applications of community thresholds include supporting the development of numerical environmental criteria to prevent loss of biodiversity and ecosystem function (King & Richardson 2003) or identification of biological reference conditions to better characterize community dynamics in the absence of disturbance (Wang, Robertson, & Garrison 2007; Brenden et al. 2008; Utz, Hildebrand, & Boward 2009). Current statistical methods used for identifying thresholds were not developed for multivariate species abundance data (Brenden et al. 2008; Andersen et al. 2009). The vast majority of taxa in community data sets have low occurrence frequencies (i.e., do not occur in a large proportion of the sample units) and highly variable abundances (McCune & Grace 2002). Consequently, most investigators aggregate community data into univariate responses, selecting a priori attributes that presumably represent an important facet of community structure, such as the number of taxa or deriving synthetic variables from multivariate analysis of taxa composition among sites (e.g., dissimilarity metrics, ordination axes; King and Richardson 2003, Walsh et al. 2005). While aggregating taxa into one or more response variables may, in some instances, increase community signal in response to anthropogenic gradients, it also likely obscures nonlinear changes in one or more taxa, potentially underestimating or misrepresenting the effect of an anthropogenic gradient on ecological communities. Thus, evaluating ecological community thresholds with existing approaches often involves undesirable generalities, loss of information or assumptions regarding taxon-specific responses. We introduce a new analytical approach, Threshold Indicator Taxa ANalysis (TITAN), with the goals of (i) exploring and identifying abrupt changes in both the occurrence frequency and relative abundance of individual taxa along an environmental, spatial or temporal gradient; (ii) quantifying uncertainty around locations of abrupt change; and (iii) estimating the relative synchrony and uncertainty of those changes as a nonparametric indicator of a community threshold. We contend that a flexible, taxon-specific approach can yield insight informing and augmenting understanding obtained from existing analytical methods. We illustrate accuracy and diagnostic advantages of TITAN using simulated community data and retrospective analysis of two aquatic macroinvertebrate data sets spanning different types of anthropogenic environmental gradients. Materials and methods Background and component calculations TITAN combines and extends change-point analysis and indicator species analysis. Change-point analysis (nCPA; King & Richardson 2003; Qian, King, & Richardson 2003) is a nonparametric technique that orders and partitions observations along an environmental gradient, identical to a single-split, multivariate regression tree analysis (MRT; De’Ath 2002). MRT replaces the univariate response of typical regression tree (RT; Breiman et al. 1984; De’Ath & Fabricius 2000) with a measure of multivariate dissimilarity among sample pairs. In both nCPA and MRT, optimal partitioning is obtained by maximizing a deviance reduction statistic that compares within-group vs. between-group dissimilarity determined by a user-selected ecological distance metric. As optimal partitioning can be sensitive to sample distribution along the environmental gradient, nCPA adds a bootstrap resampling procedure to assess uncertainty associated with the observed change-point value (King & Richardson 2003; King et al. 2005). TITAN replaces the aggregate, community-level, dissimilarity response of nCPA with taxon-specific, indicator value (IndVal) scores from indicator species analysis (Dufrêne & Legendre 1997). Indicator species analysis is a widely accepted method for identifying indicator taxa in noisy biological data, pruning dendrograms from a hierarchical cluster analysis to an optimal number of groups, or evaluating how sampling unit groupings explain species distributions. IndVal scores are a simple and intuitive measure designed specifically to capture strength-of-association between any particular taxon and an external sample grouping (e.g., an a priori cluster analysis typology). Indicator species analysis produces an IndVal score estimating the association of each taxon with each group. Thus, two IndVal scores are computed for a single taxon in a two-group classification. IndVal scores are the product of cross-group relative abundance (proportion of abundance among all sample units belonging to group i) and within-group occurrence frequency (proportion of sample units in group i with a positive abundance value). IndVal uses occurrence frequency within each group to weight a taxon’s relative abundance by how consistently it is observed in each group. A large abundance within one sample group results in a greater IndVal score only if the taxon also occurs with great regularity in that same group. IndVal scores are superior to simple abundance as a measure of association because integration of occurrence frequency and abundance leads to a measure of association that is unbiased by group size (Dufrêne & Legendre 1997). IndVal scores are scaled from 0% to 100% with a value of 100 indicating that a taxon was collected in every sample within a group and not in any other group. Conversely, in a two-group classification, a value of 50 can mean that a taxon occurred in half the samples within only one group, or in equal abundances in all samples of both groups (e.g., Fig. 2). The probability of obtaining an equal or larger IndVal score from random data (P) is estimated by comparing the magnitude of each observed IndVal score with those generated when group membership is assigned via ≥250 randomized permutations (Dufrêne & Legendre 1997). Figure 2Open in figure viewerPowerPoint Response of indicator value (IndVal) and z scores for six hypothetical taxa abundances along a uniform environmental gradient. Dashed vertical arrows illustrate the maximum z score and corresponding environmental change point. Change-point identification and uncertainty for a single taxon TITAN uses IndVal scores instead of deviance reduction (as in nCPA or MRT) to identify change points across a continuous environmental gradient (x; Fig. 1, Step 1). Midpoints between observed values of x are candidate change points (xi) used to iteratively split observations into two groups, and thus produce two IndVal scores at each split (Fig. 1, Step 2·1). The relative magnitude of IndVal scores for groups on each side of a candidate change point reflects whether a taxon shows greater association with the left (negative response with respect to x) or the right (positive response) side of each split (Fig. 1, Step 2·1C, D). The greater the difference in taxon fidelity (association) created by a specific split, the greater the IndVal score for one of the two groups. The greatest IndVal score at each split and the side of the split on which it occurs are retained for comparison with those at other candidate change points (Fig. 2). In practice, we use a minimum group size of five observations (De’Ath & Fabricius 2000), so for any sample of n observations and depending on the number of unique x-values observed, TITAN compares up to 2n−20 IndVal scores for each taxon (i.e., 2n for IndVals at each split, less the 10 samples and splits needed to satisfy the minimum group size). Any value of x resulting in an IndVal maximum among candidate splits is identified as the observed change point or the optimal partition for that taxon (Fig. 2). Figure 1Open in figure viewerPowerPoint Flow chart of Threshold Indicator Taxa ANalysis (TITAN). TITAN estimates uncertainty surrounding taxon-specific responses using the distribution of change-point values across a series of bootstrap replicates of the entire data set (Fig. 1, Step 3; Manly 1997; Toms & Lesperance 2003). The bootstrap procedure is necessary because unlike a priori group classification required by indicator species analysis, optimal group partitioning along x is initially unknown in TITAN, and is in fact the objective of the analysis. Whereas the permutation procedure is used to estimate the probability that an equal or larger IndVal could be obtained from random data, the bootstrap procedure estimates uncertainty around change-point locations (optimal partitioning along x), as well as consistency in the response direction of each taxon (negative or positive). Variability in change-point location, directionality (positive or negative with respect to x) and magnitude (relative to the absence of structure along x) constitute the information content of the indicator response for each taxon in TITAN. Two important diagnostic indices measuring the quality of the indicator response for any taxon are obtained from bootstrap resampling: purity and reliability. Indicator purity is the proportion of change-point response directions (positive or negative) among bootstrap replicates that agree with the observed response. Pure indicators (e.g., purity ≥0·95) are consistently assigned the same response direction, regardless of abundance and frequency distributions generated by resampling the original data. If bootstrap resampling substantially alters the probability of obtaining an equal or larger IndVal based on 250 random permutations of the data, then that particular taxon is not a reliable indicator. Indicator reliability is estimated by the proportion of bootstrap change points whose IndVal scores consistently result in P-values below one or more user-determined probability levels (e.g., P ≤ 0·05). Reliable indicators (e.g., ≥0·95 of the bootstrap replicates achieving P ≤ 0·05, or some other user-defined proportion of replicates) are those with repeatable and consistently large IndVal maxima. For each pure indicator taxon, TITAN uses bootstrap replicates to estimate empirical quantiles of the change-point distribution. Variation in change-point estimates highlights uncertainty in the location of the maximum IndVal with respect to x. Sharp, nonlinear responses in taxon abundances are reflected by relatively narrow intervals between upper and lower change-point quantiles (e.g., 5%, 95%), whereas taxa with linear or more gradual responses will have broad quantile intervals spanning most of the range of x. If the gradient is long enough to produce a Gaussian-like response, the interval will likely encompass the mode because of reliably strong, but impure bootstrap change points. However, we discourage strict-interpretation of these quantiles as confidence limits because any method of computing confidence limits will be unreliable for taxa with low occurrence frequencies (Manly 1997). Indentifying ecological community thresholds from multiple change points Once IndVals for each candidate change point and taxon are classified according to response direction, the aggregate response of all indicator taxa at each candidate change point may be used as evidence for a community-level threshold. As a dendrogram assessment tool, Dufrêne & Legendre (1997) recommended that the number of groups resulting in the largest sum of IndVal scores (scores significant at P < 0·05 or other criteria) be considered optimal for dendrogram pruning. This approach uses the greatest cumulative IndVal signal to distinguish groups, whether from a single taxon or many taxa, because a primary goal is to facilitate accurate classification of new observations using taxa most characteristic of each group. In contrast, evidence for a community threshold in TITAN requires substantial change across more than the most predominant taxa. As rare taxa are often highly sensitive to environmental alterations and the focus of biodiversity conservation, changes in their distribution are of great interest, although often more difficult to detect. Because the absolute value of IndVal scores is influenced by a taxon’s overall abundance, it is less important in TITAN than the magnitude of change relative to each taxon’s abundance distribution (see below). All taxa are not required to have identical IndVal maxima to produce a community threshold, rather only that large IndVal scores for many taxa occur close together along x. Integrating information about taxon-specific changes across groups, we apply the additive indicator score concept in a novel way, using rescaled indicator responses to partition observations while preserving the information provided by rare taxa. IndVal scores are rescaled as z scores within TITAN according to degree of departure from expected values by subtracting the mean of randomized permutations from the observed IndVal, and dividing by its permuted SD. Exploratory simulations suggest that permuted IndVal z scores are similar to nonparametric alternatives (i.e., substituting median and interquartile range for mean and standard deviation). We use z scores because the central focus of TITAN is somewhat distinct from the original purpose of IndVal. IndVal was developed to interpret a pre-existing sample typology (Dufrêne & Legendre 1997), whereas TITAN seeks to use IndVal scores to select among candidate groupings. Rather than raw IndVal magnitudes, which would favour the most widely distributed or abundant taxa, standardization facilitates cross-taxa comparison by emphasizing change in IndVals across candidate splits given a specific pattern of abundance and occurrence (Fig. 2). Rescaling makes little difference during determination of taxon-specific change points, but it can make a substantial difference during interpretation of their relative information content. Rare or infrequently occurring taxa with smaller IndVal magnitudes can have a very strong z score if their response to environmental change is dramatic. Standardized taxa responses increasing at the change point (z+) are distinguished from those decreasing (z−) and those showing no response. Evidence for community-level thresholds among negative and positive taxa is assessed separately by tabulating and summing all z− and z+ scores for each value of x. The value(s) of x resulting in the largest cumulative z scores for negative [sum(z−)] and positive [sum(z+)] responses correspond to the maximum aggregate change in the frequency and abundance of their respective taxa. Large values of sum(z) scores occur when many taxa have strong responses at a similar value of the environmental gradient, whereas weak or variable responses result in lower sum(z) values without a distinctive maximum. TITAN community-level change points may be assessed by plotting sum(z) scores vs. x, and are easily interpreted via tabular and graphical summaries of change-point distributions of individual taxa. If sum(z) maxima involve synchronous change in many taxa, including overlap from bootstrap distributions of taxon-specific change points, then these values of x may be interpreted as evidence for observed community thresholds. Bootstrap replicates used to evaluate taxon-specific change points are also summarized to develop distributions of sum(z) responses. Variation in the bootstrapped values of x that produce the greatest sum(z−) or sum(z+) values is used to estimate uncertainty associated with community change points, and quantiles (e.g., 0·05, 0·95) of these distributions serve as empirical confidence limits (e.g., Qian et al. 2003; Toms & Lesperance 2003). Narrow confidence limits represent further evidence for a community threshold, whereas wide confidence bands suggest other response (e.g., linear, modal or random) dynamics are more likely. Case studies Simulated community data We developed two simulations to evaluate how TITAN deals with taxa with distinct distributions, classifies responses, and detects change points. By controlling the statistical properties of the response and predictor variables across 500 unique data sets generated by each simulation, we demonstrate TITAN’s efficacy and flexibility. The first simulation involved a shift in respective taxa at two distinct values of x. Eight taxa abundances were simulated along a uniformly distributed environmental gradient (runif in R 2·9·2; n = 100, range = 0–100). Abundances were generated using a negative binomial distribution (rnbinom in R 2·9·2) to simulate noisy, heteroscedastic and sparse site-by-taxa matrices typical of community data (McCune & Grace 2002). Two taxa (Sp1–2) were thus assigned frequency and abundance values that differed below and above a value of 40 along the environmental gradient to simulate samples from populations with a threshold decline (Fig. 3). Likewise, three taxa (Sp6–8) were assigned increasing frequency and abundance values above 60 to assess TITAN’s ability to distinguish disparate change points and different response directions. One taxon (Sp3) was assigned an increasing change point at 40 and a decreasing change point at 60 to approximate a unimodal distribution. Finally, two taxa (Sp4–5) differed in frequency and abundance but varied randomly with respect to the environment. The entire simulation was repeated 500 times to generate different data sets, obtain diagnostic statistics, and evaluate the ability of nCPA (Bray-Curtis distance among sample units) and TITAN to correctly identify and assess thresholds in the data. Figure 3Open in figure viewerPowerPoint Threshold Indicator Taxa ANalysis (TITAN) and change-point analysis (nCPA, Bray-Curtis distance) of a two-threshold community response to a simulated environmental gradient (Simulation 1). (a) Simulated abundances of eight taxa in response to the environmental gradient (x-axis). According to negative binomial probability distributions used in simulating frequency and abundance, Sp1–2 should decline at 40 (black vertical line), Sp3 should increase at 40 (red vertical line) and decrease at 60 (black vertical line), Sp4–5 should not vary predictably with the environment, and Sp6−8 should increase at 60 (red vertical line). (b) Pure (≥0·95) indicator taxa are plotted in increasing order with respect to their observed environmental change point. Black symbols correspond to negative (z−) indicator taxa, whereas red corresponds to positive (z+) indicator taxa. Symbols are sized in proportion to z scores. Horizontal lines overlapping each symbol represent 5th and 95th percentiles among 500 bootstrap replicates. Vertical lines indicate the simulated true values for negative (black vertical line) and positive (red vertical line) underlying thresholds. (c) TITAN sum(z−) and sum(z+) values corresponding to all candidate change points (xi) along the environmental gradient. Black and red vertical lines represent the cumulative frequency distribution of change points (xcp, or thresholds) among 500 bootstrap replicates for sum(z−) and sum(z+), respectively. (d) Deviance reduction in Bray-Curtis distance values for each candidate change point (xi) along the environmental gradient. The dashed blue line represents the cumulative frequency distribution of change points (thresholds) among 500 bootstrap replicates. The second scenario involved similar distributions of the environmental gradient to contrast threshold responses with noisy data and generalized, wedge-shaped distributions typical of complex taxon responses to multiple limiting factors (e.g., Cade & Noon 2003; Brenden et al. 2008; Konrad et al. 2008). In this scenario, two taxa (Sp1–2) were assigned decreasing frequency and abundance values to simulate samples from populations with a threshold decline at 20 (Fig. 4). Three taxa (Sp6–8) were assigned probabilities of frequency and abundance that increased in proportion to x (i.e., generating wedge-shaped distributions). Sp3 was assigned an increasing change point at 20 and a decreasing change point at 60 to approximate a broad, unimodal distribution, whereas Sp4–5 varied randomly with respect to the environment. This simulation was also repeated 500 times to obtain diagnostic statistics. Figure 4Open in figure viewerPowerPoint Threshold Indicator Taxa ANalysis (TITAN) and change-point analysis (nCPA, Bray-Curtis distance) of a negative threshold and positive wedge-shaped community response to a simulated environmental gradient (Simulation 2). (a) Sp1–2 should decline at 20 (black vertical line), Sp3 should increase at 20 (red vertical line) and decrease at 60 (black vertical line), Sp4–5 should not vary predictably with the environment, and Sp6–8 should increase with a wedge-shaped distribution (red dotted lines) instead of a threshold response. See Fig. 3 for additional details. Everglades data These data were taken from a previous study designed to identify a concentration of surface-water total phosphorus (TP) that corresponded to abrupt changes in macroinvertebrate species composition in the Florida Everglades, USA (King & Richardson 2003). Macroinvertebrate densities (no/m2, 164 taxa, species or morphospecies-level taxonomy) were measured from 126 marsh sampling stations along a 10-km TP gradient. Concentrations of TP in the data set ranged from <10 μg/L to >100 μg/L. The authors used several community variables and estimated TP change points using univariate nCPA. The resulting change points ranged from c. 10 μg/L to 25 μg/L TP, and authors concluded that TP >12–15 μg/L likely corresponded to ecologically significant changes in taxonomic composition. Maryland stream data These data were the subject of a previous study on analytical considerations for linking watershed land cover to ecological communities in streams (King et al. 2005). In the previous analysis, we used axis scores from non-metric multidimensional scaling (nMDS) of Bray-Curtis dissimilarity in nCPA to identify a level of watershed percent developed land corresponding to an abrupt change in macroinvertebrate community composition in wadeable streams (295 sites, 177 taxa abundances, mostly genus-level identification). Our previous analysis identified a relatively sharp change in community composition (inferred from nMDS axis 1 scores) between 20% and 30% watershed developed land (5th–95th bootstrap percentiles). Data analyses We performed TITAN analysis on all four data sets in R (R Development Core Team 2009, version 2·9·2) using a custom package TITAN (see Appendix S3) written by MEB and RSK. nCPA was performed using a custom function within TITAN based on the db-MRT method of De’Ath (2002) in the package mvpart. We log10(x + 1) transformed taxa abundances to reduce the influence of highly variable taxa on indicator score calculations in each data set, which was particularly important for taxa with low occurrence frequencies. Taxa with <5 occurrences were deleted (following previous analyses of these data) and we used Bray-Curtis distance as the dissimilarity metric for all nCPA assessments (King et al. 2005). We compared TITAN and nCPA by plotting sum(z) and deviance reduction values as a function of increasing values of x and identified community change points as the x value that resulted in the maximum sum(z) or deviance reduction, respectively. We computed cumulative change-point distributions by finding the maximum IndVal (individual t