Integrated mRNA and miRNA transcriptome analysis on the regulational of mechanism tuber expansion in Yam

Background Tuber is the storage organ of yam derived from modified stems. The development of tube expansion is a complex process and depends on the expression of genes that are involved by environmental and endogenous factors. However, little is known about the mechanism that regulates tuber expansion. In this study, in order to identify genes and miRNAs involvd in tuber enlargement, we examined the transcriptomes and small RNA in tuber initiation and expansion stages. Results A total of 14238 transcripts expressed differentially in the expansion stage were firstly identified using transcriptome technology. Of these, 5723 and 8515 genes were up and down regulated, respectively. The functional analysis revealed the coordination of tuber plant involved in processes of cell events, metabolism, biosynthesis and signalling transduction pathways at the transcriptional level, suggesting these DEGs are somehow involved in the responses of tuber expansion. In additional, 536 transcription factor genes were found differently expression in expansion stage at the transcript level. REVEILLE 6-like is identified to be up regulated in circadian rhythm pathway. REVEILLE 6-like LHY protein, zinc finger CCCH domain-containing protein 14, and DELLA genes were involved up regulated in expansion stage. Moreover, these genes were respectively involved in circadian rhythm pathway, starch and sucrose metabolism pathway, and GA pathway by KEGG analysis. Noteworthy, the analysis of the data showed that there were 23 known tuber miRNAs belonged to 11 miRNA families, and 50 novel miRNAs, miRNA160 miRNA396, and miR535 may be involved in complexity network to regulate cell division and differentiation in the expansion stage in yam. The integrated analysis of miRNA-mRNA were identified to be preferentially expressed in hormone signalling in expansion stage, like the ARFs by miRNA160 highlighting the involvement of miRNA-mRNA in the regulation hormone of tuber expansion in yam. The transcriptome and miRNA datasets presented here identified a subset of candidate genes and detect differential miRNA expressions (DEMs). The results showed that miRNAs expression was dynamically regulated in expansion stages. A total of 44 differentially expressed miRNAs were identified including 11 known and 33 novel miRNAs, It showed that, 40 were down regulated, and 4 were up regulated in the expansion stage(Table 3). It is interesting that miR160

miRNAs putatively associated with tuber expansion in yam, thus a hypothetical model of genetic regulatory network associated with tuber expansion in yam was put forward, which may provide a foundation for molecular regulatory mechanism researching on tuber expansion in Dioscorea species.

