Differential responses of community-level functional traits to mid- and late-season experimental drought in a temperate grassland

Extreme precipitation events are becoming more intense and frequent due to climate change. This climatic shift is impacting the structure and dynamics of natural communities and the key ecosystem services they provide. Changes in species abundance under these conditions are thought to be mediated by functional traits, morpho-physiological characteristics of an organism that impact its fitness. Future environmental conditions may, therefore, favour different traits to those in present-day communities. After six years of manipulated precipitation levels, including drought (−50% of ambient precipitation), irrigation (+50% of ambient precipitation), and control (ambient precipitation), we measured five key functional traits (plant height, leaf dry matter content [LDMC], leaf thickness, specific leaf area [SLA], and leaf phosphorus concentration) in 586 individual vascular plants to study the effects of precipitation changes on community weighted functional traits. Additionally, we tested whether the precipitation change effects on the traits depend on the time of the growing season. As expected, reduced precipitation impacted community composition only for the late season timing, after the seasonal field mowing, but led to a significant change in all community-level plant traits between season timings. Under drought, communities shifted towards shorter individuals with thicker but small leaves and lower phosphorous content. Overall, a combination of community reassembly and intraspecific variation contributed to community-weighted differences between control and drought plots for plant height, SLA, and LDMC traits. Species turnover was the main driver of community-weighted means shifts in all traits in the late-season but SLA. Whereas all traits showed variations at the community level with drought, SLA and LDMC were the most responsive traits at the species level. Nevertheless, our results suggest underestimation of intraspecific variation due to sensitive species lower abundance under stress. No differences in community-weighted means of functional traits were observed between control and irrigated plots.

Synthesis. Our findings suggest that functional trait composition of grassland communities may shift under climate change-induced drought, depending on the growing season timings . Trait-based attempts to predict ecosystem functioning must account for such temporal variation in community trait values.

Methods
Study site

