Slug Feeding Triggers Dynamic Metabolomic and Transcriptomic Responses Leading to Induced Resistance in Solanum dulcamara

Induced plant responses to insect herbivores are well studied, but we know very little about responses to gastropod feeding. We aim to identify the temporal dynamics of signaling- and defense-related plant responses after slug feeding in relation to induced resistance. We exposed Solanum dulcamara plants to feeding by the gray field slug (GFS; Deroceras reticulatum) for different periods and tested disks of local and systemic leaves in preference assays. Induced responses were analyzed using metabolomics and transcriptomics. GFS feeding induced local and systemic responses. Slug feeding for 72 h more strongly affected the plant metabolome than 24 h feeding. It increased the levels of a glycoalkaloid (solasonine), phenolamides, anthocyanins, and trypsin protease inhibitors as well as polyphenol oxidase activity. Phytohormone and transcriptome analyses revealed that jasmonic acid, abscisic acid and salicylic acid signaling were activated. GFS feeding upregulated more genes than that it downregulated. The response directly after feeding was more than five times higher than after an additional 24 h without feeding. Our research showed that GFS, like most chewing insects, triggers anti-herbivore defenses by activating defense signaling pathways, resulting in increased resistance to further slug feeding. Slug herbivory may therefore impact other herbivores in the community.

Induced plant responses to insect herbivores are well studied, but we know very little about responses to gastropod feeding. We aim to identify the temporal dynamics of signaling-and defense-related plant responses after slug feeding in relation to induced resistance. We exposed Solanum dulcamara plants to feeding by the gray field slug (GFS; Deroceras reticulatum) for different periods and tested disks of local and systemic leaves in preference assays. Induced responses were analyzed using metabolomics and transcriptomics. GFS feeding induced local and systemic responses. Slug feeding for 72 h more strongly affected the plant metabolome than 24 h feeding. It increased the levels of a glycoalkaloid (solasonine), phenolamides, anthocyanins, and trypsin protease inhibitors as well as polyphenol oxidase activity. Phytohormone and transcriptome analyses revealed that jasmonic acid, abscisic acid and salicylic acid signaling were activated. GFS feeding upregulated more genes than that it downregulated. The response directly after feeding was more than five times higher than after an additional 24 h without feeding. Our research showed that GFS, like most chewing insects, triggers anti-herbivore defenses by activating defense signaling pathways, resulting in increased resistance to further slug feeding. Slug herbivory may therefore impact other herbivores in the community.

