Tree species diversity in forest ecosystems could reduce their vulnerability to extreme droughts through improved microclimate and belowground water source partitioning driven by contrasting species-specific water use patterns. However, little is known about the seasonal dynamics of belowground water uptake that determine whether diversity positively or negatively impacts tree carbon assimilation and water exchange.

Using a network of 30 permanent plots in Mediterranean forests with increasing tree species diversity (from monospecific to four-species mixtures), we examined the seasonal patterns of in-situ aboveground carbon and water relations and belowground water sources on 265 trees from four pine and oak species over two years using hydraulic and stable isotope approaches.

We found that increasing species diversity in broadleaf and conifer mixtures induced strong soil water source partitioning between oak and pine species. As conditions became drier during the summer in mixed stands, oak species took up water from deeper soil sources, while pines were systematically limited to shallow ones. Despite significant belowground moisture partitioning, stronger drought-induced reductions in photosynthesis, stomatal conductance, and leaf water potential were still observed in diverse compared to monospecific stands for pines but with some benefits for oaks.

Synthesis:

Our findings reveal that tree species diversity promoted belowground water source partitioning in mixed oak and pine stands, potentially reducing competition for water in more diverse ecosystems. Yet, our results show that it is insufficient to buffer the adverse impacts of severe droughts on aboveground tree carbon and water use, leading to higher water stress, especially for pines.

Site description

The study was conducted in Mediterranean forests in the Alto Tajo natural park (Guadalajara, Castilla La Mancha, 40.66°N, 02.27°W) in central Spain, where 30 plots (30m x 30m) within a 20 km² area were selected from the FunDivEUROPE Exploratory Platform (Fig. S1, Baeten et al., 2013). To assess the effect of tree diversity, the plots were established in 2011 in non-managed mature even-aged (i.e., more than 50 years) forests with limited variation in altitude (i.e., from 980 to 1300 m a.s.l.), topography, soil type, and density (see Table S1 and Baeten et al. (2013) for more details on plot selection). The soils in all plots were shallow (from 20 to 40 cm) calcic cambisol soils (FAO/UNESCO soil classification) on a cracked limestone bedrock but with plant roots that may extend down to several meters through the fractured bedrock (Peñuelas & Filella, 2003). This area has a continental Mediterranean climate with hot and dry summers and cold winters. The long-term annual precipitation sum was 416 mm (2011-2022), with 516 and 367 mm in 2021 and 2022, respectively. The long-term mean annual temperature was 11°C (2011-2022), with 12°C and 11°C in 2021 and 2022, respectively. Maximum daily temperatures were 31°C and 33°C for 2021 and 2022, respectively (Fig.1). The soil aridity index (P/PET) was calculated monthly using meteorological data from the nearest station (Molina De Aragon, Castilla La Mancha, ES; 20 km away, Fig. S1). We first computed the monthly potential evapotranspiration (PET) using the Thornthwaite equation (Yates & Strzepek, 1994) based on the average monthly temperature, daylight length, and heat index. P/PET was calculated by dividing the monthly precipitation sum by the monthly PET. P/PET varied from 0.34 and 0.38 in May (i.e., corresponding to wet soil conditions) to 0.11 and 0.05 in July (i.e., the driest and hottest period) for 2021 and 2022, respectively.

The area is characterized by the natural dominance of four tree species, i.e., two coniferous species (Pinus nigra subsp. salzmannii (Dunal) Franco and Pinus sylvestris L.) and two broadleaved ones (Quercus faginea Lam. and Quercus ilex L.), which were selected for this study. We selected plots with increasing tree species richness, including monospecific (where the target species represents more than 90% of the total basal area by itself), monofunctional (i.e., two-species mixtures of either coniferous or broadleaved species), multifunctional (i.e., two-species mixtures of coniferous and broadleaved species) and mixtures of the four species. Each richness level was replicated three times for all species with all possible species combinations, except for the two-species mixture between P. sylvestris and Q. ilex, resulting in 30 plots (Table S1). In mixed plots, the target species had similar abundances with a lower limit of 60% of maximum evenness in the basal area (Baeten et al., 2013). The understory vegetation (representing less than 10% of the total basal area) was mainly composed of shrub species (Arctostaphylos uva-ursi, Buxus sempervirens, and Genista scorpius) and juveniles of the dominant tree species.

