Figure 2. Schematic representation of differences between direct interactions and higher-order interactions (HOIs). Direct interactions involve two plants, a neighbor and a target, while HOIs (in their simplest form), include a neighbor pair and a target plant. For example, although species A (neighbor) may outcompete species B (target) in the absence of species C (negative direct interaction), a mixture also including species C (neighbor), reducing the intensity of competition between A and B (less negative interaction) by supplying a scarce resource to both, may actually benefit (facilitate) the growth of species B (positive HOI).
Exposure to more extreme dry and wet spells may lead to fundamental changes in species interaction networks (Grant et al., 2014; Olsen et al., 2016). The stress-dominance hypothesis predicts that the relative importance of environmental filtering increases while competition decreases along a gradient of environmental stress (Coyle et al., 2014; Weiher & Keddy, 1995). By extension, the stress-gradient hypothesis presumes an increase in the relative importance of facilitation compared to competition along such gradients (Bertness & Callaway, 1994). Despite recent evidence questioning the generality of these hypotheses (De Boeck et al., 2018; Grant et al., 2014; Maestre & Cortina, 2004; Metz & Tielborger, 2016; Singh & Baruah, 2021), directional shifts in mean interactions along environmental gradients have been observed in many ecosystems, indicating important changes in ecosystem functioning (Coyle et al., 2014; De Boeck et al., 2018; Olsen et al., 2016).
We explored if and how increasing persistence in summer precipitation regimes (i.e. precipitation variability) changes the number, strength and directionality of species interactions that influence plant biomass production, including both direct interactions and HOIs. We subjected experimental grassland mesocosms with 12 common C3 species (six grasses, three non-N-fixing forbs, three legumes) to an 8-level gradient of precipitation regimes. These regimes varied the duration of consecutive dry and wet periods from 1 to 60 days while keeping total precipitation constant (Reynaert et al., 2021). Our research questions were whether species interactions, either direct ones or HOIs, determine species biomass under altered precipitation persistence, and whether these interactions change as drought stress intensifies (i.e., the main driver of plant responses to increasing summer precipitation persistence; Reynaert et al. (2021); Reynaert et al. (2022)).
We hypothesized that (1) HOIs would generally counteract direct competitive interactions, weakening net species interactions and stabilizing community dynamics; (2) species interactions would overall become more positive (or less negative) across the community with increasing drought severity, indicating either a decline in competitive interactions and potentially increased facilitation, or a more prominent role of environmental filtering, leading to increased competitive release; (3) species identity would determine the nature of interactions, with more facilitation in interactions with N-fixing legumes and more competition when considering dominant grasses.
Methods & materials
Data collection
From July to November 2019, a 120-day open air precipitation manipulation experiment was conducted in Flanders, Belgium (51°09′41″N, 04°24′9″E). In total, 256 mesocosms (29.5 cm inner diameter and 50 cm depth), each containing the same synthesized grassland plant community with three individuals of 12 common C3 species (36 plants per community; Fig. S1), were evenly distributed across eight identical plots of 3 m diameter and subjected to a gradient of increasingly longer dry-wet period alternations (Reynaert et al., 2021). The precipitation persistence gradient contained eight different regimes of dry-wet period alternations, ranging from 1 to 60 consecutive days with or without precipitation and starting with either a dry or a wet period. This yielded 16 regimes in total, including 1, 3, 6, 10, 15, 20, 30 and 60 consecutive days dry/wet. All regimes had a total of 60 days with and 60 days without precipitation. Historically, Belgian weather patterns following 1-day dry-wet alternations are most common (KMI, 2019; Reynaert et al., 2021). The precipitation regimes were created by blocking out ambient precipitation with rainout screens which were automatically deployed only during rain, and applying dripper irrigation (Reynaert et al., 2021). Precipitation was equal in all regimes on each irrigation day (6.87 L/m2, being 1.5 times the Belgian daily average of a rainy day to compensate for excess evapotranspiration in mesocosms) as well as in total across the experimental period. Mesocosms had holes in the bottom, allowing free drainage and preventing flooding (Reynaert et al., 2021). Per regime, a total of 16 replicate mesocosms was evenly distributed across 4 plots. The position of these replicates was different in every plot to account for potential edge effects.
The 12 grassland species represented a broad gradient in plant functional traits (Table S1) and included six grasses (Agrostis capillaris  L. (AC), Anthoxanthum odoratum  L. (AO), Deschampsia cespitosa  (L.) P. Beauv. (DC), Phleum pratense  L. (PHP), Poa pratensis  L. (POP), Holcus lanatus  L. (HL)); three N‐fixing forbs (Lotus corniculatus  L. (LC), Trifolium pratense  L. (TP), Trifolium medium  L. (TM)); and three non‐N‐fixing forbs (Centaurea jacea  L. (CJ), Lychnis flos‐cuculi  L. (SF), Plantago lanceolata  L. (PL)). Seeds were obtained from a seed company in the Netherlands (Cruydt-Hoeck) and sown in separate seedling containers at the start of April 2019. By May, viable seedlings of similar size per species were transplanted directly into the mesocosms (Reynaert et al., 2021). We maximized interspecific interactions by planting them in a hexagonal grid where each individual had different neighbors at 4.5 cm interspace, and avoiding clumping (Fig. S1). The spatial configuration in our containers was aimed at representing each species pair at the most comparable extent possible, with 92% of the 66 possible pairs occurring either one or two times (Table S2). Species positions and north-south orientation were equal in all mesocosms. Mesocosms were weeded regularly over the course of the experiment to prevent (re)colonization.
Average volumetric soil water content (SWC) over 0-30 cm depth was recorded half-hourly in two to four mesocosms per precipitation regime by a CS650-DS Reflectometer (Campbell® Scientific INC., Logan, Utah, USA). Because the impact of increasing weather persistence was shown to be mainly related to drought in the 2019 growing season (Reynaert et al., 2021), the daily averages of mean SWC values per regime were utilized to calculate different drought indices which represented the persistency of the precipitation regime. These included drought stress intensity (at relative extractable water (REW) of 0.4, 0.2 and 0.1), total time below permanent wilting point, and the length of the longest consecutive period below permanent wilting point (Reynaert et al., 2021; Reynaert et al., 2022; Vicca et al., 2012). Following Vicca et al. (2012), drought stress intensity was calculated by determining a soil moisture threshold, and integrating the duration and extent to which soil moisture declined below it over the experimental period.
At the end of the experiment, all standing dead and live biomass above 4.5 cm height was harvested per species per mesocosm, in 12 out of the 16 replicate mesocosms per regime (all available individuals of the same species per mesocosm per paper bag). Dead biomass was included since death at the time of measurement, does not indicate a lack of previous interact with its neighbors. In half of these 12 mesocosms, one third of the individuals per species was randomly removed for another study. We corrected for this by multiplying the plant mass in these bags by 1.5 (per regime, every species was discarded twice in every unique position; Fig. S1). Combining the biomass of all harvested individuals per species in a mesocosm accounted for variation caused by differences in local neighbor setting. Biomass bags were oven-dried at 70 °C for >72 h and immediately weighed with 0.01 g accuracy. To achieve dry plant mass per species per mesocosm, we subtracted the average mass of 58 oven-dried empty bags (SD = 0.08 g). Plants that had little or no biomass above 4.5 cm height because they remained too small during the experiment (< 17 %, n = 2281) were given a mass of 0.01 to account for being below our detection limit (Lubbe et al., 2021).
Statistical analysis
All statistical analyses were performed in R version 4.0.2 (R Core Team, 2019). Significance was assumed for p-values < 0.05. Graphs were constructed utilizing the packages dplyr (Wickham et al., 2015), igraph (Csardi & Nepusz, 2006) and ggplot2 (Wickham, 2016).
To model direct interactions and HOIs, we adapted the approach from Mayfield and Stouffer (2017), which allows the use of biomass data. Their methodology suggests to construct models with the characteristic of interest of a target species in the community (here: biomass) as response variable and (the biomass of) all neighboring species and their statistical interactions as explanatory variables, after which all significant terms from the most parsimonious equivalent models (e.g., Δ AIC < 2) per target species are of interest. The estimates of significant model terms can then be utilized to obtain information about species interactions, direct ones and HOIs, in the studied plant community. Because biomass data was pooled per species per mesocosm in our experiment, we did not focus on the role of individual plant position but rather on species identity. Our central question was if and how the growth performance of a target species, expressed as accumulated dry biomass over the growing season, was affected by the growth performance of neighboring species under increasingly persistent precipitation regimes in diverse grassland communities.
As described earlier, we opted to represent our treatments in the models by a drought index. Out of all calculated drought indices (section 3.1), the length of the longest consecutive period below permanent wilting point captured species biomass effects caused by increasingly persistent summer precipitation regimes the best (lowest AIC, adj. R² = 0.56; Table S3), allowing us to more accurately linearize the biomass responses (Fig. S2 & S3) and model the wet and dry start regimes along a common dimension (Reynaert et al., 2021).
First, we normalized all biomass data by taking the natural logarithm of each biomass observation divided by the mean biomass per species across all regimes. This approach is similar to the log-response ratio, a technique often applied in meta-analysis, where values >0 indicate above average plant performance, and values <0 below average (Van Sundert et al., 2021). By doing so, effect sizes become comparable among species with distinct intrinsic characteristics, i.e. the coefficients of models with different target species are standardized and directly comparable. Since all effects become relative, they can be interpreted as the extent to which a change in mass of a neighbor species (relative to its mean size) affects the mass of a target species (relative to its mean size). We then created multiple linear models of all possible species combinations including every species as response variable (target) and all others (neighbors), the drought index and their two-way interactions as explanatory variables (Fig. 2; Eqn. 1).
𝔼 (BMA,i) = Intercept + α1BMB,i + α2BMC,i + α3BMB,iBMC,i + γIndex + β1BMB,iIndex + β2BMC,iIndex + β3BMB,iBMC,iIndex + … + ε i (Eqn 1)
With 𝔼 (BMA,i) the expected standardized total dry biomass (BM) of target species A in a specified mesocosm i, BMB,i and BMC,i the standardized total dry biomass of neighbor species B and C, respectively, in that mesocosm, Index the length of the longest consecutive period below permanent wilting point, α and β the species interaction (model) coefficients and ε i the errors. Model terms including only one species (BMB,i or BMC,i) represent direct interactions. Model terms including both species (BMB,iBMC,i) indicate HOIs.
To avoid overfitting and other issues related to insufficient data (Martyn et al., 2021), we only included the simplest, i.e. the three-way, HOIs (Fig. 2). Intraspecies interactions were excluded from the models because we pooled the species biomass for the three individuals per mesocosm, which accounts for variation caused by individual plant position. Following Grant et al. (2014), we interpret the general variation in interspecies interactions knowing that distant (i.e. non-direct neighbor; Fig. S1) intraspecies effects may also contribute to the estimated interaction coefficients. Because plant density and positions are the same in all mesocosms we assume that found model coefficients primarily reflect variation in plant performance related to changes in interspecies interactions and drought effects. Per response variable (target species), we then compared full models with direct-interaction only and null models (only including the drought index as explanatory variable), retaining only the most parsimonious model per target species (Δ AIC < 2). In the remnant models, the strength and direction of the interaction(s) was determined by the model estimates of the explanatory variable(s). Significant model coefficients excluding the drought index (α coefficients; Eqn 1) represent interaction strengths at a drought length of zero. All other significant model coefficients (β coefficients; Eqn 1) indicate a change of interaction strength along the gradient of increasing drought length. To calculate interaction strengths at maximum drought index, the maximum index value (39), was used to fill in Eqn 1. If the model coefficient was not significant, the interaction strength was assumed to be zero. By doing so, we were able to isolate the most important direct interactions and three-way HOIs per target species in the plant community.
Although our models suggest the capacity to detect unidirectional effects of species on each other and thus causation (influenced species as dependent variable or target, influencing species as independent variables or neighbors in Eq. 1), they can actually only detect the balance of these unidirectional effects, and thus correlation. We therefore considered every found interaction coefficient as the bidirectional outcome of the (potential) unidirectional effects, regardless of whether the species of interest was included as target or neighbor in the model.
Utilizing the interaction coefficients obtained in this way, we then calculated the total number of interactions that involved this species (either as target or neighbor in the models). This was done separately for direct interactions and HOIs. Direct interaction coefficients between the same species pair originating from different models were only counted once per species. Other coefficients represented unique and equally relevant 3-way HOIs and were all counted per species. Next, we calculated the mean direct interaction and HOI strengths by averaging all positive and negative interactions per species, again utilizing the mean values of coefficients between the same species pairs for direct interactions and all coefficients for HOIs. We then visualized species interaction networks for direct interactions, HOIs, net interactions (sum of direct interactions and HOIs per species pair) and their changes along the drought gradient. Finally, to test if species characteristics determined any of the observed patterns, we conducted various non-parametric Spearman correlation tests, utilizing mean interaction strengths and total number as response variables and average species biomass (at zero or maximum drought index) as well as Ellenberg nutrient and moisture indices as explanatory variables.
Results
Species interactions (including HOIs) played an important role in determining species performance, expressed as aboveground biomass production, under changing precipitation regimes. Out of the 12 most parsimonious models (ΔAIC < 2), 75% contained significant direct interactions and 92% significant HOIs (Table S4). Null models, which only included the drought index as explanatory variable, explained 10% of the variance in species biomass on average (Table S4). Inclusion of direct interactions increased the variance explained to 28% on average (Table S4). Full models, also including HOIs, explained 43% of variance in target species biomass on average (Table S4).
At zero drought index, the majority of direct interactions (75%, n = 24) (stacked bars above zero; Fig 3a) were positive compared to only half of the HOIs (50%, n=108) (Fig. 3b). Although a larger proportion of significant direct interactions became more negative (50%, n = 12) (stacked bars below zero; Fig. 3c) with increasing drought length (i.e., index change), there were still more positive mean direct interactions at maximum drought index overall (72%, n = 11; Fig. 4c). For HOIs, the total number of positive and negative interactions stayed similar (52% positive, n = 75) (Fig. 3d), though on average, more species were involved in facilitative HOIs at maximum drought index (58%, n = 12; Fig. 4d). Whereas opposite interaction effects generally neutralized each other at zero drought index, overall, competition and facilitation became more pronounced with intensifying drought stress (Fig. 4). As such, average direct competition increased in strength and prevalence along the drought gradient while more HOIs became facilitative (Fig. 4). Despite these opposite changes, mean direct interaction and HOI strength did not inversely correlate at either side of the drought gradient (Fig. S4).