The experiment was conducted at the Upper Seeds field site (51°46'16.8"N 1°19'59.1"W, 155 m a.s.l) in Wytham woods, Oxfordshire, UK. The study site is a calcareous temperate grassland characterised by a shallow soil depth (300-500 mm depth) alkaline soils (Gibson & Brown, 1991), a daily average temperature ranging between -5 °C and 26 °C, a mean annual temperature of 11.5 ºC (2016-2021), daily total precipitation of 0-40 mm (2016-2021) and an annual total precipitation of 686 mm (2016-2021) (Table S1). According to the Köppen-Geiger climate classification, the study site constitutes a maritime temperate climate (Cfb) (Peel et al., 2007). The site is currently managed for maintenance by mowing twice per year, a common practice across most European grasslands (Török et al., 2018). The first mowing takes place mid-growing season (late July), and the second at the end of the growing season (late September).

To test our hypotheses, we experimentally manipulated precipitation levels using the RainDrop (rainfall and drought platform) long-term ecological experiment for six years, and in July and September 2021 we determined species abundance and sampled the most abundant species using functional traits measurements. RainDrop is integrated in the DroughtNet global coordinated research network (https://droughtnet.weebly.com/). The RainDrop experiment has been running since 2016 and consists of twenty-five 5 m × 5 m permanent plots distributed across the study site in five randomised blocks (A:E) (Figure S1). Each block includes five plots that receive one of the following treatments during the growing season (March-September): drought (−50% ambient precipitation), irrigation (+50% ambient precipitation), two ambient control plots (no manipulation), and procedural control (inverted rain shelters). Rain shelters intercept 50% of precipitation for the drought plots. In each drought plot, precipitation is intercepted by gutters and collected in containers situated next to each rain shelter. Pipes connect these deposits to sprinklers that spray the water onto the adjacent irrigation plots. This design is based on the proposal of Yahdjian & Sala (2012) and Gherardi & Sala (2013), which has been applied across >100 nodes of the DroughtNet network worldwide for exploring climate change-induced drought and deluge effects on community composition (Fischer et al., 2013; Smith et al., 2024). Furthermore, specific to our study site, climatic projections forecast shifts in annual precipitation from -20% to +20% under the 8.5RCP for 2080-2099 compared to 1981-2000 (UKCP18 Project, Met Office). Hence, the plots that experienced no precipitation manipulation served as ambient precipitation treatment (control) plots. Additionally, to control for shelter effects, each block has one procedural control plot. These consisted of rain shelters with inverted gutters that allow 100% of precipitation to pass through.

Ongoing work at this field site has revealed no differences in community composition between the procedural and ambient controls (Jackson et al., 2024). Therefore, we did not measure traits from the procedural control plots, and instead sampled one additional ambient control plot per block. To minimise edge effects, we split the 5 m × 5 m plot into four quarters and marked out a 1 m × 1 m quadrat in the centre of the study quarter from which we took all trait and abundance measurements. Data were collected at two different times along the growing season: mid-growing season (July 2021) and late-growing season (September 2021), just before each seasonal mowing.

Data and metrics

Species abundance

We collected species abundance data to assess community composition dissimilarity between the precipitation treatments (H1) and possible differences in community-weighted means between treatments (H2). To do so, we quantified species percentage cover for all vascular plant species and bare soil percentage cover using a 1 m × 1 m quadrat (with 10 cm grid) at each examined permanent plot. We first estimated the percentage cover independently for every species in each quadrat (one per plot), and next transformed these estimates into relative abundances that sum to 100%. Because the mid-season mowing removed the inflorescence from all grasses, species identity was difficult to know for many graminoid species during the late-season, which may impact the observed community effects. In the late season, only two of the graminoid species (Brachypodium pinnatum and Brachypodium sylvaticum) were identifiable at the species level because of their distinctive leaves. For these two Brachypodium species, we recorded percentage cover as described above. Separately, we recorded the pooled abundance of all other graminoid species. Although we identified graminoid species in the mid-season, to ensure sound comparisons between mid-season and late-season abundance data (NMDS analysis), we combined the mid-season abundance of non-Brachypodium graminoids in analyses to contrast community-weighted means among precipitation treatments at the different times.

Trait measurement and community weighted means

We quantified key functional traits of the grassland community to test whether CWM trait values differ among precipitation treatments and whether the differences are driven by different factors (interspecific vs. intraspecific trait variation) at the two contrasted growing seasons timings (H2). Trait measurement was performed at the same time species abundance was measured (July and September 2021). We measured height, specific leaf area, leaf thickness, leaf dry matter content and leaf phosphorus on the most abundant species in each quadrat, totalling 586 individual plants. We measured traits using a standardised protocol (Pérez-Harguindeguy et al., 2013), briefly summarised in Table 1. For the selected species to be representative of the community, we focused on species that contributed to a cumulative abundance of at least 80% within each quadrat, following Garnier et al. (2004) and Pakeman & Quested (2007). Next, in each quadrat, we randomly selected three mature, healthy individuals per species for trait measurement. Leaf traits were measured on one young but fully developed leaf per individual with the exception of leaf phosphorus, which require pooling three leaves per species per plot to obtain the 50 mg of dry weight required for the analysis (Esslemont et al. 2000). After pooling, 69 samples (out of a total of 195 potential samples) reached the dry weight threshold and thus were eligible for analysis. Leaf phosphorus concentration was estimated using inductively coupled plasma mass spectrometry (ICP-MS) following Esslemont et al. (2000). Briefly, samples were digested with 1 mL of concentrated nitric acid and 0.7 mL of hydrogen peroxide at 50ºC overnight and diluted 25 times with MiliQ water prior analysis.

Community-weighted means of each trait were calculated in each quadrat for all treatments. CWM are commonly used in trait-based ecology to quantify shifts in community mean trait values due to macro and micro-environmental selection (Garnier et al., 2004; Bruelheide et al., 2018; Griffin-Nolan et al., 2019; Kambach et al., 2023). We calculated CWM of the five assessed traits using the framework proposed by Lepš et al. (2011), and the trait.transform and trait.CWM functions from Lepš et al. (2011) and  Götzenberger et al. (2020), which are included in cati R package (Taudiere and Violle, 2016). This framework disentangles the effect of interspecific variability (species turnover) and the combination of species turnover and their intraspecific trait variability as fixed community-weighted means (fixed CWM) and specific community-weighted means (specific CWM), respectively. Briefly, to calculate the specific CWMs, we multiplied the mean trait value per dominant species in each treatment by each species’ relative abundance in the quadrat. The resulting species products were summed, and their abundances rescaled following de Bello et al. (2021). Our target of sampling species with a cumulative abundance of 80% was achieved for all traits in all quadrats, with the only exception of leaf phosphorus, for which the threshold was achieved in 25% of quadrats only (Table S2). Because we observed that the measured traits did not vary significantly across the growing season for any of the studied species, we used the same trait values for both growing season timings, with the exception of plant height. Due to the graminoids identification difficulties at the late-season, and because all graminoids (Agrostis capillaris, Arrhenatherum elatius, Brachypodium pinnatum, Brachypodium sylvaticum, Dactylis glomerata, Holcus lanatus, and Trisetum flavescens) showed similar height in comparison to other functional groups in the mid-season, we used the height of Brachypodium pinnatum, the most abundant identifiable graminoid (ca. 34%, whereas the non-Brachypodium graminoids showed relative abundances <18% at the mid-season), as the height of late-season non-Brachypodium graminoids. To calculate fixed CWM, a single mean trait value for individual species was used for all quadrats. While changes in fixed averages across treatments would reflect species turnover, changes in specific averages reflect both between and within-species variability in traits (Lepš et al., 2011).

When trying to test for differences at the species level (intraspecific variation), not all species presented a consistent number of samples, as their abundance might depend on the sensitivity to the treatment. We therefore restricted analysis of each species responses to different treatments (i.e., precipitation-driven intraspecific variation) to the species that had enough replication to do so, considering the recommended replicate number to account for natural trait variation (n= 30 at each precipitation level, Pérez-Harguindeguy et al., 2013). These species consisted of three graminoids (Brachypodium pinnatum, Trisetum flavescens, and Arrhenatherum elatius), three legumes (Medicago lupulina, Trifolium repens, and Trifolium pratense), and one non-leguminous forb (Crepis capillaris).

Statistical analysis

H1: Do precipitation treatments affect community composition?

We performed non-metric multidimensional scaling (NMDS) to contrast community composition between treatments (H1) using the vegan R package (Oksanen et al., 2020). NMDS is a form of dimension reduction that allows for differences in communities to be quantified. In its application to plant trait-based ecology, NMDS is based on the rank-order of species abundances and maximises the correlation between real-world distance and Bray Curtis distance in the ordination space. We assessed the significance of any differences in community composition using analysis of similarity (ANOSIM). To determine which species were responsible for any dissimilarities among treatments, we used similarity percentage (SIMPER). Differences between groups when performing multidimensional analysis (ANOSIM and SIMPER) were considered significant when the p-value was less than 0.05.

H2: How do precipitation manipulations alter CWM of traits in different periods of a season?

We fitted hierarchical linear mixed-effects models to test our expectations that both the specific and fixed CWM of each of the five assessed traits would differ between the three precipitation treatments in mid-season and late-season periods, but the pattern of the differences for each trait would differ between the two growing season timings (H2). We used the R package lme4 to account for the blocked experimental design, where models were fit using maximum likelihood (ML) (Bates et al., 2015). Precipitation treatment (three levels: drought, control, irrigated) and timing in the growing season (two levels: mid- vs. late-season) were fixed effects, and block was considered as a random effect (five levels: A:E, where all treatments are represented). A significant interaction effect of precipitation treatment and seasonal timing indicated that the impact of the precipitation manipulation differed between the two timings. To meet the assumptions of gaussian distributions in our models, we log-transformed plant height, while other traits’ data were not transformed. When interpreting the output of mixed-effects models, we focused on differences based on the 95% confidence intervals overlap rather than relying on p-values following Flechner & Tseng (2011) and Bates et al. (2015), who discouraged the use of p-values in mixed-effect models.

 Finally, to explore the link between trait and species abundance precipitation sensitivity, we contrasted the observed relative abundance decrease under the drought treatment (as we observed significant differences, see results) with the position of all species along the two principal components of species traits variation using principal component analysis (PCA) from the PCAtools R package (Blighe & Lun, 2022). In addition, to further explore how specific CWM differ between treatments in each species, we fitted linear mixed-effects models, with only precipitation treatment as a fixed effect for individual species with enough replication across treatments.

All analysis were performed in R (R Core Team, 2022, v. 4.2.1).