Introduction

MicroRNAs (miRNAs) are small noncoding RNAs that exist in most eukaryotic genomes (Carthew and Sontheimer 2009; Bartel 2004; Baulcombe 2004). They play a fundamental role in almost all biological processes, including development, hormone signaling and response to biotic and abiotic stresses (Bartel 2004; Rubio-Somoza et al. 2009; Jin 2008; Sunkar et al. 2007; Padmanabhan et al. 2009). miRNAs function through partial complementary base pairing to messenger RNAs (mRNAs), resulting mainly in post-transcriptional regulation of the targeted genes. Information of miRNA abundance under various conditions can thus provide insight to their regulatory functions.

Several studies have shown that small RNAs play a critical role in disease resistance responses (Jin 2008; Padmanabhan et al. 2009; Navarro et al. 2006; Katiyar-Agarwal et al. 2006, 2007; Gibbings and Voinnet 2010). miR393 can be induced by bacterial elicitor flg22 and positively contributes to pathogen associated molecular pattern (PAMP)-triggered immunity (PTI) by silencing auxin receptors and subsequently suppressing auxin signaling (Navarro et al. 2006). Auxin promotes plant growth, which provides carbon and nitrogen sources for biotrophic pathogens and increase their virulence and susceptibility. In addition, auxin can also promote pathogenesis by suppressing salicylic acid (SA)-mediated defense responses (Chen et al. 2007; Achard et al. 2007; Grant and Jones 2009; Spoel and Dong 2008). Small RNA profiling analysis using deep sequencing revealed that in addition to miR393, miR167 and miR160 were also induced by a non-pathogenic Pseudomonas syringae pv. tomato (Pst) DC3000 strain with a mutated type III secretion system hrcC (Fahlgren et al. 2007). Profiling of AGO1-bound small RNA after flg22 treatment revealed the role of miR160a, miR398b and miR773 in callose deposition, indicating their involvement in PAMP triggered immunity (Li et al. 2010). miR398 is involved in both biotic and abiotic responses (Sunkar et al. 2006; Jagadeeswaran et al. 2009). It is down-regulated by bacterial pathogen Pst DC3000 (avrRpt2) or Pst (avrRpm1). In A. tumefaciens, introduced T-DNA genes can trigger RNA silencing, which was then suppressed by an oncogenic protein. Accumulation of miR393 and miR167 was greatly reduced in C58-induced tumors in plants (Dunoyer et al. 2006). Furthermore, induced host pathogenesis-related (PR)-gene expression and TMV resistance in A. tumefaciens were correlated with an increased expression of miR393 in the infiltrated area with C58 but not those with a disarmed control (Pruss et al. 2008). In Pinus taeda, the expression of 10 tested miRNA (pta-miRNA) families were significantly repressed in the fusiform rust fungus (Cronartium quercuum) infected galled stem compared to healthy stem (Lu et al. 2007). Furthermore, the expression of ~82 plant disease-related transcripts was detected to be altered in response to miRNA regulation in pine (Lu et al. 2007). In Brassica, miR1885, which was predicted to target TIR-NBS-LRR class genes, was induced in Turnip mosaic virus infected plants (He et al. 2008).

In this paper, we studied the expression of the miRNAs in Arabidopsis in response to infection of non-pathogenic, virulent and avirulent bacterial Pst strains using Illumina SBS high-throughput sequencing. Our results indicate that miRNAs contribute to plant immunity by regulating multiple hormone-signaling pathways and some stress-responsive genes. We also compared the results from sequencing-based profiling with those from small RNA Northern blot and miRNA quantitative RT-PCR, revealing some inconsistency between these general techniques. We discussed the procedural differences between these two methods that may lead to the inconsistency.

Results

Identification of Pst-regulated miRNAs by deep sequencing

To examine the expression of Arabidopsis miRNAs in response to the infection of various strains of bacterial pathogen Pst, we performed a small-RNA profiling analysis using Illumina SBS deep sequencing. RNAs were isolated from 4-week-old short-day grown Arabidopsis Col-0 plants after 6 and 14 h post-inoculation (hpi) of mock (10 mM MgCl2), a non-pathogenic strain Pst DC3000 hrcC, a virulent strain Pst DC3000 carrying an empty vector (EV), and an avirulent strain Pst DC3000 carrying an effector avrRpt2 (eight libraries). We included biological replicates for mock 6-hpi, mock 14-hpi, Pst DC3000 hrcC 14-hpi and Pst DC3000 (avrRpt2) 14-hpi samples; we also sequenced one of the Pst DC3000 (avrRpt2) 14-hpi samples twice to ensure the reproducibility of the sequencing results (five libraries in total). Pst DC3000 hrcC induces PTI triggered by the recognition of PAMPs by host PAMP recognition receptors. Pst DC3000 (EV) induces effector–triggered susceptibility (ETS) by injecting effector proteins into the plant cell to suppress PTI (Chisholm et al. 2006; Jones and Dangl 2006). Pst DC3000 (avrRpt2) induces effector-triggered immunity (ETI) and hypersensitive response (HR) in Col-0 carrying cognate resistance gene RPS2 and leads to cell death symptom usually at 15-hpi under the condition and concentration we used. Our samples were collected at 14-hpi before HR could be visualized. Sequencing reads with recognizable adaptor sequences removed were compared with Arabidopsis nuclear, chloroplast, mitochondrial genome sequences and cDNAs. From a total of more than 24.6 million sequence reads from all libraries, we obtained 13,985,938 total reads that have high sequencing quality (i.e., with no N’s), contain a 3′ sequencing adaptor sequence, and can map perfectly to the Arabidopsis genome and cDNA sequences with zero mismatches. This large collection of deep-sequencing data has been deposited to the NCBI/GEO databases (GSE19694). The data have been recently analyzed to identify Arabidopsis miRNA precursors that can harbor more than one authentic miRNAs (Zhang et al. 2010) and miRNA precursors that produce, at the same sites, miRNAs as well as siRNAs, some of which guide DNA methylation at their target loci (Chellappan et al. 2010). In the current study, we further analyze the deep sequencing data to understand the function of miRNAs in response to the challenge of bacterial infection.

