Increasing instability of a rocky intertidal meta-ecosystem.

Climate change threatens to destabilize ecological communities, potentially moving them from persistently occupied "basins of attraction" to different states. Increasing variation in key ecological processes can signal impending state shifts in ecosystems. In a rocky intertidal meta-ecosystem consisting of three distinct regions spread across 260 km of the Oregon coast, we show that annually cleared sites are characterized by communities that exhibit signs of increasing destabilization (loss of resilience) over the past decade despite persistent community states. In all cases, recovery rates slowed and became more variable over time. The conditions underlying these shifts appear to be external to the system, with thermal disruptions (e.g., marine heat waves, El Niño-Southern Oscillation) and shifts in ocean currents (e.g., upwelling) being the likely proximate drivers. Although this iconic ecosystem has long appeared resistant to stress, the evidence suggests that subtle destabilization has occurred over at least the last decade.

Overpopulation and industrialization threaten human society and global ecosystems alike via the existential challenge of climate change (1). Ecologists are among those at the forefront in understanding and meeting this challenge. Crucial tools include long-term datasets, field experimentation, and modeling, enabling deep insight into patterns of change, potential causal links, and possible future outcomes. In the marine environment, novel occurrences such as marine heatwaves (MHWs) (2) and disease epidemics (3) add new, more acute, and more episodic threats to the more gradual increases in warming and ocean acidification typifying climate change. Increasing evidence indicates that marine ecosystems are already responding to climate change and MHWs (1, 4–6).

Ecological change in response to a changing environment can vary from gradual, monotonic responses to sudden shifts that result in radically different states (7, 8). Theory and evidence have focused on the latter concept, i.e., alternative stable states, in which a system can exist in two or more configurations (9–11). For example, extreme wildfires show how gradual changes in temperature or rainfall can lead to sudden catastrophic and often irreversible ecosystem changes (12, 13). Despite such examples, however, a recent summary of meta-analyses suggests that overall empirical evidence for thresholds is scarce (14), possibly because sudden changes are hard to detect. This situation seems likely to change as climate change accelerates.

While predicting sudden changes is difficult, systems nearing a tipping point may exhibit “early warning signals” that potentially enable anticipation of these shifts (e.g., refs. 15 and 16). One postulated sign of a system nearing a tipping point is increasing variability in community structure (16–18). Another is “critical slowing down,” in which a system recovers increasingly slowly from small perturbations (16). These ideas are conceptualized in Fig. 1, portraying possible trajectories of communities that are resistant, resilient, or diminishingly resilient to perturbations (Fig. 1, black, red, and blue lines, respectively). In response to perturbations, resistant or resilient systems can maintain or recover to their original configurations (Fig. 1, red and black lines, respectively). Systems with diminishing resilience may show signs of instability and eventually reconfigure to a future new state (Fig. 1, blue lines).

Fig. 1.

Conceptual models showing community responses to perturbations (e.g., annual disturbances). Variation of community states of resistant (black lines), resilient (red lines), or diminishingly resilient (blue lines) communities. (A) The resistant community (black line) is slightly affected by the disturbance and varies in structure around a long-term “original” state, while two other communities respond to disturbance by shifting to a different state. One (red line) is resilient, recovering within some time period (in this example, a year), while the other (blue line) exhibits diminished resilience, with increasingly slow recovery. Such slowing recovery rates might lead to a shift to a different future state. These dynamics can be differentiated by measuring trends in (B) mean community state and (C) variability of community state or community processes over time.

To test the resilience of low intertidal communities across multiple sites spanning a range of environmental conditions, we conducted annual disturbance–recovery experiments. At the continental scale, these communities vary in state in relation to gradients in air and sea temperatures, upwelling, precipitation, shelf width, and tidal range along the US West Coast (19–22). Multiple ecological processes in these communities are responsive to changing environmental conditions. Examples include predation rates by sea stars on mussels increased with rising temperature (23). Prey recruitment rates fluctuated with the North Pacific Gyre Oscillation (NPGO) and upwelling (24, 25). Growth and elemental composition of macrophytes were sensitive to fluctuations in upwelling, El Niño, and the Pacific Decadal Oscillation (PDO) (26).

