Statistical analyses
First, we tested the overall effects of plant species richness and plant
functional group richness on the response variables using linear mixed
models. Total values were calculated as the sum of each response during
the entire sampling duration. We then tested the overall relationship
between plant diversity and phenological synchrony. As a measure of
phenological synchrony, we calculated Pearson’s correlation index
(r ) between all possible pairs of activities by plot (6 pairs in
total). Correlations were calculated based on data collected during the
maximum temporal extent possible; that is, for aboveground-aboveground
and aboveground-belowground correlations for the entire growing season,
and belowground-belowground correlations, for the entire year.
Significant correlations are referred to here as positively or
negatively coupled, when the direction of the correlation is positive
and negative, respectively. A decoupled phenology here implies
independence, i.e. non-significant relationships. We fitted individual
linear mixed-effect models with the correlation coefficients as response
and plant species richness or plant functional group richness as
predictors. We treated block and plot as nested random factors for both
analyses, and species richness was log-transformed.
To evaluate whether time of year altered plant diversity effects on the
strength, direction, and predictability of relationships among
aboveground and belowground response variables, we fit a series of
mixed-effect models into a structural equation model (SEM; Grace, 2006;
Eisenhauer et al., 2015), following the conceptual framework depicted in
Fig. S1. Given that the response variables were sampled in different
periods of the year, and differed in pattern (aboveground shoot
responses were unimodal, while belowground responses were bi- or
multimodal), we subset the dataset into three parts. The “spring”
dataset encompasses the beginning of the aboveground measurements until
the first mowing (14/Jul). The “summer” dataset encompasses the
sampling after the first mowing until the second mowing, when the
aboveground measurements finished (Fig. S1a). Finally, the “winter”
dataset encompasses only belowground phenology, from the end of
September until February (Fig. S1b). To simplify models and to avoid
multicollinearity, we ran a stepwise selection of variables for each
mixed-effect model within the SEM, using the ‘step’ function. The
resulting models were then used to build the initial SEM, using the
piecewiseSEM package (Lefcheck, 2015). We inspected the initial SEM
model results according to the goodness-of-fit tests for both the SEM
and individual causal relationships and selected the final model by
excluding the insignificant factors and adding missing relationships
that significantly improved the model’s global fit. We treated block and
plot as nested random factors, allowing responses to vary randomly
between blocks and plots. Given that samples taken from closer sampling
point times are more alike, we also accounted for temporal
autocorrelation fitting a ‘corCAR1’ term in each model. The variables
mean plant community height, root production, and species richness were
log-transformed. Due to the relationship between species richness and
functional group richness, we have also incorporated correlated errors
between those variables. We assessed the homogeneity of residuals with
residuals vs. fitted values plots and Q-Q plots for data normality using
‘Pearson’ correlation (Zuur et al., 2009) for each of the mixed-effect
equations used in the SEM. Statistical analyses were performed with R v.
4.2.2 (R Core Team, 2022).