Within the small-RNA population, 21- and 24-nt small RNA species are the two most abundant among all size classes (Fig. 1a). Within the reads that can map to the known miRNAs, 80% have length of 21-nt (Fig. 1c). More than 90% of these miRNA reads start with a uracil (Fig. 1d), which is in shape contrast to the first nucleotide distribution of all small RNA reads (Fig. 1b). This result is consistent with the current understanding of miRNA composition. We also examined these distributions for individual conditions, but no obvious difference among them and the results in Fig. 1 was observed (data not shown).

Fig. 1
figure 1

Distributions of lengths and the first nucleotides of small RNAs and miRNAs from the 13 small-RNA sequencing data. a Length distribution of all sequencing reads (“all”) that have a 3′ sequencing adaptor, are of high quality and can perfectly map to the Arabidopsis genome or cDNA sequences, and distribution of those reads in “all” that exclude miscellaneous noncoding RNAs, including tRNAs, rRNAs, snRNAs and snoRNAs (“all no misc”). Horizontal axis shows the sequence lengths and the vertical axis is the percent of sequences of a particular length. b Distribution of the first nucleotides of “all” and “all no misc” reads in a, respectively. c Length distribution of sequence reads that perfectly map to the known miRNAs. c Distribution of the first nucleotides of reads mapped to the known miRNAs

Full size image

Although a detailed characterization of miRNA expression, particularly the absolute expression level of a miRNA, is useful, identification of differentially expressed (DE) miRNAs at the whole genome level in response to environmental stimuli or endogenous cues is often desirable for elucidating miRNA function in particular cellular processes. To identify DE miRNAs, we need to consider the reproducibility of sequencing data. We first sequenced one of the libraries, Pst DC3000 (avrRpt2) 14-hpi, twice and compared the numbers of sequence reads that map to miRNAs. The counts of reads in each library were normalized to the numbers of reads per million (RPM) to adjust for the differences in library sizes or sequencing depths (see “Methods”). As shown in Fig. 2a, the normalized numbers of reads for miRNAs in two technical duplicates for the Pst DC3000 (avrRpt2) 14-hpi sample were strongly correlated, with the square of the correlation coefficient exceeded 0.94 and the slope of linear regression being 0.9766. We also sequenced biological replicates under 4 conditions. The normalized result on the biological replicates for mock 14-hpi was also highly correlated, with the square of the correlation coefficient of 0.838 and the slope of linear regression of 0.809 (Fig. 2b). The results on the other three pairs of biological replicates, i.e., that for mock 6-hpi, DC3000 hrcC 14-hpi and Pst DC3000 (avrRpt2) 14-hpi, were also well correlated (data not shown), similar to that in Fig. 2b. These results showed that the read counts of miRNAs from Illumina high-throughput sequencing are highly reproducible and the read counts of a miRNA in different small-RNA libraries are good indicators of the abundance of the miRNA in the corresponding biological samples.

Fig. 2
figure 2

Correlation of miRNA abundances between a two technical duplicates for the Pst (avrRpt2) 14-hpi, and b two biological replicates for mock 14-hpi from Illumina high-throughput deep sequencing

Full size image

We mapped the adaptor-trimmed sequence reads to the known miRNA loci to detect the expressed miRNAs. In this mapping process, we only considered perfect matches and those reads whose first nucleotide to be within 3-nt upstream or down-stream of the starting position of the annotated starting position of a miRNA to compensate the imperfect cleavage activities of Dicer-like proteins (DCLs). We detected the expression of Arabidopsis miRNAs in 86 families (Supplemental Table 1).