Prior research suggested that sensitivities of community pattern may be decoupled from process such that apparent relative stasis in structure obscures variable dynamics (27, 28). In the low intertidal zone, despite climate-sensitive process rates (e.g., predation, competition, recruitment, and growth; see citations in previous paragraph), community structure has been remarkably stable over at least 14 y () and, based on observations by the lead author, back to the early 1980s. This persistence occurred despite interannual fluctuations in environmental conditions, including climate (PDO, NPGO, El Niño), upwelling, temperature, sea star wasting, and orders-of-magnitude variation in prey recruitment [ref. 29; ].

Our experiments tested whether the apparent stasis of community structure hid underlying changes in dynamical processes, using repeated annual disturbances to test resilience. Studies were conducted at six sites paired within each of three regions or capes ranging across 260 km of the Oregon coast from 2011 to 19 in the low intertidal zone (). Disturbed plots (n = 5 per site) were cleared each July and were paired with intact uncleared plots (all 0.5 × 0.5 m). Plots were photographed monthly to seasonally, and the annual percent cover of all taxa in each treatment was estimated from these photos. These experiments were used to ask whether the kinetics underlying the apparent stasis of community structure were also temporally unchanging (please see and for more details).

Results

Spatial Variation in Community Structure.

As documented previously (19), communities within capes were similar, but cape-scale differences were large, with the central cape differing greatly from northern and southern capes (Fig. 2; ). Although communities in intact and disturbed treatments varied temporally, each cape remained in distinct areas of 2D ordination space. Intact and disturbed communities at northern (sites Fogarty Creek [FC], Boiler Bay [BB]) and southern (sites Cape Blanco [CB], Rocky Point [RP]) capes were dominated by macrophytes (algae and surfgrass), while those in the central cape (sites Yachats Beach [YB], Strawberry Hill [SH]) were dominated by sessile invertebrates (mussels, barnacles, anemones) and bare space. Intact communities in the northern cape had higher kelp and surfgrass (Phyllospadix spp.) abundances, while those at the southern cape had more red algae (e.g., Mazzaella spp., Polysiphonia spp., Microcladia spp., Odonthalia spp.) (Fig. 2; see for functional group compositions).

Fig. 2.

nMDS analysis of the recovery experiment. Time series trajectories by year (2012–2019) and sites (symbols) in intact control plots (solid lines) and clearance recovery plots (dashed lines). Site codes are listed in . Cape Foulweather sites are shown in light and dark blue shading, Cape Perpetua in red and purple shading, and Cape Blanco region sites in light and dark green shading. Vector overlays depict taxon functional groups driving community separation.

Recovery Experiments.

Ordinations distinguished results on both taxon-composition and treatment scales. Thus, nonmetric multidimensional scaling (nMDS) 1 (x axis) separated algal-dominated (Fig. 2, Left) and invertebrate-dominated (Fig. 2, Right) communities, while except for RP, nMDS2 (y axis) reflected treatments (controls higher, removals lower). In most cases, disturbance moved communities toward structures dominated by bare space and “weedier” taxa such as barnacles and filamentous algae (Fig. 2; removal communities moved downward and rightward in nMDS space).

Temporal Trends.

In disturbed communities, overall average recovery rates declined (slope = −0.499, P = 0.0015, variance explained = 18%) and overall mean variability increased (slope = 0.505, P < 0.0001, variance explained = 34%) (Fig. 3 ). Site-scale regressions indicate that consistent with the overall trend, mean community structure decreased at most sites, and especially at southern sites. The variability (SD) of community structure increased at all sites, most strongly at the northern sites ().

Fig. 3.

Rates of change and variability of disturbed and intact experimental plots. Data are final community structure after 1 y of recovery (see ) in treatments across all six sites from 2011 to 2019. (A) Annual community structure in disturbed plots. (B) Annual variation (SD from the mean) in recovered community structure in disturbed plots. (C) Annual community structure of intact communities. (D) Annual variation (SD from the mean) in intact community structure. Linear regression equations: (A) Mean Community Structure = 1,024 − 0.05 × Year, F = 0.0015, n = 48; (B) Community Variability = −991.5 + 0.505 × Year, F < 0.0001, n = 48; (C) Mean Community Structure = 139.3 − 0.06 × Year, F = 0.29, n = 48; (D) Community Variability = −750.8 + 0.38 × Year, F = 0.021, n = 48.

