Transcriptional Variation in Glucosinolate Biosynthetic Genes and Inducible Responses to Aphid Herbivory on Field-Grown Arabidopsis thaliana

Recently, increasing attempts have been made to understand how plant genes function in natura. In this context, transcriptional profiles represent plant physiological status in response to environmental stimuli. Herein, we combined high-throughput RNA-Seq with insect survey data on 19 accessions of Arabidopsis thaliana grown at a field site in Switzerland. We found that genes with the gene ontology (GO) annotations of “glucosinolate biosynthetic process” and “response to insects” were most significantly enriched, and the expression of these genes was highly variable among plant accessions. Nearly half of the total expression variation in the glucosinolate biosynthetic genes (AOPs, ESM1, ESP, and TGG1) was explained by among-accession variation. Of these genes, the expression level of AOP3 differed among Col-0 accession individuals depending on the abundance of the mustard aphid (Lipaphis erysimi). We also found that the expression of the major cis-jasmone activated gene CYP81D11 was positively correlated with the number of flea beetles (Phyllotreta striolata and Phyllotreta atra). Combined with the field RNA-Seq data, bioassays confirmed that AOP3 was up-regulated in response to attack by mustard aphids. The combined results from RNA-Seq and our ecological survey illustrate the feasibility of using field transcriptomics to detect an inducible defense, providing a first step towards an in natura understanding of biotic interactions involving phenotypic plasticity.

Different defense responses of a host plant species are elicited depending on the feeding habits and host specializations of insect herbivores attacking the plant. For example, leafchewing herbivores crush plant tissues, which activates the GSL-myrosinase system (Barth and Jander, 2006;Martinoia et al., 2018;Shirakawa and Hara-Nishimura, 2018). Although the hydrolysis products of GSLs are toxic to generalist chewers (Lambrix et al., 2001;Kliebenstein et al., 2002;Barth and Jander, 2006), specialist herbivores exploit GSL and its hydrolysis products as a host plant signal (Ratzka et al., 2002;Renwick et al., 2006). Sapsuckers, such as aphids and thrips, consume plant fluids and very rarely crush plant tissues (Mewis et al., 2005;Kempema et al., 2007). In addition, damaged plants may emit volatile chemicals that elicit defenses of other individual plants or alter feeding behaviors of other insect species (Bruce et al., 2008;Matthes et al., 2011;Yazaki et al., 2017). However, only a few field studies have conducted a genome-wide analysis to address which inducible defenses may accurately function under simultaneous attacks by insect herbivores with different feeding modes .
The purpose of this study was to first reveal to what extent variations in gene expression could be explained by plant genotypes under field conditions and then to identify which herbivores could modulate plant defense responses. To address these issues, we combined our previously established protocol of cost-effective RNA-Seq (Nagano et al., 2015;Kamitani et al., 2016;Ishikawa et al., 2017) with insect monitoring data on 19 accessions of A. thaliana individuals (Table 1; Figure 1). This joint approach using transcriptome analysis and insect surveys provides an overall picture of how A. thaliana responds to multiple attackers under field conditions.

Field Experiment
In a field experiment, we used 17 natural accessions and two glabrous mutants (Table 1) of A. thaliana, thus covering a range of phenotypic variation in both trichomes (Larkin et al., 1999;Atwell et al., 2010;Bloomer et al., 2012) and glucosinolates (Chan et al., 2010). We initially prepared 10 replicates of the 19 accessions (190 plants in total) in an environmental chamber, and then we transferred them to the outdoor garden of the University of Zurich at the Irchel campus (Zurich, Switzerland: 47°23′N, 8°33′E, altitude ca. 500 m) ( Figure 1A). Plants were cultivated using mixed soils of agricultural composts (Profi Substrat Classic CL ED73, Einheitserde Co.) and perlite, with a compost to perlite ratio of 3:1 by volume. No additional fertilizers were supplied because the agricultural soils contained fertilizers. Seeds were sown on the soil and stratified under constant dark conditions at an ambient temperature of 4°C for 1 week. Plants were grown under short-day conditions (8:16-h light:dark (L:D) at 20°C and a relative humidity of 60%) for 1 month. The tray positions were rotated every week to minimize growth bias due to light conditions. Individual plants were moved to plastic pots (6.0 × 6.0 × 6.0 cm) and acclimated for 3 days in a shaded area outdoors prior to field experiments. The potted plants were randomly placed in a checkered pattern between three blocks, each containing 68, 69, or 53 plants. The potted plants were set on water-permeable plastic sheets without being embedded in the ground ( Figure 1A). Blocks were more than 1.0 m apart from each other, and the plants were watered every morning and dawn. These experiments were conducted from July 13 to August 3, 2016. Several accessions of field-grown A. thaliana were small, and it was difficult to obtain sufficient amounts of leaf tissues for metabolome analyses. Because of this limitation, we used transcriptomes as proxy to plant defense responses.
Insects and the leaf holes made by flea beetles were counted every 2-3 days on individual plants, and the final observation data were used for statistical analyses. The initial plant size (evaluated by the length in mm of the largest leaf at the start of field experiment) and presence or absence of flowering stems (2 weeks after the start of experiment) were also recorded so that these phenotypes could be included as covariates in statistical analyses. All monitoring was conducted by a single observer during the daytime (08:00-17:00) for 3 weeks after the beginning of the field experiment to minimize variation. Details of insect abundance and diversity are reported in our previous publication (Sato et al., 2019a). To avoid unintentional activation of plant defenses by mechanical damage, we did not sample any leaves until the end of the field experiment.

