Yasuhiro Sato1,2
Ayumi Tezuka2Makoto Kashima2Ayumi Deguchi2,3
Rie Shimizu-Inatsugi4Misako Yamazaki4
Kentaro K. Shimizu4,5*
Atsushi J. Nagano6*- 1PRESTO, Japan Science and Technology Agency, Kawaguchi, Japan
- 2Research Institute for Food and Agriculture, Ryukoku University, Otsu, Japan
- 3Graduate School of Horticulture, Chiba University, Matsudo, Japan
- 4Department of Evolutionary Biology and Environmental Studies, University of Zurich, Zurich, Switzerland
- 5Kihara Institute for Biological Research, Yokohama City University, Yokohama, Japan
- 6Department of Plant Life Sciences, Faculty of Agriculture, Ryukoku University, Otsu, Japan
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.
Introduction
As sessile organisms, plants are exposed to multiple stresses under naturally fluctuating environments (Wilczek et al., 2009; Carrera et al., 2017; Mishra et al., 2017). Recently, increased efforts have been made to understand how plants cope with complex field conditions (Kerwin et al., 2015; Carrera et al., 2017; Kono et al., 2017; Sugiyama et al., 2017; Taylor et al., 2017; Hiraki et al., 2018). Such in natura studies are important to gain a comprehensive understanding of gene functions from the laboratory to the field (Shimizu et al., 2011; Kudoh, 2016; Yamasaki et al., 2018; Zaidem et al., 2019; Nagano et al., 2019). Insect herbivores are the most diverse group of organisms that impose biotic stresses on land plants (Schoonhoven et al., 2005; Ahuja et al., 2010; Escobar-Bravo et al., 2018). To deal with various threats, plants activate some defense mechanisms only when necessary. Such inducible defenses are triggered through jasmonate (JA) signaling after wounding or insect attacks (Mewis et al., 2005; Escobar-Bravo et al., 2017; Tsuda, 2017; Nakano and Mukaihara, 2018; Zhou et al., 2018; Zhu et al., 2018), whereas constitutive defenses are continuously expressed. Notably, the magnitude of inducible defense varies among plant genotypes (Agrawal et al., 2002; Kuśnierczyk et al., 2008; Snoeren et al., 2010). Furthermore, spatiotemporal variation in herbivory and insect abundance in the field could modulate defense metabolism through phenotypic plasticity or gene-by-environment interactions (Agrawal, 2001; Kerwin et al., 2015; Kerwin et al., 2017).
In the glucosinolate (GSL)–myrosinase system of Arabidopsis thaliana and related Brassicales, methionine-derived or aliphatic GSLs confer plant defenses against herbivory (Kliebenstein et al., 2002; Brachi et al., 2015; Kerwin et al., 2015). The accumulation and profiles of these GSLs are variable among A. thaliana accessions worldwide (Kliebenstein et al., 2001b; Kroymann et al., 2003; Chan et al., 2010; Brachi et al., 2015). Production of aliphatic GSLs is initiated by MYB28 and MYB29 transcription factors (Hirai et al., 2007). Double mutants of both of these genes accumulate few aliphatic GSLs (Sønderby et al., 2007). During the accumulation of aliphatic GSLs, amino acids and side chain structures are modified by methylthioalkylmalate synthase (MAM), 2-oxoglutarate-dependent dioxygenase (encoded in alkenyl hydroxalkyl producing (AOP) loci), and flavin-monooxygenase glucosinolate S-oxygenase (Kliebenstein et al., 2001a; Kroymann et al., 2003; Hansen et al., 2007). When insect herbivores bite plant tissues, the enzyme myrosinase [also known as thioglucoside glucohydrolase (TGG)] catalyzes breakdown of GSLs, resulting in emission of isothiocyanates, nitriles, or other hydrolysis products (Lambrix et al., 2001; Barth and Jander, 2006; Zhang et al., 2006; Shirakawa and Hara-Nishimura, 2018). Epithiospecifier proteins (ESPs, also known as TASTY; Jander et al., 2001; Lambrix et al., 2001) promote the hydrolysis of GSL with some modification by the EPITHIOSPECIFIER MODIFIER1 (ESM1) locus (Zhang et al., 2006), resulting in different defense activities against insect herbivores (Ratzka et al., 2002).
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, leaf-chewing 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 (Broekgaarden et al., 2010).
In wild populations, A. thaliana has multiple generations per year (Thompson, 1994; Wilczek et al., 2009; Taylor et al., 2017) and is attacked by various herbivores (Arany et al., 2005; Harvey et al., 2007; Sato et al., 2019a). Previously, Sato et al. (2019a) found that 12 insect species, including the mustard aphid (Lipaphis erysimi; Homoptera), flea beetles (Phyllotretastriolata and Phyllotreta atra; Coleoptera), the diamondback moth (Plutella xylostella; Lepidoptera), and the western flower thrips (Frankliniella occidentalis; Thysanoptera), colonized an experimental summer population of A. thaliana near Zurich, Switzerland. In addition, the insect community composition significantly varied among A. thaliana accessions (Sato et al., 2019a); however, there was no correlation between herbivore abundance in the field and GSL profiles of Arabidopsis accessions grown in growth chambers (Sato et al., 2019a). To develop a more comprehensive picture of plant defense expression in the field, transcriptomics techniques have been applied to assess the environmental responses of many plant species (e.g., Kempema et al., 2007; Broekgaarden et al., 2010; Kamitani et al., 2016; Sun et al., 2017; Lin et al., 2017; Wang et al., 2018; Xu et al., 2018).
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.
Table 1 List of Arabidopsis thaliana accessions used in this study.
Figure 1 Outline of the field study on A. thaliana. (A) Procedure of the field experiment. (B) Observed variation in insect abundance among plant accessions. (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.
Materials and Methods
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 (NH4)2SO4, 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 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 log2(rpm + 1) for each gene between all plant samples and eliminated genes with an average log2(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 <105 total reads. Overall trends did not change when we set the threshold at 105, 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 log2(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 (NH4)2SO4, 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.
Results
Insect Herbivores Observed on Field-Grown A. thaliana
The major insect herbivores observed included L. erysimi, P.striolata, P. atra, P. xylostella, and F. occidentalis (Sato et al., 2019a). Of these herbivores, L. erysimi and F. occidentalis are sucking insects, while P. xylostella, P. striolata, and P. atra are leaf chewers (Ahuja et al., 2010; Escobar-Bravo et al., 2018). L. erysimi, P. xylostella, P. striolata, and P. atra are specialists of Brassicaceae (Ahuja et al., 2010), whereas F. occidentalis is a generalist that feeds on various plant families (Escobar-Bravo et al., 2018) (Figure 1B).
Gene Expression Variation Among 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 PFDR < 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).
Figure 2 Natural variation in the expression levels of genes involved in glucosinolate biosynthesis and hydrolysis. (A) Histogram showing the proportion of variation explained by plant accessions, (B) the expression of ESM1, (C)AOP3, (D)AOP2, (E)TGG1, (F)TGG2, (G)ESP, and (H)MYB29. Grey bars and vertical lines indicate mean ± SE. The list of the top 5% most variable genes is available in the Supporting Information (Table S1).
Table 2 Gene ontology (GO) enrichment analysis of genes with the top 5% of expression variation explained by plant accessions. P-values are corrected by the false discovery rate (PFDR). Shown are the significant GOs within the biological process terms (PFDR < 0.05).
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 field-grown 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 (PFDR < 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 (PFDR = 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, PFDR = 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, PFDR < 10−6), AT1G34844 (plant accession × holes, PFDR < 10−13), and AT2G47570 (plant accession × holes, PFDR = 0.007).
Figure 3 Correlations between candidate gene expression and insect abundance in the field. (A) Relationship between the number of specialist mustard aphids (L. erysimi) and the expression of AOP3 in Ler-1 (an accession with constitutive expression) and Col(gl1-2) (an accession with potentially induced expression). (B) Relationship between the number of Phyllotreta beetles and expression of CYP81D11; residuals unexplained by plant accessions are shown in the Y-axis.
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).
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).
Figure 4 Inducible response of AOP3 to the specialist mustard aphid L. erysimi in the laboratory-grown A. thaliana Col-0 accession. (A) Photograph of a laboratory-reared colony of L. erysimi. (B) Aphid-infested seedlings on Col-0 accession plants. (C) Relative expression of AOP3 in RT-qPCR analysis of aphid-infested and control plants; seedlings were infested by aphids for 7 days (Aphid) or not infested (Cont.).
Discussion
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 (Martinez-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 (Broekgaarden et al., 2010), 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.
Funding
This study was funded by the Japan Society for the Promotion of Science (JSPS) Postdoctoral Fellowship (grant number, 16J30005) and Japan Science and Technology Agency (JST) PRESTO (JPMJPR17Q4) to YS, JST CREST (JPMJCR15O2) to AN, and JST CREST (JPMJCR16O3), MEXT KAKENHI (18H04785), and the Swiss National Science Foundation to KS. The field experiment was supported by the University Research Priority Program of Global Change and Biodiversity at the University of Zurich.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
The authors thank F. Kobayashi for assistance in RNA extraction and Dynacom Co. Ltd. for help with the RNA-Seq data analyses. This manuscript has been released as a bioRxiv Pre-Print at https://www.biorxiv.org/content/10.1101/563486v1 (Sato et al., 2019b).
Abbreviations
ANOVA, analysis of variance; AOP, alkenyl hydroxalkyl producing; ESM1, epithiospecifier modifier 1; ESP, epithiospecifier protein; FDR, false discovery rate; GO, gene ontology; GSL, glucosinolate; MAM, methylthioalkylmalate synthase; rpm, reads per million; TGG, thioglucoside glucohydrolase.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00787/full#supplementary-material
References
Agrawal, A. A. (2001). Phenotypic plasticity in the interactions and evolution of species. Science 294, 321–326. doi: 10.1126/science.1060701
Agrawal, A. A., Conner, J. K., Johnson, M. T., Wallsgrove, R., Poulin, R. (2002). Ecological genetics of an induced plant defense against herbivores: additive genetic variance and costs of phenotypic plasticity. Evolution 56, 2206–2213. doi: 10.1554/0014-3820(2002)056[2206:EGOAIP]2.0.CO;2
Ahuja, I., Rohloff, J., Bones, A. M. (2010). Defence mechanisms of Brassicaceae: implications for plant-insect interactions and potential for integrated pest management: a review. Agron. Sustain. Dev. 30, 311–348. doi: 10.1051/agro/2009025
Arany, A. M., de Jong, T. J., van der Meijden, E. (2005). Herbivory and abiotic factors affect population dynamics of Arabidopsis thaliana in a sand dune area. Plant Biol. 7, 549–556. doi: 10.1055/s-2005-865831
Atwell, S., Huang, Y. S., Vilhjálmsson, B. J., Willems, G., Horton, M., Li, Y., et al. (2010). Genome-wide association study of 107 phenotypes in Arabidopsis thaliana inbred lines. Nature 465, 627–631. doi: 10.1038/nature08800
Barth, C., Jander, G. (2006). Arabidopsis myrosinases TGG1. TGG2 have redundant function in glucosinolate breakdown and insect defense. Plant J. 46, 549–562. doi: 10.1111/j.1365-313X.2006.02716.x
Benjamini, Y., Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289–300. doi: 10.2307/2346101
Bloomer, R. H., Juenger, T. E., Symonds, V. V. (2012). Natural variation in GL1 and its effects on trichome density in Arabidopsis thaliana. Mol. Ecol. 21, 3501–3515. doi: 10.1111/j.1365-294X.2012.05630.x
Bolger, A. M., Lohse, M., Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Brachi, B., Meyer, C. G., Villoutreix, R., Platt, A., Morton, T. C., Roux, F., et al. (2015). Coselected genes determine adaptive variation in herbivore resistance throughout the native range of Arabidopsis thaliana. Proc. Natl. Acad. Sci. U.S.A. 112, 4032–4037. doi: 10.1073/pnas.1421416112
Broekgaarden, C., Poelman, E. H., Voorrips, R. E., Dicke, M., Vosman, B. (2010). Intraspecific variation in herbivore community composition and transcriptional profiles in field-grown Brassica oleracea cultivars. J. Exp. Bot. 61, 807–819. doi: 10.1093/jxb/erp347
Bruce, T. J., Matthes, M. C., Chamberlain, K., Woodcock, C. M., Mohib, A., Webster, B., et al. (2008). cis. Arabidopsis genes that affect the chemical ecology of multitrophic interactions with aphids and their parasitoids. Proc. Natl. Acad. Sci. U.S.A. 105, 4553–4558. doi: 10.1073/pnas.0710305105
Carlson, M. (2017). GO.db: a set of annotation maps describing the entire gene ontology. R package version 3.5.0. doi: 10.18129/B9.bioc.GO.db
Carrera, D.Á., Oddsson, S., Grossmann, J., Trachsel, C., Streb, S. (2017). Comparative proteomic analysis of plant acclimation to six different long-term environmental changes. Plant Cell Physiol. 59, 510–526. doi: 10.1371/journal.pone.0075705
Chan, E. K. F., Rowe, H. C., Kliebenstein, D. J. (2010). Understanding the evolution of defense metabolites in Arabidopsis thaliana using genome-wide association mapping. Genetics 185, 991–1007. doi: 10.1534/genetics.109.108522
Czechowski, T., Stitt, M., Altmann, T., Udvardi, M. K., Scheible, W. (2005). Genome-wide identification and testing of superior reference genes for transcript normalization in Arabidopsis. Plant Physiol. 139, 5–17. doi: 10.1104/pp.105.063743
Escobar-Bravo, R., Klinkhamer, P. G., Leiss, K. A. (2017). Induction of jasmonic acid-associated defenses by thrips alters host suitability for conspecifics and correlates with increased trichome densities in tomato. Plant Cell Physiol. 58, 622–634. doi: 10.1093/pcp/pcx014
Escobar-Bravo, R., Ruijgrok, J., Kim, H. K., Grosser, K., Van Dam, N. M., Klinkhamer, P. G., et al. (2018). Light intensity-mediated induction of trichome-associated allelochemicals increases resistance against thrips in tomato. Plant Cell Physiol. 59, 2462–2475. doi: 10.1093/pcp/pcy166
Gonzalez, A., Zhao, M., Leavitt, J. M., Lloyd, A. M. (2008). Regulation of the anthocyanin biosynthetic pathway by the TTG1/bHLH/Myb transcriptional complex in Arabidopsis seedlings. Plant J. 53, 814–827. doi: 10.1111/j.1365-313X.2007.03373.x
Gudiño, M. E., Blanco-Touriñán, N., Arbona, V., Gómez-Cadenas, A., Blázquez, M. A., Navarro-García, F. (2018). β-Lactam antibiotics modify root architecture and indole glucosinolate metabolism in Arabidopsis thaliana. Plant Cell Physiol. 59, 2086–2098. doi: 10.1093/pcp/pcy128
Harvey, J. A., Witjes, L. M. A., Benkirane, M., Duyts, H., Wagenaar, R. (2007). Nutritional suitability and ecological relevance of Arabidopsis thaliana and Brassica oleracea as foodplants for the cabbage butterfly, Pieris rapae. Plant Ecol. 189, 117–126. doi: 10.1007/s11258-006-9204-6
Hansen, B. G., Kliebenstein, D. J., Halkier, B. A. (2007). Identification of a flavin-monooxygenase as the S-oxygenating enzyme in aliphatic glucosinolate biosynthesis in Arabidopsis. Plant J. 50, 902–910. doi: 10.1111/j.1365-313X.2007.03101.x
Hauser, M.-T., Harr, B., Schlötterer, C. (2001). Trichome distribution in Arabidopsis thaliana and its close relative Arabidopsis lyrata: molecular analysis of the candidate gene GLABROUS1. Mol. Biol. Evol. 18, 1754–1763. doi: 10.1093/oxfordjournals.molbev.a003963
Hirai, M. Y., Sugiyama, K., Sawada, Y., Tohge, T., Obayashi, T., Suzuki, A., et al. (2007). Omics-based identification of Arabidopsis Myb transcription factors regulating aliphatic glucosinolate biosynthesis. Proc. Natl. Acad. Sci. U.S.A. 104, 6478–6483. doi: 10.1073/pnas.0611629104
Hiraki, H., Uemura, M., Kawamura, Y. (2018). Calcium signaling-linked CBF/DREB1 gene expression was induced depending on the temperature fluctuation in the field: views from the natural condition of cold acclimation. Plant Cell Physiol. 60, 303–317. doi: 10.1093/pcp/pcy210
Ishikawa, T., Kashima, M., Nagano, A. J., Ishikawa-Fujiwara, T., Kamei, Y., Todo, T., et al. (2017). Unfolded protein response transducer IRE1-mediated signaling independent of XBP1 mRNA splicing is not required for growth and development of medaka fish. eLife 6, e26845. doi: 10.7554/eLife.26845
Jander, G., Cui, J., Nhan, B., Pierce, N. E., Ausubel, F. M. (2001). The TASTY locus on chromosome 1 of Arabidopsis affects feeding of the insect herbivore Trichoplusia ni. Plant Physiol. 126, 890–898. doi: 10.1104/pp.126.2.890
Kamitani, M., Nagano, A. J., Honjo, M. N., Kudoh, H. (2016). RNA-Seq reveals virus–virus and virus–plant interactions in nature. FEMS Microbiol. Ecol. 92, fiw176. doi: 10.1093/femsec/fiw176
Kempema, L. A., Cui, X., Holzer, F. M., Walling, L. L. (2007). Arabidopsis transcriptome changes in response to phloem-feeding silverleaf whitefly nymphs: similarities and distinctions in responses to aphids. Plant Physiol. 143, 849–865. doi: 10.1104/pp.106.090662
Kerwin, R. E., Feusier, J., Corwin, J., Rubin, M., Lin, C., Muok, A., et al. (2015). Natural genetic variation in Arabidopsis thaliana defense metabolism genes modulates field fitness. eLife 4, e05604. doi: 10.7554/eLife.05604.001
Kerwin, R. E., Feusier, J., Muok, A., Lin, C., Larson, B., Copeland, D., et al. (2017). Epistasis × environment interactions among Arabidopsis thaliana glucosinolate genes impact complex traits and fitness in the field. New Phytol. 215, 1249–1263. doi: 10.1111/nph.14646
Kliebenstein, D. J., Lambrix, V. M., Reichelt, M., Gershenzon, J., Mitchell-Olds, T. (2001a). Gene duplication in the diversification of secondary metabolism: tandem 2-oxoglutarate-dependent dioxygenases control glucosinolate biosynthesis in Arabidopsis. Plant Cell 13, 681–693. doi: 10.1105/tpc.13.3.681
Kliebenstein, D. J., Kroymann, J., Brown, P., Figuth, A., Pedersen, D., Gershenzon, J., et al. (2001b). Genetic control of natural variation in Arabidopsis glucosinolate accumulation. Plant Physiol. 126, 811–825. doi: 10.1104/pp.126.2.811
Kliebenstein, D., Pedersen, D., Barker, B., Mitchell-Olds, T. (2002). Comparative analysis of quantitative trait loci controlling glucosinolates, myrosinase and insect resistance in Arabidopsis thaliana. Genetics 161, 325–332.
Kliebenstein, D. J., Rowe, H. C., Denby, K. J. (2005). Secondary metabolites influence Arabidopsis/Botrytis interactions: variation in host production and pathogen sensitivity. Plant J. 44, 25–36. doi: 10.1111/j.1365-313X.2005.02508.x
Kroymann, J., Donnerhacke, S., Schnabelrauch, D., Mitchell-Olds, T. (2003). Evolutionary dynamics of an Arabidopsis insect resistance quantitative trait locus. Proc. Natl. Acad. Sci. U.S.A. 100, 14587–14592. doi: 10.1073/pnas.1734046100
Kono, M., Yamori, W., Suzuki, Y., Terashima, I. (2017). Photoprotection of PSI by far-red light against the fluctuating light-induced photoinhibition in Arabidopsis thaliana and field-grown plants. Plant Cell Physiol. 58, 35–45. doi: 10.1093/pcp/pcw215
Kudoh, H. (2016). Molecular phenology in plants: in natura systems biology for the comprehensive understanding of seasonal responses under natural environments. New Phytol. 210, 399–412. doi: 10.1111/nph.13733
Kuśnierczyk, A., Winge, P., Jrstad, T. S., Troczyska, J., Rossiter, J. T., Bones, A. M. (2008). Towards global understanding of plant defence against aphids timing and dynamics of early Arabidopsis defence responses to cabbage aphid (Brevicoryne brassicae) attack. Plant Cell Environ. 31, 1097–1115. doi: 10.1111/j.1365-3040.2008.01823.x
Langmead, B., Trapnell, C., Pop, M., Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25. doi: 10.1186/gb-2009-10-3-r25
Lambrix, V., Reichelt, M., Mitchell-Olds, T., Kliebenstein, D. J., Gershenzon, J. (2001). The Arabidopsis epithiospecifier protein promotes the hydrolysis of glucosinolates to nitriles and influences Trichoplusia ni herbivory. Plant Cell 13, 2793–2807. doi: 10.1105/tpc.13.12.2793
Larkin, J. C., Walker, J. D., Bolognesi-Winfield, A. C., Gray, J. C., Walker, A. R. (1999). Allele-specific interactions between ttg and gl1 during trichome development in Arabidopsis thaliana. Genetics 151, 1591–1604.
Li, B., Dewey, C. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323. doi: 10.1186/1471-2105-12-323
Lin, C. W., Huang, L. Y., Huang, C. L., Wang, Y. C., Lai, P. H., Wang, H. V., et al. (2017). Common stress transcriptome analysis reveals functional and genomic architecture differences between early and delayed response genes. Plant Cell Physiol. 58, 546–559. doi: 10.1093/pcp/pcx002
Lin, Y., Jiang, L., Chen, Q., Li, Y., Zhang, Y., Luo, Y., et al. (2018). Comparative transcriptome profiling analysis of red-and white-fleshed strawberry (Fragaria × ananassa) provides new insight into the regulation of anthocyanins pathway. Plant Cell Physiol. 59, 1844–1859. doi: 10.1093/pcp/pcy098
Martínez-Ballesta, M., Moreno-Fernández, D. A., Castejón, D., Ochando, C., Morandini, P. A., Carvajal, M. (2015). The impact of the absence of aliphatic glucosinolates on water transport under salt stress in Arabidopsis thaliana. Front. Plant Sci. 6, 524. doi: 10.3389/fpls.2015.00524
Martinoia, E., Mimura, T., Hara-Nishimura, I., Shiratake, K. (2018). The multifaceted roles of plant vacuoles. Plant Cell Physiol. 59, 1285–1287. doi: 10.1093/pcp/pcy113
Matthes, M. C., Bruce, T. J., Ton, J., Verrier, P. J., Pickett, J. A., Napier, J. A. (2010). The transcriptome of cis-jasmone-induced resistance in Arabidopsis thaliana and its role in indirect defence. Planta 232, 1163–1180. doi: 10.1007/s00425-010-1244-4
Matthes, M. C., Bruce, T. J., Chamberlain, K., Pickett, J. A., Napier, J. A. (2011). Emerging roles in plant defense for cis-jasmone-induced cytochrome P450 CYP81D11. Plant Signal. Behav. 6, 563–565. doi: 10.4161/psb.6.4.14915
Mewis, I., Appel, H. M., Hom, A., Raina., R., Schultz, J. C. (2005). Major signaling pathways modulate Arabidopsis glucosinolate accumulation and response to both phloem-feeding and chewing insects. Plant Physiol. 138, 1149–1162. doi: 10.1104/pp.104.053389
Mewis, I., Khan, M. A. M., Glawischnig, E., Schreiner, M., Ulrichs, C. (2012). Water stress and aphid feeding differentially influence metabolite composition in Arabidopsis thaliana (L.). PLoS ONE 7, e48661. doi: 10.1371/journal.pone.0048661
Mishra, N., Sun, L., Zhu, X., Smith, J., Prakash Srivastava, A., Yang, X., et al. (2017). Overexpression of the rice SUMO E3 ligase gene OsSIZ1 in cotton enhances drought and heat tolerance, and substantially improves fiber yields in the field under reduced irrigation and rainfed conditions. Plant Cell Physiol. 58, 735–746. doi: 10.1093/pcp/pcx032
Nagano, A. J., Honjo, M. N., Mihara, M., Sato, M., Kudoh, H., (2015). “Detection of plant viruses in natural environments by using RNA-Seq,” in Plant Virology Protocols (New York, NY: Humana Press), 89–98. doi: 10.1007/978-1-4939-1743-3_8
Nagano, A. J., Kawagoe, T., Sugisaka, J., Honjo, M. N., Iwayama, K., Kudoh, H. (2019). Annual transcriptome dynamics in natural environments reveals plant seasonal adaptation. Nat. Plants. 5, 74–83. doi: 10.1038/s41477-018-0338-z
Nakano, M., Mukaihara, T. (2018). Ralstonia solanacearum type III effector RipAL targets chloroplasts and induces jasmonic acid production to suppress salicylic acid-mediated defense responses in plants. Plant Cell Physiol. 59, 2576–2589. doi: 10.1093/pcp/pcy177
R Core Team. (2015). R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.
Ratzka, A., Vogel, H., Kliebenstein, D. J., Mitchell-Olds, T., Kroymann, J. (2002). Disarming the mustard oil bomb. Proc. Natl. Acad. Sci. U.S.A. 99, 11223–11228. doi: 10.1073/pnas.172112899
Renwick, J. A. A., Haribal, M., Gouinguené, S., Städler, E. (2006). Isothiocyanates stimulating oviposition by the diamondback moth, Plutella xylostella. J. Chem. Ecol. 32, 755–766. doi: 10.1007/s10886-006-9036-9
Sato, Y. (2018). Associational effects and the maintenance of polymorphism in plant defense against herbivores: review and evidence. Plant Spec. Biol. 33, 91–108. doi: 10.1111/1442-1984.12201
Sato, Y., Kudoh, H. (2017). Fine-scale frequency differentiation along a herbivory gradient in the trichome dimorphism of a wild Arabidopsis. Ecol. Evol. 7, 2133–2141. doi: 10.1002/ece3.2830
Sato, Y., Shimizu-Inatsugi, R., Yamazaki, M., Shimizu, K. K., Nagano, A. J. (2019a). Plant trichomes and a single gene GLABRA1 contribute to insect community composition on field-grown Arabidopsis thaliana. BMC Plant Biol. 19, 163. doi: 10.1186/s12870-019-1705-2
Sato, Y., Tezuka, A., Kashima, M., Deguchi, A., Shimizu-Inatsugi, R., Yamazaki, M., et al. (2019b). Transcriptional variation in glucosinolate biosynthetic genes and inducible responses to aphid herbivory on field-grown Arabidopsis thaliana. bioRxiv 563486. doi: 10.1101/563486
Schoonhoven, L. M., van Loon, J. J. A., Dicke, M., (2005). Insect-plant biology. 2nd ed. Oxford, UK: Oxford University Press.
Schmid, M., Davison, T. S., Henz, S. R., Pape, U. J., Demar, M., Vingron, M., et al. (2005). A gene expression map of Arabidopsis thaliana development. Nat. Genet. 37, 501–506. doi: 10.1038/ng1543
Shimizu, K. K., Kudoh, H., Kobayashi, M. J. (2011). Plant sexual reproduction during climate change: gene function in natura studied by ecological and evolutionary systems biology. Ann. Bot. 108, 777–787. doi: 10.1093/aob/mcr180
Shirakawa, M., Hara-Nishimura, I. (2018). Specialized vacuoles of myrosin cells: chemical defense strategy in Brassicales plants. Plant Cell Physiol. 59, 1309–1316. doi: 10.1093/pcp/pcy082
Snoeren, T. A. L., Kappers, I. F., Broekgaarden, C., Mumm, R., Dicke, M., Bouwmeester, H. J. (2010). Natural variation in herbivore-induced volatiles in Arabidopsis thaliana. J. Exp. Bot. 61, 3041–3056. doi: 10.1093/jxb/erq127
Sønderby, I. E., Hansen, B. G., Bjarnholt, N., Ticconi, C., Halkier, B. A., Kliebenstein, D. J. (2007). A systems biology approach identifies a R2R3 MYB gene subfamily with distinct and overlapping functions in regulation of aliphatic glucosinolates. PLoS ONE 2, e1322. doi: 10.1371/journal.pone.0001322
Sugiyama, A., Yamazaki, Y., Hamamoto, S., Takase, H., Yazaki, K. (2017). Synthesis and secretion of isoflavones by field-grown soybean. Plant Cell Physiol. 58, 1594–1600. doi: 10.1093/pcp/pcx084
Sun, H., Chen, L., Li, J., Hu, M., Ullah, A., He, X., et al. (2017). The JASMONATE ZIM-Domain gene family mediates JA signaling and stress response in cotton. Plant Cell Physiol. 58, 2139–2154. doi: 10.1093/pcp/pcx148
Taylor, M. A., Cooper, M. D., Sellamuthu, R., Braun, P., Migneault, A., Browning, A., et al. (2017). Interacting effects of genetic variation for seed dormancy and flowering time on phenology, life history, and fitness of experimental Arabidopsis thaliana populations over multiple generations in the field. New Phytol. 216, 291–302. doi: 10.1111/nph.14712
Thompson, L. (1994). The spatiotemporal effects of nitrogen and litter on the population dynamics of Arabidopsis thaliana. J. Ecol. 82, 63–88. doi: 10.2307/2261386
Tsuda, K. (2017). Division of tasks: defense by the spatial separation of antagonistic hormone activities. Plant Cell Physiol. 59, 3–4. doi: 10.1093/pcp/pcx208
Wang, M., Gu, Z., Wang, R., Guo, J., Ling, N., Firbank, L. G., et al. (2018). Plant primary metabolism regulated by nitrogen contributes to plant–pathogen interactions. Plant Cell Physiol. 60, 329–342. doi: 10.1093/pcp/pcy211
Wentzell, A. M., Rowe, H. C., Hansen, B. G., Ticconi, C., Halkier, B. A., Kliebenstein, D. J. (2007). Linking metabolic QTLs with network and cis-eQTLs controlling biosynthetic pathways. PLoS Genet. 3, e162. doi: 10.1371/journal.pgen.0030162
Wilczek, A. M., Roe, J. L., Knapp, M. C., Cooper, M. D., Lopez-Gallego, C., Martin, L. J., et al. (2009). Effects of genetic perturbation on seasonal life history plasticity. Science 323, 930–934. doi: 10.1126/science.1165826
Wu, S., Tohge, T., Cuadros-Inostroza, Á., Tong, H., Tenenboim, H., Kooke, R., et al. (2018). Mapping the Arabidopsis metabolic landscape by untargeted metabolomics at different environmental conditions. Mol. Plant 11, 118–134. doi: 10.1016/j.molp.2017.08.012
Xu, B., Yu, G., Li, H., Xie, Z., Wen, W., Zhang, J., et al. (2018). Knockdown of STAYGREEN in perennial ryegrass (Lolium perenne L.) leads to transcriptomic alterations related to suppressed leaf senescence and improved forage quality. Plant Cell Physiol. 60, 202–212. doi: 10.1093/pcp/pcy203
Yamaji, Y., Maejima, K., Komatsu, K., Shiraishi, T., Okano, Y., Himeno, M., et al. (2012). Lectin-mediated resistance impairs plant virus infection at the cellular level. Plant Cell 24, 778–793. doi: 10.1105/tpc.111.093658
Yamasaki, E., Altermatt, F., Cavender-Bares, J., Schuman, M. C., Zuppinger-Dingley, D., Garonna, I., et al. (2018). Genomics meets remote sensing in global change studies: monitoring and predicting phenology, evolution and biodiversity. Curr. Opin. Environ. Sustain. 29, 177–186. doi: 10.1016/j.cosust.2018.03.005
Yano, S. (1994). Ecological and evolutionary interactions between wild crucifers and their herbivorous insects. Plant Spec. Biol. 9, 137–143. doi: 10.1111/j.1442-1984.1994.tb00094.x
Yazaki, K., Arimura, G. I., Ohnishi, T. (2017). ‘Hidden’ terpenoids in plants: their biosynthesis, localization and ecological roles. Plant Cell Physiol. 58, 1615–1621. doi: 10.1093/pcp/pcx123
Yoshida, Y., Sano, R., Wada, T., Takabayashi, J., Okada, K. (2009). Jasmonic acid control of GLABRA3 links inducible defense and trichome patterning in Arabidopsis. Development 136, 1039–1048. doi: 10.1242/dev.030585
Zaidem, M., Groen, S. C., Purugganan, M. D. (2019). Evolutionary and ecological functional genomics, from lab to the wild. Plant J. 97, 40–55. doi: 10.1111/tpj.14167
Zhang, Z., Ober, J. A., Kliebenstein, D. J. (2006). The gene controlling the quantitative trait locus EPITHIOSPECIFIER MODIFIER1 alters glucosinolate hydrolysis and insect resistance in Arabidopsis. Plant Cell 18, 1524–1536. doi: 10.1105/tpc.105.039602
Zhang, X., Ivanova, A., Vandepoele, K., Radomiljac, J., de Velde, J., Berkowitz, O., et al. (2017). The transcription factor MYB29 is a regulator of ALTERNATIVE OXIDASE1a. Plant Physiol. 173, 1824–1843. doi: 10.1104/pp.16.01494
Zhou, S., Richter, A., Jander, G. (2018). Beyond defense: multiple functions of benzoxazinoids in maize metabolism. Plant Cell Physiol. 59, 1528–1537. doi: 10.1093/pcp/pcy064
Zhu, J., Wang, X., Guo, L., Xu, Q., Zhao, S., Li, F., et al. (2018). Characterization and alternative splicing profiles of the lipoxygenase gene family in tea plant (Camellia sinensis). Plant Cell Physiol. 59, 1765–1781. doi: 10.1093/pcp/pcy091
Keywords: AOP3, in natura, Lipaphis erysimi, RNA-Seq, plant–insect interaction
Citation: Sato Y, Tezuka A, Kashima M, Deguchi A, Shimizu-Inatsugi R, Yamazaki M, Shimizu KK. and Nagano AJ. (2019) Transcriptional Variation in Glucosinolate Biosynthetic Genes and Inducible Responses to Aphid Herbivory on Field-Grown Arabidopsis thaliana. Front. Genet. 10:787. doi: 10.3389/fgene.2019.00787
Received: 03 May 2019; Accepted: 25 July 2019;
Published: 11 September 2019.
Edited by:
James J. Cai, Texas A&M University, United StatesReviewed by:
Philippe Reymond, Université de Lausanne, SwitzerlandGeorg Jander, Boyce Thompson Institute, United States
Copyright © 2019 Sato, Tezuka, Kashima, Deguchi, Shimizu-Inatsugi, Yamazaki, Shimizu and Nagano. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kentaro K. Shimizu, kentaro.shimizu@ieu.uzh.ch; Atsushi J. Nagano, anagano@agr.ryukoku.ac.jp