In intact communities, mean community structure was persistent and stable throughout the study period (Fig. 3, P = 0.29, variance explained = 0.3%). However, as in recovery plots, variability of intact community structure increased (Fig. 3, P = 0.02, variance explained = 9.1%). Site-scale regressions on variation in community structure were generally weakly positive ().

Analysis of temporal differences (vector lengths) between intact and recovery plots yielded a similar result (Fig. 4; see for an example visualization). Mean vector length and vector length variability (SD of replicate means) increased through time, indicating that recovery rates were slowing and becoming more variable.

Fig. 4.

Analysis of vector-based rate of community recovery. Vector lengths are the Euclidean distance between intact and disturbed plot final community structures based on x and y coordinates obtained in ordinations. (A) Mean vector length (linear regression: length = −55.83 + 0.028 × year; P = 0.039, adj. R2 = 0.08, n = 42) and (B) SD of mean vector lengths (linear regression: SD = −86.37 + 0.043 × year, P = 0.0009, adj. R2 = 0.223, n = 42). See for an example of between-treatment vector lengths through time.

Detailed analysis revealed that communities differed most strongly by region (cape) but also varied with year, site within cape, and treatment (permutational analysis of variance [PERMANOVA]; ). Pairwise comparisons suggested that divergence of disturbed communities from the original state started in ∼2015 and generally persisted (). This trend occurred across all sites within capes but was strongest at the northern and southern capes (). Ordination was consistent with these results; in each cape, early-2010s (2011 to ∼2014) community states shifted to a different late-2010s (2015–2019) configuration ().

Communities on each cape varied similarly through time, suggesting responsiveness to a common driver or set of drivers. For example, most significant differences in pairwise tests between years occurred in comparisons between early years to years after 2015–2016 (, Lower left 16 cells). Similarly, differences between years after ∼2015 were minimal (, lower right six cells). With the exception of Cape Perpetua, trajectory divergence was weaker in intact than in disturbance plots [; 8 versus 29 lower left cells (of 64 total) differed from earlier years]. Inspection of functional group-scale versus year time series suggests these patterns were driven by changes in abundance of several taxa, both invertebrates and algae, in the late 2010s (). For example, at CP a massive gooseneck barnacle recruitment event occurred in 2014, leading to exceptional abundances in intact plots at YB and SH and in disturbed plots at SH in 2015, declining relatively slowly thereafter (). Similar changes in the 2015–2016 timeframe occurred for other taxa but differed more clearly between treatments and among sites ().

What drove these changes? Community structure varies with population changes, species interactions, disturbance, environmental stress, and propagule supply. Prior research revealed among-cape (spatial) differences in environmental factors such as upwelling, shelf width, temperature, nutrients, pH, and phytoplankton [(19, 30); ]. Furthermore, community processes like recruitment, predation rates, and growth also responded to longer-term climatic fluctuations such as ENSO (El Niño-Southern Oscillation), NPGO, and PDO (29). However, these factors do not necessarily vary in synchrony (). Analysis of environmental regime [upwelling (cross-shelf upwelling transport index) (31), climate oscillations, intertidal sea surface temperature (ISST), and years since sea star wasting (YSW), a disease perturbation] revealed cape-scale differences, explaining 19% of the variation (PERMANOVA; ). PCO (principal coordinates ordination) shows that environmental regimes shifted strikingly during the experiment, with large changes occurring beginning in 2014 to even more extreme changes in 2015–2016, thereafter cycling in the new region of 2D space (). Although environmental trends followed similar trajectories at each cape, those at Cape Blanco sites (green lines, ) were distinct from those at the other capes. As reflected by ISST, ENSO, PDO, and possibly YSW, the temporal shift was related to temperature (, axis PCO1 explaining 54.8% of total variation), while the north–south cape-scale difference was driven by stronger upwelling southward (, axis PCO2 explaining 22.1% of the variance).

We assessed the level of association between experimental community structure and the environmental regime using distance-based linear models (DistLM) and distance-based redundancy analysis (). As expected from among-cape environmental differences, the southern cape results for intact and disturbed plots were distinct from those for the central and northern cape. The interannual trajectory analysis showed that while environmental trajectories for each site within each cape followed a similar pattern, the northern and central site trajectories were strongly overlapping and occupied 2D space that was quite different from that occupied by southern cape sites.

