Faunal communities mediate the effects of plant richness, drought, and invasion on ecosystem multifunctional stability
DesignPlant richness. Sixteen locally frequent native plant species in the barren mountain areas (around Taizhou University, Zhejiang, China) invaded by the exotic plant Symphyotrichum subulatum60 were selected as the native species pool. These species were chosen because they spanned the dicotyledon plant taxonomy (including 7 Orders, 10 Families, and 14 Genus, in the Class Magnoliopsida), differed widely in their functional traits (related to height, life form, dominance in local communities, and leaf habit) (Supplementary Table 3), and were occasionally found to be associated with the invasive species Symphyotrichum subulatum60 in the local secondary-succession communities. With this species pool, we were able to imitate the locally natural, spatially stochastic, compositionally ruderal, and functionally varied plant community61, which is a typical attribute of the secondary-succession communities in the local barren mountains invaded by the exotic plant Symphyotrichum subulatum. Based on this native species pool, monocultures of each species (16 total), and random mixtures of 2, 4 or 8 species (with 10, 10, or 9 distinct assemblages, respectively) were designed, creating a complete set (Fig. 1d) of 45 different plant assemblages (pots) in total. Each plant assemblage was replicated 6 times, for a total of 270 pots. To eliminate the non-random effects during the 1-year development of the 270 pots, their distributions were randomized, such that not all replicates of an assemblage were next to each other (Fig. 1d–f).DroughtAfter 1-year development of the native plant assemblages, three drought treatments (non-, moderate-, and intensive-drought) were manipulated by adjusting irrigation using automatic drip irrigation systems, with 100%, 50%, and 25% of the equivalent to the amount received in the areas where native species were collected, respectively. Two random complete sets were selected for each drought treatment, each complete set being composed of 45 different plant assemblages (Fig. 1d–f).Exotic plant invasionNine months after drought treatment, the two complete sets (Fig. 1d) of each drought treatment were randomly exposed (invasion) or not exposed to (non-invasion) the invasive species Symphyotrichum subulatum (Michx.) G. L. Nesom (Fig. 1e, f). S. subulatum, an annual herbaceous plant native to North America, is a common invasive species in the subtropical and tropical regions of China18,60, and tends to interact with the native species via, for example, competing for space and resources62,63, enriching for pathogens or herbivores, and changing soil faunal, bacterial or fungal microbiomes18,64,65.ExperimentThe experiment based on the design mentioned above was conducted at Taizhou University, Zhejiang province, China (28.66°N, 121.39°E). The seeds of the 16 native plant species (Supplementary Table 3) and the soil were collected from nearby mountain areas (Wugui, 28.65°N, 121.38°E; Baiyun, 28.67°N, 121.42°E; Beigu, 28.86°N, 121.11°E). The seed-mixtures were obtained by mixing seeds of the 16 species pro rata, in proportion to germination rates. The soil (fine-loamy, mixed, semiative, mosic, Humic Hapludults) was sieved to pass a 2-mm mesh, and thoroughly mixed. 270 plastic pots (72 cm length × 64 cm width × 42 cm depth) were prepared, and each was filled with a 27-cm soil layer, followed by a 10-cm mixture of soil and vermiculite-compost to provide water-, air- and fertility-support for germination, seedling establishment, and plant growth (Supplementary Table 4).Native plant assemblagesAll the 270 pots were placed inside a plastic shelter, which allowed for both air ventilation and protection from rain. Each pot was sown with a seed-mixture of ca. 800 seeds. One month after germination, for each pot, the undesired seedlings were removed manually according to the plant richness design (Fig. 1d–f), and thus 32 vigorous seedlings (with the same number of seedlings per species, e.g., 4 seedlings for each species of the 8-species mixtures) were spatial-evenly retained. In this manner, the plant richness was manipulated for each plant assemblage. During the development of the 270 plant assemblages, the soil volumetric water content was controlled at ca. 20%, which was similar to that of the nearby mountainous soil, using the automatic drip irrigation systems. Weeds and undesired species were removed monthly (Fig. 1f).Drought treatmentAfter 1-year development of native plant assemblages, the drought treatments (non-, moderate-, and intensive-drought) were manipulated according to the experimental design mentioned above (Fig. 1d, e). Two complete sets (Fig. 1d) of different plant assemblages (2 × 45 pots) were selected for each drought treatment. Every other week, 40 pots each drought treatment were randomly selected for measuring soil water content and soil temperature at the depth of 0–20 cm, using the ProCheck analyzer (Decagon, Pullman, Washington, USA), and irrigation was adjusted accordingly using automatic drip irrigation systems. The irrigation for non-, moderate-, or intensive-drought was adjusted to accomplish an irrigation level amounts to 100%, 50%, or 25% that of the mountain areas where seeds were collected. Because of the distinct seasonal temperature and evaporation conditions, the irrigation frequencies were approximately daily in May-September, every other day in March–April and October–December, and weekly in January–February. With this manipulation, the volumetric soil water contents of non-, moderate-, and intensive-drought were controlled within ranges of 13.8–23.4%, 6.8–13.7%, and 1.4–7.4%, respectively, throughout the manipulation of drought treatment (Fig. 1e, f). Eight months after drought introduction, fresh litter was collected form the two replicate pots of each drought treatment, and then oven-dried at 40 °C, cut into ca. 2-cm pieces, and filled into litterbags (2-g litter in each litterbag).Invasion treatmentNine months after drought introduction, one complete set (45 pots) of the plant assemblages (Fig. 1d) from each drought treatment, was chosen and exposed to invasion disturbance by sowing 50 seeds of S. subulatum in each pot, and the other was specified as the non-invasion treatment (Fig. 1e, f). The prepared litterbags were embedded under the litter-layer of each pot (5 litterbags in each pot), correspondingly.SamplingSix months after invasion introduction, one litterbag was collected for litter-fauna extraction. Nine months after invasion, five soil cores (20-cm depth) were collected with augers (6.4 cm in diameter) and mixed for extraction of soil-fauna, and measurement of soil property and enzyme activity (Fig. 1f). The aboveground biomass of both native and invasive plants in each pot was harvested, sorted to species, oven-dried to a constant mass at 80 °C, and weighed. The belowground plant biomass was also sampled, sorted to native and invasive groups, oven-dried, and weighed (Fig. 1f).Plant, litter-, and soil-faunal communitiesPlant communitySince exotic plant invasion was treated as a disturbance factor, the biomass of the invasive species S. subulatum was not included for analyses concerning plant community and ecosystem (multi)functionality. The aboveground biomasses of native plant species in each of the 270 pots were collected for plant community analysis.Litter- and soil-faunal communitiesOne litterbag or fifty grams of mixed-soil samples were used for litter- or soil-fauna extraction using a Tullgren funnel apparatus (dry funnel method)66. The obtained microarthropods were stored in 70% alcohol, identified with double-tube anatomical lens, and classified to Family level. For both litter and soil samples, the numbers (abundances) of all faunal taxa were counted for litter/soil-faunal community analysis.Phylogenetic information of plant, litter-, and soil-faunal communitiesSimilar procedures were used to construct the plant and faunal phylogenetic trees. First, protein sequences of 12 faunal mitochondrial coding genes and 16 plant plastid coding genes (Supplementary Data 1) were obtained by searching plant or faunal taxonomies from NCBI protein database (https://www.ncbi.nlm.nih.gov/protein/) with Edirect software (https://www.ncbi.nlm.nih.gov/books/NBK179288/). All available sequences at plant species level or faunal Family level were fetched. If unavailable, the missing sequences were sampled from plant genus or faunal Order level. Sequoiadendron giganteum and Echinococcus were specified as out-group references for plant and faunal trees, respectively. Then, the sequences of each plant or faunal taxon were clustered at 97% or 90% identity independently, and the centroids were used as representative markers. The markers were aligned with MUSCLE67, followed by concatenation. Finally, using MEGA X68, the maximum likelihood trees were constructed based on BioNJ initial trees69 and 500 bootstrap checking nodal support. The parameters for plant tree construction were specified as follow: 70% partial deletion (with 4824 positions retained) and the best-fit substitution model JTT + G + I + F70,71; parameters for faunal tree: 90% partial deletion (2778 positions) and LG + G + I + F model71,72. The Linux codes for processing the protein sequences were submitted to GitHub (https://github.com/YuanGe-Lab/JZW_2022/tree/main/linux)The plant and faunal taxonomies, representative markers, and marker accessions are provided as Supplementary Data 1.Ecosystem function-related variablesA total of 14 individual function-related variables were collected. These variables belonged to three functional groups: (1) biomass production, including aboveground and belowground biomass of native plants, light interception efficiency, litter-fauna abundance, and soil-fauna abundance; (2) soil properties, including contents of soil organic carbon, soil nitrogen, soil phosphorus, and GRSP (relating to soil physical properties and stocks of carbon and nutrient73); and (3) processes, including rate of litter decomposition, and activities of β-glucosidase, protease, nitrate reductase and dehydrogenase.Light interception efficiency, the fraction of incident photosynthetically active radiation (PAR) intercepted by each plant community canopy, was determined between 12:00 and 14:00 on clear days using LI-191R line PAR sensors (LI-COR Inc., NE, USA), and the mean of 4 measurements (monthly from May to August the third year; Fig. 1f) was used. Total soil organic carbon and nitrogen were measured with an elemental analyzer (vario Max; Elementar, Germany). Total soil phosphorus was determined using the molybdenum blue method with a UV–visible spectrophotometer (Shimadzu, Kyoto, Japan). GRSP was determined using the method described by Shen et al.18. Litter decomposition rate was assessed by embedding litterbags and fitting litter mass loss against decomposition time (Fig. 1f). Enzyme activities were analyzed by the spectrophotometric method using the substrates, p-Nitrophenyl-β-d-glucopyranoside (pNPG; for β-glucosidase), caseinate (protease), nitrate (nitrate reductase) and triphenyltetrazolium chloride (TTC; dehydrogenase)18.Quantifying community stability and multifunctional stabilityCommunity data was comprised of native plant biomasses or faunal abundances, and the associated phylogenetic information. Multifunctionality data was comprised of 14 function-related variables, each variable (V) being transformed (V’) using the formula ({V}^{{prime} }=frac{V-{{{{{rm{min }}}}}}left(Vright)}{{{{{{rm{sd}}}}}}left(Vright)}) to guarantee even contribution to global variance. We calculated community similarity (1 minus Weighted-UniFrac distance) and multifunctional similarity (1 minus Bray–Curtis distance), based on the community data and the multifunctionality data, respectively. The specific subsets of each symmetric similarity matrix were used to assess three different aspects of stability: (1) Invariability (against stochastic fluctuations), reflected as the pairwise similarities (1476 pairs) within treatment groups, at same plant richness*drought*invasion condition; (2) Drought resistance, the similarities (2148 pairs) between drought (moderate- and intensive-drought) and non-drought treatments, at same plant richness*invasion condition; and (3) Invasion resistance, the similarities (n = 1611 pairs) between invasion and non-invasion treatments, at same plant richness*drought condition (Supplementary Fig. 1).We also assessed the three aspects of stability of each individual function in a similar way, but by calculating the similarity using the formula ({{{{{{{mathrm{SIM}}}}}}}}_{{ij}}=1-frac{|{V}_{i}-{V}_{j}|}{{V}_{i}+{V}_{j}}) (Vi and Vj are ith and jth elements in a function vector; SIMij is the similarity between Vi and Vj).Statistics and reproducibilityPERMANOVA (10,000 randomizations) was conducted to test the influences of the manipulated factors on ecosystem multifunctionality or communities of plant, litter- and soil-fauna, using “vegan::adonis” in R74. Mantel test (10,000 randomizations; Spearman’s R) was conducted to test the community-community or the community-multifunctionality relationships, using “vegan::mantel” in R74.As each similarity-pair of each aspect of community or multifunctional stability mentioned above was in strict correspondence to single level of each manipulated factor (plant richness, drought, and invasion) (Supplementary Fig. 1), the direct/indirect effects of treatments on the community or multifunctional stability can be assessed using SEM. To test direct and indirect effects (by modulating community stability) of the manipulated factors on multifunctional stability, we built three SEMs (Fig. 1a–c) based on three different aspects of stability (i.e., invariability, drought resistance, and invasion resistance) under the conditions of corresponding parings of manipulated factors (Supplementary Fig. 1), with the LAVAAN package75. The standardized paths (direct effects) in SEMs can be conceived as the partial correlations after teasing all side effects away. Bootstrapping with 10,000 randomizations was conducted to generate the unbiased mean effect. The significance of effect was tested using a Mantel-like permutation (10,000 randomizations) test76, where the null hypotheses (H0) were that the independent factors plant richness, drought, and invasion, had no direct/indirect effects (effect = 0) on multifunctional stability. Based on H0, permutation procedure was conducted by permuting the index of dependent factors (both columns and rows of a symmetric matrix; Supplementary Fig. 1) simultaneously to gain null models and null effects. p-values (probability of H0 acceptance) were calculated as the percentage of observed positive (or negative) effect that was greater (or less) than the null effects. We also assessed the direct and indirect effects of factors on the stability of each individual function based on the same SEMs, to consolidate our findings on multifunctional stability. The R codes and examples solving the permutation test for the significance of effects derived from SEMs that based on multidimensional similarity (or distance) were submitted to GitHub (https://github.com/YuanGe-Lab/JZW_2022/tree/main/R). All the analyses were conducted using R (https://www.r-project.org).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