INTRODUCTION
In addition to producing constitutively expressed defenses, plants can increase their defenses upon herbivory (Agrawal, 1998;Walling, 2000;Wittstock and Gershenzon, 2002;Howe and Jander, 2008). Herbivore-induced responses are generally observed in local, damaged, as well as in systemic, undamaged, organs (van Dam and Raaijmakers, 2006;Chung et al., 2013). Deployment of inducible defense mechanisms, in concert with reorganizing primary metabolism, allows the plant to optimize resource investment under variable herbivore pressure (Schwachtje and Baldwin, 2008;Orians et al., 2011;Steinbrenner et al., 2011;Zhou et al., 2015). Induced responses show temporal dynamics; they commonly start with a build-up phase after the first damage and may decline after the herbivore has left (relaxation; de Vos et al., 2005;Mathur et al., 2011;Mason et al., 2017). The time-lag and attenuation of induced responses are cost-saving strategies to limit resource investment in defense when the chances of recurring feeding damage or the benefits of defense production are low (Karban, 2011;Underwood, 2012;Backmann et al., 2019). Moreover, plants may tolerate transient herbivory, e.g., by reallocating resources to unattacked organs for later regrowth (van der Meijden et al., 1988;Schwachtje et al., 2006;Orians et al., 2011). Tolerance and relaxation often act in synchrony with resistance mechanisms to attain an optimal response (Strauss and Agrawal, 1999;Stowe et al., 2000;Núñez-Farfán et al., 2007;Stout, 2013).
Induced defense responses can be tailored to specific herbivores (Agrawal, 2000;van Zandt and Agrawal, 2004;Chung and Felton, 2011). Plants perceive herbivory by chemical signals associated with feeding damage and with the herbivore itself, such as the herbivore's saliva (Mithöfer and Boland, 2008;Bonaventure et al., 2011;Heil et al., 2012). These signals differ among herbivores, allowing plants to mount responses specific to the attacker (Voelckel and Baldwin, 2004;Basu et al., 2018). Herbivore-induced responses that follow lead to extensive molecular reprogramming (Leon et al., 2001;Maffei et al., 2007). This process is regulated by phytohormones, most importantly jasmonic acid (JA) and salicylic acid (SA; Erb et al., 2012). Crosstalk between the JA and SA signaling pathways results in specific responses to herbivores and pathogens (Lorenzo and Solano, 2005;Beckers and Spoel, 2006;Thaler et al., 2012;Schweiger et al., 2014). Other phytohormones, in particular abscisic acid (ABA) and ethylene (ET) are also involved in herbivore-induced defense responses (Nguyen et al., 2018).
While the mechanisms and consequences of insect-induced responses are intensively studied, similar studies addressing responses to gastropod feeding damage are scarce. It is likely that slug or snail feeding also elicits inducible defenses, because of the damage they inflict and the selective pressure they exert on plant populations (Strauss et al., 2009;Korell et al., 2016). The few studies analyzing gastropod-induced defenses show that they can increase defensive plant metabolites, including volatiles (Khan and Harborne, 1990;Viswanathan et al., 2005;Desurmont et al., 2016;Danner et al., 2017). Surprisingly, studies on the molecular regulation of gastropod-induced responses focussed on the effect of locomotion mucus, instead of actual feeding damage. These studies provided evidence that gastropod mucus applied to-artificially wounded-leaves can induce local and systemic responses (Orrock, 2013;Falk et al., 2014;Kästner et al., 2014;Meldau et al., 2014). Mucus-induced responses involve both JA-and SA-related defense signaling (Falk et al., 2014;Meldau et al., 2014). The mucus of the gray field slug (GFS; Deroceras reticulatum) contains SA , which may suppress JA-dependent defense responses via negative cross-talk (Beckers and Spoel, 2006;Thaler et al., 2012;Schweiger et al., 2014). In addition to gastropod-specific elicitors in locomotion mucus, also the feeding mode of gastropods is different from that of chewing insects. Gastropods possess a radula, a tongue-like chitinous structure covered with minute teeth, which serves to scrape off leaf material. Artificial damage mimicking the damage inflicted by chewing herbivores followed by mucus application probably would not yield the full array of responses elicited by natural gastropod feeding . It is therefore essential to subject plants to real feeding damage to understand the full spectrum of responses induced by gastropod feeding.
We recently showed that 72 h feeding by GFS induces resistance to slugs in damaged leaves of bittersweet nightshade (Solanum dulcamara; Calf et al., 2019). The magnitude of the induced resistance varied among plant genotypes. However, slugs are voracious and opportunistic generalist herbivores, which leave their host plant after feeding at night. During the day, they commonly hide in the litter layer, after which they move to the same or another host-plant the following night (see Supplementary Video in Kästner et al., 2014). These temporal dynamics of slug-induced responses are of special interest, because plants may relax the induction of defenses if they experience only one feeding bout (Baldwin, 1996). This strategy would be reflected in the expression profiles of sluginduced responses. More specifically, we expect that long-term exposure to feeding could elicit stronger responses than shortterm feeding, and that the induced response may relax after short-term feeding has stopped.
Here we experimentally assess the temporal and spatial dynamics of slug-induced responses in S. dulcamara, a common host for gastropods (Viswanathan et al., 2005;Lortzing et al., 2017). We specifically addressed the following questions: (1) Does GFS feeding induce GFS resistance in both local, damaged and systemic, undamaged leaves?; (2) What are the chemical, physiological and molecular mechanisms underlying induced responses and resistance to GFS feeding?; (3) Do induced responses and resistance depend on the duration of, and the time elapsed since, GFS feeding? In our first experiment, we exposed S. dulcamara plants to GFS feeding for 24 or 72 h and sampled local as well as systemic leaves at different time points thereafter. We used an eco-metabolomic approach  in which we combined untargeted metabolomic profiling of leaf samples with slug preference assays. This allowed us to assess locally and systemically induced metabolomic responses and their consequences for later arriving slugs. We further quantified the levels of polyphenol oxidase (PPO) activity and trypsin protease inhibitors (TPI). These are putative defenses that are commonly induced by insect herbivory on S. dulcamara (Viswanathan et al., 2005;Nguyen et al., 2016a;Lortzing et al., 2017). In addition to an increased production of defense compounds, herbivore feeding often reduces photosynthetic activity and induces leaf senescence (Wingler and Roitsch, 2008;Mellway et al., 2009;Zhu et al., 2017). For this reason, we measured chlorophyll and anthocyanin levels. In a second experiment, we exposed plants to GFS feeding for 24 h and performed microarray analyses on the transcriptome of local leaves directly at the end of the treatment period, as well as after an additional 24 h without feeding. In this experiment we also quantified the phytohormones JA, its bioactive isoleucine conjugate (JA-Ile), SA and ABA. This allowed us to assess the signaling events and induction processes directly after feeding by GFS and compare them to those after a 24 h relaxation period.

General Experimental Design
Two independent experiments were performed to assess (1) induced resistance and metabolomic responses in local and systemic leaves, and (2) transcriptomic and phytohormone responses after feeding by GFS in S. dulcamara. To grow plants we used seeds collected in native S. dulcamara populations (locations provided below) following the procedure in Zhang et al. (2016). All plants and slugs were cultured under the same conditions as in Calf et al. (2018). Plants in the feeding treatments were fitted with clip cages on the tip of one or two leaves. This prevented that the whole leaf was consumed. Clip cages were left empty (control plants) or received one adult GFS. Depending on the experiment, GFS were left on the plant for 24 or 72 h (Figure 1). When harvesting, we assessed that all plants in the slug treatments showed substantial feeding damage (estimated by eye, >1 cm 2 ).