We randomly selected five dominant or co-dominant trees per species in each plot, leading to 265 trees in total. To assess the seasonal dynamics of aboveground and belowground tree water and carbon use, we conducted in-situ measurements (detailed below) repetitively on each tree in 2021 and 2022 at the beginning (i.e., in May when soil moisture was high, Fig.1), middle (i.e., in July, corresponding to the driest and hottest period), and the end of the growing season (i.e., in September, representing the recovery transition from dry to wet soils) (n=6 measurements per physiological traits and tree in total). The two years presented contrasting climatic conditions with a rather mild summer drought in 2021 and an extremely long drought in 2022 (Fig. 1), allowing us to assess tree carbon and water dynamics under different environmental conditions. Indeed, the late rain events in 2022 led to lower soil moisture during the measurement campaign of September, while conditions were wetter during the last campaign in 2021, representing full recovery after the summer drought (Fig.1).

Leaf-level gas exchange and water potential

We measured the leaf-level light-saturated net photosynthesis (A_net, µmol m^-2 s^-1) and stomatal conductance (g_s, mol m^-2 s^-1) on one fully developed leaf (or multiple needles for conifers) per tree from a 50cm to 1m-long (for pine and oak species, respectively) sun-exposed branch. The branches were sampled using an extension pole, directly placed in a water bucket, and recutted under water to restore water flow (Lange et al., 1986). Within 15 min after sampling, A_net and g_s were measured using one infrared gas exchange analyzer (LI-6800 or LI-6400, LICOR Biosciences, USA), according to the method described by Bayar & Özçelik (2024). The measurements were done between 9 am and 1 pm (local time). The relative humidity inside the chambers was set between 30 to 50% (to match the average daily ambient environmental conditions during the measurements), the CO₂ concentration to 400 ppm, the photosynthetic photon flux density (PPFD) to 1500 μmol m^-2 s^-1 (to ensure saturating light conditions), and the air temperature inside the cuvette from 20 to 30°C depending of the sampling dates (to fit the mean midday air temperature during the measurements; Fig. S2). Each leaf was measured when the gas exchange values were stable (i.e., after max. 5 minutes).  On the same day as gas exchange measurements, one twig per tree was collected before sunrise (Ψ_pd) and at midday (Ψ_md) to measure the leaf water potential (MPa) with a Scholander-type pressure chamber (M1505D, PMS Instruments, USA). The difference between Ψ_pd and Ψ_md (Δ_Ψ) was calculated to describe the tree’s stomatal regulation strategy.

Water uptake patterns

On each sampling date, in the monospecific and four-species mixtures (n=15 plots), we collected three 10 cm-long twig samples across the canopy from each tree between 9 am and 3 pm (local time). After removing the bark, the samples were immediately sealed in airtight vials (Exetainer, Labco Limited, UK). The vial lid was wrapped with parafilm and placed in cool conditions to avoid evaporation. On the same day as the twig sampling, soil samples were collected every 10 cm at four depths (0–10, 10–20, 20–30, and 30-40 cm) and three random positions in each plot, using a manual soil corer or/and a pickaxe when the soil was too rocky. For each stand, the soil samples per depth were immediately pooled together in 50mL vials and stored like twigs. As the depth of the limestone bedrock varied within and between plots (between 20 and 70 cm; Table S1), the maximum depth of soil sample collection varied by date, plot, and position. Nevertheless, we could extract all the soil layers from 0 to 40 cm in every plot at each campaign, except for the Q.ilex monocultures in spring 2021, due to technical limitations. Precipitation water (used as a proxy of the groundwater isotopic values, see below) was collected by two Tube-dip-in-water collector types with pressure equilibration (RS1, Palmex, HR), spread into the study area (Fig. S1). Due to the unique design of the rain samplers avoiding evaporation for up to one year (Gröning et al., 2012), we collected the rainwater every two months during the growing season and once during the winter each year. The water was placed in vials and sealed with a lid and parafilm.

