A comparative analysis of small RNA sequencing data in tubers of purple potato and its red mutant reveals small RNA regulation in anthocyanin biosynthesis

Plant materials

SD92, commonly known as Hei Jingang, was a tetraploid potato with purple skin and purple flesh. SD140 is a mutant of SD92. Its skin and flesh colors were red (Liu et al., 2018; Liu et al., 2015). Two materials were planted in a greenhouse for two months at 20  ± 2 °C with a photoperiod of 16 h light/8 h dark.

Sample library construction and sequencing

Fresh tubers (diameter 4–5 cm) from three individual plants were harvested for three biological replicates, cleaned with sterilized water, frozen in liquid nitrogen and finally stored at −80 °C. Total RNA extraction of the samples was performed with a modified Trizol reagent (Liu et al., 2018) for library construction and validation of miRNA sequencing data.

Small RNA was isolated and the library was constructed in accordance with the protocol of Preparing Samples for Analysis of Small RNA (Illumina, San Diego, CA, USA). The 18-30 nt RNA segments were separated from total RNA by polyacrylamide gel electrophoresis, then ligated with 3′ adaptor (GAACGACATGGCTACGATCCGACTT) and 5′ adaptor (AGTCGGAGGCCAAGCGGTCTTAGGAAGACAA). The resulting segments were employed to synthesize first-strand cDNA. The cDNA was amplified and only cDNA with both 3′ and 5′ adaptors was enriched. Finally, the fragments of 100–120 bp were separated to construct the library. After library quantification and single-stranded DNA cyclization, the library was sequenced by BGISEQ-500 technology. The raw data was deposited into NCBI BioProject database (PRJNA824931).

miRNA identification and prediction

