Changes in plant community composition can have long-lasting consequences for ecosystem functioning. However, how the duration of plant growth of functionally distinct grassland plant communities influences abiotic and biotic soil properties and thus ecosystem functions is poorly known. In a field experiment, we established identical experimental subplots in two successive years comprising of fast- or slow-growing grass and forb community mixtures with different forb:grass ratios. After one and two years of plant growth, we measured above- and belowground biomass, soil abiotic characteristics (pH, organic matter, soil nutrients), soil microbial properties (respiration, biomass, community composition), and nematode abundance. Fast- and slow-growing plant communities did not differ in above- and belowground biomass. However, fast-and slow-growing plant communities created distinct soil bacterial communities, whereas soil fungal communities differed most in 100% forb communities compared to other forb:grass ratio mixtures. Moreover, soil nitrate availability was higher after two years of plant growth, whereas the opposite was true for soil ammonium concentrations. Furthermore, total nematodes and especially bacterial-feeding nematodes were more abundant after two years of plant growth. Our results show that plant community composition is a driving factor in soil microbial community assembly and that the duration of plant growth plays a crucial role in the establishment of plant community and functional group composition effects on abiotic and biotic soil ecosystem functioning under natural field conditions.

The experimental site was established in spring 2015 in a restored grassland site (abandoned from agricultural use in 1996), “De Mossel” (Natuurmonumenten, Ede, The Netherlands, 52° 04´ N, 5° 45´ E). The area around “De Mossel” is characterized by a mean daily temperature of 16.7 °C in summer months and 1.7 °C in winter months and monthly precipitation ranges from 48 to 76 mm (based on open source data from long-term climate models; www.climate-data.org). The soil of the field site is described as a holtpodzol and soil texture is characterized as sandy loam (94% sand, 4% silt, 2% clay, ~5% organic matter, 5.2 pH, 2.5 mg kg^-1 N, 4.0 mg kg^-1 P, 16.5 mg kg^-1 K) (Jeffery et al. 2017).

In total, 100 experimental plots (1.66 x 2.50 m) were installed and each plot was divided into two subplots (each 0.83 ´ 2.50 m) resulting in 200 subplots (De Long et al. 2019). Experimental plots were randomly arranged in four blocks. To study the effects of the duration of plant growth of distinct plant communities, the top-soil (about 4 cm depth) including the previously existing plant communities of 100 subplots was removed in May 2015 (later referred to as “two-year subplots”) whereupon experimental plant communities were sown immediately. In May 2016 (later referred to as “one-year subplots”) the top-soil including the previously existing plant communities was removed of the remaining 100 subplots and experimental plant communities were then sown. In total, 24 predominantly perennial plant communities based on their economic spectrum (fast- versus slow-growing species) and differences in forb:grass ratio were chosen. Fast- or slow-growing grasses or forbs were selected from a pool of 24 grassland species that all co-occur locally at this site. Plants were assigned to fast- versus slow-growing species based on known growth rates (Fitter and Peat 1994, Fry et al. 2014) or after consultation with botanists (Jasper van Ruijven, Henrik Poorter, personal communication). Three fast- (Fast 1, Fast 2, Fast 3) and three slow- (Slow 1, Slow 2, Slow 3) growing plant communities were used in the experiment, each mixture consisted of 3 or 6 plant species (three grasses and/or three forbs) (Table 1). Specifically, sown plant communities differed in forb:grass ratios as per seeding such as a) three fast- or slow-growing forb species (100%); b) three fast- or slow-growing grass species (100%); or c) three fast- or slow-growing forbs and three fast- or slow-growing grasses (25% grass and 75% forb or 75% forb and 25% grass). Sowing density amounted to 12,000 seeds per subplot, representing each plant species in equal amounts of seeds. Seeds were obtained from specialized suppliers that provide seeds collected from wild plants (Cruydt Hoeck, Nijeberkoop, The Netherlands and MediGran, Hoorn, The Netherlands) in 2015. Additionally, after top-soil removal, in each block two experimental subplots were kept bare one from May 2015 onwards, and one from May 2016. These plots served as a control to permit the comparison of ecosystem processes in the absence of vegetation. In 2015 and 2016, during the growth season (May through September) all sown subplots and bare control subplots were regularly weeded to maintain the sown plant community composition. In total, this resulted in 2 temporal treatments (two- and one-year; represented by the level of subplots) and 25 plant community treatment combinations (2 community growth rates (fast, slow) ´ 4 forb:grass ratios (0:100; 25:75; 75:25 or 100:0 % forb:grass; represented by the level of plots) ´ 3 species combinations + 1 bare control), which were replicated across four blocks (200 subplots). The efficiency of the establishment of experimental plant communities was reported previously in De Long et al. (2019). Here, it was shown that the six different plant communities (Fast 1, Fast 2, Fast 3, Slow 1, Slow 2, Slow 3) significantly differed in their composition whereas fast-growing plant communities and slow-growing plant communities were clustering separately. Furthermore, it was reported that the actual percentage cover of both forbs and grasses corresponded well to the treatments.