Background
Yams are monocotyledonous plant belonging to the family of Dioscoreaceae, and the tuber is its harvested organ. Tuber originates from the expansion of underground stem, which is suitable for storage of nutrients, with lots of large parenchyma cells. Tuber morphogenesis, which involves three universal process: induction, initiation and formation; starch and stored protein of accumulation are two main processes in tuber growth [ 1]. The tuber morphogenesis of yam could be divided into three periods: initiation stage, expansion stage, and maturation stage, and the expansion stage could be divided into three periods: early expansion stage, middle expansion stage, late expansion stage [ 2,3].Tuber morphogenesis is a complex physiological process regulated by hereditary, environment, hormones,etc [ 4]. Great efforts have been made to the explore the physiological factors affecting tuber morphogenesis of yam, especially the initiation stage developing into expansion stage, short days promoting tuber expansion at the initiation tuber stage of water yam (D. alata), whereas long days inhibiting tuber expansion [ 5,6]. Endogenous hormones played an important role in tuber morphogenesis, for example gibberellins(GA), acetic acid(IAA) and abscisic acid(ABA) playing a key role at the beginning of tuber expansion stage, and trans-zeatin(tZ), jasmonic acid(JA) were also involved in tuber expansion [ 2,7,8].To date, exogenous hormones have been used to study tuber expansion, GAs could promote tuber expansion and yield by treatment in vitro and in vivo [ 9,10].Exogenous GA combined with ABA promoted microtuber growth and expansion [ 11].Exogenous JA was found to be essential for yam tuberization and induced an increase of tuber number in vitro and in vivo [ 12,13]. However, fundamental knowledge about endogenous metabolic networks is utterly lacking in tuber expansion.
The induction and growth of microtubers in vitro in yam have been found to be under the control of nutrients, sucrose concentration was the most important factor affecting tuberization and frequency of proliferation [ 7,14]. Yam tuber morphology were significantly correlated with nutrient accumulation and enzymatic activity, sucrose, soluble sugars,proteins are significantly increased in tuber expansion stage, subsequently decreased in maturity stage, starch content was increased in the whole tuber morphogenesis, and sucrose synthase, sucrose phosphate synthase, and AGPase were significantly correlated with these nutrient accumulation [ 15]. Although many DNA molecular markers were used to uncover the genetic diversity and relationship among yam germplasms [ [16][17][18], it is little known about specific genes involving tuber morphogenesis. The sucrose synthase 4 and sucrose-phosphate synthase 1 were associated with the earliest stages of starch biosynthesis and storage, a SCARECROW-LIKE gene was involved in formation of adventitious roots [ 19]. PE2.1 and PE53 are the members of pectinesterase(PE) superfamily, which are likely to be associated with the regulation of starch and sucrose metabolism and signalling pathways. Therefore, they may play important roles in microtuber formation [ 20]. The tuber morphogenesis is a complex biological process involving many specific genes and proteins, especially yam tuber expansion. Transcriptome techniques can efficiently find and detect these genes and proteins. Potato is a tuber crop, many transcriptome analysis revealed that there were many genes were regulated in early stages of stolon to tuber transitions or tuberization by nutrients, photoperiodic conditions, exogenous hormones and stress [21][22][23][24]. Former transcriptome study revealed that some putative gens were involved in dioscin biosynthesis [ 25 ], and chalcone isomerase (CHS), flavanone 3-hydroxylase (F3H), flavonoid 3'monooxygenase(F3'H), dihydroflavonol 4-reductase(DFR), leucoanthocyanidin dioxygenase(LDOX), and flavonol 3-O-glucosyltransferase(UF3GT) were significantly expressed in flavonoid biosynthesis [ 26 ]. However, there are not the reports on the transcriptome study on tuber expansion. Thus, to investigate the key genes involved in the metabolism of tuber expansion in yam will provide a clear network to elucidate the molecular mechanism of tuber growth by transcriptome. microRNAs (miRNAs) are endogenous small non-coding RNA with important functions in many biological processes, such as the regulations of growth and development, stress response, and metabolism. Many studies have indicated that miRNAs play important roles in root and tuber formation or development [27][28][29].MiR165/166 regulated root growth through determining the fate of root cells combined with phytohormone crosstalk in Arabidopsis by negative regulating its target genes auxin response factor ARF10, ARF16 and ARF17 [ 30]. miRNA 172 and miR156 were involved in tuberization process, acting either as a component or a regulator of long distance gibberellin signaling pathways [ 31,32]. Potato specific miRNA193, miRNA152, and the conserved miR172-1, miRNA172-5 showed significant expression during developmental stages of tuberization [ 28 ]. Although a number of studies have been identified that miRNAs are involved in tuber and root development, the miRNA-mediated regulatory network during tuber expansion is still unclear.
Although the whole genome sequencing of the heterozygous diploid Guinea yam(D. rotundata) had been done for sex determination at the seedling stage [ 33 ], a more detailed comparative transcriptome and miRNA analysis in the expansion stage should be detected. In this study, both transcriptome and small RNA sequencing were carried out to thoroughly investigate the crucial genes and miRNA in tuber expansion stages in yam.

Results
Overview of transcriptome dynamics and small RNA sequencing In order to identify transcriptome and miRNAs co-regulated regulation in tuber expansion, we examined the transcriptomes and small RNA in tuber initiation and expansion stages (Fig.1).
Meanwhile, the transcriptome library was constructed from mixed RNA pools consisting of initiation and expansion stages in order to construct small RNA and transcriptomes. Approximately 6.67Gb bases in total on BGISEQ-500 sequencing platform were gained by transcriptome de novo analysis(Additional file 1: Table S1). After filtering low-quality reads and adaptor sequences, 7.7% of the total unigenes longer than 2000bp in the length of distribution of all the assembled yams transcripts were gained, which indicated that it is difficult to calculate the number of genes in yams due to the deficiency information of its genome dates (Fig.2) (Fig.4), the species distributions of the best matches for each sequence matched the genes of Elaeis guineensis (33.21%), followed by Phoenix dactylifera (23.35%) , Ananas comosus (7.47%), Musa acuminate subsp.malaccensis (7.31%).
A total of 32026 genes were detected, accounting for approximately 6.6 G of clean reads of each mRNA library by RNA-seq analysis(Additional file 1: Table S1). The average mapping rate of transcriptome library was 82.57%(Additional file 1: Table S1). Both the tuber initiation and expansion stages samples showed good correlations in each other based on the Pheatmap cluster analysis (Fig.5).