Water from xylem and soil samples was extracted using a custom-made cryogenic vacuum distillation system at the Swiss Federal Institute for Forest, Snow, and Landscape Research (WSL, Birmensdorf, CH) (Diao et al., 2022). The extraction system consisted of 20 tubes connected to 20 U-shaped collection tubes specifically designed for this system. A frozen sample was placed in the extraction tube and submerged in water at 80°C, while the associated collection tube was submerged in liquid nitrogen. The system was then evacuated to 5.10^-2 mBar. The extraction was maintained for 2 h for both xylem and soil samples to achieve a complete extraction following the recommendations of West et al. (2006) (i.e., a minimum of 60 min extraction time for a broad range of plant and soil materials). This process led to an extraction of 99.96% of the water in the samples, with more than 1 mL extracted for each sample, limiting the uncertainties in plant water isotopic composition due to cryogenic vacuum distillation (Diao et al., 2022). After the extraction, water samples were transferred into cap-crimp 2-ml vials and stored at −20 ^◦C until the isotopic analysis. Uncertainties associated with bulk water extractions using cryogenic distillation systems could occur that would underestimate the contributions of soil water and overestimate the ones from groundwater (Barbeta et al., 2021). Yet, as all samples were treated similarly, the errors would only affect the actual values, not the comparison between mixtures and seasons.

The δ²H and δ¹⁸O of all water samples (i.e., twig, soil, and precipitation) were measured with a high-temperature conversion elemental analyzer coupled to a DeltaPlus XP isotope ratio mass spectrometer (TC/EA-IRMS; Thermo, DE). Isotope ratios were reported in per mil (‰) relative to Vienna Standard Mean Ocean Water (VSMOW). Calibration versus the international standards was achieved by analysis of a range of certified water of different isotope ratios, resulting in a precision of 2‰ for δ²H and 0.3‰ for δ¹⁸O. The combined use of δ²H and δ¹⁸O provides a more comprehensive and nuanced understanding of the water movement through ecosystems as they can differently account for potential fractionation effects (Meißner et al., 2014).

Water source contribution

From the natural abundance of δ²H and δ¹⁸O in plant xylem and soil water, we used a Bayesian stable isotope mixing model to quantify the contribution of potential tree water sources for each species in the monospecific and four-species mixtures for each sampling date. As the rainwater isotope ratio differs throughout the season (i.e., isotopically more depleted rainwater in winter compared to summer) and the water evaporation decreases with soil depth, each water source has a significantly different stable isotopic composition (Fig. S3). These distinct soil isotopic profiles allow us to determine the contribution of each water source to the tree xylem water under the assumption that there is no isotopic fractionation during water uptake by the roots (Dawson & Ehleringer, 1991). Therefore, the natural isotopic abundance of xylem sap should reflect the water sources used by the plant. We used the package simmr in R (Parnell, 2019), where the isotopic composition (δ ¹⁸O and δ ²H) for each potential source (i.e., 0-10 cm, 10-20 cm, 20-30 cm, 30-40 cm, rainwater) and each target tree were assigned into the model. We set the TEF (trophic enrichment factor) and the concentration dependence to 0 due to the presumed absence of isotopic fractionation by the roots (but see Barbeta et al., 2020; Ellsworth & Williams, 2007). Using the isotopic values from each plant and the soil water source of the corresponding plot for each date, we ran the model where 3600 iterations out of 10000 runs were produced over 4 Markov chain Monte Carlo (MCMC) (Sun et al., 2022). To increase the clarity of presentation, the contributions from the water sources were grouped a posteriori into three layers: shallow (i.e., 0–20 cm), deep (i.e., 20–40 cm), and water stored in the fractured bedrock (i.e., rainwater). Indeed, the winter precipitation that penetrates deep soil layers and bedrock cracks could be a substantial water source for trees in Mediterranean forests during summer droughts (Eliades et al., 2018). As the bedrock water is not subjected to evaporation (Ehleringer & Dawson, 1992), we used the precipitation collected at our site as a proxy, similar to Grossiord et al., (2017). The cumulative rainwater collected during the winter of 2021 (i.e., January-April) and 2022 (i.e., October-April) was considered bedrock water source for the May campaigns of 2021 and 2022, respectively. During the remaining growing season, the stable isotope composition of the rainwater accumulated until the sampling date was added to the winter precipitation (e.g., values from winter, May, and June rainwater were used for the campaigns in July).