Experiment 1
Seeds of a single S. dulcamara population in the Netherlands (Goeree: 51 • 49 23.8 N, 3 • 53 19.2 E, provided by the Radboud University Genebank, Nijmegen, Netherlands) were used to grow plants. At the start of the experiment, the plants were 26 days old and 30-51 cm tall. Each plant was randomly assigned to one of five treatments (n = 6 plants/treatment; Figure 1, top half). Leaf number 10 and 11 from the apex were selected to receive GFS. Three groups of plants were exposed to GFS feeding for 24 h, and harvested 24, 48, or 72 h after the slugs had been removed (Figure 1, upper panel). An additional group of plants was exposed to 72 h GFS feeding and harvested 24 h after slug removal. Control plants were fitted with empty clip cages. These were left on the plant until the last clip cages were removed from the other plants (Figure 1). The leaves of the different treatment groups were harvested within 1 h (Figure 1).
The part of the leaf outside the clip cage was used to produce leaf disks using a cork-borer (1.5 cm Ø). One leaf disk was punched from each leaf (n = 12 disks/treatment/position). Leaf veins were avoided while punching out the disks. The remaining leaf material was sampled in liquid nitrogen and stored at −80 • C until further processing. Disks and samples of leaf number 5 and 6 from the apex, which vascular tissue is fully connected to the treated leaves (Viswanathan and Thaler, 2004), were used to assess systemically induced responses.

Slug Preference Assays
The leaf disks were used in multiple choice slug preference assays to test the effect of GFS feeding on slug preference. One leaf disk of each treatment was offered to GFS in Petri dishes (n = 12; 5 disks per Petri dish) following the procedure described in Calf et al. (2019). Local and systemic leaves were tested separately.

FIGURE 1 | Treatment overview of two independent experiments in which
Solanum dulcamara plants were exposed to an undamaged control treatment (C; white box) or to feeding by the gray field slug (GFS; Deroceras reticulatum) for 24 h (gray box) or 72 h (black box). Depending on the experiment and treatment, leaf harvest (dashed line) took place 0, 24, 48, or 72 h after the end of the treatment period.

Untargeted Metabolomic Profiling
An untargeted metabolomics approach was used to assess the effect of GFS feeding on the metabolomic profiles of S. dulcamara leaves. Leaf samples obtained from experiment 1 (n = 6 samples/treatment/leaf position) were extracted and analyzed by liquid chromatography -Time of Flight -Mass Spectrometry (LC-qToF-MS) following the method described in Calf et al. (2019). For methodological details see Document S1. Local and systemic leaf samples were analyzed in different sample runs. This approach resulted in a dataset with 142 and 170 mass signal clusters for the local and systemic sample sets, respectively. These groups potentially represent single metabolites. The signal with the highest intensity in the majority of the samples was selected to represent the respective metabolite in later analyses.

Targeted Analyses of Induced Responses
The levels of total proteins, PPO activity, TPI activity, chlorophyll and anthocyanin in local and systemic leaves (n = 6 samples/treatment/leaf position) were spectrophotometrically quantified from leaf extracts following methods by Bradford (1976) total proteins, Bode et al. (2013) PPO, Thaler et al. (1996) TPI, Wintermans and Demots (1965) chlorophyll, Nakata and Ohme-Takagi (2014) Figure 1, bottom half). Leaf number 10 from the apex was exposed to feeding by GFS for 24 h or fitted with an empty clip cage for the same time (undamaged control). Leaf harvest took place either directly (0 h) after removing the slugs or after an additional 24 h without GFS (Figure 1, lower panel). The total leaf area that was covered by the clip cage was collected for gene expression analyses. The leaf material outside the clip cage was sampled separately for phytohormone quantification. Thus we avoided contamination with locomotion mucus, which may contain SA . All samples were flash frozen in liquid nitrogen and stored at −80 • C until further processing.

Transcriptomic Analyses
The transcriptomic response upon GFS feeding was assessed using microarray analyses. Total RNA extraction and purification were performed on ground fresh leaf material using the RNeasy R Plant Mini kit (Qiagen) and DNAse I (Fermentas, RNAse-free) following the manufacturer's instructions. Quality and quantity of the RNA were estimated using a NanoDrop 1000 device (Thermo Fisher Scientific, Waltham, MA, United States) and 1.5% agarose gel electrophoresis. Equal quantities of total RNA from four individuals of a single plant population at each time point were pooled, resulting in four total RNA samples per treatment (control or induced). The RNA concentration, integrity and purity were assessed by electrophoretic analysis with the RNA 6000 Pico Kit using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States) 1 . All samples had an RNA integrity number (RIN) between 6.6 and 7.2 and were hybridized on 8x60K Agilent microarrays. cDNA labeling, microarray hybridization, design were performed as described in Lortzing et al. (2017). The design and the experimental data of the microarray are available at the Gene Expression Omnibus of the National Centre for Biotechnology Information (NCBI; GEO Accession: GSE131208). RT-qPCR analyses were performed for technical validation of the microarrays (see Supplementary Data Sheet S1 and "Statistical Analyses" section).

