Chapter 3- SPECIES-LEVEL DIFFERENCES IN PLANT TRAITS AND THEIR EFFECT ON ARTHROPOD COMMUNITY COMPOSITION

1- Load packages and tables
2- Plant trait ~ Species
- 2.1- Linear mixed models on all plant traits
- 2.2- Plant trait distributions
3- Arthropod guild ~ Species
- 3.1- Linear mixed models
- 3.2- Visualization of variance in arthropod guilds across species
4- Arthropod orders ~ Species
- 4.1- Permanova
- 4.2- Ordination plots
5- Plant traits ~ Species
- 5.1- Permanova
- 5.2- Ordination plots
6- Arthropod community composition ~ Species constrained by plant traits
- 6.1- Distance-based Canonical Correspondence Analysis (dbCCA)
- 6.2- Ordination plots

Description of research

This study seeks to describe the interspecific variation in plant traits and arthropod community commposition and the degree to which differences in plant traits explain the differences in arthropod community composition across 3 years in a common garden of Coastal Sage Scrub plants.

The metadata used for our analyses was collected annually in April-May for three consecutive years from 2021-2023 in a common garden at Concordia University Irvine. The common garden consisted of 13 species of woody shrubs native to the Coastal Sage Scrub ecosystem, arranged in 20 plots, each with 1 individual plant per species.

Data collected, abbreviations used and variable type assigned for the metadata can be found in the table below.

Data	Description	Column Name	Variable Type
Plant #	Individual plant	Plant	Factor
Plant Species	Plant Species	Species	Factor
Plot	Individual Plot	Plot	Factor
Year	Year data collected	Year	Factor
Specific leaf area	Leaf area : dry mass	SLA	Numeric
Percent water content	Leaf water content	PWC	Numeric
Relative growth rate	ln(final mass) - ln(initial mass)/years	RGR	Numeric
Plant quality	Average mass of surviving spodoptera	Avg_spod_mass	Numeric
Carbon to Nitrogen ratio	C:N ratio	CN	Numeric
Percent Carbon	Percent Carbon	Perc_C	Numeric
Percent Nitrogen	Percent Nitrogen	Perc_N	Numeric
Total arthropod density	Total number of arthropods	Abundance_Total	Numeric
Herbivore density	Total number of herbivores	Abundance_Herb	Numeric
Predator density	Total number of predators	Abundance_Pred	Numeric
Spider density	Total number of spiders	Aranea	Numeric
Mite density	Total number of mites	Acarina	Numeric
Spider density	Total number of spiders	Aranea	Numeric
True bug density	Total number of true bugs.	Heteroptera	Numeric
Wasp density	Total number of wasps	Hymenoptera_Vespidae	Numeric
Mantis density	Total number of mantis	Mantodea	Numeric
Leafhopper density	Total number of leafhoppers	Auchenorrhyncha	Numeric
Aphid density	Total number of aphids	Sternorrhyncha	Numeric
Thrip density	Total number of thrips	Thysanoptera	Numeric
Butterfly/moth density	Total number of adult butterflies/moths	Lepidoptera	Numeric
Caterpillar density	Total number of caterpillars	Lepidoptera_juvenile	Numeric
Honeybee density	Total number of honeybees	Hymenoptera_Apoidea	Numeric
Fly density	Total number of flies	Diptera	Numeric
Beetle density	Total number of beetles	Coleoptera	Numeric
Ant density	Total number of ants	Hymenoptera_Formicidae	Numeric
Springtail density	Total number of springtails	Collembola	Numeric
Barklice density	Total number of barklice	Psocoptera	Numeric

Methods

Plants trait and arthropod data was collected from 2021-2023 during peak growth season in May. Spodoptera mass was collected after a 10-day laboratory feeding assay. Briefly, Spodoptera exigua, a generalist herbivore, were hatched from eggs and immediately fed leaves from each individual plant. At the end of 10 days, the number of surviving larvae and final mass of each larvae were recorded. See You, An, and Li (2020) for details. Plant estimated biomass was collected using a branch-biomass estimation method as previously described in Nell and Mooney (2019). Arthropods were collected from plants using a vacuum method as previously described in Pratt et al. (2017), sorted from plant material and identified to taxonomic order.