Environmental regime effects on temporal trends in community structure differed between treatments (). Intact plot trajectory shifts reflected variation in temperature and temperature-related metrics (ENSO, NPGO, PDO, YSW, and ISST) and upwelling (). Specifically, using the proportion of variance explained column in DistLM sequential tests of , 45% (0.135/0.298) of cumulative variance was explained by thermal variation, and 54% (0.163/0.298) was explained by upwelling. Disturbed plots, in contrast, first responded mostly to thermal changes (2012/2013–2014/2015), then underwent changes more related to upwelling (). Temperature-related metrics explained 32% (0.072/0.223) of cumulative variation and upwelling explained 68% (0.151/0.223).

Discussion

Our experiments suggest that despite relative stasis of community structure over past decades, these low intertidal communities have been destabilizing since at least ∼2011. Annual recovery rates have slowed, and recovery variability has increased. These shifts were associated with variation in environmental regime, particularly ocean warming (; 2015–2016 ENSO and the 2014–2016 MHW). The onset of these changes also coincided with an unprecedented wasting event (3) that killed most individuals of the sea star Pisaster ochraceus. Based on classic experimental results (32), we expected this shock to the system to trigger invasion of the low shore by mussels and thus a coast-wide shift to mussel domination. Although small increases in mussels occurred in intact plots at the four northern sites (), this shift did not occur. Instead, we detected more subtle changes in composition and abundance of several algal and invertebrate taxa, but overall states of each cape persisted (). That is, invertebrate-dominated communities remained invertebrate-dominated and macrophyte-dominated communities remained macrophyte-dominated (e.g., ref. 29). Consistent with this, YSW accounted for only 2.5% (intact) and 1.6% (disturbed) of cumulative variance in community structure (; proportion of variance explained column in sequential tests). On the other hand, temperature-related variance accounted for 45% (intact; see previous paragraph) and 32% (disturbed) of community variation ().

Our disturbance treatments mimicked complete loss of all algae, as could occur if environmental stress were severe (e.g., in an MHW). Another form of disturbance in this system is removal by waves during storms. In the northeast Pacific, significant wave heights have increased (33), suggesting that losses from wave-caused disturbances may increase. Two factors suggest that such effects may not lead to the magnitude of losses expected from warming stress. First, excepting surfgrass, most intertidal macrophytes die back in fall and winter, overwintering as tiny alternate life stages (e.g., gametophytes) or as holdfasts. Second, as a consequence, since storms are most frequent in winter months, losses of low intertidal macrophytes from wave action are typically minimal and mostly limited to surfgrass (27). Among sessile invertebrates, barnacles and anemones are largely immune to direct wave-induced loss because of their strong attachments. Mussels, well-known to be susceptible to wave disturbance (34), are sparse in the low zone (except at the central cape). Hence, wave-driven disturbance is unlikely to be an important driver of low intertidal stability.

Basins of Attraction.

A useful stability framework is the ball-and-cup paradigm, in which balls represent local communities and cups represent basins of attraction across an environmentally driven seascape (e.g., refs. 7, 35, and 36). In our system, the seascape is defined by large-scale environmental drivers such as varying upwelling and thermal regimes generating multiple basins of attraction (here, regions or capes). Each basin can include multiple local communities or ecosystems connected by flows of energy and materials [i.e., meta-ecosystems (19, 37, 38)] within and between basins.

Our ordinations suggest that across the three capes, there were two distinct community states (Fig. 2) and two distinct physical environmental states (). Specifically, the northern and southern cape communities lie on the algal-dominated side of 2D space, and the central communities lie on the invertebrate-dominated side (Fig. 2). However, contrasting patterns in the physical environments were revealed by the environmental ordination (). That is, as driven primarily by upwelling differences, physical environments at northern and central capes were similar but distinct from that at the southern cape. With one key exception, shelf width (19), other metrics of physical conditions are also similar between the northern and central capes (39), including air temperature; wind velocity and direction; wave period, height, and power; and salinity. Thus, in this context, the northern (macrophyte-dominated) and central (sessile invertebrate-dominated) cape ecosystems might be considered alternative states since they have persisted for decades in different structural configurations despite similar physical environments. On the other hand, by altering the cross-shelf nearshore currents, the wider shelf off the central cape (19) drives two key large-scale biological differences that underlie the community differences. Specifically, retentive currents along the central cape (40) hold phytoplankton and invertebrate larvae relatively close to shore, while strong cross-shelf currents along the northern cape move plankton offshore. Through a variety of direct and indirect effects, prior research indicates that these biological consequences of bathymetric variation are the drivers of community differences between northern and central capes (19, 24–27, 29, 41).

