Forecasting in the face of ecological complexity: Number and strength of species interactions determine forecast skill in ecological communities.

The potential for forecasting the dynamics of ecological systems is currently unclear, with contrasting opinions regarding its feasibility due to ecological complexity. To investigate forecast skill within and across systems, we monitored a microbial system exposed to either constant or fluctuating temperatures in a 5-month-long laboratory experiment. We tested how forecasting of species abundances depends on the number and strength of interactions and on model size (number of predictors). We also tested how greater system complexity (i.e. the fluctuating temperatures) impacted these relations. We found that the more interactions a species had, the weaker these interactions were and the better its abundance was predicted. Forecast skill increased with model size. Greater system complexity decreased forecast skill for three out of eight species. These insights into how abundance prediction depends on the connectedness of the species within the system and on overall system complexity could improve species forecasting and monitoring.

INTRODUCTION

Over the last decades, it has become increasingly important to proficiently predict the consequences of climate change and biodiversity loss (e.g. Dietze, 2017; Godfray & May, 2014). Ecological forecasting, formally defined as the prediction of natural capital and ecosystems states and services, has advanced to be an imperative scientific and applied discipline (Clark et al., 2001; Dietze et al., 2018; Houlahan et al., 2017). Examples of its applications include predicting ecotoxicological effects on community responses (e.g. Clements & Rohr, 2009), forecasting the successes of species invasions (e.g. Romanuk et al., 2009) and predicting how communities respond to climate change (e.g. Gauzere et al., 2018; Hattab et al., 2016; McCarthy et al., 2018). However, in the context of the complexity of real‐world systems, skilful ecological forecasting remains a major challenge to the point that its feasibility has been questioned (Beckage et al., 2011; Hayes & Barry, 2008; Planque, 2016).

Generally, an ecological network or system is more complex the more biotic and abiotic variables are part of it (Bradbury & Vehrencamp, 2014; Mitchell, 2009). The number of possible indirect interaction pathways between variables rapidly increases with increasing network size (Borrett & Patten, 2003), and this is believed to hinder skilful predictions (Wootton, 2002; Yodzis, 1988). In fact, some studies have found that prediction skill deteriorated with increasing system complexity (e.g. Doak et al., 2008; Jonsson et al., 2018; Novak et al., 2011) and that species interactions can reduce community predictability (Thompson et al., 2021), resulting in the view that ecology is unpredictable due to its complexity (Beckage et al., 2011). Yet, some recent results suggest the opposite: complexity can increase rather than decrease prediction skill (Barbier et al., 2021; Iles & Novak, 2016; Mougi, 2017). For instance, it was found that the total abundance (the sum of all species abundances) was more predictable when the system consisted of more species (Dornelas et al., 2011), while another study showed that the prediction of interaction strengths improved with increasing food web size (Berlow et al., 2009). With evidence pointing in both directions, it remains unclear whether there is a general relation between system complexity and forecast skill or whether each result is specific to the system and to the quantity forecasted.

