- Research article
- Open Access
Transcriptome responses to different herbivores reveal differences in defense strategies between populations of Eruca sativa
BMC Genomics volume 20, Article number: 843 (2019)
Intraspecific variations among induced responses might lead to understanding of adaptive variations in defense strategies against insects. We employed RNA-Seq transcriptome screening to elucidate the molecular basis for phenotypic differences between two populations of Eruca sativa (Brassicaceae), in defense against larvae of the generalist and specialist insects, Spodoptera littoralis and Pieris brassicae, respectively. The E. sativa populations originated from desert and Mediterranean sites, where the plants grow in distinct habitats.
Responses to elicitation of the plants’ defenses against wounding and insect herbivory resulted in more upregulated transcripts in plants of the Mediterranean population than in those of the desert. PCA analysis differentiated between the two populations and between the elicitation treatments. Comprehensive analysis indicated that defense responses involved induction of the salicylic acid and jasmonic acid pathways in plants of the desert and Mediterranean populations, respectively. In general, the defense response involved upregulation of the aliphatic glucosinolates pathway in plants of the Mediterranean population, whereas herbivory caused downregulation of this pathway in desert plants. Further quantitative RT-PCR analysis indicated that defense response in the desert plants involved higher expression of nitrile-specifier protein (NSP) than in the Mediterranean plants, suggesting that in the desert plants glucosinolates breakdown products are directed to simple-nitriles rather than to the more toxic isothiocyanates. In addition, the defense response in plants of the desert population involved upregulation of flavonoid synthesis and sclerophylly.
The results indicated that differing defense responses in plants of the two populations are governed by different signaling cascades. We suggest that adaptive ecotypic differentiation in defense strategies could result from generalist and specialist herbivore pressures in the Mediterranean and desert populations, respectively. Moreover, the defense responses in plants of the desert habitat, which include upregulation of mechanical defenses, also could be associated with their dual role in defense against both biotic and abiotic stresses.
The term “induced defense” in plants refers to their ability to respond to herbivory by elevating their defense mechanisms, which rely mainly on the jasmonic acid signaling cascade and its interactions with other phytohormones, mainly salicylic acid and ethylene [1, 2]. Consequently, to optimize their defenses plants differ in their responses to different herbivores or to wounding in general, e.g., [3,4,5,6,7]. Such differential responses need the plant to recognise specific attacks, via mechanisms that mainly are associated with specific elicitors in the oral secretions of the chewing insect [8,9,10].
In members of the Brassicaceae, induced resistance against herbivory includes accumulation of glucosinolates, the main chemical defense metabolites, which provide effective defense against a wide range of herbivores [11, 12]. Mechanical damage to the leaf, caused by a herbivore, releases the enzyme myrosinase, which hydrolyzes intact glucosinolate molecules to their bioactive breakdown products: epithionitriles, nitriles, thiocyanate, and more-toxic isothiocyanates (ITC). However, specialist herbivores such as Pieris rapae and Plutella xylostella evolved mechanisms to suppress the release of the toxic products of glucosinolates breakdown and thereby to avoid the harmful effects of ITC [13,14,15]. Moreover, in Arabidopsis, it was shown that herbivory by the specialist P. rapae did not induce accumulation of aliphatic glucosinolates, as herbivory by the generalist Spodoptera exigua did . In support of this specialist/generalist paradigm, it also was shown that specialist and generalist chewing herbivores elicited different phytohormone responses, as in the case of Boechera divaricarpa (Brassicaceae) . However, contrary to this generalist/specialist paradigm, a microarray analysis of A. thaliana showed that plant responses to various specialist and generalist lepidopteran species could not be attributed to the insects’ degree of specialization .
It is assumed that the level of herbivory can lead to spatial intraspecific variation in defense against herbivores; more strongly defended plants are common in habitats where herbivores are more dominant, and vice versa . In consistent with this general concept, chemotypic variation in glucosinolates profiles among 75 populations of A. thaliana in Europe was strongly correlated with the abundance of the specialist aphids Brevicoryne brassicae and Lipaphis erysimi . Other studies have shown that genetic variations in defenses were associated with plants’ competitive ability [20, 21], which suggests that response to herbivores may also mediate allelopathic interactions [21, 22]. These studies strengthen the general view of the optimality theory  in that they stress the role of the tradeoff between growth and defense in shaping genetic variation of induced compared with constitutive defenses . In addition, because herbivores are expected to prosper in competitive plant communities characterized by high vegetation cover, variation in herbivore communities may represent a selective pressure that leads to intraspecific variations of defenses .
Previously we showed that in Israel, populations of Eruca sativa (Brassicaceae) originating from Mediterranean and desert habitats, differed in their induced defenses against specialist and generalist herbivores, i.e., P. brassicae and S. littoralis, respectively . Herbivory by larvae of the generalist and the specialist insects induced accumulation of glucosinolates in plants of the Mediterranean population but not in those of the desert population. Furthermore, our previous results suggested that trypsin proteinase inhibitor activity was involved in the induced defense response in plants of the desert population . To further test our hypothesis that populations of E. sativa exhibit ecotypic differentiation in defense strategies against herbivores, in the present study we applied RNA-Seq technology to analyse molecular patterns associated with herbivory in plants of the two populations. Most previous studies used transcriptome screening to examine the molecular basis of the responses of cruciferous species to various herbivores [3, 5,6,7]. In the present study a molecular comparison between plants of the two populations provided a more comprehensive understanding of the global transcript patterns and specific genes relevant to responses of plants to specific treatments (i.e., exposure to larvae of the generalist S. littoralis and specialist P. brassicae). It also enabled comparison between the effects of different environments which, in turn, can be linked to counter-adaptation of plant populations to different herbivores.
Plant growth conditions
Seeds of two populations, characteristic of desert (32° 04′ 49“ N, 35° 29’ 46” E; ≤ 200 mm annual rainfall) and Mediterranean (32° 46′ 39“ N, 35° 39’ 29” E; ≥430 mm annual rainfall) habitats, created under uniform conditions , were germinated on moistened Whatman No. 1 filter paper in 9-cm Petri dishes. The seeds were set to germinate in a growth chamber at 25 °C with an 8/16-h day/night photoperiod. Four-day-old seedlings were transferred to germination trays and placed in an insect-free net house; after 2 weeks the plants were transferred to 1-L pots containing a mixture of 50% peat, 30% tuff, and 20% perlite (Shacham, Israel) and were irrigated with an automatic irrigation system at 150 mL/day per pot. The experiment was conducted in February and March, with average max/min temperatures of 20/14 °C.
Elicitation and sampling
Defense mechanisms in the plants of the two investigated populations were induced by: 1) wounding with a pattern wheel; 2) wounding and application of 20 μL of oral secretions (OS) of S. littoralis and P. brassicae, diluted 1:5 (v:v) with distilled water. In the first group of plants, distilled water was applied to the wounded leaves (designated as “wounded”). Fivefold-diluted OS of S. littoralis was shown to cause transcript change in plants of A. thaliana . Similarly, preliminary experiments (not shown) and further quantitative RT-PCR tests [see Results section] showed that the diluted OS of both the generalist and the specialist larvae caused transcriptional changes as compared with control and wounded leaves. Treatment with OS enabled collection of leaf samples up to 48 h after elicitation without reduction of tissue area, as found when larvae were reared directly on the plants. To collect the OS, larvae of the generalist (S. littoralis) and specialist (P. brassicae) herbivores were reared from hatching to the third-instar stage on E. sativa plants, and their OS were collected with a vacuum system.
Five leaves on each plant were elicited, and were harvested at several time points after elicitation, from 0 (control) to 48 h. They were collected in chronological order starting at the first elicitation time point, and at each time point one leaf was harvested from every plant, all at the same development stage; an additional leaf was collected from each non-induced plant, as a control. The samples were immediately frozen in liquid nitrogen and were further used for RNA isolation.
RNA isolation and RNA-Seq analysis
Total RNA of each collected sample was extracted by using the TRI reagent (Sigma-Aldrich, Israel); 1 μL (2 u) of DNAse (Ambion, Thermo Fisher Scientific) was added to each sample to remove traces of DNA. RNA quality and integrity were verified with the Tapestation 2200 system (Agilent Technologies, USA). A 3-μg sample of RNA from each of the three leaves of a single plant was pooled to create one sample of the early (i.e., taken after 0.5–2 h) elicitation response (ER) to each treatment. Late elicitation responses (LR) for each plant included a pooled samples of RNA extracted from leaves collected 24 and 48 h after elicitation (4.5-μg for each time point). Samples from each population comprised a total of six biological replicates per treatment — three for each of the early- and late-induced responses; a single plant was used as a biological replicate. The control samples comprised three replicates. Thus, there was a total of 42 analysed samples. cDNA libraries were then prepared with the TruSeq Library Prep Kit V (Ilumina, Inc.). The samples were sequenced with the Illumina Hi-Seq 2000 system at the Technion Genome Center (Haifa, Israel).
De novo transcriptome assembly
Raw reads were subjected to a filtering and cleaning procedure as follow: First, the SortMeRNA tool was used to filter out rRNA ; then the FASTX Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html, version 0.0.13.2) was used for: (i) trimming read-end nucleotides with quality scores < 30 by using the fastq_quality_trimmer; (ii) removing reads with less than 70% of base pairs with quality score ≤ 30 by using the fastq quality filter. The total of ~ 1 T cleaned reads, obtained after processing and cleaning, were assembled de novo by using the Trinity software (version trinityrnaseq_r20140717 2.1.1) . The resulting de novo-assembly-generated transcriptome consisted of 80,946 contigs with N50 of 1081 bp.
The sequencing data were deposited in the NCBI Sequence Read Archive (SRA) database as bioproject PRJNA511735.
Sequence similarity and functional annotation
To assess the similarity of the transcriptome to those of other models and closely related species, sequence similarity was analyzed with the BLAST (Basic Local Alignment Search Tool) algorithm with an E-value cut-off of 10− 5 . The BLASTX algorithm was used to search protein databases by using a translated nucleotide query for comparison of the assembled contigs with sequences deposited in the Arabidopsis information resources (TAIR, http://www.arabidopsis.org). The transcriptome was used as a query for a search of the NCBI non-redundant (nr) protein database that was carried out with the DIAMOND program .
Differential expression and cluster analysis
The transcript was quantified (i.e., the number of reads per gene was determined) from RNA-Seq data by the using Bowtie aligner  and the Expectation-Maximization method (RSEM) . Differential expression was analysed with the edgeR software suit  and transcripts that were more than twofold differentially expressed with false-discovery-corrected statistical significance of ≤0.05 were considered differentially expressed . The expression patterns of the transcripts in different samples were studied by using cluster analysis of the differentially expressed transcripts in at least one pairwise sample comparison. Then, the Trinity protocol  was used to design expression normalization by using TMM (trimmed mean of M-values), following calculation of FPKM (fragments per feature kilobase per million). Hierarchical clustering of the normalized gene expression (by using centralized and log2 transformation ) and heat-map visualization were performed by using R Bioconductor . The VENNY tool  was used for construction of Venn diagrams. Principal component analysis (PCA) was applied with the FactoMineR package of R  to gain more insight on separation between samples of each treatment separately (based on the normalized expression of the average of three replications).
GO enrichment and pathways analysis
Gene ontology (GO) enrichment analysis was carried out by using the Plant MetGenMAP (http://bioinfo.bti.cornell.edu/cgi-bin/MetGenMAP/home.cgi) , based on the TAIR homology results. The tool enables integration of the functional categories between populations and among treatments with multiple testing correction of False Discovery Rate (FDR < 0.05) . Differentially expressed genes were displayed on diagrams of the Secondary Metabolism Map by using MapMan .
Quantitative RT-PCR was applied for further analysis of two myrosinase-associated proteins — the nitrile-specifier protein NSP2 and the epithiospecifier modifier ESM1 — as marker genes for glucosinolate breakdown . Reverse transcription with oligo dT (Fermentas Thermo) was used to synthesize cDNA from RNA samples (see RNA isolation and RNA-Seq analysis above). The cDNA samples were diluted to a uniform concentration (62 ng/μL) and qRT-PCR amplifications were performed with a Rotor-Gene 6000 instrument (Corbett-Qiagen, Valencia, CA, USA) by using components supplied in the KAPA SYBR FAST kit (Kapa Biosystems, Woburn, MA, USA), as previously described . The threshold cycle (Ct) was automatically determined with Rotor-Gene 6000 software and the relative expression levels of target genes were calculated with the aid of a ‘two-standard curve’ (i.e., that of the gene of interest and that of actin), implemented in the Rotor-Gene software. Each sample was analysed in two technical replicates for each target gene. Standard curves were created in each run by using a pooled cDNA sample; a reference cDNA calibration sample was used to normalize the multi-run results.
Quality trimming and filtration resulted in a total of ~ 1 T cleaned reads with an average of 22.3 M clean reads per sample. These were assembled by using Trinity, and generated 80,946 contigs for the transcriptome catalogue, with N50 of 1081 bp. Matching against the TAIR database yielded a significant hit of 54,552 contigs (67.4%). Annotating the transcriptome catalogue by aligning the contigs to the NCBI non-redundant (nr) protein database resulted in 61,403 out of 80,946 contigs (75.85%), with at least one DIAMOND hit to a protein. Of these contigs, 30,251 (~ 50%) matched sequences from the genomes of Brassica napus, followed by ~ 20% matching with Brassica oleracea var. oleracea (Additional file 1: Figure S1). A summary of the transcriptome catalogue, presenting the information of the full assembly and the number of raw reads, clean and mapped reads of each sample is provided in Additional file 4.
Principal component analysis (PCA) divided the transcriptome profiles mainly according to population and time after elicitation (Fig. 1 and Additional file 2: Figure S2). The first two axes of the PCA represent the differentiation between the two populations and between the controls and ER and LR of the three elicitation treatments (Fig. 1). Analysis of the number of differentially expressed genes (DEGs) derived for each elicitation treatment were determined according to their significance (FDR < 0.05; twofold) as compared with control non-elicited plants. Overall, the results revealed that the ER to wounding or to herbivory by S. littoralis and P. brassicae yielded more transcriptional changes than the LR (Fig. 2). Elicitation by wounding, S. littoralis and P. brassicae caused substantially more changes in the numbers of ER-upregulated DEGs in the Mediterranean plants than in those of the desert population — 1.7-, 1.8- and 1.3-fold, respectively.
Comparison of the transcript patterns between treatments: co-expressed and treatment-specific DEGs
Venn diagrams enabled clustering of the DEGs into up- and downregulated genes of overlapping and treatment-specific groups, relative to unelicited control samples (Fig. 3a). Many DEGs were co-expressed by the three elicitation treatments in plants of the desert and Mediterranean populations, both up- and downregulated DEGs. The percentage of upregulated DEGs specifically elicited by OS of S. littoralis was higher in the Mediterranean plants than in the desert ones: 18.7 and 14.0%, respectively. The responses to OS of the generalist herbivore also resulted in a higher percentage of downregulated DEGs in the Mediterranean plants than in the desert ones, at 29.5 and 12.8%, respectively. But the response of plants to elicitation by the specialist herbivore yielded higher percentages of non-overlapping up- and downregulated DEGs in the desert plants — 14.9 and 21.6%, respectively — than in the Mediterranean ones — 10.1 and 8.4%, respectively (Fig. 3a).
To improve our understanding of the differences between the responses of plants of the two populations to the elicitation treatments, we then analyzed gene ontology (GO) to cluster the core set of treatment-specific upregulated and downregulated DEGs to biological processes. Overall, in plants of both populations most of the differentially upregulated DEGs in the three elicitation treatments were categorized to metabolic or cellular processes (Fig. 3b). Nevertheless, following exposure to OS of the generalist herbivore, slightly more treatment-specific upregulated genes were associated with defense responses (categorized as ‘response to stress’) in the Mediterranean plants than in the desert ones: 8.8 and 7.9%, respectively (Fig. 3b).
Classification according to pathways categorized the total non-overlapping DEGs of plants of the Mediterranean and desert populations to a maximum of 159 upregulated pathways in Mediterranean plants and 95 in desert plants, in response to OS of S. littoralis and P. brassicae, respectively. The numbers of downregulated non-overlapping pathways were highest in the desert and Mediterranean plants following treatment with OS of the specialist and generalist herbivores, respectively: 203 and 262 pathways, respectively. The complete list of up- and downregulated pathways is presented in Additional file 5.
Figure 4 presents a heat-map of the P values of pathways that were significantly changed in at least one treatment; it shows that salicylate and phenylpropanoid biosynthesis were among the pathways in plants of the desert population that were significantly changed by all elicitation treatments. In addition, treatment with OS of both the generalist and specialist herbivores significantly induced the flavonoid (P < 0.01) and suberin (P < 0.05) biosynthesis pathways in plants of the desert population (Fig. 4a). The results also indicated that glucosinolate breakdown was significantly changed in response to OS of the specialist herbivore only in the desert plants (P = 0.05). In Mediterranean plants the mevalonate pathway was significantly upregulated by all elicitation treatments, and the jasmonic acid pathway was significantly changed in response to treatment with OS of both the generalist and the specialist herbivores. Treatment with OS of the generalist insect significantly changed the putrescine biosynthesis pathway in plants of both populations (Fig. 4a).
Most of the downregulated pathways encompassed primary metabolic processes: photosynthesis and sucrose biosynthesis, among others (Fig. 4b). Interestingly, the aliphatic glucosinolate biosynthesis pathway was significantly downregulated in plants of the desert population, in response to wounding and to elicitation by OS of S. littoralis (Fig. 4b).
Differences between populations of E. sativa, in their transcriptome profiles
Clustering in heat-maps of the significantly differentially expressed transcripts, which differentiated between the transcriptome profiles of the two populations and between ER and LR (Fig. 5), aided interpretation of the differences between the responses of plants of the two populations to the elicitation treatments. In general, the transcriptome profiles could be divided into two main groups of gene clusters: Transcripts in the first group differentiated similarly between early and late responses in plants of both populations (clade designated as ‘A’, Fig. 5a, b). The second group included transcripts that differentiated between the transcriptome profiles of plants of the two populations (i.e., clade B in Fig. 5a and b, clade C in Fig. 5b and clade D in Fig. 5c). In the second group: Clade B1 clustered transcripts that were upregulated in the desert plants but downregulated in the Mediterranean ones and vice versa for B2 (Fig. 5a, b); Clade C mostly clustered LR transcripts that were exclusively upregulated in response to elicitation by OS of S. littoralis in Mediterranean plants (Fig. 5b); Clade D differentiated between the transcriptome profiles of the ER and the LR of plants of both populations that were elicited with OS of the specialist herbivore (Fig. 5c).
Clade A, which differentiated between the upregulated ER and the downregulated LR, in plants of both populations, included 359 and 658 transcripts from leaves that were either wounded or treated with OS of the generalist herbivore, respectively — 24.1 and 32.1%, respectively, of the total (Fig. 5a, b; Additional file 6). Following wounding, more transcripts were clustered in clade B2 than in B1 (538 and 447, respectively) (Fig. 5a). Further classification according to pathways showed that following wounding, the jasmonic acid pathway (P < 0.001) and its conjugates (P = 0.024), the LOX-HPL cascade (P = 0.01) and the abscisic acid (P = 0.023) and flavonoid biosynthesis (P = 0.001) pathways were among the significantly changed pathways in clade A (Fig. 6a; Additional file 6). The salicylate biosynthesis and the mevalonate pathways were among the significantly changed pathways in both clades B1 and B2 (Fig. 6a).
The transcriptome profile of plants that were treated with OS of the generalist herbivore (Fig. 5b) included 1.6 times as many transcripts in clade B1 as in B2 — 539 and 337 genes, respectively — whereas clade C clustered 426 LR genes that were mostly upregulated in Mediterranean plants — 20.7% of the total processed hits (Additional file 6). Glucosinolate biosynthesis from phenylalanine (P = 0.05) was among the significantly changed pathways of clade A, together with several signaling pathways such as jasmonic acid (P < 0.001) and its conjugates (P = 0.037), the LOX-HPL pathway (P = 0.025), and the ethylene (P = 0.023), abscisic acid (P = 0.05) and cytokinins biosynthesis (P = 0.002) pathways. The salicylate (P < 0.001), suberin (P = 0.05) and phenylpropanoid synthesis (P = 0.005) pathways were among the significantly changed pathways of clade B1; the aliphatic glucosinolates (P = 0.04) and the mevalonate pathways (P = 0.004) were among the significantly changed pathways of clade B2 (Fig. 6b; and Additional file 6).
Similar numbers of transcripts of the transcriptome profile of plants that were treated with OS of P. brassicae were clustered in clades B1 and B2; they comprised 31.3 and 28.6%, respectively, of the total processed hits (Fig. 5c; Additional file 6). However, 3.7 times as many transcripts were clustered in clade D1 as in D2; this differentiated the ER from the LR transcriptome of plants of the two populations. These were classified into glucosinolate biosynthesis from phenyalanine (P = 0.04), flavonoid synthesis (P = 0.02), and the signaling cascades: jasmonic acid synthesis (P < 0.001) and its conjugates (P = 0.034), and the LOX-HPL (P = 0.021) and abscisic acid synthesis (P = 0.042) pathways. The salicylate biosynthesis (P = 0.04) and the mevalonate pathways (P = 0.007) were among the significantly changed pathways in clades B1 and B2, respectively (Fig. 6c; and Additional file 6).
A MapMan representation of the comparison between expression patterns of genes in plants of the two population (Additional file 3: Figure S3) illustrates the differences between responses to the three elicitation treatments, as well as differences between the responses of plants to the elicitation treatment. The results highlighted the defense roles of glucosinolates and terpenoids in Mediterranean plants, and that of flavonoids in desert plants. The results showed that OS of the specialist herbivore activated differing response in plants of the two populations; it induced genes in the alkaloids and carotenoids biosynthesis pathways in the Mediterranean plants and dihydroflavonols and anthocyanin in the desert plants (Additional file 3: Figure S3).
Expression of the genes encoding two myrosinase-associated proteins — NSP2 and ESM1 — were measured, because these reflect the division of GS-breakdown products into simple nitriles and ITC, respectively [39, 40]. The uninduced level (control) of NSP2 was significantly higher in the desert plants than in the Mediterranean ones (t-test, P < 0.01; Fig. 7). Moreover, significant induction of NSP2 expression was observed only in desert plants 24 h after treatment with the generalist OS (t-test, P < 0.05), so that there were significant differences in its expression between plants of the two populations (t-test, P < 0.05). Significantly higher expression levels in the desert plants than in the Mediterranean ones also were detected 6 and 48 h after wounding (t-test, P < 0.05 and P < 0.01, respectively), and 24 h after application of the specialist herbivore OS (t-test, P < 0.01; Fig. 7).
Wounding, as well as the application of OS of both herbivores, significantly increased the accumulation of ESM1 in plants of the two populations 6 h after treatment (t-test, P < 0.05). Nevertheless, 48 h after wounding, its expression in Mediterranean plants was significantly 2.5 times as high as that in the desert plants (t-test, P < 0.01; Fig. 7); similarly, significant (t-tests, P < 0.01) 2.2- and 2.5-fold differences were found 6 and 48 h, respectively, after application of the generalist OS to wounded leaves (Fig. 7). There were no differences between populations in the expression of ESM1 in response to treatment with the specialist OS.
One of the main findings of this study was that the jasmonic acid and salicylic acid signaling cascades were among the significantly changed pathways that differentiated between the desert and the Mediterranean populations. In the desert plants, OS of both the generalist and specialist herbivores induced the salicylic acid biosynthesis pathway, whereas in the Mediterranean plants herbivory significantly upregulated the jasmonic acid pathway (Figs. 4, 5 and 6). Thus, the differences between the two populations, in their global plant responses to herbivory, such as were demonstrated in the PCA analysis (Fig. 1), with respect to the numbers of down- and upregulated DEGs (Figs. 2 and 3) and to gene clustering (Fig. 5), included changes in the regulatory pathways that could have very important effects on plants’ physiology, growth, and defense responses [41, 42]. The negative crosstalk between the jasmonic and salicylic acid pathways can, therefore, account for the differences between the defense responses of the two populations (Figs. 4, 5, 6 and 7; discussed below).
Defense against insects in members of the Brassicaceae relies mainly on glucosinolates, among which aliphatic glucosinolates were shown to be effective against generalist herbivores . The effects of jasmonic and salicylic acids on the synthesis and accumulation of glucosinolates were shown in Arabidopsis e.g. [16, 43, 44]; they typified the antagonistic effects of these phytohormones, especially on gene expression and accumulation of aliphatic glucosinolates . Similarly, the aliphatic glucosinolates synthesis pathway was upregulated in plants of the Mediterranean population (Fig. 4a), whereas it was significantly downregulated by the three elicitation treatments in plants of the desert population (Fig. 4b). Furthermore, the RNA-Seq results (Fig. 4) and qRT-PCR analysis (Fig. 7) indicated that differences in responses to herbivory between populations of E. sativa included differences in the expression of glucosinolate breakdown genes: OS of the generalist herbivore caused higher expression of ESM1 in Mediterranean plants than in desert ones (Fig. 7). Since ESM was shown to repress nitrile formation in Arabidopsis in favor of release of more harmful ITC , the present results support our previous opinion  that accumulation of aliphatic glucosinolates and ITC breakdown products plays the primary direct defensive role against herbivores in plants of the Mediterranean population. The activation of a combination of defensive traits was evident in plants of the Mediterranean population, when all three elicitation treatments induced the mevalonate pathway (Figs. 4 and 6), which might indicate an additional indirect defense .
In general, responses of plants of the Mediterranean population to S. littoralis were substantially stronger than those of the desert plants (Figs. 3 and 5B). Defense responses of the desert plants included upregulation of synthesis of flavonoids, known for their role in direct defense against herbivores , and of mechanical defenses (suberin and lignin synthesis pathways) (Figs. 4 and 6) that strengthen the physical barrier of the cell wall. Moreover, the induction of putrescine biosynthesis in response to OS of the generalist insect (Fig. 4a) possibly indicated involvement of polyamines in wound healing , but also might suggest that the response in plants of the desert population involves sclerophylly. Thus, the role of flavonols and mechanical barriers in induced defenses in the desert plants (Fig. 5b) may be associated with their dual role in defense against biotic and abiotic stresses .
All three elicitation treatments induced higher expression of NSP2 in plants of the desert population than in Mediterranean plants (Fig. 7). It previously was shown that the release of nitriles provided crucifer plants with the ability to avoid herbivory by specialist herbivores, which employ ITC as recognition cues to detect their host plants [50, 51]. Our observations indicate that the specialist Plutella xylostela was more abundant in the desert monospecific habitat than in the heterogeneous Mediterranean plant community . Accordingly, it can be deduced that the direction of glucosinolate breakdown to producing simple nitriles is efficient against P. xylostela, whereas glucosinolates and their breakdown products ITC provide efficient defense against generalist herbivores in the Mediterranean site.
According to the generalist vs. specialist herbivore paradigm, the ability to perceive herbivory enables plants to respond differently to different herbivores. As a result, in the evolutionary battle against specialist herbivores, various plant defense mechanisms are induced to act against adapted insects . It previously was shown that accumulation of flavonols in A. thaliana increased plant resistance to larvae of the specialist P. brassicae . In the present study all the results indicate that both populations responded differently to the three defense-elicitation treatments (Figs. 1, 2, 3, 4, 5, 6 and 7; Additional file 3: Figure S3). Furthermore, elicitation with OS of P. brassicae resulted in the induced synthesis of simple coumarins in plants of both populations (clade D1; Figs. 5 and 6); activation of alkaloid and carotenoid synthesis in Mediterranean plants, and synthesis of dihydroflavones and anthocyanin in desert plants (Figs. 5 and 6; Additional file 3: Figure S3). That these pathways were not induced by OS of S. littoralis or by wounding (Additional file 3: Fig. S3) indicates the ability of plants of both populations to perceive different herbivores and to maximize their defenses and fitness in accordance with the type of attack. It is interesting to note that OS of the specialist herbivore induced genes of the aromatic glucosinolate pathway in plants of both populations (Fig. 6c). Since the metabolism of aromatic glucosinolate to benzyl cyanide increases the vulnerability of eggs to parasitic wasps , our results highlight the role of indirect defenses after an attack by specialist herbivores.
In summary, the result of the transcriptome screening indicated ecotypic differentiation between two defense strategies in populations of E. sativa originating from desert and Mediterranean habitats, respectively. The fact that the defense response is governed by differences between regulatory signaling pathways, i.e., jasmonic and salicylic acid pathways, might indicate global in-planta differentiation that can be associated with environmental adaptation. In Boechera sp., interspecific variation in induced defenses against herbivores was associated with differences in their adaptation to drought . Accordingly, the role of flavonols and mechanical defenses in induced defense in plants of the desert population may be associated with their dual role in defense against biotic and abiotic stresses . In addition, we suggest that differentiation between the defense responses in the two populations is directly linked to biotic selection agents . Thus, the direction of glucosinolate breakdown to production of simple nitriles, and activation of multiple alternative pathways provide efficient defense against both specialist and generalist insects in the desert habitat.
Availability of data and materials
All relevant data supporting our findings is provided in the article and supporting information. The sequencing data were deposited in the NCBI as detailed in the Methods section.
Differentially regulated genes
Epithio specifier modifier
Principal component analysis
quantitative real-time PCR
Nguyen D, Rieu I, Mariani C, van Dam NM. How plants handle multiple stresses: hormonal interactions underlying responses to abiotic stress and insect herbivory. Plant Mol Biol. 2016;91:727–40.
Schweiger R, Heise A-M, Persicke M, Müller C. Interactions between the jasmonic and salicylic acid pathway modulate the plant metabolome and affect herbivores of different feeding types. Plant Cell Environ. 2014;37:1574–85.
Consales F, Schweizer F, Erb M, Gouhier-Darimont C, Bodenhausen N, Bruessow F, Sobhy I, Reymond P. Insect oral secretions suppress wound-induced responses in Arabidopsis. J Exp Bot. 2011;63:727–37.
Duceppe M-O, Cloutier C, Michaud D: Wounding, insect chewing and phloem sap feeding differentially alter the leaf proteome of potato, Solanum tuberosum L Proteome Science 2012, 10:73.
Ogran A, Landau N, Hanin N, Levy M, Gafni Y, Barazani O. Intraspecific variation in defense against a generalist lepidopteran herbivore in populations of Eruca sativa (mill.). Ecology and Evolution. 2016;6:363–74.
Truong D-H, Nguyen HC, Bauwens J, Mazzucchelli G, Lognay G, Francis F. Plant defense in response to chewing insects: proteome analysis of Arabidopsis thaliana damaged by Plutella xylostella. J Plant Interact. 2018;13:30–6.
Vogel H, Kroymann J, Mitchell-Olds T. Different transcript patterns in response to specialist and generalist herbivores in the wild Arabidopsis relative Boechera divaricarpa. PLoS One. 2007;2.
Duran-Flores D, Heil M. Sources of specificity in plant damaged-self recognition. Curr Opin Plant Biol. 2016;32:77–87.
Felton GW, Tumlinson JH. Plant–insect dialogs: complex interactions at the plant–insect interface. Curr Opin Plant Biol. 2008;11:457–63.
Wu J, Baldwin IT. Herbivory-induced signalling in plants: perception and action. Plant Cell Environ. 2009;32:1161–74.
Burow M, Halkier BA. How does a plant orchestrate defense in time and space? Using glucosinolates in Arabidopsis as case study. Curr Opin Plant Biol. 2017;38:142–7.
Textor S, Gershenzon J. Herbivore induction of the glucosinolate-myrosinase defense system: major trends, biochemical bases and ecological significance. Phytochem Rev. 2009;8:149–70.
Wittstock U, Agerbirk N, Stauber EJ, Olsen CE, Hippler M, Mitchell-Olds T, Gershenzon J, Vogel H. Successful herbivore attack due to metabolic diversion of a plant chemical defense. Proc Natl Acad Sci U S A. 2004;101:4859–64.
Ratzka A, Vogel H, Kliebenstein DJ, Mitchell-Olds T, Kroymann J. Disarming the mustard oil bomb. Proc Natl Acad Sci U S A. 2002;99:11223–8.
Winde I, Wittstock U. Insect herbivore counteradaptations to the plant glucosinolate–myrosinase system. Phytochemistry. 2011;72:1566–75.
Mewis I, Tokuhisa JG, Schultz JC, Appel HM, Ulrichs C, Gershenzon J. Gene expression and glucosinolate accumulation in Arabidopsis thaliana in response to generalist and specialist herbivores of different feeding guilds and the role of defense signaling pathways. Phytochemistry. 2006;67:2450–62.
Bidart-Bouzat MG, Kliebenstein D. An ecological genomic approach challenging the paradigm of differential plant responses to specialist versus generalist insect herbivores. Oecologia. 2011;167:677–89.
Strauss SY, Rudgers JA, Lau JA, Irwin RE. Direct and ecological costs of resistance to herbivory. Trends Ecol Evol. 2002;17:278–85.
Zust T, Heichinger C, Grossniklaus U, Harrington R, Kliebenstein DJ, Turnbull LA. Natural enemies drive geographic variation in plant defenses. Science. 2012;338:116–9.
Kempel A, Schädler M, Chrobock T, Fischer M, van Kleunen M. Tradeoffs associated with constitutive and induced plant resistance against herbivory. Proc Natl Acad Sci U S A. 2011;108:5685–9.
Cipollini DF. Does competition magnify the fitness costs of induced responses in Arabidopsis thaliana? A manipulative approach. Oecologia. 2002;131:514–20.
Kong C, Hu F, Xu X. Allelopathic potential and chemical constituents of volatiles from Ageratum conyzoides under stress. J Chem Ecol. 2002;28:1173–82.
Karban R, Baldwin I. Induced responses to Herbivory. Chicago: University of Chicago Press; 2007.
Barazani O, Quaye M, Ohali S, Barzilai M, Kigel J. Photo-thermal regulation of seed germination in natural populations of Eruca sativa Miller (Brassicaceae). J Arid Environ. 2012;85:93–6.
Kopylova E, Noé L, Touzet H. SortMeRNA: Fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics. 2012;28:3211–7.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q. Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat Biotechnol. 2011;29:644.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2014;12:59.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Robinson MD, McCarthy DJ. Smyth GK: edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57:289–300.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5:R80.
Oliveros JC: VENNY. An interactive tool for comparing lists with Venn Diagrams. 2007, http://bioinfogp.cnb.csic.es/tools/venny/index.html.
Husson F, Lê S, Pagès J. Exploratory multivariate analysis by example using R: chapman and hall/CRC; 2017.
Joung J-G, Corbett AM, Fellman SM, Tieman DM, Klee HJ, Giovannoni JJ, Fei Z. Plant MetGenMAP: an integrative analysis system for plant systems biology. Plant Physiol. 2009;151:1758–68.
Usadel B, Poree F, Nagel A, Lohse M, Czedik-Eysenberg A, Stitt M. A guide to using MapMan to visualize and compare Omics data in plants: a case study in the crop species, maize. Plant Cell Environ. 2009;32:1211–29.
Wittstock U, Burow M. Tipping the scales - Specifier proteins in glucosinolate hydrolysis. IUBMB Life. 2007;59:744–51.
Wittstock U, Burow M: Glucosinolate breakdown in Arabidopsis: mechanism, regulation and biological significance The Arabidopsis Book, American Society of Plant Biologists 2010, 8:e0134.
Heidel-Fischer HM, Musser RO, Vogel H. Plant transcriptomic responses to herbivory. In: Voelckel C, Jander G, editors. Annual plant reviews. UK: Sussex; 2018. p. 155–96.
Heil M, Baldwin IT. Fitness costs of induced resistance: emerging experimental support for a slippery concept. Trends Plant Sci. 2002;7:61–7.
Kliebenstein DJ, Figuth A, Mitchell-Olds T. Genetic architecture of plastic methyl jasmonate responses in Arabidopsis thaliana. Genetics. 2002;161:1685–96.
Mikkelsen MD, Petersen BL, Glawischnig E, Jensen AB, Andreasson E, Halkier BA. Modulation of CYP79 genes and glucosinolate profiles in Arabidopsis by defense signaling pathways. Plant Physiol. 2003;131:298–308.
Zhang Z, Ober JA, Kliebenstein DJ. The gene controlling the quantitative trait locus EPITHIOSPECIFIER MODIFIER1 alters glucosinolate hydrolysis and insect resistance in Arabidopsis. Plant Cell. 2006;18:1524–36.
Paré PW, Tumlinson JH. Plant volatiles as a defense against insect herbivores. Plant Physiol. 1999;121:325–32.
Treutter D. Significance of flavonoids in plant resistance: a review. Environ Chem Lett. 2006;4:147.
Angelini R, Tisi A, Rea G, Chen MM, Botta M, Federico R, Cona A. Involvement of polyamine oxidase in wound healing. Plant Physiol. 2008;146:162–77.
Baldoni E, Genga A, Cominelli E. Plant MYB transcription factors: their role in drought response mechanisms. Int J Mol Sci. 2015;16:15811–51.
Mumm R, Burow M, Bukovinszkine’Kiss G, Kazantzidou E, Wittstock U, Dicke M, Gershenzon J. Formation of simple nitriles upon glucosinolate hydrolysis affects direct and indirect defense against the specialist herbivore, Pieris rapae. J Chem Ecol. 2008;34:1311–21.
Hopkins RJ, van Dam NM, van Loon JJA. Role of glucosinolates in insect-plant relationships and multitrophic interactions. Annu Rev Entomol. 2009;54:57–83.
Ogran A, Conner J, Agrawal AA, Barazani O: Evolution of phenotypic plasticity: genetic differentiation and additive genetic variation for induced plant defense in wild arugula Eruca sativa. J Evol Biol 2019, doi:https://doi.org/10.1111/jeb.13558.
Renwick JAA. The chemical world of crucivores: lures, treats and traps. Entomologia Experimentalis Et Applicata. 2002;104:35–42.
Onkokesung N, Reichelt M, van Doorn A, Schuurink RC, van Loon JJ, Dicke M. Modulation of flavonoid metabolites in Arabidopsis thaliana through overexpression of the MYB75 transcription factor: role of kaempferol-3, 7-dirhamnoside in resistance to the specialist insect herbivore Pieris brassicae. J Exp Bot. 2014;65:2203–17.
Fatouros NE, Broekgaarden C, Bukovinszkine'Kiss G, van Loon JJ, Mumm R, Huigens ME, Dicke M, Hilker M. Male-derived butterfly anti-aphrodisiac mediates induced indirect plant defense. Proc Natl Acad Sci U S A. 2008;105:10033–8.
Alsdurf JD, Ripley TJ, Matzner SL, Siemens DH. Drought-induced trans-generational tradeoff between stress tolerance and defence: consequences for range limits? AoB Plants. 2013;5.
We thank Prof. Murad Ghanim (Agricultural Research Organization, Israel) for providing larvae of S. littoralis and Mrs. Michal Barzilai for her technical assistance in the experiments.
This research was supported by the Israel Science Foundation (Grants No. 1222/10 and 2037/17). The funding body had no role in the design of the study, analysis and interpretation of data, and in writing the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Comparison of the resulting 80,946 contigs against the transcriptome catalogue.
Results of the principal component analysis (PCA) of the transcriptome profiles of the three replicates of the desert and Mediterranean (Med) populations. Control, early (ER) and late (LR) responses to the three elicitation treatments: wounding (Wo) or OS of S. littoralis (Sl) or P. brassicae (Pb).
MapMan illustrations of the upregulated defense pathways in plants of the Mediterranean (blue) and the desert (red) populations in response to elicitation by wounding (A), or OS of S. littoralis (B) and P. brassicae (C). The upper and lower rows represent the early and late responses, respectively.
Statistical summary of the transcriptome catalogue.
List of exclusive up- and downregulated pathways of the desert and Mediterranean plants. The list present the pathway and strength of P values of leaves that were elicited by wounding and by OS treatments.
Categorization to pathways of transcripts designated to different clades.
About this article
Cite this article
Ogran, A., Faigenboim, A. & Barazani, O. Transcriptome responses to different herbivores reveal differences in defense strategies between populations of Eruca sativa. BMC Genomics 20, 843 (2019). https://doi.org/10.1186/s12864-019-6217-9
- RNA Seq
- Jasmonic acid
- Salicylic acid
- Generalist vs. specialist insects