Recording SA and the SA-specific leaf transcriptome
To identify reliable SA marker genes in soybean we watered soil-grown plants at V2 stage with 150 ml of an aqueous SA solution (1 mM). We verified the uptake of SA by quantifying the level of free and total SA in the first trifoliate leaf by HPLC analysis. Figure 1a shows that the SA drench led to a ~400-fold increase in free SA within 6 hours with a maximum of ~50 µg g-1 free SA. At 6 h after treatment (hpt), the level of total SA was only slightly higher than that of free SA, suggesting that SA was mainly present in its free form (82%). The level of free SA was much lower at 24 and 48 hpt than at 6 hpt, but still 30-40-fold higher than in the water-treated controls (Fig. 1a). The decrease in free SA at later times coincided with a further increase in total SA up to ~100 µg g-1 FW at 48 hpt (Fig. 1a). At 24 and 48hpt bound SA (calculated as the difference between total and free SA) contributed the biggest part to total SA (~5% free SA and ~95 % bound SA) (Fig. 1a).
To gain comprehensive knowledge of the changes in gene expression following SA accumulation, we performed global transcriptome analyses. For read statistics of the RNA-seq analysis, please refer to Tables S1 and S2. The analysis identified 3,130 genes with induced expression upon SA treatment (Fig. 1c). Because we aimed to identify genes that may serve as a molecular marker to hormone treatment we here focused on genes that were activated upon accumulation of free SA in leaves but disregarded genes whose expression was suppressed by SA. Their identity is given in Figure S1. In close correlation with the high levels of free SA at 6 hpt and its low levels at 48 hpt (Fig. 1a), the overall response to SA was highest at 6 and 24 hpt with 1,541 and 1,568 significantly activated genes and lower at 48hpt with 1,014 significantly activated genes (Fig. 1c). While free SA levels severely dropped between 6 and 24 hpt (Fig. 1a), the number of SA-activated genes remained high during this period (Fig. 1c).
Recording JA and the JA-specific transcriptome in leaves
To identify JAs-responsive marker genes in soybean, we exposed plants to MeJA.. Subsequently we measured the accumulation of JA and its bioactive derivative JA-Ile 35,36 in leaves at various times after exposure. As seen in Figure 1b, the accumulation in leaves of JA and JA-Ile had an apparent maximum at 24 h after exposure to MeJA constituting ~40 µg g-1 dry weight (DW) for JA and 9 µg g-1 DW for JA-Ile (Fig. 1b). After 24 h, the MeJA-enhanced JA and JA-Ile levels dropped. However, the observed decrease in JA-Ile between 24 and 48 hpt was less pronounced than that observed for free JA (Figure 1b). The accumulation of mRNA transcript of GmVSP-B (Glyma.08g200100), which contains a MeJA-responsive domain in its promoter region 37, also pointed to similar kinetics of JA presence in MeJA-exposed plants. In fact, GmVSP-B expression was activated at all time points following MeJA exposure (Fig. S2).
Global transcriptome analysis of MeJA-exposed plants identified 1,459 genes whose expression was significantly activated by MeJA (Fig. 1d). Gene activation by MeJA followed JA-Ile accumulation (Fig. 1b and d). At 6 hpt, 559 genes were significantly more expressed in MeJA-treated plants than in the controls. The number of activated genes increased with time constituting 953 genes at the 48 hpt time point. The identity of genes whose expression was suppressed by MeJA treatments is given in Fig. S1.
Gene ontology (GO)-enrichment analysis
GO-enrichment analysis revealed that among the SA-induced genes those associated with systemic acquired resistance (SAR) were the most significantly enriched loci (Fig. 1e). However, various genes in some other SA signaling-related GO terms were enriched as well. They include “regulation of defense responses”, “regulation of plant hypersensitive response”, “SA-biosynthetic process”, “SA-mediated signaling pathway” and “defense response to fungus” (Fig. 1e).
GO categories that were overrepresented among the MeJA-induced genes comprised “response to wounding” and “response to chitin” (Fig. 1f), the latter also being enriched upon SA treatment (Fig. 1e). Among genes in other GO terms, MeJA activated the expression of genes in the “ethylene biosynthetic process”, “phenylpropanoid biosynthetic process”, “SA catabolic process” and “lipid oxidation” bins (Fig. 1f). In addition to the genes whose expression was specifically activated by either SA or MeJA, multiple genes in the GO categories associated with primary metabolism (e.g. photosynthesis) were repressed by any of the two hormones (Fig. S1). The overall pattern of gene activation upon treatments might reflect the reported growth-to-defense transition upon exposure to SA or JAs 38–40.
Identification of genes that are specifically activated by SA or JA
To identify genes whose expression is specifically activated by either SA or JAs we compared the gene expression profiles in trifoliate leaves of soybean plants at 6, 24 and 48 hpt with SA or MeJA and compared it to that of untreated control leaves. We then selected genes that were significantly induced (p ≤ 0.01) by only one of the two hormones. From these sets of genes, we discarded those genes that also were differentially expressed among mock-treated and untreated plants, independent of the time point of differential expression. To identify marker genes that, in addition to strong induction, would exhibit high expression upon treatment, we sorted identified genes by the difference in fragment per kilobase per million mapped reads (FPKM)-values between the untreated control and treatments (ΔFPKM). Figure 2 depicts genes whose transcript levels are stable in all control treatments but significantly increased to high levels only with exclusively one of the two hormones applied. These candidate genes were further divided into categories of top 15 genes with induced expression at all tested time points after treatment (category 1), early responsive genes with induced expression after 6 and 24 hpt (category 2), late responsive genes with increased expression after 24 and 48 hpt (category 3), and genes that were induced only at one single time point after treatment (6 or 24 or 48 hpt) (category 4) (Fig. 2). Consistent with the number of genes induced by SA or MeJA (Fig. 1c, d), the induced expression of genes by SA was more pronounced at the early than late time points. The opposite was true for the genes induced by MeJA that had their peak in expression at later time points (Fig. 2). For example, SA-activated genes were induced stronger in category 2 compared to category 3, whereas mRNA transcript abundance of MeJA-activated genes was mostly lower in category 2 than in category 3 (Fig. 2). Taken together, the results in Figure 2 revealed that exogenously applied SA and MeJA elicit distinct temporal gene expression patterns in soybean, with SA-responsive genes being activated preferably at the early and MeJA-induced genes at the later times upon treatment.
Verifying specific responsiveness of identified marker gene candidates to SA or JAs
To evaluate the suitability of identified genes as specific molecular markers for SA or MeJA, we profiled the expression of each 5 selected candidate genes by qRT-PCR analysis. We focused on genes with strongest and/or longest-lasting activation upon treatment because we considered those genes as probably the most robust marker genes for the hormone in question. As seen in Figures 2 and 3, detected expression kinetics of the selected genes were highly consistent between the RNA-seq and qRT-PCR analyses. Only SA, but neither control or MeJA treatment affected the expression of the selected SA-marker gene candidates GmNIMIN1 (Glyma.10G010100),GmNIMIN1.2 (Glyma.02G009500), GmWRKY40 (Glyma.17G222500) GmUGT (Glyma.02G029900) and GmGH3 (Glyma.17G165300) (Fig. 3a). Expression of those genes was highly activated by SA with fold changes varying from ~200 to >2000-fold (Fig. 3a) when compared to the untreated controls. In a similar manner, expression of the identified JAs-responsive genes GmBPI1 (Glyma.18G231700),GmKTI1 (Glyma.08G342000),GmAAT (Glyma.10G109500), GmCYP79B2 (Glyma.11G197300) and GmG3PA (Glyma.10G119900) was triggered by MeJA, but not SA or mock treatment, with fold-changes in expression from ~60 to >2,000-fold over the untreated controls (Fig. 3b).
Using identified marker genes to deciphersoybean defense responses
To further verify the suitability of identified marker genes for analyzing recruitment of SA and/or JAs signaling in the soybean, we next compared the responsiveness of newly identified SA-marker genes and the level of SA in plants during a compatible interaction of soybean cultivar Williams 82 (W82) with the biotrophic fungal pathogen Pp (isolate BR05). To increase the chance of detecting possible transient accumulation of the hormone in question, we focused on marker genes in category 1, because genes in this category are induced for a longer period (Figs. 2 and 3). Both, GmNIMIN1 and GmWRKY40 showed a biphasic induction of expression in inoculated leaves of susceptible W82 plants with an apparent first expression maximum at 12-24 hpi and a second expression peak from 144 hpi whereas GmNIMIN1.2 was expressed at higher levels in Pp-challenged than mock-inoculated leaves at all time points assayed (Fig. 4a). However, differences in expression of all three identified SA-marker genes among infected and mock-inoculated leaves were most distinct at 144 hpt. In contrast to the SA marker genes, expression of the top JAs-marker gene GmBPI1 was not detected at any tested time point upon Pp inoculation (Fig. 4a).
To find out whether expression of presumed SA-marker genes would correlate with endogenous accumulation of SA, we quantified SA in leaves of Pp-infected and mock-inoculated soybean W82 plants. The highest content of free SA (~ 9.5 µg/g DW) was detected at 144 hpi after Pp inoculation (Fig. 4b). Free SA also accumulated early after inoculation (12-24 hpi). However, free SA levels at these time points remained below the accumulation at 144 hpi (Fig. 4b). At 72 hpi, however, free SA did not accumulate in response to Pp inoculation compared to mock controls. This observation is consistent with the expression of newly identified soybean SA-marker genes GmNIMIN1 and GmWRKY40 (Fig. 4a, b). In line with the absence of GmBPI1 expression, accumulation of the most active JAs derivate JA-Ile was also absent at any time point tested after the Pp inoculation (Fig. 4a, b). Differential accumulation of free SA between mock and Pp-inoculated W82 leaves thus correlated with the abundance of mRNA transcript of adequate hormone marker genes (Fig. 4a, b).
Insect feeding often leads to the accumulation of JAs and to the activation of JAs-responsive defense genes in plants 41–43. Spodoptera exigua is a generalist herbivore that feeds on soybean and many other plant species 44,45. A detached leaf assay disclosed that feeding of S. exigua larvae caused activation of GmAAT, GmBPI1 and GmKTI1 gene expression in soybean (Fig. 4b). The expression of the three JAs-marker genes had an apparent maximum at 72-96 hpt (Fig. 4a). Of the three JAs-responsive genes, GmBPI1 exhibited the strongest and most prolonged expression (Fig. 4a). To assess the specificity of identified marker genes for JA and SA we also quantified the level of transcript accumulation of the top SA-marker gene GmNIMIN1 (Fig. 3). In contrast to the JAs-marker genes, GmNIMIN1 was not much responsive to feeding of soybean W82 leaflets with S. exigua larvae (Fig. 4a). In fact, expression of the gene was only slightly activated at 48 hpt upon feeding (Fig. 4a). Quantification of JA-Ile revealed a strong accumulation only after insect feeding with an apparent maximum of ~1.4 µg/g DW at 96 hpi (Fig. 4b). In contrast, free SA accumulated to higher levels in mock-treated leaves than upon insect-feeding (Fig. 4b). Thus, GmNIMIN1, GmNIMIN1.2, GmWRKY40, GmBPI1, GmKTI1 and GmAAT serve as reliable molecular markers to indirectly monitor hormonal changes upon insect feeding on soybean plants.
Identification of genes for normalizing gene expression
For the precise determination of gene expression by qRT-PCR analysis and for adequate scientific conclusions reliable reference genes are required which are stably expressed across different experimental conditions. To identify such reference genes in soybean, we searched for genes that would display minimal variation in expression in all samples of our transcriptome dataset (untreated, mock-treated, SA/MeJA-treated). By doing so we identified 6 candidate genes with minimal changes in expression over all the samples (Table S3). The stability of expression of the genes was compared to the widely used reference gene GmUBQ 46–49. Normfinder analysis 50 disclosed that the 6 candidate reference genes were expressed with lower variably than GmUBQ in trifoliate leaves (Table S3). Glyma.08G211200 was identified as the most stably expressed gene across all conditions tested in our analysis (Table S3). Normalizing the expression of the SA- and JA-marker genes GmNIMIN1 and GmBPI1 to the expression of GmUBQ and Glyma.08G211200, respectively, revealed a similar pattern of expression (Fig S3). However, relative expression values were higher when normalized to the expression of Glyma.08G211200 than GmUBQ3 due to the overall lower expression level of Glyma.08G211200 compared to GmUBQ3 (Fig. S3). Our data, therefore, recommend Glyma.08G211200 for use as a robust and lower expressed alternative to highly expressed GmUBQ for normalizing qRT-PCR-based evaluation of gene expression in soybean leaves, at least in the conditions used in this study.