DOI: 10.5061/dryad.3r2280gp2
title: Herbivory and allelopathy contribute jointly to the diversity-invasibility relationship

description: To construct experimental communities of different species richness (1, 2,
 4 and 8 species), we used eight native species (Arctium lappa L., Patrinia
 villosa (Thunb.) Juss., Achyranthes bidentata Bl., Reynoutria japonica
 Houtt., Taraxacum mongolicum Hand.-Mazz., Plantago asiatica Ledeb.,
 Platycodon grandiflorus (Jacq.) A. DC., Antenoron filiforme (Thunb.) Rob.
 et Vaut.). These species were chosen because they are common in natural
 shrub-grassland communities in the Zhejiang province, China, and two of
 them, P. villosa and A. bidentata, are usually very dominant in those
 communities. Solidago canadensis L. was selected as the invader. In May
 2017, we prepared 188 plots of 2 m × 2 m in a previous agricultural field
 that had been abandoned seven years earlier, in Luqiao (121°23′46.16″ E,
 28°31′21.62″ N). During October to December 2016, the top 10 cm of the
 soil was removed to deplete the seed bank, and the remaining soil was
 ploughed and repeatedly rotovated to produce a fine soil. In March to
 April 2017, the plots were weeded three times to further deplete the seed
 bank. The experimental area was divided into four blocks of 47 plots. In
 each plot, we constructed monocultures for each of the eight species with
 one replicate each, six two-species mixtures of different compositions and
 with three replicates each, six four-species mixtures with different
 compositions and three replicates each, and one eight-species mixture with
 three replicates. Each native species was used in an equal number of plots
 at each species-richness level. When the invasion treatment started, we
 also established three monoculture plots of the invader (S. canadensis) to
 assess potential side effects of activated carbon application