1- Load packages and tables

2- Plant trait ~ Species

Use linear mixed models to evaluate plant traits that differ significantly among species.

2.1- Linear mixed models on all plant traits

Plant trait ~ Species + (1|Plant) + (1|Plot)

SLA (log-transformed)

PWC- non-normal distribution of residuals

Spodoptera mass- problem with convergence of the model + residuals are not normally distributed

Chisq	Df	Pr(>Chisq)
220	13	8.88e-40

RGR

Percent C- residuals non-normally distributed

Percent N

C:N Ratio- residuals non-normally distributed

Table 1 Summary statistics for linear mixed effects models testing for among-species variation in plant traits

Plant trait	Chi Sq	Chi df	P value	Marginal R2	Conditional R2
SLA	355.35	13	<0.0001*	0.69	0.85
PWC	332.02	13	<0.0001*	0.66	0.83
Plant quality	220.34	13	<0.0001*	0.60	0.80
RGR	130.24	12	<0.0001*	0.43	0.72
%C	404.64	13	<0.0001*	0.69	0.70
%N	253.84	13	<0.0001*	0.58	0.80
C:N ratio	339.80	13	<0.0001*	0.66	0.79

2.2- Plant trait distributions

Among-species variation in plant traits, mean +/- SE

Species	se_SLA	se_PWC	se_Spod	se_RGR	se_C	se_N	se_CN	trait	mean_value	se_value
Ag	0.364	0.0118	0.00485	0.067	0.142	0.0632	0.443	SLA cm2 g-1	162	0.364
Ag	0.364	0.0118	0.00485	0.067	0.142	0.0632	0.443	PWC	0.711	0.0118
Ag	0.364	0.0118	0.00485	0.067	0.142	0.0632	0.443	S. exigua growth rate	0.0355	0.00485
Ag	0.364	0.0118	0.00485	0.067	0.142	0.0632	0.443	RGR	0.183	0.067
Ag	0.364	0.0118	0.00485	0.067	0.142	0.0632	0.443	% Carbon	44.5	0.142
Ag	0.364	0.0118	0.00485	0.067	0.142	0.0632	0.443	% Nitrogen	2.9	0.0632

3- Arthropod guild ~ Species

Use linear mixed models to assess among-species variation in arthropod guilds.

3.1- Linear mixed models

Equation #1- Repeated measures using year, plot and plant as random effects.

Arthropod guild ~ Species + (1|Year) + (1|Plot) + (1|Plant)

Equation #2- Averaging guilds across years for each plant with plot as random effect

Average arthropod guild ~ Species + (1|Plot)

Summary statistics for Linear mixed effect models

Arthropod.type	Random.effects	Chi.Sq	Chi.df	P.value	Marginal.R2	Conditional.R2
Predator	Year/Plant/Plot	108.81	13	<0.0001*	0.13	0.29
Avg Predator	Plot	373.47	13	<0.0001*	0.32	0.36
Herbivore	Year/Plant/Plot	94.25	13	<0.0001*	0.13	0.18
Avg Herbivore	Plot	310.13	13	<0.0001*	0.29	0.31
Total	Year/Plant/Plot	126.94	13	<0.0001*	0.16	0.21
Avg Total	Plot	409.35	13	<0.0001*	0.34	0.38

Takeaway

I used average arthropod counts across all years of the study instead of each year separately. Plant and plot were included as random effects.

3.2- Visualization of variance in arthropod guilds across species

LS mean tables for Herbivore, Predator and Total abundance (all years combined)

Herbivore

Species	lsmean	SE	df	lower.CL	upper.CL
Acmispon_glaber	125.336961	27.63643	598	71.060696	179.61323
Artemisia_californica	71.987285	22.48779	598	27.822637	116.15193
Diplacus_aurantiacus	12.320459	22.73922	598	-32.337981	56.97890
Encelia_californica	28.574674	23.01786	598	-16.631005	73.78035
Eriogonum_fasciculatum	17.892474	23.77993	598	-28.809856	64.59480
Grindelia_camporum	80.134053	27.67969	598	25.772839	134.49527
Isocoma_menziesii	46.468465	24.17745	598	-1.014573	93.95150
Malacothamnus_fasciculatus	161.047598	24.62397	598	112.687618	209.40758
Malacothrix_saxatilis	4.731602	24.79710	598	-43.968389	53.43159
Mirabilis_laevis	27.709308	22.92159	598	-17.307293	72.72591
Salvia_apiana	66.314342	23.54401	598	20.075343	112.55334
Salvia_mellifera	75.307862	23.32124	598	29.506380	121.10934
Sisyrinchium_bellum	3.826240	24.80862	598	-44.896372	52.54885
Stachys_ajugoides_var_rigida	-1.514460	25.61004	598	-51.811020	48.78210