The sequence reads mapped to miRNAs were then normalized (see “Methods”). Based on the normalized read counts of miRNAs, DE miRNAs were identified using stringent criteria combining fold change and a Chi-square test across the infection conditions of interest and the respective mock treatments. Specifically, we used strict criteria of at least 10 raw sequence reads, at least threefold change between the infected and the mock, and a P-value for a Chi-square test no greater than 0.05 (see “Methods”). We identified a substantial number of DE miRNAs responding to pathogen infections. The DE miRNAs for the three infection conditions at the two time points and their fold changes are shown in Fig. 3. In particular, we identified 15, 27 and 20 miRNA families to be differentially expressed upon Pst DC3000 hrcC, Pst DC3000 (EV) and Pst DC3000 (avrRpt2) infections at either 6- or 14-hpi, respectively (Fig. 3). Overall, Pst DC3000 hrcC had the least effect on miRNA expression, whereas Pst DC3000 (EV) affected the expression of the most miRNAs among the three strains we studied (Fig. 3). Among all the miRNAs, miR158 was highly induced by all three strains at 14-hpi, and by Pst DC3000 (EV) and Pst DC3000 (avrRpt2) at 6-hpi (Fig. 3). miR158 targets several genes encoding PENTATRICOPEPTIDE REPEAT (PPR) domain-containing proteins, as well as FUCOSYLTRANSFERASE genes (AtFUT1 and AtFUT2) that encode glycosyltransferases for cell wall xyloglucan biosynthesis. PPR-domain containing proteins are mostly mitochondria or chloroplasts localized putative RNA binding proteins involved in RNA processing, metabolism, editing or translation. Their function in plant disease resistance remains to be explored. An untypical PPR-domain-containing protein gene PPRL is a target of an endogenous nat-siRNA and negatively regulates Pst DC3000 (avrRpt2)-triggered ETI (Katiyar-Agarwal et al. 2006). Interestingly, miR403, which targets RNAi effector genes ARGONAUTE 2 (AGO2) and AGO3, was the only miRNA that was suppressed by all the three strains at 14 hpi, and weakly by Pst DC3000 (avrRpt2) and Pst DC3000 hrcC at 6 hpi (Fig. 3), this result suggests the potential role of AGO2 and AGO3 in plant immunity.

Fig. 3
figure 3

Differentially expressed miRNAs identified in three bacterial infections—a type III secretion system mutated strain Pst DC3000 hrcC, b virulent strain Pst DC3000 carrying an empty vector (EV), and c avirulent strain Pst DC3000 carrying effector avrRpt2, at 6-hpi and 14-hpi with respect to the respective mock treatments. The bars shown in the figures are log (base 2) ratios of fold changes relative to the mock. miRNAs are sorted based on the fold changes at 14-hpi. FC stands for fold change

Full size image

Notably, several miRNAs that regulate plant hormone signaling and biosynthesis were differentially expressed under bacterial infection. Firstly, the expression of a group of miRNAs involved in auxin signaling pathways was largely altered. Both miR160 and miR167 that target Auxin Response Factors (ARFs) were both induced by all of the three Pst strains at 6-hpi (Fig. 3). miR390 (miR390a), which triggers the production of TAS3 trans-acting siRNAs (ta-siRNAs) that also regulate another subset of ARFs, mainly ARF3 and ARF4, was down-regulated at 6-hpi of virulent strain Pst DC3000 (EV) (Fig. 3b). More than 30% of the 23 ARFs in Arabidopsis could be regulated by miR160, miR167 or miR390/TAS3 tasiRNAs, and function as positive or negative regulators of auxin signaling pathways (Guilfoyle and Hagen 2007). In addition, miR393, which targets auxin receptors, TIR1, AFB2 and AFB3, was also differentially regulated by Pst DC3000 (EV), Pst DC3000 hrcC and Pst DC3000 (avrRpt2) at 6- or 14-hpi, respectively (Fig. 3). Secondly, miR159 was down-regulated by Pst DC3000 (EV) and Pst DC3000 (avrRpt2) at 6-hpi, but up-regulated by Pst DC3000 (avrRpt2) at 14-hpi (Fig. 3b, c). miR159 targets transcription factors MYB33, MYB65 and MYC101, the homologous genes of the barley GAMYB that activates Gibberellin (GA)-signaling pathways (Reyes and Chua 2007; Millar and Gubler 2005). MYB33 and MYB101 act as positive regulators of ABA signaling pathways in Arabidopsis (Reyes and Chua 2007).

Furthermore, stress-related miRNAs were also differentially regulated by bacteria-challenges. miR398, which targets Copper/Zinc superoxide dismutase genes (CSDs), was repressed at 14-hpi of Pst DC3000 hrcC and Pst DC3000 (avrRpt2) (Fig. 3a, c), which is consistent with an early result (Jagadeeswaran et al. 2009). Pst DC3000 hrcC triggers PTI and Pst DC3000 (avrRpt2) triggers ETI (Chisholm et al. 2006; Jones and Dangl 2006). Both immunity responses trigger oxidative burst resulting in the rapid accumulation of reactive oxygen species (ROS) (Lamb and Dixon 1997). Down-regulation of miR398 under the same conditions will lead to elevated level of CSD that help detoxify increased ROS. miR398 is also suppressed by flg22 (Li et al. 2010). Overexpression of miR398b negatively regulates callose deposition and exhibits enhanced susceptibility to virulent and non-pathogenic strains of Pst (Li et al. 2010). miR408, which was predicted to target copper protein plantacyanin, laccase copper protein and copper ion binding protein genes, was first induced at 6-hpi (Fig. 3) and then down-regulated at 14-hpi by Pst DC3000 hrcC and Pst DC3000 (avrRpt2) (Fig. 3a, c). miR398 and miR408 are also differentially regulated by abiotic stresses (Trindade et al. 2010; Abdel-Ghany and Pilon 2008), suggesting their general roles in stress responses.

Expression validation using small RNA Northern blot analysis