Data collection

Above- and belowground biomass - In the beginning of June 2017, aboveground biomass production was assessed within two randomly selected squares (25 x 25 cm; minimum 10 cm distant to the edge) in each subplot by cutting plants just above the soil level. Aboveground biomass of all samples was determined by weighing after drying at 70°C for 72 h and later converted to g per m². After the aboveground biomass was clipped in a subplot, one soil core (diameter 3 cm, 10 cm deep) was taken from the center of the square. The fresh soil weight was determined, and the roots were then washed over a 0.425 mm sieve. Root biomass was then determined after drying at 40 °C for 72 h and later converted to gram per gram fresh weight of soil.

Soil sampling - Soil sampling was carried out in two sampling campaigns. First, soil samples were collected in March 2017. Here, 20 soil samples (metal corer: diameter 1 cm, 10 cm deep) were taken randomly per subplot, pooled and homogenized. Approximately 2 g of fresh soil was immediately frozen at -80 °C for DNA-based determination of the soil microbial community composition. In mid May 2017, an additional 20 soil samples (metal corer: diameter 3 cm, 10 cm deep) were taken randomly per subplot. The soil samples were pooled, homogenized and divided into several subsamples. About 50 g of fresh soil for nematode extraction was kept in a fridge (4 °C). Approximately 20 g was sieved (2 mm mesh size) to remove stones, roots, and invertebrates >2 mm, and then stored at -20 °C for measurements of soil microbial properties (soil microbial respiration and biomass), the remaining soil was dried at 40 °C and used to measure soil abiotic characteristics.

Soil abiotic characteristics - The soil was dried at 40 °C until the soil weight was stable and sieved (mesh size: 1.4 mm) to remove roots and stones. Three grams of dried soil were mixed with 30 ml of 0.01 M CaCl₂ and shaken for 2 h on a mechanical shaker with linear movement at 250 rpm. Samples were then centrifuged for 5 minutes at 3000 rpm. Then, 15 mL of the supernatant was filtered through a Whatman Puradisc Aqua 30 syringe filter with a cellulose acetate membrane and 130 μL HNO₃ was added to 12.87 mL of the filtrate. Soil extractable elements (Fe, K, Mg, P, S, Zn) were analyzed using an inductively coupled plasma - optical emission spectrometer (ICP-OES, Thermo Scientific iCAP 6500 Duo Instrument with axial and radial view and CID detector microwave digestion system). The remaining filtrate (2.13 mL) was used to measure soil pH and to measure nitrite (NO₂-N) + nitrate (NO₃-N) and ammonium (NH₄-N) on a QuAAtro Autoanalyzer (Seal analytical).

Soil organic matter - Soil organic matter content was estimated by the loss‐on‐ignition (LOI) method (Heiri et al. 2001). Approximately 5 g of soil was oven dried at 105 °C for 16 h and weighed. The sample was then burned at 550 °C for 5 h and weighed again. Soil organic matter was calculated as the percentage weight loss between the oven-dried and burned samples.

Soil microbial respiration and biomass - Approximately 5 g soil (fresh weight) was weighed into 50 mL centrifuge tubes to determine soil microbial respiration and soil microbial biomass. The lid of each tube was sealed gas-tight using an O-ring and a rubber septum in the middle. For basal respiration measurements, the tubes were capped and flushed with CO₂-free air to remove any CO₂ from the headspace. After 24 h hours of incubation at 20 °C, 12 ml of headspace was sampled using a gas-tight syringe. Microbial biomass was determined after addition of D-glucose-monohydrate using the substrate-induced respiration method (SIR) (Anderson and Domsch, 1978). Next, 2 ml of 75 mM D-glucose solution were added to each soil sample and placed on a horizontal shaker for 1 h. Tubes were capped, flushed with CO₂-free air and incubated for 4 h at 20 °C. Again, 12 ml of headspace was sampled. Measurements of the CO₂ concentrations were carried out on a Trace CG Ultra gas chromatograph (Thermo Fisher Scientific, Milan, Italy). Gravimetric soil water content was determined by drying the soil samples overnight at 60 °C to constant weight and calculating the difference in weight between fresh and dried soil.

Soil microbial community composition - DNA was extracted from 0.75 g of soil using the PowerSoil DNA Isolation Kit (Mo Bio Laboratories, Carlsbad, CA, USA) following the manufacturer’s protocol. The DNA quantity was measured using a Nanodrop spectrophotometer (Thermo Scientific, Hudson, NH, USA). Approximately 100 ng of DNA was used for a PCR. We used the primers ITS4ngs and ITS3mix targeting the ITS2 region of fungal genes (Tedersoo et al. 2015) and the primers 515FB and 806RB (Caporaso et al. 2012, Apprill et al. 2015, Parada et al. 2016) targeting the V4 region of the 16Sr RNA gene in bacteria. Presence of PCR product was checked using agarose gel electrophoresis. The PCR products were purified using Agencourt AMPure XP magnetic beads (Beckman Coulter). Adapters and barcodes were added to samples using Nextera XT DNA library preparation kit set A (Illumina, San Diego, CA, USA). The final PCR product was purified again with AMPure beads, checked using agarose gel electrophoresis and quantified using a Nanodrop spectrophotometer before equimolar pooling. We pooled all fungal samples (200) in one Miseq PE250 run and divided the bacterial samples in two separate runs (100 samples each; 1^st run: block 1 and 2, 2^nd run: block 3 and 4). Libraries were sequenced at McGill University and Genome Quebec Innovation Center, Canada. Extraction negatives were also sequenced. A mock community, containing 10 fungal species, was included to investigate the accuracy of the bioinformatics analysis.