Predator

Species	lsmean	SE	df	lower.CL	upper.CL
Acmispon_glaber	9.560587	2.535040	598	4.5819231	14.539251
Artemisia_californica	11.350171	2.289203	598	6.8543162	15.846026
Diplacus_aurantiacus	8.841596	2.300804	598	4.3229580	13.360235
Encelia_californica	6.943668	2.313252	598	2.4005828	11.486752
Eriogonum_fasciculatum	6.611974	2.349628	598	1.9974471	11.226500
Grindelia_camporum	8.717712	2.537419	598	3.7343748	13.701048
Isocoma_menziesii	12.737303	2.368477	598	8.0857592	17.388848
Malacothamnus_fasciculatus	8.525310	2.389095	598	3.8332734	13.217346
Malacothrix_saxatilis	4.045409	2.396645	598	-0.6614541	8.752273
Mirabilis_laevis	2.099486	2.308992	598	-2.4352335	6.634205
Salvia_apiana	7.468225	2.337749	598	2.8770282	12.059422
Salvia_mellifera	7.814532	2.327850	598	3.2427754	12.386288
Sisyrinchium_bellum	3.036841	2.397077	598	-1.6708725	7.744554
Stachys_ajugoides_var_rigida	2.130050	2.434344	598	-2.6508524	6.910952

Total

Species	lsmean	SE	df	lower.CL	upper.CL
Acmispon_glaber	146.559476	30.46565	599	86.727001	206.39195
Artemisia_californica	132.444995	24.10878	599	85.096988	179.79300
Diplacus_aurantiacus	49.663619	24.42409	599	1.696364	97.63087
Encelia_californica	66.709470	24.77291	599	18.057159	115.36178
Eriogonum_fasciculatum	60.467278	25.72370	599	9.947671	110.98688
Grindelia_camporum	102.547263	30.51580	599	42.616298	162.47823
Isocoma_menziesii	83.993921	26.21755	599	32.504419	135.48342
Malacothamnus_fasciculatus	236.559814	26.77110	599	183.983185	289.13644
Malacothrix_saxatilis	26.754160	26.98462	599	-26.241798	79.75012
Mirabilis_laevis	60.567514	24.65241	599	12.151858	108.98317
Salvia_apiana	160.678396	25.42973	599	110.736134	210.62066
Salvia_mellifera	107.911727	25.15192	599	58.515052	157.30840
Sisyrinchium_bellum	14.529304	26.99811	599	-38.493152	67.55176
Stachys_ajugoides_var_rigida	7.788107	27.98543	599	-47.173389	62.74960

4- Arthropod orders ~ Species

Here we seek to determine to what degree differences in species explain variance in arthropod orders.

4.1- Permanova

We first ran a permanova using the following model:

Arthropod order1, Arthropod order 2, ... Arthrpod order X ~ Species

Arthropod data consisted of the percent composition of each arthropod order for each plant pooled across all years of the study.

	df	Sum.Sq	R2	F	Pr…F.
Species	13	22.844	0.41919	13.158	<0.0001 *
Residual	237	31.652	0.58081	NA
Total	250	54.496	1.00000	NA

TAKEAWAY

As expected, we see a significant amount of variation in arthropod community composition is explained by the differences in species (~42%).

4.2- Ordination plots

We will now visualize the results of this permanova analysis by creating ordination plots which take the multivariate data and condense them to two dimensional space. We will use two paths in this calculation: eigenanalysis using Principle Coordinate Analysis (PCoA) and non-metric multidimensional scaling methods (NMDS).

The PCoA assumes a linear correlation which may not be appropriate for this dataset but does allow us to create a unique ordination result and order axes in terms of variance explained.