RNA-Seq Experiments and Data Filtering
Leaves were collected from field-grown plants at the end of the experiment (August 4, 2016). The leaf samples were immediately soaked in an RNA preservation buffer of pH 5.2 consisting of 5.3-M (NH 4 ) 2 SO 4 , 20-mM EDTA, and 25-mM trisodium citrate dihydrate at 4°C overnight and stored at −80°C until RNA extraction. Total RNA was extracted using the Maxwell 16 Lev Plant RNA Kit (Promega Japan, Tokyo) according to the manufacturer's protocol. Selective depletion of rRNAs and highly abundant transcripts was conducted prior to RNA-Seq library preparation as previously described (Nagano et al., 2015). Then, RNA-Seq libraries were prepared as previously described (Ishikawa et al., 2017). Sequencing using Illumina HiSeq ® 2500 was carried out by Macrogen Co. We sequenced 92 samples per lane and obtained 829,681 mapped reads per sample on average.
The fastq files generated by sequencing were preprocessed using Trimmomatic version 0.32 (Bolger et al., 2014). The preprocessed sequences were mapped onto the A. thaliana reference genome (TAIR10 cDNA) using Bowtie version 1.1.1 (Langmead et al., 2009) and then quantified using RSEM version 1.2.21 (Li and Dewey, 2011). The parameter settings of Trimmomatic, Bowtie, and RSEM were the same as those described by  Atwell et al., 2010), and GSL accumulation (Chan et al., 2010). *Obtained from the Kiyotaka Okada Laboratory of Kyoto University, Japan. † Obtained from Dr. M. Ohto. ‡ Estimated from the trichome density relative to the Col-0 accession presented in previous publications (Hauser et al., 2001 andYoshida et al., 2009 for Ms-0 andCol(gl1-2), respectively). $ Short-chain: the sum of C3-and C4-aliphatic GSLs; Long-chain: the sum of C7-and C8-aliphatic GSLs obtained from Chan et al. (2010). (C) Filtering and statistical analysis of RNA-Seq data. In the ANOVA formula, "Herbivore" represents the main effect of the number of herbivores, while the "A × H" indicates the interaction term between the plant accession and the number of herbivores. Kamitani et al. (2016). Following Kamitani et al. (2016), we calculated the raw read counts and reads per million (rpm) from the expected read counts generated with RSEM. Transposable elements were excluded prior to statistical analyses. We calculated the total raw read counts for each plant sample and discarded shallow-read samples belonging to the lower fifth percentile of the total raw read counts ( Figure 1C). Consequently, samples with >12,130 reads were subjected to statistical analyses. To exclude non-expressed genes, we then averaged log 2 (rpm + 1) for each gene between all plant samples and eliminated genes with an average log 2 (rpm + 1) of zero ( Figure 1C). Overall, we obtained a final dataset on 24,539 genes for 173 plants. In this final dataset, 53 out of 173 samples had <10 5 total reads. Overall trends did not change when we set the threshold at 10 5 , although the statistical power decreased due to lower sample size. Sequence data from our RNA-Seq were submitted to the NCBI Sequence Read Archive repository under the BioProject number, PRJNA488315 (https://www.ncbi.nlm.nih.gov/bioproject/ PRJNA488315). Read count data and source code are available via the GitHub repository (https://github.com/naganolab/ AthRNAseq2016Zurich_Sato_et_al).

Statistical Analysis
We used Type III analysis of variance (ANOVA; Sokal and Rohlf, 2012) to screen genes showing large expression variation among accessions ( Figure 1C). We formulated the linear model as: Y ~ Accession ID (factorial) + Flowering stem (binary data) + Initial leaf length (mm), where Y indicates log 2 (rpm + 1) of a focal gene. Sum of squares (SS) were calculated to partition expression variation attributable to each explanatory variable. The proportion of expression variation explained by the plant accession was evaluated as the SS of the plant accession ID divided by the total SS. Genes in the top 5% of expression variation were selected and subjected to statistical analysis, as described below. All statistical analyses were performed using R version 3.2.0 (R Core Team, 2015).
Subsequently, gene ontology (GO) enrichment analysis was applied to the genes that showed the top 5% of values of the proportion of expression variation explained by the plant accession. The "GO.db" package (Carlson, 2017) and "TAIR10" gene annotation package were used to obtain the GO terms. The statistical significance of each GO term was determined using Fisher's exact tests against the entire database. The P-values were adjusted by the false discovery rate (FDR; Benjamini and Hochberg, 1995) using the "p.adjust" function of R. When significant GO terms were detected, the "GOBPOFFSPRING" database in the "GO.db" package was used to find the most specific GO within the biological process.
We next incorporated the effects of herbivores into an ANOVA to determine whether insect herbivory altered gene expression among plant accessions. We formulated linear models as: Y ~ Accession ID (factorial) + No. of herbivores + Flowering stem (binary data) + Initial leaf length (mm) + (Accession ID × No. of herbivores). This ANOVA was repeated for five different explanatory variables for the No. of herbivores term: the number of mustard aphids (L. erysimi), flea beetles (P. striolata and P. atra), leaf holes made by P. striolata and P. atra, diamondback moths (P. xylostella), and western flower thrips (F. occidentalis). Significance of these terms in ANOVAs tested whether the expression of a focal gene might be induced by the herbivore species. The interaction term (Accession ID × No. of herbivores; A × H in Figure 1C) represented the combined effect exerted by the plant accession and the number of herbivores, and it tested whether the presence or magnitude of inducible defenses to the herbivore species depended on plant accessions having different genomic backgrounds. The number of herbivores was log-transformed to improve normality. Given that previous laboratory experiments detected inducible defenses 24-48 h after insect attacks (e.g., Mewis et al., 2005;Kuśnierczyk et al., 2008;Matthes et al., 2011), we used insect abundance data on August 3, 2016, that is, 1 day prior to RNA sampling, for explanatory variables. The P-values were calculated using F-tests corrected by FDR. The "aov" function in R was used to perform ANOVAs, and we compared models with or without a focal term by the F-tests.

Laboratory Bioassay and RT-qPCR
To determine whether AOP3 was induced after attack by the mustard aphid, L. erysimi, we released this aphid species onto plants of the Col-0 accession under controlled conditions in the laboratory, and then we quantified expression of AOP3 in infested and non-infested plants (Figure 4A, B). Lipaphis erysimi were collected from Rorippa indica growing at the Seta campus of Ryukoku University, Japan (34°58′N, 135°56′E), and maintained on leaves of Raphanus sativus "Longipinnatus" before the bioassay. Seeds were sown in plastic pots (6 cm in diameter and height) filled with moist vermiculite. After germination, seedlings were thinned to four per pot. Seedlings were grown under 16L:8D conditions at an ambient temperature of 20°C for 1 month. Liquid fertilizer (diluted 2,000×) was supplied to plants during their cultivation (Hyponex, Hyponex Japan, Osaka; N:P:K = 6:10:5). We assigned eight plants (two pots of four plants each) to the aphid treatment and eight plants (two pots of four plants each) for controls. Approximately 80 wingless aphids were released per pot for the aphid treatment, and the pots were separately covered with mesh. Leaf sampling was conducted once a week after the release of aphids. Leaves were soaked in an RNA preservation buffer of pH 5.2 consisting of 5.3-M (NH 4 ) 2 SO 4 , 20-mM EDTA, and 25-mM trisodium citrate dihydrate overnight and stored at −80°C until RNA extraction.
Total RNA was extracted using a Maxwell 16 Lev Plant RNA Kit (Promega Japan, Tokyo), and RNA concentration was measured using a Quant-iT RNA Assay Kit Broad Range (Invitrogen, Carlsbad, CA, USA). cDNA was synthesized from 300 ng of the total RNA using a reaction solution composed of 10 μL of template RNA, 4.0 μL of 5× SuperScript IV Reverse Transcriptase Buffer (Invitrogen, Carlsbad, CA), 0.5 μL of RNasin ® Plus RNase inhibitor (Promega Japan, Tokyo), 2.0 μL of 100-mM DTT (Invitrogen, Carlsbad, CA), 0.4 μL of 25-mM dNTPs (Clontech, Palo Alto, CA), 0.5 μL of SuperScript IV Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA), and 0.6 μL of 100-μM random primer (N)6 (TaKaRa, Kusatsu, Japan), with RNase-free water added up to 20 μL. For the reverse transcription step, the mixture was incubated at 25°C for 10 min, followed by 50 min at 56°C. SuperScript IV was inactivated by heating the mixture at 75°C for 15 min. The cDNA was 10× diluted with RNase-free water, and then reverse transcription-quantitative polymerase chain reaction (RT-qPCR) was performed using a Roche LightCycler ® 480 with 10 μL of reaction solution composed of 2 μL of the template, 5 μL of KAPA SYBR Fast RT-qPCR solution (Kapa Biosystems, Inc., Woburn, MA), 0.5 μL of 10-μM forward and reverse primers, and 2 μL of RNase-free water. The forward and reverse primers of the target gene AOP3 was 5′-TCAGGGGTCGGTTTTGAAGG-3′ and 5′-GTGAAAGGTTTCGGGCACAC-3′, respectively. Expressions of ACT2 and EF-1α were measured for the internal controls (see Czechowski et al., 2005 for primer sequences). The cDNAs were amplified following denaturation, with 35 cycles of 10 s at 95°C, 20 s at 63°C, and 10 s at 72°C. Three technical replicates were used for individual plants and primers. The primer of the target gene AOP3 was designed based on its full-length coding sequence using NCBI Primer-BLAST with a product length parameter of 50-150 bp. We tried six candidate primers and selected the AOP3 primer above based on the melting curve of laboratory-grown Col-0 and Ler-1 accessions. Cp values were calculated following the second derivative maximum method and averaged among three technical replicates. The geometric mean of Cp values between ACT2 and EF-1α was used as the internal control. Delta Cp values were calculated for each individual plant between the target and internal control. A Wilcoxon rank sum test was used to test differences in the delta Cp values between the intact and aphid-infested plants.

