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.