4.2.1- PCoA Analysis

PCoA coordinates by species with influential arthropod order vectors

Ellipses + centroids by species with influential arhtropod order vectors

Ellipse centroids by species with influential arthropod order vectors

Influential orders (PCoA)

	Dim1	Dim2	p_value	r_squared	Order
Heteroptera	0.5081125	-0.3003102	0.001	0.3483645	Heteroptera
Auchenorrhyncha	0.4979738	-0.6995661	0.001	0.7373707	Auchenorrhyncha
Sternorrhyncha	-0.8966145	-0.3731780	0.001	0.9431795	Sternorrhyncha
Thysanoptera	0.0269831	0.6696459	0.001	0.4491537	Thysanoptera
Hymenoptera_Formicidae	0.0452956	0.3808026	0.001	0.1470623	Hymenoptera_Formicidae

TAKEAWAY:The first two axes for ordination plots using PCoA explain 23.37% and 20.1% of the variation in arthropod orders between species. Variation in species was mainly associated with variation in Thrips, Ants, Aphids, Leafhoppers and True Bugs.

4.2.2- NMDS Analysis (Stress = 0.16)

NMDS ellipsoid by species with influential arthropod order vectors

Influential orders (NMDS)

	NMDS1	NMDS2	NMDS3	p_value	r_squared	Order
Aranea	0.1595844	0.0008817	0.4577357	0.001	0.2349899	Aranea
Heteroptera	0.5382180	-0.3419455	-0.2530034	0.001	0.4706160	Heteroptera
Hynemoptera_Vespidae	-0.0494004	-0.1731885	0.3913999	0.001	0.1856285	Hynemoptera_Vespidae
Auchenorrhyncha	0.2423074	-0.6800503	-0.3916503	0.001	0.6745713	Auchenorrhyncha
Sternorrhyncha	-0.8645202	-0.1872644	-0.0204871	0.001	0.7828829	Sternorrhyncha
Thysanoptera	-0.0190376	0.5918458	-0.6005022	0.001	0.7112467	Thysanoptera
Coleoptera	0.0459593	0.0022978	0.4587997	0.001	0.2126147	Coleoptera
Diptera	0.0697957	0.1890036	0.2531161	0.001	0.1046616	Diptera
Hymenoptera_Formicidae	0.1636215	0.4932945	0.3311416	0.001	0.3797663	Hymenoptera_Formicidae

Takeaway- Many orders are influential in explaining the variance in community composition among species as noted in the table. The stress value for the plot is ~.16 which indicates the ordination is an ok representation of the variance in community composition explained by species.

5- Plant traits ~ Species

Here we seek to determine to what degree differences in species explain variance in each of the plant traits measured.

5.1- Permanova

We first ran a permanova using the following model:

Plant trait 1, Plant trait 2, ... Plant trait X ~ Species

Plant trait data consisted of the z-transformed values for each plant trait across all species. Because some of the resulting values were negative, we used Euclidean distance measures as opposed to Bray measures that were used for the arthropod permanova.

	df	Sum.Sq	R2	F	Pr..F.
Species	12	282.55	0.55641	14.768	<0.001***
Residual	140	228.45	0.44359
Total	152	515.00	1.00000

TAKEAWAY

As expected, we see a significant amount of variation in plant trait variance is explained by the differences in species (~57%).

5.2- Ordination plots

The PCoA assumes a linear correlation which may not be appropriate for this dataset but does allow us to create a unique ordination result and order axes in terms of variance explained.

PCoA coordinates by species with influential trait vectors

Ellipse + centroids by species with influential trait vectors

Ellipse centroids by species with influential trait vectors

Influential plant traits (PCoA)

Plant.trait	Dim1	Dim2	P.value	R2
SLA	-0.7647073	0.1079398	<0.001*	0.5964283
PWC	-0.7920841	0.0346527	<0.001*	0.6285981
RGR_1	0.2891126	-0.8901602	<0.001*	0.8759712
RGR_2	0.1302221	-0.9174157	<0.001*	0.8586094
Perc_C	0.6350711	0.1261060	<0.001*	0.4192181
Perc_N	-0.7716332	-0.3643893	<0.001*	0.7281974
CN	0.7962027	0.4090267	<0.001*	0.8012416

