1. Study site
The “Dispositif Instrumenté en Agroforesterie Méditerranéenne sous contrainte hydrique” (DIAMs) is located 10 km south of Montpellier (France) under a Mediterranean climate (43.612°N; 3.976°E). The soils are classified as Skeletic Rhodic Luvisols (IUSS Working group WRB 2014), according to the high proportions of stones (up to 60%), a pronounced red color and a layer deeper than 100 cm with accumulation of clay. The layers were identified as Ap1: plow layer with a lumpy texture (0–20 cm), Ap 2: second plow layer (20–50 cm), Bt 1: illuvial layer with lattice clays with stones (50–80 cm), Bt2: illuvial layer with lattice clays without stones (80–125 cm) and IICk: bedrock with calcium carbonates (Fig. 1b).
This agroforestry site of 5 ha is divided into three blocks separated by several hundred meters and representing independent replicates. The alley-cropping system was planted in 2017 with black locust (Robinia pseudoacacia), a N-fixing species, in 2-m-wide rows covered with an understory vegetation strip (UVS). The herbaceous species sown in 2019 (Festuca arundinacea, Dactylis glomerata, Medicago sativa, Trifolium pratense), however, were only scarcely present, and the dominant species was ray grass according to vegetation surveys conducted in 2020. The trees were planted every 2 m on the rows, leading to a density of 294 trees ha− 1, and were intercropped with 17 m crop alleys made of rotations of cereal and legume (wheat/barley/pea) crops. Beginning in October 2019, the soil was plowed to 20 to 50 cm. At the end of October 2019, durum wheat (Triticum turgidum L. subsp. durum) was sown as a winter crop with a sowing density of 140 kg ha− 1. A phytosanitary treatment (Atlantis Pro) was applied on January 30, 2020. Fertilization campaigns were conducted on the 06th of February 2020 with Nexens 46 fertilizer (130.4 kg N ha− 1) and on the 11th of April 2020 with Smart N 46 fertilizer (173.9 kg N ha− 1). Crop harvest occurred on June 25th, 2020, and straw was left on the topsoil.
To evaluate the homogeneity of N fertilizer spread, 3 pseudoreplicates x 3 blocks = 9 replicated cups with a diameter of 19.5 cm were set in the crop at 1.5 and at 4 m from the tree line. Immediately after fertilizer application, the granules were collected in the cup, weighed, and reported on a per m2 basis.
2. Sampling strategy
In May 2020, corresponding to the flowering stage of wheat, one pit (4 m long × 1.7 m wide × 1.5 m deep) was dug per plot (3 replicates) with a backhoe. To focus on the interaction between the UVS and wheat in the transition zone, the pits were dug between two trees (spaced by 2 m) instead of at the base of a tree, perpendicular to the herbaceous strip and starting from 0.5 m from the middle of the tree row to the intercrop (Fig. 1a). All pits were oriented toward the south, i.e., the sunny side of the trees. In addition to the vertical axis with 3 soil depths (0–20, 20–50 and 50–100 cm of depth) established according to the identified soil layers (Fig. 1b), the sampling strategy also considered a horizontal axis with three locations depending on the distance to the tree (Fig. 1a): the UVS (between 0.5 and 1 m from the tree line), the Crop-1 m (between 1 and 2 m from the tree line) and the Crop-4 m (between 3.5 and 4.5 m from the tree line). In summary, the sampling strategy was defined by 3 independent plots × 3 locations according to the distance from the tree line × 3 soil depths.
Just before digging the pits, aboveground biomass was sampled on each subplot and location in a 1 m2 area. Samples were air-dried at 50°C for 72 hours and weighed, and values were reported on a per m2 basis.
In this study, few medium roots (diameters between 2 and 10 mm) were found and identified as being from the UVS plants. Because they were so rare (less than 4 roots considering all pits), they were not considered in this study, which focused on fine roots (smaller than 2 mm in diameter).
3. Root characterization
Root biomass
During the excavation of the pits, a backhoe bucket was used to sample a volume of 28 dm3 of soil at each location and depth. Living roots were sorted manually from this volume, carefully washed on a 0.5 mm sieve, air-dried at 50°C for 72 hours and weighed. Knowing the volume of the backhoe bucket, the root biomass density (g dm− 3) was calculated for each location and soil depth (n = 3). Root dry biomass was reported on a per m2 basis (g m− 2), and the proportion of root biomass in each soil depth (%) was calculated as the ratio of the root biomass in the soil depth to the total root biomass in the entire sampled soil profile.
Root mapping
For each pit and at each location, 4 replicated grids of size 50 × 100 cm were hung against the refreshed pit soil profile. The grids were made of elementary tiles of 10 x 10 cm (5 in width and 10 in length). The number of root intersections in each elementary tile in the grids was manually counted. At each location, the root map could then be established from 0 to 100 cm of depth with a mean of the root impact density calculated on the 3 blocks and the 4 pseudoreplicates (n = 12).
The root intersect method from Chopart and Siband (1999) was applied to analyze the root isotropy (preferential orientation of the root direction of growth). Two stainless steel cubes of 10 x 10 x 10 cm (1 dm3) with opposite orientations were placed in 2 soil profiles (pseudoreplicates) of each pit and at each location. However, because the soil compaction at 50–100 cm was too high, it was not possible to use the cubes at that depth. The cubes were then extracted as intact as possible. On each of the 3 faces (H- perpendicular to the profile and parallel to the soil surface, T- parallel to the profile and perpendicular to the soil surface and L- perpendicular to the profile and to the soil surface), the root impact density was determined as the number of root intersections with the vertical plane; see Maurice et al. (2010) for details. The A coefficient was calculated as an indicator of the root isotropy (Chopart and Siband 1999):
$$A=\sqrt{\frac{{{(RID}_{T}-\stackrel{-}{RID})}^{2}+{{(RID}_{L}-\stackrel{-}{RID})}^{2}+{{(RID}_{H}-\stackrel{-}{RID})}^{2}}{6\times {\stackrel{-}{RID}}^{2}}}$$
where \({RID}_{i}\) is the root impact density on face i and \(\stackrel{-}{RID}\) is the average root impact density on all faces of the cube. When A = 0, there is root isotropy, no specific orientation for root growth. When A = 1, the roots grow only in a preferred direction.
In total, 72 cubes were sampled with 2 cube orientations × 2 soil depths × 3 locations × 2 pit soil profiles × 3 pits. After determining the root impact density on the 3 faces of the soil cubes, the roots inside the soil cubes were carefully sorted, washed with water on a 0.5 mm sieve and stored at 4°C prior to analyses.
Root functional traits
Roots extracted from the cube sampling were spread out in a transparent bin with deionized water and scanned at 300 dpi with a scanner (Epson Expression © 10000 XL). The resulting images were processed with image analysis software (WinRHIZO v. 2005b Regent, Canada©) for each cube content, allowing to determine the total root length (cm) and the mean root diameter (mm). Knowing the volume of the cubes, the root length density (cm dm− 3) was calculated. After the scans, roots from the cube sampling were air-dried at 50°C for 72 hours and weighed. The specific root length (m g− 1) was calculated as the ratio of fresh root length to root dry mass in each cube. Furthermore, the root biomass density obtained in the cubes at 0–20 and 20–50 cm compared the root biomass density obtained from the backhoe sampling detailed above: there was no significant difference between the two methods (p value = 0.12, data not shown). Only the results from the backhoe method were used for root biomass quantification.
Root chemical composition
Dry roots from the backhoe sampling were used for the biochemical analyses because of the higher root biomass retrieved. C fractions (soluble compounds, cellulose, hemicellulose and lignin) were determined with a fiber analyzer (Fibretherm®, Gerhardt) on a 500 mg root litter subsample following the protocol of Van Soest (Goering and Van Soest 1970). Root C and N contents were determined from 3 mg subsamples of root litter with an automatic elemental analyzer (Flash 2000, ThermoFischer Scientific). For total root P content, 50 mg of litter powder was mixed with 65% HNO3 and mineralized for 15 min at 200°C in a Milestones Ethos Easy microwave with a standard and blank. The total P content was then quantified colorimetrically with the yellow vanadomolybdate assay (NF U42-246). The C stocks in living roots (kg C m3 soil) were calculated by multiplying the root C content (%) by the root biomass density (kg m− 3). Several root quality indexes were calculated: the lignin:N ratio (Melillo et al. 1982), the soluble compounds:cell wall ratio (McClaugherty et al. 1985) and the ligno-cellulosic index (LCI,Melillo et al. 1989) calculated as:
$$LCI= \frac{lignin}{lignin+cellulose+hemicellulose}$$
where lignin, cellulose and hemicellulose are the corresponding root C fractions (%).
Root anatomy
Just after digging the pits, living wheat roots were sampled from the pit’s wall at each soil depth and at 2 locations (Crop-1 m and Crop-4 m) and immersed in a mixture of formaldehyde, acetic acid and
ethanol (Kladnik 2013) kept at 4°C pending further preparation for anatomical observation. The roots were then fixed in a 25% glutaraldehyde and 10% paraformaldehyde solution in a phosphate buffer (0.2 Na, pH 7) with 1% caffeine and stored in a 70% ethanol solution at 4°C (Salma 2015). Alcohol baths with increasing concentrations of absolute alcohol allowed gradual dehydration. Embedding of samples in resin was performed with a Technovit® 7100 kit. Samples were then cut into 4.5 µm sections with a manual Leica RM2255 microtome and mounted on a glass slide. Sections were then successively stained in periodic acid-Schiff’s reagent (to reveal insoluble carbohydrates) and in naphthol blue black to reveal proteins (Buffard-Morel et al. 1992). Screening of the glass slides was performed automatically with a Hamamatsu NDP slide scanner (Hamamatsu Nanozoomer 2.0HT, MRI). Several measurements and counts were performed on 1 section × 4 pseudoreplicated roots (for each location, soil depth and plot) with ImageJ Software (version 1.53e): total root area, stele area, number of metaxylem vessels, area of each metaxylem vessel and endoderm width. For each measurement, we calculated the mean value of the pseudoreplicated roots to analyze the root anatomy for each soil depth and each location (n = 3).
4. Soil characterization
During the excavation of the pits, more than 5 kg of soil was sampled at each location and soil depth (n = 3). The soil was immediately passed through a 2 mm sieve. Subsamples were stored at 4°C keeping their field moisture prior to analyses, other subsamples were frozen, and others were dried at 65°C for 72 hours. Moisture was determined for all samples after drying for 48 hours at 105°C.
Dry soils (dried at 65°C for 72 hours) were analyzed for total C and N contents by dry combustion (Matejovic 1997). Other dry soil subsamples were analyzed by the Laboratoire d’Analyse des sols (INRAE, Arras) for soil CaCO3 content (Allison 1960), cation exchange capacity by cobalt hexamine (Ciesielski and Sterckeman 1997), soil pH in a 1:5 soil-water suspension, total P content (Ciesielski et al. 1997; Ivanov et al. 2010), available P content (Olsen 1954), organic P content and soil texture using five fractions: clay, silt (fine and coarse) and sand (fine and coarse). The organic C content was calculated as the difference between the total C content and inorganic C content as measured in CaCO3.
On wet soils, mineral N was extracted with a 1:4 soil-1 M KCl solution. NO3− and NH4+ were determined by continuous flow colorimetry (Continuous Flow Analyzer, Skalar), and the sum of NO3− and NH4+ represented the mineral soil N content. Dissolved organic C was extracted with a 1:5 soil-distilled water solution, and the dissolved organic C was measured with a TOC/TN analyzer (TOC-Vsch-TNM SHIMADZU). The soil MBC (microbial biomass C) and MBN (microbial biomass N) were quantified with the chloroform fumigation-extraction technique (Vance et al. 1987). A correction factor of 0.45 was included.
In frozen soils, the abundance of the 16S and 18S rRNA genes was estimated by real-time quantitative polymerase chain reaction (q-PCR) of ribosomal DNA (Smith and Osborn 2009) using a thermal cycler (BioRad, Hercules, California). The potential activities of 5 extracellular enzymes involved in the C, N and P cycles were measured according to the protocol of Bell et al. (2013): β-1,4-glucosidase (BG), cellulobiohydrolase (CBH), N-acetyl-glucosaminidase (NAG), leucine aminopeptidase (LAP) and alkaline phosphatase (AP). The enzymatic activities were obtained in nmol g− 1soil min− 1 with a fluorometric microplate reader (Victor 3, Perkin Elmer, 365 nm excitation and 450 nm emission).
The soil bulk density was measured in the topsoil (0–20 cm) using a water method as it is suitable for stony soils (Jolivet et al. 2018). An excavation was made from 0 to 20 cm. The soil was passed through a 2 mm sieve, dried at 105°C for 48 hours and weighed. The volume was determined by putting a plastic bag in the excavation hole and measuring the volume of water that could fit in it. At deeper soil depths (20–50 and 50–100 cm), an imaging-sensor-based measurement method was used as a convenient method for stony soils (Coulouma et al. 2021). An excavation was made in the wall of the pit at depths of 20–50 and 50–100 cm. The soil was passed through a 2 mm sieve, dried at 105°C for 48 hours and weighed. The volume of the excavation was determined by three-dimensional imaging. For all the samples from depths of 0–20, 20–50 and 50–100 cm, the gross fraction (> 2 mm) was washed with water, dried and weighed to obtain the mass of stones. The stone volume was measured on several stones immerged in a burette, and the stone density was calculated as the ratio of the mass of stones to the stone volume (2.13 g cm− 3).
8. Data analyses
In the cube sampling, the 2 cubes with opposite orientations taken next to each other were spatially correlated for each soil depth and location. Furthermore, all cubes coming from the same soil profile, location and depth were considered as pseudoreplicates. Consequently, we considered them all as one sample and used the mean value (n = 3).
Using the linear regression coefficient of the relationship between the root length density and the mean root impact density obtained in each cube, we converted the root impact density map into a root length density map according to Maurice et al. (2010).
The C stocks in the living roots were calculated by multiplying the root biomass density by the root C content for each location and depth (n = 3). For the crop (at the Crop-1 m and Crop-4 m locations), all roots entered the soil after harvest, and 100% of the C stocks in living roots were thus considered as annual root C inputs to the soil. In contrast, for the perennial herbaceous species in the UVS, the annual root C inputs to soil were considered to be 53% of the C stocks in living roots using a root turnover rate of 53% annually in the grasslands (Gill and Jackson 2000). This estimation method is conservative, as it did not take into account the release of C through rhizodeposition, which could yet represent an important part of additional annual C inputs, as C in the rhizodeposits representing approx. 3% of the plant gross primary production (Pausch and Kuzyakov 2018).
The C-enzymatic activity was calculated as the sum of BG and CBH. The N-enzymatic activity was calculated as the sum of NAG and LAP. The P-enzymatic activity was that of AP. The C-, N- and P-enzymatic activities were integrated over time (180 min) to obtain cumulative enzymatic activities (µmol g− 1soil) or divided by MBC to obtain biomass-specific enzymatic activities (µmol g− 1MBC min− 1). An eco-enzymatic stoichiometry analysis was conducted as in Fanin et al. (2016). In short, the relative proportion of C versus P acquiring enzymatic activities (C/[C + P]) and the relative proportion of C versus N acquiring activities (C/[C + N]) were calculated (x and y respectively). Then, the length of the vector quantifying relative C versus nutrient acquisition was calculated as the square root of the sum of squared values of x and y [Length = Sqrt(x2 + y2)]. The angle of the vector quantifying the relative P versus N acquisition was calculated as the arctangent of the point (x, y) [Angle (degrees) = Degrees(Atan2(x, y))].
The soil bulk density was calculated as follows:
$${BD}_{soil}=\frac{{M}_{soil}}{{V}_{excavation}}$$
where \({BD}_{soil}\) is the soil bulk density (g/cm3), \({M}_{soil}\) is the total dry mass of the sample (g) and \({V}_{excavation}\) is the volume of the excavation (cm3).
The fine soil bulk density was calculated as:
$${BD}_{fine soil}=\frac{{M}_{fine soil}}{({V}_{excavation}-\frac{{M}_{stones}}{{D}_{stones}})}$$
where \({BD}_{fine soil}\) is the fine soil bulk density (g/cm3), \({M}_{fine soil}\) is the dry mass of the fine fraction < 2 mm (g), \({V}_{excavation}\) is the volume of the excavation (cm3), \({M}_{stones}\) is the mass of stones in the sample (g) and \({D}_{stones}\) is the stone density (g cm− 3).
The soil organic C stocks were calculated at each plot, location and soil depth following the ‘M4’ method described by Poeplau et al. (2017) as recommended for stony soils:
$${C}_{stock\_i}={C}_{org}\times {BD}_{fine soil}\times (1-Stones)$$
where \({C}_{stock\_i}\) is the soil organic C stock at point \(i\) (specific location × depth × plot) (kgC msoil−3), \({C}_{org}\) is the organic C content in the fine soil measured in \(i\) (mgC gsoil−1), \({BD}_{fine soil}\) is the bulk density of the fine soil in \(i\) (g cm− 3) and \(Stones\) is the volumetric fraction of stones (%).
To analyze the effect of location as a fixed factor and the three replicated profiles as random factors on aerial and total root biomass and on the mass of spread N fertilizer, linear mixed models were fitted for each variable. For variables analyzed at different depths, the effect of soil depth, location and their interactions as fixed factors and the three replicated profiles as random factors were tested using linear mixed models: soil physical, chemical and microbiological properties, root traits (biomass density, functional traits, chemical composition, quality index and anatomy), soil C stocks and living root C stocks. Post hoc Tukey tests were used to assess differences between soil depths and locations.
In addition to the ANOVA-like mixed effect model, a complementary set of covariance analyses (ANCOVA) mixed effect models were used to test our hypotheses. To investigate the correlation between soil physico-chemical properties and root variables along the vertical gradient and if these correlations changed according to the 3 locations in the agroforestry system (Hypothesis 1), we used a model comparison approach. Only the most informative root variables associated with root quantity (biomass density), root quality (lignin content) and root anatomy (stele diameter) were investigated. The selection of the soil physico-chemicals used as covariates in these analyses was based on the PCA (Fig. 7b) so that they represented the main axes of variations: soil bulk density, mineral N and Olsen P contents (axis 1, Fig. 7b) and the stone volumetric content (axis 2, Fig. 7b). First, to test the interaction (if the effect of the soil physico-chemical covariate on root variables varies according to the location), we compared a model including the soil physico-chemical covariate, the location and their interaction (lm1) with the same model but excluding the interaction (lm2) using the anova() function in R. Second, to test if the soil physico-chemical covariate explained a part of the root variable that was not explained by the location (the vertical variation of root variables), we compared the model with both the soil physico-chemicals and the location (lm2) with a model with the location only (lm3).
We used a second set of model comparisons to investigate the correlation of root properties with microbial activity and soil C stocks along the horizontal gradient of the agroforestry system and whether these correlations changed according to the 3 depths investigated (Hypothesis 2). We selected the most informative root variables as covariates in these analyses to reveal the potential contribution of roots to soil C stocks: the living root C stocks, the root quality (informing the root decomposability): the lignin:N and C:N ratios and the ligno-cellulosic index. First, we tested if the effect of root covariate on microbial activity or soil C stock variables changed according to depth by comparing a model including the root covariate, the depth and their interaction (lm1) with the same model but excluding the interaction (lm2) using the anova() function in R. Second, we tested whether the root covariate explained a part of the soil variable that was not explained by the depth (the horizontal variation of soil variables). We compared the model with both the root covariate and the depth (lm2) with a model with the depth only (lm3).
For all the linear mixed models and analyses of variance, the lme4 and car packages were used. The normality of the residues was always verified with a Shapiro–Wilk test, and the homogeneity of the variances was verified with a Bartlett test. When necessary (p values < 5%), logarithmic, square root, Box–Cox or Yeo Johnson transformations were applied. Simple ordinations of the root variables and the soil physico-chemical and microbiological properties were conducted using principal component analyses (PCAs) with the vegan and factoextra packages.
All statistical analyses were performed with R software (version 4.0.5). For each measurement, data are presented as the mean values ± standard deviation of three independent replicates.