Planktonic functional diversity changes in synchrony with lake ecosystem state

Abstract Managing ecosystems to effectively preserve function and services requires reliable tools that can infer changes in the stability and dynamics of a system. Conceptually, functional diversity (FD) appears as a sensitive and viable monitoring metric stemming from suggestions that FD is a universally important measure of biodiversity and has a mechanistic influence on ecological processes. It is however unclear whether changes in FD consistently occur prior to state responses or vice versa, with no current work on the temporal relationship between FD and state to support a transition towards trait‐based indicators. There is consequently a knowledge gap regarding when functioning changes relative to biodiversity change and where FD change falls in that sequence. We therefore examine the lagged relationship between planktonic FD and abundance‐based metrics of system state (e.g. biomass) across five highly monitored lake communities using both correlation and cutting edge non‐linear empirical dynamic modelling approaches. Overall, phytoplankton and zooplankton FD display synchrony with lake state but each lake is idiosyncratic in the strength of relationship. It is therefore unlikely that changes in plankton FD are identifiable before changes in more easily collected abundance metrics. These results highlight the power of empirical dynamic modelling in disentangling time lagged relationships in complex multivariate ecosystems, but suggest that FD cannot be generically viable as an early indicator. Individual lakes therefore require consideration of their specific context and any interpretation of FD across systems requires caution. However, FD still retains value as an alternative state measure or a trait representation of biodiversity when considered at the system level.


