Cold stress responsive microRNAs and their targets in Musa balbisiana

Cold stress is an environmental factor affecting plant development and production. Recently, microRNAs (miRNAs) have been found to be involved in several plant processes such as growth regulation and stress responses. Although miRNAs and their targets have been identified in several banana species, their participation during cold accumulation in banana remains unknown. In this study, two small RNA libraries were generated from micropropagated plantlets of Musa balbisiana grown at normal and low temperature (5°C). A total of 69 known miRNAs and 32 putative novel miRNAs were detected in the libraries by Solexa sequencing. Sixty-four cold-inducible miRNAs were identified through differentially expressed miRNAs analysis. Among 43 miRNAs belonging to 26 conserved miRNA families with altered expression, 18 were upregulated and 25 downregulated under cold stress. Of 21 putative novel miRNAs with altered expression, four were downregulated and 17 upregulated. Furthermore, eight miRNAs were validated by stem-loop qRT-PCR and their dynamic differential expression was analyzed. In addition, 393 target genes of 58 identified cold-responsive miRNAs were predicted and categorized by function. These results provide important information for further characterization and functional analysis of cold-responsive miRNAs in banana.


Introduction
MicroRNAs (miRNAs) are a class of non-coding small RNAs, which post-transcriptionally regulate gene expression in many organisms by cleavage or translation suppression of their target mRNAs [1][2][3][4] . Increasing evidence demonstrates that miRNAs have important roles in many biological and metabolic processes including regulation of plant growth, development and response to biotic and abiotic stresses [5,6] . Cold temperature exposure is a major abiotic stress that affects plant development and growth and significantly impacts crop production. Plants have evolved various processes to cope with cold stress at the cellular and developmental levels. Not unexpectedly, the involvement of miRNAs in response of plants to cold has been evaluated in several plant species such as Arabidopsis thaliana [7][8][9] , poplar [10] , Brachypodium [11] , and rice [12] . In Arabidopsis, the expression of miR165/166, miR319c, miR169, miR393 and miR396 is induced by cold stress [7][8][9] . In Brachypodium, poplar, and Arabidopsis, miR397 and miR169 were found to be upregulated in response to cold stress [11] and miR168 was upregulated in Arabidopsis [9] and poplar [10] , although it was downregulated in rice [12] during cold stress. Recently, an analysis of deep sequencing data from a wheat thermosensitive genic male sterile line showed that 78 unique miRNAs may play a regulatory role in male sterility during cold stress [13] . miR319 has also been shown to regulate cold responses in both sugarcane [14] and poplar [10] .
Banana (Musa spp.), one of the most important staple crops widely cultivated in the tropics and subtropics, is a tropical giant perennial herb belonging to the Musaceae family of the order Zingiberales [15] . Cultivated bananas have evolved from the hybridization of wild species of Musa acuminata (AA) and Musa balbisiana (BB) [16] . Since they are mainly cultivated in tropical and subtropical regions, bananas are quite sensitive to cold. However, cold spells often occur during winter or early spring in banana planting areas in China, which can have a serious impact on banana development, fruit quality and yield. Musa balbisiana, native to South Asia with a geographical distribution ranging from India to south China, is immune to several diseases and pests, such as Fusarium wilt and Black leaf streak disease, and represent excellent sources of natural resistance [17] . The "B" genome is used in breeding programs for its traits of vigor and cold resistance and could also be interesting for tolerance to pests and pathogens [16] . Furthermore, genotypes with the "B" genome are more tolerant to abiotic stresses than those solely based on the "A" genome [18] .
Although miRNAs in bananas were identified in previous studies [15,[19][20][21][22][23] , compared with the number of identified miRNAs in other plants such as rice and maize, the number of known miRNAs in banana is still very low. Moreover, there is no report on cold-responsive miRNAs from M. balbisiana, a wild banana species. In the present study, high-throughput sequencing, which has been widely used for miRNA research, was used to identify conserved and novel miRNA of M. balbisiana. The altered expression levels of these miRNAs under cold treatment were analyzed and compared with controls, and the potential roles of their target genes were investigated. In addition, the expression patterns of several miRNAs were examined by RT-qPCR methods. This is the first study to prove that miRNAs are widely involved in regulating the coldinduced process of banana. These data provide a new insight to cold-tolerance molecular mechanisms of banana under cold stress.

Plant material and stress treatment
Micropropagated plantlets of M. balbisiana (BB) were obtained from the Bioversity International Musa Germplasm Transit Centre (ITC) in Leuven, Belgium [Musa Germplasm Information System (MGIS) accession no. ITC0545]. These plantlets were grown in a growth chamber at 25°C and 80% air humidity under a 16/8 h (light/dark) photoperiod. Five-leaf stage plantlets were used in the study. For the cold-stress treatments, leaf samples were first collected from plants at 25°C (CK), then, the temperature was lowered to 5°C and samples collected at 4, 8, and 24 h. The leaves from these three sampling times were cut into pieces and mixed well (CS). The collected leaf tissues were immediately frozen in liquid nitrogen and stored at -80°C until used.

Small RNA sequencing and bioinformatics analysis
Total RNA was extracted from banana leaves using plant total RNA isolation kit (Tiandz Inc., Beijing, China). RNA degradation and contamination was monitored on 1% agarose gels; RNA purity was checked using a Nano-Photometer spectrophotometer (Implen, Westlake Village, CA, USA) and RNA integrity was assessed using the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). A total of 3 µg RNA per sample was used for the construction of small RNA libraries. Sequencing libraries were generated using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (New England Biolabs, Ipswich, MA, USA) following the manufacture's recommendations and index codes were added to attribute sequences to each sample. Briefly, total RNA was purified and fractionated by 15% TBE urea denaturing polyacrylamide gel electrophoresis (PAGE) and small RNA regions corresponding to the 18-30 nt bands in the marker lane were excised and recovered. Then 18-30 nt small RNAs were ligated to a 5′-adaptor and a 3′adaptor sequentially by T4 RNA ligase. The adapterligated small RNAs were subsequently transcribed into cDNA by SuperScript III Reverse Transcriptase (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA) and then PCR-amplified for 17 cycles using the adaptor primers. The PCR products were purified for highthroughput sequencing. The sequencing was performed by Novogene (Beijing, China).
After sequencing, clean reads were obtained by removing reads containing the following: poly-N, poly A/T/G/C, adapter-contaminated tags, low quality reads and reads less than 18 nt in length from the raw data. All of the high quality sequences were considered significant and analyzed further. The small RNA tags were mapped to a reference sequence by Bowtie (http://banana-genome-hub. southgreen.fr/) [15] . The sequences matching non-coding RNAs (rRNA, scRNA, snoRNA, snRNA and tRNA) deposited in the Rfam 10.1 database were removed from the sequences. The remaining unique sequences that mapped to known mature rice miRNAs in miRBase 21 (http://www.mirbase.org) [24] were considered to be known miRNAs. The remaining unannotated sRNA sequences were analyzed by an integrated combination of miREvo [25] and mirDeep2 [26] software to predict potential novel miRNAs through exploring hairpin structure, Dicer cleavage sites and the minimum free energy. The criteria used for novel miRNA were based on the work of Meyers et al. [27] .

Expression analysis of miRNA under cold stress
The differential expression analyses of miRNA were based on sequence reads generated from cold treatment and control libraries. The expression of miRNAs in the two libraries was normalized to obtain the number of miRNAs per million reads [normalized expression = (number of miRNA reads/total number of clean reads) Â 10 6 ]. Normalized miRNA reads with values less than one were excluded from the differential analysis. Differential expression analysis of two samples was performed using the DEGseq [28] . P-values were adjusted using q-value [29] . q-value < 0.05 and fold change>1 or < -1 were set as the threshold for significant differential expression by default.

Quantitative real-time PCR analysis of miRNA expression
Identified M. balbisiana miRNAs were experimentally validated using quantitative real-time PCR (qRT-PCR). For the cold treatments, the plantlets were treated for 0, 4, 8 and 24 h at 5°C in a growth chamber. Untreated plantlets were used as controls. qRT-PCR were carried out as previously described [20] , with 5S rRNA used as an internal control (endogenous reference). All reactions were run in triplicate for each sample. Validated miRNAs and primers are listed in Table S1(Appendix A). qRT-PCR analysis was performed using FastStart Universal SYBR Green Master (ROX) (Roche Diagnostics GmbH., Mannheim, Germany, Cat. No. 04913850001) on a Mx3005P™ Real-Time PCR System (Stratagene, La Jolla, CA, USA). The relative expression quantification was calculated using the 2 -ΔΔC T method [30] .

Target gene prediction
Target predictions were performed by the online software psRNATarget (http://plantgrn.noble.org/psRNATarget) [31] with default parameters. Gene ontology (GO) enrichment analysis was used for the target gene candidates of differentially expressed miRNAs. GOseq based Wallenius non-central hypergeometric distribution [32] , which could adjust for gene length bias, was implemented for GO enrichment analysis.

Deep sequencing of Musa balbisiana sRNAs
A total of 13151282 (CK) and 17315104 (CS) raw reads were generated from Solexa Sequencing of the two libraries. After processing, 12198214 (CK) and 16738624 (CS) clean reads were obtained. Analysis of these sequences resulted in identification of 914848 (CK) and 12677356 (CS) sequences ranging in size from 18 to 30 nt in length. About 683275 (CK) and 6615966 (CS) sequence reads (about 53% to 60% of 18-30 nt reads), could be mapped to the reference Musa genome [15] ( Table 1). The majority of the small RNAs are 19-22 nt in size in the two libraries, with the 21 (CK) and 20 (CS) nt size classes predominant (Fig. 1).

Identification of known miRNAs in M. balbisiana
To identify the known miRNAs in banana leaves, small RNA sequences were mapped to the known rice miRNAs from the miRBase database. After Blastn search and further sequence analysis, a total of 69 known miRNAs were identified (52 and 68 miRNAs in control and cold treatments, respectively), which belonged to 28 miRNA families (Appendix A, Table S2). Among these miRNAs, 51 known miRNAs (56%) were detected in both libraries, while 17 miRNAs were only detected in the cold treatment library (Fig. 2a). MiR396 was the largest family with seven members, followed by miR166, miR159, miR169 and miR171, with 6, 5, 5 and 5 members, respectively. Of the remaining 23 miRNA families, 10 comprised two to four members, and the others had only one member (Fig. 3). Read numbers based on expression level analysis showed significant differences in known miRNA families, which ranged from one to more than 100000 reads. Five miRNA families (miR156, miR159, miR164, miR166 and miR167) were strongly expressed, with more than 10000 reads detected in the CS library (Appendix A, Table S2).

Identification of novel miRNAs in M. balbisiana
Based on the plant miRNA annotation criteria [27] , the formation of a stable hairpin structure is the primary criterion of a miRNA. In total, 32 potential novel miRNAs were predicted from the two libraries (Appendix A, Table S3). Twenty-one of the 32 novel miRNAs candidates were shared by the two libraries, and 10 predicted miRNAs candidates only existed in the cold treatment library, but only one of the 32 novel miRNA candidates was unique to the control library (Fig. 2b). In this study, the length of these 32 potential novel miRNA precursors ranged from 44 to 273 nt, which was close to the length of miRNA precursors in Cavendish banana fruits [22] . Except for mba-miR19 and mba-miR22, the minimum free energy (MFE) folding value of these miRNA precursors ranged from -92.5.1 kJ$mol -1 to -841.0 kJ$mol -1 , with an average of 282.2 kJ$mol -1 , which was consistent with the previously published low MFE characteristics of miRNA [33] and much higher than the published sequencing data for miRNA precursors from banana fruit (194.6 kJ$mol -1 ) [22] . In addition, Meyers et al. reported the existence of complementary sequences increases the authenticity of predicted novel miRNAs [27] . In our study, 15 of the 32 novel miRNAs have complementary miRNAs. Moreover, all of the novel complementary miRNAs exhibited lower expression levels than the corresponding mature miRNAs (Appendix A, Table S3), which might be due to the rapid degradation of complementary strands during the biogenesis of mature miRNAs [34] .    Tables S2 and S3 (Appendix A). Based on the selected criteria (fold change>1 or < -1 and a P-value < 0.05), 43 known miRNAs and 21 novel miRNAs showed differential expression in response to cold stress (Table 2).
Among the known miRNAs with altered expression, 18 upregulated and 25 downregulated were detected. As for novel miRNAs, 17 upregulated and 4 downregulated were detected in the novel miRNAs. To confirm the expression of identified miRNAs and detect their dynamic response to cold stress at different treatment stages, the expression of six conserved and two novel miRNAs (miR162a, miR167d, miR396g, miR444a, miR535, miR5538, miRn20 and miRn31) with or without significantly altered expression following treatment were analyzed by stem-loop qRT-PCR after 0, 4, 8 and 24 h of cold stress at 5°C. Eight potential miRNAs were expressed and exhibited different expression patterns at different treatment stages (Fig. 4). The eight miRNAs were divided into three groups according to their expression pattern. Group I included three miRNAs (miR396g, miR444a and miR535), and their expression increased rapidly and reached maximum level at 8 h, then decreased. The expression of Group II miRNAs (miR167d, miR5538, miRn20 and miRn31) decreased rapidly and reached a minimum level at 4 h, then increased at 8 h. The expressions of Group I and II were upregulated in 5°C cold stress. During the cold stress, the average expression level of these miRNAs was higher than the control (0 h). Group Ш had only one miRNA (miR162a), and its expression was reduced during the whole period of the cold treatment.
Deep sequencing results showed that miR162a, miR167d, miR535, miR5538, miRn20 and miRn31 were upregulated, whereas miR396g and miR444a were downregulated (Table 2). qRT-PCR analysis indicated that the expression of six miRNAs showed the same variation trend with deep sequencing in spite of variation in the absolute change between the two methods. However, opposite results were detected in miR396g and miR444a. Such inconsistency between deep sequencing and qRT-PCR verification has also been reported by other researchers [35,36] . For the 64 differentially expressed miRNAs, the target genes were predicted in the psRNATarget database to illustrate the putative mRNAs regulated by miRNAs under cold stress. The detailed annotation results of the prediction are shown in the Table S4 (Appendix A). Among these 64 miRNAs, six miRNAs (miR171d, miR398a, miR5083, miRn11, miRn14 and miRn36) had no predicted targets. A total of 393 target mRNAs were predicted from the 58 cold-responsive miRNAs. GO analysis showed that the predicted targets were classified into three main categories: biological processes, cellular components, and molecular functions (Fig. 5). Of these, cellular and metabolic process in biological process, cell and organelle part in cellular component, and binding and catalytic activity in molecular function were the most  established categories. According to the identified miRNA targets results, many were homologous to conserved target genes existing in other plants species: these targets included squamosa promoter binding-like protein (miR156), auxin response factor (miR160 and miRn10), MADS-box transcription factor (miRn17, miR444 and miR5179), MYB family (miRn07 and miR159), AP-2 and ethylene-responsive transcription factor (miR172), basic helix-loop-helix transcription factor (miR5179), WRKY transcription factor (miRn07), NAC transcription factor (miRn20, miRn25, miRn31 and miR164), TCP family transcription factor (miRn14) and homeobox-leucine zipper protein (miRn27). In addition, other predicted target genes are involved in diverse metabolic and physiological processes, including catalytic activities such as dehydrogenase, decarboxylase, helicase, synthase, protein kinases, oxidoreductase, binding protein, lipase, transporters, disease resistance protein, f-box proteins, zinc finger protein, heat shock protein and laccase (Appendix A, Table S4).

Discussion
In recent years, increasing evidence has demonstrated that miRNAs have critical roles in both biotic and abiotic stress responses. However, there has been no systematic identification of cold-responsive miRNAs from banana reported. Taking advantage of the sequenced genome of M. acuminata cv. Pahang (AA) and the high-throughput sequencing technology, the present study reports the first identification of cold-inducible miRNAs in M. balbisiana. We constructed two small RNA libraries and identified 69 known miRNAs and 32 novel miRNAs (Appendix A, Tables S2 and S3). The 69 known miRNAs belonged to 28 conserved miRNA families. Previous studies indicate that 25 conserved miRNA families exist in both monocot and dicotyledonous species [37] . Notably, three miRNA families (miR444, miR5179 and miR5538) are conserved in monocots but not in Arabidopsis or popular, suggesting that these are monocot-specific miRNAs. Except for the miR169 family, the other conserved miRNA families were detected in the two libraries. The 32 novel miRNAs significantly increased the number of miRNAs in M. balbisiana.
Relative expression of the conserved miRNAs varied widely, based on the sequence data in our study. Based on the bioinformatics analysis, 43 conserved miRNAs were identified as cold-responsive ( Table 2). The most striking ones with 8-fold changes were miR156c, miR169e and miR393a. The remaining changes in the expression of miRNAs ranged up to 7-fold ( Table 2). In this study, the expression of 13 conserved miRNAs, including miR164, miR167, miR169, miR171, miR390, miR393, miR394, miR395, miR5179, miR528, miR529, miR535 and miR5538, was upregulated, while miR166, miR396, miR397 and miR444 showed significant downregulation after cold stress. However, miR169 in rice [12] , wheat [38] , poplar [39] and Poncirus [40] was found to be downregulated under cold stress. MiR166 in rice [12] and miR396 in Arabidopsis [7][8][9] and Poncirus [40] were upregulated under cold stress. These results indicate these miRNAs have a species-specific role in response to cold stress in plants. Moreover, we also found miRNA members in the same miRNA family, such as miR160a and miR160e, miR397a and miR397b, and miR398a and miR398b, exhibited reverse expression patterns. This suggests that they may have disparate roles in regulating expressions of the stressresponsive genes or that their downstream target genes are diversified.
The prediction of miRNA target genes is important for understanding miRNA gene-regulation. In this study, targets of 64 cold-responsive miRNAs were predicted (Appendix A, Table S4). The targets of conserved miRNAs, such as miR156, miR171 and miR172, have been investigated in several banana cultivars, and their functions were in accordance with most of the previous studies [20][21][22][23] . For the three conserved and three novel mba-miRNAs, our results revealed that the functions of most target genes are unknown. However, some mba-miRNA could target genes that belong to MYB, bHLH and bZIP transcription factor (TF) families,which have been reported as important regulators of cold stress responses [41][42][43] . As TFs are crucial master switches for biotic and abiotic stress response through controlling a wide range of downstream genes, it is conceivable that the miRNAs may exert an extensive impact on gene expression by regulating these TFs [44] . NACs are plant-specific TFs with diverse roles in development and stress regulation [45] . In this work, NAC was predicted to be the target of miR164, miRn20, miRn25, and miRn31, which was suppressed by cold, implying that the NAC gene may be activated under cold. The enhanced NAC expression may help the plant to fight against the cold stress, as overexpression of a NAC gene, SNAC2, has been reported to confer cold tolerance in rice [46] .
In addition, miRNAs of the same family may target the same genes. For instance, miR396a, miR396e, miR396f and miR396g have the same target genes (growthregulating factor, GRF). It was reported that GRF7 functions as a repressor of a broad range of osmotic stress-responsive genes to prevent growth inhibition under normal conditions [47] . MiR396 was also detected to target a zinc ion binding protein family. Metal ion (copper, zinc, or iron) binding is important in normal respiratory function and the process of protein damage during oxidative stress in A. thaliana [48] . Zinc ion binding increases the structural complexity of abscisic acid stress ripening proteins, which is associated with DNA binding [49] . The functional analysis of these candidate target genes, combined with research on miRNA-mediated post-transcriptional regulation, will finally help in the illustration of the roles of miRNAs under cold stress in banana.

Conclusions
The present study is the first to examine the miRNA expression profile and evaluate miRNA function under cold stress in banana. We identified 69 conserved and 32 putative novel miRNAs from the banana leaf under cold treatment using high-throughput sequencing. The majority of the identified miRNAs were significantly responsive to cold stress in banana. Furthermore, eight miRNAs were validated by stem-loop qRT-PCR and their dynamic differential expressions were analyzed. In addition, 393 target genes of 58 identified cold-responsive miRNAs were predicted and categorized by function. These results provide important information for further characterization and functional analysis of cold-responsive miRNAs in banana.