To confirm the results we obtained from small-RNA deep sequencing, we performed a small-RNA Northern blot analysis to further validate several miRNAs that were detected as differentially expressed by deep sequencing. We found that the majority of small RNA Northern blot results were consistent with the deep sequencing data (Fig. 4). For example, at 6-hpi miR160 was induced in all three treatments as determined by both methods (Figs. 3, 4a). miR408 was induced at 6-hpi in Pst DC3000 hrcC and Pst DC3000 (avrRpt2) as detected by both methods (Figs. 3a, c, 4a). Up-regulation of miR822 in Pst DC3000 (avrRpt2) and Pst DC3000 (EV) treatments observed in Northern blot (Fig. 4a) was consistent with the result from deep sequencing (Fig. 3b, c). At 14-hpi, miR158 was induced by all the three strains and the inductions shown in the data from the deep sequencing and Northern blot matched each other well. Similarly, miR159 and miR393 were induced whereas miR408 and miR822 were suppressed by Pst DC3000 (avrRpt2) (Fig. 4b), which were also consistent with the results from the deep sequencing data (Fig. 3c).

Fig. 4
figure 4

Small RNA Northern blot verification of some of the differentially expressed miRNAs detected by deep sequencing from the three bacterial pathogen infections and the mock controls at a 6-hpi and b 14-hpi. Similar results were obtained from two biological duplicates. U6 RNA was used as a control for measuring the relative amount of the bands (shown below each panel). Imagequant software version 2.1 was used for relative amount quantification (GE Healthcare Life Sciences)

Full size image

On the other hand, we also observed some discrepancies between the levels of miRNA expression measured by high-throughput sequencing and Northern blot (Figs. 3, 4). For example, in the deep sequencing results of 6-hpi, miR159 was repressed by Pst DC3000 (EV) and Pst DC3000 (avrRpt2) (Fig. 3b, c), miR393 was down-regulated by all three infections (Fig. 3), and miR162 was induced by Pst DC3000 (EV) (Fig. 3b). On the contrary, none of these miRNAs showed clear change under any of these conditions in the Northern blot analysis (Fig. 4). Similarly, our deep sequencing data showed that miR398 was induced by Pst DC3000 (EV) at 6-dpi (Fig. 3b), however, Northern blot analysis by Jagadeeswaran et al. did not detect any clear difference of miR398 expression under the same condition (Jagadeeswaran et al. 2009). On the other hand, Northern blot analysis also detected some expression changes that were not detected by deep sequencing. For example, our small RNA Northern blot analysis revealed that miR319 and miR852 were clearly induced by Pst DC3000 hrcC and Pst DC3000 (avrRpt2) at 14-hpi (Fig. 4b), while the deep sequencing profiling failed to detect these differential expressions (Fig. 3). To further examine this observation, we selected three miRNAs (miR159, miR408 and miR319) and quantified them by miRNA real time quantitative RT-PCR (QRT–PCR) (Fig. 5). In consistent with the result from Northern blot and deep sequencing results, QRT–PCR results showed that miR159 was up-regulated and miR408 was down-regulated in Pst DC3000 hrcC, Pst DC3000 (EV) and Pst DC3000 (avrRpt2) as compared with mock treatment at 14-hpi (Fig. 5a, b). miR319 showed an overall up-regulation in the treatments of Pst DC3000 hrcC, Pst DC3000 (avrRpt2) and Pst DC3000 (EV) as compared with the mock (Fig. 5c), which is consistent with the Northern blot analysis (Fig. 4), although it detected more induction of miR319 by Pst DC3000 hrcC and less induction of miR319 by Pst DC3000 (avrRpt2) than that in Northern blot analysis. The possible causes of discrepancies among these three techniques are included in “Discussion” section.

Fig. 5
figure 5

Relative expression levels of miR159, miR408 and miR319 in Pst DC3000 (EV), Pst DC3000 hrcC and Pst DC3000 (avrRpt2) compared to the mock by quantitative RT-PCR. The comparative threshold cycle (Ct) method was used to determine relative transcript levels (Bulletin 5279, Real-time PCR applications Guide, Bio-Rad), using actin as an internal control. The expression level of a mock was set to 1, and the values for the other conditions are fold changes for those conditions. a miR159, b miR408, c miR319. Similar results were obtained from two biological replicates

Full size image

Expression analysis of the targets of bacteria-regulated miRNAs

To determine whether the expression change of miRNAs in response to Pst infection altered the expression of their targets, we analyzed the public microarray dataset from the AtGenExpress Consortium (ME00331) collected from the same treatments as our small RNA libraries (see “Methods”), i.e., Pst DC3000 (EV)- and Pst DC3000 hrcC-infected plants at 6-hpi. We used miRNA targets that were previously validated to analyze the functions of some of the miRNAs that had a significant expression change between a treatment and the mock (see “Methods”). We focused on those miRNAs that had at least threefolds of expression change in an infection relative to the mock (Fig. 3) and their protein-coding targets that exhibited an anti-correlated expression with respect to that of the targeting miRNAs. The results are listed separately in Supplemental Table 2. The expression of 12 mRNA target genes, including miR160 and miR167 targets ARF8, ARF10, ARF16, and ARF17, miR393 targets TIR1, AFB2 and AFB3, and miR159 targets MYB33 and MYB65, was negatively correlated with their corresponding miRNAs under the infection of Pst DC3000 hrcC at 6-hpi (Supplemental Table 2a). Under the infection of Pst DC3000 (EV) at 6-hpi, 28 mRNA target genes had expression anti-correlated to the expression of their targeting miRNAs (Supplemental Table 2b). Again, miR160, miR167, miR393 that involved in auxin signaling pathway and miR159 that involved in ABA signaling pathway were on the list.