A. thaliana Accessions
Our statistical analysis showed that, when ordered by the expression variation explained by plant accessions, the top 5% of genes had more than 20% of their variation attributable to the plant accessions (Figure 2A; Table S1). In these highly variable genes, 22 GOs were significantly enriched at P FDR < 0.05 ( Table 2). The GO annotations of "response to insect" and "glucosinolate biosynthetic process" were most and second most significantly enriched, respectively ( Table 2).
Several key genes of aliphatic GSL biosynthesis, such as AOPs (Kliebenstein et al., 2001a), had over half of their expression variation attributable to plant accession (Figure 2A, C, D). AOP2 and AOP3 are located nearby in the genome and encode 2-oxoglutarate-dependent dioxygenases involved in side chain modification of aliphatic GSLs (Kliebenstein et al., 2001a). If AOP3 is expressed, plants accumulate hydroxypropyl. If AOP2 is expressed, plants accumulate alkenyl (Kliebenstein et al., 2001a;Chan et al., 2010). Mean expression levels of AOP3 or AOP2 among 17 A. thaliana accessions transplanted in the field were highly correlated with the known levels of hydroxypropyl or alkenyl accumulated in plants grown in a growth camber (Figure S1). Genes involved in GSL hydrolysis, such as TGGs, ESM1, and ESP (Lambrix et al., 2001;Barth and Jander, 2006;Zhang et al., 2006), exhibited a similarly large variation in expression; these genes had nearly half of their variation attributable to the plant accession (Figures 2B, E,  F, G). TGG1 and TGG2 are functionally redundant, and their double mutants are known to become palatable for generalist caterpillars, but not to aphids and specialist caterpillars (Barth and Jander, 2006). ESM1 and ESP were initially screened by quantitative trait locus (QTL) mapping utilizing the natural variation between the Col and Ler accession (Jander et al., 2001;Lambrix et al., 2001;Zhang et al., 2006). Furthermore, the transcription factor gene MYB29, which is responsible for the high accumulation of aliphatic GSLs (Hirai et al., 2007), showed 32% variation in the expression level among fieldgrown accessions (Figure 2H).