Within a system of a given complexity, commonly only a few strong species interactions are present with most interactions being weak (e.g. Bascompte et al., 2005; Berlow et al., 2004; Paine, 1992; Wootton & Emmerson, 2005), though still being important for system stability (e.g. Kadoya et al., 2018; McCann et al., 1998; O'Gorman & Emmerson, 2009). Moreover, there is some evidence indicating that generalists (i.e. species with many possible interacting partners) mostly have weak interactions, while specialist species (i.e. fewer possible interacting partners) show stronger interactions (Wootton & Stouffer, 2016). As in the case of increasing system complexity, the more interactions a species has in a network the more indirect pathways exist that can influence its abundance. Hence, if prediction skill generally decreases with increasing complexity, we might hypothesise that the prediction skill for a given species will also decrease the more interactions that species has within a network, and consequently also the weaker these interactions are on average.

In this study, we investigated whether the prediction of species abundances depends on how many interactions the species have, a hitherto untested question that could help explain why some species can be forecast better than others (Harris, 1994). We also tested whether the forecast skill of species abundances depends on how many system variables are considered in the forecast.

Using a laboratory‐based aquatic microbial community as our study system, we carried out a 5‐month‐long experiment. The community consisted of algae, bacteria, ciliates, flagellates and rotifer species. These species are characterised by short generation times (e.g. Altermatt et al., 2015), which renders them convenient study organisms for our experiment and questions. We exposed half of the community replicates to either a constant temperature or a fluctuating temperature setting. The fluctuating temperatures added a layer of complexity to the system, by potentially affecting the species and their interactions in both direct and indirect ways. We chose temperature fluctuations due to their relevance and relative ease of manipulation.

We forecasted species abundances and estimated the number and the strength of species interactions using the nonparametric time‐series analysis framework empirical dynamic modelling (EDM, Ye et al., 2015). We built forecasts that included increasingly more variables as predictors. We hypothesised that the more a species is isolated (i.e. fewer and weaker interactions), the better its abundance can be predicted as it is less dependent on the system state. We expected that the fluctuating temperatures would decrease forecast skill, unless they are a strong enough driver of system dynamics to outweigh the effect of increased system complexity. Further, we hypothesised that up to a threshold forecasting improves when more system variables are included in the prediction, but that fewer variables are necessary to achieve the highest or close to the highest forecast skill for more isolated species. We also tested whether a variable that interacts strongly with a focal species is a good predictor variable of the abundance of said species.

MATERIAL AND METHODS

Experiment: design, setup and sampling

We carried out a laboratory‐based experiment to record the dynamics of microbial communities (i.e. microcosms) at constant (17.3°C) and at fluctuating temperatures over a period of 154 days. We used three different fluctuating temperature time series. One was identical to the temperature of a local small stream (Furtbach ZH, Switzerland); the other times series had the same mean, variance and autocorrelation structure as the natural time series, but differed in the temperature order (Cohen et al., 1999; Petchey, 2000). See the Section S2.1 for further information regarding the temperature time series.

The tri‐trophic microcosms were semi‐naturalistic with respect to the potential co‐occurrences of the species and the functional groups present (Table S5, Figure S2). The first trophic level of the community consisted of three bacteria (Serratia fonticola, Brevibacillus brevis and Bacillus subtilis), an autotroph alga (Chlamydomonas reinhardtii), a mixotroph alga (Euglena gracilis) and a mixotroph ciliate (Euplotes daidaleos). The latter two are mixotrophic species and their trophic level is between the first and the second level (Ward & Follows, 2016). The second level contained three bacterivore ciliate species (Colpidium striatum, Dexiostoma campylum and Spirostomum sp.), one omnivore ciliate species (Paramecium caudatum) and one omnivore rotifer species (Rotifer sp.), while one ciliate predator species (Didinium nasutum) made up the top level. Small non‐identified flagellate species present in the species stock cultures were also part of the microbial communities and we classified them into the three groups ‘small and white flagellates’, ‘green and white flagellates’ and ‘big and white flagellates’. Prior to the experiment, we kept the ciliate and algae species in stock culture jars at 20°C containing organic protozoan pellet medium (Carolina Biological Supply Company, Burlington NC; concentration of 0.55 gl−1, Altermatt et al., 2015). For the heterotrophic and mixotrophic species, we bacterised the medium with the same bacteria species. We fed D. nasutum with P. caudatum ad libitum and freshly established all stock cultures 2 weeks prior to the experiment.

We set up the microcosms in 2 L screw‐capped glass bottles filled with 250 ml of the non‐bacterised medium containing C. reinhardtii at 50 cells/ml, 750 ml of the bacterised medium, a magnetic stirrer and 20 wheat seeds for slow and continuous nutrient release. We added the remaining species (except D. nasutum and Spirostomum) at a density of 0.1 cells/ml. As Spirostomum sp. only reached low abundances in the stock cultures we inoculated it at a density of 0.005 cells/ml. We added the predator D. nasutum at a density of 0.02 cells/ml 9 days after the start of the experiment. We reintroduced all species (except bacteria and C. reinhardtii) at low densities (<0.01 cells/ml) once a week. This rate was high enough that extinct species could potentially re‐establish in the long term and low enough to not influence population dynamics in the short term.

We kept the experimental bottles in temperature‐controlled incubators with a 14/10 h light–dark cycle. We distributed 18 replicates across six incubators. We set three of the incubators to the constant temperature and assigned one incubator to each of the three fluctuating temperature time series. Thus, nine replicates were in the constant and nine in the fluctuating temperature environment.

We sampled the microcosms three times per week (Mondays, Wednesdays and Fridays) for 22 weeks (66 data points per microcosm). We measured dissolved oxygen concentration using a non‐invasive oxygen recorder (Precision Sensing GmbH, Germany) with oxygen sensing optodes attached to the inside of the bottles. Before sampling, we homogenised the microcosms on a magnetic plate. We sampled 65 ml per replicate and added the same amount of bacterised medium to them afterwards. We measured the abundances of the grouped bacteria species and of the small, intermediate and large species by, respectively, using flow cytometry, FlowCAM imaging, video microscopy and manual counts (Table S5). Video microscopy involved the R‐package bemovi (Pennekamp et al., 2015). For the video and the FlowCAM data, we used automated species classifications. For more details regarding measurements and classifications, see supplementary Section S2.2.

Processing of recorded time series

Preceding analyses, we processed the recorded time series as follows (Benincà et al., 2008, see Section S2.1): We first interpolated the time series using a cubic hermite spline to obtain equally distanced time points (time step of 2.3 days). To flatten sharp changes in abundances, we carried out a fourth‐root power transformation. We then regressed the time series against time and henceforth used the resulting residuals, which are trendless, after we standardised them.

Throughout the experiment, Spirostomum sp. remained practically extinct (Figure S3m). As this species was effectively not part of the microbial community, we did not consider it in the subsequent analyses. Further, the predator D. nasutum did not show stable abundance. We used this species only as a predictor variable. Accordingly, the forecasted target species were C. reinhardtii, E. gracilis, E. daidaleos, C. striatum, D. campylum, P. caudatum, Rotifer sp. and the three bacteria species considered as one group (for simplicity henceforth referred to as a species).

Forecasting of species abundances

We forecasted the abundances of species using the empirical dynamic modelling (EDM) technique multiview embedding (Ye et al., 2015; Ye & Sugihara, 2016), as species dynamics are often nonlinear (Blonder et al., 2017; Clark & Luis, 2020). In EDM, forecasting is based on the assumption that similar system states will lead to subsequent system states that are again similar. In this method, state variables are used as predictors in both a non‐lagged and a lagged fashion, following Takens' theorem that the time series of a variable contains information of interacting variables (Takens, 1981). The lagged and non‐lagged time series re‐construct the attractor manifold and the number of time series used for this is the embedding dimension E.

Multiview embedding (Ye & Sugihara, 2016, Section S3) is an extension of this method in which for a fixed embedding dimension all possible combinations (called ‘views’) of the predictor time series are constructed, which are then ranked by in‐sample forecast skill and the best k views are used for an average out‐of‐sample forecast. We used an embedding dimension of E = 3 (sensitivity analysis in Section S8.3.1) and a maximum lag of l = 3 (i.e. we lagged the predictors by zero days, 2.3 days and 4.7 days, with species generation time ranging from hours to days, see Leary & Petchey, 2009; Altermatt et al., 2015).

For each species, we repeatedly forecasted its abundance using increasingly more predictors (i.e. we increased the forecast model size). Excluding temperature, there were 13 possible non‐lagged predictors for each target species (the eight target species, D. nasutum, the three flagellate groups and the dissolved oxygen). As the number of predictors, we used n = {1,2,3,4,6,8,10,13}. For each value of n we calculated the number of possible non‐lagged predictor combinations . Out of these combinations, we randomly selected 200 if the number of combinations exceeded this value. We used the function Multiview() (R‐package rEDM, Park et al., 2021) which adds the lags to the predictor variables. We used the first 44 time points of the time series as the in‐sample data and the last 22 points as the out‐of‐sample time points to be predicted (one‐step ahead forecasts). For each predictor combination, we evaluated up to 25 values for k logarithmically spaced between 1 and 100. We previously confirmed that values above 100 do not improve forecasts (Section S3.1). We then repeated all forecasts with the temperature added to the same predictor combinations and also for when temperature was the sole predictor. In total, we fitted more than 7.4 million multiview embedding models (Table S6). As a measure of forecast error, we calculated the root mean square error (RMSE) of each fitted model based on the out‐of‐sample data. Because we standardised the time series, an RMSE below one indicates that the used model predicts the abundance of a species better than its mean abundance does.

Estimation of the number and strength of interactions

For each target species T, we determined which and how many state variables influenced its abundance with convergent cross‐mapping (CCM, Sugihara et al., 2012), following the recommendations of Deyle et al. (2016, supplementary Section S4). CCM is a test of causation that reveals whether there is a causal link between variables. We defined the variables that showed a significant effect on a target species in a replicate as its interactors and their number as the number of interactions N .

We then estimated the interaction strength time series of the species that were causally linked using S‐map EDM (Deyle et al., 2016). The estimated pairwise species interaction strength time series are S (t) = ∂T(t + τ)/∂I(t), where t is a time point and τ = 2.3 is the smallest time step, ∂ indicates the partial derivative and T(t) and I(t) are the transformed abundance time series of the target and the interactor state variables respectively. In S‐map EDM, at each time point t the interaction Jacobian (i.e. , where and are the different target and interactor state variables) is calculated (Chang et al., 2021). This calculation is done by including information of when the system was in a similar state at other times using locally weighted multivariate linear regressions. The parameter θ determines how nearby system states are weighted in the regression. We used an intermediate value (θ = 5, sensitivity analysis in Section S8.2.1). The interactor variables were limited to those that influenced the target in a given replicate (based on the CCM analysis described above). We used the same eight target species as for the abundance forecasting.

Forecast error analyses

Relation between the number of interactions, mean interaction strength and forecast error

We calculated the mean interaction strength μ of the target species T with the N state variables it interacted with as: In Equation 1, I is the set consisting of the interactors that affected the target T, S is the interaction strength time series between target T and interactor I ∈ I , |S (t)| is its absolute value at time point t, N is the number of interactions and L is the number of time points in the time series. We then computed the sum of interaction strengths ∑ by multiplying Equation 1 with the number of interactions N : ∑  = N μ .

We investigated the relations between the three explanatory variables N , μ and ∑ and the forecast error (RMSE) of species abundances. The RMSE values were based on the forecast models in which all variables were used as predictors. We fitted three separate linear mixed models with the RMSE as the response variable and one of the three explanatory variables as the regressor. We fitted a fourth mixed model between μ (response variable) and N . In each model, we included the temperature regime (i.e. constant or fluctuating) and its interaction with the other explanatory variable.

Forecast error as a function of the number of predictors and the number of interactors

We investigated the relation between the median forecast error and the number of predictors and the temperature regime (constant or fluctuating) with a linear mixed model conjointly for the eight different target species. We used the median RMSE as the response variable, while the log10‐transformed number of used predictors, the temperature regime, a binary variable indicating whether temperature was used as a predictor and the target species were the explanatory variables, alongside their pairwise interactions. We included bottle ID nested in target species as a random intercept.

Using the same settings as before, we then forecasted species abundances again using as predictors only variables that influenced the target species (based on the CCM analysis). Among all forecast models, we selected the ones that predicted the target species the best. For this, from the models that yielded an RMSE within 1% of the lowest achieved RMSE (for a given species and replicate), we selected the models with the fewest predictor variables (i.e. the smallest models). We then used the number of predictors in the best forecast models as the response variable in a mixed model that included the number of interactions, the temperature regime and their interaction as explanatory variables.

Interactor strength versus predictor importance

We investigated whether stronger interactors are also better predictors. In each replicate, we calculated the mean interaction strength of each evaluated target‐interactor species pair as the mean of the absolute values of their interaction strength time series. We log10‐transformed this variable and used it in a mixed model in a three‐way interaction with the temperature regime and the target species. The response variable was the RMSE of the forecast model in which the interacting species was the only predictor.

We fitted all linear mixed models using the function lmer (R‐package lme4, Bates et al., 2015). We included bottle ID as a random intercept in all models, if not specified otherwise.

RESULTS

Relation between the number of interactions, mean interaction strength and forecast error

Both the number and the mean strength of interactions of a target species had a significant effect on the forecast error of species abundances (Figure 1, Table S1). Neither the temperature to which the microcosms were exposed nor its interaction with the other explanatory variables had significant effects on any of the following results (Table S1). Forecast error decreased (i.e. forecast skill increased) the more interactions a species had (Figure 1a). Quantitatively, with every unit increase in number of interactions the forecast error decreased by 5.2% (constant temperature) and 3.9% (fluctuating temperature) with respect to its maximum (t 140 = −6.04, p < 0.001). Meanwhile, the stronger a species interacted on average, the worse it was predicted (Figure 1b): the forecast error increased by 6.5% (constant) and 4.8% (fluctuating) towards its maximum for every 0.1 increase in mean interaction strength (t 127 = 6.45, p‐value < 0.001). The number of interactions and the mean interaction strength of a species were negatively correlated, with the mean interaction strength decreasing by 0.053 (constant) and 0.055 (fluctuating) for every unit increase in number of interactions (t 140 = −9.66, p < 0.001, Figure 1c), indicating that the more interactions a species had the weaker these were. These two measures consequently share considerable explanatory power regarding the RMSE (adjusted R 2: 0.31 for the single‐regressor models, 0.34 for the model that included both quantities). The sum of interaction strengths (the product of these two quantities) was thus unrelated to forecast error (t 140 = 1.40, p = 0.163, Figure 1d). Each target species had a comparable number of interactions across replicates (Figure S8).

FIGURE 1

Species with more but weaker interactions are predicted better. (a, b) Forecast error (RMSE), respectively, as a function of the number of interactions (a) and the mean interaction strength (b). (c) The relation between the number of interactions and the mean interaction strength. (d) Forecast error (RMSE) as a function of the sum of interaction strengths. The black line represents the fit of the respective regression models, and the shaded regions indicate the corresponding 95% confidence intervals. In each sub‐figure, the left and the right panel show the data for the constant and the fluctuating temperature case respectively. The respective intercepts and slopes do not significantly differ between the two temperature settings.

The described patterns persisted in a sensitivity analysis for the parameters θ and E (Figure S13 and Figure S15) and across several robustness analyses reported in detail in the Section S8. For instance, we repeated the forecasting by using recurrent neural networks and ARIMA (Figure S22), and estimated interaction strengths with regularised S‐map (Cenci et al., 2019; Figure S14) and with a multivariate auto‐regressive system state model (Holmes et al., 2012; Figure S23). These analyses show that our results are robust to the forecasting method used as well as how interaction strengths are quantified. They also provide evidence that the negative correlation between the number and the mean strength of interactions is not a statistical artefact of the analysis methods.

Forecast error as a function of the number of predictors and temperature

In general, the median forecast error of species abundances decreased the more variables were used as predictors (F = 1847.81, p < 0.001, Figure 2, Table S2): the respective slopes ranged from −0.579 to −0.055 across target species (F = 175.86, p < 0.001), with 0.022 subtracted to these slopes in the case of fluctuating temperatures (F = 5.76, p = 0.016). Fluctuating temperatures significantly affected forecast errors (F = 25.09, p < 0.001), but this was not the case for the forecasting of all target species as the difference in forecast error ranged from −0.034 to 0.373 across them (F = 4.79, p < 0.001, Tables S2 and S7). Specifically, the fluctuating temperatures increased the forecasting error of C. reinhardtii, Rotifer sp. and P. caudatum (respective differences in estimated marginal mean forecast errors of 0.373, 0.196 and 0.183, with respective p‐values of <0.001, 0.003 and 0.005, see Section S6). The inclusion of temperature as a predictor decreased the forecast error (F = 26.62, p < 0.001), with the change ranging from −0.079 to −0.005 across targets (F = 7.53, p < 0.001), with no significant difference between the two temperature regimes (F = 0.31, p = 0.580). The value of using temperature as a predictor decreased the more other predictors were used in the forecasting, with temperature reducing forecast errors by 0.048 less for every 10‐fold increase in the number of predictors (F = 27.53, p < 0.001).

FIGURE 2

Median forecast error of species abundance decreases the more state variables are used as predictors. Median forecast error as a function of the number of predictors used in the forecast model, of the predicted target species (sub‐figures), of the temperature regime and of whether temperature was used as a predictor (colour). In each sub‐figure, the left panel shows the data as a scatter plot with the points connected with solid lines if they come from the same replicate, while in the right panel the solid lines represent the fit of the fitted regression model, and the shaded regions indicate the corresponding 95% confidence intervals. (a) Bacteria. (b) Chlamydomonas reinhardtii. (c) Colpidium striatum. (d) Dexiostoma campylum. (e) Euglena gracilis. (f) Euplotes daidaleos. (g) Paramecium caudatum. (h) Rotifer sp.

The number of predictors used in the best forecast model of the abundance of a species was independent of the number of interactions of said species (t 130 = 0.47, p = 0.637, Figure 3a, Table S3) regardless of the temperature regime (t 137 = −0.26, p = 0.793). Across replicates, when the temperature was constant 38.9% of the best models had three predictors, 27.8% had two predictors and 22.2% had four predictors, and similarly, when the temperature varied in most cases the best models had three (43.1%), two (30.6%) and four (15.3%) predictors (Figure 3b). The sensitivity analysis for the embedding dimension confirmed these results (Figure S16). Adding predictors to the best forecast models did either worsen or not improve them (Figure S6).

FIGURE 3

The best forecast models consistently contained only few predictors, regardless of the number of interactions. (a) The relation between the number of predictors in the best forecast models and the number of interactions of the predicted species. The black lines represent the fit of the regression model, and the shaded regions indicate the corresponding 95% confidence intervals. (b) Bar‐plot of the number of predictors in the best forecast models. The left panels are for the constant temperature case and the right ones for the fluctuating temperature case.

Interactor strength versus predictor importance

Overall, we found no relation between the interaction strength of a state variable with a target species and the forecast error of the abundance of the target species with the state variable as the sole predictor (F = 2.30, p = 0.129, Figure 4, Table S4), regardless of temperature regime (F = 0.05, p = 0.822). An alternative analysis confirmed these results (Figure S12). The estimated slopes varied across target species from −0.523 to 0.165 (F = 7.42, p < 0.001), but only for C. reinhardtii and E. gracilis there was evidence that the stronger interactors predicted these targets better (Table S8) and that the CCM skills (used to determine whether the species interacted, see above) positively correlated with interaction strengths (Section S7.1).

FIGURE 4

How strongly a species interacts with a target does not influence how well the species predicts the abundance of the target. The forecast error (RMSE) of the abundance of a target species with a model having only one species as predictor as a function of the interaction strength between the predictor and the predicted target. The black lines show the fit of the regression model, and the shaded regions indicate the corresponding 95% confidence intervals. The left panels are for the constant temperature case and the right ones for the fluctuating temperature case. (a) Bacteria. (b) Chlamydomonas reinhardtii. (c) Colpidium striatum. (d) Dexiostoma campylum. (e) Euglena gracilis. (f) Euplotes daidaleos. (g) Paramecium caudatum. (h) Rotifer sp.

DISCUSSION

We found that the forecast skill of the abundance of a species increases the more interactions the species has within the system (e.g. with other species) but also that it increases the weaker these interactions are on average. We found that these two measures—the number of interactions and their mean strength—are negatively correlated resulting in the abundance of species with many but on average weak interactions to be predicted the most skilfully. While the fluctuating temperatures did not influence these findings, they lowered the median skill of forecasting the abundances of three out of eight target species.

Previous studies reported contrasting results as in some cases predictions improved with increasing system complexity (Berlow et al., 2009; Dornelas et al., 2011; Iles & Novak, 2016; Mougi, 2017), while in others the opposite was the case (e.g. Doak et al., 2008; Jonsson et al., 2018; Novak et al., 2011). The latter led to the statement that ecological forecasting is limited by the low intrinsic predictability of real‐world systems due to their great complexity (Beckage et al., 2011). In our study, the addition of complexity (i.e. the fluctuating temperature) to the system lowered forecast skill for some but not all species. Thus, a universal association between an increasing system complexity and the predictability of the abundance of its components is less likely to exist. Future studies could investigate the effects of other complexifying factors.

In contrast, the more connected a species was in the system, the better it was forecasted. The negative correlation between the number of interactions and the mean interaction strength could explain this as it might indicate that species with many but weak interactions were less dependent on the state of individual system components and more dependent on the state of the whole system. Regardless, this result provides a first insight into why certain aspects of ecological systems are more predictable than others (e.g. some species abundances more than others; Harris, 1994). It suggests that species with few, strong interactions should be sampled more frequently than those with many weak interactions to achieve a comparable forecast skill. Thus, it has the potential of improving the monitoring of species in real‐world ecosystem, which can be a costly endeavour (e.g. Jones, 2011; Manley et al., 2004).

Yet, it remains unclear why a species with few strong interactions is not predicted more skilfully than a species with many weaker interactions. However, this result is corroborated by the finding that interaction strength is not a good indicator of how well an interacting variable predicts the abundance of the target species. Knowledge about good predictors of species abundances could help our understanding of ecological forecasting and our skill in carrying it out (Petchey et al., 2015). Based on these results, it is likely that interaction strength can be excluded as an indicator of good species abundance predictors.

In this context, we found that the median forecast skill increased the more system variables we included as predictors in the forecast models. However, in most cases, we achieved the best forecast skill already with few predictors included (i.e. between two and four predictors in approximately 90% of forecast models), regardless of how many interactions the forecasted species had. This suggests that if it is known a priori which system components are good predictors of the abundance of a specific species, then data collection can be streamlined by focusing on these variables rather than on the whole system.

In our experiment, we compared fluctuating and constant temperatures as the former are more truthful to natural conditions. The lower predictability of the abundance of some species in the fluctuating temperature setting when compared to the constant temperature setting suggests that prediction skill might be overestimated in experiments in which temperatures are constant. In fact, in laboratory‐based or simulated time series experiments, the temperature is usually kept at one or more constant levels (e.g. in Yeo et al., 2003; Ferguson & Ponciano, 2014; Daugaard et al., 2019) and only rarely fluctuating temperatures are used (e.g. Descamps‐Julien & Gonzalez, 2005; Jiang & Morin, 2007). Given that temperature is a strong driver of species metabolic rates (Brown et al., 2016) and thus also of their dynamics (e.g. Bernhardt et al., 2018; Lee et al., 2007), fluctuating temperatures should more frequently be considered to better reconcile results from laboratory or simulation experiments with real‐world insights.

The distribution of interaction strengths in a system is known to be right‐skewed, with the bulk of the interactions being weak and only comparably few interactions being strong (e.g. Paine, 1992; Wootton, 1997; Wootton & Emmerson, 2005). This was also the case in our study (Figure S9). Moreover, our finding that the number of interactions and the average interaction strength are strongly negatively correlated strengthens the recent empirical evidence (Ratzke et al., 2020; Ushio, 2022) of the theoretical finding that generalists have predominately weak interactions while specialist are responsible for the right‐skew of the interaction strength distribution (Wootton & Stouffer, 2016). Given that weak interactions have been identified as system stabilising (e.g. McCann & Hastings, 1997; Neutel et al., 2002; Otto et al., 2007), our results support previous observations stating that generalist species have a stabilising function due to the weak interactions they engage in (e.g. Mougi & Nishimura, 2007; Chakraborty, 2015; Brechtel et al., 2019, note, however, that we did not carry out a stability analysis of the system in this study).

Several robustness analyses confirmed that the results are not sensitive to the specifications of the experimentation and analyses (e.g. the choice of forecast method did not affect the results, Section S8.5.1). Noticeably, the analyses of the potential influence of the different measurement methods on the results revealed that the main results most often still occurred within measurement method (Section S8.4.6). However, the 95% confidence intervals often overlapped zero due to small sample sizes. The analyses also showed that any effect of measurement method was not due to anything as simple as differences in measurement error across methods (Section S8.4.5). Since measurement method is confounded with species identity, we cannot tease apart their possible influences on the main results, and therefore cannot completely rule out that measurement method has, for some reason, some power in explaining the main results reported.

In conclusion, we provide novel insights into why the abundance of some species is better predictable than others in the same system. The dependency of forecast skill on the number and the strength of species interactions not only improves our knowledge of ecological forecasting. It has also the potential of improving the resource allocation for the sampling and monitoring of species, as comparable forecast skill across species likely requires varying amounts of data per predicted species based on how much and how strongly this species interacts. We also shed further light on the relationship between elements of system complexity and forecast skill, showing that the relationship can be both species‐specific and of different signs within and across systems. Thus, forecasting skill may deteriorate with increasing complexity, but this cannot be taken for granted and can depend on whether one is comparing across or within systems.

AUTHORSHIP

DI, FP and OLP conceived the experiment. UD, SBM and OLP conceived the research questions and designed the methodology. DI collected the data. UD analysed the data and led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/ele.14070.

OPEN RESEARCH BADGES

This article has earned an Open Data badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at [https://doi.org/10.5281/zenodo.6524695].

Appendix S1

Click here for additional data file.