Identification of microRNAs and target genes in the fruit and shoot tip of Lycium chinense: a traditional Chinese medicinal plant.

Although Lycium chinense (goji berry) is an important traditional Chinese medicinal plant, little genome information is available for this plant, particularly at the small-RNA level. Recent findings indicate that the evolutionary role of miRNAs is very important for a better understanding of gene regulation in different plant species. To elucidate small RNAs and their potential target genes in fruit and shoot tissues, high-throughput RNA sequencing technology was used followed by qRT-PCR and RLM 5'-RACE experiments. A total of 60 conserved miRNAs belonging to 31 families and 30 putative novel miRNAs were identified. A total of 62 significantly differentially expressed miRNAs were identified, of which 15 (14 known and 1 novel) were shoot-specific, and 12 (7 known and 5 novel) were fruit-specific. Additionally, 28 differentially expressed miRNAs were recorded as up-regulated in fruit tissues. The predicted potential targets were involved in a wide range of metabolic and regulatory pathways. GO (Gene Ontology) enrichment analysis and the KEGG (Kyoto Encyclopedia of Genes and Genomes) database revealed that "metabolic pathways" is the most significant pathway with respect to the rich factor and gene numbers. Moreover, five miRNAs were related to fruit maturation, lycopene biosynthesis and signaling pathways, which might be important for the further study of fruit molecular biology. This study is the first, to detect known and novel miRNAs, and their potential targets, of L. chinense. The data and findings that are presented here might be a good source for the functional genomic study of medicinal plants and for understanding the links among diversified biological pathways.

Introduction

Lycium chinense (Mill.), popularly known as Chinese boxthorn or goji berry and originating from China, is widely distributed in the temperate and subtropical regions of Japan, Korea and southeastern parts of Asia and is even found in European countries [1]. This species belongs to the nightshade (Solanaceae) family. The dry root bark and fruit of this plant have been used in traditional medicine, and young leaves have been widely consumed as a new functional vegetable. Numerous biochemicals, nutrients, minerals and vitamins in remarkable quantities have been reported [2]. The presence of protein, fat, carbohydrates, crude fiber, minerals, calcium, phosphorous, iron, carotene, thiamine, riboflavin, nicotinic acid and ascorbic acid in the dry fruits have ranked this species as one of the highest valued functional foods [1]. Although the mechanism is unknown, according to Chinese traditional medicine, this plant improves visual acuity [3]. The dry fruits of Lycium species also possess a wide range of clinical importance, such as anti-inflammatory effects [4], anti-aging qualities [5], anticancer and immunomodulating activities [6, 7], and blood glucose- and lipid level-reducing capabilities [8]. Although this species has many fascinating aspects, there is limited genome and molecular information. The gene regulation and functions of miRNAs need to be elucidated in this species.

A group of RNAs with a low molecular weight, termed small RNAs (sRNAs), possess tremendous regulatory functions. Small RNAs are classified into the following two categories depending on their biogenesis and action: short-interfering RNAs (siRNAs) and microRNAs (miRNAs) [9, 10]. Among these sRNAs, miRNAs are non-coding entities that are endogenous, 20–24 nt long and widely found in eukaryotes [11, 12]. These miRNAs regulate gene expression at the transcriptional and post-transcriptional levels [13], and play a crucial role in plant growth and development [10, 14, 15]. MicroRNAs are processed from the stem-loop position of long primary transcripts by a Dicer-like (DCL) enzyme in plants, and loaded into silencing complexes, where they cleave the complementary sequence of targets [16]. Plant miRNAs have only been known for approximately a decade [17], but are already known to play pivotal roles at every major developmental stage, especially in core gene regulatory networks, and are also known to regulate transcription factors (TF). Therefore, miRNAs are considered as “master regulators” [12, 18].

