Rapid and independent evolution of ancestral and novel defenses in a genus of toxic plants (Erysimum, Brassicaceae)

Phytochemical diversity is thought to result from coevolutionary cycles as specialization in herbivores imposes diversifying selection on plant chemical defenses. Plants in the speciose genus Erysimum (Brassicaceae) produce both ancestral glucosinolates and evolutionarily novel cardenolides as defenses. Here we test macroevolutionary hypotheses on co-expression, co-regulation, and diversification of these potentially redundant defenses across this genus. We sequenced and assembled the genome of E. cheiranthoides and foliar transcriptomes of 47 additional Erysimum species to construct a highly resolved phylogeny, revealing that cardenolide diversity increased rapidly rather than gradually over evolutionary time. Concentrations, inducibility, and diversity of the two defenses varied independently among species, with no evidence for trade-offs. Closely related species shared similar cardenolide traits, but not glucosinolate traits, likely as a result of specific selective pressures acting on distinct molecular diversification mechanisms. Ancestral and novel chemical defenses in Erysimum thus appear to provide complementary rather than redundant functions.

included in all three experiments, others could only be grown once due to limited seed availability or 159 germination. To maximize germination success, seeds were placed on water agar (1%) in Petri dishes 160 and cold-stratified for two weeks. After stratification, Petri dishes were moved to a growth chamber 161 set to 24 °C day / 22 °C night at a 16:8 h photoperiod. Viable seeds germinated within 3-10 days of 162 placement in the growth chamber. As soon as cotyledons had fully extended, we transplanted the 163 seedlings into 10 x 10 cm plastic pots filled with a mixture of peat-based germination soil 164 (Seedlingsubstrat, Klasmann-Deilmann GmbH, Geeste, Germany), field soil, sand, and vermiculite at 165 a ratio of 6:3:1:5. Plants were moved to a climate-controlled greenhouse set to 24 °C day / 16 °C night 166 and 60 % RH with natural light and supplemented artificial light set to a 14:10 h photoperiod. Plants Hombrechtikon, Switzerland), including a step for on-column DNase digestion, and following the 287 manufacturer's instructions. The purified total RNA was dissolved in 50 µL RNase-free water, split 288 into three aliquots, and stored at -80 °C until further processing. Assessment of RNA quality, library 289 preparation, and sequencing were all performed by the Next Generation Sequencing Platform of the 290 University of Bern (Bern, Switzerland). RNA quality was assessed in one aliquot per extract using a 291
Finally, coalescent species tree inference was performed with these 9,868 gene trees as input, using 322 ASTRAL-III v5.6.3 (Zhang et al. 2018). 323 In the second approach, we used transcriptome sequences translated by TransDecoder to 324 predict protein sequences, after which the longest predicted protein for each gene was retained. For E. 325 Characteristic fragmentation allowed us to identify candidate cardenolide compounds in a 432 genin-guided approach, where the presence of characteristic genin fragments in a chromatographic 433 peak indicated the likely presence of a cardenolide molecule. We then identified the parental mass of 434 cardenolides, additional fragments corresponding to the loss of the outer sugar moiety allowed us to 437 determine the mass and order of sugar moieties in the linear glycoside chain of the molecule. We 438 screened our data for the presence of glycosides of strophanthidin, digitoxigenin, and cannogenol, and 439 additional genins known to occur in Erysimum species (Makarevich et al. 1994 Aldrich) to each sample as an internal standard, but between-sample variation (technical noise) was 448 negligible compared to between-species variation. 449 For glucosinolate and cardenolide data separately, raw ion counts for each compound were 450 averaged across experiments to yield robust chemotype data. Raw ion counts were standardized by the 451 dry sample weight, possible dilution of samples, and internal standard concentrations (where 452 available). For pooled samples, ion counts were standardized by the average dry weight calculated 453 from all samples that contributed to a pool. The full set of standardized compound ion counts was 454 then analyzed using linear mixed effects models (package nlme v3.1-137 in R v3.5.3). Because 455 standardized ion counts still had a heavily skewed distribution, we applied a log(+0.1) transformation 456 to all values. Log-transformed ion counts were modelled treating experiment as a fixed effect, and a 457 species-by-compound identifier as the main random effect. Nested within the main random effect, we 458 fitted a species-by-compound-by-experiment identifier as a second random effect to account for the 459 difference of pooled or individual samples among experiments. The fixed effect of this model thus 460 captures the overall differences in compound ion counts between experiments, while the main random 461 effect captures the average deviation from an overall compound mean for each compound in each 462 species. We extracted the overall compound mean and the main random effects from these models, 463 providing us with average ion counts for each compound in each species on the log-scale. Negative 464 values on the log-scale were set to zero as they would correspond to values below the limit of reliable 465 detection of the LC-MS on the normal scale. 466 467