Inducible Responses to Leaf-Chewing and Sap-Sucking Herbivores
The results of ANOVA indicated that 27, 25, 20, and 19 candidate genes were significantly related to inducible responses to mustard aphids, flea beetles, diamondback moths, and western flower thrips, respectively (Table S2).
In response to mustard aphid feeding (L. erysimi), AOP3 had a significant interaction between its expression with the number of aphids (P FDR < 0.001: Figure 3A, Table S2), indicating that its induction by aphids depended on background genomic variation. This statistical interaction explained 7% of variation in AOP3 expression. Similar interactions with aphid herbivory were observed for MYB113 and JAX1. MYB113 is known to be involved in anthocyanin biosynthesis and induced via JA signaling (Gonzalez et al., 2008), and JAX1 encodes jacalin-type lectin resistance to potexvirus and exhibits varying resistance among natural accessions (Yamaji et al., 2012).
The number of flea beetles (P. striolata and P. atra) was positively correlated with the expression level of CYP81D11 (P FDR = 0.027: Figure 3B, Table S2), which is known to be a major cis-jasmone activated gene (Matthes et al., 2010;Matthes et al., 2011). However, its expression was not affected by plant accession (plant accession × beetles, P FDR = 0.99), indicating limited effects of background genomic variation. Additionally, the number of leaf holes made by the flea beetles was only related to the expression of three loci, all of unknown function, AT2G41590 (plant accession × holes, P FDR < 10 −6 ), AT1G34844 (plant accession × holes, P FDR < 10 −13 ), and AT2G47570 (plant accession × holes, P FDR = 0.007).
In addition, AT5G48770, which encodes disease resistance proteins of the TIR-NBS-LRR family and has a GO annotation of "defense response, " was expressed in response to the diamondback moth (P. xylostella) ( Table S2).
Finally, the presence of the western flower thrips (F. occidentalis) was correlated with the expression of one locus, AT2G15130. This locus encodes a basic secretory protein family protein in plants and has the GO annotation "defense response" (Table S2).  (Table S1). September 2019 | Volume 10 | Article 787 Frontiers in Genetics | www.frontiersin.org

