Effects of Dam Removal on Fish Community Interactions and Stability in the Eightmile River System, Connecticut, USA

Abstract

New multivariate time-series methods have the potential to provide important insights into the effects of ecosystem restoration activities. To this end, we examined the temporal effects of dam removal on fish community interactions using multivariate autoregressive models to understand changes in fish community structure in the Eightmile River System, Connecticut, USA. We sampled fish for 6 years during the growing season; 1 year prior to, 2 years during, and for 3 years after a small dam removal event. The multivariate autoregressive analysis revealed that the site above the dam was the most reactive and least resilient sample site, followed in order by the below-dam and nearby reference site. Even 3 years after the dam removal event, the stream was still in a recovery stage that had failed to approximate the community structure of the reference site. This suggests that the reorganization of fish communities following dam removals, with the goal of ecological restoration, may take decades to centuries for the restored sites to approximate the community structure of nearby undisturbed sites. Results from this study also highlight the utility of multivariate autoregressive modeling for examining temporal interactions among species in response to adaptive management activities both in aquatic systems and elsewhere.

Appendices

Appendix 1

Mean (+SE) environmental characteristics for each sample site-year combination during the Zemko Dam removal project in the Eightmile River in Connecticut. Asterisks (*) indicate significant differences in environmental characteristics over the time series for each sample site location according to a mixed model analysis. Fig. 3

Fig. 3

Mean (+ S.E.) environmental characteristics for each sample site-year combination during the Zemko Dam removal project in the Eightmile River in Connecticut. Asterisks (*) indicate significant differences in environmental characteristics over the time series for each sample site location according to a mixed model analysis

Appendix 2

Detailed description of the MAR modeling framework

MAR models are stochastic and they examine temporal shifts in community stability, species interactions, disturbances, and environmental characteristics, while accounting for temporal autocorrelation in species abundances. Multivariate state–space models (MARSS) were also considered prior to choosing MAR because of their ability to account for process error and enhanced variance partitioning (Hampton et al. 2013). MARSS model interaction strengths and off-diagonal terms of the B matrices were similar to the MAR model results, and the MARSS models did not result in better model fits (i.e., lower AIC/BIC) than the MAR models. Therefore, we felt confident in using the MAR analytical framework in this study. Hampton et al. (2013) define MAR as a system of p linear equations describing the abundances for each species in the community, which is written as:

$${{\bf{X}}_t} = {\bf{B}}{{\bf{X}}_{t - 1}} + {\bf{A}} + {\bf{C}}{{\bf{U}}_{t - 1}} + {{\bf{E}}_t}$$

where X _t is the p × 1 vector of log abundances for each of the p species at time t, A is the p × 1 vector of a values for each species, B is a p × p interaction matrix, whose elements b _ij describe the effect of the density of species j on the per capita growth rate of species i, U _t−1 is the q × 1 vector of covariate values at time t−1, and C is the p × i matrix whose elements c _ij describe the effect of covariate j on species i. The vector of process errors E _t is assumed to be drawn from a multivariate normal distribution with a mean of 0 and covariance matrix ∑.

The MAR framework is a set of multiple regressions solved simultaneously for interacting taxa and environmental covariates, but also accounting for the serial autocorrelation in time series data (see Beisner et al. 2003) through the calculation of autoregression coefficients which depend on the correlated response of one variable to the others in the time series of interest. Therefore, autoregression coefficients depend upon patterns of change in the data. For example, if a given variable does not change, it is not influencing the changes in abundance of taxa in the dataset and the autoregression coefficients will be zero. Furthermore, the MAR model serves as a linear approximation for nonlinear stochastic mutispecies processes (Ives et al. 2003; Hampton et al. 2013).

Stability metrics of community resilience and reactivity can be estimated by MAR using the B species interaction matrix, as defined by Ives et al. (2003). The metrics included in our study were three resilience metrics and two reactivity metrics. Resilience metrics included: (1) det(B)^2/p: the variance of the stationary distribution of community states relative to the variance of the process error that drives the stochasticity as quantified by the B and the ∑ matrices together (lower values are more stable); (2) max (λ _B ): the dominant eigenvalue of B (lower values are more stable); and (3) max(λ _B⊗B): the maximum eigenvalue of the Kronecker product or direct matrix product, B⊗B (higher values = slower return rates). Two reactivity metrics were also described by Ives et al. (2003). The first is −tr(Σ)/tr(V_∞), the trace of the process error over the trace of the p × p covariance matrix, V _∞. This metric describes the degree to which species interactions increase the variance of the stationary distribution relative to the variance of the process (environmental) error; the greater the variance of the individual species abundances relative to the variance of the process error, the higher the reactivity of the system, and hence the lower the stability (less negative values are more reactive). The second is max(λ _B′B ) ₋ 1: the dominant eigenvalue of the matrix B′B minus one, where B′B is the transpose of B multiplied by B. This second metric characterizes the “worst-case” reactivity of a system that depends only on the matrix B (higher values are more reactive).

Appendix 3

Plot of the best-fit model AIC relative to the other top models retained in the analysis for (A) the reference site, (B) the above-dam site, and (C) the below-dam site. The AIC value for the selected best-fit model is denoted with an asterisk along the bottom of the plot (Beisner et al. 2003; Hampton et al. 2013; Figs. 4–6).