bs_bs_banner
Plant, Cell and Environment (2014) 37, 2086–2101
doi: 10.1111/pce.12296
Original Article
Transcriptome analysis of soybean lines reveals transcript diversity and genes involved in the response to common cutworm (Spodoptera litura Fabricius) feeding Yongli Wang1,2, Hui Wang2, Rui Fan3, Qing Yang1 & Deyue Yu2 1 College of Life Sciences, 2National Center for Soybean Improvement, National Key Laboratory of Crop Genetics and Germplasm Enhancement, Nanjing Agricultural University, Nanjing, Jiangsu 210095, China and 3College of Chemical and Material Engineering, Quzhou University, Quzhou, Zhejiang 324000, China
ABSTRACT The interaction between soybeans and the destructive common cutworm insect is complicated. In this paper, the time course of induced responses to common cutworm was characterized in two soybean lines, and the results showed that the induced resistance peaked at different times in the resistant (WX) and susceptible (NN) soybean lines. Two sets of transcriptome profiles from the WX and NN lines at the peak of their induced resistance were compared using microarray analysis. In total, 827 and 349 transcripts were differentially expressed in the WX and NN lines, respectively, with 80 probes common regulated and seven regulated in the opposite direction. All common- and unique-regulated genes were grouped into 10 functional categories based on sequence similarity searches, which showed that most of the genes were related to stress and defence responses. qRTPCR analysis of 22 genes confirmed the results of the microarray analysis. The spatiotemporal expression patterns of the six genes revealed the consistency of systemic expression levels with the timing of the resistance response observed in the bioassay experiments. In summary, we described the conceptual model of induced resistance in two soybean lines and provided the first large-scale survey of common cutworm-induced defence transcripts in soybean. Key-words: induced resistance; microarray; spatiotemporal analysis.
INTRODUCTION Soybean (Glycine max L. Merrill) is a main oil-seed crop worldwide. The demand for soybean has grown rapidly because of its use in a wide range of applications in food, feed and industrial products (Wilson 2008). However, insect pests greatly increase the cost of soybean production. Meanwhile, concerns over insecticide residues in food and in the environment are becoming increasingly prevalent (McPherson 2004). Common cutworms can feed on members of several Correspondence: Q. Yang. e-mail: [email protected]; and D. Yu. Fax: +86 25 843 96410; e-mail: [email protected] 2086
taxa, and this pest has become one of the primary defoliators of soybean plants in southern China in the recent years (Cui & Gai 1997; Zhan et al. 2001). As generalist herbivore insects, common cutworms are sensitive to induced plant defences, and therefore, they serve as a useful model for research into insect-induced resistance (Stotz et al. 2002). Plant resistance to herbivores can be classified as constitutive resistance or induced resistance. Constitutive resistance always exists in plants, whether the attacker present or absent, whereas induced resistance is activated or enhanced only after a plant is attacked (Zhang et al. 2008). The benefits of induced defences against herbivores are their reduced cost and higher efficacy compared with constitutive resistance (Karban et al. 1997). Induced resistance to different attackers can be variable in different species or genotypes. Some plants only slightly respond to herbivory, whereas others produce high levels of toxins. For example, there are remarkable genotypic differences in the induced resistance of Populus (Havill & Raffa 1999; Osier & Lindroth 2001). Additionally, the level and effect of a herbivore-induced response may exhibit great temporal and spatial variation. Thus, plants may respond to injury within hours or days, with relatively short-lasting effects (Edwards et al. 1992; Baldwin & Schmelz 1996), whereas long-term induced responses can affect herbivores in subsequent seasons (Haukioja & Niemelä 1979; Hunter 1987; Bryant et al. 1991). Within a damaged plant, elevated levels of defence chemicals are commonly detected close to site of the initial injury (Croxford et al. 1989; Edwards et al. 1992), but systemic responses at sites far from the initial injury have also been demonstrated (Karban & Myers 1989; Alborn et al. 1996). The induced defence mechanisms of plants can be classified as direct defences that result in the inhibition of herbivore growth and/or development, or indirect defences, which include the release of plant volatiles that will attract the herbivore’s parasitoids and predators (Paré & Tumlinson 1999; Walling 2000). Complex cross-talk networks have been uncovered and have been shown to connect various signalling pathways in the regulation of defence induction (Walling 2000). Although methyl jasmonate (MeJA) signalling plays a © 2014 John Wiley & Sons Ltd
Transcriptional response to common cutworm feeding in soybeans primary role in chewing insect defence (Reymond et al. 2004; De Vos et al. 2006), ethylene-mediated expression is also involved (Stotz et al. 2000; Kessler & Baldwin 2002; von Dahl & Baldwin 2007). Additionally, salicylic acid (SA) is important to activate plant defence responses against pathogen attack during biotrophic pathogen interactions (McDowell & Dangl 2000; Glazebrook 2005). Accumulation of these hormones results in the activation of defence-related genes, followed by the production of various metabolic defences. Both positive and negative interactions between these pathways have been reported (Rojo et al. 2003; Bostock 2005; Beckers & Spoel 2006; Pieterse et al. 2009). Ecologists have long understood that plants exhibit multimechanistic resistance towards herbivores, but the molecular mechanisms underlying these complicated responses have remained elusive (Baldwin et al. 2001). Because the tools to unravel some of the signalling pathways involved in the regulation of gene expression in response to insect attacks have become available, the results of these studies have opened up new and unexpected areas of research. For example, a single microarray-based study revealed that Arabidopsis undergoes changes in levels of over 700 mRNAs during the defence response (Schenk et al. 2000). In contrast, only 100 mRNAs were up-regulated by spider mite infestation in lima beans, although a further 200 mRNAs were up-regulated in an indirect response mediated by feeding-induced volatile signal molecules (Arimura et al. 2000). However, many of these genes are of unknown function, and many changes in gene expression do not represent the induction of defence-related proteins. Insect-induced responses in soybean have only been studied for a small number of genes that are involved in secondary metabolism and octadecanoid signalling (Jouanin et al. 1998; Marchetti et al. 2000; Zhu-Salzman et al. 2003; Wu et al. 2008). We previously identified an induced resistance response caused by the feeding of common cutworms in a susceptible soybean line, and we detected 11 differently expressed proteins that may be involved in this induced defence response using proteomic approaches (Fan et al. 2012). However, the biochemical mechanisms and signalling pathways of the induced response in soybeans that leads to the decreased herbivore preference and/or performance are still poorly understood. In this study, two different local soybean lines, Wanxianbaidongdou (WX) and Nannong99-10 (NN), which are constitutively resistant and susceptible to common cutworm, respectively (Wang et al. 2011) were chosen to test the progression of induced resistance levels after insect feeding. Additionally, two sets of induced transcriptome profiles from WX and NN lines at their peak induction were compared using microarray analysis. This comparison allowed us to reveal molecular signatures underlying the different induced resistance phenotypes of these two soybean lines, and these results may facilitate a better understanding of the resistance mechanisms of soybean in response to common cutworm. qRT-PCR validation of 22 randomly selected genes from the two transcription profiles revealed a significantly strong correlation with the microarray results in both lines and confirmed that findings of the microarray data
2087
were reliable and rigorous. Furthermore, the spatial and temporal expression analysis of several insect defence-related genes provided us with a more revealing detailed kinetic analysis and revealed that genes expressed in different parts of the plant respond differently to insect stress across various time points. Our study provided broad insight into the herbivore-induced resistance in soybeans and created an excellent dataset of potential candidate genes involved in dynamic induced resistance in resistant and susceptible soybean lines which can then be tested further to determine their roles in this biological process using functional genomics approaches.
MATERIALS AND METHODS Plant growth, treatment and harvest Seeds from the WX and NN soybean lines were sown in 15 cm diameter plastic pots that contained a sterile soil mixture (top soil : sand : vermiculite, 3:2:1). The seedlings were grown in a greenhouse maintained at 30 ± 5 °C, with 70 ± 10% room humidity, and a 14 h/10 h (day/night) photoperiod with supplementary metal-halide illumination. To prevent the interference of other insects, all potted plants were covered with gauze. The potted plants were watered once per day with the same volume of water. Common cutworms (provided by Institute of Plant Protection, Jiangsu Academy of Agricultural Sciences, Nanjing, China) were applied to the soybean plants at the V4 stage (vegetative stage with four nodes) as defined by Fehr et al. (1971); control plants were left without exposure to common cutworms. Defoliation (approximately 20–30%) was induced by placing three third or fourth instar common cutworm larvae on the first and second trifoliate leaf for 48 h, after which the larvae were removed. The leaves of the control and treated plants were excised at six sampling times [0, 0.5, 1, 3, 5 and 7 days after induction (dai)] for the identification of induced resistance. For the transcriptome analysis and the first qRT-PCR validation experiments, the damaged and undamaged leaves mixed together in the treated plants and leaves at the same places in the control plants were taken at 5 and 1 dai from WX and NN lines, respectively. For spatiotemporal expression analysis, RNA samples were separately isolated from the upper undamaged leaves and from the middle damaged leaves in the treated and control plants at the same two places at 0, 0.5, 1, 3, 5, 7 and 15 dai.
Identification of induced resistance To evaluate the levels of soybean resistance induced by common cutworm, a dual-choice test and a force-feeding experiment were conducted. In the dual-choice test, three to five leaves from each pair of experimental plants were arranged opposite each other at the margins of a porcelain dish (30 cm × 20 cm). Ten common cutworm larvae, which had been starved for 24 h, were released into the middle area of the porcelain dish and were allowed to feed for 12 h. The leaf area was measured using a LICOR-3000 area meter
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
2088
Y. Wang et al.
(LI-COR, Lincoln, NE, USA). The consumed area in the control (C) and treated (T) was used to calculate the feeding preference index (PI), where PI = 2C/(T + C) (Kogan 1972). The PI values ranged from 0 to 2, with a PI = 1 indicating no feeding preference for either the control or treated leaves, a PI > 1 indicating a preference for the control leaves, and a PI < 1 indicating a preference for the treated leaves. In the force-feeding experiment, the control and treated leaves were placed in 15 cm diameter Petri dishes, along with five common cutworm larvae. In addition to the PI values, the relative growth rate (RGR) was also calculated using the formula RGR = (W1 −W2)/W1, where W1 and W2 are the weights of the five common cutworm larvae before being placed in the dish and after feeding for 24 h, respectively. Statistical analysis of the data was performed using Microsoft Excel 2007. Each experiment was repeated three times, and the mean values and standard deviation (SD) were calculated. Analysis of variance (GLM procedure of SAS; SAS 1989) was used to determine whether the genotypes differed significantly in the timing of their induced response, with sample time, genotype and methods included as factors.
Total RNA isolation RNA from frozen soybean leaves sampled for the transcriptome and qRT-PCR experiments was isolated using the RNeasy plant mini kit (QIAGEN Inc., Valencia, CA, USA) according to the manufacturer’s protocol. One volume of plant RNA isolation aid (Ambion, Austin, TX, USA) was added for each unit of mass of frozen tissue (mL g−1) before the tissue homogenization step to remove common contaminants, such as polysaccharides and polyphenolics. The RNA was treated with RNase-Free DNase set (QIAGEN Inc.) to digest any genomic DNA which might be present. RNA was quantified using an ultraviolet (UV) spectrophotometer and the RNA quality and integrity were examined using an RNA Lab-On-A-Chip (Caliper Technologies Corp., Mountain View, CA, USA), which was evaluated on an Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA). The samples were stored at −80 °C.
CapitalBio Corporation (http://www.capitalbio.com) according to the Affymetrix standard protocols. The data from these 12 chips are publicly available in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) under accession number GSE51823. Expression signals were first analysed using the GeneChip operating software 1.4 (GCOS, Affymetrix Inc., Santa Clara, CA, USA) to determine the ‘present’ probe set list. The MAS5.0 Algorithm in the Affymetrix’s GCOS was used to scale and normalize the data to provide signal value intensities, as well as to logarithm base-2 transform the data. Pairwise comparisons were performed between the groups (the common cutworm-induced sample versus the control non-induced sample), and genes with at least twofold differential expression and a P value of less than 0.05 were selected for processing using a forward step-wise false discovery rate (FDR) method in the Significance Analysis of Microarrays (SAM 3.0) software, an FDR of 5% was used. Genes passing the P value after FDR analysis were categorized as statistically significant differentially expressed genes. Annotation in the SoyBase for the Affymetrix Soybean Genome Array is limited, and most of the probes represent unknown genes. Therefore, we used the accession numbers of the probes to BLAST search the NCBI mRNA reference database to predict the putative identifications of the genes, an e-values cut-off of 1 × e−20 was used. The 1096 probes and corresponding accession numbers with putative identifications are listed in Supporting Information Table S1. Some identical genes were found to be represented by more than one probe. Identical probes that represent the same or homologous genes were merged into one putative identification. We then consolidated each list of identifications or genes into 10 broad categories based on putative function, including signalling, binding and transport, defence/stress, transcriptional regulation, secondary metabolism, cell wall modification, growth and development, primary metabolism, others and proteins of unknown function (Fig. 4). It should be noted that some of the genes’ functional categorizations might be arbitrary, and there may be some overlaps. All together, 64 probes from WX and 26 probes from NN failed to find putative identifications in the NCBI reference database with e-values less than 1 × e−20.
Microarray experimental design and data analysis
Real-time qRT-PCR
For the 5 and 1 dai samples of WX and NN, three biological replicates were used for each of the four treatments (WX treated, WX control, NN treated, NN control), requiring 12 soybean Affymetrix GeneChips, which contained 37 744 G. max probe sets (35 611 transcripts). Thus, some redundancy was present. This high-density array consists of an 11-probe pair (25 bp per oligonucleotide) with an 11-μm feature size. The details of the GeneChip® soybean genome array are available at the Affymetrix® website (http://affymetrix.com/ index.affx). Methods for the preparation of cDNA from mRNA, as well as the subsequent steps leading to the hybridization and scanning of the soybean GeneChip arrays, were performed as described by Ko & Han (2004) and Panthee et al. (2007). Microarray hybridization was conducted by the
A total of 22 genes were selected for qRT-PCR validation from the genes of interest. The housekeeping soybean b-tubulin gene (GenBank accession no. AY907703) was used as a reference gene for relative quantification (Chen et al. 2009). Target gene-specific primers were designed for qRTPCR using the online software Primer3 version 0.4.0 (Rozen & Skaletsky 2000). Reverse transcription reactions were performed with an oligo (dT) 18 primer and ReverTra Ace Moloney murine leukemia virus reverse transcriptase (Toyobo, Osaka, Japan) according to the manufacturer’s instructions. All primers were designed based on the sequences obtained using the BLASTN search in the NCBI database or soybean genome sequences database (http:// www.phytozome.net). Primer sequences can be found in
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
Transcriptional response to common cutworm feeding in soybeans
2089
Figure 1. Dynamics of induced resistance based on feeding preference indices of leaves area in both dual-choice (dotted line) and force-feeding (full line) tests across time after common cutworm feeding for WX (a) and NN (b) soybean lines. A preference index (PI) of 1 indicates no preference (no induced response). PI values > 1 indicate induced resistance and PI values < 1 indicate induced susceptibility. Relative growth rate of common cutworms fed with control and treated leaves in force-feeding tests across time after common cutworm feeding for WX (c) and NN (d) soybean lines. Error bars represent the standard deviation. *Significant (P < 0.05). **Highly significant (P < 0.01).
Supporting Information Table S5. qRT-PCR was performed on an ABI 7500 real-time PCR system (Applied Biosystems, Forster City, CA, USA) using the SYBR Green Realtime Master Mix (Toyobo). The data were analysed with the SDS 2.0 software (Applied Biosystems). The concordance of microarray expression data in the resistant and susceptible lines with the qRT-PCR analysis was measured using Pearson’s correlation based on log2 fold changes. Pearson’s correlation analyses were performed using the PASW Statistic 17.0 software suite (SPSS Inc., Chicago, IL, USA).
experiment at six different time points (Fig. 1c,d).An analysis of variance test based on PIs indicated that induced resistance did change over time (Table 1; ‘time’). The induced responses between the WX and NN soybean lines were highly different across the various time points (Table 1; ‘Varieties’ or ‘Varieties*time’). However, there was no difference in timing of responses between the two methods (Table 1; ‘Methods’). The induced resistance response was relatively rapid in the NN line, with the highest level of resistance observed at 1 dai; the response tapered off only by 3 dai (Fig. 1b). Between 5 and
RESULTS
Table 1. Analysis of variance for induction as a function of time since damage and varieties by two methods, dual-choice test and force-feeding experiment
Dynamic analysis of induced resistance against common cutworm in WX and NN To characterize the progression of the soybean-induced response to common cutworm herbivory, the levels of induced resistance were measured in dual-choice test and the forcefeeding experiment at six time points in one discrete event following damage.PIs of leaf area in the two experiments were calculated at six different times after induction (Fig. 1a,b).The RGR of the common cutworm larvae that were fed on treated and control leaves was also calculated in the force-feeding
Source
d.f.
Type III SS
F
F0.05
F0.01
Varieties Methods Time Varieties*time Error Total
1 1 5 5 46 71
0.45 0.01 1.39 2.54 0.31 5.20
66.85** 1.82 41.53** 75.93**
4.05 4.05 2.42 2.42
7.22 7.22 3.44 3.44
Highly significant (P < 0.01).
**
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
2090
Y. Wang et al.
7 dai,there was a period of increased susceptibility to common cutworm (Fig. 1b). In WX, however, the level of induced resistance appeared unstable during the first 24 h, followed by a subsequent increase in induced resistance from 3 to 5 dai. By 7 dai, the induced resistance had decreased, but it still existed (Fig. 1a). Furthermore, in the force-feeding experiment, the RGR of the common cutworm feeding on the control leaves was significantly (P < 0.01) higher than those feeding on the treated leaves at 5 and 1 dai in the WX and NN lines, respectively (Fig. 1c,d). Additionally, the induced resistance level of the WX and NN soybean lines at 5 and 1 dai from the dualchoice and force-feeding tests were independently compared. The PI of the leaf area was higher in the WX plants than in the NN plants, and the RGRs of the common cutworm larvae were lower in the resistant lines, regardless of whether the common cutworms were fed on control or treated leaves. The weight of larvae that were fed on the WX treated leaves did not increase, but rather slightly declined (Supporting Information Fig. S1).
Transcriptome changes in the resistant and susceptible lines in response to common cutworm feeding
Figure 2. Hierarchical clustering of 1096 differentially expressed
The common phenomenon was that no matter in resistant or susceptible line, there always exist a peak level of induced resistance at some time point. Thus, we focused on the transcriptome changes induced by insect feeding at the peak time points in the WX and NN lines, which could facilitate a better understanding of the common and special induced resistance mechanisms at the molecular level. Of the 37 744 probe sets in the Affymetrix Soybean Genome Array, 26 004 and 26 559 probe sets were separately detected in the WX and NN samples. We focused on probe sets that were at least twofold differentially expressed with P values less than 0.05 after 5% FDR processing. In total, 827 and 349 transcripts were determined to be differentially expressed in the WX and NN lines, respectively (Fig. 2).A total of 1096 probe sets were determined to be differentially expressed genes either between treated and control WX plants or between treated and control NN plants with three biological replicates. The hierarchical clustering analysis classified the resistant and susceptible soybean lines into two clusters, regardless of whether they were induced by common cutworms. The three biological replicates of the four treatments all clustered together, indicated that there was no difference among the biological replicates (Fig. 2). This result confirmed the robustness of the experimental design. Overall, 827 probes were up-regulated and no probes were downregulated in the WX line. Only 249 probes were up-regulated in the NN line, of which 73 probes showed the same expression trends in the WX line. Moreover, there were an additional 100 probes that were down-regulated in the NN line, of which seven probes were up-regulated in the WX line (Fig. 3).
probe sets of 12 chips representing four treatments [treated resistant soybean line (WX-T) and control (WX-CK), and treated susceptible soybean line (NN-T) and control (NN-CK)] with three biological replicates. Red and green indicate high and low expression levels.
after insect feeding. We then consolidated the 95 identifications into 10 broad categories based on their putative function (Fig. 4a). Most of the identifications were up-regulated in both soybean lines, except for the probes that represented several genes, which were mainly involved in signalling (Nos. 10, 11); binding and transport (Nos. 14, 15, 18) and primary metabolism (Nos. 59, 61, 63) were down-regulated in the NN line (Table 2).The major category was genes whose functions have not yet been ascertained, which were categorized as ‘unknown functions’ (30%). The fold changes of this category ranged from 2.19 to 62.33 for the WX line and from 2.17 to 38.23 for the NN line, respectively. Genes associated with pest/ pathogen defence, abiotic stress responses and detoxification/ redox processes were grouped under the defence/stress category (19%), which formed the second largest group, and
Functional classification of common- and unique-regulated genes in the two soybean lines
Figure 3. Venn diagram showing the number of
A total of 182 and 143 probes in the WX and NN lines, representing 95 identifications, were common response genes
up-regulated/down-regulated/oppositely regulated genes in WX and NN soybean lines in treated plants compared with control plants identified by microarray data.
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
Transcriptional response to common cutworm feeding in soybeans
2091
Figure 4. Functional classification of gene lists common or unique regulated in WX and NN. Pie charts represent functional classification according to putative identifications of (a) common regulated genes in WX and NN; (b) unique up-regulated genes in WX; (c) unique up-regulated genes in NN; (d) unique down-regulated genes in NN.
the fold changes in this group ranged from 2.1 to 24.42 up-regulated in the WX line and 0.41 to 19.46 up- or downregulated in the NN line.The third most significant fraction of genes that were common to both lines were genes involved in signalling category (13%), and the fold changes of these genes ranged from 2 to 7.07 and 0.38 to 6.36 for the WX and NN lines, respectively. Except for the signalling category, the other categories related to regulatory functions were the binding and transport category (6%) and transcriptional regulation category (4%). Several transcripts were associated with secondary metabolism (9%), including metabolites such as protease inhibitor, polyphenol oxide (PPO), isoflavone and phytoalexin, which may serve as both antimicrobials and antioxidants. The rest of the overlapping identifications were related to cell wall modification (7%), growth and development (3%), primary metabolism (5%) and others (4%). Our comparative transcriptome analysis of the two lines overlapped to a certain extent, but more differences were observed in gene regulation between the two lines. Of all the regulated probes, 581 probes representing 426 identifications were particularly up-regulated in WX, 115 probes representing 94 identifications were uniquely up-regulated in NN and 65 probes representing 60 identifications were uniquely down-regulated in NN. The identifications uniquely up- or
down-regulated in the WX and NN lines are listed in Supporting Information Tables S2, S3 and S4. We further grouped the unique identifications into the same 10 functional categories (Fig. 4b,c,d). The unique identifications that were up-regulated in WX were enriched for the following categories: unknown function (45%), transcriptional regulation (14%) and binding and transport (9%). These categories included numerous transcription factors (TFs), components of molecular complexes and metal ion transporters. Other complementary categories were defence/stress (6%), primary metabolism (6%), signalling (5%), secondary metabolism (4%), cell wall modification (3%), growth and development (3%) and others (5%). Pathogenesis-related proteins (PR1, PR2 and PR3) and genes related to programmed cell death and water stress were uniquely responded to the insect feeding in the WX line. Additionally, a wide variety of enzymes involved in the phenylpropanoid pathway and the biosynthesis of secondary metabolites, such as anthocyanin, alkaloid, lignin and volatile compounds, were also noteworthy. However, in the susceptible line, uniquely up-regulated probes were more corresponded to other categories:unknown function (46%), defence/stress (12%), secondary metabolism (10%) and binding and transport (10%). Examples of genes
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
2092
Y. Wang et al.
Table 2. Functional classification of identified genes overlapping between the resistant and susceptible lines after common cutworm induction compared with control plants
Categorizations Signalling
Binding and transport
Defence/stress
Transcriptional regulation
Secondary metabolism
Cell wall modification
Growth and development
Primary metabolism
Others
Unknown function aWe
1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68–95.
Resistant line (WX)
Susceptible line (NN)
(Putative) identificationsa
Nos. of probes
Fold change range
Nos. of probes
Fold change range
Lipoxygenase Allene oxide synthase Allene oxide cyclase TIFY 10A 1-aminocyclopropane-1-carboxylate oxidase Tryptophan aminotransferase-related protein Gibberellin receptor GID1B-like Phytosulfokine Two-component response regulator ARR5 CBL-interacting protein kinase S-receptor-like kinase Serine/threonine protein kinase Small ubiquitin-related modifier 2-like Peroxisomal small heat shock protein ABC transporter family members Bidirectional sugar transporters Nitrate transporters UDP-glycosyltransferase PR5 MLP-like protein 34-like SRG1-like protein LURP-one-related 15-like Disease resistance protein Patatin group A-3-like Stress-induced protein SAM22 Carboxylesterase 31 kDa protein /vegetative storage protein B 24 kDa seed coat protein Wound-induced protein Differentially expressed in soybean Senescence/dehydration-associated protein Peroxidases Glutathione S-transferase Gamma-glutamyl transpeptidase 1-like Glutaredoxin Heme oxygenase WRKY 40 WRKY56 BTB/POZ domain-containing protein Ethylene-responsive transcription factor Polyphenol oxidase Protease inhibitor 1 Bowman-Birk type proteinase inhibitor-like Flavonol synthase Isoflavone synthase Isoflavone 2′-hydroxylase-like Isoflavone reductase Momilactone A synthase-like Xyloglucan endotransglucosylase Laccases Hydroxyproline-rich glycoprotein Periaxin-like Repetitive proline-rich cell wall protein Glucuronoxylan glucuronosyltransferases Polyamine oxidases O-acyltransferase WSD1-like Microtubule-associated protein LOB domain-containing protein Sucrose synthase GDSL esterase/lipases Esterase-like Basic 7S globulin-like Aspartic proteinase nepenthesin 21 kDa protein-like CASP-like protein NAD Cytochrome P450 78A3-like Glycine max cDNA homologues
1 1 1 3 5 1 1 1 3 6 3 3 1 1 7 1 1 8 2 1 1 1 3 2 5 3 1 2 2 2 1 11 1 2 2 1 5 1 2 3 1 1 1 1 1 1 2 2 2 3 1 1 3 1 1 3 1 2 2 1 2 3 5 2 5 2 2 28
5.36 2.3 2.43 2.34–5.46 2.21–4.53 2.02 2.53 2.13 2–2.04 2.16–2.37 2.02–2.7 2.04–7.07 2.5 3.95 2.39–7.62 4.27 2.79 2.34–10 3.44–6.69 6.51 2.29 3.32 2.1–3.57 3.17–3.7 18.46–22.42 3.54–5.29 4.67 6.48–6.67 6.56–6.61 8.5–12.6 2.54 2.36–24.42 2.2 2.22–4.21 2.53–4.06 14.22 4.89–5.27 5.27 2.07–2.19 2.28–9.27 3.45 3.48 11.62 2.09 2.42 7.17 3.41–4.91 3.02–4.04 2.1–6.83 3.2–40.12 19.72 19.72 5.06–8.55 4.19 8.09 2.97–16.25 3.22 2.76–4.89 2.64–3.26 20 2.12 2.19–4.72 2.43–6.58 2.36–2.49 2.01–2.96 3.26–4.78 2.58–3.69 2.19–62.33
1 2+2 1 3 1 1 1 1 5 2 1 1+2 1 1 1 2 3 1 1 1 1 1 1 1 1 1 1 2 1 2 1 10 3 2 2 1 2 1 1 1 5 1 1 1+1 1 1 1 1 1+1 1 1 1 3 1 2 1 2 1 3 2 2 3 1 1 1+1 1 1 28
6.09 2.23–3.12/0.45–0.44 2.4 2.03–3.51 3.61 6.36 2.44 2.44 2.02–8.07 0.38–0.42 0.49 2.82/0.5 2.07 0.28 0.45 4.55–5.81 2.03–2.21 0.39 3.03 0.41 3.79 0.49 2.8 2.4 2.54 2.31 2.03 2.7–2.73 2.38 6.87–13.48 3.4 2.64–19.46 2.13–3 2.27–2.52 2.52–2.63 5.94 2.86–3.81 3.18 0.46 3.36 3.47–16.55 38.23 236.11 3.58/0.49 0.4 16.06 2.99 17.6 2.01/0.31 3.82 2.2 9.04 2.9–12.09 0.31 2.08–2.1 3.87 2.58–2.59 3.28 0.32–0.33/2.45 2.35–2.56 0.34–0.47 3.54–4.16 0.4 0.43 2.32/0.44 3.98 2.04 2.17–38.23
identified genes whose e-values exhibited at least less than 1 × e−20 with the homology searches in DNA databases.
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
Transcriptional response to common cutworm feeding in soybeans found in these function categories are mainly encoded for PR8 proteins (chitinases), cysteine proteinase inhibitor, geranyl-diphosphate synthase, calnexin and the chaperonebinding proteins, BiP C and BiP B. The remaining less prominent functional categories were cell wall modification (5%), growth and development (5%), signalling (4%), others (4%), primary metabolism (3%) and transcriptional regulation (1%). Differences in the quantities of the various groups were particularly noted in the identifications of the genes that were suppressed in the susceptible line. Some functional categories were absent in this set of genes, including transcriptional regulation (0%), cell wall modification (0%) and secondary metabolism (0%). Three functional categories were noticeably enriched: unknown function (63%), primary metabolism (13%) and defence/stress (8%). More than half of the unique identifications in the NN line were proteins with unknown functions. Interestingly, the second largest category was the primary metabolism category, which mainly includes enzymes involved in catalyzing electrogenic translocation, carbohydrate synthesis and partitioning in the photosynthetic pathway.
Assessing the accuracy of mRNA measurements with qRT-PCR To validate the microarray expression profiles and to obtain more refined gene expression data, we conducted qRT-PCR for a subset of 22 genes selected from the probe sets that were uniquely induced or suppressed, as well as from the probes that overlapped between theWX and NN lines.These 22 genes were from various functional categories that may have potential biological significance in plant defence against insects. Collectively, the fold changes observed for 20 genes in our qRT-PCR expression analysis of the treated samples and control samples of the resistant and susceptible lines were statistically significant (P < 0.05) or highly significant (P < 0.01). The P values of the qRT-PCR fold changes of calnexin and salicylic acid methyl transferase-like protein (SAMT) genes in the WX line were 0.17 and 0.38, respectively. Similarly, two other genes that encode dehydrationresponsive element binding protein and the MYB73 transcription factor were also insignificant (P > 0.05) between the treated and control samples in the NN line. Furthermore, the magnitude of the expression levels was slightly different between the microarray and real-time RT-PCR analyses (Supporting Information Table S6). However, as indicated by previous studies (Yauk & Berndt 2007), we concluded that more emphasis should be placed on the direction of the gene expression change indicated by the microarray than the magnitude of change. Moreover, we conducted a correlation analysis of the microarray and qRT-PCR expression data in the WX and NN lines. According to our results, there was a slightly better correlation between the microarray and qRT-PCR values for the NN samples (R = 0.898, P = 0.000, Fig. 5a) than for the WX samples (R = 0.843, P = 0.000, Supporting Information Fig. S2). Overall, the microarray and qRT-PCR data were strongly correlated in a highly significant (P < 0.01) manner for both the resistant and susceptible lines.
2093
Spatiotemporal expression analysis To obtain a detailed spatiotemporal expression pattern of selected genes and to further research the spatial and temporal extent of the common cutworm-induced resistance in the resistant and susceptible soybean lines, we selected the following genes for spatiotemporal expression analysis: (1) insect defence-related genes: soybean vegetative lectin (SVL, EU070415.1) and putative cysteine proteinase inhibitor gene (CysPl3, XM_003545396.1); (2) phenylpropanoid pathway-related genes: isoflavone synthase 1 gene (IFS1, AF195798.1) and NADPH:isoflavone reductase gene (N:IFR, AJ003245.1); (3) wounding-induced genes: vegetative storage protein β gene (VSPβ, M20038.1) and peroxisomal small heat shock protein gene (pHSP, NM_001253036.1). The expression pattern of SVL in the WX line demonstrated a strongly enhanced expression level at 0 dai in the damaged leaves (local) after treatment; however, in the undamaged leaves (systemic), the expression alterations of SVL were approximately threefold less than in the local response (Fig. 5a). The maximum up-regulation of SVL in the WX plants was observed at 0 and 5 dai for the local and systemic responses, respectively. By contrast, in the NN line, there were three- and 4.6-fold increases in SVL expression beginning at 0 dai, and these levels quickly reached six- to sevenfold by 0.5 dai for the local and systemic responses, respectively, this increase was subsequently followed by a general decrease in expression. Unlike SVL, a slightly higher up-regulation of the CysPl3 expression levels was observed in the damaged leaves compared with the undamaged leaves in the WX line at the beginning 1 dai (Fig. 5b), after which both the local and systemic responses were weakened. The expression of CysPl3 in response to herbivory treatment in the NN line was much stronger than in the WX line during the first 24 h. Similarly, both local and systemic responses were weakened after 1 dai. The expression of IFS1 was induced in the WX line but repressed in the NN line after insect treatment. In WX, there was a characteristic delay in IFS1 expression before the maximal levels were reached at 1 dai, expression was maintained at a high level until 5 dai, after which the expression decreased.The magnitude of the local IFS1 expression was, in most cases, stronger than systemic IFS1 expression, suggesting that the accumulation of IFS1 mRNAs is primarily a local response to herbivore damage in the WX line. On the contrary, in NN, the expression of IFS1 was suppressed at a similar magnitude in both the local and systemic leaves from 0 to 7 dai, however, expression did recover to its normal level at 15 dai (Fig. 5c). The putative phytoalexin biosynthesisrelated gene N:IFR responded slowly in the WX line by increasing linearly during the first 3 dai, with a higher level of local expression compared with systemic expression. After reaching a peak at 3 dai, N:IFR expression slightly decreased each day, especially in the local leaves. However, in the NN line, N:IFR was quickly and dramatically induced at the early time points, reaching maximal levels of 5.3- and 5.5-fold increases compared with the control plants at 1 dai in the local and systemic leaves. The high expression level of N:IFR
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
2094
Y. Wang et al.
Figure 5. Spatial and temporal analysis of six selected genes in seven time points after induction. The values are expressed as log2 of fold changes compared with control at the same time point of each one, >0 means up-regulation,
Plant, Cell and Environment (2014) 37, 2086–2101
doi: 10.1111/pce.12296
Original Article
Transcriptome analysis of soybean lines reveals transcript diversity and genes involved in the response to common cutworm (Spodoptera litura Fabricius) feeding Yongli Wang1,2, Hui Wang2, Rui Fan3, Qing Yang1 & Deyue Yu2 1 College of Life Sciences, 2National Center for Soybean Improvement, National Key Laboratory of Crop Genetics and Germplasm Enhancement, Nanjing Agricultural University, Nanjing, Jiangsu 210095, China and 3College of Chemical and Material Engineering, Quzhou University, Quzhou, Zhejiang 324000, China
ABSTRACT The interaction between soybeans and the destructive common cutworm insect is complicated. In this paper, the time course of induced responses to common cutworm was characterized in two soybean lines, and the results showed that the induced resistance peaked at different times in the resistant (WX) and susceptible (NN) soybean lines. Two sets of transcriptome profiles from the WX and NN lines at the peak of their induced resistance were compared using microarray analysis. In total, 827 and 349 transcripts were differentially expressed in the WX and NN lines, respectively, with 80 probes common regulated and seven regulated in the opposite direction. All common- and unique-regulated genes were grouped into 10 functional categories based on sequence similarity searches, which showed that most of the genes were related to stress and defence responses. qRTPCR analysis of 22 genes confirmed the results of the microarray analysis. The spatiotemporal expression patterns of the six genes revealed the consistency of systemic expression levels with the timing of the resistance response observed in the bioassay experiments. In summary, we described the conceptual model of induced resistance in two soybean lines and provided the first large-scale survey of common cutworm-induced defence transcripts in soybean. Key-words: induced resistance; microarray; spatiotemporal analysis.
INTRODUCTION Soybean (Glycine max L. Merrill) is a main oil-seed crop worldwide. The demand for soybean has grown rapidly because of its use in a wide range of applications in food, feed and industrial products (Wilson 2008). However, insect pests greatly increase the cost of soybean production. Meanwhile, concerns over insecticide residues in food and in the environment are becoming increasingly prevalent (McPherson 2004). Common cutworms can feed on members of several Correspondence: Q. Yang. e-mail: [email protected]; and D. Yu. Fax: +86 25 843 96410; e-mail: [email protected] 2086
taxa, and this pest has become one of the primary defoliators of soybean plants in southern China in the recent years (Cui & Gai 1997; Zhan et al. 2001). As generalist herbivore insects, common cutworms are sensitive to induced plant defences, and therefore, they serve as a useful model for research into insect-induced resistance (Stotz et al. 2002). Plant resistance to herbivores can be classified as constitutive resistance or induced resistance. Constitutive resistance always exists in plants, whether the attacker present or absent, whereas induced resistance is activated or enhanced only after a plant is attacked (Zhang et al. 2008). The benefits of induced defences against herbivores are their reduced cost and higher efficacy compared with constitutive resistance (Karban et al. 1997). Induced resistance to different attackers can be variable in different species or genotypes. Some plants only slightly respond to herbivory, whereas others produce high levels of toxins. For example, there are remarkable genotypic differences in the induced resistance of Populus (Havill & Raffa 1999; Osier & Lindroth 2001). Additionally, the level and effect of a herbivore-induced response may exhibit great temporal and spatial variation. Thus, plants may respond to injury within hours or days, with relatively short-lasting effects (Edwards et al. 1992; Baldwin & Schmelz 1996), whereas long-term induced responses can affect herbivores in subsequent seasons (Haukioja & Niemelä 1979; Hunter 1987; Bryant et al. 1991). Within a damaged plant, elevated levels of defence chemicals are commonly detected close to site of the initial injury (Croxford et al. 1989; Edwards et al. 1992), but systemic responses at sites far from the initial injury have also been demonstrated (Karban & Myers 1989; Alborn et al. 1996). The induced defence mechanisms of plants can be classified as direct defences that result in the inhibition of herbivore growth and/or development, or indirect defences, which include the release of plant volatiles that will attract the herbivore’s parasitoids and predators (Paré & Tumlinson 1999; Walling 2000). Complex cross-talk networks have been uncovered and have been shown to connect various signalling pathways in the regulation of defence induction (Walling 2000). Although methyl jasmonate (MeJA) signalling plays a © 2014 John Wiley & Sons Ltd
Transcriptional response to common cutworm feeding in soybeans primary role in chewing insect defence (Reymond et al. 2004; De Vos et al. 2006), ethylene-mediated expression is also involved (Stotz et al. 2000; Kessler & Baldwin 2002; von Dahl & Baldwin 2007). Additionally, salicylic acid (SA) is important to activate plant defence responses against pathogen attack during biotrophic pathogen interactions (McDowell & Dangl 2000; Glazebrook 2005). Accumulation of these hormones results in the activation of defence-related genes, followed by the production of various metabolic defences. Both positive and negative interactions between these pathways have been reported (Rojo et al. 2003; Bostock 2005; Beckers & Spoel 2006; Pieterse et al. 2009). Ecologists have long understood that plants exhibit multimechanistic resistance towards herbivores, but the molecular mechanisms underlying these complicated responses have remained elusive (Baldwin et al. 2001). Because the tools to unravel some of the signalling pathways involved in the regulation of gene expression in response to insect attacks have become available, the results of these studies have opened up new and unexpected areas of research. For example, a single microarray-based study revealed that Arabidopsis undergoes changes in levels of over 700 mRNAs during the defence response (Schenk et al. 2000). In contrast, only 100 mRNAs were up-regulated by spider mite infestation in lima beans, although a further 200 mRNAs were up-regulated in an indirect response mediated by feeding-induced volatile signal molecules (Arimura et al. 2000). However, many of these genes are of unknown function, and many changes in gene expression do not represent the induction of defence-related proteins. Insect-induced responses in soybean have only been studied for a small number of genes that are involved in secondary metabolism and octadecanoid signalling (Jouanin et al. 1998; Marchetti et al. 2000; Zhu-Salzman et al. 2003; Wu et al. 2008). We previously identified an induced resistance response caused by the feeding of common cutworms in a susceptible soybean line, and we detected 11 differently expressed proteins that may be involved in this induced defence response using proteomic approaches (Fan et al. 2012). However, the biochemical mechanisms and signalling pathways of the induced response in soybeans that leads to the decreased herbivore preference and/or performance are still poorly understood. In this study, two different local soybean lines, Wanxianbaidongdou (WX) and Nannong99-10 (NN), which are constitutively resistant and susceptible to common cutworm, respectively (Wang et al. 2011) were chosen to test the progression of induced resistance levels after insect feeding. Additionally, two sets of induced transcriptome profiles from WX and NN lines at their peak induction were compared using microarray analysis. This comparison allowed us to reveal molecular signatures underlying the different induced resistance phenotypes of these two soybean lines, and these results may facilitate a better understanding of the resistance mechanisms of soybean in response to common cutworm. qRT-PCR validation of 22 randomly selected genes from the two transcription profiles revealed a significantly strong correlation with the microarray results in both lines and confirmed that findings of the microarray data
2087
were reliable and rigorous. Furthermore, the spatial and temporal expression analysis of several insect defence-related genes provided us with a more revealing detailed kinetic analysis and revealed that genes expressed in different parts of the plant respond differently to insect stress across various time points. Our study provided broad insight into the herbivore-induced resistance in soybeans and created an excellent dataset of potential candidate genes involved in dynamic induced resistance in resistant and susceptible soybean lines which can then be tested further to determine their roles in this biological process using functional genomics approaches.
MATERIALS AND METHODS Plant growth, treatment and harvest Seeds from the WX and NN soybean lines were sown in 15 cm diameter plastic pots that contained a sterile soil mixture (top soil : sand : vermiculite, 3:2:1). The seedlings were grown in a greenhouse maintained at 30 ± 5 °C, with 70 ± 10% room humidity, and a 14 h/10 h (day/night) photoperiod with supplementary metal-halide illumination. To prevent the interference of other insects, all potted plants were covered with gauze. The potted plants were watered once per day with the same volume of water. Common cutworms (provided by Institute of Plant Protection, Jiangsu Academy of Agricultural Sciences, Nanjing, China) were applied to the soybean plants at the V4 stage (vegetative stage with four nodes) as defined by Fehr et al. (1971); control plants were left without exposure to common cutworms. Defoliation (approximately 20–30%) was induced by placing three third or fourth instar common cutworm larvae on the first and second trifoliate leaf for 48 h, after which the larvae were removed. The leaves of the control and treated plants were excised at six sampling times [0, 0.5, 1, 3, 5 and 7 days after induction (dai)] for the identification of induced resistance. For the transcriptome analysis and the first qRT-PCR validation experiments, the damaged and undamaged leaves mixed together in the treated plants and leaves at the same places in the control plants were taken at 5 and 1 dai from WX and NN lines, respectively. For spatiotemporal expression analysis, RNA samples were separately isolated from the upper undamaged leaves and from the middle damaged leaves in the treated and control plants at the same two places at 0, 0.5, 1, 3, 5, 7 and 15 dai.
Identification of induced resistance To evaluate the levels of soybean resistance induced by common cutworm, a dual-choice test and a force-feeding experiment were conducted. In the dual-choice test, three to five leaves from each pair of experimental plants were arranged opposite each other at the margins of a porcelain dish (30 cm × 20 cm). Ten common cutworm larvae, which had been starved for 24 h, were released into the middle area of the porcelain dish and were allowed to feed for 12 h. The leaf area was measured using a LICOR-3000 area meter
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
2088
Y. Wang et al.
(LI-COR, Lincoln, NE, USA). The consumed area in the control (C) and treated (T) was used to calculate the feeding preference index (PI), where PI = 2C/(T + C) (Kogan 1972). The PI values ranged from 0 to 2, with a PI = 1 indicating no feeding preference for either the control or treated leaves, a PI > 1 indicating a preference for the control leaves, and a PI < 1 indicating a preference for the treated leaves. In the force-feeding experiment, the control and treated leaves were placed in 15 cm diameter Petri dishes, along with five common cutworm larvae. In addition to the PI values, the relative growth rate (RGR) was also calculated using the formula RGR = (W1 −W2)/W1, where W1 and W2 are the weights of the five common cutworm larvae before being placed in the dish and after feeding for 24 h, respectively. Statistical analysis of the data was performed using Microsoft Excel 2007. Each experiment was repeated three times, and the mean values and standard deviation (SD) were calculated. Analysis of variance (GLM procedure of SAS; SAS 1989) was used to determine whether the genotypes differed significantly in the timing of their induced response, with sample time, genotype and methods included as factors.
Total RNA isolation RNA from frozen soybean leaves sampled for the transcriptome and qRT-PCR experiments was isolated using the RNeasy plant mini kit (QIAGEN Inc., Valencia, CA, USA) according to the manufacturer’s protocol. One volume of plant RNA isolation aid (Ambion, Austin, TX, USA) was added for each unit of mass of frozen tissue (mL g−1) before the tissue homogenization step to remove common contaminants, such as polysaccharides and polyphenolics. The RNA was treated with RNase-Free DNase set (QIAGEN Inc.) to digest any genomic DNA which might be present. RNA was quantified using an ultraviolet (UV) spectrophotometer and the RNA quality and integrity were examined using an RNA Lab-On-A-Chip (Caliper Technologies Corp., Mountain View, CA, USA), which was evaluated on an Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA). The samples were stored at −80 °C.
CapitalBio Corporation (http://www.capitalbio.com) according to the Affymetrix standard protocols. The data from these 12 chips are publicly available in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) under accession number GSE51823. Expression signals were first analysed using the GeneChip operating software 1.4 (GCOS, Affymetrix Inc., Santa Clara, CA, USA) to determine the ‘present’ probe set list. The MAS5.0 Algorithm in the Affymetrix’s GCOS was used to scale and normalize the data to provide signal value intensities, as well as to logarithm base-2 transform the data. Pairwise comparisons were performed between the groups (the common cutworm-induced sample versus the control non-induced sample), and genes with at least twofold differential expression and a P value of less than 0.05 were selected for processing using a forward step-wise false discovery rate (FDR) method in the Significance Analysis of Microarrays (SAM 3.0) software, an FDR of 5% was used. Genes passing the P value after FDR analysis were categorized as statistically significant differentially expressed genes. Annotation in the SoyBase for the Affymetrix Soybean Genome Array is limited, and most of the probes represent unknown genes. Therefore, we used the accession numbers of the probes to BLAST search the NCBI mRNA reference database to predict the putative identifications of the genes, an e-values cut-off of 1 × e−20 was used. The 1096 probes and corresponding accession numbers with putative identifications are listed in Supporting Information Table S1. Some identical genes were found to be represented by more than one probe. Identical probes that represent the same or homologous genes were merged into one putative identification. We then consolidated each list of identifications or genes into 10 broad categories based on putative function, including signalling, binding and transport, defence/stress, transcriptional regulation, secondary metabolism, cell wall modification, growth and development, primary metabolism, others and proteins of unknown function (Fig. 4). It should be noted that some of the genes’ functional categorizations might be arbitrary, and there may be some overlaps. All together, 64 probes from WX and 26 probes from NN failed to find putative identifications in the NCBI reference database with e-values less than 1 × e−20.
Microarray experimental design and data analysis
Real-time qRT-PCR
For the 5 and 1 dai samples of WX and NN, three biological replicates were used for each of the four treatments (WX treated, WX control, NN treated, NN control), requiring 12 soybean Affymetrix GeneChips, which contained 37 744 G. max probe sets (35 611 transcripts). Thus, some redundancy was present. This high-density array consists of an 11-probe pair (25 bp per oligonucleotide) with an 11-μm feature size. The details of the GeneChip® soybean genome array are available at the Affymetrix® website (http://affymetrix.com/ index.affx). Methods for the preparation of cDNA from mRNA, as well as the subsequent steps leading to the hybridization and scanning of the soybean GeneChip arrays, were performed as described by Ko & Han (2004) and Panthee et al. (2007). Microarray hybridization was conducted by the
A total of 22 genes were selected for qRT-PCR validation from the genes of interest. The housekeeping soybean b-tubulin gene (GenBank accession no. AY907703) was used as a reference gene for relative quantification (Chen et al. 2009). Target gene-specific primers were designed for qRTPCR using the online software Primer3 version 0.4.0 (Rozen & Skaletsky 2000). Reverse transcription reactions were performed with an oligo (dT) 18 primer and ReverTra Ace Moloney murine leukemia virus reverse transcriptase (Toyobo, Osaka, Japan) according to the manufacturer’s instructions. All primers were designed based on the sequences obtained using the BLASTN search in the NCBI database or soybean genome sequences database (http:// www.phytozome.net). Primer sequences can be found in
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
Transcriptional response to common cutworm feeding in soybeans
2089
Figure 1. Dynamics of induced resistance based on feeding preference indices of leaves area in both dual-choice (dotted line) and force-feeding (full line) tests across time after common cutworm feeding for WX (a) and NN (b) soybean lines. A preference index (PI) of 1 indicates no preference (no induced response). PI values > 1 indicate induced resistance and PI values < 1 indicate induced susceptibility. Relative growth rate of common cutworms fed with control and treated leaves in force-feeding tests across time after common cutworm feeding for WX (c) and NN (d) soybean lines. Error bars represent the standard deviation. *Significant (P < 0.05). **Highly significant (P < 0.01).
Supporting Information Table S5. qRT-PCR was performed on an ABI 7500 real-time PCR system (Applied Biosystems, Forster City, CA, USA) using the SYBR Green Realtime Master Mix (Toyobo). The data were analysed with the SDS 2.0 software (Applied Biosystems). The concordance of microarray expression data in the resistant and susceptible lines with the qRT-PCR analysis was measured using Pearson’s correlation based on log2 fold changes. Pearson’s correlation analyses were performed using the PASW Statistic 17.0 software suite (SPSS Inc., Chicago, IL, USA).
experiment at six different time points (Fig. 1c,d).An analysis of variance test based on PIs indicated that induced resistance did change over time (Table 1; ‘time’). The induced responses between the WX and NN soybean lines were highly different across the various time points (Table 1; ‘Varieties’ or ‘Varieties*time’). However, there was no difference in timing of responses between the two methods (Table 1; ‘Methods’). The induced resistance response was relatively rapid in the NN line, with the highest level of resistance observed at 1 dai; the response tapered off only by 3 dai (Fig. 1b). Between 5 and
RESULTS
Table 1. Analysis of variance for induction as a function of time since damage and varieties by two methods, dual-choice test and force-feeding experiment
Dynamic analysis of induced resistance against common cutworm in WX and NN To characterize the progression of the soybean-induced response to common cutworm herbivory, the levels of induced resistance were measured in dual-choice test and the forcefeeding experiment at six time points in one discrete event following damage.PIs of leaf area in the two experiments were calculated at six different times after induction (Fig. 1a,b).The RGR of the common cutworm larvae that were fed on treated and control leaves was also calculated in the force-feeding
Source
d.f.
Type III SS
F
F0.05
F0.01
Varieties Methods Time Varieties*time Error Total
1 1 5 5 46 71
0.45 0.01 1.39 2.54 0.31 5.20
66.85** 1.82 41.53** 75.93**
4.05 4.05 2.42 2.42
7.22 7.22 3.44 3.44
Highly significant (P < 0.01).
**
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
2090
Y. Wang et al.
7 dai,there was a period of increased susceptibility to common cutworm (Fig. 1b). In WX, however, the level of induced resistance appeared unstable during the first 24 h, followed by a subsequent increase in induced resistance from 3 to 5 dai. By 7 dai, the induced resistance had decreased, but it still existed (Fig. 1a). Furthermore, in the force-feeding experiment, the RGR of the common cutworm feeding on the control leaves was significantly (P < 0.01) higher than those feeding on the treated leaves at 5 and 1 dai in the WX and NN lines, respectively (Fig. 1c,d). Additionally, the induced resistance level of the WX and NN soybean lines at 5 and 1 dai from the dualchoice and force-feeding tests were independently compared. The PI of the leaf area was higher in the WX plants than in the NN plants, and the RGRs of the common cutworm larvae were lower in the resistant lines, regardless of whether the common cutworms were fed on control or treated leaves. The weight of larvae that were fed on the WX treated leaves did not increase, but rather slightly declined (Supporting Information Fig. S1).
Transcriptome changes in the resistant and susceptible lines in response to common cutworm feeding
Figure 2. Hierarchical clustering of 1096 differentially expressed
The common phenomenon was that no matter in resistant or susceptible line, there always exist a peak level of induced resistance at some time point. Thus, we focused on the transcriptome changes induced by insect feeding at the peak time points in the WX and NN lines, which could facilitate a better understanding of the common and special induced resistance mechanisms at the molecular level. Of the 37 744 probe sets in the Affymetrix Soybean Genome Array, 26 004 and 26 559 probe sets were separately detected in the WX and NN samples. We focused on probe sets that were at least twofold differentially expressed with P values less than 0.05 after 5% FDR processing. In total, 827 and 349 transcripts were determined to be differentially expressed in the WX and NN lines, respectively (Fig. 2).A total of 1096 probe sets were determined to be differentially expressed genes either between treated and control WX plants or between treated and control NN plants with three biological replicates. The hierarchical clustering analysis classified the resistant and susceptible soybean lines into two clusters, regardless of whether they were induced by common cutworms. The three biological replicates of the four treatments all clustered together, indicated that there was no difference among the biological replicates (Fig. 2). This result confirmed the robustness of the experimental design. Overall, 827 probes were up-regulated and no probes were downregulated in the WX line. Only 249 probes were up-regulated in the NN line, of which 73 probes showed the same expression trends in the WX line. Moreover, there were an additional 100 probes that were down-regulated in the NN line, of which seven probes were up-regulated in the WX line (Fig. 3).
probe sets of 12 chips representing four treatments [treated resistant soybean line (WX-T) and control (WX-CK), and treated susceptible soybean line (NN-T) and control (NN-CK)] with three biological replicates. Red and green indicate high and low expression levels.
after insect feeding. We then consolidated the 95 identifications into 10 broad categories based on their putative function (Fig. 4a). Most of the identifications were up-regulated in both soybean lines, except for the probes that represented several genes, which were mainly involved in signalling (Nos. 10, 11); binding and transport (Nos. 14, 15, 18) and primary metabolism (Nos. 59, 61, 63) were down-regulated in the NN line (Table 2).The major category was genes whose functions have not yet been ascertained, which were categorized as ‘unknown functions’ (30%). The fold changes of this category ranged from 2.19 to 62.33 for the WX line and from 2.17 to 38.23 for the NN line, respectively. Genes associated with pest/ pathogen defence, abiotic stress responses and detoxification/ redox processes were grouped under the defence/stress category (19%), which formed the second largest group, and
Functional classification of common- and unique-regulated genes in the two soybean lines
Figure 3. Venn diagram showing the number of
A total of 182 and 143 probes in the WX and NN lines, representing 95 identifications, were common response genes
up-regulated/down-regulated/oppositely regulated genes in WX and NN soybean lines in treated plants compared with control plants identified by microarray data.
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
Transcriptional response to common cutworm feeding in soybeans
2091
Figure 4. Functional classification of gene lists common or unique regulated in WX and NN. Pie charts represent functional classification according to putative identifications of (a) common regulated genes in WX and NN; (b) unique up-regulated genes in WX; (c) unique up-regulated genes in NN; (d) unique down-regulated genes in NN.
the fold changes in this group ranged from 2.1 to 24.42 up-regulated in the WX line and 0.41 to 19.46 up- or downregulated in the NN line.The third most significant fraction of genes that were common to both lines were genes involved in signalling category (13%), and the fold changes of these genes ranged from 2 to 7.07 and 0.38 to 6.36 for the WX and NN lines, respectively. Except for the signalling category, the other categories related to regulatory functions were the binding and transport category (6%) and transcriptional regulation category (4%). Several transcripts were associated with secondary metabolism (9%), including metabolites such as protease inhibitor, polyphenol oxide (PPO), isoflavone and phytoalexin, which may serve as both antimicrobials and antioxidants. The rest of the overlapping identifications were related to cell wall modification (7%), growth and development (3%), primary metabolism (5%) and others (4%). Our comparative transcriptome analysis of the two lines overlapped to a certain extent, but more differences were observed in gene regulation between the two lines. Of all the regulated probes, 581 probes representing 426 identifications were particularly up-regulated in WX, 115 probes representing 94 identifications were uniquely up-regulated in NN and 65 probes representing 60 identifications were uniquely down-regulated in NN. The identifications uniquely up- or
down-regulated in the WX and NN lines are listed in Supporting Information Tables S2, S3 and S4. We further grouped the unique identifications into the same 10 functional categories (Fig. 4b,c,d). The unique identifications that were up-regulated in WX were enriched for the following categories: unknown function (45%), transcriptional regulation (14%) and binding and transport (9%). These categories included numerous transcription factors (TFs), components of molecular complexes and metal ion transporters. Other complementary categories were defence/stress (6%), primary metabolism (6%), signalling (5%), secondary metabolism (4%), cell wall modification (3%), growth and development (3%) and others (5%). Pathogenesis-related proteins (PR1, PR2 and PR3) and genes related to programmed cell death and water stress were uniquely responded to the insect feeding in the WX line. Additionally, a wide variety of enzymes involved in the phenylpropanoid pathway and the biosynthesis of secondary metabolites, such as anthocyanin, alkaloid, lignin and volatile compounds, were also noteworthy. However, in the susceptible line, uniquely up-regulated probes were more corresponded to other categories:unknown function (46%), defence/stress (12%), secondary metabolism (10%) and binding and transport (10%). Examples of genes
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
2092
Y. Wang et al.
Table 2. Functional classification of identified genes overlapping between the resistant and susceptible lines after common cutworm induction compared with control plants
Categorizations Signalling
Binding and transport
Defence/stress
Transcriptional regulation
Secondary metabolism
Cell wall modification
Growth and development
Primary metabolism
Others
Unknown function aWe
1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68–95.
Resistant line (WX)
Susceptible line (NN)
(Putative) identificationsa
Nos. of probes
Fold change range
Nos. of probes
Fold change range
Lipoxygenase Allene oxide synthase Allene oxide cyclase TIFY 10A 1-aminocyclopropane-1-carboxylate oxidase Tryptophan aminotransferase-related protein Gibberellin receptor GID1B-like Phytosulfokine Two-component response regulator ARR5 CBL-interacting protein kinase S-receptor-like kinase Serine/threonine protein kinase Small ubiquitin-related modifier 2-like Peroxisomal small heat shock protein ABC transporter family members Bidirectional sugar transporters Nitrate transporters UDP-glycosyltransferase PR5 MLP-like protein 34-like SRG1-like protein LURP-one-related 15-like Disease resistance protein Patatin group A-3-like Stress-induced protein SAM22 Carboxylesterase 31 kDa protein /vegetative storage protein B 24 kDa seed coat protein Wound-induced protein Differentially expressed in soybean Senescence/dehydration-associated protein Peroxidases Glutathione S-transferase Gamma-glutamyl transpeptidase 1-like Glutaredoxin Heme oxygenase WRKY 40 WRKY56 BTB/POZ domain-containing protein Ethylene-responsive transcription factor Polyphenol oxidase Protease inhibitor 1 Bowman-Birk type proteinase inhibitor-like Flavonol synthase Isoflavone synthase Isoflavone 2′-hydroxylase-like Isoflavone reductase Momilactone A synthase-like Xyloglucan endotransglucosylase Laccases Hydroxyproline-rich glycoprotein Periaxin-like Repetitive proline-rich cell wall protein Glucuronoxylan glucuronosyltransferases Polyamine oxidases O-acyltransferase WSD1-like Microtubule-associated protein LOB domain-containing protein Sucrose synthase GDSL esterase/lipases Esterase-like Basic 7S globulin-like Aspartic proteinase nepenthesin 21 kDa protein-like CASP-like protein NAD Cytochrome P450 78A3-like Glycine max cDNA homologues
1 1 1 3 5 1 1 1 3 6 3 3 1 1 7 1 1 8 2 1 1 1 3 2 5 3 1 2 2 2 1 11 1 2 2 1 5 1 2 3 1 1 1 1 1 1 2 2 2 3 1 1 3 1 1 3 1 2 2 1 2 3 5 2 5 2 2 28
5.36 2.3 2.43 2.34–5.46 2.21–4.53 2.02 2.53 2.13 2–2.04 2.16–2.37 2.02–2.7 2.04–7.07 2.5 3.95 2.39–7.62 4.27 2.79 2.34–10 3.44–6.69 6.51 2.29 3.32 2.1–3.57 3.17–3.7 18.46–22.42 3.54–5.29 4.67 6.48–6.67 6.56–6.61 8.5–12.6 2.54 2.36–24.42 2.2 2.22–4.21 2.53–4.06 14.22 4.89–5.27 5.27 2.07–2.19 2.28–9.27 3.45 3.48 11.62 2.09 2.42 7.17 3.41–4.91 3.02–4.04 2.1–6.83 3.2–40.12 19.72 19.72 5.06–8.55 4.19 8.09 2.97–16.25 3.22 2.76–4.89 2.64–3.26 20 2.12 2.19–4.72 2.43–6.58 2.36–2.49 2.01–2.96 3.26–4.78 2.58–3.69 2.19–62.33
1 2+2 1 3 1 1 1 1 5 2 1 1+2 1 1 1 2 3 1 1 1 1 1 1 1 1 1 1 2 1 2 1 10 3 2 2 1 2 1 1 1 5 1 1 1+1 1 1 1 1 1+1 1 1 1 3 1 2 1 2 1 3 2 2 3 1 1 1+1 1 1 28
6.09 2.23–3.12/0.45–0.44 2.4 2.03–3.51 3.61 6.36 2.44 2.44 2.02–8.07 0.38–0.42 0.49 2.82/0.5 2.07 0.28 0.45 4.55–5.81 2.03–2.21 0.39 3.03 0.41 3.79 0.49 2.8 2.4 2.54 2.31 2.03 2.7–2.73 2.38 6.87–13.48 3.4 2.64–19.46 2.13–3 2.27–2.52 2.52–2.63 5.94 2.86–3.81 3.18 0.46 3.36 3.47–16.55 38.23 236.11 3.58/0.49 0.4 16.06 2.99 17.6 2.01/0.31 3.82 2.2 9.04 2.9–12.09 0.31 2.08–2.1 3.87 2.58–2.59 3.28 0.32–0.33/2.45 2.35–2.56 0.34–0.47 3.54–4.16 0.4 0.43 2.32/0.44 3.98 2.04 2.17–38.23
identified genes whose e-values exhibited at least less than 1 × e−20 with the homology searches in DNA databases.
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
Transcriptional response to common cutworm feeding in soybeans found in these function categories are mainly encoded for PR8 proteins (chitinases), cysteine proteinase inhibitor, geranyl-diphosphate synthase, calnexin and the chaperonebinding proteins, BiP C and BiP B. The remaining less prominent functional categories were cell wall modification (5%), growth and development (5%), signalling (4%), others (4%), primary metabolism (3%) and transcriptional regulation (1%). Differences in the quantities of the various groups were particularly noted in the identifications of the genes that were suppressed in the susceptible line. Some functional categories were absent in this set of genes, including transcriptional regulation (0%), cell wall modification (0%) and secondary metabolism (0%). Three functional categories were noticeably enriched: unknown function (63%), primary metabolism (13%) and defence/stress (8%). More than half of the unique identifications in the NN line were proteins with unknown functions. Interestingly, the second largest category was the primary metabolism category, which mainly includes enzymes involved in catalyzing electrogenic translocation, carbohydrate synthesis and partitioning in the photosynthetic pathway.
Assessing the accuracy of mRNA measurements with qRT-PCR To validate the microarray expression profiles and to obtain more refined gene expression data, we conducted qRT-PCR for a subset of 22 genes selected from the probe sets that were uniquely induced or suppressed, as well as from the probes that overlapped between theWX and NN lines.These 22 genes were from various functional categories that may have potential biological significance in plant defence against insects. Collectively, the fold changes observed for 20 genes in our qRT-PCR expression analysis of the treated samples and control samples of the resistant and susceptible lines were statistically significant (P < 0.05) or highly significant (P < 0.01). The P values of the qRT-PCR fold changes of calnexin and salicylic acid methyl transferase-like protein (SAMT) genes in the WX line were 0.17 and 0.38, respectively. Similarly, two other genes that encode dehydrationresponsive element binding protein and the MYB73 transcription factor were also insignificant (P > 0.05) between the treated and control samples in the NN line. Furthermore, the magnitude of the expression levels was slightly different between the microarray and real-time RT-PCR analyses (Supporting Information Table S6). However, as indicated by previous studies (Yauk & Berndt 2007), we concluded that more emphasis should be placed on the direction of the gene expression change indicated by the microarray than the magnitude of change. Moreover, we conducted a correlation analysis of the microarray and qRT-PCR expression data in the WX and NN lines. According to our results, there was a slightly better correlation between the microarray and qRT-PCR values for the NN samples (R = 0.898, P = 0.000, Fig. 5a) than for the WX samples (R = 0.843, P = 0.000, Supporting Information Fig. S2). Overall, the microarray and qRT-PCR data were strongly correlated in a highly significant (P < 0.01) manner for both the resistant and susceptible lines.
2093
Spatiotemporal expression analysis To obtain a detailed spatiotemporal expression pattern of selected genes and to further research the spatial and temporal extent of the common cutworm-induced resistance in the resistant and susceptible soybean lines, we selected the following genes for spatiotemporal expression analysis: (1) insect defence-related genes: soybean vegetative lectin (SVL, EU070415.1) and putative cysteine proteinase inhibitor gene (CysPl3, XM_003545396.1); (2) phenylpropanoid pathway-related genes: isoflavone synthase 1 gene (IFS1, AF195798.1) and NADPH:isoflavone reductase gene (N:IFR, AJ003245.1); (3) wounding-induced genes: vegetative storage protein β gene (VSPβ, M20038.1) and peroxisomal small heat shock protein gene (pHSP, NM_001253036.1). The expression pattern of SVL in the WX line demonstrated a strongly enhanced expression level at 0 dai in the damaged leaves (local) after treatment; however, in the undamaged leaves (systemic), the expression alterations of SVL were approximately threefold less than in the local response (Fig. 5a). The maximum up-regulation of SVL in the WX plants was observed at 0 and 5 dai for the local and systemic responses, respectively. By contrast, in the NN line, there were three- and 4.6-fold increases in SVL expression beginning at 0 dai, and these levels quickly reached six- to sevenfold by 0.5 dai for the local and systemic responses, respectively, this increase was subsequently followed by a general decrease in expression. Unlike SVL, a slightly higher up-regulation of the CysPl3 expression levels was observed in the damaged leaves compared with the undamaged leaves in the WX line at the beginning 1 dai (Fig. 5b), after which both the local and systemic responses were weakened. The expression of CysPl3 in response to herbivory treatment in the NN line was much stronger than in the WX line during the first 24 h. Similarly, both local and systemic responses were weakened after 1 dai. The expression of IFS1 was induced in the WX line but repressed in the NN line after insect treatment. In WX, there was a characteristic delay in IFS1 expression before the maximal levels were reached at 1 dai, expression was maintained at a high level until 5 dai, after which the expression decreased.The magnitude of the local IFS1 expression was, in most cases, stronger than systemic IFS1 expression, suggesting that the accumulation of IFS1 mRNAs is primarily a local response to herbivore damage in the WX line. On the contrary, in NN, the expression of IFS1 was suppressed at a similar magnitude in both the local and systemic leaves from 0 to 7 dai, however, expression did recover to its normal level at 15 dai (Fig. 5c). The putative phytoalexin biosynthesisrelated gene N:IFR responded slowly in the WX line by increasing linearly during the first 3 dai, with a higher level of local expression compared with systemic expression. After reaching a peak at 3 dai, N:IFR expression slightly decreased each day, especially in the local leaves. However, in the NN line, N:IFR was quickly and dramatically induced at the early time points, reaching maximal levels of 5.3- and 5.5-fold increases compared with the control plants at 1 dai in the local and systemic leaves. The high expression level of N:IFR
© 2014 John Wiley & Sons Ltd, Plant, Cell and Environment, 37, 2086–2101
2094
Y. Wang et al.
Figure 5. Spatial and temporal analysis of six selected genes in seven time points after induction. The values are expressed as log2 of fold changes compared with control at the same time point of each one, >0 means up-regulation,