Bacterial sequences and fungal sequences were analyzed using the PIPITS pipeline and the Hydra pipeline, respectively (Gweon et al. 2015, Hollander 2017). In short, fungal sequences were paired using VSEARCH and quality was filtered using standard parameters. The ITS2 region was extracted using ITSx (Bengtsson-Palme et al. 2013). Short reads were removed, and sequences were clustered based on a 97% similarity threshold using VSEARCH and chimeric sequences were removed by comparing with UNITE uchime database. The representative sequences were identified using the RDP classifier against the UNITE database (Kõljalg et al. 2005). For bacterial sequences VSEARCH was used to pair sequences and cluster them or classification, SINA classification was used with the SILVA database.

Nematode extraction - Decantation (Cobb 1918) and centrifugal flotation methods (Van Bezooijen 2006) were used to extract nematodes from 50 g fresh soil collected from block 1 and 2 (100 subplots). Briefly, soil samples were weighed, suspended in 3 L water and stirred until a homogenous suspension was obtained. The suspension was settled for 15 sec and the supernatant decanted into a plastic bowl. This procedure was repeated three times and the obtained suspension was passed through one 75 µm and three 45 µm sieves and collected into a 50 ml centrifuge tube. The nematode suspension was centrifuged for 5 minutes at 3000 rpm and the supernatant was carefully poured off. Sugar solution (484 g/L) was added, fully mixed with the remaining sediment in the tube and centrifuged again for 5 min at 1000 rpm. The supernatant was poured over a 30 µm sieve. Nematodes on the sieve were carefully rinsed with water and collected in a beaker for microscopic analysis). All nematodes were counted and classified by microscopy to one of five feeding types (root feeders, fungivores, omnivores, bacterivores, or predators) according to Yeates et al. (1993). However, predatory nematodes were excluded from further analysis because of low abundance.

Data analysis

We used linear mixed-effects models to test the effects of duration of plant growth (one- and two-year subplots), community growth rates (fast vs. slow), and forb:grass ratios on above- and root biomass, soil abiotic characteristics (pH, Fe, K, Mg, P, S, Zn, NO₂-N +NO₃-N, NH₄-N), soil organic matter, soil microbial properties, and feeding types of nematodes. Thereby, plot and plant community identity (Fast 1, Fast 2, Fast 3, Slow 1, Slow 2, Slow 3) were included as random factors. Above- and belowground biomass, chemical abiotic characteristics and soil microbial properties (respiration and biomass) data were ln- transformed to meet assumptions of ANOVA. Linear mixed- effects models were performed using lme4 package (Bates et al. 2015), whereas p-values and degrees of freedom were estimated with Type III Kenward-Roger approximation using lmerTest (Kuznetsova et al. 2017). We further used comparisons of means for treatment- specific effects (Tukey’s HSD test; α < 0.05). Tukey’s tests were performed using and multcomp package (Hothorn et al. 2016). Bare plots were excluded from analysis, but their values are displayed in graphs.

To test the effects of duration of plant growth, community growth rates, and forb:grass ratios on bacterial and fungal community compositions, we first filtered out taxa that were present in less than ten samples and had an abundance of less than 0.01%. ITS sequences derived from other organisms than fungi were further removed and for 16S rRNA data mitochondria and chloroplast sequences were removed. For both bacteria and fungi, samples with less than 1000 reads of more than 80 000 reads remaining were removed from the dataset and read numbers were further normalized using total sum scaling (TSS). Mock communities consisting of 10 fungal species were used to inspect the filtering done for fungi. After filtering, we detected 13 fungal OTUs which shows that we might be slightly overestimating the diversity. Afterwards, we ran permutational multivariate analysis of variance (PERMANOVA) (based on Bray–Curtis dissimilarities, 999 permutations) on both bacterial and fungal OTUs using the “adonis” function in the vegan package (Oksanen et al. 2016). However, here we only used one- and two-year subplots to test for the differences in plant community effect on soil microbial community composition. Since the basic assumption of PERMANOVA analysis is a balanced design and homogenization of samples, bare subplots were excluded due to lower number of samples than the one-and two-year subplots. For visualization we applied a nonmetric multidimensional (NMDS) analysis of the dissimilarities (based on Bray–Curtis dissimilarities) in microbial community composition using ggplot2 package (Wickham 2016). All statistics were performed within the R statistical environment (Version 3.5.1; R Core Team, 2018).