Although many miRNAs and targets exhibited similar anti-correlated expression patterns in both Pst DC3000 hrcC and Pst DC3000 (EV), a substantial number of miRNA/target pairs were infection specific. Specifically, miR166/ATHB-15, miR167/ARF8, miR393/AFB2, and miR393/At3G23690 (a bHLH protein) were specific to Pst DC3000 hrcC at 6-hpi, whereas ten miRNAs (i.e., miR158, miR161, miR162, miR167, miR171, miR394, miR396, miR402, miR775, and miR847) and their twenty targets listed in Supplemental Table 2b were anti-correlated specifically in Pst DC3000 (EV). These results suggest that there was specific regulation mediated by miRNAs in the non-pathogenic, compatible and incompatible interactions between bacteria and host, although largely overlapped patterns of miRNA expression were observed. In particular, miR162, which targets DCL1, and miR775, which targets a galactosyltransferase family protein, may play a rather specific role in compatible interaction. Future in-depth functional analysis on these miRNAs and their targets in plant defense responses may help address this question. It is already known that bacteria pathogens can suppress miRNA pathway by injecting certain effector proteins into the plant cell (Navarro et al. 2008). It is likely that the virulent strain Pst DC3000 (EV) utilizes an additional approach to suppress miRNA pathway by inducing miR162 and subsequently down-regulating DCL1, and thus promotes disease. Furthermore, not all the target genes were negatively regulated by corresponding DE miRNAs in a given condition. For example, miR167 was negatively correlated with 3 targets, including ARF8, under the infection of Pst DC3000 hrcC at 6-hpi (Supplemental Table 2a), but regulated ARF6 instead in response to Pst DC3000 (EV) (Supplemental Table 2b), suggesting that DE miRNAs may selectively regulate their targets under certain conditions and this type of selective regulation is prevalent. Indeed, this phenomenon of selective regulation of miRNAs was observed previously in other studies, e.g., that in (Rubio-Somoza et al. 2009; Jagadeeswaran et al. 2009). For example, miR398 was suppressed by Pst DC3000 (avrRpt2) or Pst (avrRpm1). Only one of the two targets, CSD1, was derepressed, but the other target, CSD2, was slightly repressed (Jagadeeswaran et al. 2009). Both CSD1 and CSD2 were shown to be regulated by miR398 under abiotic stresses (Sunkar et al. 2006). These results illustrated that a single miRNA can respond to both biotic and abiotic stresses and selectively regulate different targets under different conditions.

Furthermore, we quantified target mRNA expression levels by real-time QRT–PCR for the 14-hpi treatments on selected differentially regulated miRNAs from Northern and miRNA QRT–PCR (Fig. 6). The results revealed that miR158 target AtFUT1 (At2g03220) and miR166 target REVOLUTA (At5g60690) were down-regulated in both Pst DC3000 (EV) and Pst DC3000 (avrRpt2) treated samples at 14-hpi (Fig. 6a, c). MYB33 (At5g06100), a miR159 target, was down regulated in all three treatments at both 6- and at 14-hpi (Fig. 6b). PLANTACYANIN (ARPN, At2g02850), a target of miR408, was induced by all the three conditions at 14-hpi (Fig. 6e). These data showed a clear negative correlation between the expression levels of the mRNA targets and their corresponding miRNAs. In addition, TCP4, the target of miR319 was strongly down-regulated by Pst DC3000 (avrRpt2), but weakly repressed by Pst DC3000 (EV) (Fig. 6d). This result is consistent with the results from miRNA QRT–PCR and Northern blot where miR319 was induced by Pst DC3000 (avrRpt2). In summary, these results provided insightful information of miRNA-mediated gene regulation in bacterial pathogen infection.

Fig. 6
figure 6

Expression analysis of miRNA target genes using quantitative real time RT-PCR. Relative expression levels of a At2g03220 (miR158 target), b At5g06100 (miR159 target), c At5g60690 (miR166 target) and d At3g15030 (miR319 target), e At2g02850 (miR408 target) at 14-hpi of Pst DC3000 hrcC, Pst DC3000 (EV), and Pst DC3000 (avrRpt2) infections as compared with the mock treatment. The comparative threshold cycle (Ct) method was used (Bulletin 5279, Real-time PCR applications Guide, Bio-Rad), with actin as an internal control. Similar results were obtained from two biological replicates

Full size image

Discussion

miRNAs orchestrate hormone signaling pathways for defense response in plant

