- Research article
- Open Access
Integrated mRNA and small RNA sequencing reveals microRNA regulatory network associated with internode elongation in sugarcane (Saccharum officinarum L.)
BMC Genomics volume 20, Article number: 817 (2019)
Internode elongation is one of the most important traits in sugarcane because of its relation to crop productivity. Understanding the microRNA (miRNA) and mRNA expression profiles related to sugarcane internode elongation would help develop molecular improvement strategies but they are not yet well-investigated. To identify genes and miRNAs involved in internode elongation, the cDNA and small RNA libraries from the pre-elongation stage (EI), early elongation stage (EII) and rapid elongation stage (EIII) were sequenced and their expression were studied.
Based on the sequencing results, 499,495,518 reads and 80,745 unigenes were identified from stem internodes of sugarcane. The comparisons of EI vs. EII, EI vs. EIII, and EII vs. EIII identified 493, 5035 and 3041 differentially expressed genes, respectively. Further analysis revealed that the differentially expressed genes were enriched in the GO terms oxidoreductase activity and tetrapyrrole binding. KEGG pathway annotation showed significant enrichment in “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction”, which might be participating in internode elongation. miRNA identification showed 241 known miRNAs and 245 novel candidate miRNAs. By pairwise comparison, 11, 42 and 26 differentially expressed miRNAs were identified from EI and EII, EI and EIII, and EII and EIII comparisons, respectively. The target prediction revealed that the genes involved in “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways are targets of the miRNAs. We found that the known miRNAs miR2592-y, miR1520-x, miR390-x, miR5658-x, miR6169-x and miR8154-x were likely regulators of genes with internode elongation in sugarcane.
The results of this study provided a global view of mRNA and miRNA regulation during sugarcane internode elongation. A genetic network of miRNA-mRNA was identified with miRNA-mediated gene expression as a mechanism in sugarcane internode elongation. Such evidence will be valuable for further investigations of the molecular regulatory mechanisms underpinning sugarcane growth and development.
Internode elongation is a major feature that affects plant growth, errectness, biomass and ultimately yield . Thus, the genetics and the regulatory mechanisms of internode elongation in crop plants have been extensively investigated. Genetic and environmental factors such as gene expression [2,3,4], genomic variation [5, 6], hormonal regulation [7, 8], nutrients [9, 10], light , water  and temperature  control internode elongation. Of these regulatory factors, hormonal manipulation is an effective and efficient approach to promote crop growth to promote productivity [14,15,16].
Complex hormonal mechanisms are associated with internode elongation. For example, auxin , gibberellin  and brassinosteroids  induce internode elongation. In contrast, abscisic acid , ethylene  and jasmonic acid  suppress internode elongation in plants. Further, different species and growing condition add additional complexity to growth and developmental regulation. Understanding sugarcane crop-specific mechanisms of hormonal regulation stem growth would help improve crop productivity. Alternatively, endogenous hormones can be manipulated by genome editing [22, 23] and RNA interference technologies [24, 25]. microRNAs (miRNAs), being an effective RNA interference mechanisms, show the prospect of regulating hormone production and action [26,27,28]. TIR1 and AFB, part of auxin signaling, are targets of miR393, and the suppressive effects of miR393 on auxin are indicated in Arabidopsis . GAMYB, a gene in gibberellin signal pathway, is regulated by miR159 . Also, hormones regulate miRNA expression in plants. For instance, with deep sequencing of abscisic acid-treated tomato (Solanum lycopersicum), 269 differentially expressed miRNAs were identified .
Development of sequencing technology has facilitated transcriptome studies that provide unprecedented detail about the molecular biological processes in plants [32, 33]. Transcriptome sequencing approaches promise increased understanding of the expression patterns and molecular regulatory mechanisms in gene expression . By transcriptome sequencing, the genes associated with culm elongation in bamboo (Dendrocalamus sinicus) were identified . In another study, transcriptome sequencing showed the changes in gene expression via induction of ethephon in maize (Zea mays) plants . These studies provide basic information about the functional genes involved in internode elongation. In cotton (Gossypium hirsutum), 64 differentially expressed miRNAs were identified during the fiber elongation process . The miRNA profiles during tissue differentiation and growth revealed by small RNA sequencing may provide new insight for epigenetic regulation, which might determine a starting point toward important questions regarding plant growth.
Sugarcane (Saccharum officinarum L.) is an economically important crop that is widely planted in tropical and subtropical regions . Sugarcane is used for producing ethanol and raw sugar; thus, this valuable crop is grown around the world . Understanding the genetic control of sugarcane growth, particularly the biological process of internode elongation, would accelerate the industrial development of sugarcane cultivation. Investigation of miRNA-mRNA networks in sugarcane could reinforce further crop gains. Although several studies demonstrate changes in gene expression or miRNAs during internode elongation [40, 41], the present study focused on the integrated analysis of miRNA and mRNA interactions, which should produce an image of the mRNA-miRNA networks that occur between transcriptional and posttranscriptional regulation in this biological process.
To better understand the molecular changes and regulation of gene expression by miRNA, we sequenced mRNA and small RNA libraries from internode tissues at different stages including the pre-elongation stage, early elongation stage and rapid elongation stage. The libraries from these tissues were sequenced by an Illumina Hiseq 4000 platform. By comparing the differential expression, the candidate genes and miRNAs involved in internode elongation were identified. Furthermore, the mRNA and miRNA interaction network was built by target prediction using a bioinformatic approach. These integrated mRNA and small RNA sequencing results provide pioneering evidence for a view of candidate internode elongation-associated miRNAs in sugarcane and may be useful for development of potential functional markers to develop molecular breeding.
mRNA expression profiles in different stages from internodes
To understand the molecular mechanism of internode elongation in sugarcane, nine cDNA libraries from the pre-elongation stage (EI), early elongation stage (EII) and rapid elongation stage (EIII) were sequenced. Three biological replicates (individuals) from the stages were included from which a total of 499,495,518 reads were obtained. After filtering, 484,324,322 clean reads were identified (Table 1). All the clean reads were used to perform de novo assembly, which generated 80,745 unigenes. The average length of unigenes was 900 bp, and N50 was 1600 bp.
The gene expression in the present study was calculated by the RPKM method. Hierarchical clustering of all the unigenes showed that the biological replicates from each stage were separately clustered together (Fig. 1a). After normalizing the gene expression of all the unigenes, the results of principal component analysis showed a distinct position of the EI, EII and EIII stages, indicating significant changes among the stages in the transcriptome (Additional file 1). There was greater separation on PC1 for samples from the EII groups, indicating a stronger transcriptional differentiation during early elongation stage (Fig. 1b).
Differentially expressed genes in different stages of internode development
To compare the gene expression in different stages of internode elongation, the RSEM package was used to identify differentially expressed genes. In the comparisons between EI and EII, EI and EIII, and EII and EIII, 493, 5035 and 3041 differentially expressed genes were identified, respectively (Fig. 2a). Between the EI and EII stages, 234 genes were up-regulated and 259 genes were down-regulated. The most differentially expressed genes were in the comparison between EI and EIII stages, which included 1601 up-regulated and 3434 down-regulated genes at EIII stage. 1,061 and 1,980 genes increased and decreased at EII stage, respectively, by compared with EII and EIII stages (Fig. 2b). The Venn plot showed that 17 differentially expressed genes overlapped among the three comparisons (Fig. 2b, Table 2). The overlap between EI vs. EIII and EII vs. EIII had the maximum number of genes (1,483 differentially expressed genes), whereas the overlap of EI vs. EII and EI vs. EIII had the fewest number of genes (221 differentially expressed genes). The overlap between EI vs. EII and EII vs. EIII had 226 differentially expressed genes (Fig. 2c, Additional file 1).
Functional annotation of the differentially expressed genes
To reveal the function of differentially expressed genes in different stages from internodes, GO and KEGG enrichment analyses were performed. Only 3 GO terms were enriched (q-value< 0.05) in the comparison EI vs. EII: two oxidoreductase activity terms and tetrapyrrole binding. It was found that 11 GO terms were enriched (q-value< 0.05) in the EI vs. EIII comparison, which included DNA binding, carboxypeptidase activity, hydrolase activity, oxidoreductase activity, nitrate reductase activity, nucleic acid binding transcription factor activity, tetrapyrrole binding, transmembrane transporter activity, and transporter activity. The comparison between EII and EIII indicated that 29 GO terms were enriched (Fig. 3a, Additional file 2).
From the KEGG enrichment results, it was found that 16 KEGG pathways were enriched. From the EI vs. EII comparison, 9 enriched pathways were identified, including “zeatin biosynthesis” and “nitrogen metabolism”, which involved tissue growth. The EI vs. EIII comparison contained 11 enriched pathways, including “plant hormone signal transduction” and “nitrogen metabolism” associated with growth. For EII vs. EIII, 7 enriched pathways were found, including “nitrogen metabolism” (Fig. 3b, Additional file 3).
The expression profiles of ten candidate genes from “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways were investigated by qPCR. The expression profiles by qPCR were similar to the results of transcriptome analysis (Fig. 4).
Sequencing of small RNAs in internodes
The sequencing of small RNAs was investigated to unveil the dynamic regulation of miRNAs on gene expression during internode elongation in sugarcane. Nine small RNA libraries were sequenced, and a total of 137,610,370 clean reads (Table 3) were generated. The length distribution showed that most of these reads were in the range of miRNA from 18 to 24 nt. After removing the rRNA, snRNA, snoRNA, and tRNA by BLAST against the GenBank and Rfam databases, the remaining small RNAs were retained for the following analysis.
Identification of known and novel miRNAs
The known miRNAs were conducted by blastn to hit the miRBase. A total of 241 known miRNAs were detected in the internode tissues from sugarcane. From all the sequenced libraries, 118 known miRNAs were overlapped in all the groups. miR168-x, miR319-y, miR168-y, miR396-x and miR166-y were the most abundant known miRNAs. The novel miRNAs were predicted by the mireap v0.2 package. A total of 245 novel candidate miRNAs were found in the internodes from sugarcane (Table 3, Additional file 4). Novel-miR0183-5p, novel-miR0209-5p and novel-miR0183-3p were the most abundant novel miRNAs.
Differentially expressed miRNAs and their targets
To understand the miRNA regulatory mechanism in internode elongation, the differentially expressed miRNAs were identified by the TPM method. In the comparisons between EI and EII, EI and EIII, and EII and EIII, 11, 42 and 26 differentially expressed miRNAs were found, respectively. Between stages EI and EII, 4 up regulated and 7 down regulated miRNAs were found. The most extensive differentially expressed miRNAs were found between stages EI and EIII, with 12 and 30 that were up-regulated and down-regulated, respectively. Finally, 5 up-regulated and 21 down-regulated miRNAs were found in the comparison EII vs. EIII (Fig. 5a). In the EII vs. EIII comparison, 14 miRNAs with their 75 targets were identified. The principal component analysis indicated a distinct position of the EI, EII and EIII stages (Fig. 5b, Additional file 5).
Targets of the differentially expressed miRNAs were detected. For 8 of the total differentially expressed miRNAs in the EI and EII comparison, the target unigenes were identified (78 in total), whereas no targets were identified for the other 3 miRNAs. In the EI and EIII comparison, 204 targets for 31 miRNAs were identified (Additional file 6).
Differentially expressed mRNA and miRNA pairs related to internode elongation
Internode elongation is a tissue growth process, and therefore, the mRNA and miRNA network related to tissue growth were identified. As found in the present analysis, “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways were involved in internode elongation. For identification of mRNAs and their corresponding miRNAs in these pathways, 2, 3 and 37 pairs of miRNA-mRNA were found from “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways, respectively (Additional file 7). Negative correlations (rho between − 0.87 and − 1) between these miRNA and mRNA pairs were observed. The expression of these genes and miRNAs was shown in Fig. 6. In the “zeatin biosynthesis” pathway, novel-m0140-5p and novel-m0139-3p targeted cytokinin dehydrogenase 5 precursor (CKX5) and cytokinin hydroxylase-like isoform X1 (CYP735A1), respectively. miR2592-y targeted glutamate dehydrogenase (GDH1) isoform X2 and TPA, and glutamic dehydrogenase1 (GDH1) and novel-m0204-5p targeted ferredoxin-dependent glutamate synthase (GLSF), the chloroplastic precursor in “nitrogen metabolism” pathways. “Plant hormone signal transduction” had the most miRNA and mRNA pairs. Among all miRNA and mRNA pairs in “plant hormone signal transduction” pathways, 11 genes associated with auxin were found under regulation by miRNAs.
To confirm the results of small RNA-seq, qPCR was used to analyze miR2592-y, novel-m0204-5p, miR1520-x and miR6169-x. The expression profiles from qPCR were consistent with the results of the sequencing analysis (Fig. 7).
Stem is the primary valuable tissue for sugarcane, which grows within the internode elongation process. To understand the miRNA-mRNA network during this important process, an integrated mRNA and small RNA sequencing study was done using a next-generation sequencing approach for different internode elongation stages in sugarcane in the present study. This pioneer study provided a molecular basis underlying the elongated internodes in posttranscription regulation.
To investigate the gene profiles in a wide range during sugarcane internode elongation, next-generation sequencing technology was used to analyze the transcriptome changes in the stems from the pre-elongation stage, early elongation stage and rapid elongation stage. The de novo assembly generated 80,745 unigenes. The number of differentially expressed genes in the EI vs. EII, EI vs. EIII and EII vs. EIII comparisons were 493, 5035 and 3041, respectively. Compared with the differentially expressed genes, a significantly higher ratio of unigenes showed similar expression among all the groups, strongly suggesting that the differentially expressed genes were involved in sugarcane internode elongation.
The enriched GO terms identified from the differentially expressed genes were within our expectations. Oxidoreductase activity and tetrapyrrole binding were highly enriched in EI compared with those in the other two groups, indicating the function of these genes in sugarcane internode elongation. Contributing to oxidoreductase activity, flavonoid 3-monooxygenase was significantly higher in the EI pre-elongation stage than in the other groups. Flavonoid 3-monooxygenase is responsible for catalyzing flavonoid to 3′-hydroxyflavonoid. Flavonoids affect plant resistance , and flavonoid accumulation may inhibit lignin synthesis . Tetrapyrrole plays various roles in biological processes in plants, including plant growth . For the tetrapyrrole binding term, higher expression of ent-kaurenoic acid oxidase 1 was identified at EII and EIII than EI. This enzyme catalyzes the gibberellin biosynthesis pathway and regulates lignin synthesis .
To understand the key pathways in sugarcane internode elongation, the enriched KEGG pathways with differentially expressed genes were identified, and “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways were enriched with differentially expressed genes. The differentially expressed genes in “zeatin biosynthesis” pathways included CKX3, CKX5, CKX9, CKX10 and CKO2 at EII and EIII compared with EI. It was noted that with the exception of CKX5 at EII, all the other CKX and CKO were significantly down-regulated at EII and EIII stages. CKX/CKO is responsible for degradation of cytokinin [46, 47]. Low expression of the genes in these pathways might be related to high growth activities at EII and EIII stages during internode elongation via cytokinin accumulation . Glutamate dehydrogenase and nitrate reductase were the two enzyme families in “nitrogen metabolism” pathways that might be involved in internode elongation. Transgenic studies in plants found that glutamate dehydrogenase is important for plant growth and productivity . Nitrate as a nutrient for plant growth is regulated by nitrate reductase, indicating changes in nutrient requirements during sugarcane internode elongation . The changing of “plant hormone signal transduction” during internode elongation was also expected. Auxin-related genes including auxin-responsive proteins were decreased at EII and EIII compared with those at EI. Auxin-responsive proteins affect auxin expression by acting on its promoter [51, 52], leading to regulation of sugar accumulation, ethylene biosynthesis and promoting cell elongation.
The observations that several differentially expressed genes were related to internode elongation in the present study indicated that the sugarcane internode elongation was induced by the changes in gene expression in key pathways, suggesting a complex network of the relationships between hormone secretion and internode elongation in sugarcane.
Previous studies demonstrated that plant miRNAs regulate gene expression at posttranscriptional levels related to plant development, including germination, root development, flowering and internode elongation. Here, small RNA sequencing was performed in the stems from the pre-elongation stage, early elongation stage and rapid elongation stage in the present study. A total of 241 known miRNAs and 245 novel candidate miRNAs were identified in these small RNA libraries. Most of these miRNAs were expressed widely among the three different stages. miR168-x, miR319-y, miR168-y, miR396-x and miR166-y were the most abundant known miRNAs, which were consistent with previous miRNA studies in sugarcane [53,54,55].
Subsequently, the miRNAs involved in control of the gene expression that was related to internode elongation was analyzed. As found in transcriptome analysis, several genes in “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways participated in sugarcane internode elongation. A total of 43 miRNA-mRNA pairs were found included. Only 6 miRNAs were decreased while other 37 miRNAs increased in EII and EIII stages. Thus, the overall expression trend of miRNAs in these pathways were downregulated and their targeted mRNA were upregulated during elongation process. Among them, CKX5 precursor and CYP735A1 were targeted by two novel miRNAs. This evidence showed that numbers of the candidate miRNA-mRNA remained unidentified to date. Further studies are required to demonstrate the substantive regulatory roles of the novel miRNAs. miR2592 is proposed to participate in the adventitious shoot organogenesis in Acacia crassicarpa . In the present study, it was found that miR2592-y was up regulated at EIII and EII stages, whereas its targets were GDH1 isoform X2 and TPA, and GDH1 was down-regulated. In the “plant hormone signal transduction” pathway, 11 miRNA-mRNA pairs were found related to auxin. In these pairs, 5 known miRNAs and 4 novel miRNAs were included. The 5 known miRNAs were miR1520-x, miR390-x, miR5658-x, miR6169-x and miR8154-x. miR1520 is unstable in sugarcane buds under cold stress  and only increased at EII stage, which would be expected to lead to an increase in repressing IAA4-like gene expression. The miR390 and its target as auxin response factor were previously reported , and therefore, the up-regulation at EII and EIII stages to affect the internode elongation was understandable. miR5658 targeted lateral suppressor-like protein (HaLS-L) and participated in growth and development, as reported previously in carrot . Overexpression of miR8154 in Taxus cell lines indicates that miR8154 regulates Taxol, phenylpropanoid, and flavonoid biosynthesis pathways . In this study, high expression of miR8154 at EI stage suggested that the activities of miR8154 were increasing at pre-elongation stage of internode in sugarcane.
With the integrated data from mRNA and small RNA sequencing at different stages of internode elongation, a comprehensive vision could be developed for understanding the miRNA function in this important biological process of sugarcane. Several differentially expressed genes involved in internode elongation related to “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways were identified. The miRNA data were then used to predict the target mRNA and putative mRNA-miRNA interactions that might participate in internode elongation were found (Fig. 8). Therefore, the present findings demonstrated that miRNAs mediated gene expression in sugarcane internode elongation, which indicates a regulatory network between miRNAs and genes at first glance.
In summary, mRNA and small RNA profiles were first revealed in sugarcane internode elongation, and comprehensive analysis was performed to identify the regulatory network during this biological process. The results suggested that potential miRNA-mRNA pairs involved in “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways controlled stem growth and development. This evidence provides valuable information for further functional characterization of the miRNAs and their targets in sugarcane internode elongation.
Plant cultivation and tissue collection
The sugarcane variety GT42, which was bred through sexual hybridization , was used as the material in the present study. The sugarcane was grown in an intelligent greenhouse at the Sugarcane Research Institute of Guangxi Academy of Agricultural Sciences, Nanning, China. The bud body from the upper middle stem was selected on March 15, 2016. The seedcane was cut into single-bud setts, and cultivated in a sand table. When the seedlings grew to a height of 8–10 cm, those with consistent size were transplanted to an intelligent greenhouse. After two more month’s growth, when 9–10 true leaves emerged, the plants entered the pre-elongation stage (EI), and the tissue of the stem internode covered by the second true leaf was collected. When 12–13 true leaves emerged, the plants entered the early elongation stage, and the tissue from the stem internode covered by the second true leaf was collected (EII). The rapid elongation stage was determined when 15–16 true leaves emerged, and the stem internode tissue covered by the second true leaf was collected (EIII). Three biological replicates (individuals) were included in the present study for each stage. These tissues were stored at − 80 °C after liquid nitrogen treatment.
Construction of RNA and small RNA sequencing libraries
According to the supplier’s instructions, total RNA was extracted with RNA Trizol (Invitrogen, Carlsbad, CA, USA). The RNA integrity and quantity were detected by an Agilent 2100 bioanalyzer (Agilent, Santa Clara, CA, USA). Nine mRNA-seq libraries were constructed. First, the mRNAs were enriched by oligo (dT) beads and fragmented into short fragments by fragmentation buffer; second, the first-strand cDNA was reversely transcribed by random primers followed by synthesis of the second-strand cDNA using DNA polymerase I and RNase H; finally, the cDNA fragments were purified with a Qiaquick PCR extraction kit (Qiagen, Hilden, Germany) and ligated with Illumina sequencing adapters (Illumina, San Diego, CA, USA). The mRNA-seq libraries were prepared using pair-end methods (read length was 150 bp) and sequenced using a HiSeq™ 4000 (Illumina, San Diego, CA, USA).
Small RNAs from total RNAs were isolated with 3.5% agarose gel electrophoresis by isolating the 18–30 nt fragments. After ligation with 3′ adapters, the small RNAs were extracted by denaturing urea polyacrylamide gel electrophoresis by isolating the 36–44 nt desired band. Then, 5′ adapters were ligated with the isolated small RNAs following reverse transcription PCR, and the products were isolated again using 3.5% agarose gel electrophoresis at 62–75 bp. The small RNA libraries were sequenced using a HiSeq™ 2000 (Illumina, San Diego, CA, USA).
Raw data processing and de novo assembly of mRNA-seq reads
Raw reads were stored as fastq and filtered to obtain high-quality clean reads. Briefly, the adapters from raw reads were first removed; then, the reads containing > 10% unknown nucleotides (N) were excluded; and finally, low-quality reads that had 50% low quality (q-value lower than 10) bases were removed. The remaining reads were clean reads.
De novo assembly was conducted using the Trinity program (Version: v2.1.1). The inchworm assembled reads using a K-mer-based approach to generate contigs. Chryalis built a de Bruijn graph for clusters of contigs. Finally, based on de Bruijn graphs, Butterfly analyzed the read pairings from the contigs to generate transcripts and unigenes.
Identification of differentially expressed genes
Gene expression was calculated based on Reads Per kb per Million reads (RPKM) . The RPKM formula was the following: RPKM = (1 × 106 × C)/(N × L/1000), where C represents the number of reads mapped to the target unigenes, N represents the number of reads mapped to all the unigenes, and L is the length of the target unigenes. The calculation was performed by the RSEM package. The adjustment of the P-value for multiple testing by FDR correction was performed for identification of differentially expressed genes with FDR cutoff< 0.05. The differentially expressed genes were identified using the standard as |log2Ratio| ≥ 1 and q < 0.05.
The unigenes from De novo assembly were annotated by four public databases, including the NCBI non-redundant protein database (www.ncbi.nlm.nih.gov), Swiss-Prot protein database (www.expasy.ch/sprot), KOG database (www.ncbi.nlm.nih.gov/KOG) and Kyoto Encyclopedia of Genes and Genomes database (www.genome.jp/kegg), using the Blastx package from NCBI (www.ncbi.nlm.nih.gov/BLAST/). The cutoff e-value was 1e-5, and the best-aligned sequence was determined as the annotation for unigenes. Gene Ontology (GO) was used to reflect the gene function of the large unigene set. The GO annotation of the unigenes was performed using Blast2GO software based on the annotation results from the NCBI non-redundant protein database.
Processing of small RNA sequencing data
The clean small RNA reads were obtained from raw reads after removing the adapters and low-quality reads (the reads with more than one base with a Q-value lower than 20). Only small RNAs ranging from 18 to 35 nt in length were included in the further analysis. The clean small RNA reads were matched to the GenBank database (www.ncbi.nlm.nih.gov/genbank/, Version 209.0) and Rfam database (rfam.xfam.org/, Version 11) to identify and remove rRNA, scRNA, snRNA and tRNA by blastn (the cutoff e-value was 1e-5).
Identification of known miRNAs and novel miRNAs
Known miRNAs were identified by hitting the candidate miRNAs miRBase (www.mirbase.org/, Version 21) using blastn (the cutoff e-value was 1e-5). The annotated miRNAs were determined as known miRNAs. Novel miRNAs were predicted by mireap v0.2 (https://sourceforge.net/projects/mireap/). The parameters of mireap were -A 18, −B 25, −a 20, −b 23, −u 20, −e 18, −d 300, −p 16, −v 4, −s 4, and -f 20.
Prediction of the targets of miRNAs
The software patmatch (ftp://ftp.arabidopsis.org/home/tair/Software/Patmatch/, Version1.2) with default parameters was used to predict the targets of miRNAs. Then, based on the patmatch prediction, the gene-miRNA pairs with Spearman’s correlation coefficient for ranked data lower than − 0.5 were further confirmed the targets of miRNAs. The target gene functions of miRNAs were clarified by GO and KEGG annotation. GO enrichment was calculated by FDR correction, with FDR lower than 0.05 as significantly enriched. The KEGG pathways were performed with FDR correction, with FDR lower than 0.05 as a threshold.
Identification of differentially expressed miRNAs
The expression of miRNAs was assessed by transcripts per million (TPM), and TPM = miRNA counts/(Total miRNA counts× 106). The fold-changes were calculated using the formula fold-changes = log2(group1/group2). The miRNAs with fold change greater than 2 and P-value less than 0.05 were identified as differentially expressed miRNAs.
Quantitative PCR (qPCR)
Sugarcane tissues were collected as described previously. The total RNA from the different stages EI, EII and EIII was extracted with RNA Trizol (Invitrogen, Carlsbad, CA, USA). Ten of the differentially expressed genes from “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways were selected for confirmation by qRT-PCR. The primers for qPCR were designed by Primer Express v3.0 (Applied Biosystems, Waltham, MA, USA) (Additional file 8). β-actin was used as the internal control. After RNA extraction, the first-strand cDNA synthesis was obtained using total RNA and a PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China) following the manufacturer’s instructions. The cDNA was diluted to 1:20 and used as template for qPCR, which was performed using SYBR®Premix Ex Taq™ II (TliRNaseH Plus) (TakaRa, Dalian, China) according to the manufacturer’s instructions. The PCR was started with initial denaturation at 95 °C for 30 s, followed by 40 amplification cycles at 95 °C for 10 s and 60 °C for 30 s using an ABI 7500 Fast Real Time PCR System (Applied Biosystems, Waltham, MA, USA). The threshold cycle (CT) values were used to calculate relative expression by the 2-ΔΔCT method , which was performed using 7500 software version 2.0.4 (Applied Biosystems, Waltham, MA, USA). For qPCR of miRNA, the reverse transcriptase reaction was performed using a microRNA RT kit (TakaRa, Dalian, China). The cycles were the same as with mRNA analysis. All the reactions were performed with three biological repeats. The data were presented as the mean ± SD. The significant differences among the groups were assessed using one-factor analysis of variance (one-way ANOVA) followed by a Bonferroni post hoc test. P-value< 0.05 was considered statistically significant.
Availability of data and materials
All the raw data are available in the Sequencing Read Archive (SRA) of NCBI under the BioProject number PRJNA429786.
Cytokinin dehydrogenase 5 precursor
Cytokinin hydroxylase-like isoform X1
Early elongation stage
Rapid elongation stage
False discovery rate
Ferredoxin-dependent glutamate synthase
Lateral suppressor-like protein
Kyoto encyclopedia of genes and genomegenomes
Eukaryotic orthologous group
RNA sequencing or the transcriptome sequencing
Reads Per kb per Million reads
Ribosomal ribonucleic acid
Small nucleolar RNA
Small nuclear ribonucleic acid
Transcripts per million
Lingle SE, Thomson JL. Sugarcane internode composition during crop development. Bioenergy Res. 2012;5:168–78.
Li J, Sima W, Ouyang B, Wang T, Ziaf K, Luo Z, et al. Tomato SlDREB gene restricts leaf expansion and internode elongation by downregulating key genes for gibberellin biosynthesis. J Exp Bot. 2012;63:6407–20.
Zhou H-L, He S-J, Cao Y-R, Chen T, Du B-X, Chu C-C, et al. OsGLU1, a putative membrane-bound endo-1,4-ß-D-glucanase from rice, affects plant internode elongation. Plant Mol Biol. 2006;60:137–51.
Cui K, He C-Y, Zhang J-G, Duan A-G, Zeng Y-F. Temporal and spatial profiling of internode elongation-associated protein expression in rapidly growing culms of bamboo. J Proteome Res. 2012;11:2492–507.
Zhao K, Tung C-W, Eizenga GC, Wright MH, Ali ML, Price AH, et al. Genome-wide association mapping reveals a rich genetic architecture of complex traits in Oryza sativa. Nat Commun. 2011;2:467.
Sripongpangkul K, Posa GBT, Senadhira DW, Brar D, Huang N, Khush GS, et al. Genes/QTLs affecting flood tolerance in rice. Theor Appl Genet. 2000;101:1074–81.
Bishop GJ, Koncz C. Brassinosteroids and plant steroid hormone signaling. Plant Cell. 2002;14:S97–S110.
Hager A. Role of the plasma membrane H +-ATPase in auxin-induced elongation growth: historical and new aspects. J Plant Res. 2003;116:483–505.
La Nafie YA, Santos CBD, Brun FG, Van Katwijk MM, Bouma TJ. Waves and high nutrient loads jointly decrease survival and separately affect morphological and biomechanical properties in the seagrass Zostera noltii. Limnol Oceanogr. 2012;57:1664–72.
Hwang SJ, Hamayun M, Kim HY, Na CI, Kim KU, Shin DH, et al. Effect of nitrogen and silicon nutrition on bioactive gibberellin and growth of rice under field conditions. J Crop Sci Biotechnol. 2008;10:281–6.
Peralta G, Pérez-Lloréns JL, Hernández I, Vergara JJ. Effects of light availability on growth, architecture and nutrient content of the seagrass Zostera noltii Hornem. J Exp Mar Biol Ecol. 2002;269:9–26.
Hattori Y, Miura K, Asano K, Yamamoto E, Mori H, Kitano H, et al. A major QTL confers rapid internode elongation in response to water rise in Deepwater rice. Breed Sci. 2007;57:305–14.
Mazzella MA, Bertero D, Casal JJ. Temperature-dependent internode elongation in vegetative plants of Arabidopsis thaliana lacking phytochrome B and cryptochrome 1. Planta. 2000;210:497–501.
Cleland RE. Auxin and cell elongation. In: Davies PJ, editor. Plant hormones: biosynthesis, signal transduction, action! London: Springer; 2010. p. 204–20.
Davies PJ. The plant hormones: their nature, occurrence, and function. In: Davies PJ, editor. Plant hormones: biosynthesis, signal transduction, action! Dordrecht: Springer; 2010. p. 1–15.
Davies PJ. Plant hormones and their role in plant growth and development. Berlin/Heidelberg: Springer Science & Business Media Press; 2012.
Ingram TJ, Reid JB, MacMillan J. The quantitative relationship between gibberellin A1 and internode growth in Pisum sativum L. Planta. 1986;168:414–20.
Clouse SD, Sasse JM. BRASSINOSTEROIDS: essential regulators of plant growth and development. Annu Rev Plant Physiol Plant Mol Biol. 1998;49:427–51.
Azuma T, Hirano T, Deki Y, Uchida N, Yasuda T, Yamaguchi T. Involvement of the decrease in levels of abscisic acid in the internodal elongation of submerged floating rice. J Plant Physiol. 1995;146:323–8.
Metraux JP, Kende H. The role of ethylene in the growth response of submerged deep water rice. Plant Physiol. 1983;72:441–6.
Heinrich M, Hettenhausen C, Lange T, Wünsche H, Fang J, Baldwin IT, et al. High levels of jasmonic acid antagonize the biosynthesis of gibberellins and inhibit the growth of Nicotiana attenuatastems. Plant J. 2013;73:591–606.
Pieterse CMJ, Van der Does D, Zamioudis C, Leon-Reyes A, Van Wees SCM. Hormonal modulation of plant immunity. Annu Rev Cell Dev Biol. 2012;28:489–521.
Costantini E, Landi L, Silvestroni O, Pandolfini T, Spena A, Mezzetti B. Auxin synthesis-encoding transgene enhances grape fecundity. Plant Physiol. 2007;143:1689–94.
Schijlen EGWM, de Vos CHR, Martens S, Jonker HH, Rosin FM, Molthoff JW, et al. RNA interference silencing of chalcone synthase, the first step in the flavonoid biosynthesis pathway, leads to parthenocarpic tomato fruits. Plant Physiol. 2007;144:1520–30.
Wasson AP, Pellerone FI, Mathesius U. Silencing the flavonoid pathway in Medicago truncatula inhibits root nodule formation and prevents auxin transport regulation by rhizobia. Plant Cell. 2006;18:1617–29.
Denance N, Sanchez-Vallet A, Goffner D, Molina A. Disease resistance or growth: the role of plant hormones in balancing immune responses and fitness costs. Front Plant Sci. 2013;4:155.
Guo HS, Xie Q, Fei J-F, Chua N-H. MicroRNA directs mRNA cleavage of the transcription factor NAC1 to downregulate auxin signals for Arabidopsis lateral root development. Plant Cell. 2005;17:1376–86.
Jones-Rhoades MW, Bartel DP. Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol Cell. 2004;14:787–99.
Si-Ammour A, Windels D, Arn-Bouldoires E, Kutter C, Ailhas J, Meins F, et al. miR393 and secondary siRNAs regulate expression of the TIR1/AFB2 auxin receptor clade and auxin-related development of Arabidopsis leaves. Plant Physiol. 2011;157:683–91.
Tsuji H, Aya K, Ueguchi-Tanaka M, Shimada Y, Nakazono M, Watanabe R, et al. GAMYB controls different sets of genes and is differentially regulated by microRNA in aleurone cells and anthers. Plant J. 2006;47:427–44.
Cheng HY, Wang Y, Tao X, Fan YF, Dai Y, Yang H, et al. Genomic profiling of exogenous abscisic acid-responsive microRNAs in tomato (Solanum lycopersicum). BMC Genomics. 2016;17:423.
Edwards D, Batley J. Plant genome sequencing: applications for crop improvement. Plant Biotechnol J. 2010;8:2–9.
Varshney RK, Nayak SN, May GD, Jackson SA. Next-generation sequencing technologies and their implications for crop genetics and breeding. Trends Biotechnol. 2009;27:522–30.
Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10:57–63.
Cui K, Wang H, Liao S, Tang Q, Li L, Cui Y, et al. Transcriptome sequencing and analysis for culm elongation of the world’s largest bamboo (Dendrocalamus sinicus). PLoS One. 2016;11:e0157362.
Wei X, Zhang W, Zhang Q, Sun P, Li Z, Zhang M, et al. Analysis of differential expression of genes induced by ethephon in elongating internodes of maize plants. Front Agric Sci Eng. 2016;3:263–82.
Wang Y, Ding Y, Liu J-Y. Identification and profiling of microRNAs expressed in elongating cotton fibers using small RNA deep sequencing. Front Plant Sci. 2016;7:1722.
Tew TL, Cobill RM. Genetic improvement of sugarcane (Saccharum spp.) as an energy crop. In: Vermerris W, editor. Genetic improvement of bioenergy crops. New York: Springer; 2008. p. 273–94.
Nonato RV, Mantelatto PE, Rossell CE. Integrated production of biodegradable plastic, sugar and ethanol. Appl Microbiol Biotechnol. 2001;57:1–5.
Peng Z, Zhang C, Zhang Y, Hu T, Mu S, Li X, et al. Transcriptome sequencing and analysis of the fast growing shoots of moso bamboo (Phyllostachys edulis). PLoS One. 2013;8:e78944.
Bosch M, Mayer C-D, Cookson A, Donnison IS. Identification of genes involved in cell wall biogenesis in grasses by differential gene expression profiling of elongating and non-elongating maize internodes. J Exp Bot. 2011;62:3545–61.
Havsteen BH. The biochemistry and medical significance of the flavonoids. Pharmacol Ther. 2002;96:67–202.
Besseau S, Hoffmann L, Geoffroy P, Lapierre C, Pollet B, Legrand M. Flavonoid accumulation in Arabidopsis repressed in lignin synthesis affects auxin transport and plant growth. Plant Cell. 2007;19:148–62.
Tanaka R, Tanaka A. Tetrapyrrole biosynthesis in higher plants. Annu Rev Plant Biol. 2007;58:321–46.
Helliwell CA, Chandler PM, Poole A, Dennis ES, Peacock WJ. The CYP88A cytochrome P450, ent-kaurenoic acid oxidase, catalyzes three steps of the gibberellin biosynthesis pathway. Proc Natl Acad Sci U S A. 2001;98:2065–70.
Schmülling T, Werner T, Riefler M, Krupková E, Manns IB. Structure and function of cytokinin oxidase/dehydrogenase genes of maize, rice, Arabidopsis and other species. J Plant Res. 2003;116:241–52.
Frebort I, Kowalska M, Hluska T, Frebortova J, Galuszka P. Evolution of cytokinin biosynthesis and degradation. J Exp Bot. 2011;62:2431–52.
Caboni E, D'Angeli S, Chiappetta A, Innocenti AM, van Onckelen H, Damiano C. Adventitious shoot regeneration from vegetative shoot apices in pear and putative role of cytokinin accumulation in the morphogenetic process. Plant Cell Tissue Organ Cult. 2002;70:199–206.
Dubois F, Tercé-Laforgue T, Gonzalez-Moro M-B, Estavillo J-M, Sangwan R, Gallais A, et al. Glutamate dehydrogenase in plants: is there a new story for an old enzyme? Plant Physiol Biochem. 2003;41:565–76.
Crawford NM. Nitrate: nutrient and signal for plant growth. Plant Cell. 1995;7:859–68.
Hagen G, Guilfoyle T. Auxin-responsive gene expression: genes, promoters and regulatory factors. Plant Mol Biol. 2002;49:373–85.
Tiwari SB, Hagen G, Guilfoyle T. The roles of auxin response factor domains in auxin-responsive transcription. Plant Cell. 2003;15:533–43.
Bottino MC, Rosario S, Grativol C, Thiebaut F, Rojas CA, Farrineli L, et al. High-throughput sequencing of small RNA transcriptome reveals salt stress regulated microRNAs in sugarcane. PLoS One. 2013;8:e59423.
Gentile A, Ferreira TH, Mattos RS, Dias LI, Hoshino AA, Carneiro MS, et al. Effects of drought on the microtranscriptome of field-grown sugarcane plants. Planta. 2013;237:783–98.
Ortiz-Morea FA, Vicentini R, Silva GFF, Silva EM, Carrer H, Rodrigues AP, 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.
Liu W, Yu W, Hou L, Wang X, Zheng F, Wang W, et al. Analysis of miRNAs and their targets during adventitious shoot organogenesis of Acacia crassicarpa. PLoS One. 2014;9:e93438.
Yang Y, Zhang X, Chen Y, Guo J, Ling H, Gao S, et al. Selection of reference genes for normalization of microRNA expression by RT-qPCR in sugarcane buds under cold stress. Front Plant Sci. 2016;7:86.
Marin E, Jouannet V, Herz A, Lokerse AS, Weijers D, Vaucheret H, et al. miR390, Arabidopsis TAS3 tasiRNAs, and their auxin response factor targets define an autoregulatory network quantitatively regulating lateral root growth. Plant Cell. 2010;22:1104–17.
Sarangzai AM. Profiling the carrot (Daucus carota L.) microRNAs and their targets. Pak J Bot. 2013;45:353–8.
Zhang M, Dong Y, Nie L, Lu M, Fu C, Yu L. High-throughput sequencing reveals miRNA effects on the primary and secondary production properties in long-term subcultured Taxus cells. Front Plant Sci. 2015;6:604.
Wang LW, Liao JX, Tan F, Tang SY, Huang JY, Li X, et al. Breeding of new high-yield, high-sugar and lodging-resistant sugarcane variety Guitang 42 and its high-yield cultivation technique. J Southern Agric. 2015;46:1361–6.
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5:621–8.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25:402–8.
The authors are grateful to Dr. Prakash Lakshmanan (E-Mail: firstname.lastname@example.org) for his critical revision of this manuscript. We thank Pu Li from Guangxi Pfomic Info Technology Co., Ltd. for the technical support of bioinformatics analysis.
National Natural Science Foundation of China (31701363 and 31360312), Guangxi Natural Science Foundation Project (2018GXNSFAA138149, 2017GXNSFBA198050, 2016GXNSFBA380034, 2015GXNSFDA139011 and 2015GXNSFBA139095), Guangxi Key Laboratory Construction Project (15–140-13 and 16–380-18), Guangxi Science and Technology Project (Guike AA17202042–13), Guangxi Bangui Scholars and Special Experts Special Funds (2013), Guangxi Academy of Agricultural Sciences Fund Project (2015YM13 2018YM02 and 2018YT01), and National Modern Agricultural History Technology System, Guangxi Sugarcane Innovation Team Project (nycytxgxcxtd-03-01). The funding bodies had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Ethics approval and consent to participate
No special permits were required to collect tissues from sugarcane.
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.
The RPKM of genes by RNA-Seq.
GO analysis of differentially expressed genes.
KEGG pathway of differentially expressed genes.
Novel miRNAs predicted by mireap v0.2.
miRNA TPM values.
miRNA target predictions.
The differentially expressed mRNA and miRNA pairs in “zeatin biosynthesis”, “nitrogen metabolism” and “plant hormone signal transduction” pathways.
Primers used in this study.
About this article
Cite this article
Qiu, L., Chen, R., Fan, Y. et al. Integrated mRNA and small RNA sequencing reveals microRNA regulatory network associated with internode elongation in sugarcane (Saccharum officinarum L.). BMC Genomics 20, 817 (2019). https://doi.org/10.1186/s12864-019-6201-4
- Next-generation sequencing
- Zeatin biosynthesis
- Nitrogen metabolism
- Plant hormone signal transduction