Genome-wide miRNA expression profiling in potato (Solanum tuberosum L.) reveals TOR-dependent post-transcriptional gene regulatory networks in diverse metabolic pathway

Plant material and growth condition

In this study, the potato seedling of BP12-OE line was used for the sRNA sequencing (Deng et al., 2017). All plants were grown in Murashige & Skoog (MS) medium under 16 h light/8 h dark in growth chambers at 22 °C. Four-week-old potato seedlings growing on the MS media for 48 h with TOR inhibitors (RAP + KU) and DMSO were collected for the sRNA sequencing.

RNA extraction, library preparation and sequencing

The extraction method and quality control of total RNA were according to Deng et al. (2017). The total amount of 3 ug RNA for each sample was used for small RNA library construction. The NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB, USA) was used to generate the sequencing libraries and the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips was used to assess the library quality. The cBot Cluster Generation System using TruSeq SR Cluster kit v3-cBot-HS (Illumina) was used to cluster the index-coded samples. And then, the libraries were sequenced on an Illumina Hiseq 2500/2000 platform. At last, total of 50 bp single-end reads were generated.

Data analysis

Quality control and mapping

Raw reads of fastq format were processed through python scripts and custom perl to remove the low quality reads. The GC content, Q20 and Q30 of the raw reads were calculated. All downstream analyses were according to a range of length of clean reads, such as 18–30 nt for plant. The small RNA tags were mapped to the reference potato genome sequence (http://plants.ensembl.org/Solanum_tuberosum/Search/New?db=core) by Bowtie without mismatch to analyze their distribution and expression on the reference (Xu et al., 2011; Langmead et al., 2009).

Bioinformatic identification of known and novel miRNAs, and their differential expression under TOR inhibitors treatment

The mapped small RNA tags were aligned to the miRNA precursor/mature miRNA of plants and animals in the miRBase20.0 database for looking for known miRNA. The potential miRNAs and their secondary structures were obtained using modified software mirdeep2  (Friedlander et al., 2012) and srna-tools-cli (http://srna-tools.cmp.uea.ac.uk/). The Custom scripts were used to analyze the known miRNA counts and their base bias either on each position or on the first position with certain length.

The rRNA, tRNA, snRNA, snoRNA, repeat and protein-coding genes were deleted from the small RNA libraries based on RepeatMasker and Rfam database. And then the remaining sequences were used for predicting the novel miRNAs with the software miREvo  (Wen et al., 2012) and mirdeep2  (Friedlander et al., 2012) through exploring the secondary structure, the minimum free energy and the Dicer cleavage site. The analysis of novel miRNA counts and their base bias was the same as that of known miRNA.

All miRNA abundances were evaluated and normalized using the tags per million reads (TPM) method through the following criteria  (Zhou et al., 2010): TPM = number of mapped miRNA reads ×10⁶/ number of clean sample reads. Differential expression analysis between the treatment and control was performed using the DESeq R package (1.8.3). The P-values were adjusted using q-value  (Storey, 2003). /log₂fold-change/ >1 and the corrected P-value (q-value) <0.01 were set as the threshold for significantly differential expression by default.

Prediction and annotation of target genes of miRNAs

The putative target genes of miRNAs were predicted by psRobot-tar in psRobot  (Wu et al., 2012) for plants. All the target genes were annotated in the NCBI nr database and potato genome database (http://plants.ensembl.org/Solanum_tuberosum/Search/New?db=core). For further annotating the target gene candidates of differentially expressed miRNAs, Gene Ontology (GO) enrichment and the statistical enrichment of KEGG pathways were implemented using GOseq based Wallenius non-central hyper-geometric distribution and KOBAS, respectively (Mao et al., 2005; Young et al., 2010).

M RNA and miRNA validation by qRT-PCR

Four-week-old potato seedlings growing on the MS media were treated with TOR inhibitors (RAP + KU) and DMSO for 48 h and were collected for qRT-PCR. The mRNA validation was according to  Deng et al. (2017). The primers for qRT-PCR designed by Primer premier 5 software were listed in Table S1. To validate the expression of miRNA, 1 µg of total RNA was reverse transcribed using Mir-X™ miRNA qRT-PCR TB Green™ Kit (Takara). The RT-qPCR analysis was performed according to the manufacturer’s protocol (Mir-X™ miRNA qRT-PCR TB Green™ Kit). The miRNA-specific primers were listed in Table S1 and reverse primers was provided by kit. Potato actin gene was used as constitutive references (Nicot et al., 2005; Payyavula, Singh & Navarre, 2013; Liu et al., 2015). All reactions were conducted with 3 biological replicates, and the relative expression of genes was conducted using the 2^−ΔΔc(t) method.

sRNA sequencing

The basic information of sRNA sequencing data from potato (BP12-OE line) with or without TOR inhibitor (RAP + KU) was shown in Fig. 1 and Table S2. The values of Q30 (bases correct recognition rate >99.9%) of the raw data were more than 97%. After removal of adapter contaminants, low quality reads, small reads and polyN, about 92.53% and 95.76% reads in total raw reads were clean reads in two groups, respectively (Fig. 1A). Finally, about 79.18% and 80.61% of the total clean reads with 18–30 nt read length could be mapped and matched to the potato genome sequence with bowtie, respectively (Langmead et al., 2009). For avoiding some sRNAs mapped to more than one annotation, the mapped clean total sRNAs were annotated according to the following priority rule: known miRNA >rRNA >tRNA >snRNA >snoRNA >repeat >NAT-siRNA >gene >novel miRNA >ta-siRNA (Fig. 1B). Most clean sRNA reads ranged from 21 to 24 nt, with 24 nt being the most abundant group of small RNAs, representing the typical length of Dicer-like protein 3 (DCL3)-derived products (Fig. 1C). The dominance of the 24 nt read length sRNA in potato is consistent with previous studies of other species, such as rice (Morin et al., 2008), cucumber (Cucumis sativus L.) (Mao et al., 2012), tomato (Solanum lycopersicum L.) (Cao et al., 2014), tobacco (Nicotiana tabacum L.) (Yin et al., 2015) and radish (Liu et al., 2018). The second abundant group is 21-nt miRNAs, which have the canonical size derivedfrom DCL1 processing. Most of the first nucleotide from the 5′end of the known and novel miRNAs was uridine (U), which was a bias for AGO1 (Mi et al., 2008), suggesting the important character of miRNAs was not changed under the treatment of TOR specific inhibitors (Fig. S1). In addition, in all the libraries, the total rRNA proportion was less than 60% which was used as a marker as the plant samples quality control.

Identification of known and potential novel miRNAs in potato

In order for identifying the known miRNAs from sRNA libraries in potato, we compared the sRNA information with the known, up-dated and mature plant miRNAs deposited in the miRBase without any mismatching (Kozomara & Griffiths-Jones, 2014). Based on the blastn searches and further sequence analysis, the total of 193 known miRNA were observed and 174 known miRNA exist both in the RAP + KU and control libraries (Table S3). The read counts of the known miRNAs ranged from 0 to 24,000 in the RAP + KU and control libraries, which showed high diversity. About 17 and 22 known miRNAs were found to have more than 10,000 redundancies in the above groups, respectively (Table 1). Among these miRNAs, the miR166a-3p and miR319a-3p had the most reads.

The available softwares miREvo (Wen et al., 2012) and mirdeep2 (Friedlander et al., 2012) were integrated to predict potential novel miRNAs. Only those with stable hairpin structures were considered, because this was an essential characteristic for identification of novel miRNAs. Additionally, the binding locations of Dicer enzymes and free energies were used to evaluate these candidate miRNAs. Ultimately, 79 novel miRNAs were discovered from two potato libraries (Table S3). The novel miRNAs displayed lower expression levels, ranging from 0 to 11,652, when compared with known miRNAs. In the RAP + KU and control group, only 17 and 19 novel miRNAs were found to have more than 100 redundancies, respectively (Table 1). In general, many novel miRNAs discoveries in potato enriched the plant miRNA dataset.

TOR regulates potato miRNA expression

After obtaining the readcounts of all the miRNAs, the quantification and normalization of them were conducted by TPM (transcript per million) (Zhou et al., 2010) (Table S4). The two treatments have the similar distributions of the expression levels of all the miRNAs, and the high correlation (R² = 0.873) of miRNA expression level showed the experimental reliability and reasonable sample selection (Figs. 2A–2B). MiRNAs with TPMs over 60 were regarded as expressing at a very high level and miRNAs with TPMs in the interval 0–1 were deemed not to be expressed at very low levels. About 86% of the total miRNAs (272) were expressed (TPM ≥ 1) and more than 146 miRNAs were highly expressed (TPM >  60) in the two treatments. Among them, one novel and 19 known miRNA were only observed in the sRNA libraries from TOR inhibition (RAP + KU), while 2 novel and 19 known miRNA were only detected in the control DMSO library, indicating TOR may completely inhibit or induce these miRNAs expression in potato.

In order to identify the differentially expressed miRNAs, we compared the expression of the known and novel miRNAs between the TOR inhibitors and control treatment samples using the DEGseq (2010) R package. Total of 86 differentially expressed miRNAs were detected in the RAP + KU vs DMSO group, including 41 up-regulated and 45 down-regulated miRNAs (Fig. 2C). To obtain the miRNA expression patterns, we performed the hierarchical clustering of all the differentially expressed miRNAs based on the log₁₀ (TPMs + 1) for the two sRNA libraries (Fig. 2D). The expression levels of all the differentially expressed miRNAs were listed in Table 2 and these miRNAs (stu-miR5303c, stu-miR8010, stu-miR8024a-3p and stu-miR8028-5p) had remarkable up-regulated differences, while these miRNAs (stu-miR3627-5p, stu-miR7985, stu-miR156f-3p and stu-miR8013) had remarkable down-regulated differences. The differential expression of about 1/3 of total miRNAs suggested that TOR may play a crucial role in the regulation of miRNAs expression in potato.

The prediction of miRNA target genes

In order to investigate the function of all the previously identified miRNAs in different biological processes, their potential target genes were predicted using psRobot_tar in psRobot for plants (Wu et al., 2012). Overall, 192 miRNAs got total 4,127 predicted target genes in 272 identified miRNAs (Table S5). Except the numbered EPlSTUG genes, all the rest (4081) numbered PGSC0003DMG genes were annotated in the website (http://plants.ensembl.org/Solanum_tuberosum/Search/New?db=core) and listed in Table S6. The number of potential targets in the potato is variable for each miRNA from 1 to 772, and most of them have multiple target genes, which are consistent with other reports (Curaba et al., 2012; Xie, Frazier & Zhang, 2011). The top three miRNAs are stu-miR7797a (772), stu-miR5303h (751) and stu-miR5303g (743). The results indicated the single miRNA might possess wide-ranging functions and involve in different kinds of signal pathways in potato. The products of the target genes include functional protein, transcriptional factors and enzyme, etc.

The putative function of differentially expressed miRNAs

Some differential expression miRNAs had been reported in other species, such as miR156 and miR169 etc., (Xiong et al., 2013, Jiao et al., 2010; Meng et al., 2011; Wu et al., 2009). In addition to these miRNAs, there are a large number of miRNA functions that have not been studied. In order to understand the regulatory functions of differentially expressed miRNAs, we carried out the miRNA target prediction analysis and annotated the target genes according to the website (http://plants.ensembl.org/Solanum_tuberosum/Search/New?db=core) (Table S7). The putative miRNA target genes participated in many biological processes such as cell wall modification, the cell cycle, photosynthesis, nutrient transport, carbon and nitrogen utilization, autophagy, ubiquitination, senescence, protein and lipid metabolism, chromatin structure, hormone metabolism, signaling and stress-related processes. This result was also consistent with differentially expressed genes related to the TOR-regulated pathways reported by previous transcriptome studies (Park et al., 2002; Payyavula, Singh & Navarre, 2013; Xiong & Sheen, 2014; Caldana et al., 2013). And then Numerous GO functional and KEGG pathway enrichment of differentially expressed miRNAs were detected, suggesting TOR regulated miRNAs participated in many biological metabolic processes in potato (Tables S8A–S8B). These results suggest that miRNA may play a vital role in the regulation of metabolic processes in TOR signaling.

Integrative miRNA-mRNA expression and function analysis

In our previous study, our focus is on the relationship between auxin, TOR and adventitious root formation (Deng et al., 2017). Our transcriptome data showed that a large number of genes associated with root development were differentially expressed in potato seedling under TOR inhibitor treatment. Loss of auxin signaling after TOR suppression results in significant down-regulation of these genes (Morin et al., 2008; Deng et al., 2017). In addition, a large number of genes related to biosynthesis and metabolism are differentially expressed. How TOR affects the expression of these genes is still unclear. To further determine the regulatory function of TOR at post-transcriptional level, we performed miRNA-mRNA reverse correlation analysis to identify reliable miRNA and their targets. In present study, thirteen up-regulated miRNAs and 14 down-regulated miRNAs were corresponding to 47 down-regulated 82 up-regulated mRNAs, respectively (Tables S9A–S9D). GO functional (Table S10A–S10B) and KEGG pathway (Fig. 3) analysis showed that the negatively correlated miRNA/mRNA interaction pairs participate in some key metabolic processes such as lipid synthesis, amino acid metabolism, nitrogen metabolism, flavonoid biosynthesis and so on. To validate the results of high-throughput sRNA sequencing, the expression patterns of 8 miRNAs were analyzed by qRT-PCR. Meanwhile, to detect the expression patterns of the targets identified by transcriptome analysis, 6 target genes were selected for qRT-PCR analysis. The results showed that all the selected miRNAs and targets shared similar expression tendency with the original results (Fig. 4).

Differentially expressed miRNAs were involved in some anabolic and biosynthetic pathways

TOR, as a central regulator of anabolism, participates in many biosynthetic processes, such as ribosome biosynthesis, and protein translation and synthesis. In the integrative miRNA-mRNA data, DE miRNAs target genes also participate in many anabolism processes. As shown in Fig. 3 and Table 3, differentially expressed miRNAs were involved in phenylpropanoid biosynthesis, cutin, suberine and wax biosynthesis. On the other hand, TOR also negatively regulates the synthesis of some metabolites, such as glucosinolates and flavonoids. Suppression of TOR could up-regulated genes responsible glucosinolates and flavonoids synthesis (Caldana et al., 2013). It was found that the differentially expressed miRNAs target genes were also involved in the biosynthesis of flavonoids, suggesting that TOR may regulate the biosynthesis of flavonoids relying on miRNAs. Interestingly, TOR, as the core regulatory element of ribosomal biosynthesis, large number genes of ribosome was differentially expressed in transcriptomics. But only few differentially expressed miRNAs target genes were involved in ribosomal biosynthesis. These results indicated that miRNAs are only part of TOR regulation networks in transcription and post-transcriptional level.

Differentially expressed miRNAs related to some catabolic processes

Catabolic pathways is another set of metabolism that breaks down large molecules (such as lipids, polysaccharides, proteins and nucleic acids) into smaller unites (such as fatty acids, monosaccharides, amino acids and nucleotides, respectively) which are either oxidized to release energy or used in other anabolic reactions. Autophagy is a natural regulated project of the cell that disassembles unnecessary or dysfunctional components for turnover and recycling of intracellular macro molecules and whole organelles. In the mammals, many miRNAs (such as miR-145, miR-181, miR-93-5p, miR-33 and miR-21) (Zhang et al., 2017b) as the upstream regulator of TOR signaling pathway were shown to modulate disease biogenesis through inhibiting autophagy. In this study, the KEGG pathway “regulation of autophagy” was affected under the TOR suppression, and a total of three down-regulated miRNAs as the downstream of the TOR signaling pathway were targeting the autophagy-related genes (Table 3). Ubiquitination is a key posttranslational modification carried out by a set of three enzymes: ubiquitin-activating enzyme (E1), ubiquitin-conjugating enzymes (E2), and ubiquitin-protein ligase (E3). The “ubiquitin mediated proteolysis” KEGG pathway was one of the most enriched pathways in the miRNA data. A total of four differentially expressed ubiquitination-related miRNAs were observed and three were significantly down-regulated (Table 3). These observations strongly support the previous results in which TOR negatively regulated the catabolic process (autophagy and ubiquitination) in yeast, mammals and plants (Xiong & Sheen, 2014; Dong et al., 2015). In previously studies, it has showed that suppression of TOR could lead to accumulation of sugars and lipids in Arabidopsis (Caldana et al., 2013; Dobrenel et al., 2013). In this study, we found that many miRNA targeted genes were involved in these metabolic processes. TOR may influence the metabolic processes of sugars and lipids by regulating corresponding miRNAs. Thus, our results showed that miRNAs may act as a bridge between TOR and catabolism.