We identified a group of miRNAs that were differentially expressed after the infection of bacterial pathogen Pst using small RNA profiling by Illumina SBS high-throughput deep sequencing. Notably, many of these bacteria-regulated miRNAs manipulate the components of plant hormone signaling pathways, including auxin, ABA and jasmonic acid (JA) signaling pathways. Among them, miR160, miR167, miR390 and miR393 all regulate genes involved in the auxin signaling pathway, including different ARFs and auxin receptors. Auxin is an essential phytohormone for growth and development. Many biotrophic pathogens could synthesize auxin or auxin-like molecules to promote virulence (Navarro et al. 2006; Chen et al. 2007; Wang et al. 2007). Pst DC3000 (avrRpt2) (when expressed in a rps2 mutant background) promotes auxin signaling/biosynthesis (Chen et al. 2007). As a result, host plants have developed several counter measures, including miRNA-mediated gene regulation, to suppress auxin signaling and subsequently inhibit pathogen growth (Navarro et al. 2006; Wang et al. 2007). Our data suggested that miR160, miR167, miR390 and miR393 regulated the expression of a group of signaling genes that were involved in various steps of auxin signaling, and thus contributed to the inhibition of pathogen growth. miR393 was shown to be induced by bacterial elicitor flg22 to repress the expression of auxin receptor genes (Navarro et al. 2006). Pst DC3000 hrcC could induce miR160, miR167 and miR393 at 3-hpi (Fahlgren et al. 2007). We observed similar induction patterns of miR160 and miR167 at 6-hpi in our study. However, we did not observe clear changes of miR393 expression at 6- and 14-hpi after Pst treatment in our Northern blot analysis, suggesting that miR393 may regulate auxin signaling at an early stage of bacterial infection.

Our small RNA Northern blot analysis revealed that miR319 was induced by Pst DC3000 hrcC and Pst DC3000 (avrRpt2) at 14-hpi (Fig. 4). miR319 targets TCP (TEOSINTE BRANCHED/CYCLOIDEA/PCF) transcription factor family genes which directly regulate LIPOXYGENASE2 (LOX2). LOX2 encodes a chloroplast-localized enzyme that is responsible for the first step in the JA biosynthesis pathway (Schommer et al. 2008). Because JA signaling is usually antagonistic to salicylic acid (SA) signaling, while SA signaling is important for plant defense against biotrophic pathogens, including Pst. Down-regulation of JA biosynthesis by miR319 through TCP genes in response to Pst DC3000 hrcC and Pst DC3000 (avrRpt2) would promote SA-mediated plant resistance responses.

Furthermore, miR159 was induced by Pst DC3000 (avrRpt2) at 14-hpi (Figs. 3c, 4). miR159 targets transcription factors MYB33 and MYB101, both of which are positive regulators of ABA signaling pathways (Reyes and Chua 2007). Although ABA plays a positive role in antibacterial defense at early infection stage by regulating the pre-invasive stomata-based responses (Melotto et al. 2006). ABA was shown to attenuate SA-mediated resistance at later infection stages and can also suppress callose deposition in response to flg22 (de Torres-Zabala et al. 2007). Bacterial effectors, such as AvrPtoB, promote ABA biosynthesis to suppress SA signaling pathway (de Torres-Zabala et al. 2007). Induction of miR159 suppresses ABA signaling via MYB33 and MYB101, which activates SA signaling pathway and subsequently promotes SA-mediated defense responses.

It is worth noting that, in many cases, the avirulent strain Pst DC3000 (avrRpt2) induces miRNAs to a higher level than by non-pathogenic Pst DC3000 (hrcC) and virulent strains Pst DC3000 (EV). Figure 4 shows several such examples, including miR158, miR159, miR166 and miR319. miR159 and miR319 suppress components involved in ABA and JA signaling pathways and subsequently promote SA-mediated defense signaling pathway. These results suggest that plants have more miRNAs involved in multiple hormone signaling pathways to trigger stronger and faster defense responses in ETI than in non-pathogenic strain and virulent strain-induced PTI. Similar observations were made in early mRNA profiling studies (Tao et al. 2003).

In summary, almost all the plant hormones, including SA, JA, ethylene, auxin, ABA, Gibberellin and cytokinin, have been implicated in playing a direct or indirect role in plant defenses (Grant and Jones 2009; Bari and Jones 2009; Spoel et al. 2007). Bacterial pathogens often produce hormones or hormone like molecules, such as auxin or coronatine, to interfere with host hormone homeostasis and signaling, which leads to the suppression of plant immune responses. Plants have developed successful defense mechanisms against bacterial-triggered hormone perturbation. Our results suggested that a group of miRNAs work coordinately to promote defense responses by regulating and fine-tuning multiple host hormone signaling pathways in response to pathogen invasion.

Discrepancy between deep sequencing and Northern blot for miRNA profiling

Comparing with hybridization-based profiling methods, such as microarray, Northern blot, and PCR, high-throughput sequencing-based profiling methods have at least two advantages. First, they can not only profile the known genes and known small RNAs, but also detect and characterize novel transcripts, while hybridization-based methods can only provide information of the annotated genes. This advantage is important for identification of novel transcripts that are expressed under rare conditions, such as under a peculiar pathogen infection, or at low abundance. Second, sequencing-based methods can profile transcripts of single copy, and more importantly have a resolution of single nucleotide to detect even single nucleotide polymorphisms. This feature is essential for distinguishing individual members of a miRNA family that often differ by one or two nucleotides. In contrast, hybridization-based methods, such as microarray and Northern blot, could often detect a combined signal of individual members of the same miRNA families, due to the possible cross-hybridization that may not be able to discriminate a small number of mismatches between a probe and a sequence annealed. Thus hybridization-based methods measure the abundance of the intended miRNA families, not necessarily individual miRNAs.