System state metrics
Five common measures of ecosystem state were used to provide a representation of the development of each lake's plankton community through time. The five metrics used were community composition, total plankton density, Fisher information (FI), the multivariate index of variability (MVI), and trophic ratio. Total plankton density is simply the sum of plankton densities at each time point (Kraemer et al. 2017) while trophic ratio is the division of total zooplankton density by total phytoplankton density (Jeppesen et al. 2011). Both of these metrics have been used previously to represent overall community functioning.
Fisher information (FI) is derived from information theory and attempts to quantify the amount of information that observable data can provide on an unmeasured parameter (Fisher & Russell 1922). It consequently has been co-opted to act as an indicator of system dynamics and stability, where decreasing FI indicates decreasing system stability (Cabezas et al. 2010).
Here, we used a simplified discrete time equation of Fisher and Russel (1922)'s mathematic proof following Karunanithi et al. (2008): where qi 2 is the amplitude of the probability of observing states of the system at time window i and m is the number of possible 'states'. Possible states are defining by comparing the difference between temporally adjacent data windows to a reference 'uncertainty', typically suggested to be the standard deviation of each variable (in our case each plankton species' density) across the entire time series multiplied by 2 (Karunanithi et al. 2008). If the absolute difference in density is less than the reference deviation for all variables (i.e. species densities), then the windows are binned in to the same 'state'. In this study, we used a rolling window of five time points to ensure maximum coverage. While Fisher information has long been used in statistics, it has recently been suggested as an appropriate stability measure for ecosystems specifically identifying regime shifts (Ahmad et al. 2016;Konig et al. 2019

Fuzzily-coded trait dissimilarity
Fuzzy coding is an approach that has been successfully applied to multiple freshwater macroinvertebrate communities (Brown et al. 2018;Múrria et al. 2020), and has been suggested as an appropriate tool to circumvent some of the challenges non-species level data can bring to trait-based approaches in plankton (Martini et al. 2021). In practice, a fuzzy trait is sub-categorised and 'affinity' is the proportion of the taxon's species which fall within each of the fuzzily coded sub-categories (Gower 1971;de Bello et al. 2021). During this study, each of the lakes have sufficient monitoring documentation that genus level data identifies the species that have had their abundances pooled to form a taxon count. Therefore, we were able to code only species known to contribute to the count rather than coding all possible taxon members present in the trait databases. A fuzzily coded trait matrix was consequently uniquely constructed for each lake and plankton guild (phytoplankton vs zooplankton), and from it a dissimilarity matrix based on a modified Gower index was derived using the 'gawdis' package in R (R Core Team 2020). This dissimilarity matrix was then Cailliez transformed (Cailliez 1983) to improve suitability for Euclidean-based analysis during construction of the 'trait space'. Our choice of a dissimilarity-based functional diversity methodology stems from our aim to quantify changes in trait space size through time (Mammola et al. 2021).

Functional diversity metrics
Functional diversity is usually estimated using three metrics (functional richness -FRic, functional evenness -FEve, and functional dispersion -FDis) which together attempt to quantify different aspects of the trait space available to a community at any one point in time (Laliberté & Legendre 2010). To achieve this, principal coordinate analysis is performed on the trait dissimilarity matrix to estimate species coordinates within the resulting multidimensional space (Magneville et al. 2022). From these coordinates, the three functional diversity metrics can be calculated. Functional richness characterises the area encompassed by the community as measured by a convex hull, functional evenness is the length of the minimum spanning tree between species positions, and functional dispersion quantifies the average distance to the point of community 'gravity' (the centroid of community values weighted towards abundant species). Figure S2 gives a visual representation of these metrics.

Convergent cross mapping
Convergent cross mapping invokes Takens' embedding theorem (Takens 1981) which suggests that an underlying latent system/attractor manifold can be reconstructed from one or Causality is identified by comparing the quality of prediction of one time series from the reconstructed system built upon the time-delay embedding of another (Runge et al. 2019); if X can be predicted using the reconstruction from Y, then X had a causal effect on Y. The reciprocal relationship is then also tested to establish whether bi-directional causation is present.
The appropriate embedding dimension for CCM was selected using nearest neighbour forecasting, specifically simplex projection (Sugihara & May 1990). Simplex projections generate forecasts of the manifold reconstruction and identifies the embedding dimension with the highest prediction skill (represented as the correlation between observed and predicted values). In this analysis, we explored dimensions from 0 to 10 with an increment of 1, using * = 1 (a month) for embedding. Using the identified embedding dimension, CCM analysis was performed for each time-to-prediction lag/delay with a library size set to the maximum possible to improve prediction error and the likelihood of convergence (Sugihara et al. 2012).  Figure S2. Diagrammatic representation of the three functional diversity metrics based upon trait dissimilarity. Points denote a species' trait value/coordinate within the trait space, with the size of the point proportional to the abundance of that species at time t. FRic is functional richness, FEve is functional evenness, and FDis is functional dispersion. Figure S3. Density of the permuted instantaneous (Lag0) cross-correlation coefficients following 10,000 simulations between functional diversity and each of the system state metrics. Points represent the observed correlation coefficient and stars represent the 'significant' relationships (where the observed value transgresses the 2.5-97.5% confidence intervals). Permuted CCF t0 between detrended FD and system state Figure S4. Density of permuted cross-correlation coefficients following 10,000 simulations between functional diversity and each of the system state metrics (LagX). Points represent the absolute strongest observed correlation coefficient and stars represent the 'significant' relationships (where the observed value transgresses the 2.5-97.5% confidence intervals). The reported number is the optimal lag (in months) that the observed value was identified. Positive lags indicate functional diversity lagging state, whereas negative lags indicate functional diversity leading state.  Figure S5. Density plots of the optimal lags identified by cross-correlations pooled across all functional diversity:system state combinations (blue area) or subset to significant combinations only (yellow area). Dashed, vertical lines indicate 10 th and 90 th quartiles. Figure S6. Density of instantaneous (Lag0) permuted cross mapping skills where the system state metrics map functional diversity following 1000 simulations. Points represent the observed cross map skill and stars represent the 'significant' relationships (where the observed value exceeds the 95% quantile). These points therefore represent the predicted strength of 'forward' causality where functional diversity 'causes' system state.  Figure S7. Density of instantaneous (Lag0) permuted cross mapping skills where functional diversity maps system state following 1000 simulations. Points represent the observed cross map skill and stars represent the 'significant' relationships (where the observed value exceeds the 95% quantile). These points therefore represent the predicted strength of 'reverse' causality where system state 'causes' functional diversity. Permuted cross skill of FD mapping system state (reverse) Figure S8. Density of permuted cross mapping skills (LagX) where the system state metrics map functional diversity following 10,000 simulations. Points represent the observed cross map skill and stars represent the 'significant' relationships (where the observed value exceeds the 95% quantile). These points therefore represent the predicted strength of 'forward' causality where functional diversity 'causes' system state. Optimal lags (in months) are indicated by the numbers above each distribution. Negative lags indicate functional diversity leading state whereas positive lags indicate generalised synchrony.  Figure S10. Density of permuted cross mapping skills (LagX) where functional diversity maps system state following 10,000 simulations. Points represent the observed cross map skill and stars represent the 'significant' relationships (where the observed value exceeds the 95% quantile). These points therefore represent the predicted strength of 'reverse' causality where system state 'causes' functional diversity. Optimal lags (in months) are indicated by the numbers above each distribution. Negative lags indicate functional diversity leading state whereas positive lags indicate generalised synchrony.          Table S7. Summary statistics for each functional diversity-system state cross mapping across lags (LagX). This is the reverse relationship (functional diversity maps system state) where a significant relationship suggests diversity is caused by state. Variation is reported as standard errors.  Table S8. Causality comparison for cross mappings between functional diversity and system state, across lags. The proportion of lakes with forward, reverse, bi-directional and zero causality are reported. Variation is reported as standard errors.