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