The impurities of raw data, including 5′ primer contaminants, no-insert tags, oversized insertion tags, low quality tags, poly-A tags and the tags without 3′ primer, were excluded from the raw data to obtain clean tags. Low-quality tags were tags whose base quality values were less than 20, accounting for more than 50% of the total bases. The clean tags were mapped to potato reference genome PGSC_DM v4.03 (http://solanaceae.plantbiology.msu.edu/data) by Bowtie2 (Langmead et al., 2009) and small RNA databases miRBase (Kozomara & Griffiths-Jones, 2014), snoRNA (Yoshihama, Nakao & Kenmochi, 2013) and Rfam (Nawrocki et al., 2015). If a small RNA could be mapped to more than one database, the small RNA annotation followed the searching priority of miRBase > snoRNA > Rfam. One small RNA was only mapped to one category. The small RNAs mapped to Rfam database were validated by cmsearch (Nawrocki & Eddy, 2013). The novel miRNA was determined by miRA (Evers et al., 2015) according to the characteristic hairpin structure of miRNA precursor. Small interfering RNA (siRNA), a 22–24 nt double-strand RNA, was identified by the characteristic of one strand 2 nt shorter than the other (Jagla et al., 2005).

miRNA expression and screening of differentially expressed miRNAs (DEMs)

The expression level of miRNA was estimated by the transcripts per kilobase million (TPM) (’t Hoen et al., 2008). The differential expression was calculated by DEGseq (Wang et al., 2010) based on MA-plot method (Yang et al., 2002). The P-value calculated for each gene was adjusted to Q-value for multiple testing corrections by two alternative strategies. The miRNAs with expression fold change > 2 and Q-value < 0.001 were defined as differentially expressed miRNAs. The volcano plot and heatmap of differentially expressed miRNAs were obtained by Excel 2016 and MeV (Saeed et al., 2003), respectively.

Target gene prediction, Gene Ontology (GO) and KEGG pathway enrichment analyses

TargetFinder (Fahlgren & Carrington, 2010) and psRobot (Wu et al., 2012) were used to predict the target genes of miRNAs. All target genes were mapped to GO-terms in the database (http://www.geneontology.org/) and KEGG Orthology (Kanehisa et al., 2008) pathways. The enrichment analyses of GO terms and KEGG pathways were performed by the hypergeometric test based on GO::TermFinder (Boyle et al., 2004). The P-values were adjusted by Bonferroni method (Abdi, 2007). The adjusted P-value was defined as Q-value. The terms with Q-value < 0.05 were defined as significantly enriched terms.

Expression validation of miRNAs

RNAs were digested by DNaseI (Thermo, USA) to remove genome DNA. First-strand cDNA was synthesized by miRNA First Strand cDNA Synthesis Kit (Sangon Biotech, China) using tailing reaction method. Real-time quantitative PCR (RT-qPCR) was performed with UltraSYBR Mixture Kit (CWBIO, China) by using 18S rRNA (GenBank: X67238.1) as a reference gene. The primers of 18S rRNA and miRNAs were listed in Table 1. The universal reverse primer for miRNAs was supplied from miRNA First Strand cDNA Synthesis Kit. Three biological replicates were performed. Significant difference of miRNA expression between SD92 and SD140 was identified by Student’s t-test (P < 0.05).

Sequencing and classification of potato small RNAs

To identify the miRNAs regulating potato flesh color, six small RNA libraries were constructed and sequenced. The counts of raw tags of six libraries ranged from 28,058,311 to 30,152,601 (Table 2). Low quality tags, invalid adapter tags, poly-A tags and short valid length tags (shorter than 18 nt) were removed to obtain clean tags. The percentages of clean tags of six libraries ranged from 92.10% to 95.22%, which indicated the sequencing data could be used for subsequent analyses. Of the six libraries, 19-25 nt length tags accounted for 87.9%–96.4% of the total tags, and the 24 nt tags were the most abundant (Table S1). More than 85.04% of the total clean tags from six libraries were mapped to the reference genome (Table S2). Therefore, the sequencing data should accurately reflect small RNA expression and could be used for differential expression analysis of small RNA. To classify and annotate small RNAs, the clean tags were mapped to small RNA databases miRBase, snoRNA and Rfam. The types and proportion of identified small RNAs were similar within six libraries. The intergenic miRNAs were the most abundant (Table S3).

Identification of known and novel miRNAs

There were about 300 known miRNAs and 160 novel miRNAs identified in every library (Table 3). In total, 356 known miRNAs belonging to 121 miRNA families were identified (Table S4), and miR172 family was the most abundant family where 21 members were identified. The nucleotide bias analyses on these non-redundant known miRNAs (Fig. S1A) showed that the first and 24th nucleotide preferred to be uracil (U), and adenine (A) was the dominant nucleotide in the 10th nucleotide position. Meanwhile, several nucleotide positions were conserved. For example, the proportions of four kinds of nucleotides were nearly equal in the 4th, 9th and 16th nucleotide position (Fig. S1A).

Unmapped tags were further used to predict novel small RNAs. Totally, 171 novel miRNAs were identified for six libraries. The mature sequences, star sequences and precursor sequences of 171 novel miRNAs were listed in Table S5. The length of the novel miRNAs ranged from 19 to 30 nucleotides. Most of the nucleotide positions preferred to be uracil (U) or adenine (A) (Fig. S1B). Two exceptions were the 9th and 11th nucleotide where the dominant nucleotides were guanine (G) and cytosine (C), respectively.

Differentially expressed miRNAs between SD92 and SD140

To further validate the expression changes of miRNAs between SD92 and SD140, 15 miRNAs from 11 different miRNA families were randomly selected to be tested by RT-qPCR (Fig. 1). The results of RT-qPCR showed the same expression regulation pattern with miRNA sequencing data, which suggested that the miRNA sequencing result was reliable. What’s more, the results showed 6 miRNAs were differentially expressed between SD92 and SD140 (P < 0.05). Different miRNAs from the same miRNA family displayed the same regulation pattern. For example, both miR166a-3p and miR166d-5p_2 were from miR166 family and exhibited higher expression levels in SD140 than in SD92.

A total of 179 differentially expressed miRNAs were identified in this study, including 107 known miRNAs and 72 novel miRNAs (Fig. 2A, Table S6). Among the differentially expressed miRNAs, 65 and 114 were confirmed to be up- and down-regulated in SD140, respectively. These miRNAs belonged to 49 miRNA families. Of the 49 miRNA families, miR399 and miR172 family were the two largest families, which contained 10 and 9 miRNA members, respectively. Interestingly, the members of miR399 and miR172 families were significantly down-regulated in SD140, respectively.

Target gene prediction of miRNAs

To further explore the function of miRNAs, the target genes (mRNAs) of all miRNAs were predicted by psRobot and TargetFinder. Totally, 7,416 target genes were identified for 450 miRNAs where 897 target genes were confirmed as targets of 116 miRNAs by both softwares. Among these 897 target genes, 305 genes were regulated by 31 differentially expressed miRNAs (Table S7).

GO and KEGG pathway enrichment analysis of target genes

GO enrichment analysis of the above 305 target genes showed that the biological process ontology included 47 GO terms. “Cellular macromolecule metabolic process” and “macromolecule metabolic process” were the most abundant GO terms, containing 77 genes, respectively.

The cellular component ontology included 16 GO terms, and the most abundant terms were “cell” and “cell part”, which contained 115 genes, respectively. The molecular function ontology included 10 GO terms. The GO term “binding” contained 126 genes, which was the most abundant term in molecular function (Fig. 3).

To explore the possible function of target genes, KEGG pathway enrichment analysis was performed. The 305 target genes of 31 DEMs were distributed in 6 first-level and 33 second-level KEGG pathways, respectively. The first-level KEGG pathway term “metabolism” was the most abundant, including 10 second-level KEGG pathway terms. Thirty-eight target genes were assigned in the second-level KEGG pathway term “signal transduction”, which was the most abundant second-level KEGG pathway term (Fig. 4).

Among the enriched top 20 pathways, only two pathways, “plant hormone signal transduction” and “plant-pathogen interaction”, were defined as significantly enriched pathways (P < 0.05), which comprised 24 target genes, respectively (Fig. 5 and Table S8). This indicated that the DEMs between SD92 and SD140 might be involved in plant-pathogen interaction and hormone signal transduction.

Target genes of miRNAs involved in regulation of anthocyanin biosynthesis

Generally, plant miRNAs regulate target mRNAs through two major mechanisms, transcript cleavage and translational inhibition (Chen, 2009), thus there are negative regulation relationship in the expressions of miRNA and corresponding target genes. In our previous study, a comparative transcriptome analysis was performed between SD92 and SD140 (Liu et al., 2018). By combining transcriptome sequencing data (SRA accession number: SRP125987) and miRNA sequencing data of present study, 31 differentially expressed miRNAs and corresponding target mRNAs were identified and listed in Table S9. Among them, the differentially expressed miRNAs negatively regulating target mRNAs were screened, and 140 miRNA-mRNA pairs were confirmed. In these miRNAs-mRNAs pairs, miRNAs contained 5 known miRNA families and 12 novel miRNAs. These mRNAs corresponded to 71 genes (Table 4). These genes mainly encoded transcription factors, quamosa promoter binding protein, hormone response factors, protein kinases and disease resistance protein.

Transcription factors affect anthocyanin biosynthesis by regulating the expressions of structural genes (D’Amelia et al., 2014; Liu et al., 2016). In this study, we focused on the regulation of miRNA on transcription factors in anthocyanin biosynthesis (Table 4). PGSC0003DMG400006604, PGSC0003DMG400011046 and PGSC0003DMG400012038, which were regulated by miR172b, encoded AP2 transcription factor SlAP2e, RAP2-7-like and RAP2-7, respectively. The target gene of miR530b_4, PGSC0003DMG400025479, encoded AP2-like transcription factor TOE3. PGSC0003DMG400011457 encoded WRKY transcription factor 48 and was regulated by miR172e-5p. Both PGSC0003DMG400004826 and PGSC0003DMG400018279, which were regulated by novel_mir170, encoded transcription factor ERF039-like and MYB35-like, respectively.

Hormones improve the biosynthesis of anthocyanins (Zhang et al., 2011; Palma-Silva et al., 2016), so we did research on miRNA regulating hormones in this experiment in order to throw light on miRNA regulation mechanism on anthocyanins biosynthesis. In this study, RAP2-7 and RAP2-7-like, which were regulated by miR172b, were ethylene-responsive transcription factors. TOE3 transcription factor, which was regulated by miR172b and miR530b_4, was also responsive to ethylene (Table 4). The target gene of miR171b-3p, PGSC0003DMG400012683, encoded the DELLA protein that was an inhibitor of GA signal transduction.

Protein kinases are involved in anthocyanin biosynthesis (Li et al., 2016). Protein kinases regulated by miRNA were investigated in this study. Both PGSC0003DMG400018811 and PGSC0003DMG400024795, which were regulated by novel_mir170, encoded LRR receptor-like serine/threonine protein kinase ERECTA and RCH1, respectively. PGSC0003DMG400026383 encoded receptor-like protein kinase and was regulated by novel_mir117.

There were also significant changes in the expression of target genes regulated by other miRNAs, such as PGSC0003DMG402007414, which was target gene of novel_mir105 and novel_mir143, but the gene function was unknown.