At the larger scale of the Oregon coast, our experiments and observations suggest that although these low intertidal communities have persisted in similar cape-specific states for years, thus demonstrating resistance and resilience to major perturbations such as the 1982–1983 and 1997–1998 El Niño events, system resilience may be diminishing (e.g., Fig. 1). Across all sites and capes, rates of recovery of community structure have slowed and the variability in community structure has increased (Figs. 3 and 4), consistent with the ideas of “critical slowing down” and “increasing variance” as a system destabilizes. However, predicting the future trajectories of these communities is presently impossible, and in fact, the trends shown here may reverse. For example, site-scale regressions of mean community state and SD versus year in disturbed and intact communities show a few shifts from negative trends to positive trends through the 2010s (; quadratic regressions). Specifically, at FC mean community state declined, then increased in recovery plots, and at SH variability in intact communities increased, then decreased from 2011 to 2019.

Nonetheless, coastal environments clearly are worsening. For example, ENSO and MHW events are predicted to increase in frequency, if not magnitude (2, 42). However, some indication of a possible future might be inferred from recent collapses of subtidal macrophyte-dominated (kelp) systems (43–44). In macrophyte-dominated low intertidal zones, kelps such as Hedophyllum sessile, Lessoniopsis littoralis, Egregia menziesii, and Laminaria spp. can be dominant canopy-forming foundation species (45, 46). Many other species of algae and invertebrates are harbored beneath these species, so kelp loss can trigger a cascade of additional losses as well as create habitat for new species (45, 47). A recent study reports that these intertidal kelps are sensitive to the warmer conditions created by ENSO and MHW, with populations suffering huge losses during the 1997–1998 El Niño and suppressed growth during the warm waters of 2014–2016 (46). We infer that continued climate warming may cause intertidal macrophyte populations to decline, leading to drastic shifts in community state. Impacts on sessile invertebrate–dominated low intertidal communities are less clear, but our evidence suggests that mussel performance is also declining, while as yet barnacles seem unaffected by warming.

Mechanisms.

How do warmer conditions affect taxon performance? The literature indicates that direct thermal stress affects algal performance and survival of, and colonization by, sessile animals, while indirect negative effects of ENSO- and MHW-driven nutrient stress likely affect algal growth and phytoplankton abundance. For example, culturing experiments with kelp indicate that juvenile life history stages fail to survive at temperatures ≥18 °C (48). A recent study suggested inhibition of kelp growth by warmer water temperatures (49), and another indicated that surfgrass growth is suppressed by high (desiccating) air temperatures (50). Many intertidal invertebrates are known to be near their upper thermal tolerance limits, even in the temperate habitats along the Pacific coast of North America (51, 52), which are often exposed to air temperatures ≥25 °C (e.g., ref. 53). During our years of observation, periodic mass mortalities of adult mussels and barnacle recruits have occurred, as have near-annual algal and surfgrass bleaching (desiccation) events. Other stressors, such as changes in carbonate chemistry, low oxygen levels, and variable upwelling, can also negatively affect marine biota in this region (19, 30, 54). Collectively, through trophic linkages, these direct and indirect effects are likely to underlie more subtle reductions in performance of the sessile biota, thus affecting recovery rates by modulating growth rates and interaction strengths.

Stability.

The issue of stability in ecological systems has a long and contentious history, with foci on the stability-diversity relationship (e.g., Ref. 55) and on the meanings of components of stability [e.g., resilience, resistance, persistence (9, 56, 57)]. Connell and Sousa (56) assessed evidence needed to document system stability and argued that structures maintained as alternative states through press manipulations were not true alternative states. They noted that by definition, true alternative states must occur without intervention, under similar physical environmental conditions, and that they must persist longer than the longest-lived species in the system. Others (58) suggested that persistence through multiple generations is “less strict” and “more biologically practical.” The Connell-Sousa (56) stipulations thus preempted many early studies cited as showing alternative states, since the required long-term studies in ecology then were exceedingly rare. Since that time, examples apparently meeting the Connell-Sousa (56) long-term and nonintervention criteria have appeared (reviews in refs. 13, 59, and 60), indicating that alternative states have empirical reality (but see Ref. 14). As noted previously, one interpretation of our results is that the northern and central capes are alternative states. The community differences between sites on these capes appear to meet one of the Connell-Sousa criteria: they have persisted for decades, which doubtless exceeds the length of life of most if not all intertidal taxa. Furthermore, the physical environments are very similar. However, biological environments differ, and these differences explain the contrasting community states (19, 24–27, 29, 41), so whether or not the northern and central communities are truly alternative states is controversial.