The discrepancy between the results from Illumina deep sequencing and Northern blot and miRNA QRT–PCR as observed in the current study is perhaps due to the differences between the procedures of these two profiling techniques. A deep-sequencing method often uses small RNA adaptor ligation that may introduce a bias that some small RNA species are more likely to be ligated than others (Chellappan and Jin 2009; Katiyar-Agarwal and Jin 2007). Some miRNA modifications may also reduce the efficiency of ligation or even block ligation reactions, whereas small-RNA Northern blot is based on nucleotide hybridization, which could reveal all small RNAs that are reverse complementary to the probes for hybridization regardless the modification. However, Northern blot will not be able to measure the exact size and sequence distribution at a locus. Furthermore, PCR amplification is the last step of the Illumina small RNA library preparation before sequencing. PCR is an inherently biased process because it tends to amplify highly abundant sequences more effectively than low abundant sequences. Therefore, PCR may exaggerate the differences between the high abundant and low abundant sequences. This discrepancy between the current sequencing-based approaches and hybridization-based approaches has been observed in a recent comparison of microarray-based miRNA profiling methods and Illumina deep sequencing (Git et al. 2010). The observation we made in the current study is also in agreement with the result from a comparison of two deep sequencing protocols for Illumina deep sequencing platform, one used RNA amplification and the other did not (Kapranov et al. 2010). The result in (Kapranov et al. 2010) showed that, although the amplification-based protocol in general preserves the general trend of abundances of small RNAs, the normalized abundances of individual small RNAs can vary up to 1–2 orders of magnitude in the amplification-based protocol versus the non-amplification-based protocol, revealing the biases introduced in amplification.

In short, the sequencing-based and hybridization-based small-RNA profiling techniques have their own advantages and drawbacks. They may provide different measures of the same transcript in the same biological sample. Likewise, two miRNAs (or siRNAs) with similar numbers of sequencing reads may in fact differ substantially in their absolute abundances in a sample. Therefore, we must exercise extra caution when directly compare small RNA expression information using deep sequencing data and PCR-based assays or Northern blot hybridization.

Methods

Plant materials, small-RNA library construction and deep sequencing

Arabidopsis Col-0 was used as the plant material in this study. Profiling experiments were performed on 4-week-old Arabidopsis plants grown at 23°C with 12-h light. Infiltration of the following bacteria Pst DC3000 EV, Pst DC3000 hrcC, Pst DC3000 (avrRpt2) and mock (10 mM MgCl2) was carried out as described previously (Katiyar-Agarwal et al. 2006). The small RNA libraries for deep sequencing were constructed by 5′phosphate-dependent method as described (Chellappan and Jin 2009), and the libraries were sequenced at the Illumina Inc. The raw and processed sequencing reads are available at NCBI/GEO under the accession number GSE19694.

Small-RNA Northern blot analysis

For Northern blot analysis 80 μg total RNA was resolved on 15% denaturing polyacrylamide gel. We used the following oligo probes to detect miRNAs. Blots were exposed to phosphorscreens and scanned using PhosphorImager (Molecular dynamics).

  • miR160-anti: tggcatacagggagccaggca

  • miR148-anti: gccagggaagaggcagtgcat

  • miR822-anti: catgtgcaaatgcttcccgca

  • miR159-anti: tagagctcccttcaatccaaa

  • miR162-anti: ctggatgcagaggtttatcga

  • miR393-anti: ggatcaatgcgatccctttgga

  • miR158-anti: tgctttgtctacatttggga

  • miR166-anti: ggggaatgaagcctggtccga

  • miR408-anti: gccagggaagaggcagtgcat

  • miR852-anti: cagaactaaggcgcttatctt

Quantitative RT-PCR to quantify miRNA expression levels

The expression levels of miRNAs were quantified by a modified quantitative RT-PCR from described (Varkonyi-Gasic et al. 2007). About 2 μg of DNase I (Invitrogen) treated total RNA was used for reverse transcription with artificially designed stem-loop specific RT-primers and oligo dT primer separately using Superscript III (Invitrogen). To capture individual miRNAs, specific stem-loop RT primers that were complementary to each miRNA by a six-nucleotide extension at the 3′ end was used. The specificity of stem-loop RT primers to individual miRNAs is conferred by a six-nucleotide extension at the 3′ end. The reaction was incubated at 16°C for 30 min followed by 30 s at 30°C and 42°C and 1 s at 50°C of 60 cycles and finally 5 min at 85°C. Following RT, quantitative PCR in real time was performed using SensiMix SYBR (Quantace) using miRNA specific forward primer and universal reverse primers as listed below. The comparative threshold cycle (Ct) method was used for determining relative transcript levels (Bulletin 5279, Real-time PCR applications Guide, Bio-Rad). Actin was used as an internal control.

  • QmiR159:RT-primer

  • GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTAGAGC

  • miR159Fp:CGGCGGTTTGGATTGAAGGGA

  • qmiR319:RT-primer.

  • GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGGGAGC

  • qmiR319Fp:TTGCGGTTGGACTGAAGGGAG

  • qmiR408:RT-primer.

  • GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGCCAGG

  • qmiR408Fp:TTGCGGATGCACTGCCTCTTC

  • qmiRUniversal Rp: GTGCAGGGTCCGAGGT

  • ActinF: AGTGGTCGTACAACCGGTATT

  • ActinR: GATGGCATGAGGAAGAGAGAA