Water source partitioning

To estimate the plot-level vertical water source distribution in the four-species mixture plots, we calculated the belowground water source partitioning at the community level (P_W, unitless) from the sum of the differences in either the natural abundance of δ²H or δ ¹⁸O in plant water between all interacting species in every four-species mixtures for each sampling date. As different water uptake depths between trees correspond to greater differences in the xylem water stable isotope within a tree cluster, higher P_W indicates that the trees are taking up water from more distant water sources from each other. Therefore, to calculate P_W based on δ²H, we used the following equation developed by Grossiord et al. (2018):

P_W = |δ²H_PN - δ²H_PS| + |δ²H_PN - δ²H_QF| + |δ²H_PN - δ²H_QI| +

|δ²H_PS - δ²H_QF| + |δ²H_PS - δ²H_QI| + |δ²H_QF - δ²H_QI|                  (Equ. 1)

where δ²H_PN, δ²H_PS, δ²H_QF, and δ²H_QI correspond to δ²H of P. nigra, P. sylvestris, Q. faginea, and Q. ilex, respectively. For the monospecific plots, we applied the same formula but with the sum of the differences in δ²H or δ ¹⁸O in plant water between the tree individuals for each sampling date. P_W resulting from δ¹⁸O and δ²H weighted similarly (Fig. S4), supporting the absence of significant differences in using δ¹⁸O and δ²H to determine water source partitioning found by Wang et al. (2019). Hence, as δ²H had higher range of values increasing the vizualisation of the data, only the P_W based on δ²H was further used in this paper.

Statistical analyses

The effects of tree species diversity on Ψ_pd, Ψ_md, Δ_Ψ, A_net, g_s, xylem δ²H and δ¹⁸O were determined through linear mixed-effects models for each species using the package lmer. The effect of the season (i.e., spring, summer, fall), year (i.e., 2021, 2022), and species diversity (i.e. monospecific/monofunctional/ multifunctional/four-species mixtures) were used as fixed effects, and the individual plot and tree was treated as a random intercepts. Similar models were used to determine differences in soil water δ²H and δ¹⁸O. Sampling dates, species diversity, species, and soil depth were used as explanatory factors in the fixed part of the model. Significant differences in soil water δ²H and δ¹⁸O between depths for each species' diversity and sampling dates were found, allowing us to use the Bayesian isotope mixing model to determine the water source contribution of trees as described above (Fig. S3). The output of this model was analyzed similarly with linear mixed-effects models for each species. First, the effect of soil depth (i.e., shallow, deep, bedrock), season, year, and species diversity (i.e., monospecific and four-species mixture) were set as fixed effects, and the individual plot as a random effect.  Then, we ran similar models for each soil depth and species where the season, year, and species diversity were used as fixed effects and the plots as random effects. To reveal significant differences between species richness for each measurement at each sampling date and each species, post hoc analyses were performed with Tukey's HSD test, with FDR correction for multiple testing. Linear regressions were used to test the relationships between Δ_Ψ, Ψ_pd, g_s, xylem δ¹⁸O, P/PET, and P_W. All analyses were performed using the R v.4.2.2 statistical software (R Development Core Team, Vienna, Austria, 2022). Before performing each model, the homogeneity of variances and the normality of residuals were assessed, and data were log-transformed if necessary (Table S2).