In atredennick/community_synchrony: Calculates syunchrony and stability metrics from time series data
\renewcommand\linenumberfont{\normalfont\tiny\sffamily\color{gray}}
library("methods")
library("knitr")
basename <- "ms_communitySynchrony"
opts_chunk$set(fig.path = paste("components/figure/", basename, "-", sep=""),
cache.path = paste("components/cache/", basename, "/", sep=""))
opts_chunk$set(cache = 2)
opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE,
comment = NA, verbose = TRUE, echo=FALSE)
# PDF-based figures
opts_chunk$set(dev='pdf')
\begin{singlespace}
\begin{center}
\large{\textbf{Environmental responses, not species interactions, determine \textcolor{blue}{synchrony of dominant species in semiarid grasslands}}}
\renewcommand*{\thefootnote}{\fnsymbol{footnote}}
\vspace{1em}
\normalsize{Andrew T. Tredennick\textsuperscript{1}, Claire de Mazancourt\textsuperscript{2}, Michel Loreau\textsuperscript{2}, and Peter B. Adler\textsuperscript{1}}
\vspace{1em}
\textit{\small{\textsuperscript{1}Department of Wildland Resources and the Ecology Center, 5230 Old Main Hill, Utah State University, Logan, Utah 84322 USA}}
\textit{\small{\textsuperscript{2}Centre for Biodiversity Theory and Modelling, Theoretical and Experimental Ecology Station, CNRS and Paul Sabatier University, Moulis, 09200, France}}
\end{center}
\vspace{2em}
\textbf{Running Head}: Drivers of species synchrony in \textcolor{blue}{grasslands}
\textbf{Corresponding Author}: \
\hspace{2em}Andrew Tredennick \
\hspace{2em}Department of Wildland Resources and the Ecology Center\
\hspace{2em}Utah State University \
\hspace{2em}5230 Old Main Hill \
\hspace{2em}Logan, Utah 84322 USA \
\hspace{2em}Phone: +1-970-443-1599 \
\hspace{2em}Fax: +1-435-797-3796 \
\hspace{2em}Email: atredenn@gmail.com
\end{singlespace}
\renewcommand*{\thefootnote}{\arabic{footnote}}
\setcounter{footnote}{0}
\newpage{}
\begin{abstract}
Temporal asynchrony among species helps diversity to stabilize ecosystem functioning, but identifying the mechanisms that determine synchrony remains a challenge.
Here, we refine and test theory showing that synchrony depends on three factors: species responses to environmental variation, interspecific interactions, and demographic stochasticity. We then conduct simulation experiments with empirical population models to quantify the relative \textcolor{blue}{influence} of these factors \textcolor{blue}{on the synchrony of dominant species in five semiarid grasslands.}
We found that the average synchrony of per capita growth rates, which can range from 0 (perfect asynchrony) to 1 (perfect synchrony), was higher when environmental variation was present (0.62) rather than absent (0.43).
Removing interspecific interactions and demographic stochasticity had small effects on synchrony.
\textcolor{blue}{For the dominant species in these plant communities}, where species interactions and demographic stochasticity have little influence, synchrony reflects the covariance in species responses to the environment.
\vspace{1em}
\emph{Key words: synchrony, compensatory dynamics, environmental stochasticity, demographic stochasticity, interspecific competition, stability, grassland}
\end{abstract}
\setlength{\parindent}{5ex}
\section{Introduction}
Ecosystems are being transformed by species extinctions [@Cardinale2012], changes in community composition [@Dornelas2014a; @Vellend2013], and anthropogenic environmental change [@Vitousek1997], impacting the provisioning and stability of ecosystem services [@Loreau2001a; @Hooper2005; @Rockstrom2009].
Experiments have provided compelling evidence that decreases in species richness will decrease productivity [@Tilman2001] and the temporal stability of productivity [@Tilman2006; @Hector2010].
\textcolor{blue}{The stabilizing effect of species richness arises from a combination of selection effects and complementarity} [@Loreau2001c].
\textcolor{blue}{Selection effects occur when a dominant species has lower than average temporal variability, which generates a positive effect on ecosystem stability} [e.g., @Grman2010].
\textcolor{blue}{Complementarity occurs when species have unique temporal dynamics, causing their abundances to fluctuate asynchronously, increasing ecosystem stability.}
\textcolor{blue}{The premise of this paper is that understanding the mechanisms driving species' temporal dynamics, and resulting (a)synchrony, is necessary to predict the impacts of global change on ecosystem stability.}
\textcolor{blue}{Asynchronous dynamics, also known as compensatory dynamics} [@Gonzalez2009], \textcolor{blue}{occur whenever species synchrony is not perfect} and result from individual species responding in different ways to environmental fluctuations, random chance events, and/or competitive interactions [@Hector2010; @Isbell2009; @DeMazancourt2013a; @Gross2014].
Species richness affects \textcolor{blue}{the degree of synchrony in a community} because larger species pools are more likely to contain species that respond dissimilarly to environmental conditions, reducing synchrony and increasing stability [@Yachi1999].
\textcolor{blue}{Species richness can also affect synchrony if the strength of species interactions varies systematically with richness, because competition generally increases synchrony and reduces stability} [@Loreau2013; but see @Gross2014].
\textcolor{blue}{The effects of environmental change and species losses on ecosystem stability will depend on whether synchrony is driven by species-specific responses to environmental conditions or interspecific competition} [@Hautier2014].
\textcolor{blue}{If responses to environment are important, then environmental change could alter synchrony and stability.
If competition is important, then the direct effects of environmental change may not affect synchrony and stability, but species gains/losses will.}
\textcolor{blue}{The relative role of environmental responses and competition in driving synchrony in natural plant communities} remains controversial [reviewed in @Gonzalez2009].
One source of the controversy is that quantifying the relative strengths of each driver based on the covariance matrix of species abundances [e.g., @Houlahan2007] is impossible.
This is because observed synchrony can arise from non-unique combinations of factors [@Ranta2008].
For example, weak synchrony of population abundances could reflect positive environmental correlations (synchronizing effect) offset by strong competition (desynchronizing effect), or negative environmental correlations and weak competition.
\textcolor{blue}{Theory can help us resolve this empirical question.}
Recent theoretical work has identified three determinants of species synchrony: environmental stochasticity, interspecific interactions, and demographic stochasticity [@Loreau2008a; @Gonzalez2009; @Loreau2013].
\textcolor{blue}{This theory has been developed by focusing on simple limiting cases in which only one of these three drivers operates.}
For example, in a community composed of large populations (no demographic stochasticity) with weak interspecific interations, community-wide species synchrony should be determined by the covariance of species' responses to the environment [@Loreau2008a].
However, this prediction relies on a relatively simple population model and \textcolor{blue}{depends on two restrictive} assumptions: (i) species' responses to the environment are similar in magnitude and (ii) all species have similar growth rates.
Whether such theoretical predictions hold in natural communities where species differences are unlikely to be symmetrical is unkown because few studies have explicitly tested theory on the drivers of species synchrony in natural communities [@Mutshinda2009; @Thibaut2012], and they did not consider demographic stochasiticity.
An intuitive way to quantify the effects of environmental stochasticity, demographic stochasticity, and interspecific interactions is to remove them one-by-one, and in combination.
\textcolor{blue}{If synchrony changes more when we remove environmental stochasticity than when we remove interspecific interactions, we would conclude that environmental fluctuations are the more important driver.}
In principle, this could be done in an extremely controlled laboratory setting [e.g., @Venail2013], but empirically-based models of interacting populations, fit with data sets from natural communities, offer a practical alternative.
For example, @Mutshinda2009 fit a dynamic population model to several community time series of insect and bird abundances.
They used a statistical technique to decompose temporal variation into competition and environmental components, and found that positively correlated environmental responses among species determined community dynamics.
@Thibaut2012 used a similar approach for reef fish and came to a similar conclusion: environmental responses determine synchrony.
While a major step forward, Mutshinda et al.'s (2009) and Thibaut et al.'s (2012) modeling technique does have some limitations.
\textcolor{blue}{First, although both studies quantified the relative importance of environmental stochasticity and interspecific interactions to explain the observed variation of species synchrony, they did not use the model to quantify how much synchrony would change if each factor were removed.
Second, they relied on popluation abundance data that may or may not reliably capture competitive interactions occuring at the individual level.
Third, fluctuations in abundance may mask the mechanisms that underpin species synchrony.
The synchrony of species' abundances ultimately determines the stability of total community biomass, but the processes that drive species synchrony are most tightly linked to each species' immediate response to environmental conditions and competition} [@Loreau2008a].
\textcolor{blue}{Therefore, we focus on per capita growth rates, which represent a species' immediate response to interannual fluctuations.}
Here, we use multispecies population models fit to long-term demographic data from five semi-arid plant communities to test theory on the drivers of species synchrony (Fig. 1).
Our objectives are to (1) derive and test theoretical predictions of species synchrony and (2) determine the relative influence of environmental stochasticity, interspecific interactions, and demographic stochasticity on \textcolor{blue}{the synchrony of dominant species} in natural plant communities.
\textcolor{blue}{While our focus is limited to dominant species due to data constraints, previous work indicates that the dynamics of dominant species determine ecosystem functioning in grasslands} [@Smith2003; @Bai2004; @Grman2010; @Sasaki2011].
To achieve our objectives, we first refine theory that has been used to predict the effects of species richness on ecosystem stability [@DeMazancourt2013a] and species synchrony [@Loreau2008a] to generate predictions of community-wide species synchrony under two limiting cases.
We then confront our theoretical predictions with simulations from the empirically-based population models.
Second, we use the multi-species population models to perform simulation experiments that isolate the effects of environmental stochasticity, demographic stochasticity, and interspecific interactions on community-wide species synchrony.
Given that our population models capture the essential features of community dynamics important to synchrony (density-dependence, and demographic and environmental stochasticity), and that these models successfully reproduce observed community dynamics [@Chu2015], perturbing the models can reveal the processes that determine synchrony \textcolor{blue}{of dominant species} in our focal grassland communities.
\section{Theoretical Model}
\subsection{The model}
While existing theory has identified \textcolor{blue}{environmental responses, species interactions, and demographic stochasticity as the drivers of the temporal synchrony,} we do not have a simple expression to predict synchrony in a particular community with all factors operating simultaneously.
However, we can derive analytical predictions for species synchrony under special limiting cases.
The limiting case predictions we derive serve as baselines to help us interpret results from empirically-based simulations (described below).
We focus on synchrony of per capita growth rates, rather than abundances, because growth rates represent the instantaneous response of species to the environment and competition, and are less susceptible to the legacy effects of drift and disturbance [@Loreau2008a].
We present equivalent results for synchrony of species abundances in the Appendix S1, and show that they lead to the same overall conclusions as synchrony of per capita growth rates.
Following @Loreau2008a and @DeMazancourt2013a, we define population growth, ignoring observation error, as
\begin{align}
r_i(t) &= \text{ln}N_i(t+1) - \text{ln}N_i(t) \
&= r_{mi} \left[ 1- \frac{N_i(t)+\sum_{j \neq i} \alpha_{ij}N_j(t)} {K_i} + \sigma_{ei}u_{ei}(t) + \frac{\sigma_{di}u_{di}(t)}{\sqrt{N_i(t)}} \right]
\end{align}
\noindent where $N_i(t)$ is the biomass of species i in year t, and $r_i(t)$ is its population growth rate in year t.
$r_{mi}$ is species i's intrinsic rate of increase, $K_i$ is its carrying capacity, and $\alpha_{ij}$ is the interspecific competition coefficient representing the effect of species j on species i.
Environmental stochasticity is incorporated as $\sigma_{ei}u_{ei}(t)$, where $\sigma_{ei}^2$ is the \textcolor{blue}{temporal variance of species \emph{i}'s response to the environment and $u_{ei}(t)$ is a normal variable with unit variance that is independent among species but may be correlated.}
\textcolor{blue}{The product, $\sigma_{ei}u_{ei}(t)$, is species \emph{i}'s environmental response.}
Demographic stochasticity arises from variations in births and deaths among individuals (e.g., same states, different fates), and is included in the model as a first-order, normal approximation [@Lande2003; @DeMazancourt2013a].
$\sigma_{di}^2$ is the demographic variance (\textcolor{blue}{i.e., the intrinsic demographic stochasticity of species \emph{i}}) and $u_{di}(t)$ are independent normal variables with zero mean and unit variance \textcolor{blue}{that allow demographic stochasiticity to vary through time.}
To derive analytical predictions we solved a first-order approximation of Equation 2 [@DeMazancourt2013a and Appendix S1].
Due to the linear approximation approach, our analytical predictions will likely fail in communities where species exhibit large fluctuations due to limit cycles and chaos [@Loreau2008a].
Indeed, one of the advantages of focusing on growth rates rather than abundances is that growth rates are more likely to be well-regulated around an equilibrium value, if the long-term average of a species' growth rate is relatively small (e.g., $r < 2$).
\subsection{Predictions}
Our first prediction assumes no interspecific interactions, no environmental stochasticity, identical intrinsic growth rates, and that demographic stochasticity is operating but all species have identical demographic variances.
This limiting case, $\mathcal{M}{D}$, represents a community where dynamics are driven by demographic stochasticity alone.
Our prediction for the synchrony of per capita growth rates for $\mathcal{M}{D}$, $\phi_{R,\mathcal{M}_{D}}$, is
\begin{equation}
\phi_{R,\mathcal{M}_{D}} = \frac{\sum_i p_i^{-1}}{\left(\sum_i p_i^{-1/2} \right)^2},
\end{equation}
\noindent where $p_i$ is the average frequency of species i, $p_i = N_i/N_T$.
When all species have identical abundances and $p_i = 1/S$, where S is species richness, synchrony equal 1/S [@Loreau2008a].
Our second limiting case assumes only environmental stochasticity is operating ($\mathcal{M}{E}$). Thus, we assume there are no interspecific interactions, demographic stochasticity is absent, intrinsic growth rates are identical, and environmental variance is identical for all species.
Our prediction for the synchrony of per capita growth rates for $\mathcal{M}{E}$, $\phi_{R,\mathcal{M}_{E}}$, is
\begin{equation}
\phi_{R,\mathcal{M}{E}} = \frac{\sum{i,j}\text{cov}(u_{ei},u_{ej})}{S^2},
\end{equation}
\noindent{} where $\text{cov}(u_{ei},u_{ej})$ is the standardized covariance of environmental responses between species i and species j.
Confronting our theoretical predictions with data requires estimates of species dynamics of large populations (no demographic stochasticity) growing in isolation (no interspecific interactions) to calculate the covariance of species' environmental responses.
To estimate environmental responses in natural communities, we turn to our population models built using long-term demographic data.
\section{Empirical Analysis}
\subsection{Materials and methods}
Data.--
We use long-term demographic data from five semiarid grasslands in the western United States [described in detail by @Chu2015].
Each site includes a set of 1-$\text{m}^2$ permanent quadrats within which all individual plants were identified and mapped annually using a pantograph (Hill 1920).
The resulting mapped polygons represent basal cover for grasses and canopy cover for shrubs.
Data come from the Sonoran desert in Arizona [@Anderson2012], sagebrush steppe in Idaho [@Zachmann2010], southern mixed prairie in Kansas [@Adler2007], northern mixed prairie in Montana [@Anderson2011], and Chihuahuan desert in New Mexico [Anderson et al. in preparation, @Chu2015] (Table 1).
Calculating observed synchrony.--
The data consist of records for individual plant size in quadrats for each year.
To obtain estimates of percent cover for each focal species in each year, we summed the individual-level data within quadrats and then averaged percent cover, by species, over all quadrats.
We calculated per capita growth rates as $\text{log}(x_t) - \text{log}(x_{t-1})$, where x is species' percent cover in year t.
Using the community time series of per capita growth rates or percent cover, we calculated community synchrony using the metric of @Loreau2008a in the `synchrony' package [@Gouhier2014] in R [@R2013].
Specifically, we calculated synchrony as
\begin{equation}
\phi_r = \frac{\sigma^{2}{T}}{(\sum{i}\sigma_{r_{i}})^{2}}
\end{equation}
\noindent where $\sigma_{r_{i}}$ is the temporal \textcolor{blue}{standard deviation} of species i's per capita population growth rate ($r_i$) and $\sigma^{2}_{T}$ is the temporal variance of the aggregate community-level growth rate. $\phi$ ranges from 0 at perfect asynchrony to 1 at perfect synchrony [@Loreau2008a].
We use the same equation to calculate observed synchrony of species' percent cover, which we present to relate our results to previous findings, even though we focus on synchrony of growth rates in our model simulations (see below).
Fitting statistical models.--
Vital rate regressions are the building blocks of our dynamic models: an integral projection model (IPM) and an individual-based model (IBM).
We followed the approach of @Chu2015 to fit statistical models for survival, growth, and recruitment (see Appendix S1 for full details).
We modeled survival probability of each genet as function of genet size, temporal variation among years, permanent spatial variation among groups of quadrats, and local neighborhood crowding from conspecific and heterospecific genets.
Regression coefficients for the effect of crowding by each species can be considered a matrix of interaction coefficients whose diagonals represent intraspecific interactions and whose off-diagonals represent interspecific interactions [@Adler2010]. These interaction coefficients can take positive (facilitative) or negative (competitive) values.
We modeled growth as the change in size of a genet from one year to the next, which depends on the same factors as the survival model.
We fit the survival and growth regressions using INLA [@Rue2014], a statistical package for fitting generalized linear mixed effects models via approximate Bayesian inference [@Rue2009], in R [@R2013].
Spatial (quadrat groupings) variation was treated as a random effect on the intercept and temporal (interannual) variation was treated as random effects on the intercept and the effect of genet size in the previous year (Appendix S1).
\textcolor{blue}{Interspecific and intraspecific crowding, which represent species interactions, can be included as fixed effects or as random effects that vary each year.
Strong year-to-year variation in these crowding effects would indicate a statistical interaction between environmental conditions and species interactions.
We tested for such a dynamic by comparing models with and without random year effects on crowding.
Based on the results, we decided to treat crowding as a fixed effect without a temporal component because most 95\% credible intervals for random year effects on crowding broadly overlapped zero and, in a test case, including interannual variation in crowding did not change our results.}
We modeled recruitment at the quadrat scale, rather than the individual scale, because the original data do not attribute new genets to specific parents [@Chu2015].
Our recruitment model assumes that the number of recruits produced in each year follows a negative binomial distribution with the mean dependent on the cover of the parent species, permanent spatial variation among groups, temporal variation among years, and inter- and intraspecific interactions as a function of total species' cover in the quadrat.
We fit the recruitment model using a hierarchical Bayesian approach implemented in JAGS [@Plummer2003] using the `rjags' package [@Plummer2014] in R [@R2013].
Again, temporal and spatial variation were treated as random effects.
Building dynamic multi-species models.--
Once we have fit the vital rate statistical models, building the population models is straightforward.
For the IBM, we initialize the model by randomly assigning plants spatial coordinates, sizes, and species identities until each species achieves a density representative of that observed in the data.
We then project the model forward by using the survival regression to determine whether a genet lives or dies, the growth regression to calculate changes in genet size, and the recruitment regression to add new individuals that are distributed randomly in space.
Crowding is directly calculated at each time step since each genet is spatially referenced.
Environmental stochasticity is not an inherent feature of IBMs, but is easily included since we fit year-specific temporal random effects for each vital rate regression.
To include temporal environmental variation, at each time step we randomly choose a set of estimated survival, growth, and recruitment parameters specific to one observation year.
For all simulations, we ignore the spatial random effect associated with variation among quadrat groups, so our simulations represent an average quadrat for each site.
The IPM uses the same vital rate regressions as the IBM, but it is spatially implicit and does not include demographic stochasticity.
Following @Chu2015, we use a mean field approximation that captures the essential features of spatial patterning to define the crowding index at each time step (Supporting Online Information).
Temporal variation is included in exactly the same way as for the IBM.
For full details on the IPMs we use, see @Chu2015.
Simulation experiments.--
We performed simulation experiments where drivers (environmental stochasticity, demographic stochasticity, or interspecific interactions) were removed one-by-one and in combination.
To remove interspecific interactions, we set the off-diagonals of the interaction matrix for each vital rate regression to zero.
This retains intraspecific interactions, and thus density-dependence, and results in simulations where species are growing in isolation.
We cannot definitively rule out the effects of species interactions on all parameters, meaning that a true monoculture could behave differently than our simulations of a population growing without interspecific competitors.
To remove the effect of a fluctuating environment, we removed the temporal (interannual) random effects from the regression equations.
To remove the effect of demographic stochasticity, we use the IPM rather than the IBM because the IPM does not include demographic stochasticity (demographic stochasticity cannot be removed from the IBM).
Since the effect of demographic stochasticity on population dynamics depends on population size [@Lande2003], we can control the strength of demographic stochasticity by simulating the IBM on areas (e.g. plots) of different size.
Results from an IBM with infinite population size would converge on results from the IPM.
Given computational constraints, the largest landscape we simulate is a 25 $\text{m}^2$ plot.
We conducted the following six simulation experiments:
(1) IBM: All drivers (environmental stochasticity, demographic stochasticity, or interspecific interactions) present;
(2) IPM: Demographic stochasticity removed;
(3) IBM: Environmental stochasticity removed;
(4) IBM: Interspecific interactions removed;
(5) IPM: Interspecific interactions and demographic stochasticity removed;
(6) IBM: Interspecific interactions and environmental stochasticity removed.
\textcolor{blue}{We did not include a simulation with only interspecific interactions because our population models run to deterministic equilibriums in the absence of environmental or demographic stochasticity.}
We ran IPM simulations for 2,000 time steps, after an initial 500 iteration burn-in period.
This allowed species to reach their stable size distribution.
We then calculated yearly per capita growth rates from the simlated time series, and then caluclated the synchrony of species' per capita growth rates over 100 randomly selected contiguous 50 time-step sections.
We ran IBM simulations for 100 time steps, and repeated the simulations 100 times for each simulation experiment.
From those, we retained only the simulations in which no species went extinct due to demographic stochasticity.
Synchrony of per capita growth rates was calculated over the 100 time steps for each no extinction run within a model experiment.
\textcolor{blue}{In the IBM simulations, the strength of demographic stochasticity should weaken as population size increases, meaning that synchrony should be less influenced by demographic stochasticity in large populations compared to small populations.
To explore this effect}, we ran simulations \textcolor{blue}{(4)} and (6) on plot sizes of 1, 4, 9, 16, and 25 $\text{m}^2$.
All other IBM simulations were run on a 25 $\text{m}^2$ landscape \textcolor{blue}{to most closely match the implicit, large spatial scale of the IPM simulations}.
\textcolor{blue}{Our simulations allow us to quantify the relative importance of environmental responses, species interactions, and demographic stochasticity by comparing the simulated values of community-wide species synchrony.}
The simulation experiments also allow us to test our theoretical predictions.
First, in the absence of interspecific interactions and demographic stochasticity, populations can only fluctuate in response to the environment.
Therefore, we can use results from simulation (5) to estimate the covariance of species' responses to the environment ($\text{cov}(u_{ie}, u_{je})$) and parameterize Equation 4.
Parameterizing Equation 3 does not require simulation output because the only parameters are the species' relative abundances.
Second, simulations (5) and (6) represent the simulated version of our limiting case theoretical predictions.
\textcolor{blue}{This approach to testing theoretical predictions may seem circular, but recall we derived the predictions using strict assumptions about equivalence in species' growth rates, environmental response variances, and demographic variances.
Our empirically-based models do not make these assumptions.
Comparisons between parameterized predictions and simulated synchrony reveal whether the assumptions we must make to derive analytical predictions hold in the natural community our model represents.}
\subsection{Results}
\textcolor{blue}{Observed} synchrony of species' per capita growth rates at our study sites range from 0.36 to 0.89 and synchrony of percent cover ranged from 0.15 to 0.92 (Table 2).
Synchrony of per capita growth rates and CV of percent cover are positively correlated (Pearson’s $\rho$ = 0.72, N = 5).
For all five communities, species synchrony from IPM and IBM simulations closely approximated observed synchrony (Fig. S1).
IBM-simulated synchrony is consistently, but only slightly, lower than IPM-simulated synchrony (Fig. S1), likely due to the desynchronizing effect of demographic stochasticity.
Across the five communities, our limiting case predictions closely matched synchrony from the corresponding simulation experiment (Fig. 2 and Table S1).
The correlation between our analytical predictions and simulated synchrony was 0.97 for $\phi_{R,\mathcal{M}D}$ (N = 5) and 0.997 for $\phi{R,\mathcal{M}E}$ (N = 5).
The largest difference between predicted and simulated synchrony was 0.05 in New Mexico for $\phi{R,\mathcal{M}_D}$ (Table S1).
Simulation experiments revealed that removing environmental fluctuations has the largest impact on synchrony, leading to a reduction in synchrony of species growth rates in four out of five communities (Fig. 2).
Removing environmental fluctuations ("No E.S" simulations) decreased synchrony by 33% in Arizona, 48% in Kansas, 39% in Montana, and 40% in New Mexico.
Only in Idaho did removing environmental fluctuations cause an increase in synchrony (Fig. 2), but the effect was small (9% increase; Table S2).
Overall, species' temporal random effects in the statistical vital rate models are positively, but not perfectly, correlated (Table S3).
\textcolor{blue}{These temporal random effects represent environmental responses, meaning that positively correlated temporal random effects indicate positively correlated environmental responses.}
Species interactions are weak in these communities [Table S4 and @Chu2015], and removing interspecific interactions had little effect on synchrony (Fig. 2; "No Comp." simulations).
Removing interspecific interactions caused, at most, a 5% change in synchrony (Fig. 2 and Table S2).
Removing demographic stochasticity ("No D.S." simulations) caused synchrony to increase slightly in all communities (Fig. 2), with an average 6% increase over synchrony from IBM simulations on a 25$\text{m}^2$ area.
In IBM simulations, the desynchronizing effect of demographic stochasticity, which increases as population size decreases, modestly counteracted the synchronizing force of the environment, but not enough to lower synchrony to the level observed when only demographic stochasticity is operating (Fig. 3).
In the largest, 25 $\text{m}^2$ plots, synchrony was driven by environmental stochasticity (e.g., $\mathcal{M}_E$).
At 1 $\text{m}^2$, synchrony reflected \textcolor{blue}{the combined effects of} demographic stochasticity and environmental stochasticity (e.g., \textcolor{blue}{simulated synhcrony fell} between $\mathcal{M}_E$ and $\mathcal{M}_D$).
For context, population sizes increased from an average of 17 individuals per community in 1 $\text{m}^2$ IBM simulations to an average of 357 individuals per community in 25 $\text{m}^2$ IBM simulations.
Results for synchrony of percent cover are qualitatively similar, but more variable and less consistent with analytical predictions and observed synchrony (Appendix S1, Figs. S2-S3).
\section{Discussion}
Our study produced four main findings that were generally consistent across five natural plant communities: (1) limiting-case predictions from the theoretical model were well-supported by simulations from the empirical models; (2) demographic stochasticity decreased community synchrony, as expected by theory, and its effect was largest in small populations; (3) environmental fluctuations increased community synchrony relative to simulations in constant environments because species-specific responses to the environment were positively, though not perfectly, correlated; and (4) interspecific interactions were weak and therefore had little impact on community synchrony.
We also found that analyses based on synchrony of species' percent cover, rather than growth rates, were uninformative (Figs. S2-S3) since the linear approximation required for analytical predictions is a stronger assumption for abundance than growth rates, especially given relatively short time-series (Appendix S1).
Thus, our results provide further evidence that it is difficult to decipher mechanisms of species synchrony from abundance time series, as expected by theory [@Loreau2008a].
Observed synchrony of per capita growth rates was positively correlated with the variability of percent cover across our focal communities, which confirms that we are investigating an important process underlying ecosystem stability.
\subsection{Simulations support theoretical predictions}
Our theoretical predictions were derived from a simple model of population dynamics and required several restrictive assumptions, raising questions about their relevance to natural communities.
For example, the species in our communities do not have equivalent environmental and demographic variances (Figs. S4-S7), as required by our predictions.
However, the theoretical predictions closely matched results from simulations of population models fit to long-term data from natural plant communities (Table S1).
Strong agreement between our analytical predictions and the simulation results should inspire confidence in the ability of simple models to inform our understanding of species synchrony even in complex natural communities, and allows us to place our simulation results within the context of contemporary theory.
\textcolor{blue}{However, whether the theoretical model adequately represents more complex communities remains unknown because our analysis was restricted to dominant species.}
\subsection{Demographic stochasticity decreases synchrony \textcolor{blue}{in small populations}}
In large populations, removing demographic stochasticity had no effect on species synchrony (Fig. 2).
In small populations, demographic stochasticity partially counteracted the synchronizing effects of environmental fluctuations and interspecific interactions on per capita growth rates, in agreement with theory [@Loreau2008a].
\textcolor{blue}{This is shown in Fig. 3, where IBM simulations with environmental forcing and demographic stochasticity have higher synchrony than simulations with only demographic stochasticity.
The differences between simulations are smaller at lower population sizes because as the number of individuals decreases with area, the strength of demographic stochasticity increases, reducing the relative effect of environmental forcing.}
Even in small populations \textcolor{blue}{(e.g., 1 $\text{m}^2$ lanscapes)}, however, demographic stochasticity was not strong enough to compensate the synchronizing effects of environmental fluctuations and match the analytical prediction where only demographic stochasticity is operating ($\mathcal{M}_D$ in Fig. 3).
These results confirm the theoretical argument of @Loreau2008a that independent fluctuations among interacting species in a non-constant environment should be rare.
Only in the Idaho community does synchrony of per capita growth rates approach $\mathcal{M}_D$ in a non-constant environment (Fig. 3).
This is most likely due to the strong effect of demographic stochasticity on the shrub Artemisia tripartita since even a 25 $\text{m}^2$ quadrat would only contain a few individuals of that species.
Our analysis of how demographic stochasticity affects synchrony demonstrates that synchrony depends on the observation area.
As the observation area increases, population size increases and the desynchronizing effect of demographic stochasticity lessens (Fig. 3).
Thus, our results suggest that community-wide species synchrony will increase as the observation area increases, rising from $\mathcal{M}_D$ to $\mathcal{M}_E$.
Such a conclusion assumes, however, that species richness remains constant as observation area increases, which is unlikely [@Taylor1961].
Recent theoretical work has begun to explore the linkage between ecosystem stability and spatial scale [@Wang2014; @Wang2016], and our results suggest that including demographic stochasticity in theoretical models of metacommunity dynamics may be important for understanding the role of species synchrony in determining ecosystem stability across spatial scales.
\subsection{Environmental fluctuations drive community synchrony}
In large populations where interspecific interactions are weak, synchrony is expected to be driven exclusively by environmental fluctuations (Equation 4).
Under such conditions community synchrony should approximately equal the synchrony of species' responses to the environment [@Loreau2008a].
Two lines of evidence lead us to conclude that environmental fluctuations drive species synchrony in our focal plant communities.
First, in our simulation experiments, removing interspecific interactions resulted in no discernible change in community-wide species synchrony of per capita growth rates (Fig. 2).
Second, removing environmental fluctuations from simulations consistently reduced synchrony (Fig. 2).
Our results lead us to conclude that environmental fluctuations, not species interactions, are the primary driver of community-wide species synchrony \textcolor{blue}{among the dominant species} we studied.
Given accumulating evidence that niche differences in natural communities are large [reviewed in @Chu2015], and thus species interactions are likely to be weak, our results may be general in natural plant communities.
In the Idaho community, removing environmental fluctuations did not cause a large decrease in synchrony.
However, that result appears to be an artifact.
Removing environmental variation results in a negative invasion growth rate for A. tripartita.
Although we only analyzed IBM runs in which A. tripartita had not yet gone extinct, it was at much lower abundance than in the other simulation runs.
When we removed A. tripartita from all simulations, the Idaho results conformed with results from all other sites: removing environmental stochasticity caused a significant reduction in species synchrony (Fig. S8).
Our main results for Idaho (Fig. 2), with A. tripartita included, demonstrate how the processes that determine species synchrony interact in complex ways.
A. tripartita has a facilitative effect on each grass species across all vital rates, except for a small competitive effect on H. comata's survival probability (Tables S8-S10).
At the same time, all the perennial grasses have negative effects on each other for each vital rate (Tables S8-S10).
We know synchrony is affected by interspecific competition [@Loreau2008a], but how facilitative effects manifest themselves is unknown.
The interaction of facilitation and competition is clearly capable of having a large effect on species synchrony, and future theoretical efforts should aim to include a wider range of species interactions.
\textcolor{blue}{Environmental responses synchronized dynamics relative to a null expectation of independent species interactions (e.g., "No Comp. + No E.S." simulations in Fig. 2), but observed and simulated synchrony was still less than one in all cases (Fig. 2).
Synchrony was far from complete because of differences in species' responses to interannual environmental variation.
Many studies of ecosystem stability in semiarid grasslands focus on trad-offs among dominant functional types} [@Bai2004; @Sasaki2011].
\textcolor{blue}{Such groupings are based on the idea that ecologically-similar species will have similar responses to environmental fluctuations.
At first glance our results may appear to support the grouping of perennial grasses in one functional type because their environmental responses were positively correlated.
However, even though environmental responses among the dominant species we studied were similar, they were dissimilar enough to cause synchrony to be less than perfect (Fig. 2).
The subtle differences among ecologically-similar dominant species do impact species synchrony and, ultimately, ecosystem stability.
Ignoring such differences could mask important dynamics that underpin ecosystem functioning.}
\subsection{Interspecific interactions had little impact on community synchrony}
We expected community synchrony of per capita growth rates to decrease when we removed interspecific interactions [@Loreau2008a].
We found that community synchrony was virtually indistinguishable between simulations with and without interspecific interactions (Fig. 2).
The lack of an effect of interspecific interactions on synchrony is in contrast to a large body of theoretical work that predicts a strong role for competition in creating compensatory dynamics [@Tilman1988] and a recent empirical analysis [@Gross2014].
Our results do not contradict the idea that competition can lead to compensatory dynamics, but they do highlight the fact that interspecific competition must be relatively strong to influence species synchrony.
The communities we analyzed are composed of species with very little niche overlap [@Chu2015] and weak interspecific interactions (Tables S1 and S3-S17).
Mechanistic consumer-resource models [@Lehman2000] and phenomenological Lotka-Volterra models [@Lehman2000; @Loreau2013] both confirm that the effect of competition on species synchrony diminishes as niche overlap decreases.
In that sense, our results are not surprising: interspecific interactions are weak, so of course removing them does not affect synchrony.
Our conclusion that species interactions have little impact on synchrony only applies to single trophic level communities.
Species interactions almost certainly play a strong role in multi-trophic communities where factors such as resource overlap [@Vasseur2007], dispersal [@Gouhier2010], and the strength of top-down control [@Bauer2014] are all likely to affect community synchrony.
\section{Conclusions}
Species-specific responses to temporally fluctuating environmental conditions is an important mechanism underlying asynchronous population dynamics and, in turn, ecosystem stability [@Loreau2013].
When we removed environmental variation, we found that synchrony decreased in four out of the five grassland communities we studied (Fig. 2).
A tempting conclusion is that our study confirms that compensatory dynamics are rare in natural communities, and that ecologically-similar species will exhibit synchronous dynamics [e.g., @Houlahan2007].
Such a conclusion misses an important subtlety.
The perennial grasses we studied do have similar responses to the environment (Table S2), which will tend to synchronize dynamics.
However, if community-wide synchrony \textcolor{blue}{among dominant species} is less than perfect, as it is in all our focal communities, some degree of compensatory dynamics must be present [@Loreau2008a].
\textcolor{blue}{Even ecologically-similar species, which are sometimes aggregated into functional groups, have environmental responses that are dissimilar enough to limit synchrony.
Subtle differences among dominant species ultimately determine ecosystem stability and should not be ignored.}
\section{Acknowledgements}
This work was funded by the National Science Foundation through a Postdoctoral Research Fellowship in Biology to ATT (DBI-1400370) and a CAREER award to PBA (DEB-1054040). CdM and ML were supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41). We thank the original mappers of the permanent quadrats at each site and the digitizers in the Adler lab, without whom this work would not have been possible. Compute, storage, and other resources from the Division of Research Computing in the Office of Research and Graduate Studies at Utah State University are gratefully acknowledged. We also thank Patrick Venail and several anonymous reviewers who provided critical feedback that improved the manuscript.
\setlength{\parindent}{0ex}
\newpage{}
\section{Tables}
\singlespacing
Table 1: Site descriptions and focal species.
\footnotesize
| Site Name | Biome | Location (Lat, Lon) | Obs. Years | Species |
| --------- | ----- | ------------------- | ---------- | ------- |
| New Mexico | Chihuahuan Desert | 32.62° N, 106.67° W | 1915-1950 | Bouteloua eriopoda |
| | | | | Sporobolus flexuosus |
| Arizona | Sonoran Desert | 31°50' N, 110°53' W | 1915-1933 | Bouteloua eriopoda |
| | | | | Bouteloua rothrockii |
| Kansas | Southern mixed prairie | 38.8° N, 99.3° W | 1932-1972 | Bouteloua curtipendula |
| | | | | Bouteloua hirsuta |
| | | | | Schizachyrium scoparium |
| Montana | Northern mixed prairie | 46°19' N, 105°48' W | 1926-1957 | Bouteloua gracilis |
| | | | | Hesperostipa comata |
| | | | | Pascopyrum smithii |
| | | | | Poa secunda |
| Idaho | Sagebrush steppe | 44.2° N, 112.1° W | 1926-1957 | Artemisia tripartita |
| | | | | Pseudoroegneria spicata |
| | | | | Hesperostipa comata |
| | | | | Poa secunda |
\normalsize
\pagebreak{}
library(communitySynchrony)
library(plyr)
library(reshape2)
library(synchrony)
library(ggplot2)
library(xtable)
library(gridExtra)
library(RColorBrewer)
site <- "Kansas"
spp_list <- c("BOCU","BOHI","SCSC")
num_spp <- length(spp_list)
ks_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
ks_data <- rbind(ks_data, spp_data)
} #end species looping for raw data
ks_data <- ks_data[2:nrow(ks_data),] #remove first NA row
tmp1<-which(ks_data$quad_data=="q25" & (ks_data$year<35 | ks_data$year>62))
tmp2<-which(ks_data$quad_data=="q27")
tmp3<-which(ks_data$quad=="q28")
tmp4<-which(ks_data$quad=="q30")
tmp5<-which(ks_data$quad=="q31" & (ks_data$year<35 | ks_data$year>39))
tmp6<-which(ks_data$quad=="q32" & (ks_data$year<35 | ks_data$year>41))
tmp<-c(tmp1,tmp2,tmp3,tmp4,tmp5,tmp6)
ks_data<-ks_data[-tmp,]
# exclude the records later than 1968, to keep the same random year effect...
ks_data<-subset(ks_data,year<68)
site <- "Idaho"
spp_list <- sort(c("PSSP","HECO","POSE","ARTR"))
num_spp <- length(spp_list)
id_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
id_data <- rbind(id_data, spp_data)
} #end species looping for raw data
id_data <- id_data[2:nrow(id_data),] #remove first NA row
# id_data <- subset(id_data, species!="ARTR") #take out the shrub
site <- "Montana"
spp_list <- sort(c("BOGR","HECO","PASM","POSE"))
num_spp <- length(spp_list)
mt_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
mt_data <- rbind(mt_data, spp_data)
} #end species looping for raw data
mt_data <- mt_data[2:nrow(mt_data),] #remove first NA row
site <- "NewMexico"
spp_list <- sort(c("BOER","SPFL"))
num_spp <- length(spp_list)
nm_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
nm_data <- rbind(nm_data, spp_data)
} #end species looping for raw data
nm_data <- nm_data[2:nrow(nm_data),] #remove first NA row
site <- "Arizona"
spp_list <- sort(c("BOER","BORO"))
num_spp <- length(spp_list)
az_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
az_data <- rbind(az_data, spp_data)
} #end species looping for raw data
az_data <- az_data[2:nrow(az_data),] #remove first NA row
out_ks <- get_comm_synchrony(ts_data = ks_data)
out_id <- get_comm_synchrony(ts_data = id_data)
out_mt <- get_comm_synchrony(ts_data = mt_data)
out_nm <- get_comm_synchrony(ts_data = nm_data)
out_az <- get_comm_synchrony(ts_data = az_data)
get_table_metrics <- function(output){
tmp <- output
exp_tmp <- round(tmp$pgr_expected_synch_ind_flucts, 2)
obs_tmp <- round(as.numeric(tmp$pgr_synchrony[1]), 2)
demo_tmp <- obs_tmp - exp_tmp
exp_cover <- round(as.numeric(tmp$cover_expected_synch_ind_flucts), 2)
obs_cover <- round(as.numeric(tmp$abund_synchrony[1]), 2)
demo_cover <- obs_cover - exp_cover
plant_size <- tmp$avg_size
cover_var <- tmp$variability
return(c(obs_tmp, exp_tmp, demo_tmp, obs_cover, exp_cover,
demo_cover, cover_var))
}
table_data <- as.data.frame(rbind(get_table_metrics(out_az),
get_table_metrics(out_id),
get_table_metrics(out_ks),
get_table_metrics(out_mt),
get_table_metrics(out_nm)))
table_data$Site <- c("Arizona", "Idaho", "Kansas", "Montana", "New Mexico")
alltable_data <- table_data
colnames(table_data) <- c("obs_gr_synch", "demo_only_exp_grsynch",
"diff_noneed", "obs_cover_synch",
"demo_only_exp_coversynch", "diffcover_noneed",
"cover_variance", "site")
saveRDS(table_data, "../results/calculate_site_metrics.RDS")
table_data <- alltable_data[,c("Site", "V1", "V4", "V7")]
table_data$spp <- c(length(unique(az_data$species)),
length(unique(id_data$species)),
length(unique(ks_data$species)),
length(unique(mt_data$species)),
length(unique(nm_data$species)))
colnames(table_data) <- c("Site",
"$\\phi_{R}$",
"$\\phi_{C}$",
"CV of Total Cover",
"Species richness")
myorder <- c(5,1,3,4,2)
xtable_data <- table_data[myorder,]
# colnames(table_data) <- c("Site",
# "$\\phi_{r}$",
# "$\\text{E}(\\phi_{r})$",
# "$\\text{E}(\\sigma_{r,c} + \\sigma_{r,e})$",
# "$\\phi_{C}$",
# "$\\text{E}(\\phi_{C})$",
# "$\\text{E}(\\sigma_{C,c} + \\sigma_{C,e})$",
# "CV of Total Cover",
# "Species richness")
# table_data[,3] <- round(1/table_data[,ncol(table_data)],2) # reset expected sybchrony to 1/S
# table_data[,4] <- table_data[,2] - table_data[,3] # calculate obs - exp synchrony
synch_cap <- "Observed synchrony among species' per capita growth rates ($\\phi_{R}$), observed synchrony among species' percent cover ($\\phi_{C}$), the coefficient of variation of total community cover, and species richness for each community. Species richness values reflect the number of species analyzed from the community, not the actual richness."
print(xtable(xtable_data, caption = synch_cap),
caption.placement="top",
include.rownames = F,
sanitize.colnames.function = identity,
comment=FALSE,
size="normalsize")
sink(file = "../results/correlation_Rsynch_cv.txt")
cat("correlation between synchrony of growth rates and cover CV is:")
cor(table_data$`CV of Total Cover`, table_data[,2])
Sys.Date()
sink()
\pagebreak{}
\section{Figures}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/synchrony_flowchart_bw.png}
\caption{Diagram of our coupled theoretical-empirical approach. We followed this workflow for each of our five focal communities.}
\end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/formatted_figure1.png}
\caption{Community-wide species synchrony of per capita growth rates from model simulation experiments. Synchrony of species' growth rates for each study area are from simulation experiments with demographic stochasticity, environmental stochasticity, and interspecific interactions present (All Drivers''), demographic stochasticity removed (No D.S.''), environmental stochasticity removed (No E.S.''), interspecific interactions removed (No Comp.''), interspecific interactions and demographic stochasticity removed (No Comp. + No D.S.''), and interspecific interactions and environmental stochasticity removed (No Comp. + No E.S.''). Abbreviations within the bars for the New Mexico site indicate whether the IBM or IPM was used for a particular simulation. Error bars represent the 2.5\% and 97.5\% quantiles from model simulations. All IBM simulations shown in this figure were run on a 25 $\text{m}^2$ virtual landscape. Points show observed and predicted synchrony aligned with the model simulation that corresponds with each observation or analytical predicion.}
\end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/formatted_figure2.png}
\caption{Synchrony of species' growth rates for each study area from IBM simulations across different landscape sizes when only demographic stochasticity is present (No Comp. + No E.S. (D.S. Only)'') and when environmental stochasticity is also present (No Comp. (D.S. + E.S.)''). The horizontal lines show the analytical predictions $\mathcal{M}_D$ (dashed line) and $\mathcal{M}_E$ (dotted line). The strength of demographic stochasticity decreases as landscape size increases because population sizes also increase. Theoretically, No Comp. + No E.S. (D.S. Only)'' simulations should remain constant across landscape size, whereasNo Comp. (D.S. + E.S.)'' simulations should shift from the $\mathcal{M}_D$ prediction to the $\mathcal{M}_E$ prediction as landscape size, and thus population size, increases, but only if demographic stochasticity it strong enough to counteract environmental forcing. Error bars represent the 2.5\% and 97.5\% quantiles from model simulations.}
\end{figure}
\pagebreak{}
References
atredennick/community_synchrony documentation built on May 10, 2019, 2:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
\renewcommand\linenumberfont{\normalfont\tiny\sffamily\color{gray}}
library("methods") library("knitr") basename <- "ms_communitySynchrony" opts_chunk$set(fig.path = paste("components/figure/", basename, "-", sep=""), cache.path = paste("components/cache/", basename, "/", sep="")) opts_chunk$set(cache = 2) opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, comment = NA, verbose = TRUE, echo=FALSE) # PDF-based figures opts_chunk$set(dev='pdf')
\begin{singlespace}
\begin{center} \large{\textbf{Environmental responses, not species interactions, determine \textcolor{blue}{synchrony of dominant species in semiarid grasslands}}}
\renewcommand*{\thefootnote}{\fnsymbol{footnote}}
\vspace{1em}
\normalsize{Andrew T. Tredennick\textsuperscript{1}, Claire de Mazancourt\textsuperscript{2}, Michel Loreau\textsuperscript{2}, and Peter B. Adler\textsuperscript{1}}
\vspace{1em}
\textit{\small{\textsuperscript{1}Department of Wildland Resources and the Ecology Center, 5230 Old Main Hill, Utah State University, Logan, Utah 84322 USA}}
\textit{\small{\textsuperscript{2}Centre for Biodiversity Theory and Modelling, Theoretical and Experimental Ecology Station, CNRS and Paul Sabatier University, Moulis, 09200, France}}
\end{center}
\vspace{2em}
\textbf{Running Head}: Drivers of species synchrony in \textcolor{blue}{grasslands}
\textbf{Corresponding Author}: \
\hspace{2em}Andrew Tredennick \
\hspace{2em}Department of Wildland Resources and the Ecology Center\
\hspace{2em}Utah State University \
\hspace{2em}5230 Old Main Hill \
\hspace{2em}Logan, Utah 84322 USA \
\hspace{2em}Phone: +1-970-443-1599 \
\hspace{2em}Fax: +1-435-797-3796 \
\hspace{2em}Email: atredenn@gmail.com
\end{singlespace} \renewcommand*{\thefootnote}{\arabic{footnote}} \setcounter{footnote}{0}
\newpage{}
\begin{abstract} Temporal asynchrony among species helps diversity to stabilize ecosystem functioning, but identifying the mechanisms that determine synchrony remains a challenge. Here, we refine and test theory showing that synchrony depends on three factors: species responses to environmental variation, interspecific interactions, and demographic stochasticity. We then conduct simulation experiments with empirical population models to quantify the relative \textcolor{blue}{influence} of these factors \textcolor{blue}{on the synchrony of dominant species in five semiarid grasslands.} We found that the average synchrony of per capita growth rates, which can range from 0 (perfect asynchrony) to 1 (perfect synchrony), was higher when environmental variation was present (0.62) rather than absent (0.43). Removing interspecific interactions and demographic stochasticity had small effects on synchrony. \textcolor{blue}{For the dominant species in these plant communities}, where species interactions and demographic stochasticity have little influence, synchrony reflects the covariance in species responses to the environment.
\vspace{1em}
\emph{Key words: synchrony, compensatory dynamics, environmental stochasticity, demographic stochasticity, interspecific competition, stability, grassland}
\end{abstract}
\setlength{\parindent}{5ex}
\section{Introduction} Ecosystems are being transformed by species extinctions [@Cardinale2012], changes in community composition [@Dornelas2014a; @Vellend2013], and anthropogenic environmental change [@Vitousek1997], impacting the provisioning and stability of ecosystem services [@Loreau2001a; @Hooper2005; @Rockstrom2009]. Experiments have provided compelling evidence that decreases in species richness will decrease productivity [@Tilman2001] and the temporal stability of productivity [@Tilman2006; @Hector2010]. \textcolor{blue}{The stabilizing effect of species richness arises from a combination of selection effects and complementarity} [@Loreau2001c]. \textcolor{blue}{Selection effects occur when a dominant species has lower than average temporal variability, which generates a positive effect on ecosystem stability} [e.g., @Grman2010]. \textcolor{blue}{Complementarity occurs when species have unique temporal dynamics, causing their abundances to fluctuate asynchronously, increasing ecosystem stability.} \textcolor{blue}{The premise of this paper is that understanding the mechanisms driving species' temporal dynamics, and resulting (a)synchrony, is necessary to predict the impacts of global change on ecosystem stability.}
\textcolor{blue}{Asynchronous dynamics, also known as compensatory dynamics} [@Gonzalez2009], \textcolor{blue}{occur whenever species synchrony is not perfect} and result from individual species responding in different ways to environmental fluctuations, random chance events, and/or competitive interactions [@Hector2010; @Isbell2009; @DeMazancourt2013a; @Gross2014]. Species richness affects \textcolor{blue}{the degree of synchrony in a community} because larger species pools are more likely to contain species that respond dissimilarly to environmental conditions, reducing synchrony and increasing stability [@Yachi1999]. \textcolor{blue}{Species richness can also affect synchrony if the strength of species interactions varies systematically with richness, because competition generally increases synchrony and reduces stability} [@Loreau2013; but see @Gross2014].
\textcolor{blue}{The effects of environmental change and species losses on ecosystem stability will depend on whether synchrony is driven by species-specific responses to environmental conditions or interspecific competition} [@Hautier2014]. \textcolor{blue}{If responses to environment are important, then environmental change could alter synchrony and stability. If competition is important, then the direct effects of environmental change may not affect synchrony and stability, but species gains/losses will.}
\textcolor{blue}{The relative role of environmental responses and competition in driving synchrony in natural plant communities} remains controversial [reviewed in @Gonzalez2009]. One source of the controversy is that quantifying the relative strengths of each driver based on the covariance matrix of species abundances [e.g., @Houlahan2007] is impossible. This is because observed synchrony can arise from non-unique combinations of factors [@Ranta2008]. For example, weak synchrony of population abundances could reflect positive environmental correlations (synchronizing effect) offset by strong competition (desynchronizing effect), or negative environmental correlations and weak competition.
\textcolor{blue}{Theory can help us resolve this empirical question.} Recent theoretical work has identified three determinants of species synchrony: environmental stochasticity, interspecific interactions, and demographic stochasticity [@Loreau2008a; @Gonzalez2009; @Loreau2013]. \textcolor{blue}{This theory has been developed by focusing on simple limiting cases in which only one of these three drivers operates.} For example, in a community composed of large populations (no demographic stochasticity) with weak interspecific interations, community-wide species synchrony should be determined by the covariance of species' responses to the environment [@Loreau2008a]. However, this prediction relies on a relatively simple population model and \textcolor{blue}{depends on two restrictive} assumptions: (i) species' responses to the environment are similar in magnitude and (ii) all species have similar growth rates. Whether such theoretical predictions hold in natural communities where species differences are unlikely to be symmetrical is unkown because few studies have explicitly tested theory on the drivers of species synchrony in natural communities [@Mutshinda2009; @Thibaut2012], and they did not consider demographic stochasiticity.
An intuitive way to quantify the effects of environmental stochasticity, demographic stochasticity, and interspecific interactions is to remove them one-by-one, and in combination. \textcolor{blue}{If synchrony changes more when we remove environmental stochasticity than when we remove interspecific interactions, we would conclude that environmental fluctuations are the more important driver.} In principle, this could be done in an extremely controlled laboratory setting [e.g., @Venail2013], but empirically-based models of interacting populations, fit with data sets from natural communities, offer a practical alternative. For example, @Mutshinda2009 fit a dynamic population model to several community time series of insect and bird abundances. They used a statistical technique to decompose temporal variation into competition and environmental components, and found that positively correlated environmental responses among species determined community dynamics. @Thibaut2012 used a similar approach for reef fish and came to a similar conclusion: environmental responses determine synchrony.
While a major step forward, Mutshinda et al.'s (2009) and Thibaut et al.'s (2012) modeling technique does have some limitations. \textcolor{blue}{First, although both studies quantified the relative importance of environmental stochasticity and interspecific interactions to explain the observed variation of species synchrony, they did not use the model to quantify how much synchrony would change if each factor were removed. Second, they relied on popluation abundance data that may or may not reliably capture competitive interactions occuring at the individual level. Third, fluctuations in abundance may mask the mechanisms that underpin species synchrony. The synchrony of species' abundances ultimately determines the stability of total community biomass, but the processes that drive species synchrony are most tightly linked to each species' immediate response to environmental conditions and competition} [@Loreau2008a]. \textcolor{blue}{Therefore, we focus on per capita growth rates, which represent a species' immediate response to interannual fluctuations.}
Here, we use multispecies population models fit to long-term demographic data from five semi-arid plant communities to test theory on the drivers of species synchrony (Fig. 1). Our objectives are to (1) derive and test theoretical predictions of species synchrony and (2) determine the relative influence of environmental stochasticity, interspecific interactions, and demographic stochasticity on \textcolor{blue}{the synchrony of dominant species} in natural plant communities. \textcolor{blue}{While our focus is limited to dominant species due to data constraints, previous work indicates that the dynamics of dominant species determine ecosystem functioning in grasslands} [@Smith2003; @Bai2004; @Grman2010; @Sasaki2011].
To achieve our objectives, we first refine theory that has been used to predict the effects of species richness on ecosystem stability [@DeMazancourt2013a] and species synchrony [@Loreau2008a] to generate predictions of community-wide species synchrony under two limiting cases. We then confront our theoretical predictions with simulations from the empirically-based population models. Second, we use the multi-species population models to perform simulation experiments that isolate the effects of environmental stochasticity, demographic stochasticity, and interspecific interactions on community-wide species synchrony. Given that our population models capture the essential features of community dynamics important to synchrony (density-dependence, and demographic and environmental stochasticity), and that these models successfully reproduce observed community dynamics [@Chu2015], perturbing the models can reveal the processes that determine synchrony \textcolor{blue}{of dominant species} in our focal grassland communities.
\section{Theoretical Model} \subsection{The model} While existing theory has identified \textcolor{blue}{environmental responses, species interactions, and demographic stochasticity as the drivers of the temporal synchrony,} we do not have a simple expression to predict synchrony in a particular community with all factors operating simultaneously. However, we can derive analytical predictions for species synchrony under special limiting cases. The limiting case predictions we derive serve as baselines to help us interpret results from empirically-based simulations (described below). We focus on synchrony of per capita growth rates, rather than abundances, because growth rates represent the instantaneous response of species to the environment and competition, and are less susceptible to the legacy effects of drift and disturbance [@Loreau2008a]. We present equivalent results for synchrony of species abundances in the Appendix S1, and show that they lead to the same overall conclusions as synchrony of per capita growth rates.
Following @Loreau2008a and @DeMazancourt2013a, we define population growth, ignoring observation error, as
\begin{align} r_i(t) &= \text{ln}N_i(t+1) - \text{ln}N_i(t) \ &= r_{mi} \left[ 1- \frac{N_i(t)+\sum_{j \neq i} \alpha_{ij}N_j(t)} {K_i} + \sigma_{ei}u_{ei}(t) + \frac{\sigma_{di}u_{di}(t)}{\sqrt{N_i(t)}} \right] \end{align}
\noindent where $N_i(t)$ is the biomass of species i in year t, and $r_i(t)$ is its population growth rate in year t. $r_{mi}$ is species i's intrinsic rate of increase, $K_i$ is its carrying capacity, and $\alpha_{ij}$ is the interspecific competition coefficient representing the effect of species j on species i. Environmental stochasticity is incorporated as $\sigma_{ei}u_{ei}(t)$, where $\sigma_{ei}^2$ is the \textcolor{blue}{temporal variance of species \emph{i}'s response to the environment and $u_{ei}(t)$ is a normal variable with unit variance that is independent among species but may be correlated.} \textcolor{blue}{The product, $\sigma_{ei}u_{ei}(t)$, is species \emph{i}'s environmental response.} Demographic stochasticity arises from variations in births and deaths among individuals (e.g., same states, different fates), and is included in the model as a first-order, normal approximation [@Lande2003; @DeMazancourt2013a]. $\sigma_{di}^2$ is the demographic variance (\textcolor{blue}{i.e., the intrinsic demographic stochasticity of species \emph{i}}) and $u_{di}(t)$ are independent normal variables with zero mean and unit variance \textcolor{blue}{that allow demographic stochasiticity to vary through time.}
To derive analytical predictions we solved a first-order approximation of Equation 2 [@DeMazancourt2013a and Appendix S1]. Due to the linear approximation approach, our analytical predictions will likely fail in communities where species exhibit large fluctuations due to limit cycles and chaos [@Loreau2008a]. Indeed, one of the advantages of focusing on growth rates rather than abundances is that growth rates are more likely to be well-regulated around an equilibrium value, if the long-term average of a species' growth rate is relatively small (e.g., $r < 2$).
\subsection{Predictions} Our first prediction assumes no interspecific interactions, no environmental stochasticity, identical intrinsic growth rates, and that demographic stochasticity is operating but all species have identical demographic variances. This limiting case, $\mathcal{M}{D}$, represents a community where dynamics are driven by demographic stochasticity alone. Our prediction for the synchrony of per capita growth rates for $\mathcal{M}{D}$, $\phi_{R,\mathcal{M}_{D}}$, is
\begin{equation} \phi_{R,\mathcal{M}_{D}} = \frac{\sum_i p_i^{-1}}{\left(\sum_i p_i^{-1/2} \right)^2}, \end{equation}
\noindent where $p_i$ is the average frequency of species i, $p_i = N_i/N_T$. When all species have identical abundances and $p_i = 1/S$, where S is species richness, synchrony equal 1/S [@Loreau2008a].
Our second limiting case assumes only environmental stochasticity is operating ($\mathcal{M}{E}$). Thus, we assume there are no interspecific interactions, demographic stochasticity is absent, intrinsic growth rates are identical, and environmental variance is identical for all species. Our prediction for the synchrony of per capita growth rates for $\mathcal{M}{E}$, $\phi_{R,\mathcal{M}_{E}}$, is
\begin{equation} \phi_{R,\mathcal{M}{E}} = \frac{\sum{i,j}\text{cov}(u_{ei},u_{ej})}{S^2}, \end{equation}
\noindent{} where $\text{cov}(u_{ei},u_{ej})$ is the standardized covariance of environmental responses between species i and species j.
Confronting our theoretical predictions with data requires estimates of species dynamics of large populations (no demographic stochasticity) growing in isolation (no interspecific interactions) to calculate the covariance of species' environmental responses. To estimate environmental responses in natural communities, we turn to our population models built using long-term demographic data.
\section{Empirical Analysis} \subsection{Materials and methods}
Data.--
We use long-term demographic data from five semiarid grasslands in the western United States [described in detail by @Chu2015]. Each site includes a set of 1-$\text{m}^2$ permanent quadrats within which all individual plants were identified and mapped annually using a pantograph (Hill 1920). The resulting mapped polygons represent basal cover for grasses and canopy cover for shrubs. Data come from the Sonoran desert in Arizona [@Anderson2012], sagebrush steppe in Idaho [@Zachmann2010], southern mixed prairie in Kansas [@Adler2007], northern mixed prairie in Montana [@Anderson2011], and Chihuahuan desert in New Mexico [Anderson et al. in preparation, @Chu2015] (Table 1).
Calculating observed synchrony.--
The data consist of records for individual plant size in quadrats for each year. To obtain estimates of percent cover for each focal species in each year, we summed the individual-level data within quadrats and then averaged percent cover, by species, over all quadrats. We calculated per capita growth rates as $\text{log}(x_t) - \text{log}(x_{t-1})$, where x is species' percent cover in year t. Using the community time series of per capita growth rates or percent cover, we calculated community synchrony using the metric of @Loreau2008a in the `synchrony' package [@Gouhier2014] in R [@R2013]. Specifically, we calculated synchrony as
\begin{equation} \phi_r = \frac{\sigma^{2}{T}}{(\sum{i}\sigma_{r_{i}})^{2}} \end{equation}
\noindent where $\sigma_{r_{i}}$ is the temporal \textcolor{blue}{standard deviation} of species i's per capita population growth rate ($r_i$) and $\sigma^{2}_{T}$ is the temporal variance of the aggregate community-level growth rate. $\phi$ ranges from 0 at perfect asynchrony to 1 at perfect synchrony [@Loreau2008a]. We use the same equation to calculate observed synchrony of species' percent cover, which we present to relate our results to previous findings, even though we focus on synchrony of growth rates in our model simulations (see below).
Fitting statistical models.--
Vital rate regressions are the building blocks of our dynamic models: an integral projection model (IPM) and an individual-based model (IBM). We followed the approach of @Chu2015 to fit statistical models for survival, growth, and recruitment (see Appendix S1 for full details). We modeled survival probability of each genet as function of genet size, temporal variation among years, permanent spatial variation among groups of quadrats, and local neighborhood crowding from conspecific and heterospecific genets. Regression coefficients for the effect of crowding by each species can be considered a matrix of interaction coefficients whose diagonals represent intraspecific interactions and whose off-diagonals represent interspecific interactions [@Adler2010]. These interaction coefficients can take positive (facilitative) or negative (competitive) values. We modeled growth as the change in size of a genet from one year to the next, which depends on the same factors as the survival model. We fit the survival and growth regressions using INLA [@Rue2014], a statistical package for fitting generalized linear mixed effects models via approximate Bayesian inference [@Rue2009], in R [@R2013]. Spatial (quadrat groupings) variation was treated as a random effect on the intercept and temporal (interannual) variation was treated as random effects on the intercept and the effect of genet size in the previous year (Appendix S1).
\textcolor{blue}{Interspecific and intraspecific crowding, which represent species interactions, can be included as fixed effects or as random effects that vary each year. Strong year-to-year variation in these crowding effects would indicate a statistical interaction between environmental conditions and species interactions. We tested for such a dynamic by comparing models with and without random year effects on crowding. Based on the results, we decided to treat crowding as a fixed effect without a temporal component because most 95\% credible intervals for random year effects on crowding broadly overlapped zero and, in a test case, including interannual variation in crowding did not change our results.}
We modeled recruitment at the quadrat scale, rather than the individual scale, because the original data do not attribute new genets to specific parents [@Chu2015]. Our recruitment model assumes that the number of recruits produced in each year follows a negative binomial distribution with the mean dependent on the cover of the parent species, permanent spatial variation among groups, temporal variation among years, and inter- and intraspecific interactions as a function of total species' cover in the quadrat. We fit the recruitment model using a hierarchical Bayesian approach implemented in JAGS [@Plummer2003] using the `rjags' package [@Plummer2014] in R [@R2013]. Again, temporal and spatial variation were treated as random effects.
Building dynamic multi-species models.--
Once we have fit the vital rate statistical models, building the population models is straightforward. For the IBM, we initialize the model by randomly assigning plants spatial coordinates, sizes, and species identities until each species achieves a density representative of that observed in the data. We then project the model forward by using the survival regression to determine whether a genet lives or dies, the growth regression to calculate changes in genet size, and the recruitment regression to add new individuals that are distributed randomly in space. Crowding is directly calculated at each time step since each genet is spatially referenced. Environmental stochasticity is not an inherent feature of IBMs, but is easily included since we fit year-specific temporal random effects for each vital rate regression. To include temporal environmental variation, at each time step we randomly choose a set of estimated survival, growth, and recruitment parameters specific to one observation year. For all simulations, we ignore the spatial random effect associated with variation among quadrat groups, so our simulations represent an average quadrat for each site.
The IPM uses the same vital rate regressions as the IBM, but it is spatially implicit and does not include demographic stochasticity. Following @Chu2015, we use a mean field approximation that captures the essential features of spatial patterning to define the crowding index at each time step (Supporting Online Information). Temporal variation is included in exactly the same way as for the IBM. For full details on the IPMs we use, see @Chu2015.
Simulation experiments.--
We performed simulation experiments where drivers (environmental stochasticity, demographic stochasticity, or interspecific interactions) were removed one-by-one and in combination. To remove interspecific interactions, we set the off-diagonals of the interaction matrix for each vital rate regression to zero. This retains intraspecific interactions, and thus density-dependence, and results in simulations where species are growing in isolation. We cannot definitively rule out the effects of species interactions on all parameters, meaning that a true monoculture could behave differently than our simulations of a population growing without interspecific competitors. To remove the effect of a fluctuating environment, we removed the temporal (interannual) random effects from the regression equations. To remove the effect of demographic stochasticity, we use the IPM rather than the IBM because the IPM does not include demographic stochasticity (demographic stochasticity cannot be removed from the IBM). Since the effect of demographic stochasticity on population dynamics depends on population size [@Lande2003], we can control the strength of demographic stochasticity by simulating the IBM on areas (e.g. plots) of different size. Results from an IBM with infinite population size would converge on results from the IPM. Given computational constraints, the largest landscape we simulate is a 25 $\text{m}^2$ plot.
We conducted the following six simulation experiments: (1) IBM: All drivers (environmental stochasticity, demographic stochasticity, or interspecific interactions) present; (2) IPM: Demographic stochasticity removed; (3) IBM: Environmental stochasticity removed; (4) IBM: Interspecific interactions removed; (5) IPM: Interspecific interactions and demographic stochasticity removed; (6) IBM: Interspecific interactions and environmental stochasticity removed. \textcolor{blue}{We did not include a simulation with only interspecific interactions because our population models run to deterministic equilibriums in the absence of environmental or demographic stochasticity.} We ran IPM simulations for 2,000 time steps, after an initial 500 iteration burn-in period. This allowed species to reach their stable size distribution. We then calculated yearly per capita growth rates from the simlated time series, and then caluclated the synchrony of species' per capita growth rates over 100 randomly selected contiguous 50 time-step sections.
We ran IBM simulations for 100 time steps, and repeated the simulations 100 times for each simulation experiment. From those, we retained only the simulations in which no species went extinct due to demographic stochasticity. Synchrony of per capita growth rates was calculated over the 100 time steps for each no extinction run within a model experiment. \textcolor{blue}{In the IBM simulations, the strength of demographic stochasticity should weaken as population size increases, meaning that synchrony should be less influenced by demographic stochasticity in large populations compared to small populations. To explore this effect}, we ran simulations \textcolor{blue}{(4)} and (6) on plot sizes of 1, 4, 9, 16, and 25 $\text{m}^2$. All other IBM simulations were run on a 25 $\text{m}^2$ landscape \textcolor{blue}{to most closely match the implicit, large spatial scale of the IPM simulations}.
\textcolor{blue}{Our simulations allow us to quantify the relative importance of environmental responses, species interactions, and demographic stochasticity by comparing the simulated values of community-wide species synchrony.} The simulation experiments also allow us to test our theoretical predictions. First, in the absence of interspecific interactions and demographic stochasticity, populations can only fluctuate in response to the environment. Therefore, we can use results from simulation (5) to estimate the covariance of species' responses to the environment ($\text{cov}(u_{ie}, u_{je})$) and parameterize Equation 4. Parameterizing Equation 3 does not require simulation output because the only parameters are the species' relative abundances. Second, simulations (5) and (6) represent the simulated version of our limiting case theoretical predictions. \textcolor{blue}{This approach to testing theoretical predictions may seem circular, but recall we derived the predictions using strict assumptions about equivalence in species' growth rates, environmental response variances, and demographic variances. Our empirically-based models do not make these assumptions. Comparisons between parameterized predictions and simulated synchrony reveal whether the assumptions we must make to derive analytical predictions hold in the natural community our model represents.}
\subsection{Results} \textcolor{blue}{Observed} synchrony of species' per capita growth rates at our study sites range from 0.36 to 0.89 and synchrony of percent cover ranged from 0.15 to 0.92 (Table 2). Synchrony of per capita growth rates and CV of percent cover are positively correlated (Pearson’s $\rho$ = 0.72, N = 5). For all five communities, species synchrony from IPM and IBM simulations closely approximated observed synchrony (Fig. S1). IBM-simulated synchrony is consistently, but only slightly, lower than IPM-simulated synchrony (Fig. S1), likely due to the desynchronizing effect of demographic stochasticity.
Across the five communities, our limiting case predictions closely matched synchrony from the corresponding simulation experiment (Fig. 2 and Table S1). The correlation between our analytical predictions and simulated synchrony was 0.97 for $\phi_{R,\mathcal{M}D}$ (N = 5) and 0.997 for $\phi{R,\mathcal{M}E}$ (N = 5). The largest difference between predicted and simulated synchrony was 0.05 in New Mexico for $\phi{R,\mathcal{M}_D}$ (Table S1).
Simulation experiments revealed that removing environmental fluctuations has the largest impact on synchrony, leading to a reduction in synchrony of species growth rates in four out of five communities (Fig. 2). Removing environmental fluctuations ("No E.S" simulations) decreased synchrony by 33% in Arizona, 48% in Kansas, 39% in Montana, and 40% in New Mexico. Only in Idaho did removing environmental fluctuations cause an increase in synchrony (Fig. 2), but the effect was small (9% increase; Table S2). Overall, species' temporal random effects in the statistical vital rate models are positively, but not perfectly, correlated (Table S3). \textcolor{blue}{These temporal random effects represent environmental responses, meaning that positively correlated temporal random effects indicate positively correlated environmental responses.}
Species interactions are weak in these communities [Table S4 and @Chu2015], and removing interspecific interactions had little effect on synchrony (Fig. 2; "No Comp." simulations). Removing interspecific interactions caused, at most, a 5% change in synchrony (Fig. 2 and Table S2). Removing demographic stochasticity ("No D.S." simulations) caused synchrony to increase slightly in all communities (Fig. 2), with an average 6% increase over synchrony from IBM simulations on a 25$\text{m}^2$ area.
In IBM simulations, the desynchronizing effect of demographic stochasticity, which increases as population size decreases, modestly counteracted the synchronizing force of the environment, but not enough to lower synchrony to the level observed when only demographic stochasticity is operating (Fig. 3). In the largest, 25 $\text{m}^2$ plots, synchrony was driven by environmental stochasticity (e.g., $\mathcal{M}_E$). At 1 $\text{m}^2$, synchrony reflected \textcolor{blue}{the combined effects of} demographic stochasticity and environmental stochasticity (e.g., \textcolor{blue}{simulated synhcrony fell} between $\mathcal{M}_E$ and $\mathcal{M}_D$). For context, population sizes increased from an average of 17 individuals per community in 1 $\text{m}^2$ IBM simulations to an average of 357 individuals per community in 25 $\text{m}^2$ IBM simulations.
Results for synchrony of percent cover are qualitatively similar, but more variable and less consistent with analytical predictions and observed synchrony (Appendix S1, Figs. S2-S3).
\section{Discussion} Our study produced four main findings that were generally consistent across five natural plant communities: (1) limiting-case predictions from the theoretical model were well-supported by simulations from the empirical models; (2) demographic stochasticity decreased community synchrony, as expected by theory, and its effect was largest in small populations; (3) environmental fluctuations increased community synchrony relative to simulations in constant environments because species-specific responses to the environment were positively, though not perfectly, correlated; and (4) interspecific interactions were weak and therefore had little impact on community synchrony. We also found that analyses based on synchrony of species' percent cover, rather than growth rates, were uninformative (Figs. S2-S3) since the linear approximation required for analytical predictions is a stronger assumption for abundance than growth rates, especially given relatively short time-series (Appendix S1). Thus, our results provide further evidence that it is difficult to decipher mechanisms of species synchrony from abundance time series, as expected by theory [@Loreau2008a]. Observed synchrony of per capita growth rates was positively correlated with the variability of percent cover across our focal communities, which confirms that we are investigating an important process underlying ecosystem stability.
\subsection{Simulations support theoretical predictions} Our theoretical predictions were derived from a simple model of population dynamics and required several restrictive assumptions, raising questions about their relevance to natural communities. For example, the species in our communities do not have equivalent environmental and demographic variances (Figs. S4-S7), as required by our predictions. However, the theoretical predictions closely matched results from simulations of population models fit to long-term data from natural plant communities (Table S1). Strong agreement between our analytical predictions and the simulation results should inspire confidence in the ability of simple models to inform our understanding of species synchrony even in complex natural communities, and allows us to place our simulation results within the context of contemporary theory. \textcolor{blue}{However, whether the theoretical model adequately represents more complex communities remains unknown because our analysis was restricted to dominant species.}
\subsection{Demographic stochasticity decreases synchrony \textcolor{blue}{in small populations}} In large populations, removing demographic stochasticity had no effect on species synchrony (Fig. 2). In small populations, demographic stochasticity partially counteracted the synchronizing effects of environmental fluctuations and interspecific interactions on per capita growth rates, in agreement with theory [@Loreau2008a]. \textcolor{blue}{This is shown in Fig. 3, where IBM simulations with environmental forcing and demographic stochasticity have higher synchrony than simulations with only demographic stochasticity. The differences between simulations are smaller at lower population sizes because as the number of individuals decreases with area, the strength of demographic stochasticity increases, reducing the relative effect of environmental forcing.} Even in small populations \textcolor{blue}{(e.g., 1 $\text{m}^2$ lanscapes)}, however, demographic stochasticity was not strong enough to compensate the synchronizing effects of environmental fluctuations and match the analytical prediction where only demographic stochasticity is operating ($\mathcal{M}_D$ in Fig. 3). These results confirm the theoretical argument of @Loreau2008a that independent fluctuations among interacting species in a non-constant environment should be rare. Only in the Idaho community does synchrony of per capita growth rates approach $\mathcal{M}_D$ in a non-constant environment (Fig. 3). This is most likely due to the strong effect of demographic stochasticity on the shrub Artemisia tripartita since even a 25 $\text{m}^2$ quadrat would only contain a few individuals of that species.
Our analysis of how demographic stochasticity affects synchrony demonstrates that synchrony depends on the observation area. As the observation area increases, population size increases and the desynchronizing effect of demographic stochasticity lessens (Fig. 3). Thus, our results suggest that community-wide species synchrony will increase as the observation area increases, rising from $\mathcal{M}_D$ to $\mathcal{M}_E$. Such a conclusion assumes, however, that species richness remains constant as observation area increases, which is unlikely [@Taylor1961]. Recent theoretical work has begun to explore the linkage between ecosystem stability and spatial scale [@Wang2014; @Wang2016], and our results suggest that including demographic stochasticity in theoretical models of metacommunity dynamics may be important for understanding the role of species synchrony in determining ecosystem stability across spatial scales.
\subsection{Environmental fluctuations drive community synchrony} In large populations where interspecific interactions are weak, synchrony is expected to be driven exclusively by environmental fluctuations (Equation 4). Under such conditions community synchrony should approximately equal the synchrony of species' responses to the environment [@Loreau2008a]. Two lines of evidence lead us to conclude that environmental fluctuations drive species synchrony in our focal plant communities. First, in our simulation experiments, removing interspecific interactions resulted in no discernible change in community-wide species synchrony of per capita growth rates (Fig. 2). Second, removing environmental fluctuations from simulations consistently reduced synchrony (Fig. 2). Our results lead us to conclude that environmental fluctuations, not species interactions, are the primary driver of community-wide species synchrony \textcolor{blue}{among the dominant species} we studied. Given accumulating evidence that niche differences in natural communities are large [reviewed in @Chu2015], and thus species interactions are likely to be weak, our results may be general in natural plant communities.
In the Idaho community, removing environmental fluctuations did not cause a large decrease in synchrony. However, that result appears to be an artifact. Removing environmental variation results in a negative invasion growth rate for A. tripartita. Although we only analyzed IBM runs in which A. tripartita had not yet gone extinct, it was at much lower abundance than in the other simulation runs. When we removed A. tripartita from all simulations, the Idaho results conformed with results from all other sites: removing environmental stochasticity caused a significant reduction in species synchrony (Fig. S8). Our main results for Idaho (Fig. 2), with A. tripartita included, demonstrate how the processes that determine species synchrony interact in complex ways. A. tripartita has a facilitative effect on each grass species across all vital rates, except for a small competitive effect on H. comata's survival probability (Tables S8-S10). At the same time, all the perennial grasses have negative effects on each other for each vital rate (Tables S8-S10). We know synchrony is affected by interspecific competition [@Loreau2008a], but how facilitative effects manifest themselves is unknown. The interaction of facilitation and competition is clearly capable of having a large effect on species synchrony, and future theoretical efforts should aim to include a wider range of species interactions.
\textcolor{blue}{Environmental responses synchronized dynamics relative to a null expectation of independent species interactions (e.g., "No Comp. + No E.S." simulations in Fig. 2), but observed and simulated synchrony was still less than one in all cases (Fig. 2). Synchrony was far from complete because of differences in species' responses to interannual environmental variation. Many studies of ecosystem stability in semiarid grasslands focus on trad-offs among dominant functional types} [@Bai2004; @Sasaki2011]. \textcolor{blue}{Such groupings are based on the idea that ecologically-similar species will have similar responses to environmental fluctuations. At first glance our results may appear to support the grouping of perennial grasses in one functional type because their environmental responses were positively correlated. However, even though environmental responses among the dominant species we studied were similar, they were dissimilar enough to cause synchrony to be less than perfect (Fig. 2). The subtle differences among ecologically-similar dominant species do impact species synchrony and, ultimately, ecosystem stability. Ignoring such differences could mask important dynamics that underpin ecosystem functioning.}
\subsection{Interspecific interactions had little impact on community synchrony} We expected community synchrony of per capita growth rates to decrease when we removed interspecific interactions [@Loreau2008a]. We found that community synchrony was virtually indistinguishable between simulations with and without interspecific interactions (Fig. 2). The lack of an effect of interspecific interactions on synchrony is in contrast to a large body of theoretical work that predicts a strong role for competition in creating compensatory dynamics [@Tilman1988] and a recent empirical analysis [@Gross2014].
Our results do not contradict the idea that competition can lead to compensatory dynamics, but they do highlight the fact that interspecific competition must be relatively strong to influence species synchrony. The communities we analyzed are composed of species with very little niche overlap [@Chu2015] and weak interspecific interactions (Tables S1 and S3-S17). Mechanistic consumer-resource models [@Lehman2000] and phenomenological Lotka-Volterra models [@Lehman2000; @Loreau2013] both confirm that the effect of competition on species synchrony diminishes as niche overlap decreases. In that sense, our results are not surprising: interspecific interactions are weak, so of course removing them does not affect synchrony.
Our conclusion that species interactions have little impact on synchrony only applies to single trophic level communities. Species interactions almost certainly play a strong role in multi-trophic communities where factors such as resource overlap [@Vasseur2007], dispersal [@Gouhier2010], and the strength of top-down control [@Bauer2014] are all likely to affect community synchrony.
\section{Conclusions} Species-specific responses to temporally fluctuating environmental conditions is an important mechanism underlying asynchronous population dynamics and, in turn, ecosystem stability [@Loreau2013]. When we removed environmental variation, we found that synchrony decreased in four out of the five grassland communities we studied (Fig. 2). A tempting conclusion is that our study confirms that compensatory dynamics are rare in natural communities, and that ecologically-similar species will exhibit synchronous dynamics [e.g., @Houlahan2007]. Such a conclusion misses an important subtlety. The perennial grasses we studied do have similar responses to the environment (Table S2), which will tend to synchronize dynamics. However, if community-wide synchrony \textcolor{blue}{among dominant species} is less than perfect, as it is in all our focal communities, some degree of compensatory dynamics must be present [@Loreau2008a]. \textcolor{blue}{Even ecologically-similar species, which are sometimes aggregated into functional groups, have environmental responses that are dissimilar enough to limit synchrony. Subtle differences among dominant species ultimately determine ecosystem stability and should not be ignored.}
\section{Acknowledgements} This work was funded by the National Science Foundation through a Postdoctoral Research Fellowship in Biology to ATT (DBI-1400370) and a CAREER award to PBA (DEB-1054040). CdM and ML were supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41). We thank the original mappers of the permanent quadrats at each site and the digitizers in the Adler lab, without whom this work would not have been possible. Compute, storage, and other resources from the Division of Research Computing in the Office of Research and Graduate Studies at Utah State University are gratefully acknowledged. We also thank Patrick Venail and several anonymous reviewers who provided critical feedback that improved the manuscript.
\setlength{\parindent}{0ex} \newpage{}
\section{Tables}
\singlespacing
Table 1: Site descriptions and focal species. \footnotesize
| Site Name | Biome | Location (Lat, Lon) | Obs. Years | Species | | --------- | ----- | ------------------- | ---------- | ------- | | New Mexico | Chihuahuan Desert | 32.62° N, 106.67° W | 1915-1950 | Bouteloua eriopoda | | | | | | Sporobolus flexuosus | | Arizona | Sonoran Desert | 31°50' N, 110°53' W | 1915-1933 | Bouteloua eriopoda | | | | | | Bouteloua rothrockii | | Kansas | Southern mixed prairie | 38.8° N, 99.3° W | 1932-1972 | Bouteloua curtipendula | | | | | | Bouteloua hirsuta | | | | | | Schizachyrium scoparium | | Montana | Northern mixed prairie | 46°19' N, 105°48' W | 1926-1957 | Bouteloua gracilis | | | | | | Hesperostipa comata | | | | | | Pascopyrum smithii | | | | | | Poa secunda | | Idaho | Sagebrush steppe | 44.2° N, 112.1° W | 1926-1957 | Artemisia tripartita | | | | | | Pseudoroegneria spicata | | | | | | Hesperostipa comata | | | | | | Poa secunda |
\normalsize
\pagebreak{}
library(communitySynchrony) library(plyr) library(reshape2) library(synchrony) library(ggplot2) library(xtable) library(gridExtra) library(RColorBrewer)
site <- "Kansas" spp_list <- c("BOCU","BOHI","SCSC") num_spp <- length(spp_list) ks_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now ks_data <- rbind(ks_data, spp_data) } #end species looping for raw data ks_data <- ks_data[2:nrow(ks_data),] #remove first NA row tmp1<-which(ks_data$quad_data=="q25" & (ks_data$year<35 | ks_data$year>62)) tmp2<-which(ks_data$quad_data=="q27") tmp3<-which(ks_data$quad=="q28") tmp4<-which(ks_data$quad=="q30") tmp5<-which(ks_data$quad=="q31" & (ks_data$year<35 | ks_data$year>39)) tmp6<-which(ks_data$quad=="q32" & (ks_data$year<35 | ks_data$year>41)) tmp<-c(tmp1,tmp2,tmp3,tmp4,tmp5,tmp6) ks_data<-ks_data[-tmp,] # exclude the records later than 1968, to keep the same random year effect... ks_data<-subset(ks_data,year<68) site <- "Idaho" spp_list <- sort(c("PSSP","HECO","POSE","ARTR")) num_spp <- length(spp_list) id_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now id_data <- rbind(id_data, spp_data) } #end species looping for raw data id_data <- id_data[2:nrow(id_data),] #remove first NA row # id_data <- subset(id_data, species!="ARTR") #take out the shrub site <- "Montana" spp_list <- sort(c("BOGR","HECO","PASM","POSE")) num_spp <- length(spp_list) mt_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now mt_data <- rbind(mt_data, spp_data) } #end species looping for raw data mt_data <- mt_data[2:nrow(mt_data),] #remove first NA row site <- "NewMexico" spp_list <- sort(c("BOER","SPFL")) num_spp <- length(spp_list) nm_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now nm_data <- rbind(nm_data, spp_data) } #end species looping for raw data nm_data <- nm_data[2:nrow(nm_data),] #remove first NA row site <- "Arizona" spp_list <- sort(c("BOER","BORO")) num_spp <- length(spp_list) az_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now az_data <- rbind(az_data, spp_data) } #end species looping for raw data az_data <- az_data[2:nrow(az_data),] #remove first NA row
out_ks <- get_comm_synchrony(ts_data = ks_data) out_id <- get_comm_synchrony(ts_data = id_data) out_mt <- get_comm_synchrony(ts_data = mt_data) out_nm <- get_comm_synchrony(ts_data = nm_data) out_az <- get_comm_synchrony(ts_data = az_data)
get_table_metrics <- function(output){ tmp <- output exp_tmp <- round(tmp$pgr_expected_synch_ind_flucts, 2) obs_tmp <- round(as.numeric(tmp$pgr_synchrony[1]), 2) demo_tmp <- obs_tmp - exp_tmp exp_cover <- round(as.numeric(tmp$cover_expected_synch_ind_flucts), 2) obs_cover <- round(as.numeric(tmp$abund_synchrony[1]), 2) demo_cover <- obs_cover - exp_cover plant_size <- tmp$avg_size cover_var <- tmp$variability return(c(obs_tmp, exp_tmp, demo_tmp, obs_cover, exp_cover, demo_cover, cover_var)) } table_data <- as.data.frame(rbind(get_table_metrics(out_az), get_table_metrics(out_id), get_table_metrics(out_ks), get_table_metrics(out_mt), get_table_metrics(out_nm))) table_data$Site <- c("Arizona", "Idaho", "Kansas", "Montana", "New Mexico") alltable_data <- table_data colnames(table_data) <- c("obs_gr_synch", "demo_only_exp_grsynch", "diff_noneed", "obs_cover_synch", "demo_only_exp_coversynch", "diffcover_noneed", "cover_variance", "site") saveRDS(table_data, "../results/calculate_site_metrics.RDS") table_data <- alltable_data[,c("Site", "V1", "V4", "V7")] table_data$spp <- c(length(unique(az_data$species)), length(unique(id_data$species)), length(unique(ks_data$species)), length(unique(mt_data$species)), length(unique(nm_data$species))) colnames(table_data) <- c("Site", "$\\phi_{R}$", "$\\phi_{C}$", "CV of Total Cover", "Species richness") myorder <- c(5,1,3,4,2) xtable_data <- table_data[myorder,] # colnames(table_data) <- c("Site", # "$\\phi_{r}$", # "$\\text{E}(\\phi_{r})$", # "$\\text{E}(\\sigma_{r,c} + \\sigma_{r,e})$", # "$\\phi_{C}$", # "$\\text{E}(\\phi_{C})$", # "$\\text{E}(\\sigma_{C,c} + \\sigma_{C,e})$", # "CV of Total Cover", # "Species richness") # table_data[,3] <- round(1/table_data[,ncol(table_data)],2) # reset expected sybchrony to 1/S # table_data[,4] <- table_data[,2] - table_data[,3] # calculate obs - exp synchrony synch_cap <- "Observed synchrony among species' per capita growth rates ($\\phi_{R}$), observed synchrony among species' percent cover ($\\phi_{C}$), the coefficient of variation of total community cover, and species richness for each community. Species richness values reflect the number of species analyzed from the community, not the actual richness." print(xtable(xtable_data, caption = synch_cap), caption.placement="top", include.rownames = F, sanitize.colnames.function = identity, comment=FALSE, size="normalsize")
sink(file = "../results/correlation_Rsynch_cv.txt") cat("correlation between synchrony of growth rates and cover CV is:") cor(table_data$`CV of Total Cover`, table_data[,2]) Sys.Date() sink()
\pagebreak{}
\section{Figures}
\begin{figure}[!ht] \centering \includegraphics[width=6in]{./components/formatted_figures/synchrony_flowchart_bw.png} \caption{Diagram of our coupled theoretical-empirical approach. We followed this workflow for each of our five focal communities.} \end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/formatted_figure1.png}
\caption{Community-wide species synchrony of per capita growth rates from model simulation experiments. Synchrony of species' growth rates for each study area are from simulation experiments with demographic stochasticity, environmental stochasticity, and interspecific interactions present (All Drivers''), demographic stochasticity removed (No D.S.''), environmental stochasticity removed (No E.S.''), interspecific interactions removed (No Comp.''), interspecific interactions and demographic stochasticity removed (No Comp. + No D.S.''), and interspecific interactions and environmental stochasticity removed (No Comp. + No E.S.''). Abbreviations within the bars for the New Mexico site indicate whether the IBM or IPM was used for a particular simulation. Error bars represent the 2.5\% and 97.5\% quantiles from model simulations. All IBM simulations shown in this figure were run on a 25 $\text{m}^2$ virtual landscape. Points show observed and predicted synchrony aligned with the model simulation that corresponds with each observation or analytical predicion.}
\end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/formatted_figure2.png}
\caption{Synchrony of species' growth rates for each study area from IBM simulations across different landscape sizes when only demographic stochasticity is present (No Comp. + No E.S. (D.S. Only)'') and when environmental stochasticity is also present (No Comp. (D.S. + E.S.)''). The horizontal lines show the analytical predictions $\mathcal{M}_D$ (dashed line) and $\mathcal{M}_E$ (dotted line). The strength of demographic stochasticity decreases as landscape size increases because population sizes also increase. Theoretically, No Comp. + No E.S. (D.S. Only)'' simulations should remain constant across landscape size, whereasNo Comp. (D.S. + E.S.)'' simulations should shift from the $\mathcal{M}_D$ prediction to the $\mathcal{M}_E$ prediction as landscape size, and thus population size, increases, but only if demographic stochasticity it strong enough to counteract environmental forcing. Error bars represent the 2.5\% and 97.5\% quantiles from model simulations.}
\end{figure}
\pagebreak{}
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.