Laboratory Bioassay Using the Specialist Aphids
The statistical interaction between plant accession and number of aphids (Section 3.3) could be explained by the positive correlation we found between the expression level of AOP3 and the number of mustard aphids (L. erysimi) in Col(gl1-2) plants ( Figure 3A). In our laboratory bioassay, the expression of AOP3 was up-regulated in aphid-infested Col-0 plants (Wilcoxon rank sum test, W = 0, n = 16, P = 0.0002: Figure 4C).

Expression Variation in Glucosinolate Biosynthetic Genes
Glucosinolate profiles, leaf damage, and plant fitness significantly vary among A. thaliana accessions, as demonstrated by similar field experiments using GSL mutants (Kerwin et al., 2015;Kerwin et al., 2017) or a number of natural accessions (Brachi et al., 2015). Previous laboratory studies show that GSL biosynthesis is affected by many factors other than herbivory, including nutrient deficiency (Hirai et al., 2007), light conditions (Wu et al., 2018), drought (Mewis et al., 2012), and pathogen infection (Kliebenstein et al., 2005;Gudiño et al., 2018;Shirakawa and Hara-Nishimura, 2018). In our field study, where these factors were not controlled for and subject to natural variation, GO enrichment of "GSL biosynthesis" and "response to insects" genes corresponded to the top 5% most variably expressed genes among A. thaliana accessions ( Table 2). These two GOs were the most and second most significantly enriched terms, while "systemic acquired resistance" and "response to water deprivation" were also noted. This result suggests that anti-herbivore defense may be one of primary functions of GSL biosynthesis in the field. Notably, nearly half of the expression variation in AOPs, ESM1, and TGG1 was explained by plant accessions (Figure 2), showing a comparable magnitude of variation with the heritability reported by a laboratory eQTL study (Wentzell et al., 2007). A third of the expression variation in the transcription factor gene MYB29 was attributable to plant accessions, even though this gene is known to respond to water stress and other abiotic stimuli   -Ballesta et al., 2015;Zhang et al., 2017). Overall, our genome-wide analysis using RNA-Seq indicated that GSL biosynthesis and anti-herbivore functions were among the most genetically variable functions in field-grown A. thaliana. Among the GSL biosynthetic genes, AOPs showed remarkably high variation in expression among natural accessions. More specifically, the Ler-1 accession expressed AOP3 and not AOP2, whereas Cvi expressed AOP2 and not AOP3 (Figures 2C, D). These findings are similar to those of Kliebenstein et al. (2001a). In the Col accession, AOP2 encodes non-functional proteins (Kliebenstein et al., 2001a), and AOP3 is not expressed in intact leaves (Schmid et al., 2005). Strong genome-wide associations between the AOP loci and GSL profiles have been repeatedly detected among natural accessions cultivated under laboratory (Chan et al., 2010) or controlled greenhouse (Brachi et al., 2015) conditions. Based on a genome scan, Brachi et al. (2015) also detected an adaptive differentiation in the AOP loci within European A. thaliana. In an evolutionary context, the present study provided observational evidence for a link between genomic and functional variations in AOPs in the field.