description: <em>Treatments</em>
 In the herbivory-reduction
 treatment, we sprayed the plots with an insecticide. To avoid that the
 insecticide spray would drift into plots of the non-insecticide treatment,
 we randomly assigned the four blocks, each with 47 plots, to one of the
 four experimental treatments (i.e. herbivory reduction, allelopathy
 reduction, herbivory and allelopathy reduction, no herbivory and
 allelopathy reduction). Within each block, the monocultures and species
 mixtures of the native plants as well as the three monocultures of
 <em>S. canadensis</em> were randomly assigned to different
 plots. In the blocks with herbivory reduction, we sprayed each plot every
 10 days with a 0.5-L solution of Abamectin (1:500, v:v; Hongke Biochemical
 Co., Ltd, Shijiazhuang city, Hebei province in China), which is a
 broad-spectrum insecticide. The other plots were sprayed with tap water
 instead, and the plants in those plots were considered to be fully exposed
 to herbivory. In the blocks with allelopathy reduction, we had, prior to
 planting, mixed 8 L (~4 kg) of activated carbon (HG/T 3491-1999, Wuxi
 Yatai United Chemical Co., Ltd) into the top 20 cm of soil in each plot,
 resulting in a concentration of 1: 100 (v:v). Activated carbon can adsorb
 chemicals, including allelopathic substances produced by plants, and
 thereby supposedly neutralizes or reduces allelopathic interactions
 (Inderjit and Callaway, 2003; Prati & Bossdorf, 2004; Ridenour
 & Callaway, 2001; Yuan et al., 2022). However, as activated carbon
 may have undesired side-effects on plant growth (Lau et al., 2008; Kabouw
 et al., 2010; Weißhuhn & Prati, 2009; Wurst et al., 2010), we
 tested for side effects in the additional S. canadensis monoculture plots.
 We found no significant effect of activated carbon on biomass production
 of S. canadensis (Appendix S1: Fig. S1), which suggests that side effects
 of activated carbon were absent or minimal in our study. Therefore, plots
 without activated carbon were considered to have the full strength of
 allelopathic interactions among the plants.
 <em>Measurements</em>
 To quantify the likelihood of
 resource competition in each plot, we measured light interception by the
 canopy and soil nutrient contents. Light interception in each plot was
 measured twice on cloud-free days (September 25–28, 2018 and September
 21–24, 2019). Photosynthetically active radiation (PAR) was measured at
 three randomly selected points within the central 1.6 m × 1.6 m area of
 each plot using a PAR ceptometer (GLZ-C, Zhejiang Top Instrument Co., Ltd,
 China). This was done at mid-day (between 11:00 and 14:00), when the solar
 angle was maximal. Light-interception proportion was calculated as (PAR
 above the canopy - PAR at ground level)/PAR above the canopy.
 To
 determine nutrient contents, we took three soil cores (6.4 cm in diameter
 and 20 cm in depth) from the central 1.6 m × 1.6 m area of each plot, and
 thoroughly mixed the three cores to obtain one composite sample per plot.
 This was done on the days that we took the light measurments. Samples were
 air-dried for two weeks and then sieved through a 2-mm mesh. After that,
 20 g of each soil sample was ground into powder and used to measure soil
 total nitrogen and total phosphorus contents with an autoanalyzer
 (Autoanalyzer 3, BRAN+LUEBBE, Germany). We also determined the total
 organic carbon content using the method of Nelson and Sommers
 (1982). Before biomass harvest in 2018 and 2019 years,
 in each plot of herbivory and allelopathy and that of herbivory and
 allelopathy reduction treatment, 30 leaves (10 leaves was randomly
 selected in each of upper, middle and lower layers of canopy) of each
 species living in each plot were surveyed for damage by herbivores. If the
 total number of leaves of a species was less than 30, and the leaves were
 all surveyed. We counted the number of leaves with herbivore damage (e.g.
 holes). We then quantified the herbivory ratio as the number of leaves
 with herbivory damage divided by 30. In addition, we calculated for the
 damaged leaves on each plant species a herbivory intensity (an
 herbivory-severity index) as the average proportion of leaf area lost due
 to herbivory. To measure foliar flavonoid content of native speies and
 <em>S. canadensis</em>, five plots that having the species (if
 the species have died and another plot was reselected) were randomly
 selected in each of four treatments. Leaf samples (5-10 number) of the
 species were randomly gathered in each of five plots. The method of Shen
 et al. (2005) was used to measure the foliar flavonoid
 content. To determine the biomass production of plants,
 and whether this changed over time, we did two harvests. The first one was
 done in October 2018, one year after the invaders had been planted, and
 the second one was done in October 2019. At the first harvest, we
 collected aboveground biomass of all living plants in one half (1 m × 2 m)
 of each plot, and at the second harvest, we did this for the remaining
 halves. At both harvests, the plants were sorted to species, dried at 80°C
 for 48 h and weighed.