We conclude that although rocky intertidal systems are highly dynamic, this meta-ecosystem in recent years has begun to show signs of slowing recovery and increasing variability, presumably in response to unprecedented stresses related to acute warming events. Although the ecosystems have not shifted to new states, our experiments suggest that they are exhibiting diminished resilience to disturbance, indicative of “critical slowing down.” Even intact communities are becoming more variable, which may be a harbinger of instability and be an early warning signal of community state change. We plan to assess whether they do indeed shift (or hopefully not) through continued research.

Materials and Methods

Estimating Abundances.

We initiated the experiment series in 2011 by haphazardly locating, then permanently marking five pairs of 0.5 × 0.5 m low intertidal plots. In July 2011, we photographed, then cleared one plot of each pair of all macrobiota, including macrophytes and sessile invertebrates. Adjacent uncleared reference plots were left intact. Clearances (“disturbed plots”) were allowed to recover without further intrusion for 12 mo, when they were photographed and recleared. Experiments were conducted at each of two sites (separated by 0.5–10 km) nested within three regions or capes (separated by 61–197 km) along ∼260 km of the central and southern Oregon coast ().

Percent cover of each species was estimated by inspecting annual plot photographs taken after 1 y of recovery, prior to reclearing. The 0.5 × 0.5 m quadrats were subdivided into 0.1 × 0.1 subquadrats, each consisting of 4% cover. Abundances of each taxon were estimated by eye for each subquadrat, and totals were obtained by adding across all 25 subquadrats. In the ordinations, we grouped species into functional groups (see for group composition).

We used two methods to estimate temporal patterns of recovery rates. First, after estimating the final abundance of each taxon (i.e., after 12 mo) in each replicate of each site × year × treatment combination, data were square-root transformed, then standardized to the most abundant taxon in each replicate. Means and SDs were calculated across taxa for the five replicates for each site, year, and treatment. Second, using x and y coordinates obtained from ordinations, we calculated annual final vector length (Euclidean distance) between intact (control) and recovery (removal) plots by year and site (see for an example visualization). We reasoned that if the rate of recovery in cleared plots was slowing, the ordinal distance between community states in control versus removal plots should increase through time. We averaged across replicates in each site × year × treatment combination to obtain mean and SD of mean vector lengths through time.

Because communities and environments are multivariate, we used ordination approaches to quantify community and environment state. Specifically, by definition, communities consist of multiple response variables such as species, taxon, or functional group abundances and environments consist of multiple metrics such as temperature, salinity, humidity, pH, O2 and CO2 levels, and precipitation. Climate oscillations such as upwelling, ENSO, NPGO, and PDO represent integrated environmental measures driven by variation in multiple variables such as temperatures, currents, and winds, which generate environmental variability. For data analyses, we used JMP v14.0.1 for univariate analyses (primarily regressions, means, and variance measures) and PRIMER v7 with PERMANOVA+ for ordination and multivariate analyses. Community data were square-root transformed and standardized to the maximum value for each variable prior to analysis. We quantified community structure using Bray-Curtis similarity matrices. We analyzed variation in community structure using PERMANOVA, with year, cape, and treatment as fixed factors and site nested within cape as a random factor. Fixed effects sum to zero for mixed terms. The permutation method was permutation of residuals under a reduced model with 9,999 permutations.

We used DistLM to analyze the association between environmental variables and community structure by treatment and distance-based redundancy analysis to visualize the results. Environmental indices (already anomalies) and water temperature were not transformed. We included the YSW disease event as an environmental factor because it represented the occurrence of an external coast-wide disease epidemic. To visualize among-site variation in environment over time, we plotted ordinations of temporal changes in intact and disturbed plots and in the environmental regime by site using nMDS and PCO.