Inhibition of mammal Na + /K + -ATPase by leaf extracts
Although all cardenolides target the same enzyme in animal cells, structural variation among different 469 cardenolides can significantly influence binding affinity and thus affect toxicity (Dzimiri et al. 1987, 470 Petschenka et al. 2018). Cardenolide quantification from LC-MS mass signal intensity does not 471 capture such differences in biological activity, and furthermore may be challenging due to compound-472 specific response factors and narrow ranges of signal linearity. To evaluate whether total ion counts 473 are an appropriate and biologically relevant measure for between-species comparisons of defense 474 levels, we therefore quantified cardenolide concentrations by a separate method (Züst et  and further diluted 1:5, 1:50, and 1:500 using 10% DMSO. To quantify potential non-specific 486 enzymatic inhibition that could occur at high concentrations of plant extracts, we also included control 487 extracts from Sinapis arvensis leaves (a non-cardenolide producing species of the Brassicaceae) in 488 these assays. 489 Assays were carried out in 96-well microplate format. Reactions were started by adding 80 μL 490 of a reaction mix containing 0.0015 units of porcine Na + /K + -ATPase to 20 μL of leaf extracts in 10% 491 DMSO, to achieve final well concentrations (in 100 μL) of 100 mM NaCl, 20 mM KCl, 4 mM MgCl2, 492 50 mM imidazol, and 2.5 mM ATP at pH 7.4. To control for coloration of leaf extracts, we replicated 493 each reaction on the same 96-well plate using a buffered background mix with identical composition 494 as the reaction mix but lacking KCl, resulting in inactive Na + /K + -ATPases. Plates were incubated at 495 37 °C for 20 minutes, after which enzymatic reactions were stopped by addition of 100 μL sodium 496 dodecyl sulfate (SDS, 10% plus 0.05% Antifoam A) to each well. Inorganic phosphate released from 497 enzymatically hydrolyzed ATP was quantified photometrically at 700 nm following the method 498 described by Taussky and Shorr (1953). The absorbance values at four dilutions x are thus used to estimate the upper (A, fully active enzyme) 505 and lower (B, fully inhibited enzyme) asymptotes, the dilution value xmid at which 50% inhibition is 506 achieved, and a shape parameter scal. In order to estimate four parameters from four absorbance 507 values per extract, the scal parameter was fixed for all extracts and changed iteratively to optimize 508 overall model fit, judged by AIC. Individual plant extracts were treated as random effects to account 509 for lack of independence within extract dilution series. For each extract we estimated xmid from the 510 average model fit and the extract-specific random deviate. Using a calibration curve made with 511 ouabain ranging from 10 -3 to 10 -8 M that was included on each 96-well plate, we then estimated the 512 concentration of the undiluted sample in ouabain equivalents, i.e., the amount of ouabain required to 513 achieve equivalent inhibition.

Similarity in defense profiles between Erysimum species 534
To quantify chemical similarity among species, we performed separate cluster analyses on the 535 glucosinolate and cardenolide profile data averaged across the three experiments. For each species, 536 the log-transformed average ion counts of all compounds were converted to proportions (all 537 compounds produced by a species summing to 1). From this proportional data we then calculated 538 pairwise Bray-Curtis dissimilarities for all species pairs using function vegdist in the R package vegan 539 v2.5-4. We incorporated vegdist as a custom distance function for pvclust in the R package pvclust 540 v2.0 (Suzuki and Shimodaira 2014), which performs multiscale bootstrap resampling for cluster 541 analyses. We constructed dendrograms of glucosinolate and cardenolide profile similarities by fitting 542 hierarchical clustering models (Ward's D) and estimated support for individual species clusters from 543 10,000 permutations. To compare chemical similarity to phylogenetic relatedness we performed 544 principal coordinate analyses (PCoA) on Bray-Curtis dissimilarity matrices of glucosinolate and 545 cardenolide data using function pcoa in R package ape v5.0 (Paradis and Schliep 2019), and extracted 546 the first two principal coordinates for each defense trait to test for phylogenetic signal. 547 548

Relationship between plant traits and phylogenetic signal 549
We evaluated a prevalence of phylogenetic signal in chemical defense traits, myrosinase activity, and 550 principal coordinates for both chemical similarity matrices using Blomberg's K (Blomberg et al. 551 2003). K is close to zero for traits lacking phylogenetic signal; it approaches 1 if trait similarity among 552 related species matches a Brownian motion model of evolution, and it can be >1 if similarity is even 553 higher than expected under a Brownian motion model. We estimated K for all traits using function 554 phylosig in the R package phytools v0.6-60 (Revell 2012). Additionally, we used the geographic 555 coordinates for all species with known collection locations to construct pairwise geographic distances, 556 calculated pairwise geographic dissimilarities, performed principal coordinate analyses on the 557 geographic dissimilarity matrix, and estimated K for the first two components. 558 To test for directional effects in the evolution of compound number and abundance for both 559 glucosinolates and cardenolides, we applied Pagel's method (Pagel 1999). Specifically, we compared 560 a Brownian motion model of trait evolution to a model in which additionally a directional trend is 561 assessed by regressing the path length (i.e., molecular branch length from root to tip) against trait 562 values. For this analysis we used the concatenated 2,306-gene tree for which branch lengths are an 563 estimate of substitutions per site. Models were fit using function fitContinuous in R package geiger 564 v2.0.6.2 (Harmon et al. 2008), where the default setting fits a Brownian motion model, whereas the 565 additional argument 'model=drift' specifies a directional trend model. Support for directional trends 566 in defense traits was evaluated using likelihood-ratio tests between the two models. 567 568

E. cheiranthoides genome assembly 570
A total of 39.5 Gb of PacBio sequences with an average read length of 10,603 bp were assembled into 571 1,087 contigs with an N50 of 1.5 Mbp (Table 1). Hi-C scaffolding oriented 98.5 % of the assembly 572 into eight large scaffolds representing pseudomolecules (Table 1, Figure S1), while 216 small contigs 573 remained unanchored. The final assembly (v1.2) had a total length of 174.5 Mbp, representing 86% of 574 the estimated genome size of E. cheiranthoides and capturing 99% of the BUSCO gene set (Table 1, 575 Figure S2). Sequences were deposited under GenBank project ID PRJNA563696 and additionally are 576 provided at www.erysimum.org. A total of 29,947 gene models were predicted and captured 98% of 577 the BUSCO gene set ( Figure S3). In the presumed centromere regions of each chromosome, genic sequences were less abundant, whereas repeat sequences were more common (Fig 2A). Repetitive 579 sequences constituted approximately 29% of the genome (Table S2). Long terminal repeat 580 retrotransposons (LTR-RT) made up the largest proportion of the repeats identified ( Figure S4). 581 Among these, repeats in the Gypsy superfamily constituted the largest fraction of the genome (Table  582 S2). The majority of the LTR elements appeared to be relatively young, with most having estimated 583 insertion times of less than 1 MYA ( Figure S5). Synteny analysis showed evidence of several 584 chromosomal fusions and fissions between the eight chromosomes of E. cheiranthoides and the five 585 chromosomes of Arabidopsis ( Figure 2B). 586 587

Phylogenetic relationship of 48 Erysimum species 649
Assemblies of transcriptomes from 48 Erysimum species (including E. cheiranthoides) had N50 650 values ranging from 574 -2,160 bp (Table S3). Transcriptome assemblies contained completed genes 651 from 54% -94% of the BUSCO set and coding sequence lengths were generally shorter on average 652 than the E. cheiranthoides coding sequence lengths (Table S3). Transcriptome sequences were 653 deposited under GenBank project ID PRJNA563696 and at www.erysimum.org. The large number of 654 orthologous gene sequences identified among the E. cheiranthoides genome and the 48 transcriptomes 655 allowed us to infer phylogenetic relatedness with high confidence. We assume that the coalescent 656 phylogeny generated by the first approach of tree construction using 9,868 syntenic gene sequences 657 represents our best estimate of species relationships ( Figure 5). However, both phylogenies generated 658 by the second approach using 2,306 orthologous genes had highly similar tree topologies ( Figures S7,  659 S8), suggesting overall high reliability of our results regardless of the phylogenetic inference method.

Glucosinolate diversity and myrosinase activity 685
Across the 48 Erysimum species, we identified 25 candidate glucosinolate compounds with distinct 686 molecular masses and HPLC retention times (Table S4). Of these, 24 compounds could be assigned to 687 known glucosinolate structures with high certainty. The remaining compound appeared to be an 688 unknown isomer of glucocheirolin. Individual Erysimum species produced between 5 and 18 689 glucosinolates ( Figure 6A), and total glucosinolate concentrations were highly variable among species 690 ( Figure 6B). The ploidy level of species explained a significant fraction of total variation in the 691 number of glucosinolates produced (F4,38 = 4.63, p = 0.004), with hexaploid species producing the 692 highest number of compounds ( Figure S9). However, neither the numbers of distinct glucosinolates 693 nor the total concentrations exhibited a phylogenetic signal (Table 2). Similarly, there was no 694 evidence for a directional trend in either glucosinolate trait (Table 3)

712
Clustering species according to similarity in glucosinolate profiles mostly resulted in chemotype 713 groups corresponding to known underlying biosynthetic genes, although support for individual 714 species clusters was variable (Figure 7). The majority of species produced glucoiberin as the primary 715 glucosinolate. Of these, approximately half also produced sinigrin as a second dominant glucosinolate 716 compound. Further chemotypic subdivision, related to the production of glucocheirolin and 2-717 hydroxypropyl glucosinolate, appeared to be present but only had relatively weak statistical support. 718 However, eight species clearly differed from these general patterns. The species E. allionii (ALI), E.  (Table 2). 727 As glucosinolates require activation by myrosinase enzymes upon tissue damage by 728 herbivores, myrosinase activity in leaf tissue determines the rate at which toxins are released. We 729 quantified myrosinase activity of Erysimum leaf extracts and found it to be highly variable among 730 species ( Figure 6C). After grouping species into nine chemotypes defined by chemical similarity and 731 the production of characteristic glucosinolate compounds ( Figure 7C), we found that myrosinase 732 activity significantly differed among these chemotypes (Figure 8, F8,33 = 7.06, p < 0.001). Chemotypes 733 that predominantly accumulated methylsulfonyl glucosinolates, hydroxy glucosinolates, or indole 734 glucosinolates had low to negligible activity against the assayed glucosinolate sinigrin. It is important 735 to note that sinigrin is an alkenyl glucosinolate and activity with other, structurally dissimilar 736 glucosinolates may differ. After chemotype differences were accounted for, myrosinase activity was 737 related positively to total glucosinolate concentrations (F1,33 = 5.92, p = 0.021). Surprisingly, uncorrected myrosinase activity was the only glucosinolate-related trait that showed a significant 739 phylogenetic signal (Table 2).     Figure S9). To obtain an 766 estimate of biological activity and evaluate quantification from total MS ion counts, we used an 767 established assay that quantifies cardenolide concentrations from specific inhibition of animal Na + /K + -768 ATPase by crude Erysimum leaf extracts. We found generally strong enzymatic inhibition, with leaves 769 of Erysimum species containing an equivalent of 5.72 ± 0.12 µg mg -1 (± 1 SE) of the reference 770 cardenolide ouabain on average. Despite only producing trace amounts of cardenolides, E. collinum 771 (COL) extracts caused significantly stronger inhibition than the Brassicaceae control, S. arvensis 772 ( Figure S10). Overall, quantification of cardenolide concentrations by Na + /K + -ATPase inhibition was 773 highly correlated with the total MS ion count (Fig S8, r = 0.95, p < 0.001). Thus, the use of ion count 774 data for cross-species comparisons was appropriate for this purpose. Both the total numbers of 775 compounds and the total abundances exhibited a strong phylogenetic signal (Table 2), indicating that 776 closely-related species were more similar in their cardenolide traits than expected by chance. 777 However, there was again no evidence for a directional trend in the evolution of either number or 778 abundance of cardenolides (Table 3), suggesting a rapid rather than a gradual gain of cardenolide diversity, which is also evident from the considerable number of cardenolide compounds present in 780 the earliest-diverging species ( Figure 6D). distinguishable candidate cardenolide compounds identified across the 48 Erysimum species (Table  795   S5). Of these, 46 compounds had distinct molecular masses and mass fragments, while the remaining 796 compounds likely were isomers, sharing a molecular mass with another compound but having a 797 distinct HPLC retention time. The 95 putative cardenolides comprised nine distinct genins (Figures 9,  798 S11), the majority of which were glycosylated with digitoxose, deoxy hexoses, xylose, or glucose 799 moieties. Only digitoxigenin and cannogenol accumulated as free genins, while all other compounds 800 occurred as either mono-or di-glycosides. A major source of isomeric cardenolide compounds was 801 thus likely the incorporation of different deoxy hexoses of equivalent mass, such as rhamnose, fucose, 802 or gulomethylose. A subset of compounds had molecular masses that were heavier by 42.011 m/z than 803 known mono-or di-glycoside cardenolides. Such a gain in mass corresponds to the gain of an acetyl-804 group, and mass fragmentation patterns indicated that these compounds were acetylated on the first 805 sugar moiety (Table S5). Out of the nine detected genins, six had previously been described from 806 Erysimum species (Makarevich et al. 1994). In addition, we identified three previously undescribed 807 mass features with fragmentation patterns characteristic of cardenolide genins ( Figure S11). Of these 808 three, one matched an acetylated digitoxigenin (also known as oleandrigenin), a common cardenolide 809 in Nerium oleander. The other two matched molecular structures of digitoxigenin-formate (also 810 known as gitaloxigenin) and strophanthidin-formate. Formate adducts can sometimes be formed 811 during LC-MS due to the addition of formic acid to solvents, although this is less common with 812 positive ionization. To exclude to possibility that these were technical artefacts, we analyzed a subset 813 of extracts by LC-MS without the addition of formic acid and found both genin-formates at 814 comparable concentrations ( Figure S12). We therefore assume that all three novel structures are 815 natural variants of cardenolides produced by Erysimum plants, even though we currently lack final 816 structural elucidation. 817 Clustering of species based on similarity in cardenolide profiles revealed fewer obvious 818 species clusters than for glucosinolates, and particularly higher-level species clusters had only weak 819 statistical support (Figure 10). A clear exception to this was a species cluster that included E. 820 cheiranthoides (ECE) and E. sylvestre (SYL), which lacked several otherwise common cannogenol-821 and strophanthidin-glycosides, while accumulating unique digitoxigenin-glycosides. A second major 822 cluster that was visually apparent -yet not statistically significant -separated groups of species that 823 did or did not produce glycosides of the newly discovered putative strophanthidin-formate ( Figure  824 10). Similarity in cardenolide profiles among species quantified as the first and second principal 825 coordinate of the Bray-Curtis dissimilarity matrix exhibited a very strong phylogenetic signal (Table  826 2), suggesting that closely-related species not only were more similar in their total cardenolide 827 concentrations, but also had more similar cardenolide profiles than expected by chance. 828  Table S5 for additional compound information.

Macroevolutionary patterns in defense and inducibility 835
Given the very distinct patterns for glucosinolate and cardenolide diversity among Erysimum species, 836 it is unsurprising that concentrations of the two defense traits were not correlated (Pearson's 837 correlation: r = -0.09, p = 0.534). Foliar application of JA was expected to stimulate defense levels in 838 plant leaves, and among the 30 tested species, glucosinolate levels responded positively to JA, with 839 the majority of species increasing their foliar glucosinolate concentration (Figure 11). However, the 840 glucosinolate inducibility of a species was independent of constitutive glucosinolate levels (r = -0.26, 841 p = 0.169). By contrast, the majority of species exhibited lower cardenolide levels in response to JA, 842 resulting in lack of inducibility across species (Figure 11). The species E. crepidifolium (CRE) 843 heavily influenced inducibility patterns, as it not only had three times higher constitutive 844 concentrations of cardenolides than any other Erysimum species, but in addition pronouncedly 845 increased both glucosinolate and cardenolide concentrations in response to JA treatment ( Figure 11).  glucosinolate-specialized herbivores (Chew 1975, Wiklund and Åhrberg 1978, Renwick et al. 885 1989, Dimock et al. 1991. 886 As further evidence for the distinct roles of glucosinolates and cardenolides, the two defenses 887 responded differently to exogenous JA application. Glucosinolate concentrations were upregulated in 888 response to JA in the majority of species, with an average 52% increase relative to untreated controls. bridges. However, the node connecting the two geographic clades had by far the weakest support 926 across the whole phylogeny (local posterior probability = 0.5) and, as our phylogeny did not include 927 East Asian Erysimum species, it is possible that the addition of further species would change this 928 grouping. The remaining species belonged to one of three distinct central European clades with no 929 obvious geographic separation. As a prominent exception, the central European E. crepidifolium 930 (CRE) was most closely related to Greek Erysimum species rather than to other central European 931

species.
Despite vast morphological differences among sampled Erysimum species, the diversity in 933 glucosinolate profiles was relatively limited compared to the diversity that is present within 934 Arabidopsis (Kliebenstein et al. 2001). However, broader comparative studies of glucosinolate 935 diversity in other Brassicaceae species would be needed to provide a more natural 'baseline' for 936 glucosinolate diversity. The majority of Erysimum species produced glucoiberin as their main 937 glucosinolate. Aliphatic glucosinolates such as glucoiberin are derived from methionine in a process 938 that involves elongation and modification of a variable side-chain (Halkier and Gershenzon 2006), 939 and in this context the 3-carbon glucosinolate glucoiberin is one of the least biosynthetically complex 940 glucosinolates. However, the potential to produce additional aliphatic glucosinolates with longer side 941 chains clearly exists in the genus, as 4-, 5-, and 6-carbon glucosinolates with more complex 942 modifications were scattered across the phylogeny. A few species produced glucosinolates that are not 943 found in Arabidopsis, including a sub-class of aliphatic glucosinolates, the methylsulfonyl 944 glucosinolates. The homolog of 3-butenyl glucosinolate 2-hydroxylase (GS-OH), which in 945 Arabidopsis forms 2-hydroxy-but-3-enyl glucosinolate from 3-butenyl glucosinolate, does not have a 946 clear function in E. cheiranthoides due to the lack of alkenyl glucosinolates. However, it is possible 947 that the GS-OH homolog in E. cheiranthoides may code for the unknown enzyme that hydroxylates 4-948 methylsulfonylbutyl glucosinolate to form 3-hydroxy-4-methylsulfonylbutyl glucosinolate (Figure 4). 949 Methylsulfonyl glucosinolates are found in several Brassicaceae genera (Fahey et al. 2001 possibly 2-amino-5-methoxycarbonyl-pentanoic acid (Chisholm 1973). Modification of the amino 960 acid side chain during methionine-derived aliphatic glucosinolate biosynthesis as a pathway to 961 glucoerypestrin is less likely, due to the lower specific incorporation of 14 C-labeled methionine 962 compared to 14 C-labeled dicarboxylic acids into this compound (Chisholm 1973). In any case, the gain 963 of glucoerypestrin represents yet another evolutionary novelty in the Erysimum genus, but its relative 964 toxicity and the adaptive benefits of its production have yet to be elucidated. 965 We found no phylogenetic signal of glucosinolate chemotype, as more closely related species 966 were not more likely to share the same glucosinolate profile. The pattern of more complex 967 glucosinolates scattered across the phylogenetic tree may be generated by horizontal gene transfer, 968 during hybridization, or by repeated gains and losses of biosynthesis genes as species diverge. The latter may also be facilitated by changes in ploidy, as hexaploid species in particular accumulated 970 large numbers of both glucosinolate and cardenolide compounds. Alternatively, a full complement of 971 synthesis genes may be maintained in species' gene pools at low frequencies until they are favored by 972 a new environment, or they might be maintained in the genome but not expressed in leaves. More 973 extensive sampling within each species will be required to conclusively address this question, 974 although a preliminary screening of multiple E. cheiranthoides accessions suggests little to no 975 variation in glucosinolate profiles within this species (T. Züst, unpublished data). 976 Myrosinase activity levels differed among glucosinolate chemotypes, and activity was 977 positively correlated with glucosinolate abundance in plants when controlling for glucosinolate 978 chemotype. Erysimum species that predominantly produced indole glucosinolates or 4-methylsulfinyl 979 glucosinolates had negligible myrosinase activity against the assayed aliphatic glucosinolate sinigrin. 980 Indole glucosinolates can be activated by PEN2 -a thioglucosidase that is more specific for indole glucosinolates that are produced. In contrast to glucosinolate defenses, myrosinase activity was more 985 similar among related species, suggesting that the two defense components are subject to different 986 selective regimes, with the potential for maladaptive combinations between glucosinolate defense and 987 myrosinase activity. In addition, myrosinase activity is highly dependent on the presence of other 988 proteins and cofactors (Halkier and Gershenzon 2006), which may also differ between Erysimum 989

species. 990
We detected considerable amounts of the evolutionarily novel cardenolide defense in 47 out 991 of 48 Erysimum species or accessions. Among the 95 likely cardenolide compounds, there were 992 several structures that had not been described previously in Erysimum. This metabolic diversity had 993 three main sources: modification of the genin core structure, variation of the glycoside chain, or 994 isomeric variation (e.g., through the incorporation of different isomeric sugars). Structural variation in 995 cardenolides affects the relative inhibition of Na + /K + -ATPase (Dzimiri et al. 1987, Petschenka et al. 996 2018) and physiochemical properties such as lipophilicity, which play an important role in uptake and 997 metabolism of plant metabolites by insects (Duffey 1980). Individual Erysimum species produced 998 between 15 and 50 different cardenolide compounds, and the comparison of quantification by total 999 mass ion counts vs. quantification by inhibition of Na + /K + -ATPase revealed highly similar results.