Differentially expressed genes annotation by GO term and KEGG pathway
To identify the genes regulated differentially in tuber expansion stage, DEseq software was used to compare the gene expression between initiation and expansion stages. Of these, 5765 and 8515 genes were up and down regulated, respectively(Additional file 1: Table S2).
To better understand the function of the DEGs, 44 significant GO categories were identified. For cellular component, 15 GO categories including cell, membrane and membrane part, organelle, and organelle part were enriched among DEGs (Fig.6). For biological processes, DEGs related to cellular process, metabolic process and biological regulation were enriched in expansion stage (Fig.6). The molecular functions of the DEGs were mainly associated with catalytic activity, binding, transporter activity, structural molecular activity. Among the significant GO term anaylsis, 15 genes were enriched in cell wall polysaccharide metabolic process, 15 genes were involved in hemicellulose metabolic process, and 13 genes were related to xyloglucan metabolic process related to cell wall formation in expansion stage (Table 1). Moreover, the results also showed a number of significantly expression genes involved in tissue development, root morphogenesis, root system development, root development (Table 2). KEGG pathways analysis revealed that plant hormone signal transduction,biosynthesis of amino acids were enriched with DEGs in expansion stage (Fig .7,Fig.8).
The other pathways like MAPK signalling pathway, starch and sucrose metabolism and carbon metabolism were also expressed in KEGG pathway. These metabolic pathways may be closely related to the development of tuber expansion and bioactive compound synthesis.  Table S3). Among them, the majority of SPS were up-regulated, and SuSy were down-regulated in expansion stage. Interestingly, dioscorins, the major storage proteins in yam tubers, were significantly up-regulated in expansion stage(Additional file 1: Table S3). These results indicated that many functional genes were involved in expansion stage in yam.

Transcription factor analysis
Transcription factors are widely involved in various biological processes. 1254 TF-encoding genes were found and classified into 55 different families in yam tuber, of which MYB, MYB-related and AP2-EREP were enriched in tuber (Fig.10). Further analysis revealed that 541TF-encoding genes belong to 48 TF families were differentially expressed between initiation and expansion stage. There were 286 and 255 TF-encoding genes were up and down regulated, respectively(Additional file 1: Table S4). The TF gene families with the highest number were express by heat map in expansion stage (Fig.11).

Detection of known and novel miRNAs expressed in tuber initiation and expansion stages
The investigation of both known miRNA and novel putative miRNAs was performed by miRDeep2 program. This program combined the position and frequency of small RNAs with the secondary structure of miRNA precursors to provide novel miRNA which may be specifically in tuber. Data analysis showed that there were 23 known miRNAs(21 and 20 in tuber initiation and expansion stage, respectively ) and 50 novel miRNAs in yam tuber (Additional file 1: Table S5).
Further analysis revealed that 23 known miRNAs belonged to 11 miRNA families, miRNA168, miRNA159 and novel miRNA16 were the most extensively represented families. All miRNAs were analyzed to detect differential miRNA expressions (DEMs). The results showed that miRNAs expression was dynamically regulated in expansion stages. A total of 44 differentially expressed miRNAs were identified including 11 known and 33 novel miRNAs, It showed that, 40 were down regulated, and 4 were up regulated in the expansion stage (Table 3). It is interesting that miR160 family(miRNA160,miRNA160a-5p, miRNA160b-1, miRNA160h-1) and miR535a_1 were down-regulated in expansion stage, miRNA396b and miRNA168 were up-regulated compared to their miRNA families( Table 3).