In fact, the large number of miRNAs that are evolutionarily conserved from mosses and ferns to higher flowering plants have been used as helpful indicators for the prediction of miRNAs using homology searches across various species [19, 20]. Recently, next-generation high-throughput sequencing technologies have permitted the efficient identification and quantification of miRNAs with high precision, which dramatically facilitates the understanding of miRNAs. Moreover, miRBase (The miRBase database is a searchable database of published miRNA sequences and annotation, http://mirbase.org/) has released a new version (v20, June 2013) with 24,521 miRNA loci, leading to the production of 30,424 mature miRNAs products belonging to 206 species [21]. Among the plant species, the majority of miRNAs have been identified in the monocot Oryza sativa and the eudicots Populus trichocarpa and Arabidopsis thaliana. After deep sequencing, miR167, miR1857 and miR172a were identified in sweet orange and have regulatory effects in lycopene accumulation [22]. There are also reports elucidating the molecular biology of different miRNAs related to fruit ripening, softening, ethylene biosynthesis and the signal transduction of different phytohormones [23-26].

Although miRNAs have been thoroughly studied in many species, no miRNAs have been identified so far from goji, and none are listed in the miRBase. Many studies have confirmed that both species-specific and conserved miRNAs are very important for plants in different biological processes. Therefore, profiling miRNAs in a non-model plant species is very important for understanding the evolutionary trend of miRNAs and the regulation of different biological processes. In this study, high-throughput sequencing was performed for the RNAs of young shoot tips with adjacent leaves and the bulk mixing of three developmental stages of fruits of L. chinense. We depict an overall scenario of miRNAs of the vegetative tissues (shoot and leaf) and reproductive tissues (three developmental stages of fruit) of this species, and then target genes that are related to fruit development as identified through significant differentially expressed miRNAs. Comparative studies and expression differences have revealed an interesting pattern and distribution of miRNAs between vegetative (shoot and leaves) and reproductive (fruit) tissues. Therefore, the miRNAs and targeted genes that were identified in this study are unique resources for biological pathway analysis and functional characterization. These data could be an exclusive source of comprehensive understanding of miRNA expression profiling in L. chinense, as it is an important representative of medicinal plants.

Materials and Methods

Plant materials

Healthy and non-lignified shoot tips (approximately 1 cm from the tip end), including adjacent leaves of L. chinense (cv. Large leaf goji), were collected from two-month-old plants that were grown at the Wuhan Botanical Garden of Chinese Academy of Sciences (CAS), Wuhan, China. To collect the optimum stages of fruits, three developmental stages were chosen, green, color breaking and red mature (S2 Fig). The shoot and leaves were mixed together for total RNA extraction, but individual stages of fruit were extracted separately followed by equal bulk mixing for sequencing and subsequent qRT-PCR and RLM 5’-RACE experiments. The tissues collected in sealer bags and snap-frozen immediately in liquid nitrogen. The target tissues were stored in a -70°C ultra-refrigerator for later and/or further uses.

Total RNA isolation, purity check, sRNA library construction, and Illumina sequencing

The total RNAs were isolated from the tissue samples via the Trizol method with a slight modification following the manufacturer’s instructions of the RNAiso Plus (Takara) extraction kit. RNA degradation and contamination were monitored on 1% agarose gels, and purity was checked using the Nano photometer (Implen, CA, USA). The RNA concentration was measured with a Qubit 2.0 Fluorometer (Life Technologies, CA, USA), and the integrity was assessed by Agilent bioanalyzer 2100 (Agilent Technologies, CA, USA). The purified RNAs were used to construct small RNA libraries and were sequenced on an Illumina HiSeq 2500/2000 platform at the Novogene Company, Beijing, China.

Bioinformatics analyses of sequencing data

Raw sequence data were obtained via the RNA high-throughput sequencing process. Raw reads were screened out to remove the contaminating reads, sequences containing ‘adapters’, without insert tags, and reads with poly-A tails. Sequences from 18–40 nt in length were used for further analysis. Clean reads of sRNA tags were mapped to the L. barbarum transcriptome database with less than 2 bp mismatch (provided by South China Botanical Garden, CAS, Guangzhou, China) using Bowtie [60]. Tags matching non-coding RNAs (rRNAs, tRNAs, snRNAs, and snoRNAs) were BLAST-searched with the sequences of Rfam (Rfam: http://www.sanger.ac.uk/software/Rfam) and the NCBI GenBank databases (GenBank: http://www.ncbi.nlm.nih.gov/blast/Blast.cgi) and deducted for downstream analyses [61].

The remaining sRNA tags were used to search the conserved or known miRNAs in miRBase 20.0 allowing a maximum of two mismatches [62]. The modified software miRDeep2 [63] and srna-tools-cli were used to obtain the potential miRNA and draw the secondary structures. The characteristics regarding the hairpin structure of the miRNA precursor were used to predict the potential novel miRNAs. With the integrated use of the software miREvo [64] and miRDeep2, we predicted the potentially novel miRNAs by exploring the secondary structures, the DL1 cleavage sites and the minimum free energy (MFE) of the sRNA tags that were unannotated in the former steps.

To explore the occurrence of miRNA families, we used the known miRNAs in miFam.dat (http://www.mirbase.org/ftp.shtml), and the novel miRNA precursors were aligned to Rfam (http://www.rfam.sanger.ac.uk/search/). To predict the target genes of miRNAs, psRobot.tar was used in psRobot software [65].

After predicting the target genes, the expression levels of miRNAs (quantification of miRNA) were estimated using TPM (transcripts per million) according to the normalization formula that was described by Zhou et al., 2010 [66]. The DEGseq R package (to minimize the positive false discovery rate for the sequencing data without biological replicates) was used (2010) to identify differentially expressed miRNAs between the two samples. The corresponding P-values were adjusted to a q-value [67]. The threshold level of the P-value was set as P-value<0.01 and ׀log2(foldchange)׀>1 for the significant differential expression.

The dataset that was studied in this article is available in the NCBI (SRA) public repository under the accession number SRP043345.

Validation of miRNA expression using qRT-PCR

Poly (A)-tailed qRT-PCR was used to measure the gene expression variation and validate the deep sequencing results of miRNAs. The reverse-transcription reaction and real-time PCR primer design were conducted according to Rui Shi [68] and the manufacturer’s instructions (Takara-Prime script miRNA qPCR starter kit, 2.0), where each PCR reaction was performed in a volume of 25 µl containing 12.5 µl of SYBR premix ex Taq II (2×), 1 µl of each forward primer, 1 of µl universal reverse primer and reverse-transcribed cDNA from ~100 pg of total RNA. The PCR protocol was 5 seconds at 95°C, 40 cycles of 95°C for 5 seconds, 60°C for 20 seconds and 72°C for 1 minute. The primers that were used in this study are listed in Additional file 7. Real- time PCR was performed on the Applied Biosystem 7500 detection system using the SYBR Green I method, and all of the reactions were run in biological triplicates. The melting curve was used to determine the specificity of PCR products (primer amplicons).

To validate the reliability of Illumina throughput sequencing technology, 14 miRNAs were selected, 7 of which were randomly selected from known miRNAs and 7 from novel miRNAs after determining the appropriate reference genes. The delta-CT (corresponding cycle threshold) method was used to calculate the relative expressional levels of miRNAs [69]. In this method, the CT values of the desired miRNAs and reference genes were first transformed for measurement using delta-CT followed by dividing the quantities of desired miRNAs by the geometric mean of the reference genes. The standard deviation and mean were calculated using triplicate qRT-PCR assays. The CT was calculated using the machine accessory software and converted into relative copy numbers using a standard curve as previously described [30]. The gene 5.8 S rRNA was used as reference gene in the qPCR detection of miRNAs. Student’s t-test was used for the statistical analysis of the qRT-PCR.

Functional assignments for potential target genes

There are three distinct segment of GO: biological process (BP), cellular component (CC) and molecular function (MF). For the GO enrichment analysis, the over-represented P-value and the corrected P-value identified the highly significant gene categories in the GO enrichment analysis followed by KEGG database searching to determine the functional diversity of the ‘target gene candidates’ of differentially expressed miRNAs [70]. Moreover, KOBAS [71] software was used to test the ‘target gene candidates’ in the KEGG pathways.

Validation of predicted targets using RLM 5’-RACE

For the validation of candidate targets, modified RNA Ligase-Mediated 5’-RACE (RLM 5’-RACE) was performed. 5’-Full RACE Kit (Takara) was used according to the manufacturer’s instructions, with slight modifications. Briefly, total RNA was directly ligated to the 5’ adaptor followed by reverse transcription with the oligo (dT) primers. PCR was performed with 5’ primers and 3’ gene-specific primers using the cDNA as the template (S7 Table). The 5’-RACE PCR products were purified using Takara PCR product recovery kit, cloned, and sequenced.

Results

Summary of small-RNA library dataset by deep sequencing in L. chinense

To explore the regulatory networks of miRNAs in the fruit and shoot of L. chinense, we constructed two libraries from the young shoot tip with leaves and fruit tissues. To identify the target genes and regulatory pathways, Illumina sequencing technology was used. A total of 19,014,652 raw reads were obtained, of which 9,128,886 and 9,885,766 raw reads were produced in shoot and fruit sRNA libraries, respectively (Table 1). Moreover, 97.27% and 97.62% raw reads of high quality for the shoot and fruit libraries, respectively. For further study, 18 to 40 nt total small RNA sequences were selected (Fig. 1), with 21–24 nt long sequences with high frequency. The frequency of shoot sRNA was 40.19% (24 nt), followed by 10.14% (23 nt), 8.63% (21 nt) and 7.92% (22 nt), while this frequency was 40.33% (24 nt) followed by 9.83% (21 nt), 9.11% (23 nt) and 8.09% (22 nt) in the fruit sRNA library. The percentage of small RNAs of 24 nt was significantly higher than that of other sRNAs.

Table 1

Statistical summary of the data that were generated by high-throughput small-RNA sequencing in goji.

Read type	sRNA libraries		Total
	Shoot	Fruit
Total raw reads	9,128,886 (100%)	9,885,766 (100%)	19,014,652
High quality	8,879,712 (97.27%)	9,650,746 (97.62%)	18,530,458
Low quality reads	3007 (0.03%)	6160 (0.06%)	9167
N%> 10%	3612 (0.04%)	0	3612
5’ adapter contaminants	3124 (0.03%)	6104 (0.06%)	9228
3’ adapter null	187,402 (2.05%)	189,023 (1.91%)	376,425
Poly (A/T/G/C)	52,029 (0.57%)	33,733 (0.34%)	85,762
Total bases (G)	0.456	0.494	0.95

N%, undetectable gaps at base calling.

Figure 1

Length distribution of small RNAs in the shoot and fruit libraries of L. chinense.

X-axis, length of sRNA distribution; Y-axis, corresponding percentage of raw reads.

Based on a trans-species alignment, the small RNAs were mapped against the Goji berry (Lycium barbarum L.) transcriptome database as provided by the South China Botanical Garden (SCBG, CAS, Guangzhou, China). A total of 3,826,551 (44%) and 3,278,589 (43.94%) unique small RNA reads were mapped for the shoot and fruit libraries, respectively. These mapped small RNAs comprised known miRNAs, putative novel miRNAs, rRNAs, tRNAs, snRNAs, snoRNAs, repeat-associated RNAs, TAS (trans-acting small interfering RNAs) and unannotated fragments (Table 2). Additionally, unannotated reads accounted for 81.70% and 82.76% in the shoot and fruit libraries, respectively.

Table 2

Number of reads for each small RNA classification as identified in goji.

Read type	sRNA libraries		Total
	Shoot	Fruit
Total	3,826,551 (100%)	3,278,589 (100%)	7,105,140
Known miRNAs	22,966 (0.60%)	31,058 (0.95%)	54,024
Novel miRNAs	12,519 (0.33%)	28,187 (0.86%)	40,706
rRNAs	566,742 (14.81%)	425,905 (12.99%)	992,647
tRNAs	20,599 (0.54%)	11,522 (0.35%)	32,121
snRNAs	2886 (0.08%)	1715 (0.05%)	4601
snoRNAs	11,651 (0.30%)	3832 (0.12%)	15,483
Repeat associated RNAs	51,953 (1.36%)	52,636 (1.61%)	104,589
TAS	10,907 (0.29%)	10,243 (0.31%)	21,150
Unannotated	3,126,328 (81.70%)	2,713,491 (82.76%)	5,839,819

Known miRNAs, perfect matching mature miRNA sequences in the miRBase miRNA repository; Novel miRNAs, new mature miRNA sequences that were detected by characteristic stem-loop hairpins and minimum free energy (MFE); rRNAs, tRNAs, snRNAs, and snoRNAs are non-coding sRNAs; TAS, transacting small interfering RNAs; Unannotated, reads that had no known information.

Identification of conserved and evolutionarily known miRNAs

To identify the known miRNAs, the sRNA libraries were BLASTN-searched for known mature plant miRNAs that were deposited in miRBase 20.0 (v20, June 2013). Following BLASTN search and subsequent sequence analyses, 94 known hairpin structures were identified corresponding to 60 unique mature miRNA sequences that were orthologs of conserved miRNAs from other plant species in the shoot and fruit libraries of L. chinense, representing 31 diversified miRNA families (S1 Table). In addition, among these 60 known miRNAs, 53 and 46 miRNAs were expressed in the shoot and fruit libraries, respectively.

The member numbers of different miRNA families were also analyzed (Fig. 2), which showed that approximately half of the families contained more than one member, while the other half had a single member. The families, MIR162, MIR166, MIR167, MIR168, MIR171, MIR319, MIR397, MIR398, MIR482, MIR6024, MIR6025 and MIR7997 possessed 2–3 members, while MIR156, MIR172 and MIR5303 contained 5, 5 and 6 members, respectively. Sixteen families, namely MIR159, MIR160, MIR164, MIR1918, MIR1919, MIR394, MIR399, MIR530, MIR5300, MIR5301, MIR6020, MIR6022, MIR8007, MIR8011, MIR8021 and MIR8031, contained only one member. Gao et al. [27] has reported that miRNA families MIR156, MIR159, MIR160, MIR166, MIR167, miRNA171 and MIR398, are ancient, as their orthologs are distributed in Coniferophyta and Embryophyta. On the contrary, six families (MIR162, MIR164, MIR168, MIR172, MIR394 and MIR399) had homologs in angiosperms, idicating that these L. chinense miRNA families might have recently evolved.

Figure 2

The distribution of the member numbers (A) and read counts of each known miRNA families (B) in L. chinense.

X-axis, name of the miRNAs family; Y-axis, representing miRNA numbers.

Prediction of potentially novel miRNAs and nucleotide bias

The software programs miREvo and miRDeep2 were used to predict potentially novel miRNAs from the unannotated small-RNA tags by exploring the secondary structure and dicer cleavage site and measuring the minimum free energy (MFE). Custom scripts were used to obtain the identified miRNA counts and base bias on the first position. The number of novel miRNAs was summarized along with the nucleotide bias of the first position of each small RNA of a certain length (Fig. 3). The majority of these novel miRNAs had a length of 22 nt in both of the libraries, and started with a 5’ U. In contrast, for the known miRNAs, the majority had 21 nt in both of the libraries with a 5’ U. A total of 30 novel mature miRNA sequences were obtained, of which 25 and 29 were shoot- and fruit-specific, respectively. In this study, the miRNA precursor length of the putative novel miRNAs ranged from 46 to 295 nt.

Figure 3

First nucleotide bias of novel miRNA candidates in the L. chinense shoot (A) and fruit (B) libraries.

The number on top of the bars indicates the number of sequences corresponding to the miRNA length; “0” indicates the absence of sequences of that corresponding length. X-axis, length of miRNAs; Y-axis, representing percent of each nucleotide.

Differentially expressed known and novel miRNAs between the shoot and fruit

To identify miRNAs regulating fruit development, the expression of each miRNA was normalized to transcripts per million (TPM) and compared between the shoot and fruit. Of 60 known mature miRNAs, 14 were only found in the shoot, and 7 were only found in the fruit, with 39 being commonly expressed in both of the tissues. The number of expressed miRNAs varied between the shoot and fruit, with 53 known miRNAs expressing in the shoot and 46 known miRNAs in the fruit (S1 Table). However, this pattern was opposite in the case of novel miRNAs, with 29 novel miRNAs being expressed in the fruit and 25 in the shoot, with 5 fruit-specific and only 1 shoot-specific novel miRNA (Table 3). Moreover, the expression level of miRNAs varied greatly, i.e., the abundance ranged from 24.55 TPM to 63,523 TPM in the shoot and from 27.1 TPM to 693,556 TPM in the fruit. In shoot library, miR166b, miR5301 and miR162 were the most represented families with abundances of 63,523 TPM, 18,269 TPM and 14,610 TPM, respectively. These three miRNA families were also found to be highly represented in fruit, but the top three were miR167, miR166b and miR162 with 693,356 TPM, 160,252 TPM and 75,796 TPM, respectively. Interestingly miR166b was expressed most frequently, ranking first in the shoot and second in the fruit, covering nearly 67% of the differentially expressed known miRNA transcripts in the shoot and fruit libraries (S1 Table). A previous report has identified the rice homeobox-leucine zipper protein HOX32 as the target of miR166. These types of transcription factors are involved in various tissues or different plant developmental stages [28].

Table 3

Description of 30 potentially mature novel miRNAs.

miRNA	ML (nt)	Sequence	PL (nt)	miRNA reads		Ratio (Shoot/Fruit)	miRNA*
				Shoot	Fruit
LC1	22	uagggcguucggauccuucugc	125	2887	18,833	0.29	Yes
LC11	21	aauccuucugcaauccauaac	144	269	405	1.27	Yes
LC14	21	uagaaagaguuuguaggcgag	286	56	438	0.24	Yes
LC16	22	ucuuaccaauaccucccauucc	94	346	182	3.64	Yes
LC17	22	uccaaucuccucgcccauauuu	76	0	95	-	Yes
LC2	22	uugccaauuccccccauuccga	76	2662	468	10.88	Yes
LC20	21	cgggguauuguaaguggcaga	92	56	84	1.28	No
LC22	24	acaagcucugauaccauguaggaa	104	102	73	2.67	Yes
LC23	21	uccaaauaaaccuccaucgug	82	3	39	0.15	No
LC24	21	uuguugugugauaaauccagg	129	9	9	1.91	No
LC25	21	cauaggauucuugggcaugcu	240	0	18	-	Yes
LC26	22	ucucuuucggucuagaucgucu	295	4	21	0.36	No
LC27	24	acaggcucugauaccaugaaggaa	102	25	15	3.19	No
LC28	21	uuauuaucucagaaagucacu	80	3	0	-	Yes
LC29	21	uuucugguguauagugauaau	93	10	3	6.38	No
LC3	22	uuaaaaaggguuauguaguggc	279	0	1788	-	No
LC30	21	cagcagcgguaggaaugaagc	81	0	6	-	Yes
LC31	21	caaaggccacaagauucacuu	148	5	1	9.56	Yes
LC32	21	auaauacuuggaauaugcccu	222	330	127	4.97	Yes
LC33	24	aaucccgggauuguaguguuauuu	46	11	16	1.32	No
LC34	24	auguugcucggacucuuugaaaau	107	3	5	1.15	No
LC35	21	aaacccucagcgauccauaac	144	0	3	-	Yes
LC36	21	uuugcgugaacacuacagagu	66	17	18	1.81	No
LC37	21	ugucgcaggugacuuucgccc	85	673	404	3.19	No
LC4	20	gaucaugugguagcuucacc	73	85	704	0.23	No
LC40	21	agaaagcaagguuucaggugu	235	8	3	5.10	No
LC41	21	uuaaggcguguagaugugcau	143	245	41	11.43	No
LC46	21	ucggacucuucaaaaauguug	107	1	4	0.48	No
LC6	21	ugggucgcugaaggauugaug	172	465	357	2.49	Yes
LC9	21	uuagauucacgcacaaacuug	75	409	316	2.48	No

ML, mature sequence length; PL, Precursor length; miRNA*, the sequences that can form miRNA::miRNA* with miRNA variants.

Our high-throughput sequencing discovered 73.33% (44 of 60) conserved miRNAs and 60% (18 of 30) novel miRNAs that were significantly differentially expressed in the fruit, either as up- or down-regulating patterns. A total of 12 significantly differentially expressed miRNAs (7 known and 5 novels) were specific to the fruit-tissues, while 15 (14 known and 1 novel) were specific to the shoot tissues. These results suggested that the miRNAome, as well as the small-RNA transcriptome [29], might be very complex, although miRNA sequences are distinctly conserved across plant species [15].

Among the 44 differentially expressed known miRNAs, 17 were up-regulated, and 27 were down-regulated in the fruit significantly at the P<0.01 level. In the case of novel miRNAs, the pattern was opposite, as11 were up-regulated, and 7 were down-regulated in the fruit tissues. Of these 62 (44 known and 18 novel) differentially expressed miRNAs, 43 (34 known and 9 novel) demonstrated 3-fold or greater expression changes between the shoot and fruit tissues (Table 4).

Table 4

Significantly differentially expressed miRNAs showing 3 or more log2-fold expression changes.

sRNA	Shoot Reads (TPM)	Fruit Reads (TPM)	log2-fold change	P-value	Signature	Regulation	Sig level
miR156a	122.78	21,864.83	-7.48	0	TRUE	Down	**
miR156f-5p	24.56	20,160.03	-9.68	0	TRUE	Down	**
miR156h-3p	0	27.1	-5.76	5.28E-06	TRUE	Down	**
miR164a	491.10	19,104.67	-5.28	0	TRUE	Down	**
miR167	4468.99	69,355.9	-3.96	0	TRUE	Down	**
miR167a	2602.82	36,450.41	-3.80	0	TRUE	Down	**
miR171a	3486.80	351.79	3.31	0	TRUE	Up	**
miR171b	49.11	0	6.62	3.80E-18	TRUE	Up	**
miR171d	245.55	0	8.94	4.08E-64	TRUE	Up	**
miR172a	220.99	27.1	3.03	7.09E-75	TRUE	Up	**
miR172c-3p	98.22	0	7.62	1.81E-31	TRUE	Up	**
miR1918	0	54.12	-6.76	4.98E-11	TRUE	Down	**
miR1919–5p	24.56	0	5.62	2.10E-10	TRUE	Up	**
miR319	4223.44	270.6	3.96	0	TRUE	Up	**
miR319a	4223.44	270.60	3.96	0	TRUE	Up	**
miR319b	4223.44	270.60	3.96	0	TRUE	Up	**
miR394	6310.60	595.33	3.41	0	TRUE	Up	**
miR397a	0	27.06	-5.76	5.28E-06	TRUE	Down	**
miR397b	0	135.30	-8.08	2.30E-24	TRUE	Down	**
miR398	49.11	4329.67	-6.46	0	TRUE	Down	**
miR398a-3p	638.43	6359.20	-3.32	0	TRUE	Down	**
miR398b-3p	24.56	0	5.62	2.10E-10	TRUE	Up	**
miR530	73.67	0	7.20	4.56E-25	TRUE	Up	**
miR5300	0	108.24	-7.76	3.92E-20	TRUE	Down	**
miR5303	0	27.06	-5.76	5.28E-06	TRUE	Down	**
miR5303c	0	27.06	-5.76	5.28E-06	TRUE	Down	**
miR5303f	24.55	0	5.62	2.10E-10	TRUE	Up	**
miR6022	73.66	0	7.20	4.56E-25	TRUE	Up	**
miR6024	73.66	0	7.20	4.56E-25	TRUE	Up	**
miR6024–3p	73.66	0	7.20	4.56E-25	TRUE	Up	**
miR6025a	98.22	0	7.62	1.81E-31	TRUE	Up	**
miR6025d	98.22	0	7.62	1.81E-31	TRUE	Up	**
miR8007a-5p	24.55	0	5.62	2.10E-10	TRUE	Up	**
miR8011a-3p	24.55	0	5.62	2.10E-10	TRUE	Up	**
LC14	1375.07	11,852.47	-3.11	0	TRUE	Down	**
LC17	0	2570.74	-12.33	8.35E-256	TRUE	Down	**
LC23	73.66	1055.36	-3.84	3.22E-121	TRUE	Down	**
LC25	0	487.09	-9.93	2.36E-70	TRUE	Down	**
LC28	73.66	0	7.20	4.56E-25	TRUE	Up	**
LC3	0	48,384.06	-16.57	0	TRUE	Down	**
LC30	0	162.36	-8.34	1.94E-28	TRUE	Down	**
LC35	0	81.18	-7.34	1.05E-15	TRUE	Down	**
LC4	2087.16	19,050.55	-3.19	0	TRUE	Down	**

TPM, normalized transcripts per million; Signature-TRUE, differentially expressed miRNAs;

*, significant at P <0.05;

**, significant at P < 0.01.

Dramatic expression differences were also observed in the members of the same miRNA family. For example, miR156a, miR156f-5p and miR156h were expressed at significantly higher levels in the fruit than the shoot with 7.48-, 9.68- and 5.76-fold changes, respectively. Additionally, miR156h-3p was not expressed in the shoot tissues at all. The target genes of these miRNAs are responsible for DNA binding, and have some specific targets that are responsible for the biosynthesis of biochemicals and for cellular functions (S2 Table). Additionally, these 14 up-regulated known miRNAs with a more than 3-fold expression change varied dramatically in L. chinense fruit.

Prediction and annotation of potential target genes of differentially expressed miRNAs

The sequence complementarity between miRNAs and their corresponding target genes offers a straightforward process with which to efficiently search miRNAs target genes [30]. A total of 1139 target genes were found for 90 known and novel miRNAs. A total of 703 target genes were identified for 40 out of the 62 significant differently expressed miRNAs.

The annotation of these target genes was conducted based on GO enrichment and KEGG analyses. The number of target genes for each differentially expressed miRNAs ranged from one to more than one hundred. The highest number of target genes was predicted for miR5303 (142), followed by miR5303g (87) and miR5303 (47). However, the novel miRNA LC35 and miR5300 had only one target gene. Notably, there were 22 significant differentially expressed miRNAs that did not match with any target genes, making it impossible to understand the function of these miRNAs (S6 Table). This phenomenon could be explained by the low expression of target genes or the inefficiency of the cross-species mapping of these miRNAs with the L. barbarum transcriptome.

GO (Gene Ontology) categories were assigned for target genes, and 975 putative targets were identified for the top 30 enriched GO categories in terms of ‘biological process’, ‘cellular component’ and ‘molecular function’ (Fig. 4). Based on the biological process, the genes were classified into 20 categories of which the top three over-represented GO terms were “response to organic substances”, “endogenous stimuli” and “hormone stimuli”. The successive analysis revealed that target genes were highly diversified with respect to function. In the case of molecular function, the genes were classified into 6 groups, of which they are mostly involved in “ADP binding”, “helicase activity” and “four-way junction helicase activity”. Based on cellular components, the genes were classified into 4 categories, of which the top three were related to the “molybdopterin synthase complex”, “holiday junction helicase complex” and “DNA helicase complex”. The greatest numbers of genes were involved in “heterocyclic compound binding” (109 genes) and “organic cyclic compound binding” (109 genes) under the molecular function terms. The most over-represented GO term was “ADP binding”, involving 30 genes, followed by “response to organic substances” and “molybdopterin synthase complex”, with 7 and 3 genes, respectively. The significant differences in the categories for cell component occured in the nucleus, chloroplast and plastid, and a striking difference in the biological process was found for transcription regulation, protein modification and photosynthesis. The high frequency of the terms “protein modification” and “transcripts regulation” suggest that miRNAs are closely involved in transcription and translation [31]. Analyses of gene sets for statistically enriched GO terms have shown in S1 Fig.

Figure 4

Gene ontology terms and numbers of the predicted target genes for differentially expressed miRNAs.

BP, Biological Process; CC, Cellular Component; MF, Molecular Function. Right-hand-side scale, targeted gene numbers corresponding to the GO terms; left-hand-side scale, percent of targeted gene numbers corresponding to the GO terms.

Furthermore, a KEGG pathway analysis was performed to illuminate the biological interpretation of these target genes of differentially expressed miRNAs. A total of 75 highly diversified biochemical pathways were involved with identified miRNAs target genes, considering a corrected P-value and the number of involved genes. Moreover, the top 20 enriched pathways were discovered with 68 involved genes (Fig. 5). “Metabolic pathways” was the most significantly enriched with respect to the rich factor and gene number (19 genes), followed by “biosynthesis of secondary metabolites” (12 genes) and “spliceosome” (6 genes). S4 Table presented details of these pathways, id name, gene name and hyperlinks.

Figure 5

KEGG analysis with the 20 most enriched pathways.

The coloring of the q-values indicates the significance of the rich factor; the circle indicates the target genes that are involved, and the size is proportional to the gene number.

Expression validation of sequencing data with qRT-PCR and RLM 5’-RACE

Quantitative RT-PCR was used to experimentally validate the sequencing data. A total of 14 (7 conserved and 7 novel) miRNAs with different expression levels were randomly selected (Fig. 6). The CT values and RQ values of these miRNAs were obtained, and 12 miRNAs had a similar expression pattern in the fruit and shoot compared with the sequencing results. For example, miR156a was highly expressed in the fruit with 7.48-fold increases compared with the shoot, and a similar result was reflected in the qRT-PCR results, where approximately 7 times greater expression was found in the fruit. However, miR166b and LC14 showed acceptable inconsistency in the expression profile between the sequencing data and qRT–PCR. The possible reasons for this inconsistency might be the difference in the amplification performance of primers or other unknown reasons. The primer designing and their corresponding expression has listed in S5 Table. Additionally, four important targets were validated by modified RLM 5’-RACE with their miRNA cleavage sites (Fig. 7).

Figure 6

Expression ratios (Fruit/Shoot) of miRNAs in L. chinense based on the qRT-PCR (A) and deep sequencing (B) results.

X-axis, name of the miRNAs that were selected for qRT-PCR; column above the X-axis, miRNAs that were up-regulated in the fruit; column below the X-axis, miRNAs that were up-regulated in the shoot.

Figure 7

Target validation of known L. chinense miRNAs.

Arrows showing the 5’-ends of cleavage products.

Discussion

The exclusive application of high-throughput sequencing is rapidly replacing the traditional Sanger sequencing platforms for the identification of miRNAs. By now, the most conserved miRNAs have been identified in Arabidopsis, rice and poplar [32]. Recent studies, of non-model plants [33], have reported that the non-conserved and species-specific miRNAs are normally found at a lower abundance than conserved miRNAs, and such a low abundance is difficult to determine efficiently and cost-effectively. However, high-throughput sequencing has been successfully applied for both model plants [29, 34–36] and non-model plants [33, 37, 38], and has also been recognized as a blessing of miRNA identification [39].

Overview of the deep-sequencing datasets

Length distribution analysis is important to determine the compositional profiling of small RNAs. Normally, miRNA length is 21 nt or 22 nt, whereas, 24 nt in length is considered typical for siRNAs [14]. For both of the libraries, 24-nt small RNA sequences were the most abundant, suggesting that sRNAs of this size class play important roles in L. chinense. The same observation was previously reported on Arabidopsis thaliana [36], Citrus sinensis [40], Medicago truncatula [41] etc. However, an opposite scenario was reported in Pinus contorta, where 24-nt RNAs constitute only 2.5%, and 21 nt miRNAs are highly enriched up to 50% [27]. Additionally, Lippman et.al [50] and Lu et.al [29, 42] reported the 24-nt sRNAs were predominantly comprised of siRNAs associated with repeats and transposons. In the monocot model plant rice, 24-nt sRNAs comprises approximately 20% while 21-nt sRNAs comprises approximately 10% [37], indicating that a distinct difference exists but, may be not very severe. In this study, 1.36% of the sRNA reads were repeat-associated in the shoot and 1.61% in the fruit libraries. These repeat-associated reads were remarkably more frequent than any other class of sRNA in the two sRNA libraries (Table 2). These results also suggested that the siRNAs that are related to repeat and transposons might have a dominating affiliation both in the shoot and fruit tissues of L. chinense. The variation in the length distribution confirmed that the sRNA transcriptomes are complicated and might be dramatically different depending on the evolutionary status of the target species. The results of this study indicate that the MIR166, MIR159, MIR167, MIR162, MIR156 and MIR5301 miRNA families have relatively higher abundance (Fig. 2). Among these families, the highly conserved miRNAs, such as MIR156, MIR159, and MIR166, tended to have greater abundance. However, the non-conserved miRNAs, such as MIR8011, MIR8021, MIR807, MIR1919, had a relatively lower abundance.

There are minimum annotating criteria regarding the biogenesis and expression in which miRNAs can be folded into a stem-loop structure, and mature miRNAs sequences should be detected by northern blotting, qRT-PCR or sequencing analysis [43]. Furthermore, the sequence detection or cloning data of miRNA* biogenesis has been required as a proof of biogenesis because miRNA*s are complementary to miRNAs sequences in the precursor [16, 44]. Among the 30 identified mature novel miRNA candidates, 16 were found with complementary miRNA*s (Table 3). Therefore, these novel mRNAs with miRNA* are likely to be recently evolved (young miRNAs) and species-specific to L. chinense. The secondary structures of the identified novel miRNAs are listed in S3 Table.

The length of small RNAs in the shoot and fruit libraries was enriched from 21–24 nt, a typical length range in plants (Fig. 1). The small RNAs of 24 nt had a higher expression in both the shoot and fruit sRNA libraries. In addition to the length distribution of small RNAs, there was also a distinct difference between the proportions of various small RNA categories (Table 2). The proportion of rRNA, tRNA, snRNA and snoRNA were higher in the shoot than in the fruit. However, a proportion of the known miRNAs was higher in the shoot, while that of novel miRNAs and that of unannotated reads were higher (82.76%) in the fruit (Table 2). Therefore, the two small RNA libraries exhibited remarkable differences in the small RNA transcriptome in both length distribution pattern and composition, indicating that fruit-related miRNAs are worth more attention in future studies.

Target gene annotation of fruit-related miRNAs in L. chinense

It is possible to identify lowly expressed miRNAs to elucidate the regulation network through high-throughput sequencing. High-throughput technology was first used to identify small RNAs in the fruit, including miRNAs in young green tomatoes, and then extended to different tomato fruit developmental stages [33, 45, 46]. Not surprisingly, a single miRNA can regulate multi-genes, and similarly multiple miRNAs can be used to regulate a single gene, indicating that the miRNA-gene regulation network might be extremely complicated. Li et al. [47] recently identified a set of nine new miRNAs in Solanaceae (tobacco, potato, and tomato), that commonly target R (resistance) genes. Five miRNAs (miR6022–miR6024, miR6026, and miR6027) are expressed during tomato fruit development. In this study, miR5303, miR5303h, miR5303g and miR5303j regulated the same gene, comp66360_c1, which is responsible for aminopeptidase activity (S6 Table).

In this study, 7 differentially expressed known miRNAs and 5 novel miRNAs were fruit-specific. A total of 227 target genes were identified for fruit-specific known miRNAs, and one target was found for only one fruit-specific novel miRNAs. Xu et al. [22] reported that miRNA transcriptome regulated lycopene biosynthesis in a red-flesh sweet orange mutant, and three groups of miRNA-target genes were involved in lycopene accumulation. A key gene of the carotenoid biosynthesis pathway is EY752486, which is a potential target of csi-miR167. EY752486 encodes geranylgeranyl pyrophosphate synthase (GGPS) and lycopene β-cyclase (LYCb), which are key enzymes for the production of prerequisite molecules and the conversion of lycopene to downstream cyclic carotenes, respectively [22]. A similar observation was reported in this study, where homologs of miR167, such as miR167a and miR167b, showed a significantly higher expression level in the fruit compared with the shoot, suggesting that miR167 also plays important roles in L. chinense fruit. Another interesting finding was that the target gene APETALA2 (AP2), designated as RAP2.2 in Arabidopsis, interacts with the phyteone synthase (PSY) gene encoding a rate-limiting enzyme for the biosynthesis of carotenoids [48]. The RAP2.2 transcription factor was also a target gene of the csi-miR172 homologs, miR172 and miR172d-3p in L. chinense, and these two miRNAs had a remarkable low expression level in the fruit. Moreover, our GO results identified genes that are involved in photosynthesis or carbon fixation in plastids or chloroplasts. A previous study of the tomato hp-2 mutant with a high level of lycopene reported that photosynthesis-related genes are consistently up-regulated throughout fruit ripening [49]. Taken together, miRNAs that regulate biological processes, including photosynthesis, transcription regulation and carotenoid biosynthesis, should be given more attention in the study of secondary metabolite accumulation in L. chinense fruits.

Further functional analyses of fruit-related miRNAs are also required to elucidate the function of these miRNAs during fruit development. miR156, during ripening regulation in tomato, targets a member of the squamosa-promoter binding protein (SPB) family called CNR [24]. In this study, miR156 was also significantly highly expressed in the fruit. Ethylene is one of the most important phytohormones and plays important roles in the plant life cycle, including seed germination, fruit ripening, and various abiotic stresses [25]. miR159 and miR394 are associated with ethylene response in rice [50]. In L. chinense, miR394 exhibited an up-regulated expression in the shoot, and miR159 was highly abundant in both the shoot and the fruit, indicating that ethylene plays regulatory roles in both the vegetative and reproductive phases. miR159 had another target, designated as 1-aminocyclopropane-1-carboxylate synthase, which also plays important roles in ethylene biosynthesis [26]. miR394-targeted F-box proteins are also important in the signal transduction pathways of various plant hormones [51, 52].

To date, fruit-specific miRNA identification has been reported on a very limited scale and includes cloning approaches in tomato [33] and grape [53] and deep sequencing approaches in peach [54], tomato [46] and orange [22]. The results of these reports strongly suggested that the expression and function of fruit-specific miRNAs throughout plant species are lacking, and fruit ripening miRNAs might functions in different organs besides fruits [22].

Target genes annotation of shoot-specific miRNAs in L. chinense

In this study, 15 miRNAs (14 known and 1 novel) were identified as shoot-specific, and a total of 179 targets were discovered for 14 known miRNAs, while the novel miRNA did not have any target. In addition, six known miRNAs were identified with a high expression level in the shoot tissues: miR160a, miR171d, miR172d-3p, miR319, miR394 and miR5301. miR160 regulates auxin response factors (ARFs), which play a central role in the auxin-regulated gene expression of primary response genes [55]. miR319 is a regulator of Lanceolate (La), which is required for compound leaf development in tomato [56]. miR394 has a regulatory effect in abiotic stress, and miR172 could regulates the expression of the AP2-like transcription factor [57]. Moreover, miR172 is a regulator of the developmental phases and specific floral organs in Arabidopsis thaliana [57, 58]. Karlova et al. (2013) reported that the MYB class of transcription factors is the target of miR5301. Additionally, miR5301 and some of its reported homologs might function in signal transduction or disease resistance [59]. The target genes of shoot-related miRNAs suggest that the miRNAs that are specific to shoot tissues might have broad regulatory functions in transcription factors, abiotic stresses, disease resistance and vegetative development in L. chinense.

Conclusion

This study reported the first set of miRNAs in L. chinense and identified 90 conserved and potentially novel miRNAs in shoot and fruit small-RNA libraries, of which 28 had a significantly higher expression level in the fruit than in the shoot of L. chinense. Further annotation revealed 30 over-represented GO categories, in which 975 genes were involved in ‘molecular functions’, ‘cellular component’ and ‘biological processes’. Additionally, the 20 most enriched pathways were identified, and 68 target genes were involved. Therefore, this study provides a unique source of miRNAs and target genes for studying biological pathways, especially for the medicinal properties of L. chinense.

GO terms of statistically enriched target genes in biological system ontology.

A: Biological process; B: Cellular component; C: Molecular function (GO terms in the square indicate the top enriched terms; red color indicates a higher level of enrichment significance than the pink color).

(TIF)

Click here for additional data file.

Plant parts of Lycium chinense.

A: Fruits bearing shoot; B: Three developmental stages of fruit and shoot tip. g, green stage; cb, color breaking stage; mr, mature red stage.

(TIF)

Click here for additional data file.

All (known and novel) differentially expressed miRNAs in the shoot and fruit tissues of L. chinense.

(DOCX)

Click here for additional data file.

GO enrichment results with description of GO accessions and corresponding genes in the L. chinense shoot and fruit.

(XLSX)

Click here for additional data file.

The secondary structures of the identified novel miRNAs in L. chinense.

(PDF)

Click here for additional data file.

The most enriched pathways that were identified for target genes.

(XLSX)

Click here for additional data file.

The designed primers and candidate list of miRNAs for qRT-PCR with differential expression through deep sequencing.

(DOCX)

Click here for additional data file.

Annotation of target genes of known and novel miRNAs in L. chinense.

(XLSX)

Click here for additional data file.

The designed primers for RLM 5’-RACE and targets annotation.

(DOCX)

Click here for additional data file.