Quantitiative RT-PCR to quantify miRNA targets expression levels

For QRT-PCR analysis, 5 μg of total RNA was used for synthesizing cDNA. DNA contamination was removed by using DNase I (Invitrogen). Amplification of miRNA targets was carried out using real-time PCR machine (iQ5, Bio-Rad). The following oligos were used to perform the quantitative RT-PCR analysis.

  • miR158-ATFUT1(2g03220)QRTF: CCAAATGTTACTAACATCAATTC

  • miR158-ATFUT1(2g03220)QRTR: CGGACCACATCGCTTGTGAAGCTTT

  • miR159-MYB33(5g06100)QRTF: TAGATGGGCACGTATGGCTGCACATTTG

  • miR159-MYB33(5g06100)QRTR: TACTCTTGACTCCACTCAAGTGCCT

  • miR166-REVOLUTA(5g60690)QRTF: GGATTGCTCTCAATCGCAGAGG

  • miR166-REVOLUTA(5g60690)QRTR: CTCACAAACTGAGAAGCTGAAGC

  • miR319-TCP4(3g15030)QRTF:ACCATCACCTTCCTCCACCGGTTCACC

  • miR319-TCP4(3g15030)QRTR: AAATAGAGGAAGCAGAGGACGGCTTGTG

  • miR408-ARPN(2g02850)QRTF: TAGAGCCGGCGACGTTCTTGTGTT

  • miR408-ARPN(2g02850)QRTR: AAGTTATACGATCTTTGCCTGAAGTGTA

Preprocessing and normalization of sequence reads

The raw sequence reads were processed to remove reads with no 3′ sequencing adaptors, of low quality (having at least an N), or shorter than 17-nt. The adaptor trimming was done by an in-house program that recursively searches for the longest subsequence of the adaptor appearing within a sequencing read. If a raw sequence read did not have a subsequence of the adaptor longer than 5-nt, it was considered as carrying no adaptor. The adaptor-trimmed reads that are longer than 16nt and can map to the Arabidopsis genome (TAIR version 9.0, from http://www.arabidopsis.org) and cDNAs (ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR9_genome_release/TAIR9_sequences/TAIR9_cdna_20090619) with zero mismatches were retained for further analysis. The mapping of sequence reads to the genome and cDNAs were done using Bowtie (http://bowtie-bio.sourceforge.net/index.shtml).

To facilitate identification of differentially expressed miRNAs, we normalized the numbers of sequencing reads in a small RNA library. For normalization, we assumed the total number of small RNAs in a cell to be a constant. To obtain the total number of reads of all small RNAs produced in the cell, we identified reads that can map to the Arabidopsis nuclear, chloroplast and mitochondrial genomes as well as the transcripts in available cDNA libraries since some small RNAs can originate from mRNA splicing junctions; we called this the number of total mapped reads. The normalized abundance of a unique small RNA is the ratio of the number of sequencing reads for the small RNA and the number of total mapped reads, multiplied by one million.

Identification of differentially expressed miRNAs

To identify differentially expressed miRNAs across bacterial infected and mock treated Arabidopsis plants, we first mapped the pre-processed sequencing reads to the genomic loci of the known miRNAs (using Bowtie) based on miRBase annotation (version 15, http://www.mirbase.org). We only considered miRNAs that had at least 10 raw sequencing reads in an infection and the mock conditions of interest in order to deal with noise in sequencing. In the mapping, we allowed up to 3-nt shift upstream and down-stream from the annotated starting locus of a miRNA to compensate for possible variation in Dicer activities. We compared DC3000 EV 6-hpi, DC3000 hrcC 6-hpi and Pst DC3000 (avrRpt2) 6-hpi samples with the mock 6-hpi plant, and DC3000 EV 14-hpi, DC3000 hrcC 14-hpi and Pst DC3000 (avrRpt2) 14-hpi samples with the mock 14-hpi. Using normalized sequence reads, candidate differentially expressed (DE) miRNAs in an infection condition were chosen as those whose normalized abundance increased or decreased by at least threefolds with respect to the corresponding mock treatments. These candidate DE miRNAs were further tested using the Chi-square test with a threshold of 0.05. Only those miRNAs that passed the test were taken as DE genes. In the statistical test, the null hypothesis was that no difference for a miRNA between an infection and the mock treated conditions. The test involved comparing the observed frequency (number of reads) of a miRNA under the infection against its expected frequency in the mock condition.

Microarray data and processing, and miRNA target genes

Genome-wide mRNA expression profiling data from the international collaborative AtGenExpress project, available at TAIR (http://www.arabidopsis.org), were used in our study. Specifically, we used the microarray dataset for Pst DC3000 hrcC and Pst DC3000 (EV) collected at 6-hpi (TAIR Accession number ME00331), which perfectly matched two of the experimental conditions of our study. These datasets and the dataset of its corresponding control from the Affymetrix microarray platform were normalized using gcRMA in BioConductor (http://www.bioconductor.org).

The miRNA target genes we used were from those that were identified by Parallel Analysis of RNA Ends (PARE) and reported in the literature, which were collected by the Meyers lab at (PARE, http://mpss.udel.edu/common/web/targets.php?SITE=at_pare) and (http://mpss.udel.edu/at/target.php).