description: This Wang_2023_DATA_README.txt file was generated on 2023-11-21 by Jiang
 Wang GENERAL INFORMATION 1\. Title of Dataset: Data from: Herbivory and
 allelopathy contribute jointly to the diversity-invasibility relationship.
 \** ** 2\. Author Information Corresponding Investigator Name: Dr Xiaoyan
 Wang Institution: School of Life Science/Zhejiang Provincial Key
 Laboratory of Plant Evolutionary Ecology and Conservation, Taizhou
 University, China Email: [wxy3470117@163.com](mailto:wxy3470117@163.com)
 Co-investigator 1 Name: Dr Jiang Wang Institution: School of Life
 Science/Zhejiang Provincial Key Laboratory of Plant Evolutionary Ecology
 and Conservation, Taizhou University, China Email:
 [wangjiang@tzc.edu.cn](mailto:wangjiang@tzc.edu.cn) Co-investigator 2
 Name: Dr Song Gao Institution: School of Life Science/Zhejiang Provincial
 Key Laboratory of Plant Evolutionary Ecology and Conservation, Taizhou
 University, China Email: [lygs2008@163.com](mailto:lygs2008@163.com)
 Co-investigator 3 Name: Hefang Hong Institution: Linhai Branch of Taizhou
 Ecological Environment Bureau,Linhai 317000, China Email:
 [6598958@qq.com](mailto:6598958@qq.com) Co-investigator 4 Name: Wei Xue
 Institution: School of Life Science/Zhejiang Provincial Key Laboratory of
 Plant Evolutionary Ecology and Conservation, Taizhou University, China
 Email: [x_wei1988@163.com](mailto:x_wei1988@163.com) Co-investigator 5
 Name: Jiwei Yuan Taizhou Ecological and Environmental Monitoring Center,
 Taizhou, China Email: [46785047@qq.com](mailto:46785047@qq.com)  
 Co-investigator 6 Name: Mark van Kleunen School of Advanced Study, Taizhou
 University, Taizhou 318000, China;4 Ecology, Department of Biology,
 University of Konstanz, 78464 Konstanz, Germany Email:mark.
 [vankleunen@uni-konstanz.de](mailto:vankleunen@uni-konstanz.de)  
 Co-investigator 7 Name: Junmin Li School of Advanced Study, Taizhou
 University, Taizhou 318000, China;4 Ecology, Department of Biology,
 University of Konstanz, 78464 Konstanz, Germany Email:
 [limtzc@126.com](mailto:limtzc@126.com) 3\. Date of data collection: 2018
 and 2019 4\. Geographic location of data collection: Taizhou, Zhejiang,
 China 5\. Funding sources that supported the collection of the data:
 National Key R&D Program of China (2016YFC1201100), the Natural
 Science Foundation of Zhejiang Province (No. LY22C030001, LTY22C030004),
 the National Natural Science Foundation of China (No. 31870504), the
 Taizhou University National Funds for Distinguished Young Scientists (No.
 2017JQ005, 2019JQ005). 6\. Recommended citation for this dataset: Wang et
 al. (2023), Data from: Herbivory and allelopathy contribute jointly to the
 diversity-invasibility relationship, Dryad, Dataset \** ** DATA & FILE
 OVERVIEW 1\. Description of dataset To construct experimental communities
 of different species richness (1, 2, 4 and 8 species), we used eight
 native species (*Arctium lappa* L., *Patrinia villosa* (Thunb.) Juss.,
 *Achyranthes bidentata* Bl., *Reynoutria japonica* Houtt., *Taraxacum
 mongolicum* Hand.-Mazz., *Plantago asiatica* Ledeb., *Platycodon
 grandiflorus* (Jacq.) A. DC., *Antenoron filiforme* (Thunb.) Rob. et
 Vaut.). These species were chosen because they are common in natural
 shrub-grassland communities in the Zhejiang province, China, and two of
 them, P. villosa and A. bidentata, are usually very dominant in those
 communities. *Solidago canadensis* L. was selected as the invader. In May
 2017, we prepared 188 plots of 2 m × 2 m in a previous agricultural field
 that had been abandoned seven years earlier, in Luqiao (121°23′46.16″ E,
 28°31′21.62″ N). During October to December 2016, the top 10 cm of the
 soil was removed to deplete the seed bank, and the remaining soil was
 ploughed and repeatedly rotovated to produce a fine soil. In March to
 April 2017, the plots were weeded three times to further deplete the seed
 bank. The experimental area was divided into four blocks of 47 plots. In
 each plot, we constructed monocultures for each of the eight species with
 one replicate each, six two-species mixtures of different compositions and
 with three replicates each, six four-species mixtures with different
 compositions and three replicates each, and one eight-species mixture with
 three replicates. Each native species was used in an equal number of plots
 at each species-richness level. When the invasion treatment started, we
 also established three monoculture plots of the invader (*S. canadensis*)
 to assess potential side effects of activated carbon application 2\. File:
 Name: Wang_2023_data.xlsx Description: All data of this study.
 METHODOLOGICAL INFORMATION Treatments In the herbivory-reduction
 treatment, we sprayed the plots with an insecticide. To avoid that the
 insecticide spray would drift into plots of the non-insecticide treatment,
 we randomly assigned the four blocks, each with 47 plots, to one of the
 four experimental treatments (i.e. herbivory reduction, allelopathy
 reduction, herbivory and allelopathy reduction, no herbivory and
 allelopathy reduction). Within each block, the monocultures and species
 mixtures of the native plants as well as the three monocultures of S.
 canadensis were randomly assigned to different plots. In the blocks with
 herbivory reduction, we sprayed each plot every 10 days with a 0.5-L
 solution of Abamectin (1:500, v:v; Hongke Biochemical Co., Ltd,
 Shijiazhuang city, Hebei province in China), which is a broad-spectrum
 insecticide. The other plots were sprayed with tap water instead, and the
 plants in those plots were considered to be fully exposed to herbivory. In
 the blocks with allelopathy reduction, we had, prior to planting, mixed 8
 L (~4 kg) of activated carbon (HG/T 3491-1999, Wuxi Yatai United Chemical
 Co., Ltd) into the top 20 cm of soil in each plot, resulting in a
 concentration of 1: 100 (v:v). Activated carbon can adsorb chemicals,
 including allelopathic substances produced by plants, and thereby
 supposedly neutralizes or reduces allelopathic interactions (Inderjit and
 Callaway, 2003; Prati & Bossdorf, 2004; Ridenour & Callaway, 2001;
 Yuan et al., 2022). However, as activated carbon may have undesired
 side-effects on plant growth (Lau et al., 2008; Kabouw et al., 2010;
 Weißhuhn & Prati, 2009; Wurst et al., 2010), we tested for side
 effects in the additional S. canadensis monoculture plots. We found no
 significant effect of activated carbon on biomass production of S.
 canadensis (Appendix S1: Fig. S1), which suggests that side effects of
 activated carbon were absent or minimal in our study. Therefore, plots
 without activated carbon were considered to have the full strength of
 allelopathic interactions among the plants. Measurements To quantify the
 likelihood of resource competition in each plot, we measured light
 interception by the canopy and soil nutrient contents. Light interception
 in each plot was measured twice on cloud-free days (September 25–28, 2018
 and September 21–24, 2019). Photosynthetically active radiation (PAR) was
 measured at three randomly selected points within the central 1.6 m × 1.6
 m area of each plot using a PAR ceptometer (GLZ-C, Zhejiang Top Instrument
 Co., Ltd, China). This was done at mid-day (between 11:00 and 14:00), when
 the solar angle was maximal. Light-interception proportion was calculated
 as (PAR above the canopy - PAR at ground level)/PAR above the canopy. To
 determine nutrient contents, we took three soil cores (6.4 cm in diameter
 and 20 cm in depth) from the central 1.6 m × 1.6 m area of each plot, and
 thoroughly mixed the three cores to obtain one composite sample per plot.
 This was done on the days that we took the light measurments. Samples were
 air-dried for two weeks and then sieved through a 2-mm mesh. After that,
 20 g of each soil sample was ground into powder and used to measure soil
 total nitrogen and total phosphorus contents with an autoanalyzer
 (Autoanalyzer 3, BRAN+LUEBBE, Germany). We also determined the total
 organic carbon content using the method of Nelson and Sommers (1982).
 Before biomass harvest in 2018 and 2019 years, in each plot of herbivory
 and allelopathy and that of herbivory and allelopathy reduction treatment,
 30 leaves (10 leaves was randomly selected in each of upper, middle and
 lower layers of canopy) of each species living in each plot were surveyed
 for damage by herbivores. If the total number of leaves of a species was
 less than 30, and the leaves were all surveyed. We counted the number of
 leaves with herbivore damage (e.g. holes). We then quantified the
 herbivory ratio as the number of leaves with herbivory damage divided by
 30. In addition, we calculated for the damaged leaves on each plant
 species a herbivory intensity (an herbivory-severity index) as the average
 proportion of leaf area lost due to herbivory. To measure foliar flavonoid
 content of native speies and S. canadensis, five plots that having the
 species (if the species have died and another plot was reselected) were
 randomly selected in each of four treatments. Leaf samples (5-10 number)
 of the species were randomly gathered in each of five plots. The method of
 Shen et al. (2005) was used to measure the foliar flavonoid content. To
 determine the biomass production of plants, and whether this changed over
 time, we did two harvests. The first one was done in October 2018, one
 year after the invaders had been planted, and the second one was done in
 October 2019. At the first harvest, we collected aboveground biomass of
 all living plants in one half (1 m × 2 m) of each plot, and at the second
 harvest, we did this for the remaining halves. At both harvests, the
 plants were sorted to species, dried at 80°C for 48 h and weighed.
 DATA-SPECIFIC INFORMATION FOR: Wang_2023_data.xlsx 1\. Number of
 variables: 12 2\. Number of cases/rows: Sheet 1 'data': 377
 Sheet 2 'biomass': 1121 3\. Variable List: plot_name: Plot name.
 year: 1 2018 year, 2 2019 year. Species_richness: Species number in plot
 community_composition: Communities with different composition, 1-21
 indicate 21 types communities with different compositions. allelopathy: 0
 allelopathy reduction, 1 allelopathy. herbivory: 0 herbivory reduction, 1
 herbivory. total_organic_carbon: total organic carbon concentration in
 soil of each plot. total_nitrogen: total nitrogen concentration in soil of
 each plot.       total_phosphorus : total phosphorus concentration in soil
 of each plot. light_interception: light interception of plant communities
 of each plot. invader_biomass: the biomass of invader in each plot.
 complementarity_effect: Complementarity effect among native species. 
 selection_effect: Selection effect among native species. species name:
 latin name of native species in each plot. biomass: biomass of each native
 species originally planted in each plot. 4\. Missing data codes: N/A: not
 applicable 5 Script of R library(nlme) library(lme4) library(lmerTest)
 setwd("E:/wei.xue/Jiang Wang") HDR &lt;-
 read.csv("data_name.csv", header = TRUE, stringsAsFactors =
 TRUE) head(HDR) str(HDR) HDR$Y &lt;- as.factor(HDR$Y) HDR$H &lt;-
 as.factor(HDR$H) HDR$richness &lt;- as.numeric(HDR$richness) HDR$A &lt;-
 as.factor(HDR$A) HDR$species.composition &lt;-
 as.factor(HDR$species.composition) \####analysis of richenss effects mod1
 &lt;- lme(biomass~ scale(richness,scale = FALSE)  + A + H +              
 scale(richness,scale = FALSE):A + scale(richness,scale = FALSE):H + A:H +
               scale(richness,scale = FALSE):A:H             ,random = ~ 1
 | species.composition,             na.action = na.exclude,            
 data = subset(HDR, Y=="1"))# change "1" to
 "2" to test effects in year 2 summary(mod1) anova(mod1)
 intervals(mod1) var &lt;-mod1$apVar var par&lt;-attr(var,
 "Pars") par \#get variance explained by the random factor
 vc&lt;-exp(par)^2 vc \###analysis including diversity effects mod3 &lt;-
 lme(biomass~ H + A + complementarity.effect.transform +
 selection.effect.transform +             
 H:complementarity.effect.transform + A:complementarity.effect.transform +
 H:selection.effect.transform + A:selection.effect.transform +             
 H:A:complementarity.effect.transform + H:A:selection.effect.transform
             ,random = ~ 1 | species.composition,             na.action =
 na.exclude,             data = subset(HDR, Y=="1"))# change
 "1" to "2" to test effects in year 2 summary(mod3)
 anova(mod3) intervals(mod3) var &lt;-mod3$apVar var par&lt;-attr(var,
 "Pars") par vc&lt;-exp(par)^2 Vc 6\. Phyton code  # Import
 dataset import pandas as pd \# file can be replaced with SEMdata2019.csv
 file = r'F:\wj\20231003\SEMdata2018.csv'  data =
 pd.read_csv(file) data.head() \# correlations corr_matrix = data.corr()
 corr_matrix["invader_biomass"].sort_values(ascending=False) \#
 Spliting target variable and independent variables X =
 data.drop(['处理','lable','year','invader_biomass','invader_invader_native'], axis = 1) y = data['invader_biomass'] X.head() \# Splitting to training and testing data from sklearn.model_selection import train_test_split X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.3, random_state = 4) \# Import Random Forest Regressor from sklearn.ensemble import RandomForestRegressor \# Create a Random Forest Regressor reg = RandomForestRegressor() \# Train the model using the training sets reg.fit(X_train, y_train) \# Model prediction on train data y_pred = reg.predict(X_train) from sklearn import metrics as skmet \# Model Evaluation print('R^2:',skmet.r2_score(y_train, y_pred)) print('Adjusted R^2:',1 - (1-skmet.r2_score(y_train, y_pred))*(len(y_train)-1)/(len(y_train)-X_train.shape[1]-1)) print('MAE:',skmet.mean_absolute_error(y_train, y_pred)) print('MSE:',skmet.mean_squared_error(y_train, y_pred)) \# Predicting Test data with the model y_test_pred = reg.predict(X_test) \# Model Evaluation - Test set acc_rf = skmet.r2_score(y_test, y_test_pred) print('R^2:', acc_rf) print('Adjusted R^2:',1 - (1-skmet.r2_score(y_test, y_test_pred))*(len(y_test)-1)/(len(y_test)-X_test.shape[1]-1)) print('MAE:',skmet.mean_absolute_error(y_test, y_test_pred)) print('MSE:',skmet.mean_squared_error(y_test, y_test_pred)) num_vars = len(X.columns) num_vars names = list(X.columns) bounds = list(map(list,zip(X.min()-X.min()*0.95,X.max()+X.max()*0.95))) names \# print("min:",X.min()) \# print("max:",X.max()) \# pip install SALib from SALib.sample import saltelli from SALib.analyze import sobol from SALib.test_functions import Ishigami import numpy as np import SALib \# Define the model inputs problem = {     'num_vars': num_vars,     'names': names,     'bounds': bounds } \# Generate samples param_values = SALib.sample.saltelli.sample(problem, 2**14) \# Run model (example) Y = reg.predict(param_values) \# Perform analysis Si = sobol.analyze(problem, Y) total_Si, first_Si, second_Si = Si.to_df() total_Si, first_Si import matplotlib.pyplot as plt import matplotlib.colors as mcolors import random x=list(total_Si.index) y=list(total_Si['ST'])     colors=list(map(list,list(mcolors.cnames.items())[10:10+len(x)])) colors=[i[1] for i in colors] colors \# # for i in range(len(x)): \# #     colors.append()[random.randint(0,10)][1])   def autolabel(rects):     for rect in rects:         height = rect.get_height()         print(height)         plt.text(rect.get_x()+rect.get_width()/2.-0.2, 1.03*height, '%s' % round(height,2), size=10, family="Times new roman")   cm = plt.bar(x, y, color=colors,)   plt.xticks(rotation=45) autolabel(cm)   \# # plt.grid(True,linestyle=':',color='r',alpha=0.6) plt.show()   \################################# import shap plt.figure() shap.initjs() model = reg \# explain the model explainer = shap.TreeExplainer(model) shap_values = explainer.shap_values(X[0:120]) \# visualize the impact of each features 7\. shap.summary_plot(shap_values, X[0:120])