Identification of target genes of differentially expressed miRNA in tuber expansion compared with initiation stage
The miRNA regulates target mRNA through translational repression or mRNA degradation. To identify the correlation between the expression of DEMs and DEGs, a total of 11 DEMs were putatively targeted 34 DEGs(Additional file 1: Table S6). Furthermore, based on GO and KEGG analysis of the targets, it revealed that several key genes involved in expansion stage were completely regulated by miRNAs. GO analysis revealed differentially expressed miRNAs regulated in expansion stage by binding, cell, membrane, organelle, membrane part( Fig.12). KEGG pathways analysis were involved in transport and catabolism, signal transduction, biosynthesis of other secondary metabolites, energy metabolism, carbohydrate metabolism, lipid metabolism. It is notable that 3 genes which were targets of 2 differentailly expressed miRNAs(miR535a_1, miR160) were identified in the pathway of plant hormone signal transduction (Fig.13). miRNA160 was down regulated in expansion stage, ARF18(CL2135.Contig1_Total_1), ARF17(Unigene5660_Total_1) were up regulated in auxin signal transduction and Log2(GH16-E/GH16-I) of three miRNAs and their target genes showed not only all of them were differentially expressed significantly but also these miRNAs regulated the expression of their targets in expansion stage. miR535a_1 was involved in brassinosteroid signal transduction pathway by targeting SBP. Moreover, miR396b and miR396a-5b were involved in other metabolism processes and functional proteins by targeting genes, including IST1-like protein, growth-regulating factor 4, photosystem II oxygen-evolving enhancer protein 2, acetyl-CoA carboxylase biotin carboxyl carrier protein (Fig.13, Additional file 1: Table S6).
Some target genes and corresponding miRNA were involved in cell event, regulation of plant hormone levels, carbohydrate metabolism, and other metabolism, implying that they may play crucial roles during the tuber expansion stage (Fig .14).

Validation of the DEGs and DEMs data using RT-qPCR
To identify the accuracy and reliability of transcriptome and miRNA data, the RT-qPCR was used to measure the expressions of some DEGs and DEMs, including 17 mRNAs, and 11 miRNAs, and specific primers were designed (Additional file 1: Tables S6 and S7). All of mRNA and miRNA expressions were confirmed to the accuracy of the RNA-seq data. Overall, these results indicated that 16 out of 17 mRNAs showed similar expression patterns compared to DEGs analysis (Fig.15), and 9 out of 11 miRNAs also showed very similar patterns compared to DEMs analysis (Fig.16).

