Kira C. M. Neller1
Katalin A. Hudak- 1Department of Biology, York University, Toronto, ON, Canada
- 2Department of Electrical Engineering and Computer Science, York University, Toronto, ON, Canada
The American pokeweed plant, Phytolacca americana, displays broad-spectrum resistance to plant viruses and is a heavy metal hyperaccumulator. However, little is known about the regulation of biotic and abiotic stress responses in this non-model plant. To investigate the control of miRNAs in gene expression, we sequenced the small RNA transcriptome of pokeweed treated with jasmonic acid (JA), a hormone that mediates pathogen defense and stress tolerance. We predicted 145 miRNAs responsive to JA, most of which were unique to pokeweed. These miRNAs were low in abundance and condition-specific, with discrete expression change. Integration of paired mRNA-Seq expression data enabled us to identify correlated, novel JA-responsive targets that mediate hormone biosynthesis, signal transduction, and pathogen defense. The expression of approximately half the pairs was positively correlated, an uncommon finding that we functionally validated by mRNA cleavage. Importantly, we report that a pokeweed-specific miRNA targets the transcript of OPR3, novel evidence that a miRNA regulates a JA biosynthesis enzyme. This first large-scale small RNA study of a Phytolaccaceae family member shows that miRNA-mediated control is a significant component of the JA response, associated with widespread changes in expression of genes required for stress adaptation.
Introduction
The American pokeweed, Phytolacca americana, is a non-model plant with promising applications in agriculture. Pokeweed synthesizes pokeweed antiviral protein (PAP), a ribosome inactivating protein with RNA N-glycosidase activity (Endo et al., 1988). Several mutants of PAP have been generated and expressed heterologously in transgenic plants to impart novel antiviral and antifungal properties (Zoubenko et al., 1997, 2000; Wang et al., 1998; Dai et al., 2003). Additionally, pokeweed is a heavy metal hyperaccumulator with a potential role in phytoremediation (Peng et al., 2008; Liu et al., 2010; Zhao et al., 2011).
We recently assembled and annotated the pokeweed mRNA transcriptome and reported its regulation by jasmonic acid (JA; Neller et al., 2016). Jasmonates are well-characterized plant signal molecules of the oxylipin family that are involved in abiotic and biotic stress responses, as well as growth and development. JA mediates tolerance to drought, salt, heat and cold stresses, resistance to pathogens and insects, and induces embryogenesis, flowering, and senescence (reviewed in Wasternack and Feussner, 2017). JA is also implicated in tripartite interactions between insect vectors, viruses, and plants (Sun et al., 2017). Interestingly, JA treatment reduces heavy metal-induced oxidative stress (Singh and Shah, 2014) and conversely, heavy metal exposure is associated with biosynthesis of a JA precursor (Foroughi et al., 2014). Given that pokeweed exhibits broad-spectrum pathogen resistance, is a heavy metal hyperaccumulator, and undergoes marked transcriptional reprograming upon JA treatment, we are interested in exploring the range and regulation of defense strategies in this non-model plant.
JA biosynthesis is described by the octadecanoid pathway (reviewed in Schaller et al., 2000). Briefly, alpha-linolenic acid released from plastid membranes is oxidized by a13-lipoxygenase (LOX), followed by coupled dehydration-cyclization by allene oxide synthase (AOS) and allene oxide cyclase (AOC) to generate 12-oxophytodienoic acid (OPDA). This undergoes reduction by OPDA reductase 3 (OPR3) and is followed by three rounds of beta oxidation to form JA. The JA-amido synthetase JAR1 produces the bioactive form of JA, (+)-7-iso-JAIle (JA-Ile); this is perceived by the COI1-JAZ co-receptor and results in ubiquitin-mediated degradation of JAZ repressors and activation of several transcription factors (reviewed in Chini et al., 2016).
Here, we complement our previous transcriptome analysis with a paired RNA-Seq study investigating miRNA-associated regulation of the JA response. miRNAs are non-coding RNAs of ~20–24 nucleotides (nt) that regulate post-transcriptional gene expression (reviewed in Yu et al., 2017). miRNA biogenesis begins with transcription of the MIR gene by RNA polymerase II to yield a primary miRNA (pri-miRNA), containing a stem-loop region flanked by unstructured arms. In sequential steps, dicer-like 1 RNase (DCL1) excises the stem-loop to form the pre-miRNA, then produces a smaller duplex comprising the miRNA and its opposing strand, classically termed miRNA*. In the cytoplasm, the miRNA associates with argonaute 1 (AGO1) to establish the RNA-induced silencing complex (RISC). Plant miRNAs act primarily through target mRNA cleavage, owing to high miRNA/target sequence complementarity, although translational inhibition has also been observed. miRNAs have been implicated in numerous abiotic and biotic stress responses, such as those induced by temperature, salinity, drought, light, nutrient deficiencies, oxidative stress, and pathogens (reviewed in Khraiwesh et al., 2012; Li et al., 2017; Islam et al., 2018). Previous reports have characterized the effect of JA or its methyl ester (JA-Me) on miRNAs in Arabidopsis, Nicotiana attenuata, Taxus chinensis, Lycoris aurea, and Catharanthus roseus (Qiu et al., 2009; Bozorov et al., 2012; Zhang et al., 2012; Xu et al., 2016; Shen et al., 2017). Control of JA biosynthesis through miR319-targeted TCP transcription factors has also been established (Schommer et al., 2008; Zhao et al., 2015). However, despite the importance of JA in plant development and stress resistance, the expression of JA-responsive miRNAs and targets has not been correlated on a genome-wide scale.
Although a Phytolaccaceae family reference genome is unavailable, current bioinformatic tools support miRNA analysis based on an assembled transcriptome (Dai and Zhao, 2011; Yang and Li, 2011). Here, we identify JA-responsive miRNAs, their potential roles, and construct a miRNA/target interaction network that incorporates correlated expression data from biological replicates. Our results indicate that many defense-related genes are regulated by JA-induced miRNAs, including the JA biosynthesis gene OPR3. These novel targets have potential agricultural applications to improve resistance to pathogen and environmental stresses.
Materials and Methods
Plant Growth Conditions and Jasmonic Acid Treatment
Pokeweed plants at the 4-leaf stage of growth were sprayed with 5 mM JA suspended in 0.5% ethanol, or 0.5% ethanol alone (control), and leaf tissue was harvested 24 h post treatment (Neller et al., 2016). Through time-course analysis, we established previously that PAP is expressed maximally 24 h following JA treatment (Klenov et al., 2016). This time point was chosen for both the mRNA transcriptome-wide analysis of pokeweed (Neller et al., 2016) and the current paired miRNA analysis, as we were interested in identifying defense genes co-expressed with PAP.
Small RNA Sequencing and Data Processing
Total RNA was extracted from leaf tissue of pokeweed plants as reported in Neller et al. (2016) and the small RNA fraction was sequenced according to Klenov et al. (2016). Briefly, the 15–30 nt size RNA portion was gel purified and ligated to directional 5p and 3p primers. Following cDNA synthesis and PCR-amplification, the small RNA libraries were sequenced on a SOLiD 5500 XL machine. Six libraries were sequenced, with three biological replicates per treatment. Each replicate consisted of an equal amount of small RNA pooled from three independent plants. The raw sequencing reads were processed as described in Klenov et al. (2016). Non-coding RNA (rRNA, tRNA, snRNA, and snoRNA) and mature miRNA sequences were annotated by mapping reads to the sugar beet genome (Dohm et al., 2014) and the plant miRNA database (Zhang et al., 2010; http://bioinformatics.cau.edu.cn/PMRD/). Mapping was performed with Bowtie software (v. 1.0.1; Langmead et al., 2009) and up to three mismatches were allowed in alignments.
miRNA Prediction and Annotation
miRNA prediction was performed with miRDeep-P (v. 1.3; Yang and Li, 2011), a modified version of the miRDeep algorithm (Friedländer et al., 2008) that incorporates an adapted scoring system for plant miRNA biogenesis. Non-annotated small RNA sequences as well as those annotated as mature miRNAs were used as input for miRNA prediction. Small RNA sequences were aligned to the pokeweed mRNA transcriptome (Neller et al., 2016) and only perfect, full-length alignments were retained. For each alignment, a potential precursor sequence of 250 nt in length was excised and its secondary structure and minimum free energy were determined with RNAfold from the ViennaRNA package (v. 2.1.7; Lorenz et al., 2011). Mature miRNAs were predicted based on compatibility of the small RNA with the structure of its precursor. As per miRDeep-P notation, if both 5p and 3p miRNA reads were detected then the most abundant strand was labeled as the mature miRNA. The less abundant strand was defined as the star sequence. Precursor structures were visualized with The UEA Small RNA Workbench (Stocks et al., 2012). Following miRNA prediction, BLASTn-short (Altschul et al., 1997) was used to identify conserved sequences (E ≤ 0.001) against all mature sequences from the plant miRNA database. Remaining, non-annotated sequences were considered unique to pokeweed.
miRNA Differential Expression Analysis
Raw sequencing counts of predicted miRNAs were used for differential expression analysis with the Bioconductor package EdgeR (Robinson et al., 2010) in the statistical program R (R Development Core Team, 2011). Prior to analysis, any miRNA with an abundance of less than 1 read per million (RPM) in each group of replicates was removed. Subsequently, libraries were normalized based on their total number of counts, and the common and tag-wise dispersions were estimated with default parameters. An Exact Test was conducted to detect miRNAs with significant differences in expression between control and JA-treated groups (FDR < 0.05). Each treatment group consisted of three biological replicates.
Prediction of miRNA Targets
Predicted miRNAs were aligned to the pokeweed mRNA transcriptome and putative targets were identified using the Plant Small RNA Target Analysis Server with default parameters (Dai and Zhao, 2011; http://plantgrn.noble.org/psRNATarget/). This program predicts miRNA/target pairs based on reverse complementary matching and target site accessibility, incorporating plant-specific features of target recognition.
Correlation of miRNA/target Expression
A paired design was employed for small RNA and mRNA sequencing; that is, both small RNA and mRNA fractions were sequenced for each biological replicate of each treatment group. The Pearson correlation coefficient (PCC) of each miRNA/target pair was calculated from normalized expression data in RPM and fragments per kilobase per transcript per million mapped reads (FPKM), respectively. Transcript quantification and differential expression results of miRNA targets were obtained from Neller et al. (2016).
Functional Analysis of miRNA Targets
Annotation information of miRNA targets was obtained from Neller et al. (2016). Gene ontology (GO) term distribution and enrichment analysis were conducted with Blast2GO software (Conesa et al., 2005). KOBAS 2.0 (Xie et al., 2011) was used to identify enriched KEGG pathways. For enrichment analysis, the test set comprised targets with highly correlated miRNAs (PCC > |0.8|), for which both the target and miRNA were differentially expressed (FDR < 0.05), and the reference set consisted of all miRNA targets.
Construction of a miRNA/target Interaction Network
Cytoscape software (v. 3.3.0; Shannon et al., 2003) was used to generate a network of highly correlated miRNA/target pairs (PCC > |0.8|), for which both members were differentially expressed (FDR < 0.05).
Quantitative Reverse Transcription PCR (qRT-PCR)
Reverse transcription was performed by combining 250 ng of total pokeweed RNA with 1 μL of 10 μM miRNA-specific stem-loop reverse primer and 1 μL of 10 mM dNTPs in a total volume of 13 μL. The solution was incubated at 65°C for 10 min before cooling on ice for 3 min. To the mixture was added 4 μL of 5X First Strand Buffer (250 mM Tris-HCl pH 8.3, 375 mM KCl, 15 mM MgCl2), 1.38 μL of dH2O, 1 μL of 0.1 M DTT, 20 units of Murine RNase Inhibitor, and 25 units of Superscript III reverse transcriptase (Invitrogen). The reaction was incubated at 42°C for 1 h and heat inactivated at 70°C for 20 min.
Following reverse transcription, 6 μL of cDNA was combined with 2 μL of 10 μM forward primer, 2 μL of 10 μM reverse primer, 23 μL of dH2O, and 33 μL of 2X Sybr Green Mastermix (ABM). Each reaction was divided into three technical replicates and analyzed in a Qiagen Rotor-gene-Q real time PCR cycler. Ct values were calculated with the ΔΔCt relative quantification method. miR156 was used as an internal control, as it showed stable expression based on RNA-Seq results.
5′ Rapid Amplification of cDNA Ends (RACE)
Polyadenylated RNA was isolated with the NEBNext Poly(A) mRNA Magnetic Isolation Module. The method of 5′ RACE was based on the Generacer 5′ RACE kit, omitting phosphatase and decapping steps. To relax secondary structure, 250 ng of pokeweed mRNA and 250 ng of RNA adapter were combined in a total volume of 10 μL and incubated at 65°C for 5 min. After cooling on ice for 2 min, 2 μL of 10X T4 RNA Ligase Buffer (500 mM Tris-HCl pH 7.5, 100 mM MgCl2, 10 mM DTT), 2 μL of 10 mM ATP, and 10 units of T4 RNA ligase were added and brought to a final volume of 20 μL. The ligation mixture was incubated at 37°C for 1 h, extracted once with acidic phenol/chloroform, precipitated with 5 μg of linear acrylamide, and re-suspended in RNA storage buffer (22.5 mM DTT, 1 mM sodium citrate pH 6.4). Reverse transcription was performed as described above with poly-d(T) adapter primers, followed by one round of non-specific amplification with forward and reverse adapter primers. To isolate cleavage sites, 1 μL of cDNA was used as template for PCR with adapter forward primers and gene-specific reverse primers 3′ of the putative cleavage site. Following one round of nested PCR, the band of interest was gel-purified, cloned, and sequenced. Five cDNA clones were sequenced per validation.
Results
Overview of the Pokeweed Small RNA Transcriptome
A summary of our study is provided in Figure S1. Following quality control processing, a total of 57,171,256 clean reads were obtained from control and JA-treated libraries (Table 1). To identify conserved miRNAs in pokeweed, reads were aligned to mature sequences from the plant miRNA database, with up to three allowable mismatches in alignments. Conserved miRNAs comprised a small percentage (<0.5%) of reads, representing 1,381 unique sequences. Other classes of conserved, non-coding RNA in pokeweed (rRNA, tRNA, snRNA, and snoRNA) were annotated by mapping reads to the sugar beet genome. Sugar beet was selected as a reference since it is the closest related species with a sequenced genome. At most, the combination of these remaining non-coding RNA species accounted for 14% of reads. Mapping rates to the pokeweed mRNA transcriptome averaged 18 and 23% for control and JA libraries, respectively. Although relatively low, these rates agree with the finding that most small RNAs were 24-mer heterochromatic siRNAs (59 and 55% for control and JA-treated libraries, respectively; Figure S2). Small RNAs of this class are derived from non-polyadenylated precursors, so these reads were not expected to align to the mRNA transcriptome. Taken together, most of the pokeweed small RNA transcriptome remains un-annotated when compared with sequences from other plants. This indicates that the Phytolaccaceae family represents a source of novel small RNAs and corresponding targets.
Table 1. Annotation of small RNA sequences in pokeweed.
miRNA Prediction and Differential Expression Analysis
Pokeweed miRNAs were predicted computationally with miRDeep-P, a plant-specific program that incorporates principles of miRNA biogenesis. Briefly, unannotated reads and those annotated as mature miRNA sequences were aligned to the pokeweed mRNA transcriptome to identify likely miRNA precursors. For each alignment, a flanking region of mRNA was extracted and its secondary structure predicted. Mature miRNAs were identified based on the likelihood that a particular small RNA sequence originated from the biological processing of a miRNA precursor. This strategy led to the prediction of 582 miRNAs in pokeweed; their sequences are provided in Table S1. To identify which of these miRNAs were conserved, a BLASTn-short search was conducted against mature sequences from the plant miRNA database (E-value ≤ 0.001). A total of 24 conserved miRNAs were annotated (Table 2). Among these, the most abundant were miR156, miR169, miR535, and miR396.
Table 2. Conserved miRNAs in pokeweed.
The majority of predicted miRNAs were 21-nt, as expected (Figure 1A); however, this pool had low sequence diversity, as many 21-nt reads had the same sequence (Figure 1B). Conversely, non-canonical 24-nt miRNAs had high sequence diversity and likely targeted numerous different genes, even though they comprised a small fraction of the miRNA pool. These trends were consistent following differential expression analysis to identify JA-responsive miRNAs. That is, the 24-nt miRNAs contained the greatest number of unique sequences in both the total miRNA and JA-responsive miRNA pools.
Figure 1. The length distribution of predicted miRNAs in pokeweed. (A) Reads; (B) unique sequences. DE, differentially expressed; nt, nucleotides.
Of 582 miRNAs, 145 showed significant changes in abundance with JA treatment, relative to control (FDR < 0.05; Figure 2A and Table S2). Among these were two conserved miRNAs, miR172 and miR395, with log2 fold changes of −13.4 and 2.30, respectively. The numbers of up and down-regulated miRNAs were equal, with 73 and 72 sequences, respectively. Interestingly, although a majority, 443 of the 582 miRNAs, were present in both control and JA samples (Figure 2B), most differentially expressed miRNAs were condition-specific (Figure 2C). That is, most were detected in either all control replicates or all treated replicates; only 16 miRNAs (11% of those differentially expressed) were found in all samples. Furthermore, most JA-responsive miRNAs were of low abundance and displayed high fold changes within a narrow range (Figure 2D). We observed a symmetrical distribution of miRNA fold changes (Figure 2E). Of all differentially expressed miRNAs, those that decreased in abundance with JA had an average log2FC and standard deviation of −5.13 ± 1.38, while those that increased with JA had values of 4.92 ± 1.07. Therefore, JA treatment alters sequence diversity within the miRNA pool, and discrete fold changes in expression of JA-responsive miRNAs may be required for biological effect.
Figure 2. Identification of JA-responsive miRNAs in pokeweed. (A) Heat map of expression values (log2RPM, median-centered) of differentially expressed miRNAs (FDR < 0.05). (B,C) Depict treatment-specific expression patterns of all miRNAs and differentially expressed miRNAs, respectively. (D) MA plot indicating all differentially expressed miRNAs (red). (E) Histogram of miRNA fold change; x-axis values indicate upper bin numbers. DE, differentially expressed; FC, fold change.
Identification of JA-responsive miRNA/target Pairs
We used the Plant Small RNA Target Analysis Server to identify pokeweed mRNA transcripts having high sequence complementarity with the 582 predicted miRNAs, which resulted in a total of 25,931 miRNA/target pairs. To identify the most biologically relevant interactions, the list was filtered to include only JA-responsive miRNAs and targets (FDR < 0.05), resulting in a subset of 516 miRNA/target pairs. That is, for each pair, both miRNA and target were differentially expressed with JA. The rationale for limiting our analysis in this manner was to include only interactions that could potentially mediate a functional response in the plant through the alteration of mRNA abundance, given that plant miRNAs act predominately through target mRNA cleavage.
The expression change of each miRNA and target forming the 516 JA-responsive pairs is shown in Figure 3A. The number of negative and positive miRNA/target expression correlations was similar, with 274 pairs showing opposite JA-induced changes (i.e. miRNA ↑ and target ↓ or vice versa) and 242 pairs with parallel changes (i.e. miRNA and target both ↑ or ↓). Most strikingly, while miRNA expression change was limited to a narrow range, target change was much more variable. Interestingly, plotting all miRNA/target pairs revealed that the non-random distribution was inherent in the overall miRNA population, as shown by trimodal clustering (Figure 3B). These results indicate that JA treatment induces a certain threshold of miRNA expression change, for those miRNAs that are JA-responsive, suggesting an optimal miRNA fold change for effective JA control.
Figure 3. Comparison of miRNA and target expression changes. The fold change (FC) of each miRNA was plotted against that of its predicted target. Blue and orange indicate positive and negative correlations, respectively, of miRNA/target expression changes. (A) Differentially expressed miRNAs and targets (FDR < 0.05). (B) All miRNAs and targets.
This RNA-Seq study incorporated paired small RNA and mRNA samples; that is, each RNA sample was sequenced to quantify both small RNA and mRNA levels. This paired design, combined with the use of independent biological replicates, enabled the identification of miRNA/target pairs with highly correlated expression. The PCC was calculated for each of the 516 miRNA/target pairs described above and the distribution of correlations is shown in Figure S3. There was a clear enrichment in moderate and high correlations; this was expected as differential expression analysis and sequence complementarity filters had already been applied to retain the most likely miRNA/target candidates. Specifically, we were interested in highly correlated pairs, defined as those with PCC > |0.8|. This cut-off reduced the number of pairs from 516 to 171, comprising 97 positively correlated pairs and 74 negatively correlated ones. These 171 pairs are most biologically relevant since the expression of each miRNA is associated with the expression of its target.
Functional Analysis of Highly Correlated miRNA/target Pairs
Following the identification of highly correlated miRNA/target pairs, our next goal was to functionally characterize them. Target annotations were obtained from our previous study of the pokeweed mRNA transcriptome (Table S3). Based on KEGG pathway enrichment, the top pathways were associated with plant stress and defense responses, including “peroxisome,” “alpha-Linolenic acid metabolism,” and “plant-pathogen interaction” (Table 3). The corresponding targets mapping to each enriched pathway revealed key genes that may be regulated by miRNAs. Notably, these included two genes associated with JA biosynthesis (OPR3 and AOC3) as well as an integrator of JA and ethylene signaling (Ethylene-responsive transcription factor 1B, ERF1B). Although miRNA-mediated control of transcription factors involved in JA biosynthesis has been reported, our results indicate that miRNAs directly target genes of the JA biosynthesis pathway.
Table 3. Top pathways involving highly correlated, JA-responsive miRNA/target pairs.
Since pokeweed is a non-model plant with limited resources for annotation, some targets could not be mapped to GO terms or KEGG pathways, or their functional associations lacked substantial detail. To account for this, we manually assessed the 171 highly correlated miRNA/target pairs and cross-referenced the literature to find genes associated with plant stress and defense responses. The results are summarized in Table 4, revealing novel targets involved in hormone signaling, hypersensitive response, lignin biosynthesis, and oxidative stress response. Importantly, some miRNAs in Table 4 have more than one target, such as 110629_x73, which is predicted to regulate two genes encoding anti-microbial proteins, and 419489_x11, which targets genes involved in both JA signaling and cellular detoxification. Taken together, it is evident that targets of highly correlated miRNAs contribute to various aspects of biotic and abiotic stress responses.
Table 4. Identification of highly correlated, JA-responsive miRNA/target pairs implicated in plant stress and defense responses.
We constructed an interaction network to visualize the 171 highly correlated, differentially expressed miRNA/target pairs (Figure S4). The network includes 67 unique miRNAs, the majority of which (40) had more than one target. Major hubs are located toward the top of the network; these are cases where a single miRNA has many targets, and different hubs are often connected by a common target. Enrichment analysis was conducted to identify overrepresented GO terms associated with major hubs; however, no significant terms were found. This suggests that though individual miRNA/target pairs are involved in plant defense, these miRNAs also control different mRNAs with diverse functions, rather than all targets of a single miRNA being involved in a single response (Table 3).
Precursor structures of the 67 miRNAs are provided in Figure S5. Although variability was evident in terms of overall precursor shape and positioning of mature sequence, the mature miRNA tended to be in the middle or upper stem of a stable hairpin structure. We did note a few unexpected cases wherein the mature miRNA was mostly or completely contained within a non-hairpin loop; these cases may have arisen from limitations in the RNA folding and/or miRNA prediction programs.
Validation of miRNA Differential Expression and Functional Activity
To validate RNA-Seq results, we measured the expression of eight miRNAs (98462_x2365, 297830_x3990, 66102_x232, 110629_x73, 114897_x905, 144454_x98, 62255_x144, 50891_x48) by qRT-PCR from control and JA-treated plants (Figure 4). These miRNAs were selected to represent a range of fold changes (9.5 × 10−5 to 9.5) and abundances (0.4–44 RPM), and all had a highly correlated target that was functionally annotated through BLASTx. The R2 correlation value was 0.88, indicating high correspondence between the two methods of differential expression analysis. The fold-change of transcript levels quantified with RNA-Seq were higher than those quantified with qRT-PCR, as observed previously (Baran-Gale et al., 2015). This difference is likely due to lower amplification efficiency of stem-loop primers used in qRT-PCR, compared to adapter-specific primers used for RNA-Seq. To verify miRNA-induced cleavage of a specific mRNA, four targets were validated by 5′RACE: c56431_g2_i2 (squamosa promoter-binding-like protein 2, SPL2), c59355_g7_i1 (thymidine kinase, TK), c19572_g1_i1 (OPR3), and c17159_g1_i1 (transcription factor bHLH25, BHLH25) (Figure 5A). The SPL2/miR156 pair was selected as a positive control, as it is a well-established interaction. mRNA cleavage was detected at the expected site, between the 10th and 11th nucleotide of the miRNA binding site, in all cases except for OPR3, where cleavage occurred between the 11th and 12th nucleotide. OPR3 and BHLH25 were targeted by pokeweed-specific miRNAs (98462_x2365 and 144454_x98, respectively) and represent two confirmed interactions from Table 4; both miRNAs were also validated through qRT-PCR (above). TK was identified unexpectedly as a novel target of the well-conserved miR156 and was significantly JA-responsive based on mRNA-Seq results. Precursor structures of miR156 and the two novel miRNAs are shown in Figure 5B. Both miR156 and miR 98462_x2365 had detectable star sequences while 144454_x98 did not, likely owing to its lower abundance. All three precursors adopted stable hairpin formations.
Figure 4. Validation of JA-responsive miRNAs. The correlation of JA-induced expression changes obtained from RNA-Seq and qRT-PCR is shown for 8 miRNAs. Results for qRT-PCR represent the mean from two or three independent biological replicates for each miRNA. Bars indicate the standard error associated with qRT-PCR.
Figure 5. Validation of miRNA target cleavage. (A) miRNA-induced cleavage of mRNA targets was validated through 5′RACE. The upper sequence indicates the mRNA target and the lower sequence is its miRNA. Fractions indicate cleavage sites from individual cDNA clones. From top to bottom, miRNAs are: miR156, miR98462_x2365, miR144454_x98, and miR156. (B) Precursor structures are shown for each miRNA. The mature miRNA sequence is indicated in green. If a miRNA star strand was detected, it was indicated in blue. For each structure, the corresponding miRNA ID, precursor mRNA ID, and minimum free energy are provided.
Discussion
We have performed a large-scale miRNA analysis of pokeweed, with emphasis on JA-induced genes that contribute to biotic and abiotic stress tolerance. A comprehensive resource of 57,171,256 high-quality small RNA reads has been generated, representing the first available small RNA transcriptome of any Phytolaccaceae family member. Through miRNA prediction and annotation, we identified conserved and novel miRNAs in pokeweed. Differential expression analysis revealed that the majority of JA-responsive miRNAs were low in abundance, condition-specific, and exhibited a narrow range of expression change. Our paired RNA-Seq experimental design allowed the prediction of high-confidence miRNA targets and construction of an interaction network to visualize the global connection between miRNAs and their different targets.
Pokeweed miRNAs Have Low Sequence Conservation With Those of Other Plants
Unexpectedly, less than 0.5% of pokeweed small RNA reads represented conserved miRNAs, based on comparison with available plant sequences. Compared with other primary studies of non-model plants, this number seems low. For example, in black pepper, 11% of small RNA reads were conserved miRNAs, with up to two mismatches in alignments (Asha et al., 2016). In Vriesea carinata, 16% of small RNA reads were conserved miRNAs, with only perfect alignments permitted (Guzman et al., 2013). Since this is the first miRNA study of any Phytolaccaceae member, low sequence conservation in pokeweed may be a result of species specificity. Although the majority of plant miRNAs are species-specific (Chávez Montes et al., 2014), the low level of conservation in pokeweed is indicative of a relatively unique miRNA pool.
We performed a stringent miRNA annotation by requiring precursor evidence in addition to mature sequence conservation. Using this strategy, 24 conserved miRNAs were identified in pokeweed, and their expression levels were similar to those observed in other plants (Chávez Montes et al., 2014). Although 1,381 sequences were conserved, a miRNA precursor was predicted for only 1.6%. Others have reported similar findings. For example, in black pepper, 33,350 miRNA sequences were conserved but precursors were found for only 50, or 0.15% (Asha et al., 2016). Despite the common strategy of identifying miRNAs based on sequence conservation alone, our analysis indicates that this method may overestimate the number of true miRNAs.
We also observed non-canonical miRNAs that were 24-nt in length and had high sequence diversity relative to 21-mers. Similar findings have been reported by others (Wan et al., 2012; Guzman et al., 2013; Hackenberg et al., 2013; Shuai et al., 2013). Longer miRNAs tend to be species-specific and are reflective of newly evolved “proto-miRNAs” that undergo imprecise processing by DCL enzymes; interestingly, 24-nt miRNAs can enter the heterochromatic siRNA pathway to direct chromatin modifications of target genes (Axtell, 2013).
JA Treatment Induces Defined Changes in the miRNA Transcriptome
Differential expression analysis revealed that JA-responsive miRNAs tended to be condition-specific and low in abundance, with narrowly clustered expression patterns. Condition specificity and low abundance may indicate that, at certain times, sequence diversity of the miRNA pool is more biologically relevant than miRNA abundance. For example, low abundant miRNAs may be early responders to stress. In support of this idea, Pandey et al. (2008) reported only 2.4% overlap between the small RNA populations of untreated Nicotiana attenuata plants and those having undergone 45 minutes of insect elicitation, which strongly upregulated the JA response. Since low-abundant miRNAs tend to be species-specific (Chávez Montes et al., 2014), analysis of corresponding targets may lead to the identification of key genes contributing to distinctive traits.
Our expression results suggest that a certain threshold of miRNA change is necessary to mediate an effective response to JA. Specifically, log2FC values of JA-responsive miRNAs clustered around +5 and −5. JA-associated thresholding is further supported by the finding that log2FC values of insect-induced small RNAs in N. attenuata ranged from +4 to −4 (Pandey et al., 2008). Clustering of miRNAs has also been reported. In cotton under salt stress, Peng et al. (2014) observed four discrete groups when the expression of salt-responsive miRNAs was plotted against that of targets; such clustering existed for both positively and negatively correlated pairs. These patterns may reflect co-regulated miRNA expression owing to control by common transcription factors, which are themselves differentially expressed.
Although we identified many JA-responsive miRNAs, only two were conserved with other plants: miR172 and miR395, which decreased and increased with JA, respectively. In contrast, others have reported several conserved, JA-responsive miRNAs (Qiu et al., 2009; Bozorov et al., 2012; Zhang et al., 2012; Xu et al., 2016; Shen et al., 2017). One reason for this discrepancy may be our stringent requirements for differential expression, owing to the availability of biological replicates. Other important differences between our work and that cited above include the use of plants deficient in JA biosynthesis, as well as the investigation of cell cultures rather than whole plants. Although miR172 was reported as JA-responsive by others, its magnitude and direction of change differed across studies. In contrast, miR395 has not been reported as JA-responsive; its most well-established role to date is in sulfur assimilation (Matthewman et al., 2012). Interestingly, miR395 increased in Arabidopsis after heavy metal treatment and contributed to cadmium detoxification in Brassica napus (Jagadeeswaran et al., 2013; Zhang et al., 2013). Given that heavy metal stress is associated with biosynthesis of a JA precursor (Foroughi et al., 2014), and JA treatment reduces heavy metal stress (Singh and Shah, 2014), there is a likely association between miR395 and JA.
JA-responsive miRNAs Are Both Positively and Negatively Correlated With Targets
To identify the top miRNA/target candidates, three factors were considered: differential expression, sequence complementarity, and correlation analysis. Pearson correlation is widely used to infer miRNA/target relationships (Muniategui et al., 2013). However, since high correlation does not necessarily indicate causation, we applied this filter after first considering more biologically relevant parameters. Although there is a tendency in literature to report only pairs with negatively correlated expression, we observed a relatively equal distribution of positive and negative correlations. This finding has been observed in other plant studies, reinforcing the concept that both dynamics are prevalent and worthy of consideration (Lopez-Gomollon et al., 2012; Peng et al., 2014; Wen et al., 2016). Positive correlations between the expression of a miRNA and its target can arise from miRNA-mediated spatial restriction of a target (Kidner and Martienssen, 2004; Nikovics et al., 2006; Levine et al., 2007; Kawashima et al., 2009). In Arabidopsis roots, a positive temporal and negative spatial correlation was found for the expression of miR395 and its target, sulfur transporter SULTR2;1. Authors hypothesized that under sulfur starvation, phloem-specific co-expression of miR395/SULTR2;1 restricts expression of SULTR2;1 to the xylem (Kawashima et al., 2009). As another example, the extent of leaf serration in Arabidopsis is controlled by spatial restriction of MIR164A-targeted CUC2 (Nikovics et al., 2006).
Additionally, it is thought that miRNAs operate through two independent mechanisms that are not mutually exclusive: expression tuning, whereby they modify the mean target abundance, and expression buffering, in which they reduce the variance of target abundance around a preset mean (Wu et al., 2009). Negative miRNA/target expression correlations could arise from the first mechanism and positive correlations from the second. Expression buffering consists of negative feedback loops and incoherent feed forward loops (Tsang et al., 2007; Wu et al., 2009). Consider a miRNA, “miR” and its target gene, “T.” As summarized by Wu et al., in a negative feedback loop, protein T is a positive regulator of miR, such that increased T results in increased miR, dampening oscillations in the expression of either gene. An example of this dynamic has been shown within the miRNA biogenesis pathway itself: miR162 targeting of DCL1 (Xie et al., 2003). In an incoherent feed forward loop, consider an additional gene, “A”, whose product increases levels of both T and miR. Since A regulates the level of T, any change in the level of A would also affect that of T. However, variation in A is also passed to miR, which then buffers the fluctuation in T. An interesting example here is the maintenance of AGO1 homeostasis by (1) miR168 targeting; (2) transcriptional co-expression of miR168 and AGO1; and (3) preferential stabilization of miR168 by AGO1. This example combines two negative feedback loops (1 and 3), with an incoherent feed forward loop (2). Though positive correlations of miRNA and target gene levels are not commonly described, examples exist of their biological relevance.
We validated the miRNA-induced cleavage of a positively correlated target, BHLH25. Little is known about the biological role of this transcription factor, but it is reportedly JA-responsive (Heim et al., 2003) and transgenic Arabidopsis plants overexpressing this gene have altered root and shoot morphology and greater susceptibility to cyst nematode infection (Jin et al., 2011). Increased expression of the miRNA in tandem with BHLH25 could buffer levels of this transcription factor during JA induction. For example, an increase in BHLH25 protein could increase the expression of its own miRNA directly (negative feedback loop), and/or an upstream transcription factor could simultaneously increase expression of both BHLH25 and its miRNA (incoherent feed forward loop).
JA-responsive miRNA/target Pairs Are Associated With Multiple Levels of Plant Defense
We identified putative miRNA targets involved in JA biosynthesis and signaling. Two genes encoding JA biosynthesis enzymes were identified, AOC3 and OPR3, and the miRNA-induced cleavage of OPR3 mRNA was verified by 5′ RACE. OPR3 and more recently, OPR2, were shown to be required for JA biosynthesis in Arabidopsis (Stintzi and Browse, 2000; Chini et al., 2018). Although miR319 targets TCP transcription factors and indirectly regulates JA biosynthesis (Schommer et al., 2008), this is the first reported miRNA-induced cleavage of a transcript encoding a JA biosynthesis enzyme. Others have predicted complementary miRNA/target pairs comprising JA biosynthesis enzymes upstream in the pathway, namely LOX2/3 and AOS; however, these interactions have not been functionally validated (Pandey et al., 2008). Given that AOC3 and OPR3 mRNA abundance increased with JA, while their corresponding miRNAs decreased, such regulation would contribute to the known positive feedback loop of JA biosynthesis (Wasternack and Hause, 2013). Interestingly, although miR319 was identified in pokeweed, it was not JA-responsive. In tomato, JA treatment significantly reduced miR319 expression, which was associated with an increase in TCP4 mRNA (Zhao et al., 2015). miR319 targets several class II TCPs, which function as transcriptional activators of LOX2 (Schommer et al., 2008). We hypothesize that the existence of pokeweed-specific, JA-responsive miRNAs reflects an alternative strategy to regulate JA biosynthesis. This hypothesis is further supported by the absence of miR319 binding sites in any TCP-encoding mRNA transcript in pokeweed. We also predicted that ERF1B is targeted by a miRNA; the encoded transcription factor is a key mediator of JA signal transduction that also integrates ethylene-induced signals (Wasternack and Hause, 2013). To the best of our knowledge, this is the first report that OPR3 mRNA is cleaved by a miRNA, and that AOC3 and ERF1B are putative miRNA targets.
Apart from JA biosynthesis and signaling, we identified putative targets having more widespread roles in stress tolerance. These are summarized in Table 4 and include several factors involved in salicylic or abscisic acid (ABA) signaling, with well-established roles in pathogen and environmental stress resistance, and others associated with classic growth hormones. Analysis of the expression dynamics of these targets provides insight into the role of miRNAs in mediating JA crosstalk with other hormone pathways. For example, we identified two miRNA/target pairs involving negative regulators of ABA signaling: Probable protein phosphatase 2C 73 and U-box domain-containing protein 19. Both genes had increased expression with JA, suggesting a dampening of the ABA pathway. This agrees with the known relationship between ABA and pathogen susceptibility, mediated by the transcription factor MYC2, which activates the ABA pathway and inhibits the JA pathway (Derksen et al., 2013). The ABA pathway contributes heavily to environmental stress tolerance (e.g. drought), so these miRNA/target interactions could modulate the trade-off between abiotic and biotic stress responses. A separate trade-off exists between stress response and growth, and JA inhibits developmental processes including seed germination and root/shoot growth (Wasternack and Feussner, 2017). We identified a miRNA/target pair comprising Auxin response factor 1, which is a repressor of auxin signaling (Ellis et al., 2005). Expression of this gene increased with JA, suggesting inhibition of auxin-associated gene expression and related growth processes. The involvement of JA-responsive miRNAs in regulating Auxin response factor genes has been reported by others, indicating that these genes represent a conserved family of targets (Bozorov et al., 2012; Xu et al., 2016; Shen et al., 2017). Finally, we identified targets involved in redox activities: Glutathione S-transferase U9 and Isocitrate dehydrogenase, which may be relevant to the ability of pokeweed to hyperaccumulate heavy metals (Peng et al., 2008; Liu et al., 2010; Zhao et al., 2012). Notably, exogenous application of JA reduces heavy metal-associated oxidative stress (Singh and Shah, 2014). Taken together, our results suggest that miRNA-mediated regulation of the JA response contributes to widespread stress resistance activities in pokeweed. We chose to study pokeweed because it is a non-model plant known to resist biotic and abiotic pressure; however, the JA-responsive miRNA/target pairs identified in this study may have agricultural relevance in terms of improving resistance to pathogen and environmental stresses in other plants.
Data Availability
Raw small RNA reads have been submitted to the SRA under project # SRP069141.
Author Contributions
KN and KH conceived the design and coordination of the study. KN performed bioinformatic analyses, statistical testing, and drafted the manuscript. AK carried out qRT-PCR validations and 5′RACE cloning. JG determined miRNA/target expression correlation and generated the interaction network. KH edited the manuscript. All authors read and approved the final manuscript.
Funding
This work was funded by a Discovery Grant to KH from the Natural Sciences and Engineering Research Council of Canada, and CGS-M and OGS Scholarships to KN and AK.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
The authors are grateful to Dr. Yutaka Amemiya, Manager, Genomics Core Facility, Next Generation Sequencing Lab, Sunnybrook Research Institute, for small RNA sequencing and for his assistance and advice with sequence analysis.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00589/full#supplementary-material
Abbreviations
DE, differentially expressed; FC, fold change; FDR, false discovery rate; FPKM, fragments per kilobase per transcript per million mapped reads; GO, gene ontology; JA, jasmonic acid; KEGG, Kyoto encyclopedia of genes and genomes; PAP, pokeweed antiviral protein; PCC, Pearson correlation coefficient; qRT-PCR, quantitative reverse transcription PCR; RACE, rapid amplification of cDNA ends; RPM, reads per million.
References
Altschul, S. F., Madden, T. L., Schäffer, A. A., Zhang, J., Zhang, Z., Miller, W., et al. (1997). Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 25, 3389–3402. doi: 10.1093/nar/25.17.3389
Asha, S., Sreekumar, S., and Soniya, E. V. (2016). Unravelling the complexity of microRNA-mediated gene regulation in black pepper (Piper nigrum L.) using high-throughput small RNA profiling. Plant Cell Rep. 35, 53–63. doi: 10.1007/s00299-015-1866-x
Axtell, M. J. (2013). Classification and comparison of small RNAs from plants. Annu. Rev. Plant Biol. 64, 137–159. doi: 10.1146/annurev-arplant-050312-120043
Baran-Gale, J., Kurtz, C. L., Erdos, M. R., Sison, C., Young, A., Fannin, E. E., et al. (2015). Addressing bias in small RNA library preparation for sequencing: a new protocol recovers miRNAs that evade capture by current methods. Front. Genet. 6:352. doi: 10.3389/fgene.2015.00352
Beffa, R. S., Hofer, R. M., Thomas, M., and Meins, F. (1996). Decreased susceptibility to viral disease of [beta]-1,3-glucanase-deficient plants generated by antisense transformation. Plant Cell 8, 1001–1011. doi: 10.1105/tpc.8.6.1001
Berrocal-Lobo, M., Segura, A., Moreno, M., López, G., García-Olmedo, F., and Molina, A. (2002). Snakin-2, an antimicrobial peptide from potato whose gene is locally induced by wounding and responds to pathogen infection. Plant Physiol. 128, 951–961. doi: 10.1104/pp.010685
Bozorov, T. A., Baldwin, I. T., and Kim, S. G. (2012). Identification and profiling of miRNAs during herbivory reveals jasmonate-dependent and -independent patterns of accumulation in Nicotiana attenuata. BMC Plant Biol. 12:209. doi: 10.1186/1471-2229-12-209
Chini, A., Gimenez-Ibanez, S., Goossens, A., and Solano, R. (2016). Redundancy and specificity in jasmonate signalling. Curr. Opin. Plant Biol. 33, 147–156. doi: 10.1016/j.pbi.2016.07.005
Chini, A., Monte, I., Zamarreño, A. M., Hamberg, M., Lassueur, S., Reymond, P., et al. (2018). An OPR3-independent pathway uses 4,5-didehydrojasmonate for jasmonate synthesis. Nat. Chem. Biol. 14, 171–178. doi: 10.1038/nchembio.2540
Chávez Montes, R. A., De Fátima Rosas-Cárdenas, F., De Paoli, E., Accerbi, M., Rymarquis, L. A., Mahalingam, G., et al. (2014). Sample sequencing of vascular plants demonstrates widespread conservation and divergence of microRNAs. Nat. Commun. 5:3722. doi: 10.1038/ncomms4722
Coca, M., and San Segundo, B. (2010). AtCPK1 calcium-dependent protein kinase mediates pathogen resistance in Arabidopsis. Plant J. 63, 526–540. doi: 10.1111/j.1365-313X.2010.04255.x
Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610
Dai, X., and Zhao, P. X. (2011). psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 39, W155–W159. doi: 10.1093/nar/gkr319
Dai, W. D., Bonos, S., Guo, Z., Meyer, W. A., Day, P. R., and Belanger, F. C. (2003). Expression of pokeweed antiviral proteins in creeping bentgrass. Plant Cell Rep. 21, 497–502. doi: 10.1007/s00299-002-0534-0
Davin, L. B., and Lewis, N. G. (2000). Dirigent proteins and dirigent sites explain the mystery of specificity of radical precursor coupling in lignan and lignin biosynthesis. Plant Physiol. 123, 453–461. doi: 10.1104/pp.123.2.453
Derksen, H., Rampitsch, C., and Daayf, F. (2013). Signaling cross-talk in plant disease resistance. Plant Sci. 207, 79–87. doi: 10.1016/j.plantsci.2013.03.004
R Development Core Team, (2011). R: a language and environment for statistical computing. R Found. Stat. Comput. 1:409.
Dohm, J. C., Minoche, A. E., Holtgräwe, D., Capella-Gutiérrez, S., Zakrzewski, F., Tafer, H., et al. (2014). The genome of the recently domesticated crop plant sugar beet (Beta vulgaris). Nature 505, 546–549. doi: 10.1038/nature12817
Earley, K., Smith, M., Weber, R., Gregory, B., and Poethig, R. (2010). An endogenous F-box protein regulates ARGONAUTE1 in Arabidopsis thaliana. Silence 1:15. doi: 10.1186/1758-907X-1-15
Ellis, C. M., Nagpal, P., Young, J. C., Hagen, G., Guilfoyle, T. J., and Reed, J. W. (2005). AUXIN RESPONSE FACTOR1 and AUXIN RESPONSE FACTOR2 regulate senescence and floral organ abscission in Arabidopsis thaliana. Development 132, 4563–4574. doi: 10.1242/dev.02012
Endo, Y., Tsurugi, K., and Lambert, J. M. (1988). The site of action of six different ribosome-inactivating proteins from plants on eukaryotic ribosomes: the RNA N-glycosidase activity of the proteins. Biochem. Biophys. Res. Commun. 150, 1032–1036. doi: 10.1016/0006-291X(88)90733-4
Foroughi, S., Baker, A. J., Roessner, U., Johnson, A. A., Bacic, A., and Callahan, D. L. (2014). Hyperaccumulation of zinc by Noccaea caerulescens results in a cascade of stress responses and changes in the elemental profile. Metallomics 6, 1671–1682. doi: 10.1039/C4MT00132J
Friedländer, M. R., Chen, W., Adamidi, C., Maaskola, J., Einspanier, R., Knespel, S., et al. (2008). Discovering microRNAs from deep sequencing data using miRDeep. Nat. Biotechnol. 26, 407–415. doi: 10.1038/nbt1394
Guzman, F., Almerão, M. P., Korbes, A. P., Christoff, A. P., Zanella, C. M., Bered, F., et al. (2013). Identification of potential miRNAs and their targets in Vriesea carinata (Poales, Bromeliaceae). Plant Sci. 210, 214–223. doi: 10.1016/j.plantsci.2013.05.013
Hackenberg, M., Huang, P. J., Huang, C. Y., Shi, B. J., Gustafson, P., and Langridge, P. (2013). A comprehensive expression profile of micrornas and other classes of non-coding small RNAs in barley under phosphorous-deficient and -sufficient conditions. DNA Res. 20, 109–125. doi: 10.1093/dnares/dss037
Heim, M. A., Jakoby, M., Werber, M., Martin, C., Weisshaar, B., and Bailey, P. C. (2003). The basic helix-loop-helix transcription factor family in plants: a genome-wide study of protein structure and functional diversity. Mol. Biol. Evol. 20, 735–747. doi: 10.1093/molbev/msg088
Islam, W., Noman, A., Qasim, M., and Wang, L. (2018). Plant responses to pathogen attack: small rnas in focus. Int. J. Mol. Sci. 19:515. doi: 10.3390/ijms19020515
Jagadeeswaran, G., Li, Y. F., and Sunkar, R. (2013). Redox signaling mediates the expression of a sulfate-deprivation-inducible miR395 in Arabidopsis. Plant J. 77, 85–96. doi: 10.1111/tpj.12364
Jin, J., Hewezi, T., and Baum, T. J. (2011). The Arabidopsis bHLH25 and bHLH27 transcription factors contribute to susceptibility to the cyst nematode Heterodera schachtii. Plant J. 65, 319–328. doi: 10.1111/j.1365-313X.2010.04424.x
Kalunke, R. M., Tundo, S., Benedetti, M., Cervone, F., De Lorenzo, G., and D'Ovidio, R. (2015). An update on polygalacturonase-inhibiting protein (PGIP), a leucine-rich repeat protein that protects crop plants against pathogens. Front. Plant Sci. 6:146. doi: 10.3389/fpls.2015.00146
Kawashima, C. G., Yoshimoto, N., Maruyama-Nakashita, A., Tsuchiya, Y. N., Saito, K., Takahashi, H., et al. (2009). Sulphur starvation induces the expression of microRNA-395 and one of its target genes but in different cell types. Plant J. 57, 313–321. doi: 10.1111/j.1365-313X.2008.03690.x
Khraiwesh, B., Zhu, J.-K., and Zhu, J. (2012). Role of miRNAs and siRNAs in biotic and abiotic stress responses of plants. Biochim. Biophys. Acta 1819, 137–148. doi: 10.1016/j.bbagrm.2011.05.001
Kidner, C. A., and Martienssen, R. A. (2004). Spatially restricted microRNA directs leaf polarity through ARGONAUTE1. Nature 428, 81–84. doi: 10.1038/nature02366
Klenov, A., Neller, K. C., Burns, L. A., Krivdova, G., and Hudak, K. A. (2016). A small RNA targets pokeweed antiviral protein transcript. Physiol. Plant. 156, 241–251. doi: 10.1111/ppl.12393
Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25. doi: 10.1186/gb-2009-10-3-r25
Levine, E., McHale, P., and Levine, H. (2007). Small regulatory RNAs may sharpen spatial expression patterns. PLoS Comput. Biol. 3, 2356–2365. doi: 10.1371/journal.pcbi.0030233
Li, S., Castillo-González, C., Yu, B., and Zhang, X. (2017). The functions of plant small RNAs in development and in stress responses. Plant J. 90, 654–670. doi: 10.1111/tpj.13444
Liu, X., Peng, K., Wang, A., Lian, C., and Shen, Z. (2010). Cadmium accumulation and distribution in populations of Phytolacca americana L. and the role of transpiration. Chemosphere 78, 1136–1141. doi: 10.1016/j.chemosphere.2009.12.030
Liu, Y. C., Wu, Y. R., Huang, X. H., Sun, J., and Xie, Q. (2011). AtPUB19, a U-Box E3 ubiquitin ligase, negatively regulates abscisic acid and drought responses in Arabidopsis thaliana. Mol. Plant 4, 938–946. doi: 10.1093/mp/ssr030
Lopez-Gomollon, S., Mohorianu, I., Szittya, G., Moulton, V., and Dalmay, T. (2012). Diverse correlation patterns between microRNAs and their targets during tomato fruit development indicates different modes of microRNA actions. Planta 236, 1875–1887. doi: 10.1007/s00425-012-1734-7
Lorenz, R., Bernhart, S. H., Höner zu Siederdissen, C., Tafer, H., Flamm, C., Stadler, P. F., et al. (2011). ViennaRNA Package 2.0. Algorithms Mol. Biol. 6:26. doi: 10.1186/1748-7188-6-26
Lorenzo, O., Piqueras, R., Sánchez-Serrano, J. J., and Solano, R. (2003). ETHYLENE RESPONSE FACTOR1 integrates signals from ethylene and jasmonate pathways in plant defense. Plant Cell 15, 165–178. doi: 10.1105/tpc.007468
Matthewman, C. A., Kawashima, C. G., Húska, D., Csorba, T., Dalmay, T., and Kopriva, S. (2012). MiR395 is a general component of the sulfate assimilation regulatory network in Arabidopsis. FEBS Lett. 586, 3242–3248. doi: 10.1016/j.febslet.2012.06.044
Mittler, R., Kim, Y., Song, L., Coutu, J., Coutu, A., Ciftci-Yilmaz, S., et al. (2006). Gain- and loss-of-function mutations in Zat10 enhance the tolerance of plants to abiotic stress. FEBS Lett. 580, 6537–6542. doi: 10.1016/j.febslet.2006.11.002
Muniategui, A., Pey, J., Planes, F. J., and Rubio, A. (2013). Joint analysis of miRNA and mRNA expression data. Brief. Bioinformatics. 14, 263–278. doi: 10.1093/bib/bbs028
Neller, K. C. M., Klenov, A., and Hudak, K. A. (2016). The pokeweed leaf mRNA transcriptome and its regulation by jasmonic acid. Front. Plant Sci. 7:283. doi: 10.3389/fpls.2016.00283
Nersissian, A. M., Valentine, J. S., Immoos, C., Hill, M. G., Hart, P. J., Williams, G., et al. (1998). Uclacyanins, stellacyanins, and plantacyanins are distinct subfamilies of phytocyanins: plant-specific mononuclear blue copper proteins. Protein Sci. 7, 1915–1929. doi: 10.1002/pro.5560070907
Niggeweg, R., Thurow, C., Weigel, R., Pfitzner, U., and Gatz, C. (2000). Tobacco TGA factors differ with respect to interaction with NPR1, activation potential and DNA-binding properties. Plant Mol. Biol. 42, 775–788. doi: 10.1023/A:1006319113205
Nikovics, K., Blein, T., Peaucelle, A., Ishida, T., Morin, H., Aida, M., et al. (2006). The balance between the MIR164A and CUC2 genes controls leaf margin serration in Arabidopsis. Plant Cell 18, 2929–2945. doi: 10.1105/tpc.106.045617
Pandey, S. P., Shahi, P., Gase, K., and Baldwin, I. T. (2008). Herbivory-induced changes in the small-RNA transcriptome and phytohormone signaling in Nicotiana attenuata. Proc. Natl. Acad. Sci. U.S.A. 105, 4559–4564. doi: 10.1073/pnas.0711363105
Peng, Z., He, S., Gong, W., Sun, J., Pan, Z., Xu, F., et al. (2014). Comprehensive analysis of differentially expressed genes and transcriptional regulation induced by salt stress in two contrasting cotton genotypes. BMC Genomics 15:760. doi: 10.1186/1471-2164-15-760
Peng, K., Luo, C., You, W., Lian, C., Li, X., and Shen, Z. (2008). Manganese uptake and interactions with cadmium in the hyperaccumulator-Phytolacca Americana L. J. Hazard. Mater. 154, 674–681. doi: 10.1016/j.jhazmat.2007.10.080
Qi, Y., Tsuda, K., Nguyen, L. V., Wang, X., Lin, J., Murphy, A. S., et al. (2011). Physical association of Arabidopsis Hypersensitive Induced Reaction proteins (HIRs) with the immune receptor RPS2. J. Biol. Chem. 286, 31297–31307. doi: 10.1074/jbc.M110.211615
Qiu, D., Pan, X., Wilson, I. W., Li, F., Liu, M., Teng, W., et al. (2009). High throughput sequencing technology reveals that the taxoid elicitor methyl jasmonate regulates microRNA expression in Chinese yew (Taxus chinensis). Gene 436, 37–44. doi: 10.1016/j.gene0.2009.01.006
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Schaller, F., Biesgen, C., Müssig, C., Altmann, T., and Weiler, E. W. (2000). 12-Oxophytodienoate reductase 3 (OPR3) is the isoenzyme involved in jasmonate biosynthesis. Planta 210, 979–984. doi: 10.1007/s004250050706
Schommer, C., Palatnik, J. F., Aggarwal, P., Chételat, A., Cubas, P., Farmer, E. E., et al. (2008). Control of jasmonate biosynthesis and senescence by miR319 targets. PLoS Biol. 6:e230. doi: 10.1371/journal.pbio.0060230
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Shen, E. M., Singh, S. K., Ghosh, J. S., Patra, B., Paul, P., Yuan, L., et al. (2017). The miRNAome of Catharanthus roseus: identification, expression analysis, and potential roles of microRNAs in regulation of terpenoid indole alkaloid biosynthesis. Sci. Rep. 7:43027. doi: 10.1038/srep43027
Shuai, P., Liang, D., Zhang, Z., Yin, W., and Xia, X. (2013). Identification of drought-responsive and novel Populus trichocarpa microRNAs by high-throughput sequencing and their targets using degradome analysis. BMC Genomics 14:233. doi: 10.1186/1471-2164-14-233
Sibout, R., Eudes, A., Mouille, G., Pollet, B., Lapierre, C., Jouanin, L., et al. (2005). CINNAMYL ALCOHOL DEHYDROGENASE-C and -D are the primary genes involved in lignin biosynthesis in the floral stem of Arabidopsis. Plant Cell 17, 2059–2076. doi: 10.1105/tpc.105.030767
Silva, N. F., and Goring, D. R. (2002). The proline-rich, extensin-like receptor kinase-1 (PERK1) gene is rapidly induced by wounding. Plant Mol. Biol. 50, 667–685. doi: 10.1023/A:1019951120788
Singh, I., and Shah, K. (2014). Exogenous application of methyl jasmonate lowers the effect of cadmium-induced oxidative injury in rice seedlings. Phytochemistry 108, 57–66. doi: 10.1016/j.phytochem.2014.09.007
Stintzi, A., and Browse, J. (2000). The Arabidopsis male-sterile mutant, opr3, lacks the 12-oxophytodienoic acid reductase required for jasmonate synthesis. Proc. Natl. Acad. Sci. U.S.A. 97, 10625–10630. doi: 10.1073/pnas.190264497
Stocks, M. B., Moxon, S., Mapleson, D., Woolfenden, H. C., Mohorianu, I., Folkes, L., et al. (2012). The UEA sRNA workbench: a suite of tools for analysing and visualizing next generation sequencing microRNA and small RNA datasets. Bioinformatics 28, 2059–2061. doi: 10.1093/bioinformatics/bts311
Sun, Y. C., Pan, L. L., Ying, F. Z., Li, P., Wang, X. W., and Liu, S. S. (2017). Jasmonic acid-related resistance in tomato mediates interactions between whitefly and whitefly-transmitted virus. Sci. Rep. 7:566. doi: 10.1038/s41598-017-00692-w
To, J. P. C., Haberer, G., Ferreira, F. J., Deruère, J., Mason, M. G., Schaller, G. E., et al. (2004). Type-A Arabidopsis response regulators are partially redundant negative regulators of cytokinin signaling. Plant Cell 16, 658–671. doi: 10.1105/tpc.018978
Tran, L. S. P., Nakashima, K., Sakuma, Y., Simpson, S. D., Fujita, Y., Maruyama, K., et al. (2004). Isolation and functional analysis of Arabidopsis stress-inducible NAC transcription factors that bind to a drought-responsive cis-element in the early responsive to dehydration stress 1 promoter. Plant Cell 16, 2481–2498. doi: 10.1105/tpc.104.022699
Tsang, J., Zhu, J., and van Oudenaarden, A. (2007). MicroRNA-mediated feedback and feedforward loops are recurrent network motifs in mammals. Mol. Cell 26, 753–767. doi: 10.1016/j.molcel.2007.05.018
Umezawa, T., Sugiyama, N., Mizoguchi, M., Hayashi, S., Myouga, F., Yamaguchi-Shinozaki, K., et al. (2009). Type 2C protein phosphatases directly regulate abscisic acid activated protein kinases in Arabidopsis. Proc. Natl. Acad. Sci. U.S.A. 106, 17588–17593. doi: 10.1073/pnas.0907095106
Valderrama, R., Corpas, F. J., Carreras, A., Gómez-Rodríguez, M. V., Chaki, M., Pedrajas, J. R., et al. (2006). The dehydrogenase-mediated recycling of NADPH is a key antioxidant system against salt-induced oxidative stress in olive plants. Plant Cell Environ. 29, 1449–1459. doi: 10.1111/j.1365-3040.2006.01530.x
Wagner, U., Edwards, R., Dixon, D. P., and Mauch, F. (2002). Probing the diversity of the Arabidopsis glutathione S-transferase gene family. Plant Mol. Biol. 49, 515–532. doi: 10.1023/A:1015557300450
Wan, L. C., Wang, F., Guo, X., Lu, S., Qiu, Z., Zhao, Y., et al. (2012). Identification and characterization of small non-coding RNAs from Chinese fir by high throughput sequencing. BMC Plant Biol. 12:146. doi: 10.1186/1471-2229-12-146
Wang, P., Zoubenko, O., and Tumer, N. E. (1998). Reduced toxicity and broad spectrum resistance to viral and fungal infection in transgenic plants expressing pokeweed antiviral protein II. Plant Mol. Biol. 38, 957–964. doi: 10.1023/A:1006084925016
Wasternack, C., and Feussner, I. (2017). The oxylipin pathways: biochemistry and function. Annu. Rev. Plant Biol. 69, 1–24. doi: 10.1146/annurev-arplant-042817-040440
Wasternack, C., and Hause, B. (2013). Jasmonates: biosynthesis, perception, signal transduction and action in plant stress response, growth and development. An update to the 2007 review in Annals of Botany. Ann. Bot. 111, 1021–1058. doi: 10.1093/aob/mct067
Wen, M., Xie, M., He, L., Wang, Y., Shi, S., and Tang, T. (2016). Expression variations of miRNAs and mRNAs in rice (Oryza sativa). Genome Biol. Evol. 8, 3529–3544. doi: 10.1093/gbe/evw252
Wu, C. I., Shen, Y., and Tang, T. (2009). Evolution under canalization and the dual roles of microRNAs - A hypothesis. Genome Res. 19, 734–743. doi: 10.1101/gr.084640.108
Xie, Z., Kasschau, K. D., and Carrington, J. C. (2003). Negative feedback regulation of Dicer-Like1 in Arabidopsis by microRNA-guided mRNA degradation. Curr. Biol. 13, 784–789. doi: 10.1016/S0960-9822(03)00281-1
Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., et al. (2011). KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, 316–322. doi: 10.1093/nar/gkr483
Xu, S., Jiang, Y., Wang, N., Xia, B., Jiang, Y., Li, X., et al. (2016). Identification and differential regulation of microRNAs in response to methyl jasmonate treatment in Lycoris aurea by deep sequencing. BMC Genomics 17:789. doi: 10.1186/s12864-016-2645-y
Yang, X., and Li, L. (2011). miRDeep-P: a computational tool for analyzing the microRNA transcriptome in plants. Bioinformatics 27, 2614–2615. doi: 10.1093/bioinformatics/btr430
Yu, Y., Jia, T., and Chen, X. (2017). The “how” and “where” of plant microRNAs. New Phytol. 216, 1002–1017. doi: 10.1111/nph.14834
Zhang, L. W., Song, J. B., Shu, X. X., Zhang, Y., and Yang, Z. M. (2013). MiR395 is involved in detoxification of cadmium in Brassica napus. J. Hazard. Mater. 250–251, 204–211. doi: 10.1016/j.jhazmat.2013.01.053
Zhang, B., Xie, D., and Jin, Z. (2012). Global analysis of non-coding small RNAs in Arabidopsis in response to jasmonate treatment by deep sequencing technology. J. Integr. Plant Biol. 54, 73–86. doi: 10.1111/j.1744-7909.2012.01098.x
Zhang, Z., Yu, J., Li, D., Zhang, Z., Liu, F., Zhou, X., et al. (2010). PMRD: plant microRNA database. Nucleic Acids Res. 38, D806–D813. doi: 10.1093/nar/gkp818
Zhao, W., Li, Z., Fan, J., Hu, C., Yang, R., Qi, X., et al. (2015). Identification of jasmonic acid-associated microRNAs and characterization of the regulatory roles of the miR319/TCP4 module under root-knot nematode stress in tomato. J. Exp. Bot. 66, 4653–4667. doi: 10.1093/jxb/erv238
Zhao, L., Sun, Y., Le, Cui, S. X., Chen, M., Yang, H. M., Liu, H. M., et al. (2011). Cd-induced changes in leaf proteome of the hyperaccumulator plant Phytolacca americana. Chemosphere 85, 56–66. doi: 10.1016/j.chemosphere.2011.06.029
Zhao, H., Wu, L., Chai, T., Zhang, Y., Tan, J., and Ma, S. (2012). The effects of copper, manganese and zinc on plant growth and elemental accumulation in the manganese-hyperaccumulator Phytolacca americana. J. Plant Physiol. 169, 1243–1252. doi: 10.1016/j.jplph.2012.04.016
Ziegler, J., Stenzel, I., Hause, B., Maucher, H., Hamberg, M., Grimm, R., et al. (2000). Molecular cloning of allene oxide cyclase. The enzyme establishing the stereochemistry of octadecanoids and jasmonates. J. Biol. Chem. 275, 19132–19138. doi: 10.1074/jbc.M002133200
Zoubenko, O., Hudak, K., and Tumer, N. E. (2000). A non-toxic pokeweed antiviral protein mutant inhibits pathogen infection via a novel salicylic acid-independent pathway. Plant Mol. Biol. 44, 219–229. doi: 10.1023/A:1006443626864
Keywords: jasmonic acid, miRNA, Phytolacca americana, plant defense, pokeweed, small RNA, transcriptome
Citation: Neller KCM, Klenov A, Guzman JC and Hudak KA (2018) Integration of the Pokeweed miRNA and mRNA Transcriptomes Reveals Targeting of Jasmonic Acid-Responsive Genes. Front. Plant Sci. 9:589. doi: 10.3389/fpls.2018.00589
Received: 03 January 2018; Accepted: 16 April 2018;
Published: 03 May 2018.
Edited by:
Juan Caballero, Universidad Autónoma de Querétaro, MexicoReviewed by:
Andrea Chini, Consejo Superior de Investigaciones Científicas (CSIC), SpainMarcos Egea-Cortines, Universidad Politécnica de Cartagena, Spain
Copyright © 2018 Neller, Klenov, Guzman and Hudak. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Katalin A. Hudak, hudak@yorku.ca