Takeaway- Very strong model fit with ~72% of variance in plant traits explained by differences in species. Every quadrant is occupied by as least 2 Species centroids (abbreviated with first letter of Genus and first letter of species) suggesting that plant species differ widely in variance between plant traits. Influential plant traits included all except the insect bioassay. SLA and PWC show strong postive covariance wherease RGR exhibits negative covariance to SLA and PWC. CN and Percent N also exhibit negative covariance.

6- Arthropod community composition ~ Species constrained by plant traits

6.1- Distance-based Canonical Correspondence Analysis (dbCCA)

We will used a constrained ordination method to model how well the variance in plant species explains the variance in arthropod community composition when constrained by the variance in plant traits. This will allow us to understand how plant traits influence arthropod community composition across plant species.

Anova results

Global

## Permutation test for capscale under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = dist_matrix_arthropod_dbcca ~ SLA + PWC + RGR_1 + Perc_N + Perc_C + CN, data = plant_traits_data)
##           Df SumOfSqs      F Pr(>F)    
## Model      6   3.3601 2.7124  0.001 ***
## Residual 146  30.1435                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

By individual plant traits

Trait	df	Sum.Sq	F	Pr..F.
SLA	1	0.2161	1.0465	0.375
PWC	1	0.6424	3.1114	0.003	**
RGR_1	1	0.6761	3.2748	0.002	**
Perc_N	1	0.6701	3.2458	0.002	**
Perc_C	1	0.5978	2.8953	0.004	**
CN	1	0.5576	2.7006	0.01	**
Residual	146	30.1435

Takeaway- It seems that RGR, %N and C:N ratio are having a moderately significant constraint on the variance in arthropod community composition due to differences in plant species

6.2- Ordination plots

dbCCA coordinates with influential arthropod order and plant trait vectors

Centroids + ellipses with influential arthropod order and plant trait vectors

Ellipse centroids with influential arthropod orders and plant trait vectors

Takeaway: only 9.12% of the variation in arthropod community composition is explained by the constrained ordination. Of this 9.12%, ~75% is explained by the two axes, CAP1 and CAP2. In contrast, 43.47% of the variation in arthropod community composition is explained by the first two axes, PCoA1 and PCoA2 in the unconstrained model.

Influential arthropod orders and plant traits

	CAP1	CAP2	Pr(>F)	R2
RGR	-0.2124979	-0.2840967	0.001	0.1258663
% N	-0.2823634	0.1942418	0.001	0.1174589
CN	0.3545982	-0.2193679	0.001	0.1738621
Aranea	0.3253173	-0.2240724	0.001	0.1560398
Heteroptera	-0.7514170	-0.2636978	0.001	0.6341641
Hynemoptera_Vespidae	0.3776614	-0.2480675	0.001	0.2041656
Auchenorrhyncha	-0.7448663	0.4600808	0.001	0.7665000
Sternorrhyncha	0.3907819	0.5292231	0.001	0.4327875
Thysanoptera	-0.0898314	-0.5994583	0.001	0.3674200
Coleoptera	0.3994898	-0.2902892	0.001	0.2438599
Diptera	0.3087944	-0.3436676	0.001	0.2134614
Hymenoptera_Formicidae	0.3859647	0.1982934	0.001	0.1882890

Nell, Colleen S., and Kailen A. Mooney. 2019. “Plant Structural Complexity Mediates Trade-Off in Direct and Indirect Plant Defense by Birds.” Ecology 100 (10). https://doi.org/10.1002/ecy.2853.

Pratt, Jessica D., Andrew Datu, Thi Tran, Daniel C. Sheng, and Kailen A. Mooney. 2017. “Genetically Based Latitudinal Clines in Artemisia Californica Drive Parallel Clines in Arthropod Communities.” Ecology 98 (1): 79–91. https://doi.org/10.1002/ecy.1620.

You, Yanrong, Chunpeng An, and Chuanyou Li. 2020. “Insect Feeding Assays with Spodoptera Exigua on Arabidopsis Thaliana.” BIO-PROTOCOL 10 (5). https://doi.org/10.21769/BioProtoc.3538.