Discussion
In stem tuber plant, the formation of tuber is based on successive gene expression during  Table S3). The role of CBL in the root development displayed a strong gene expression in the cell proliferation in root, its expression in root was regulated by auxin and cytokinin [ 47,48]. These results demonstrated that calcium signal is regulated in expansion stage. In addition, MAPK signal pathways are known to play a central role in cell proliferation, differentiation and hormone, but it is not known to be involved in tuber or root formation [ 49]. In this study, MAPK signal related genes were up-regulated in expansion stage(Additional file 1: Table S3) 57]. In this study, a lot of transcription factor genes were found differently expression in expansion stage at the transcript level, including 66MYB, 64AP2-EREBP, 52bHLH, 37WRKY (Fig.10), and REVEILLE 6-like were identified to up regulated in circadian rhythm pathway. Zinc finger CCCH domain-containing protein 14 is involved up regulation in starch and sucrose metabolism pathway (Fig.10) ,which play a variety of important roles in growth and development, hormone response, and response to biotic and biotic stresses [ 58 ]. 3 DELLA genes were up regulated in GA pathway (Fig.10), which may play key regulatory functions in expansion stage, and are involved in cell division and expansion [ 57]. All the data suggest that these transcription factors may be potentially involved in the expansion stage. ]. In KEGG pathway annotation, the hormone signal pathway was the most enrich one, which was involved in 8 signal pathways (Fig.8). In radish, most Aux/IAA ,ARFs and SAUR genes were abundant expressed in the root cortex splitting and expanding stage, implying that these transcripts may be involved in cell expansion in the cambium[ 35].Among auxin signal related genes, Aux/IAA, ARFs were highly expressed in expansion stage in this study(Additional file 1: Table S2), and auxin has been identified to regulate cell division and expansion by changing genes expression [ 1,61], implying that Aux/IAA, ARFs may be provided excellent candidates in expansion stage.
Interestingly, in this study, DELLA genes were up regulated, while GIDI, like GID1 and GID2 were down regulated in expansion stage(Additional file 1: Table S3) ,implying they may play a major in the tuber expansion stage in yam. Yam tuber morphology were significantly correlated with sucrose and starch in tuber expansion stage [ 15]. Therefore, these starch and sucrose metabolism genes were required for the tuber expansion development.
The regulation miRNAs by targeting potential genes miRNA mediated gene regulation has been extensively studied in root and tuber development by transcriptional and posttranscriptional levels, which is better understanding of the molecular regulatory networks in tuber expansion stage. In this study, majority of miRNA were down regulated in expansion stage (Table 3). In maize leaf, miR160 were significantly up regulated in the meristem relative to the elongation and mature zones [ 63 ], however, miR160 was down regulated in yam tuber expansion stage (Table 3). In general, miR160 expression forms are differently in different plant. In this study, some miRNA-mRNA pairs were observed in expansion stage (Fig.13) For the small RNA libraries, RNA bands of around 18-30nt in length were isolated. Libraries were prepared according to the small sample preparation protocol, and were sequenced with BGISEQ-500 sequencing platform. The raw reads from mRNA and small RNA were subjected to quality control(QC) to harvest high quality of sequencing data. Raw reads were further filtered by a few of steps to obtain clean reads.
Processing and mapping of mRNA-Seq and small RNA libraries sequencing data After filtering, the clean reads were performed de novo assembly to obtain unigenes, and functional databases are NT, NR, GO, KOG, KEGG, SwissProt and InterPro used to annotate unigene function, Blastn, Blastx, Diamond ,Blast2GO and InterProScan5 were aligned unigenes to these functional data.
RNA-seq reads were mapped to the transcriptome de novo analysis data using bowtie2 results with default parameters after preprocessing mRNA-seq data [ 72]. Gene expression levels were presented as FPKM values [ 73], genes with expression levels >1 FPKM were retained for further analysis. Small RNA reads also were screened from raw sequencing reads by removing adaptors, poly A sequences, and low quality bases, sequences shorter than 18nt or longer than 32nt, after trimming, were removed.
The high quality clean reads were mapped to the reference genome and to other sRNA databases by Bowtie, Cmsearch and Rfam [ 74].
Identification of known and novel miRNAs, and prediction of the targets of miRNAs To identify the isoforms of known miRNA, miRProf tool was run to with clean sRNA reads, and the conserved miRNAs were identified with the known plant miRNAs registered in miRBase [ 75]. To identify novel miRNAs and their precursors, the unique reads were submitted to the RT-qPCR were performed using primers (Table S6, S7) on real-time PCR detection system (BIO-RAD).
For DEGs quantification expression detected by the TB Green ™ Premix Ex Taq ™ II (Tli RNaseH Plus) (Takara, Dalian, China), 10 μl reaction solution containing 1x Green TB Green TB Premix Ex Taq II, 10μM primer, one-third dilution of the cDNA sample, and then the reaction conditions were: 30s at 95°C followed by 40 cycles of 30 s at 95°C, and 30 s at 60°C. For DEMs quantification expression were detected by the Mir-X miRNA qRT-PCR SYBR kit, 25μl reaction solution containing 2X SYBR advantage premix, 50XROX dye, 10μM miRNA-specific Primer, 2μl cDNA sample, and then the reaction conditions were: 10s at 95°C followed by 40 cycles of 5s at 95°C, and 20s at 60°C. All mRNAs and miRNA expressions had three biological replicates with three technical replicates for each of biological replicate, ACTIN and U6 genes were used as reference genes, then relative expression levels were     Length distribution of assembled yam reads  Correlation analysis between the initiation and expansion stage  DEGs involved in IAA,ABA,GA and BR pathways of plant hormone signal transductions Figure 10 The numbers of up and down regulated transcription factors in expansion stage Figure 11 Heat map of the expression levels of highly expressed transcription factors Figure 12 Functional classification and pathway assignment of differentially expressed miRNAs by GO analysis Figure 13 Regulatory network from integrated analysis of miRNA-mRNA data Red stage represent nagatively correlated, and blue represent positively correlated in each miRNA-mRNA Verification of differential expressed genes by qRT-PCR of DEGs Figure 16 Verification of differential expressed miRNAs by qRT-PCR