- Research article
- Open Access
Genome-wide identification of leaf abscission associated microRNAs in sugarcane (Saccharum officinarum L.)
BMC Genomics volume 18, Article number: 754 (2017)
Sugarcane (Saccharum officinarum L.) is an economically important crop, mainly due to the production of sugar and biofuel (Azevedo RA, Carvalho RF, Cia MC, & Gratão PL, Trop Plant Biol 4:42-51, 2011). Grown mainly in tropical and subtropical countries, sugarcane is a highly polyploid plant with up to ten copies of each chromosome, which increases the difficulties of genome assembly and genetic, physiologic and biochemical analyses. The increasing demands of sugar and the increasing cost of sugarcane harvest require sugarcane varieties which can shed their leaves during the maturity time, so it is important to study the mechanism of leaf abscission in sugarcane.
To improve the understanding of miRNA roles in sugarcane leaf abscission, we reported the genome-wide characterization of miRNAs and their putative targets in sugarcane using deep sequencing for six small RNA libraries. In total, 93 conserved miRNAs and 454 novel miRNAs were identified in sugarcane using previously reported transcriptome as reference. Among them, 25 up-regulated and 13 down-regulated miRNAs were identified in leaf abscission sugarcane plants (LASP) compared to leaf packaging sugarcane plants (LPSP). Target prediction revealed several miRNA-mRNA modules including miR156-SPL, miR319-TPR2, miR396-GRF and miR408-LAC3 might be involved in the sugarcane leaf abscission. KEGG pathway enrichment analysis showed differentially expressed miRNAs may regulate pathways like “plant hormone signal transduction” and “plant-pathogen interaction”, which is consistent with previous transcriptome study. In addition, we identified 96 variant miRNAs with 135 single nucleotide polymorphisms (SNPs). The expression of sugarcane miRNAs and variant miRNAs were confirmed by qRT-PCR. We identified a possible poaceae specific miRNA called miR5384 for the first time in sugarcane.
We not only reported miR5384, a possible poaceae specific miRNA, for the first time in sugarcane but also presented some miRNA-mRNA modules including miR156-SPL, miR319-TPR2, miR396-GRF and miR408-LAC in sugarcane. These modules might be involved in the regulation of sugarcane leaf abscission during the maturity time. All of these findings may lay ground work for future application of sugarcane breeding program and benefit research studies of sugarcane miRNAs.
Abscission is the programmed developmental process that some of the organs such as leaves, flowers, or fruits are shed during the life of a plant , which can be divided into four major steps: (i) development of the abscission zone (AZ) tissue, (ii) acquisition of competence to respond to abscission-promoting signaling, (iii) activation of abscission and (iv) post abscission trans-differentiation . Global gene expression studies have shown many genes in multiple pathways participate in the abscission process [3, 4], including various transcription factors (TFs), cell wall hydrolysis enzymes, defense-related genes and genes involved in auxin/ethylene signal transduction [5,6,7]. However, the regulation mechanism of gene expression is not clear at present.
MicroRNAs (miRNAs) are endogenous small (21–24 nt) and single stranded noncoding RNAs which are important regulators of gene expression through mRNA degradation, translational repression and chromatin modification [8,9,10]. They are incorporated in the RNA-induced silencing complex (RISC) with AGO protein and guide the cleavage or translational repression of the target mRNAs by complementary or uncomplimentary base pairing . In plants, miRNAs are involved in multiple crucial biological processes including organ development and plant responses to environmental stresses [12,13,14,15,16]. Notably, several miRNAs have been reported to be involved in the abscission process. miR159 targeting MYB transcription factors, miR160/miR167 targeting auxin response factors (ARFs), miR172 targeting AP2-like ethylene-responsive transcription factors, and miR396 targeting glutamate decarboxylase have been found with different expression in tomato during pedicel abscission . In fruit senescence of Fragaria ananassa, NAC TFs, ARFs and MYB TFs have been validated as targets of miR164, miR160, miR167 and miR159, respectively, by small RNA sequencing and degradome sequencing .
Sugarcane (Saccharum officinarum L.) is an economically important crop, mainly due to the production of sugar and biofuel . Grown mainly in tropical and subtropical countries, sugarcane is a highly polyploid plant with up to ten copies of each chromosome, which increases the difficulties of genome assembly and genetic, physiologic and biochemical analyses. The increasing demands of sugar and the increasing cost of sugarcane harvest require sugarcane varieties which can shed their leaves during the maturity time, so it is important to study the mechanism of leaf abscission in sugarcane.
Recently, we have reported leaf abscission associated genes in sugarcane using transcriptome sequencing . To study miRNA expression changes between leaf abscission sugarcane plants (LASP) and leaf packaging sugarcane plants (LPSP), we constructed six small RNA libraries and sequenced them using an Illumina HiSeq 2500 system. In total, we characterized 93 conserved miRNAs and 454 sugarcane novel miRNAs using the transcriptome sequences as reference and identified 38 differentially expressed miRNAs in LASP compared to LPSP. Target prediction showed several miRNA-mRNA modules might be involved in sugarcane leaf abscission, such as miR156-SPL, miR319-TPR2, miR396-GRF and miR408-LAC3. KEGG pathway analysis for the target genes showed similar results as our transcriptome study and indicated “plant-pathogen interaction” and “plant hormone signal transduction” might be related with sugarcane leaf abscission during the maturity time. We obtained highly conserved sugarcane pre-miRNAs and mature miRNAs, which are valuable resources for future sugarcane miRNA studies. Moreover, our findings will provide better understanding of the complex mechanism of leaf abscission, miRNAs and their targets involved in leaf abscission.
Small RNA identification in sugarcane
Previously, we reported leaf abscission associated genes in sugarcane using transcriptome sequencing . To study the post-transcriptional regulation of sugarcane leaf abscission, six small RNA (sRNA) libraries constructed for three leaf abscission sugarcane plants (Q1, T1, T2) and three leaf packaging sugarcane plants (Q2, B1, B2) were sequenced by using an Illumina HiSeq 2500 system. Initially, a total of 71,579,415 raw reads were generated (Table 1). After low quality reads and sequencing adapters were removed, we obtained 70,082,453 clean reads longer than 18 nt for all samples with an average of 11.68 M clean reads. Length distribution of clean reads (Fig. 1a) showed the most abundant classes of sRNA were 21 and 24 nt, which is consistent with many plant miRNA deep sequencing studies [21,22,23]. To annotate sRNAs in each library, we mapped the clean reads to Rfam database and found 5.10–10.93% of the clean reads were derived from rRNA, tRNA, snRNA and snoRNA (Table 1). Next, all clean reads were aligned to recently published sugarcane transcriptome  for global miRNA characterization. Mapped reads were used to predict sugarcane miRNAs (pre-miRNAs and mature miRNAs) by MIREAP  and a total of 93 conserved miRNAs (including mature and passenger miRNAs) from 25 families and 454 sugarcane novel miRNAs were identified (Additional file 1: Supplementary Dataset). These miRNAs were used as references for sRNA mapping, miRNA expression profiling and SNP scanning. Notably, a possible poaceae specific miRNA miR5384, which has been identified only in Sorghum bicolor and Triticum aestivum in miRBase, is reported for the first time in sugarcane.
miRNA expression profile
We next profiled sugarcane miRNA expression in each sample by mapping the clean reads to sugarcane miRNA precursors using BLAST software . sRNA reads matched to sugarcane miRNAs without mismatches were counted. To compare miRNA expression in different libraries, in this study the number of clean reads was used as background for normalization and TPM (transcripts per million reads) was used to present the expression levels of miRNAs. A total of 439 miRNAs were detected ((Additional file 2: Table S1) and 340, 357, 349, 348, 368 and 313 miRNAs were detected more than 1 TPM in Q1, Q2, T1, T2, B1and B2, respectively (Table 1). Distribution of normalized miRNA expression (Fig. 1b) revealed approximate 71.47% ~ 82.75% of the total detected miRNAs were expressed no more than 5 TPM (excluding miRNA with 0 TPM), and 26 (7.65%), 18 (5.04%), 19 (5.44%), 21 (6.03%), 19 (5.16%) and 14 (4.47%) miRNAs were identified >100 TPM in Q1, Q2, T1, T2, B1 and B2, respectively. For subsequent analyses, miRNAs whose normalized expression values were lower than 5 TPM in all samples were removed. Venn diagram (Fig. 1c) revealed 32 miRNAs were common to all six samples, and 11, 3, 3, 4, 1 and 3 miRNAs were detected more than 5 TPM only in Q1, Q2, T1, T2, B1 and B2, respectively.
Analysis of differentially expressed miRNAs
In general, miRNAs differentially expressed in LASP and LPSP are associated with leaf abscission process in sugarcane. We next identified differentially expressed miRNAs between LASP and LPSP using edgeR . In this study, we used several statistical values as cut-offs: normalized expression > 5 TPM, log2 fold change (Log2FC) > 1 (up-regulated) or Log2FC <− 1(down-regulated), p-value < 0.05 and FDR < 0.05. With this criterial we obtained 25 up-regulated and 13 down-regulated miRNAs in LASP compared to LPSP (Fig. 2a, Table 2). It is clear that some miRNAs were differentially expressed significantly in LASP and LPSP because their statistical values were far away from the log2 fold change and FDR cut-offs, including sugarcane known miRNAs (miR167b-5p, miR167b-3p, miR167c-5p, miR408-3p, miR319-3p and miR159-5p) and novel miRNAs (miRN167-3p, miRN245-5p and miRN303-5p). It is interesting that both mature (also known as functional strand) and passenger (also known as the miRNA* strand) miRNAs from MIR167 family were up-regulated in LASP. Passenger miRNAs are usually thought to be degraded quickly during the miRNA biogenesis but they still can be active in silencing . Like Arabidopsis thaliana and Oryza sativa, passenger miRNAs from MIR167 family are derived from the 3′-arms of the pre-miRNAs in sugarcane. Sugarcane passenger miRNA of MIR167 (miR167-3p) has been reported with high abundance in sugarcane root and predicted to be functional in plant development . This is the first time to report leaf abscission associated miRNAs in sugarcane. A total of four pairs of mature and passenger miRNAs were differentially expressed between LASP and LPSP and they were derived from miRNA families like MIR156, MIR167 and MIR393 (Table 2).
Target prediction for differentially expressed miRNAs
The identification of miRNA targets and their regulation is a crucial step to understand the biological functions of miRNAs. Using the assembled sugarcane transcriptome we computationally predicted the target genes for up- and down-regulated miRNAs using the method recommended by Allen  and Schwab . As shown in Fig. 2b, we predicted 510 and 488 target genes for up-regulated and down-regulated miRNAs, respectively (Additional file 1: Supplementary Dataset), and there were 27 target genes shared by up- and down-regulated miRNAs. It is revealed that all differentially expressed miRNAs between LASP and LPSP were predicted to have target genes in sugarcane leaves (Fig. 2c). Among up- and down-regulated miRNAs, miR5384-3p and miRN141-3p were predicted regulate most target genes, respectively (Fig. 2c).
Annotation of miRNA target genes (Additional file 1: Supplementary Dataset) revealed several conserved miRNA-mRNA modules identified, including miR156-SPL (squamosa promoter-binding-like protein), miR319-TPR (topless-related protein), miR396-GRF (growth-regulating factor) and miR408-LAC3 (laccase-3). In this study, miR156 was predicted to target 15 SPL genes (e.g. SPL3, SPL11, SPL12, SPL14, SPL17 and SPL19), miR396a/b targeting 8 GRF genes (e.g. GRF3, GRF4, GRF9 and GRF12) and miR408 targeting one LAC3 gene (laccase-3). As shown in Fig. 2d, using both sRNA and transcriptome  data we showed the expression changes of normalized expression values of miRNAs and mRNAs from these modules in LASP and LPSP. It is interesting that all these three miRNAs (miR156, miR396 and miR408) were up-regulated while their target genes (SPL, GRF and LAC3) were down-regulated significantly (p-value < 0.05) in LASP compared to LPSP, which increased the reliability of the regulation of miRNA-mRNA in sugarcane leaf abscission. In addition, miR319-3p was down-regulated in LASP compared to LPSP, but its putative target gene TPR2 (topless-related protein 2) was up-regulated significantly (Fig. 2d). We also found opposite expression patterns of novel miRNAs and their target genes (Additional file 1: Supplementary Dataset), such as miRN245-RGA2 (disease resistance protein RGA2) and miRN375-GRF (growth-regulating factor).
Functional analysis of differentially expressed miRNAs
Function analysis for miRNA target genes showed differentially expressed miRNAs were involved in multiple pathways (Fig. 2e and Table 3). Gene Ontology analysis (Fig. 2d) revealed differentially expressed miRNAs between LASP and LPSP may function in “cellular process” (GO:0009987), “metabolic process” (GO:0008152), “response to stimulus” (GO:0050896), “cell” (GO:0044464), “membrane” (GO:0016020), “binding” (GO:0005488) and “catalytic activity” (GO:0003824) while KEGG pathway enrichment analysis (Table 3) revealed 5 significant pathways (p-value < 0.05 and q-value < 0.05) were regulated by differentially expressed miRNAs between LASP and LPSP. It is notable that 22 genes that were targets of 5 differentially expressed miRNAs (miR171c-3p, miR393-5p, miR5384-3p, miR167a-3p and miRN141-3p) were involved in the pathway of “plant hormone signal transduction” (ko04075) and both miR393-5p and miR393-3p were involved in the same pathway - “Glycosylphosphatidylinositol(GPI)-anchor biosynthesis” (ko00563). Two cell growth and death pathways were found to be regulated by two miRNAs (miRN245-5p and miR5384-3p). Detailed information of KEGG pathway enrichment analysis can be accessed in Additional file 2: Table S2. We found other not significant pathways regulated by differentially expressed miRNAs maybe important for leaf abscission in sugarcane as well. For example, “flavonoid biosynthesis” (ko00941) regulated by miR5384-3p, “apoptosis” (ko04210) regulated by miR5384-3p and miR166e-5p, and “plant-pathogen interaction” (ko04626) regulated by 11 miRNAs (miR166e-5p, miR167b-5p, miR393-3p, miR159-5p, miRN014-3p, miRN141-3p, miR396b-5p, miR396a-5p, miR167c-5p, miR167a-5p and miRN245-5p). In our transcriptome study , “plant-pathogen interaction” was the most significant pathway involved by differentially expressed genes. In current study, 20 genes from this pathway were predicted to be regulated by 11 miRNAs mentioned before, which indicated that “plant-pathogen interaction” may be important for sugarcane leaf abscission. Due to the limit annotations of the assembled transcriptome, functions of many target genes of differentially expressed miRNAs are currently unknown.
The biogenesis of miRNA, guide miRNA selection and miRNA-mRNA interaction indicate that single nucleotide polymorphisms (SNPs) in miRNA genes or mature miRNA should affect the biogenesis and function of miRNAs . It is implicated that sugarcane leaf abscission may be related to them, so we next scanned SNPs in sugarcane miRNAs in all samples. After variant miRNAs with low frequency (< 10 reads) and small rate (< 1%) were removed, we obtained 135 SNPs happened in 96 sugarcane miRNAs (Additional file 2: Table S3). Among these variant miRNAs two isoforms (sof-miR5564b-3p and sof-miRN396-3p) attracted our attention because of their high expression in the samples (Fig. 3a). An “A to T” and a “C to T” happened to sof-miR5564b-3p and sof-miRN396-3p, respectively. It is interesting that variant miRNAs were differentially expressed between LASP and LPSP. Variant sof-miR5564b-3p was expressed much higher in LPSP (B1 and B2) compared to LASP (T1 and T2), but there was no evidence of it in Q1 and Q2. In this study we identified two main variant isoforms of sof-miRN396-3p – “acgagaTgaatcttttgagcct” and “cgagaTgaatcttttgagcct”. They were found with high sequence similarity with normal sof-miRN396-3p (cgagaCgaatcttttgagcct). In addition, variant sof-miRN396-3p was expressed higher than its normal form in most samples. The numbers of reads mapping to normal sof-miRN396-3p were 1305, 1197, 1318, 1565, 1908 and 1092 reads in B1, B2, Q1, Q2, T1 and T2, respectively, but variant sof-miRN396-3p was aligned with 2249, 2364, 5713, 958, 1281 and 1567 reads in B1, B2, Q1, Q2, T1 and T2, respectively (Fig. 3a). In plants most miRNAs cleave their downstream targets depending on the highly complementary recognition sites . SNPs in miRNA regions, especially in their “seed” regions (2nd–7th nucleotides of the miRNA), has huge influence on miRNA biogenesis and function, such as miRNA guide strand selection and downstream target gene binding [27, 33]. We predicted the target genes using the transcriptome reference for normal and variant sof-miRN396-3p (Fig. 3b). Although normal and variant sof-miRN396-3p shared 76 target genes such as YLS9 (YELLOW-LEAF-SPECIFIC GENE 9), HOX27 (Homeobox-leucine zipper protein), DDM1 (ATP-dependent DNA helicase) and CDKA2 (Cyclin-dependent kinase A-2), variant sof-miRN396-3p losses the regulation of 11 genes (e.g. LRKS5 (L-type lectin-domain containing receptor kinase S.5), WEX (Werner Syndrome-like exonuclease), Y2242 (LRR receptor-like serine/threonine-protein kinase)) and gains11 new genes (e.g. POL2 (Retrovirus-related Pol polyprotein from transposon 297), E139 (Glucan endo-1,3-beta-glucosidase 9), Y5162 (Uncharacterized protein At5g41620)) due to the SNP (C - > T) in its “seed” region.
To validate the expression of normal miRNAs and variant miRNAs in LASP and LPSP, we performed stem-loop qRT-PCR experiment. Four miRNAs (sof-miR159-3p, sof-miR166a-3p, sof-miR396a-5p and sof-miR5564b-3p) and one variant miRNA (sof-miR5564bV-3p) were selected as candidates and actin was used as an internal control. Reverse transcriptase, forward and reverse primers were designed and synthesized in BGI-Shenzhen (Additional file 2: Table S4). After ΔCt value of each miRNA candidate was calculated in each sample, the expression of miRNAs was normalized using B1 as the control (Fig. 3b). A total of 13 out of 20 (65.5%) were identified with same expression patterns in all samples between small RNA-Seq and qRT-PCR, which support that miRNA expression profile by small RNA sequencing is reliable. To validate the variant isoform of sof-miR5564, we used two different forward primers with one nucleic acid difference and the results showed it was detectable in most of the samples. The overall correlation of normal and variant miRNAs between qRT-PCR and deep sequencing, calculated as 0.274, indicated that miRNAs identified by deep sequencing were accuracy and effective, and can be used for miRNA expression profiles analysis of sugarcane leaf abscission during maturity time.
In this study, we characterized 93 conserved miRNAs from 25 families and 454 novel miRNAs in sugarcane using deep sequencing technology. Unlike previous sugarcane miRNA studies [34,35,36,37,38,39], this is the first time to identify sugarcane miRNAs using the sugarcane transcriptome reference. Using sugarcane transcriptome reference, we not only characterized pre-miRNAs and mature miRNAs in sugarcane but also analyze miRNA-mRNA interactions. Pre-miRNAs identified in this study are more reliable because they have high similarities with other species. Polygenetic trees (Fig. 4 a ~ f) of six miRNA families (MIR156, MIR160, MIR166, MIR167, MIR171 and MIR396) performed by MEGA 7  revealed pre-miRNAs identified in this study were conserved across species, such as Sorghum bicolor (sbi), Zea mays (zma), Oryza sativa (osa) and other monocotyledons. In addition, multiple sequence alignment (Additional file 2: Figure S1A-F) performed by ClustalX (v 2.0)  showed mature miRNAs, step-loop regions and passenger miRNAs were conserved across species, unlike previously reported sugarcane miRNAs (marked with “D”). Considering this we submitted these sugarcane miRNAs to miRBase for updates.
Leaf abscission is a complicated process regulated by developmental, hormonal and environmental cues . In other plants multiple types of genes have been reported to relate with leaf abscission, such as auxin response proteins, ethylene-synthesis enzymes, cell wall-degrading enzymes, pathogenesis related proteins [43, 44]. In this study four miRNAs (miR156, miR319, miR396 and miR408) and their target genes (SPL, TPR2, GRF and LAC3) were identified with different expression trends (up-regulation or down-regulation) in LASP and LPSP (Fig. 2d). SPLs are trans-acting factors that bind specifically to the consensus nucleotide sequence (5′-TNCGTACAA-3′) of AP1 promoter and can promote both vegetative phase change (e.g. leaf epidermal differentiation) and flowering time [45,46,47]. There are 16, 18, 13, and 31 SPL genes in Arabidopsis thaliana, Oryza sativa, Physcomitrella patens, and Zea mays, respectively, most of which can be regulated by miR15648 . In Arabidopsis, overexpression of SPL3 that is a consequence of miR156 decrease can accelerate the production of abaxial trichome and short petioles of adult leaves . Another study also has confirmed the regulation of miR156 in SPLs and has revealed the decrease of miR156 and overexpression of SPLs (SPL3, SPL9 and SPL10) produce less lateral roots in Arabidopsis . In rice SPLs and miR156 have been reported to be involved in multiple developmental processes, especially the flower development . In our transcriptome study we have identified several SPL genes down-regulated in LASP compared to LPSP , the up-regulation of miR156 in current study strongly support that the module miR156-SPL may function in sugarcane leaf abscission.
In this study, we also found the up-regulation of miR396 and miR408 in LASP compare to LPSP. miR396 is a highly conserved miRNA (Fig. 4f) and always functions, through repressing GRFs (growth-regulating factors), in the regulation of developmental processes of leaf, flower and other organs in plants [51,52,53,54]. Kim and Lee have conducted an experiment to show that GRF4 mutation results in small leaves of Arabidopsis . Not only GRF4 but also other GRF genes have been confirmed to be regulated by miR396 and overexpression of miR396 results reduction of GRFs and narrow-leaf phenotypes in Arabidopsis . In other plants such as tobacco and citrus miR396 has also been reported to target GRFs and regulate the development of leaf, flower and fruit [54, 55]. In addition, it has been experimented that overexpression of miR396 can reduce the ability of salt and alkali stress tolerance in rice and Arabidopsis . Unlike miR396, the overexpression of miR408 can increase drought tolerance by targeting different genes in plants, such as rice , chickpea  and Medicago truncatula . In Arabidopsis miR408 is increased while its target LAC3 is decreased during senescence . Similarly, current study showed correlated expression patterns of miR396-GRF and miR408-LAC3 in LASP and LPSP.
By regulating MYB and TCP transcription factors, miR319 is involved in various developmental processes such as leaf development and senescence, organ curvature and hormone biosynthesis and signaling [61,62,63]. In this study we found miR319 can target MYB transcription factor but the expression of MYB was not significantly changed in LASP and LPSP (Additional file 1: Supplementary Dataset). However, TPR2 (topless-related protein 2) was predicted to be regulated by miR319 and its expression was correlated with miR319 (Fig. 2d). TPR2 is a transcriptional corepressor involved in the regulation of branch formation, strigolactones signaling and auxin signaling [64, 65]. At the early leaf development stage of Arabidopsis TIE1 (TCP Interactor containing EAR motif protein 1) can regulate leaf size by inhibiting TCP activities through recruiting the TOPLESS (TPL)/TOPLESS-RELATED (TPR) corepressors . The down-regulation of miR319 and up-regulation of its target TPR2 in LASP compared to LPSP indicate the module of miR319-TPR2 might be involved in sugarcane leaf abscission. Due to the limit annotations for miRNA target genes, it is hard to understand very well about the functions of all differentially expressed miRNAs between LASP and LPSP. Further experiments should be conducted to explore the functions of differentially expressed miRNAs and their targets in sugarcane leaf abscission.
In conclusion, we successfully identified 547 miRNAs in sugarcane, of which 25 were up-regulated and 13 were down-regulated in LASP compared to LPSP. We not only reported miR5384, a possible poaceae specific miRNA, for the first time in sugarcane but also presented some miRNA-mRNA modules including miR156-SPL, miR319-TPR2, miR396-GRF and miR408-LAC in sugarcane. These modules might be involved in the regulation of sugarcane leaf abscission during the maturity time. All of these findings may lay ground work for future application of sugarcane breeding program and benefit research studies of sugarcane miRNAs.
Plant material and growth conditions
New born leafs from six sugarcane cultivars (Q1, Q2, T1, T2, B1 and B2) with different phenotypes on leaf abscission during the maturity time, as described . Briefly, Q2 (ROC-26) is a precocious and productive sugarcane variety from Taiwan, but the sugar cane is wrapped tightly at the physiological maturity . Q1 (GT96–167) is a late-maturing and high-yield sugarcane variety which is bred by Guangxi Sugarcane Research Institute [68,69,70]. In contrast with Q2, Q1 can shed their leaves easily during the maturity time. T1, T2, B1 and B2 are four generation varieties of Q1 and Q2. T1 and T2 can shed their leaves as well during the maturity time, rather than B1 or B2. All the sugarcane varieties were proved to have stable agronomic characteristic on leaf shedding by 5 years of filed observation. All six sugarcane plants were planted in January of 2014 in the experimental field of Sugarcane Research Institute in Nanning, Guangxi Province of China. In December of 2014, after the sugarcane plants were validated in maturity time according to the sugar assay test, new born leaf tissues approximately 5 cm above the growing point were taken and stored in the liquid nitrogen before RNA extraction.
Total RNA extraction
Leaf tissues (100 mg) of Q1, Q2, T1, T2, B1 and B2 were used to isolate total RNA by using TRIzol® reagent (Invitrogen) according to the manufacture’s protocol . The quality of total RNA was controlled by Agilent 2100 Bioanalyzer.
Small RNA library construction and deep sequencing
Six small RNA libraries (Q1, Q2, T1, T2, B1 and B2) were constructed and sequenced with Illumina TruSeq deep sequencing technology (Sample Preparation Guide, Par #15004197 Rev.A, Illumina, San Diego, CA). Briefly, small RNAs (18–30 nt) were size-selected by gel fraction and extracted by centrifugation. After ligation of 5′ and 3′ adaptors, small RNAs were reverse transcribed into cDNA, then amplified using the sequencing primers for 14 cycles and the fragments (~ 150 bps) were isolated from a 6% TBE PAGE-gel. After the cDNA was purified, it was used for cluster generation and sequenced using an Illumina HiSeq 2000 platform. Image files generated by the sequencer were processed to nucleotide sequences (raw FASTQ files) using a base-calling pipeline (Illumina). FASTQ files for all six libraries have been submitted to Sequence Read Archive (SRA) of NCBI under the accession number SRA347706 (Q1: SRR3184695, Q2: SRR3184696, T2: SRR3184697, T1: SRR3184698, B1: SRR3184699, and B2: SRR3184700).
Reference transcriptome de novo assembly and annotation
Previously, we have demonstrated the transcriptome analysis for these six sugarcane varieties . The raw sequencing data is available with the accession number SRA291189 from SRA platform. The assembled transcriptome was used as the reference for small RNA identifications in current study . Transcriptome sequences were assembled by Trinity software . Then, the assembled transcripts were annotated by mapping them to NCBI non-redundant (NR), UniProtKB/Swiss-Prot, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. The mapping hits were filtered using a cut-off of e-value (1 × 10−5), as previously described . The assembled transcriptome sequences can be accessed in NCBI Transcriptome Shotgun Assembly (TSA) Database using the accession number GELY00000000.
Identification of sugarcane pre-miRNAs and mature miRNAs
Raw small RNA reads were cleaned by removing the adaptors, short (< 18 nt) reads and low quality reads, resulting clean reads were used for sequential analysis. To identify sugarcane miRNA precursors, we mapped all clean reads to the assembled transcriptome reference using SOAP2 software  under maximum two mismatches. Then, MIREAP  was used for global miRNA search in the assembled transcriptome of sugarcane  with parameters as follows: minimal length, 18; maximal length, 25; maximal copy number of miRNAs on reference, 20; maximal free energy allowed for a miRNA precursor, − 18 kcal/mol; maximal space between mature miRNA and passenger miRNA, 300; minima space between mature miRNA and passenger miRNA, 16; maximal bugle of mature miRNA and passenger miRNA, 4; flank sequence length of miRNA precursor, 20. Predicted pre-miRNAs and mature miRNAs were filtered with the expression and structure features by MIREAP. Two pre-miRNAs overlapped more than 80% nucleotides were considered as one. The structures of pre-miRNAs were predicted using RNAfold from ViennaRNA Package 2.0 . Next, predicted sugarcane mature miRNAs were compared to miRBase  (v21, http://www.mirbase.org) for conserved miRNA identification. Conserved miRNAs were selected according to their similarities (> 90%) and expression (> 10 reads). Apart of these conserved miRNAs, the left miRNAs were considered as sugarcane novel miRNAs. All the sugarcane miRNA precursor and mature sequences have been submitted to miRBase.
To profile miRNA expression in each sample, BLAST software  was used to align the clean reads to sugarcane pre-miRNAs in each sample under perfect match. For a particular miRNA, its expression was calculated by counting the reads which had an overlap more than 18 nt with the miRNA. Then, the miRNA expression was normalized to the total clean reads using TPM (transcripts per million reads) method:
Different expression of miRNAs
Differentially expressed miRNAs between LASP and LPSP were identified using edgeR . We used a strict criterial for selecting miRNAs with differential expression: normalized expression > 5 TPM, Log2FC > 1 or Log2FC < − 1, p-value < 0.05 and FDR < 0.05.
Target prediction and function annotation of miRNAs
The assembled sugarcane transcriptome was aligned against by the mature miRNA sequences for target prediction, according to the suggestions by Allen  and Schwab . Here, the criteria used for miRNA target prediction were as follows: i) no more than three mismatches between the miRNA and the target (G-U bases count as 0.5 mismatches); ii) no more than two adjacent mismatches in the miRNA/target duplex; iii) no adjacent mismatches in in positions 2–12 of the miRNA/target duplex (5′ of the miRNA); iv) no mismatches in positions 10–11 of miRNA/target duplex; v) no more than 2.5 mismatches in positions 1–12 of the of the miRNA/target duplex (5′ of the miRNA); and vi) the Minimum Free Energy (MFE) of the miRNA/target duplex should be ≥75% of the MFE of the miRNA bound to its perfect complement. Target genes were then used for GO and KEGG pathway enrichment analysis. We used p-value (Fisher’s exact test) and q-value  to show the significance of enrichment and control the false discovery rate. Significant GO items and KEGG pathways should satisfy the critical of p-value < 0.05 and q-value < 0.05. Detected KEGG pathways related to animal or human GO items were filtered.
To validate the expression patterns of miRNAs in sugarcane leaves, we selected four miRNAs (sof-miR159-3p, sof-miR166a-3p, sof-miR396a-5p and sof-5564b-3p) and one variant miRNA (sof-miR5564bV-3p) to perform stem-loop reverse transcription and quantitative real-time PCR (RT-qPCR), following the protocols . Initially, total RNA (2 μg) isolated from six sugarcane leaf samples was added into dNTP (2 μl, 2.5 mM, Geneland), 10 × RT Buffer (2 μl, QIAGEN), RT primer mix (2 μl, 1 μM), RNase Inhibitor (1 μl, Promega) and Quantiscript® Reverse Transcriptase (1 μl, 10 U/μl, QIAGEN) and diluted with RNase-free water to reverse transcription mix (20 μl). The reverse transcription mix was then incubated at 16 °C for 30 min, 42 °C for 40 min and at 85 °C for 5 min to finish cDNA synthesis. Next, cDNA template (1 μl) was mixed with 2× SYBR Green PCR mix (8 μl, QIAGEN), forward primer (0.2 μl, 50 pM/μl), reverse primer (0.2 μl, 50 pM/μl) and RNase-free H2O (6.6 μl) to make the PCR reaction mix. RT primers for cDNA synthesize and forward and reverse primers for RT-qPCR were designed and synthesized in BGI-Shenzhen (Additional file 2: Table S4). We used actin as an internal control and performed three replicates for each miRNA in every sample. The PCR reaction was performed and analyzed by the ABI ViiA 7 Real Time PCR System. The RT-qPCR conditions were as follows: preheating for 2 min at 95 °C; and 50 cycles of 94 °C for 10 s, 56 °C for 10 min and 72 °C for 40 s. All PCR products were denatured at 95 °C and cooled to 65 °C, and the fluorescence signals were accumulated consistently from 65 °C to 95 °C as the temperature increased at 0.2 °C per second. The expression levels of miRNAs were evaluated by the comparative threshold cycle (Ct) values. Relatively normalized expression (RNE, −ΔΔCT method) was used to show the expression change of a transcript in two samples . CT values greater than 35 were set to 35.
SNPs in miRNA
To characterize SNPs occurred in mature miRNAs of sugarcane, clean reads were aligned to sugarcane miRNA precursors by using BLAST software  with maximum one mismatch. Then, reads with perfect match to the pre-miRNAs and outside mature miRNA regions were removed. The expression of variant miRNAs was counted like normal mature miRNAs. Variant miRNAs with low expression (< 10 reads) or small ratio (< 1%) of total reads (variant miRNA reads and normal miRNA reads) were filtered. miRNAs that has two or more SNPs were not considered (Additional file 3).
Auxin response factors
Kyoto Encyclopedia of Genes and Genomes
Leaf abscission sugarcane plants
Leaf packaging sugarcane plants
Open reading frame
Quantitative real-time polymerase chain reaction
Single nucleotide polymorphism
Sequence read archive
Transcripts Per Kilobase Million
Sexton R, Roberts JA. Cell biology of abscission. Annu Rev Plant Physiol. 1982;33:133–62.
Patterson SE. Cutting loose. Abscission and dehiscence in Arabidopsis. Plant Physiol. 2001;126:494–500.
Gao Z, et al. High-throughput sequencing of small RNAs and analysis of differentially expressed microRNAs associated with pistil development in Japanese apricot. BMC Genomics. 2012;13:371.
Song QX, et al. Identification of miRNAs and their target genes in developing soybean seeds by deep sequencing. BMC Plant Biol. 2011;11:5.
Meir S, et al. Microarray analysis of the abscission-related transcriptome in the tomato flower abscission zone in response to auxin depletion. Plant Physiol. 2010;154:1929–56.
Meir S, et al. Identification of defense-related genes newly-associated with tomato flower abscission. Plant Signal Behav. 2011;6:590–3.
Cho HT, Cosgrove DJ. Altered expression of expansin modulates leaf growth and pedicel abscission in Arabidopsis Thaliana. Proc Natl Acad Sci U S A. 2000;97:9783–8.
Reinhart BJ, Weinstein EG, Rhoades MW, Bartel B, Bartel DP. MicroRNAs in plants. Genes Dev. 2002;16:1616–26.
Palatnik JF, et al. Control of leaf morphogenesis by microRNAs. Nature. 2003;425:257–63.
Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116:281–97.
Kurihara Y, Watanabe Y. Arabidopsis Micro-RNA biogenesis through dicer-like 1 protein functions. Proc Natl Acad Sci U S A. 2004;101:12753–8.
Sunkar R, Kapoor A, Zhu JK. Posttranscriptional induction of two cu/Zn superoxide dismutase genes in Arabidopsis is mediated by downregulation of miR398 and important for oxidative stress tolerance. Plant Cell. 2006;18:2051–65.
Li W, et al. Transcriptional regulation of Arabidopsis MIR168a and argonaute1 homeostasis in abscisic acid and abiotic stress responses. Plant Physiol. 2012;158:1279–92.
Song JB, et al. miR394 and LCR are involved in Arabidopsis salt and drought stress responses in an abscisic acid-dependent manner. BMC Plant Biol. 2013;13:210.
Jones-Rhoades MW, Bartel DP. Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol Cell. 2004;14:787–99.
Mallory AC, Vaucheret H. Functions of microRNAs and related small RNAs in plants. Nat Genet. 2006;38(Suppl):S31–6.
Xu T, et al. Small RNA and degradome sequencing reveals microRNAs and their targets involved in tomato pedicel abscission. Planta. 2015;242:963–84.
Xu X, et al. High-throughput sequencing and degradome analysis identify miRNAs and their targets involved in fruit senescence of Fragaria Ananassa. PLoS One. 2013;8:e70959.
Azevedo RA, Carvalho RF, Cia MC, Gratão PL. Sugarcane under pressure: an overview of biochemical and physiological studies of abiotic stress. Trop Plant Biol. 2011;4:42–51.
Li M, et al. De novo analysis of transcriptome reveals genes associated with leaf abscission in sugarcane (Saccharum Officinarum L.). BMC Genomics. 2016;17:195.
Zhang L, et al. Identification and temporal expression analysis of conserved and novel microRNAs in sorghum. Genomics. 2011;98:460–8.
Khatabi B, et al. High-resolution identification and abundance profiling of cassava (Manihot Esculenta Crantz) microRNAs. BMC Genomics. 2016;17:85.
Zhou J, et al. Identification of novel miRNAs and miRNA expression profiling in wheat hybrid necrosis. PLoS One. 2015;10:e0117507.
Qibin L, Jiang W. MIREAP: microRNA discovery by deep sequencing; 2008.
Altschul SF, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25:3389–402.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.
Ha M, Kim VN. Regulation of microRNA biogenesis. Nat Rev Mol Cell Biol. 2014;15:509–24.
Moyle RL, Sternes PR, Birch RG. Incorporating target sequences of developmentally regulated small RNAs into Transgenes to enhance tissue specificity of expression in plants. Plant Mol Biol Report. 2015;33:505–11.
Allen E, Xie Z, Gustafson AM, Carrington JC. microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell. 2005;121:207–21.
Schwab R, et al. Specific effects of microRNAs on the plant transcriptome. Dev Cell. 2005;8:517–27.
Sun G, et al. SNPs in human miRNA genes affect biogenesis and function. RNA. 2009;15:1640–51.
Voinnet O. Origin, biogenesis, and activity of plant microRNAs. Cell. 2009;136:669–87.
Gong J, et al. Genome-wide identification of SNPs in microRNA genes and the SNP effects on microRNA target binding and biogenesis. Hum Mutat. 2012;33:254–63.
Thiebaut F, et al. Differential sRNA regulation in leaves and roots of sugarcane under water depletion. PLoS One. 2014;9:e93822.
Thiebaut F, et al. Computational identification and analysis of novel sugarcane microRNAs. BMC Genomics. 2012;13:290.
Ferreira TH, et al. microRNAs associated with drought response in the bioenergy crop sugarcane (Saccharum spp.). PLoS One. 2012;7:e46703.
Carnavale Bottino M, et al. High-throughput sequencing of small RNA transcriptome reveals salt stress regulated microRNAs in sugarcane. PLoS One. 2013;8:e59423.
Viswanathan C, Anburaj J, Prabu G. Identification and validation of sugarcane streak mosaic virus-encoded microRNAs and their targets in sugarcane. Plant Cell Rep. 2014;33:265–76.
Ortiz-Morea FA, et al. Global analysis of the sugarcane microtranscriptome reveals a unique composition of small RNAs associated with axillary bud outgrowth. J Exp Bot. 2013;64:2307–20.
Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33(7)1870–4.
Larkin MA, et al. Clustal W and Clustal X version 2.0. Bioinformatics. 2007;23:2947–8.
Tadeo FR, et al. Molecular physiology of development and quality of citrus. Adv Bot Res. 2008;47:147–223.
Meir S, Hunter DA, Chen JC, Halaly V, Reid MS. Molecular changes occurring during acquisition of abscission competence following auxin depletion in Mirabilis jalapa. Plant Physiol. 2006;141:1604–16.
Eyal Y, Meller Y, Lev-Yadun S, Fluhr R. A basic-type PR-1 promoter directs ethylene responsiveness, vascular and abscission zone-specific expression. Plant J. 1993;4:225–34.
Cardon GH, Hohmann S, Nettesheim K, Saedler H, Huijser P. Functional analysis of the Arabidopsis Thaliana SBP-box gene SPL3: a novel gene involved in the floral transition. Plant J. 1997;12:367–77.
Wu G, Poethig RS. Temporal regulation of shoot development in Arabidopsis Thaliana by miR156 and its target SPL3. Development. 2006;133:3539–47.
Birkenbihl RP, Jach G, Saedler H, Huijser P. Functional dissection of the plant-specific SBP-domain: overlap of the DNA-binding and nuclear localization domains. J Mol Biol. 2005;352:585–96.
Preston JC, Hileman LC. Functional evolution in the plant SQUAMOSA-PROMOTER BINDING PROTEIN-LIKE (SPL) gene family. Front Plant Sci. 2013;4:80.
Yu N, Niu QW, Ng KH, Chua NH. The role of miR156/SPLs modules in Arabidopsis lateral root development. Plant J. 2015;83:673–85.
Xie K, Wu C, Xiong L. Genomic organization, differential expression, and interaction of SQUAMOSA promoter-binding-like transcription factors and microRNA156 in rice. Plant Physiol. 2006;142:280–93.
Debernardi JM, et al. Post-transcriptional control of GRF transcription factors by microRNA miR396 and GIF co-activator affects leaf size and longevity. Plant J. 2014;79:413–26.
Kim JH, Lee BH. GROWTH-REGULATING FACTOR4 of Arabidopsis Thaliana is required for development of leaves, cotyledons, and shoot apical meristem. Journal of Plant Biology. 2006;49:463–8.
Liu D, Song Y, Chen Z, Yu D. Ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth in Arabidopsis. Physiol Plant. 2009;136:223–36.
Yang F, Liang G, Liu D, Yu D. Arabidopsis MiR396 mediates the development of leaves and flowers in transgenic tobacco. Journal of Plant Biology. 2009;52:475–81.
Liu X, Guo LX, Jin LF, Liu YZ, Liu T, Fan YH, Peng SA. Identification and transcript profiles of citrus growth-regulating factor genes involved in the regulation of leaf and fruit development. Mol Biol Rep. 2016;43(10):1059–67.
Gao P, et al. Over-expression of osa-MIR396c decreases salt and alkali stress tolerance. Planta. 2010;231:991–1001.
Mutum RD, et al. Evolution of variety-specific regulatory schema for expression of osa-miR408 in indica rice varieties under drought stress. FEBS J. 2013;280:1717–30.
Hajyzadeh M, Turktas M, Khawar KM, Unver T. miR408 overexpression causes increased drought tolerance in chickpea. Gene. 2015;555:186–93.
Trindade I, Capitao C, Dalmay T, Fevereiro MP, Santos DM. miR398 and miR408 are up-regulated in response to water deficit in Medicago Truncatula. Planta. 2010;231:705–16.
Thatcher SR, Burd S, Wright C, Lers A, Green PJ. Differential expression of miRNAs and their target genes in senescing leaves and siliques: insights from deep sequencing of small RNAs and cleaved target RNAs. Plant Cell Environ. 2015;38:188–200.
Schommer C, Debernardi JM, Bresso EG, Rodriguez RE, Palatnik JF. Repression of cell proliferation by miR319-regulated TCP4. Mol Plant. 2014;7:1533–44.
Palatnik JF, et al. Sequence and expression differences underlie functional specialization of Arabidopsis microRNAs miR159 and miR319. Dev Cell. 2007;13:115–25.
Schommer C, Bresso EG, Spinelli SV, Palatnik JF. MicroRNAs in Plant Development and Stress Responses (ed Ramanjulu Sunkar). Berlin Heidelberg: Springer; 2012. p. 29–47.
Yoshida A, Ohmori Y, Kitano H, Taguchi-Shiobara F, Hirano HY. Aberrant spikelet and panicle1, encoding a TOPLESS-related transcriptional co-repressor, is involved in the regulation of meristem fate in rice. Plant J. 2012;70:327–39.
Jiang L, et al. DWARF 53 acts as a repressor of strigolactone signalling in rice. Nature. 2013;504:401–5.
Tao Q, et al. The TIE1 transcriptional repressor links TCP transcription factors with TOPLESS/TOPLESS-RELATED corepressors and modulates leaf development in Arabidopsis. Plant Cell. 2013;25:421–37.
Lu X-j, Li K-m, Ye J, Xu R, Huang J. Comprehensive evaluation of 13 newly introduced sugarcane varieties by gray fuzzy analysis. Guangxi Agricultural Sciences. 2008;6:009.
Zhang G-m, Liu H-b, Fang W-k, Chen Q-x. Analysis of main characteristics of new sugarcane parents and their potential utilization in breeding program. Sugar Crops China. 2007;2:002.
Zhou H. Evaluation on cold tolerance of sugarcane varieties under field conditions. Journal of Plant Genetic Resources. 2012;13:968–73.
You Q, et al. Genetic diversity analysis of sugarcane germplasm based on fluorescence-labeled simple sequence repeat markers and a capillary electrophoresis-based genotyping platform. Sugar Tech. 2015;17:1–11.
Ji H, et al. Deep sequencing of RNA from three different extracellular vesicle (EV) subtypes released from the human LIM1863 colon cancer cell line uncovers distinct miRNA-enrichment signatures. PLoS One. 2014;9:e110314.
Zhang G, et al. Deep RNA sequencing at single base-pair resolution reveals high complexity of the rice transcriptome. Genome Res. 2010;20:646–54.
Haas BJ, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.
Li R, et al. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics. 2009;25:1966–7.
Lorenz R, et al. ViennaRNA Package 2.0. Algorithms Mol Biol. 2011;6:26.
Griffiths-Jones S, Saini HK, van Dongen S, Enright AJ. miRBase: tools for microRNA genomics. Nucleic Acids Res. 2008;36:D154–8.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B-Stat Methodol. 1995;57:289–300.
Varkonyi-Gasic E, Wu R, Wood M, Walton EF, Hellens RP. Protocol: a highly sensitive RT-PCR method for detection and quantification of microRNAs. Plant Methods. 2007;3:12.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) method. Methods. 2001;25:402–8.
This work was financially supported by the National Science Foundation of China (31,460,093, 31,400,281), Guangxi Scientific Research and Technology Development Plan (GuiKeAB16380258, Guikegong1598006–1-1D), Nanning Scientific Research and Technology Development Plan (20131038) and the development foundation of Guangxi Academy of Agricultural Science (2015JZ15, 2015JZ91). The Funding “31,460,093” “31,400,281” supported materials cultivation and collection, while “GuiKeAB16380258” “GuiKeAB16380258” “20,131,038” “2015JZ15” “2015JZ91” supported experimental consumption and data analysis consumption.
Availability of data and materials
The raw sequence datasets supporting the results of this article are available in the sequence read archive repository (http://www.ncbi.nlm.nih.gov/sra/) under the accession number SRP116643 (Q1:SRR6051144, Q2: SRR6051145, T1:SRR6051146, T2:SRR6051147, B1:SRR6051142, B2:SRR6051143). Additional supporting tables are included as Additional files.
Ethics approval and consent to participate
Samples in the study (Six sugarcane cultivars) collected from the experimental field of Sugarcane Research Institute in Nanning, Guangxi Province of China, which permits unrestricted use.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Dataset. Tables of Sugarcane conserved and novel miRNAs and Target prediction of differentially expressed miRNAs. (XLS 408 kb)
Summary of supplementary tables and figures. (PDF 1132 kb)
Tables of overlapping miRNAs in six samples. (XLSX 36 kb)
About this article
Cite this article
Li, M., Liang, Z., He, S. et al. Genome-wide identification of leaf abscission associated microRNAs in sugarcane (Saccharum officinarum L.). BMC Genomics 18, 754 (2017). https://doi.org/10.1186/s12864-017-4053-3
- Sugarcane-NGS-smallRNA-leaf absicission associated miRNA