Statistical Analyses
Absolute leaf disk consumption (mm 2 ) in the preference assays (experiment 1) was analyzed using non-parametric statistical methods from the R "stats" and "PMCMR" packages (Pohlert, 2014;R Core Team, 2016). Friedman's rank sum test was applied followed by Conover's post hoc test with Bonferroni correction for multiple comparisons to evaluate overall preference differences among treatments. The Petri dish number was used as grouping factor. The minimum consumption to accurately assess preference was 100 mm 2 (11.7% of the offered leaf material). One replicate with less damage (local leaves) was excluded, leaving 11 replicates in this group.
The metabolomic datasets from experiment 1 were analyzed using the online R-based tool MetaboAnalyst 4.0 (Chong et al., 2018). The effect of GFS feeding on metabolomic profiles was assessed with discriminant analyses using orthogonal projection models to latent structures (OPLS-DA, Wiklund et al., 2008;Worley and Powers, 2013). These models separate metabolomic variation that is correlated with the treatment from random variation among samples. Individual OPLS-DA models were built to compare the generalized log-transformed (glog) intensity values of all metabolites in the control treatment to each of the GFS feeding treatments. Predictive significance of the models was assessed based on 1000 permutations using cross-validated predictive ability (Q 2 ) as performance measure (Westerhuis et al., 2008;Worley and Powers, 2013). The correlation (p corr ) of the response of each metabolite in relation to the first predictive component (p) from the model was used as a proxy for metabolite induction by GFS feeding. The explained variance by the predictive component [R 2 X (p) ] was used as a proxy for the modeled effect size of the treatment. The local and systemic responses were tested separately. Metabolites that showed a high correlation with the model (p corr > 0.7) and significant regulation according to Student's t-test after Bonferroni correction for multiple comparisons (P adj < 0.05), were selected as metabolites of interest. The relative levels of each metabolite of interest to the maximum mean value of this metabolite across all treatments was used to produce a heat map (%).
Microarray data were analyzed using the "limma" software packages from Bioconductor in "R" (Raj et al., 2006;R Core Team, 2016). Quantile normalized reads of 21261 out of 62970 contigs were included for data analysis after setting a critical level of quantification based on dark corner intensity (90% percentiles of non-labeled hairpin DNA probes) as well as removal of structural spots, repeatedly spotted probes, averaging repeatedly spotted contigs and selecting those contigs that were quantified in all samples of at least one treatment.
Average fluorescence values of contigs were log 2 transformed and visualized by principal component analyses (PCA) using the "prcomp" function. Data were fitted to a linear model using the "lmFit" function using dual contrasts (induced response at 0 or 24 h after exposure to GFS feeding versus the respective control). Contigs that showed an absolute log 2 -fold difference in expression >1 and significance (P adj ) <0.05 after correction for false discovery rate according to the Benjamini-Hochberg method were considered significantly regulated. Contig expression levels in the microarray were validated based on Pearson's correlation between fold-change in expression as estimated by microarray and RT-qPCR analyses for 10 genes related to primary and secondary metabolism (PR1a, OPR3, PI1, KPI4, LOXD, PPOA, ERF4, PPO, NIM1 and CS, R 2 = 0.97, Supplementary Figure S1 and Supplementary Data Sheet S1). Gene primers for qPCR were designed using the primer NCBI BLAST-tool 2 (see Supplementary Table S2).
Functional descriptions of regulated contigs were obtained from the Sol Genomics Network (SGN 3 ) and based on homology with tomato (D' Agostino et al., 2013). Gene ontology (GO) enrichment analyses for overrepresentation of molecular functions and biological processes were performed with the "TOPGO" package (Alexa and Rahnenfuhrer, 2010). GO annotation was based on Nguyen et al. (2016a). GO enrichment was assessed by comparing the distribution of the list of differentially expressed contigs to the distribution of all targets included in the data analysis. We used the "elim" algorithm, which corrects for the annotation of contigs in multiple parental GO terms (Alexa et al., 2006). Only GO terms that included a minimum of 50 annotated contigs were included for analyses. The P-values for the enrichment of each GO term are based on Fisher's exact tests.
The effect of GFS feeding on levels of total proteins, PPO activity, TPI activity, anthocyanin and chlorophyll from (Experiment 1) as well as phytohormone levels (Experiment 2) were statistically tested using parametric statistical methods from the "stats, " "car" and "agricolae" software packages in "R" (Fox and Weisberg, 2011;R Core Team, 2016;de Mendiburu, 2017). Data were log 10 or square root transformed when assumptions for homogeneity of variances among treatments and/or normal distribution of residuals were not met according to Levene's and Shapiro-Wilk normality test, respectively. Treatment effects on metabolites and enzymes were tested using one-way ANOVA, followed by Tukey's post hoc test for differences between treatments. Phytohormone levels were statistically tested using paired t-tests, including plant population of origin as paired factor. FIGURE 2 | Mean consumed leaf area relative to the total consumed area (% ± SE, n = 12) on Solanum dulcamara leaf disks in preference assays with the gray field slug (GFS; Deroceras reticulatum). Slugs were offered leaf disks from plants that were exposed to an undamaged control treatment (C; white bar) or to feeding by GFS for 24 h (gray bars) or 72 h (black bar). Leaf disks were collected 24, 48, or 72 h after the end of the treatment period (see also Figure 1). Disks of either local or systemic leaves were offered in independent preference assays. Different letters over the bars indicate significant preference differences (P < 0.05) according to Conover post hoc test of Friedman ANOVA (χ 2 local = 17.065, P local = 0.002; χ 2 systemic = 11.661, P systemic = 0.022) after Bonferroni correction of P-values for multiple comparisons.

Induced Effect on GFS Feeding Preference (Experiment 1)
Slugs preferred to feed on leaf disks of plants that were left undamaged (control treatment) compared with those of plants exposed to GFS feeding (Figure 2). This preference was significant for both local and systemic leaves. A 24 h feeding period was sufficient to reduce slug preference on local leaves; 72 h continuous exposure to GFS resulted in the highest resistance. In contrast to the local leaves, systemically induced leaves were equally resistant, without a significant effect of the duration or timing of feeding.

Metabolomic Profiling (Experiment 1)
OPLS-DA models were built to assess the effect of different temporal patterns of GFS feeding on leaf metabolomes ( Table 1). All models, except one (systemic 24 h/24h; Q 2 : 48.8, P < 0.1, Table 1) had a significant predictive ability (Q 2 : 60.6-83.9%). The effect sizes [R 2 X (p) ] of models testing the local response were always greater than the respective models testing the systemic Leaves were harvested at 24, 48, or 72 h past the end of the treatment period (see also Figure 1). Model codes indicate the duration of GFS feeding exposure and the time between the end of GFS exposure and harvest (duration/harvest). Separate models were built for local and systemic leaves. The cross-validated predictive ability (Q 2 ) as well as explained percentage of variance by the predictive component [R 2 X (p) ] are provided for each model. Symbols indicate significance of the crossvalidated model based on 1000 permutations (**P < 0.01, *P < 0.05, + P < 0.10).
response ( Table 1). In both leaf positions, the effect size of the model increased with time after 24 h of exposure to GFS; 72 h continuous exposure to slugs elicited the strongest response. We selected metabolites correlating strongly with the predictive component of each model (p corr > 0.7) and with significant regulation (P adj < 0.05, Supplementary Table S3). Moreover, we focused on the induction pattern of those metabolites which levels were consistently regulated in at least two treatments, or in both local and systemic leaves within a treatment ( Figure 3A).
Fifteen metabolites consistently changed in abundance in at least two treatments ( Figure 3A; 10 local, 5 systemic). In local leaves, most selected metabolites (7 out of 10) gradually changed either up or downward with time after 24 h of GFS feeding. The strongest response was observed upon 72 h of continuous feeding (Figure 3A, Local). Two of these gradually changing metabolites reduced in abundance upon GFS feeding (L119 and L088), whereas five increased (L026, L005, L053, L050, and L007). Among the gradually increased metabolites there are two putatively identified phenolamides (see Calf et al., 2019) N-caffeoylputrescine and a related metabolite (L025 and L026 in Figure 3A, see also Supplementary Figure S2). In addition to the gradually regulated metabolites, one metabolite was strongly reduced in local leaves in all feeding treatments (L097, <20% of the control level). Two other metabolites showed the strongest response upon 24 h of exposure to GFS (L099 and L024). The changes in systemic leaves were more erratic. The gradual response that was observed for most locally regulated metabolites was only seen for two metabolites in systemic leaves (Figure 3A, systemic). One of these metabolites increased (S129) and the other (S116) decreased in abundance.
Two metabolites were found to be significantly regulated both in local and systemic leaves at one single time point only (Figure 3A, Both). One of these metabolites was the glycoalkaloid solasonine (LGA3/SGA3 in Figure 3A, identified by an authentic standard in Calf et al., 2018) which increased most strongly in abundance after 72 h continuous exposure to slug feeding (see also Supplementary Figure S2). In the 24 h feeding treatments, solasonine levels gradually increase with time after slug removal. The levels of an unidentified metabolite were also upregulated in both leaf positions (L124/S147, m/z 377.15 in Figure 3A), but only 72 h after 24 h of exposure to feeding. The relative increase of this metabolite was small, as constitutive levels (C) were high to start with (∼60% of the induced level).

Targeted Analyses of Proteins and Compounds (Experiment 1)
We assessed the effect of GFS feeding on well-known insect inducible defensive proteins and compounds (Table 2 and Figure 3B). The levels of PPO activity, TPI activity and anthocyanins significantly increased in response to 72 h continuous exposure to feeding (Figure 3B), but only in local leaves ( Table 2). A 24 h feeding period resulted in intermediate levels of these compounds. Total protein and chlorophyll levels were not affected by slug feeding ( Table 2).

Transcriptomic Analyses (Experiment 2)
To assess the early molecular signaling mechanisms preceding the observed metabolomic changes, we analyzed transcriptomes of plants that were sampled directly (0 h) after 24 h of exposure to feeding as well as after an additional 24 h without feeding (Figure 1, lower half). Principal component analyses were performed on the expression values of all quantified contigs (21261) to assess differences in the transcriptomic profiles ( Figure 4A). The first principal component (PC1) captured the effect of GFS feeding, which accounted for 32.1% of the total variance among samples. Samples taken directly (0 h) after exposure to feeding separated further from the control samples than samples taken after an additional 24 h without feeding. The second principal component (PC2, 14.9%) revealed that transcriptomic profiles differed among plant populations. The German populations (Siethen and Erkner) separated from the Dutch populations (Goeree and Friesland).
A total of 1526 contigs (7.2% of the total set) was significantly up-or downregulated in response to feeding by GFS at either one of both time points (Figure 4B). The total number of regulated contigs directly (0 h) upon 24 h of exposure to feeding (1424) was >5 times greater than after an additional 24 h without feeding (262, Figure 4B). Overall, more contigs were upregulated than downregulated upon slug feeding. In the direct response treatment (0 h), the number of upregulated contigs (838) was twice as high as the number of downregulated contigs (426), whereas about an eightfold difference was observed after an additional 24 h without feeding (90 up, 12 down, Figure 4B). Two PPO enzymes (c2630 and c16387), five protease inhibitors (c460, c673, c1119, c4199, and c22651) and an anthocyanidin synthase (c10758) were among the 25 strongest regulated contigs at either one of both time points (Supplementary Table S4, 11-208 fold-change in expression). FIGURE 3 | Induced responses in Solanum dulcamara leaves after feeding by the gray field slug (GFS, Deroceras reticulatum). Plants (n = 6) were left undamaged (control treatment: C, white) or exposed to feeding by GFS for 24 h (gray) or 72 h (black). Samples were collected 24, 48, or 72 h after the end of the treatment period (see also Figure 1). (A) Mean relative abundance of metabolites of interest, the levels of which were affected by GFS feeding in at least two treatments or in both local and systemic leaves within one treatment. The metabolite ID for local (L) and systemic (S) leaves is composed of a unique number and the quantified mass (m/z). (B) Mean relative levels (±SE) of polyphenol oxidase (PPO) activity, trypsin protease inhibitor (TPI) activity and total anthocyanins in local leaves. Treatment effects were significant according to one-way ANOVA ( Table 2). Different letters over the bars indicate significant differences among treatments according to Tukey post hoc test (P < 0.05).
Enrichment analyses of gene ontology terms helped to assess the biological processes that were overrepresented at either one or both time points. The 77 (out of 339) GO terms that were most significantly enriched upon slug feeding (P adj < 0.0001) were selected and grouped according to the overall biological process they are involved in (group 1-6 in Table 3). In addition to the most strongly enriched terms, we selected six GO terms related to photosynthesis (group 7 in Table 3), which were previously found to be regulated upon insect feeding Nguyen et al., 2018). The total number of up-and downregulated contigs in each term as well as the specific numbers at each time point were further evaluated ( Figure 4C and Table 3). Most contigs in the selected GO terms were upregulated in response to slug feeding ( Figure 4C). The numbers of regulated contigs were always higher directly (0 h) after 24 h of exposure to slug feeding than after an additional 24 h without feeding ( Figure 4C; blue bar higher than red bar). Approximately 10-16% of all the contigs in terms related to defense responses (group 1) were significantly regulated. This group included various terms that involve responses to wounding, fungi, bacteria or nematodes ( Table 3). The many responses related to abiotic stimuli (group 2) indicate that there is large overlap in responses to biotic and abiotic stressors. Various terms related to defensesignaling phytohormones (group 3) were regulated upon slug feeding. These included JA (4 terms), SA (3 terms), ABA (2 terms) as well as single terms related to auxin, karrikin, gibberellin, ethylene and brassinosteroids ( Table 3). The highest percentage of phytohormone-related contigs was regulated for Symbols indicate significance at the following levels: ***P < 0.001, *P < 0.05. Significant treatment effects are shown in Figure 3B.
JA metabolic (GO:31) and biosynthetic processes (GO:33, both about 23%, Figure 4C and Table 3). The overall largest percentage of regulated contigs was observed for terms involved in organic compound metabolism (up to 35%, group 4). This group includes biosynthetic processes for many defense-related compounds, such as flavonoids, phenylpropanoids, coumarins and alkaloids. Also terms related to primary metabolism, such as phenylalanine and tyrosine metabolism, were enriched (group 4). Genes involved in cell wall biogenesis (group 5) are likely involved in wound healing. The upregulation of many basic cellular metabolic and transport processes (group 6) indicates extensive reprogramming of the plants' primary metabolism upon slug feeding. In contrast, only relatively few contigs involved in photosynthesis were regulated (group 7). Most genes in this category were upregulated, but only term 81 (chlorophyll biosynthetic process) was significantly enriched (P adj = 0.007; Figure 4C).

