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).