Statistical analyses
All statistical analyses were carried out with R software (R Development
Core Team 2019).
1) Is variation in GSL diversity across species correlated with to
species’ phylogenetic distance? First, to assess the effect of species
identity on the entire GSL matrix, we used non-metric multidimensional
scaling (NMDS) implemented in the vegan package (Oksanen et al.2013). Differences in the GSL composition were tested using a
permutational multivariate ANOVA (PERMANOVA), using
the adonis function in the vegan package (Oksanen et
al. 2013). The Bray-Curtis metric was used to calculate dissimilarity
among samples for both the NMDS and PERMANOVA, although results were
robust to other distance metrics. Second, we performed a Mantel test
(9999 iterations) using the function mantel.test in thevegan package (Oksanen et al. 2013) between the
phylogenetic distance matrix and the chemical distance matrix across all
species to test for a potential correlation between phylogenetic
distance and chemical distance.
2) Does variation in GSL diversity simultaneously converge with
plant species adaptation to their specific environments? To address
this question, we performed a coinertia analysis, at the species level,
between the plant functional traits matrix and the GSL matrix using the
function coinertia in the package ade4 (Dray & Dufour
2007). Using these analyses, we were able to visually detect clustering
of species into four groups (see Results below). In other words, we
detected four distinct growth forms-GSLs clusters, that also separated
species according to their optimal habitat (see Fig. 2). We thus next
performed discriminant analyses to determine to what extent GSL profiles
could predict group assignment for each Cardamine species. We
performed a linear discriminant (LD) analysis based on cluster groups
using the function lda from the MASS package (Venables &
Ripley 2002). The quality of the resulting model was assessed through
the classification success derived from a jackknife-based
cross-validation (i.e. leave-one-out process) using the ‘CV’ argument of
the lda function. Overall, 74% of samples were correctly
classified, and the first LD of the model (LD1) accounted for about 90%
of between-group variances. Differences in the distribution of leaf GSLs
profiles along LD1 were tested with a pairwise Wilcoxon test coupled
with a p-value adjustment based on the Benjamini and Hochberg method
(Benjamini & Hochberg 1995).
3) How are different metrics of phytochemical diversity related to
plant-herbivore interaction? To address this question, we first
calculated seven different diversity indices for production of GSLs
across Cardamine species including; 1) the total GSL abundance
(Sum), 2) the number of individual compounds (S; i.e. chemical
richness), 3) the Shannon diversity index (H), 4) the chemical evenness
(J), 5) the functional divergence (FDiv), 6) the functional richness
(FRic), and 7) the Rao’s quadratic entropy (RaoQ), using the packageFD (Laliberté et al. 2014). For calculating functional
diversity indices, we included as functional traits of each GSL compound
their chemical class (aliphatic, aromatic, indole), the class of their
breakdown products (isothiocyanates, oxazolidine-2-thione,
oxazolidine-2-thione), and their molecular weight. We here propose this
new functional approach for organizing plant secondary metabolite
diversity since we assume that high functional diversity derived from
the different chemical classes of the GSLs correlate with increased
activity. To assess the effect of species on the different chemical
diversity indices, we ran one-way ANOVAs with each of the diversity
indices separately as response variable.
For measuring the effect of each of the GSL diversity indices on each of
the four different insects’ growth rate value we ran Bayesian
phylogenetic mixed effect models (BPMMs), as implemented in the R
package MCMCglmm (function MCMCglmm ) (Hadfield 2010).
MCMCglmm analyses allow taking into account species phylogenetic
relationship as random factor in the model. Because the response
variables followed a normal distribution, we used a MCMCglmm with a
Gaussian distribution (Hadfield 2010). Finally, we assessed the effect
of the four species’ groups derived from the coinertia analyses
described above on insect resistance by performing ANOVAs with insect
growth as response variable and species nested in the corresponding
group as explanatory variable for each herbivore insect independently.
Between groups, differences were assessed by pairwise comparisons using
Tukey HSD post-hoc tests.