Genes Possessing Inducible Responses to Herbivory
Consistent with the field RNA-Seq data, our laboratory bioassay revealed that the Col-0 accession had an induced response of AOP3 to mustard aphids (L. erysimi). Besides L. erysimi, the generalist aphid Myzus persicae and the specialist aphid Brevicoryne brassicae are major natural enemies of A. thaliana (Züst et al., 2012). In an ecological context, previous studies reported unclear geographical associations between these three aphid species and AOP-related chemotypes in Europe, though the geographical distribution of the aphids was linked to that of MAM-related chemotypes (Züst et al., 2012). Our previous study also reported no significant correlations between laboratory-measured profiles of aliphatic GSL and the abundance of L. erysimi in the field (Sato et al., 2019a). Results of this study, and from microarray analyses (Kempema et al., 2007;Kuśnierczyk et al., 2008), showed that the aphid species differentially induced AOP3. This gene may be up-regulated in Col-0 by L. erysimi, down-regulated in Ler by B. brassicae (Kuśnierczyk et al., 2008), and not induced in Col by M. persicae (Kempema et al., 2007). Given this variation in response to different aphid species, the inducible response of AOPs help to explain why AOP loci are not tightly linked to higher phenotypes such as aphid resistance in wild populations.
Of the genes with the GO annotation of "response to insects, " CYP81D11 exhibited a significant positive correlation between its expression and the abundance of flea beetles. CYP81D11 is known to be up-regulated by cis-jasmone, a plant volatile emitted via wounds from insect attack or pathogen infection (Bruce et al., 2008;Matthes et al., 2010;Matthes et al., 2011). Matthes et al. (2011) used Col background A. thaliana as their standard accession, while in the present study, we included multiple natural accessions and found a positive correlation between flea beetle abundance and CYP81D11 expression ( Figure 3B). However, there was no significant correlation between CYP81D11 expression and the number of leaf holes made by these beetles. This result was probably because the leaf holes remained on the leaves for a few weeks and accumulated without reflecting the timing of wounding. Our results on CYP81D11 indicated that data on both herbivory and insect abundance were needed to detect the induced response to flea beetles, exemplifying the importance of detailed ecological observations during in natura studies.

CONCLUSION
A combination of insect surveys with field transcriptome analyses allowed us to detect an inducible defense against insect herbivores in A. thaliana. These results suggest that the molecular machinery of Arabidopsis defense can function accurately in complex environments. Whereas previous field studies on a Brassicaceae crops reported significant transcriptional changes in response to entire herbivore communities , our large-scale RNA-Seq and insect monitoring data allowed for assessment of transcriptional responses specific to individual herbivore species. Since the insect species studied here are also known as herbivores of cultivated and wild Brassicaceae worldwide (Yano, 1994;Ahuja et al., 2010;Sato and Kudoh, 2017;Sato, 2018), our findings may provide general molecular insights into Brassicaceae-herbivore interactions in natura. Further studies are needed to reveal how the inducible responses at the transcription level modulate defense metabolism and insect resistance under complex field conditions.

DATA AVAILABILITY
The datasets generated for this study can be found in NCBI Sequence Read Archive repository, under the BioProject number, PRJNA488315.

AUTHOR CONTRIBUTIONS
YS conducted plant sampling, insect monitoring, and statistical analyses. YS and AT conducted the bioassay and RT-qPCR analysis. AT, MK, and AD performed RNA-Seq experiments. YS, RS-I, and MY designed the field experiments. YS, KS, and AN conceived the study and wrote the paper with inputs from all co-authors.