Phytohormone Analyses (Experiment 2)
Jasmonic acid as well as the levels of its bioactive isoleucine conjugate (JA-Ile) were elevated at both time points, with the strongest increase measured directly (0 h) after 24 h of exposure to GFS feeding ( Figure 5). Salicylic acid levels were not affected, whereas ABA levels were elevated directly (0 h) after 24 h of exposure to slug feeding, but not after an additional 24 h without feeding.

DISCUSSION
Our study revealed that feeding by gray field slugs (GFS) elicits jasmonic acid-based defense responses in S. dulcamara resulting in local and systemic induced resistance to conspecifics (Karban and Baldwin, 1997). Twenty-four hours of slug feeding sufficed to increase the levels of JA, JA-Ile and ABA, as well as significant metabolomic and transcriptomic changes. GFS feeding triggered transcription of genes related to diverse pathways, including JA and SA phytohormone signaling. Plants responded to GFS feeding by increasing commonly known defenses, such as solasonin, TPIs, PPOs, phenolamides and anthocyanins. This indicates that GFS feeding induces similar signaling pathways and defense metabolites as chewing insects. Other than insect feeding, slug herbivory did not affect chlorophyll levels or downregulate the transcription of photosynthesis-related genes. Whereas the magnitude of the transcriptomic response was reduced within 24 h after the slugs were removed, the metabolomic response kept increasing with time. The longest feeding period of 72 h resulted in the strongest metabolomic response and the highest resistance level. Taken together, this suggests that S. dulcamara limits resource investments in GFSinduced defenses while conserving primary metabolism in the absence of further slug feeding. Slug feeding triggered several phytohormones involved in defense signaling; JA, JA-Ile and ABA levels were significantly elevated upon slug feeding, which was also reflected by the transcriptomic response in the respective GO categories. Remarkably, SA levels were unaffected, whereas genes involved in SA metabolism and signaling were upregulated. This may be explained by the earlier finding that GFS locomotion mucus contains SA . In this study, application of slug mucus to Arabidopsis thaliana leaves induced the local expression of PR1, a marker gene for SA defense signaling . SA and JA both are involved in responses to herbivore feeding whereby they often act in negative crosstalk (Beckers and Spoel, 2006;Thaler et al., 2012;Schweiger et al., 2014). This triggered the speculation that GFS may benefit from the SA in the mucus, because it suppresses JA herbivore defenses . We did not test whether the locomotion mucus of the field-collected slugs in our study contained SA, but the expression of a S. dulcamara PR1 homolog was unaffected by GFS feeding (c374 in Supplementary Figure S1). Moreover, JA and JA-Ile levels, as well as jasmonic acid related gene expression, increased strongly upon slug feeding in S. dulcamara. If the mucus of our GFS would have contained SA, the levels were not sufficient to antagonize the JA-related induced responses or prevent induced resistance to conspecific slugs.
Both in the transcriptome and the metabolome we identified different responses that may cause slug resistance. Glycoalkaloids likely play a key role in resistance to GFS. Solasonine levels were significantly higher in both local and systemic leaves upon 72 h of exposure to slug feeding, and leaf disks of these plants were also the least preferred by the GFS. Moreover, we previously showed that constitutive resistance to GFS in S. dulcamara is associated with high levels of glycoalkaloids (Calf et al., 2018). Taken together, this suggests that glycoalkaloids cause both constitutive and induced resistance to slugs. In a common garden experiment with different S. dulcamara glycoalkaloid chemotypes we found that gastropods and specialist flea beetles show opposite feeding preferences for glycoalkaloid chemotypes . This implies that slug-induced increases in glycoalkaloid levels may have opposing effects on the many members of the natural S. dulcamara herbivore community (Calf and van Dam, 2012).
In addition to glycoalkaloids, we found several other induced defenses that possibly affect slug preference. Genes related to TPI and PPO activity were among the most strongly regulated genes after 24 h exposure to GFS feeding. This is in line with their increased activity in local leaves, in particular after 72 h continuous exposure to slugs. Previous studies found that PI levels and PPO activity as well as transcription of genes coding FIGURE 4 | Transcriptomic profiles in Solanum dulcamara leaves upon 24 h of exposure to feeding by the gray field slug (Deroceras reticulatum). Samples were collected directly (0 h; blue) or 24 h (red) after the end of the treatment period (see also Figure 1). (A) Principal component analysis of overall transcriptomic differences among samples (n = 4/treatment/population). Ellipses show 95% confidence intervals. Letters indicate the plant population of origin; E, Erkner; S, Siethen; G, Goeree; F, Friesland. (B) Numbers of significantly up-and downregulated contigs that are unique for, or shared among harvest time points (fold change >2, P adj < 0.05, n = 4). The total number of regulated contigs at each time point is given in brackets. (C) Percentage of contigs that was significantly up-or downregulated in 83 gene ontology (GO) terms. GO terms were grouped according to overall involvement in biological processes and sorted based on the number of contigs in each term (group numbers correspond to the numbers in Table 3).  GO term numbers (#) correspond to the number in Figure 4C. GO terms were grouped according to overall involvement in biological processes and sorted based on the number of contigs in each term. The number (n) of regulated contigs is given for the total response as well as separate for two time points (0 and 24 h) after the end of the treatment. All terms in group 1-6 were significantly regulated (P adj < 0.0001) according to Fisher's exact test after correction for false discovery rate using the Bonferroni method. Symbols in group 7 indicate significance on the following level: **P adj < 0.01.
for these proteins increased upon insect feeding or oviposition on S. dulcamara (Boege and Marquis, 2005;Nguyen et al., 2016b;Geuss et al., 2017;Lortzing et al., 2017). Increases in PI and PPO levels correlated with increased resistance to several insect herbivores and reduced egg hatching (Geuss et al., 2018;Nguyen et al., 2018). Increased PI levels can reduce slug feeding. Leaf damage by slug feeding was reduced by 50% in Arabidopsis plants genetically modified to overexpress a cysteine protease inhibitor from rize, oryzacystatin (Walker et al., 1999). However, it is not fully clear yet whether TPIs, which are serine protease inhibitors, also confer resistance to GFS, which mainly produces cysteine proteases in its digestive glands (Walker et al., 1998). Although cystein PI activity was not measured from our leaf extracts, a contig for cystein PI8 (c4199 in Supplementary Table S4) was FIGURE 5 | Levels of phytohormones in Solanum dulcamara leaves upon 24 h of exposure to feeding by the gray field slug (Deroceras reticulatum).
(A) jasmonic acid, (B) jasmonic acid-isoleucine, (C) salicylic acid, and (D) abscisic acid. Samples were collected directly (0 h) or 24 h after the end of the treatment period (see also Figure 1). Bars represent mean levels (±SE, n = 4) in undamaged control (white bars) and damaged leaves (gray bars). The symbols indicate significant differences between treatments on a single time point according to Welch's t-test on the following levels: **P < 0.01, *P < 0.05, + P < 0.1.
among the 25 most strongly regulated contigs. Because TPIs are effective resistance factors to insect herbivores, in particular lepidopteran larvae (Glawe et al., 2003;Bode et al., 2013;Zhu-Salzman and Zeng, 2015) it is likely that slug-induced responses affect other herbivores in the community. The transcriptomic and metabolomic analyses showed that anthocyanins and phenolamides were also induced. Each of these metabolites, or genes coding for their biosynthesis, were previously found to be induced upon insect feeding or pathogen infection as well as primed by oviposition (Nguyen et al., 2016b;Geuss et al., 2018). Anthocyanins have multiple roles in plants, including resistance to insect herbivores (Gould, 2004). Phenolamines, such as caffeoylputrescine, are intermediates in lignin biosynthesis (Gaquerel et al., 2014). The genes involved in this pathway dynamically respond to biotic stress such as herbivore feeding. Increased resistance may be conferred by one of the many compounds in the pathway or by increased leaf toughness due to an increase in lignin concentration. The strong transcriptomic response in the GO related to phenylpropanoids (for example GO:0009698 phenylpropanoid metabolic process; GO:0009805 coumarin biosynthetic process and GO:0009808 lignin metabolic process) indicates that slug and insect feeding trigger similarly strong responses in the phenolamide pathway. Whether or not this is indeed a functional response leading to reduced herbivore damage, as suggested by Gaquerel et al. (2014) remains to be determined.
Interestingly, we observed that, overall, slug feeding upregulated transcriptional responses in the primary metabolism. The many transcriptomic responses related to abiotic stimuli indicate that there is large overlap in responses to biotic and abiotic stress factors. Similar results were reported in an experimental analysis on interactive responses to drought, flooding and insect feeding in S. dulcamara (Nguyen et al., 2018). Such cross-talks between signaling pathways allow the plant to orchestrate both resource acquisition and defense allocation in complex natural environments (Nguyen et al., 2016b). Chewing insect herbivore damage, for example, can downregulate photosynthesis-related processes in S. dulcamara Nguyen et al., 2018) and other plant species (Zangerl et al., 2002;Tang et al., 2006;Lawrence et al., 2008;Bilgin et al., 2010). However, local and systemic herbivore-induced responses on the level of photosynthetic pathways cover the whole range from positive to negative gene regulation (Zhou et al., 2015). This led to the contrasting hypotheses that upregulation of photosynthetic capacity is a way to increase carbon for defense production, whereas downregulation of the costly photosynthetic machinery would free up resources for defense related pathways (Zhou et al., 2015). In our current study, only a few genes involved in photosynthesis were significantly regulated and chlorophyll levels did not change. These results suggest that the photosynthetic capacity of the leaves was maintained. Taken together, GFS-induced resistance in S. dulcamara may not come at a cost in terms of reduced primary metabolism (Zhou et al., 2015).
The observed temporal and spatial patterns suggest that S. dulcamara may optimize its response to GFS feeding. The strength of the induced response depended both on the period the slug was on the plant and the time that past after the slug was removed. The transcriptomic response was strongly relaxed in the absence of further damage. This is in line with our observation that the metabolomic response after 72 h of feeding is stronger than after a 24 h feeding bout. In contrast to the transcriptomic response, the metabolomic response to 24 h feeding gained in strength with time after slug removal. The observed metabolomic induction pattern is in line with the observed effect on slug preference in local leaves. In systemic leaves, maximal resistance was already achieved after 24 h feeding. This result suggests that young systemic leaves more rapidly attain effective levels of defense metabolites, such as glycoalkaloids (see above), than older local leaves, possibly via rapid translocation of glycoalkaloids and other defenses (van Dam et al., 1995). This is in agreement with the optimal defense hypothesis, in which younger leaves are more valuable to the plant, have higher constitutive defense levels, and are more inducible than older leaves .
Finally, our experimental set-up allowed us to reveal geographic variation among populations for slug-induced responses. The transcriptome profiles of S. dulcamara plants clearly showed a core set of GFS induced responses, but German and Dutch populations separated on the second axis of the PCA (Figure 4A). Similar patterns were found in a panel of eight A. thaliana accessions subjected to flooding (van Veen et al., 2016). Despite the genotypic differences, they identified a core set of flooding specific transcripts. In order to separate general slug-induced responses from genotypic or population specific responses, a wider set of S. dulcamara accessions should be screened.
In conclusion, our results illustrate that induced responses and resistance to GFS feeding in S. dulcamara are regulated by similar pathways and metabolites as those involved in insect induced responses Geuss et al., 2018;Nguyen et al., 2018). Yet, plants did not downregulate the transcription of photosynthesis-related genes upon slug feeding, which is in contrast to what commonly happens after insect feeding. This illustrates that gastropods and insects may induce specific responses, despite of similarities in the regulatory pathways they induce. Highly standardized experiments comparing the response of S. dulcamara to various gastropod and insect herbivores are required to assess in how far induction patterns are specific to different classes of herbivores. Besides, we show that relaxation plays a role in the dynamics of GFSinduced defenses, and that sustained feeding can boost the locally induced response. The temporal and spatial dynamics of GFS-induced responses, in particular with respect to those related to primary and secondary metabolism, illustrate that S. dulcamara plants possibly balance resource investment after induction. Together with the rapidly induced local and systemic resistance, we may conclude that S. dulcamara shows a functional response to GFS, which is well-balanced with the incidence of feeding damage. Depending on the turnover rate of defense compounds in the absence of feeding, this induced response may have consequences for other herbivores in the community (Viswanathan et al., 2005). Further elucidation of the specific interactions between gastropods and plants may therefore provide novel insights to the general regulation and diversity of plant-herbivore interactions.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
OC, ND, AS, HH, and JP conceived and designed the experiment and took the lead in writing the manuscript. All other coauthors contributed to the manuscript by giving comprehensive feedback. OC performed the experiments and analyzed the data. TL performed and assisted in transcriptomic and phytohormone analyses. AW and YP performed the metabolomics analyses and supported the metabolomic data analyses.