Integrated mRNA and miRNA transcriptome analysis reveals a regulatory network for tuber expansion in Chinese yam (Dioscorea opposite)

doi:10.21203/rs.2.9777/v3

Download PDF

Research article

Integrated mRNA and miRNA transcriptome analysis reveals a regulatory network for tuber expansion in Chinese yam (Dioscorea opposite)

https://doi.org/10.21203/rs.2.9777/v3

This work is licensed under a CC BY 4.0 License

Journal Publication

published 03 Feb, 2020

Read the published version in BMC Genomics →

You are reading this older preprint version

Read the latest preprint version →

Background Yam tuber is a storage organ, derived from the modified stem. Tuber expansion is a complex process, and depends on the expressions of genes that can be influenced by environmental and endogenous factors. However, little is known about the regulatory mechanism of tuber expansion. In order to identify genes and miRNAs involved in tuber expansion, we examined the mRNAs and small RNAs in Dioscorea opposite (Chinese yam) cv. Guihuai 16 tuber during its initiation and expansion stages. Results A total of 14238 differentially expressed genes in yam tuber at its expansion stage were identified by using RNA sequencing technology. Among them, 5723 genes were up-regulated, and 8515 genes were down-regulated. Functional analysis revealed the coordination of tuber plant involved in processes of cell events, metabolism, biosynthesis, and signal transduction pathways at transcriptional level, suggesting that these differentially expressed genes are somehow involved in response to tuber expansion, including CDPK, CaM, CDL, SAUR, DELLA, SuSy, and expansin. In addition, 541 transcription factor genes showed differential expression during the expansion stage at transcriptional level. MADS, bHLH, and GRAS were involved in cell differentiation, division, and expansion, which may relate to tuber expansion. Noteworthy, data analysis revealed that 22 known tuber miRNAs belong to 10 miRNA families, and 50 novel miRNAs were identified. The integrated analysis of miRNA-mRNA showed that 4 known miRNAs and 11 genes formed 14 miRNA-target mRNA pairs were co-expressed in expansion stage. miRNA160, miRNA396, miRNA535 and miRNA5021 may be involved in complex network to regulate cell division and differentiation in yam during its expansion stage. Conclusion The mRNA and miRNA datasets presented here identified a subset of candidate genes and miRNAs that are putatively associated with tuber expansion in yam, 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.

Epigenetics & Genomics

Yam tuber

Dioscorea opposite

Expansion

mRNA

Small RNA

Yams (Dioscorea spposite) are monocotyledonous plants belonging to the family Dioscoreaceae, and tuber is its harvested organ. Tuber originates from the expansion of underground stem, is suitable for nutrients storage, with many large parenchyma cells. Tuber morphogenesis and starch, along with accumulated proteins are two main processes of tuber growth [1]. The tuber morphogenesis of yam can be divided into three stages: initiation stage, expansion stage, and maturation stage. The expansion stage can 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 heredity, environment, and hormones [4]. Great efforts have been made to explore the physiological factors affecting the morphogenesis of yam tubers. The short-day treatment tended to promote tuber growth at the primary tuber growth stage of the plant, and bulbil development at the rapid tuber growth, but the responses varied among species and cultivar [5, 6]. Endogenous hormones including gibberellins (GA), acetic acid (IAA) and abscisic acid (ABA) performed 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]. Exogenous hormones have been used to study the mechanism of tuber expansion, GAs could promote tuber expansion and yield through in vitro and in vivo treatment [9, 10]. Exogenous GA application combined with ABA has promoted microtuber growth and expansion [11]. Exogenous JA was found to be essential for yam tuberization, and induced an increase in the number of tubers in vitro and in vivo [12, 13]. However, fundamental knowledges of endogenous metabolic networks are poor in tuber expansion.

The induction and growth of microtubers in vitro were controlled by nutrients, and sucrose concentration was the most crucial factor affecting tuberization and frequency of proliferation in yam [7, 14]. Yam tuber morphology was significantly correlated with nutrient accumulation and enzymatic activity. Sucrose, soluble sugars, and proteins increased significantly during tuber expansion stage, then subsequently decreased at maturity stage. Starch content increased throughout tuber morphogenesis, and sucrose synthase, sucrose phosphate synthase, and AGPase were significantly correlated with nutrient accumulation [15]. Although many DNA molecular markers have been used to uncover the genetic diversity and relationship among yam germplasms [16-18], little is known about specific genes involved in 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 the formation of adventitious roots [19]. PE2.1 and PE53 are the members of pectinesterase (PE) superfamily, which may be involved in the regulation of starch and sucrose metabolism and signaling pathways. Therefore, they may play an essential role in microtuber formation [20]. Tuber morphogenesis is a complex biological process involving many specific genes and proteins, especially at yam tuber expansion stage. Transcriptome techniques can efficiently find and detect these genes and proteins. Potato is a tuber crop. Several transcriptome analyses revealed that numerous genes are regulated in early stages of stolon-to-tuber transitions, or tuberization by nutrients, photoperiodic conditions, exogenous hormones, and stress in potato tuber [21-24]. Former transcriptomic study revealed that some putative genes were involved in dioscin biosynthesis [25], along with this 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 no reports of transcriptome study on tuber expansion.

microRNAs (miRNAs) are small, endogenous, non-coding RNAs that have essential functions in many biological processes, such as the regulations of growth and development, stress response, and metabolism. Many studies have shown that miRNAs play essential roles in root and tuber formation or development [27-29]. miR165/166 regulated root growth by determining the fate of root cells in Arabidopsis combined with phytohormone crosstalk, by negatively regulating its target genes auxin response factor ARF10, ARF16 and ARF17 [30]. miRNA172 and miR156 were involved in tuberization process, either as a component or a regulator of long-distance gibberellin signaling pathways [31, 32]. Potato specific miRNA193, miRNA152, and conserved miR172-1, miRNA172-5 showed significant expression during developmental stages of tuberization [28]. However many studies have found that miRNAs are involved in tuber and root development, the miRNA-mediated regulatory network during tuber expansion is still unclear.

Although whole-genome sequencing of the heterozygous diploid Guinea yam (D. rotundata) had been performed for sex determination [33], a detailed comparative mRNA and miRNA analysis during yam tuber expansion stage need to be detected. In this study, to identify and analyze the global gene and miRNA expression dataset in tuber expansion, six libraries prepared from D. opposite (Chinese yam) cv. Guihuai 16 tuber of initiation stage (GH16_I) and expansion stage (GH16_E) were sequenced by using a BGISEQ-500 platform. Furthermore, the association analysis between mRNA and miRNA expression was done, and the elucidation of the regulatory relationship of miRNA and their corresponding mRNA targets was studied for understanding the expansion of tuber.

Overview of RNA-Seq dynamics and small RNA sequencing

To identify the regulation of mRNA and miRNAs co-regulatory network during tuber expansion, the RNA-Seq and small RNA were examined during tuber initiation stage (GH16_I) and expansion stage (GH16_E) (Fig. 1). Meanwhile, transcriptome library was constructed from a pool of mixed RNA consisting of initiation and expansion stages in order to construct RNA-Seq and small RNA (named Total_1). Approximately 74.71 Mb original data in total were gained from BGISEQ-500 platform at BGI-Shenzhen (Table 1). After filtering low-quality reads and adaptor sequences, 6.67 Gb clean reads were obtained and processed by de novo analysis using Trinity software. The assembly produced a total of 54781 transcripts. Then, Tgicl software was used on transcripts to remove abundance, and 32207 genes were gained. The N50 statistic was 1508, which meant that more than 50% of the genes were longer than 1508 bp. The length distribution of all the assembled yam genes shown in Fig. 2a, which indicated that 7.91% of the complete transcripts and 13.00% of the total genes were longer than 2000bp. A total of 32207 genes were functionally annotated with 7 functional database (NR, NT, GO, KOG, KEGG, SwissProt, and InterPro), 25694 (79.78%), 16891 (52.45%), 17603 (54.66%), 19472 (60.46%), 20191 (62.69%), 22159 (68.80%), 8270 (25.68%) reads were annotated functionally respectively (Fig. 2b). 13566 genes were commonly annotated in NR, KOG, KEGG, SwissProt and InterPro databases. Based on the functional annotation results of NR database, the proportions of different species in the notes of genes were calculated, 8533 (33.21%), 5999 (23.35%), 1920 (7.47%), 1879 (7.31%) and 1545 (6.01%) genes were aligned to Elaeis guineensis, Phoenix dactylifera, Ananas comosus, Musa acuminata subsp. Malaccensis and Asparagus officinalis (Fig. 2c).

After filtering low-quality reads and adaptor sequences, 6.56 Gb, 6.57 Gb, 6.56 Gb, 6.57 Gb, 6.57Gb, and 6.58Gb clean reads were obtained in six RNA-Seq analysis libraries (initiation stage named GH16_I, expansion stage named GH16_E) respectively, and a total of 32026 expressed genes were detected (Table 1). The average mapping rate of transcriptome library (named Total_1) was 82.57%. A heat map cluster showed good correlations among replicates which indicated high repeatability of the data (Fig. 2d). In order to show the number of genes in different FPKM intervals of each mRNA libraries more intuitively, three situations of FPKM (FPKM<=1, FPKM 1-10, FPKM>=10) were counted the number of genes (Fig. 2e), indicating that most genes were expressed in the FPKM 1-10 ranges in the libraries. Genes with expression levels >5 FPKM were retained for statistical analysis.

Furthermore, the corresponding six small RNA libraries at the three time-points were also constructed for deep sequencing. Initially, a total of 170,957,171 reads were generated (Table 2). After filtering low-quality reads and adaptor sequences, 157,958,048 clean reads longer than 18nt for six libraries with an average of 26.32 M clean reads were obtained, and length distribution of clean reads showed that the classes of sRNA were 21-24nt (Additional file 2: Fig S1). Subsequently, 4593044 (17.48%), 5032588 (18.58%), 4642869 (17.58%) reads in tuber initiation stage and 6388211 (25.11%), 5872589 (22.36%), 6086348 (22.98%) in tuber expansion stage were mapped to sRNA database (rRNA, tRNA, snRNA and snoRNA), respectively.

Differentially expressed genes annotation by GO term and KEGG pathway

To identify differentially regulated genes in tuber expansion stage, DESeq software was used to compare the changes of gene expression between initiation and expansion stages. Among them, 5723 genes were up-regulated, 8515 genes were down-regulated, respectively, and it were differentially expressed in expansion stage (GH16_E), compared to initiation stage (GH16_I) (Additional file 1: Table S1).

For better comprehension of DEGs functions, 44 GO categories were identified. For biological processes, DEGs associated with cellular process (33%), metabolic process (31%), and biological regulation (9%) were enriched during expansion stage (Fig.3a). For cellular component, 10 GO categories were enriched in DEGs, including cell (24%), membrane (19%), membrane part (18%), and organelle (18%) (Fig. 3b). The molecular functions of the DEGs were mainly associated with catalytic activity (44%), binding (41%), transporter activity (5%), structural molecular activity (4%). Among the significant GO term analysis, 15 genes were enriched in cell wall polysaccharide metabolic process (GO:0010383), 15 genes were involved in hemicellulose metabolic process (GO:0010410), and 13 genes were related to xyloglucan metabolic process (GO:0010411) related to cell wall formation during expansion stage (Table 3). Besides, the results also revealed several significant expression genes involved in tissue development, root morphogenesis, root system development, and root development (Table 4).

KEGG is a signal pathway database with vibrant signal pathway map, 20 pathways were identified during yam tuber expansion stage. Interestingly, KEGG pathway analysis showed that plant hormone signal transduction (ko04075), biosynthesis of amino acids (ko01230) were enriched with DEGs during expansion stage (Fig. 3d). Other pathways such as MAPK signaling pathway (ko04016), starch and sucrose metabolism (ko00500), and carbon metabolism (ko01200) were also identified as involving 283, 204, and 236 DEGs, respectively. The metabolic pathways may be closely related to the development of tuber expansion and bioactive compound synthesis.

Comprehensive analysis of differentially expressed genes in expansion stage

Compared with initiation stage, there were a large number of DEGs in tuber expansion stage using NR, GO, and KEGG annotation. Signal transduction, cell wall, cell division, starch, and sucrose metabolism were selected for profiling during the expansion of yam tuber (Additional file 1: Table S2).

Hormone signal

A total of 242 DEGs were identified to be highly similar to many plant hormone signal pathways, including 131 down-regulated and 111 up-regulated DEGs in expansion stage (Additional file 1: Table S1). Interestingly, most plant hormone-related genes in GA, IAA, and ABA signal pathways were discovered during expansion stage.

In auxin transduction pathway, the transcriptional level of auxin influx carrier /auxin-responsive protein IAA (AUX/IAA) and small auxin up RNA (SAUR) were significantly down-regulated during the expansion stage, while auxin-responsive GH3 gene family (GH3) was up-regulated. In contrast, two auxin response factor ARFs (CL2135.Contig1_Total_1, and Unigene5660_Total_1) were shown high expression level during expansion stage, while other two ARFs (CL2887.Contig2_Total_1, Unigene5486_Total_1) were low expression level during expansion stage(Additional file 1: Table S2).

In gibberellin transduction pathway, the expression of gibberellin receptor GID2 was low expression during expansion stage. In contrast, DELLA proteins were highly expressed during the expansion stage. Meanwhile, protein phosphatase 2C（PP2C）was highly expressed during expansion stage.

MAPK and calcium signaling

Regulation of genes related to MAPK and calcium signaling during the expansion stage were also investigated. Six mitogen-activated protein kinases (MAPK) genes were up-regulated during expansion stage, while MPK6 and MPK8 were down-regulated. In summary, 48 DETs were homologous with calcium signal-related genes (Additional file 1: Table S1), including calcium-dependent protein kinases (CDPKs), calcium-binding proteins (CBPs), and calreticulin (CBL). It is worth noticing that CBLs were down-regulated during expansion stage (Additional file 1: Table S2).

Cell wall and cell cycle

A total of 98 transcripts homologous to the genes associated with cell wall and cell cycle were observed as differentially regulated during expansion stage (Additional file 1: Table S1), including xyloglucan endotransglucosylase/hydrolase (XTH), expansin, extension, cyclin-dependent kinases (CKS), cell division protease (ftsHs), cell division cycle 5-like protein (CDC5), cell division control protein (CDC), cyclin-dependent kinases (CDKs), and cyclin-dependent kinase inhibitor (CDKIs). All of the expansin, extension, cell wall synthesis, and CKS genes were down-regulated during expansion stage. Meanwhile, most of the cell cytoskeleton and XTH were down-regulated during expansion stage in yam (Additional file 1: Table S2).

Starch and sucrose metabolism

The major constituents of starch and sucrose metabolism during expansion stage are sucrose synthase genes (SuSy), sucrose phosphate synthase genes (SPS), starch synthase (SS), and invertase genes (INV) (Additional file 1: Table S1). Among them, SuSy were down-regulated during expansion stage. Interestingly, dioscorins, the major storage proteins in yam tubers, were significantly up-regulated during the expansion stage (Additional file 1: Table S2). These results indicated that many functional genes were involved in expansion stage of yam tuber.

Transcription factor

A total of 541 TF-encoding genes belonging to 48 TF families were differentially expressed during expansion stage, MYB, MYB-related, and AP2-EREP were enriched (Fig. 4). 286 TF encoding genes were up-regulated and 255 TF encoding genes were down-regulated, respectively (Additional file 1: Table S3). The most abundant TF gene families with the highest number of expressions during expansion stage were depicted by heat map (Fig. 5). Moreover, these genes were involved in circadian rhythm pathway, starch and sucrose metabolism pathway, and GA pathway by KEGG analysis respectively.

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 precursor to provide novel miRNA that can accurately be in the tubers. To compare miRNA expression in six libraries, the number of clean reads was used as background for normalization, and transcripts per million reads (TPM) was used to present the expression levels of miRNAs. Data analysis showed that there were 22 known miRNAs (21 and 20 in tuber initiation and expansion stage, respectively) and 50 novel miRNAs in yam tuber (Additional file 1: Table S4) and 69, 69, 72 total miRNAs were detected in tuber initiation stage (GH16_I), and 68, 66, 70 total miRNAs were detected in tuber expansion stage (GH16_E), respectively (Table 2). Distribution of normalized miRNAs expression showed that approximately 75%-81% of the total detected miRNA expression exceeded 10 TPM in six libraries (Additional file 1: Table S4).

Further analysis revealed that 22 known miRNAs belonging to 10 miRNA families, miRNA168, miRNA396, and novel miRNA160 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 during expansion stages. A total of 44 differentially expressed miRNAs were identified including 11 known and 33 novel miRNAs; it showed that 40 miRNAs were down-regulated, and 4 miRNAs were up-regulated during expansion stage (Table 5). Interestingly, miR160 family (miRNA160, miRNA160a-5p, miRNA160b-1, miRNA160h-1) and miR535a_1 were down-regulated during expansion stage, miRNA396b and miRNA168 were up-regulated compared to other miRNA families (Table 5).

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 to 34 DEGs (Additional file 1: Table S5). Notably, 4 known miRNAs and 11genes formed 14 miRNA-target mRNA pairs with co-expressed expression in expansion stage (Fig. 6). Furthermore, based on GO and KEGG analysis of the targets, it was showed that differentially expressed miRNAs were involved in multiple pathways. Some of target genes of miRNA were annotated as transcription factors, such as two auxin response factor ARFs (ARF17 and ARF18, targeted by miRNA160), growth-regulating factor 4 (GRF, targeted by miRNA396b), B3 domain-containing protein (ABI3VP1, targeted by miRNA396b), and squamosa promoter-binding-like protein 3 (SBP3, targeted by miR535a_1).

On the other hand, some target genes were revealed to be involved in transport and catabolism processes, such as autophagy-related protein 13b (ATG13, targeted by miRNA160), IST1-like protein (IST1, targeted by miR396a-5p and miR396b). Finally, some genes were annotated as other functional proteins, such as receptor expression-enhancing protein1 (REEP1, targeted by miR396a-5p and miR396b), photosystem II oxygen-evolving enhancer protein 2 (psbP, targeted by miR396a-5p and miR396b), and acetyl-CoA carboxylase biotin carboxyl carrier protein (accB, targeted by miR396b).

Validation of the DEGs and DEMs data using RT-qPCR

To identify the accuracy and reliability of RNA-Seq and miRNA data, RT-qPCR was used to measure the expressions of some DEGs and DEMs, including 15 mRNAs, and 7 miRNAs, specific primers were designed (Additional file 1: Tables S6 and S7). All mRNA and miRNA expressions have confirmed the accuracy of RNA-Seq and small RNA data. Overall, these results showed that 15 mRNAs showed similar expression patterns compared to DEGs analysis (Fig. 7a), and 7 miRNAs also showed very similar patterns compared to DEMs analysis (Fig. 7b). These results indicated that the RNA-Seq and small RNA data were reliable.

In tuber plants, tuber development is based on successive gene expression. Recently, the gene expression analysis of tuber development has been studied by transcriptome analysis, such as Solanum tuberosum, Miscanthus lutarioriparius, Raphanus sativus, Nelumbo nucifera and Aconitum heterophyllum[28, 34-37]. Yam is one of the most commercial stem tuber crops. Yam tuber expansion stage is an attractive theoretical model for studying the development of underground organs because it is the primary storage organ as nutrients reservoir. Therefore, the understanding of candidate genes involved in tuber expansion stage is supreme importance. In this study, 14238 genes were identified that have shown significant differential expression during expansion stage.

Multiple signaling pathways regulation

KEGG enrichment analysis revealed that a large number of genes were involved in pathways related to MAPK signaling pathway, plant hormone signal transduction, starch and sucrose metabolism, and carbon metabolism (Fig. 3d). These genes are involved in cell wall biosynthesis, cell proliferation, cell expansion, nutrient accumulation, primary metabolism, and hormone signal transduction [34, 38, 39]. Cellular processes are triggered by specific stimuli and hormones that are involved in a series of signaling pathways. Calcium transduction pathway is a primary essential signaling transduction category required for normal growth and development of plants. As a second messenger, it plays an important role in many fundamental cellular processes such as cytoplasmic streaming, thigmotropism, gravitropism, cell division, cell differentiation, photomorphogenesis, plant defense, and various stress responses [40-42]. A potato Ca²⁺-dependent protein kinase, StCDPK1, which was reported to be expressed during tuberization, possibly be involved in potato tuberization through transcriptional activation of certain genes [40, 43]. In addition, StCDPK1 was affected by JA, ABA, GA treatment [43-45], suggested that CDPK possibly be a mediator of hormone-to-tuber stimulation.

Moreover, majority of studies have demonstrated that mRNA abundance of CBL, CBP, and GaM were affected by cell division during swelling of potato tubers [40, 43, 46]. In this study, expression of calcium signaling-related genes, including CDPK and CBPs, increased during tuber expansion stage. Interesting, CBLs were down-regulated (Additional file 1: Table S2). The role of CBL in root development is characterized by a robust gene expression during cell proliferation in the root, and its expression in root is regulated by auxin and cytokinin [47, 48]. These results demonstrated that calcium signal is regulated during tuber expansion stage. Besides, MAPK signal pathways are known to play a central role in cell proliferation, differentiation and hormones, but it is not known to be involved in tuber or root formation [49]. In this study, MAPK signal-related genes were up-regulated during expansion stage, and it means that MAPK signal pathways take part in yam tuber expansion. Cell number and cell size are two key determinants of plant and organ size, while involved in cell division, cell expansion and cell cycle [50, 51]. In this study, XTH, expansin, extension, CKS, ftsHs, CDC5, CDC, CDKs, and CDKIs involved in cell division, cell expansion and cell cycles were differentially regulated during expansion stage. These differentially regulated genes involved in cell wall and cell cycle metabolism were identified in Arabidopsis, Rehmannia glutinosa, and radish [35, 36, 52].

Transcription factors regulation

Transcription factors have been identified to play an essential role in the regulation of plant growth, development, and secondary metabolism. Most transcription factors have been identified to play critical roles in organogenesis, including MADS, bHLH, and GRAS. In sweet potato, MADS were preferentially expressed in root at initiation of tuberization, and in the vascular cambium region where the active cell proliferation regulated by jasmonic acid, cytokinins and stress response [53, 54]. In Arabidopsis, bHLH family, such as PIF3, PIF4, and PIF5 are key regulators for cell elongation during seedling development, antagonistically regulated by light and gibberellins [55]. GRAS transcription factor family specifically regulated the initiation of axillary meristem in Arabidopsis, such as SCARECROW (SCR) is expressed in cortex/endodermal initial cells in the root system and played a crucial role in regulating the radial organization of the root [56]. DELLA, such as RGA and GAI, controls cell expansion and cell division in hypocotyl, shoot, root, and floral induction [57]. In this study, 12 MADS, 52 bHLH, and 17 GRAS were found to be expressed differently during expansion stage at transcript level, which plays key regulatory functions in cell differentiation, division and expansion (Fig. 4). Zinc finger CCCH domain-containing protein 14 (C3H) has shown to be involved in growth and development, hormone response, and response to biotic and biotic stresses [58]. All data suggest that these transcription factors may be potentially involved in the expansion stage.

Hormonal signaling regulation

A set of transcriptome evidence suggests that hormone signal pathway-related genes involved in regulating of storage organ formation have been identified in sweet potato, potato, radish, lotus, and carrot, while may be potentially involved in cell division, differentiation, and expansion [1, 35, 37, 59, 60]. In the KEGG pathway annotation, hormone signal pathway was the most enriching one that was involved in eight signal pathways (Fig. 3d). In radish, most Aux/IAA, ARFs, and SAUR genes were abundantly expressed in the root cortex splitting and expansion stage, implying that these transcripts may be involved in cambium cell expansion [35]. Among auxin signal-related genes, Aux/IAA, ARFs were highly expressed in expansion stage in this study, and auxin has been identified to regulate cell division and expansion by altering genes expression [1, 61], suggesting that Aux/IAA, ARFs may offer excellent candidates during expansion stage.

Interestingly, in this study, DELLA genes were up-regulated, while GIDI, like GID2, were down-regulated during expansion stage (Additional file 1: Table S2). DELLA proteins are nuclear transcriptional regulators that repress GA signaling and may restrict plant growth presumably by causing transcriptional reprogramming. Binding of GA to GID1 enhances the interaction between GID1 and DELLA, resulting in rapid degradation of DELLAs via the ubiquitin-proteasome pathway. The GA-GID1-DELLA signal control hypocotyl root elongation by reducing GA levels and DELLA interacts directly with PIF4, PIF3, and SCL to mediate crosstalk between GA and light signals [55, 62]. DELLA is significantly expressed in GA-mediated rhizomes, suggesting that DELLAs can regulate rhizomatic expansion [37]. Besides, GA, IAA, and ethylene affect cell growth in roots via DELLA proteins [57]. Overall, these results indicate that hormone signaling-related genes have a complex regulatory network involved in tuber expansion.

Starch and sucrose metabolism regulation

Starch and sucrose are considered to be one of the significant carbohydrates source, in which expansion tubers are highly coordinated with starch and sucrose metabolism genes in potato, radish, and lotus identified by transcriptome analysis [22, 35, 37]. In this study, alpha-amylase, beta-amylase, and isoamylase were up-regulated (Additional file 1: Table S1), which is similar to increase in beta-amylase activity in swollen taproot in radish [35], where starch content decreased during the tuber expansion stage [15], and some sucrose metabolism genes were detected during expansion stage, including SuSy, SPS, INV, and invertase inhibitor (Additional file 1: Table S1). Evidence shows that sucrose can be converted to starch in storage root by SuSy and AGPase, which means that they play a vital role in the early stage of radish expansion [35]. The SuSy gene was down-regulated during expansion stage. The results were similar to the taproot expansion stage in radish [35], implying that they may play a significant role in tuber expansion stage of yam. Yam tuber morphology is significantly correlated with sucrose and starch [15]. Therefore, these starch and sucrose metabolism genes are necessary for tuber expansion.

miRNAs regulation by targeting the potential genes

miRNA mediated gene regulation has been extensively studied in root and tuber development through transcriptional and post-transcriptional levels, which provides a better understanding of molecular regulatory network during tuber expansion stage. In this study, most miRNAs were down-regulated during expansion stage (Table 5). In maize leaves, miR160 was significantly up-regulated in meristem relative to the elongation and mature zones [63]. However, miR160 was down-regulated during yam tuber expansion stage (Table 5). In general, miR160 expression is different in different plants. In this study, some miRNA-mRNA pairs were observed during expansion stage (Fig. 6), miR160 was involved in root cell division and differentiation by regulating auxin response factors (ARFs) affecting root development in Arabidopsis thaliana [30]. miRNA160 was down-regulated during expansion stage, auxin response factor ARF18 and ARF17 were up-regulated during expansion stage (Fig. 6), whereas, overexpression of miR160 in transgenic rice not only down-regulated the expressions of ARF10, ARF16, and ARF17, but also inhibited root-cap cells differentiation, lost control of cell division and led to ectopic expansion of the apical stem cell populations [64]. Previous studies have shown that the miRNA160-targeted ARF has been identified to be involved in cell expansion and cell differentiation in radish root and potato tuber [28]. In addition, miRNA396 plays an essential role in root growth and inhibits leaf cell division by UV-B radiation [65, 66]. miR535 is expressed in fruit development [67]. miRNA160, miRNA396, and miR535 may be involved in complex a network of regulating cell division and differentiation during expansion stage.

Regulatory networks associated with tuber expansion

Tuber and root development and environmental responses involve gene regulatory networks [21,24,35,37]. miRNA and target genes have been extensively studied in tuber development [21, 28], which significantly advanced our understanding of molecular regulatory networks underlying during yam tuber expansion. In this study, based on our RNA-Seq and small RNA analysis, previous transcriptomics analysis and results of other tuber crops, a putative model of regulatory network associated with yam tuber expansion was proposed (Fig. 8). During vegetative growth, environmental factors are the first signal to stimulate tuber growth, and cells development through some signal transduction pathways (hormone, calcium and, MAPK signaling) and metabolism processes (cell wall, starch and sucrose metabolism). Some specific genes promote cell differentiation division and expansion at tuber expansion. In detail, it was found that CBL, CBP, and GaM may affect cell division and proliferation during potato tubers expansion [40, 43, 46]. Moreover, MAPK signal pathway is known to play a central role in cell proliferation, differentiation and hormone signaling [50, 51]. XTH, expansin, extension, CKS, ftsHs, CDC5, CDC, CDKs, and CDKIs involved in cell division, cell expansion and cell cycles were differentially regulated during expansion stage. A finding similar to these differences in cell wall and cell cycle regulatory genes are involved in cell wall metabolism in Arabidopsis, Rehmannia glutinosa, radish [35, 36, 52]. DELLA controls cell expansion and cell division in hypocotyl, shoot, root, and floral induction [57]. DELLAs were significantly expressed in GA-mediated rhizome development process, suggesting that these DELLAs can regulate rhizomatic expansion [37]. In sweet potato, MADS were preferentially expressed in root at the initiation of the tuber, and in the vascular cambium region where the active cell proliferation was regulated by jasmonic acid, cytokinins, and stress response [53, 54]. Most Aux/IAA, ARFs and SAUR genes were enriched and expressed in the root cortex during splitting and expanding stage in radish, suggesting that these transcripts may be involved in cell expansion in the cambium [35]. SuSy, SPS, INV, SS, and beta-amylase, were highly correlated with tuber expansion, such as potato, radish and lotus [22, 35, 37]. In this study, some miRNA and target genes were observed during yam tuber expansion stage, bHLH (target by miR5021), SBP (target by miR535a), and ARFs (targeted by miRNA160 and miR396), were identified to be involved in cell expansion and cell differentiation in radish root and potato tuber [28, 35]. Taken together, the result suggested that these DEGs may play a vital role in the regulatory network of tuber expansion.

A total of 14238 genes expressed differentially during yam tuber expansion stage were identified first time using RNA-Seq technology. These results reveal the coordination of tuber involved in the processes of cellular events, metabolism, biosynthesis, and signal transduction pathways at transcriptional levels. It is worth noticing that the integrated analysis of miRNA-mRNA was identified to be differentially expressed in tuber expansion. The mRNA and miRNA datasets presented here identify a subset of candidate genes and miRNAs putatively associated with tuber expansion in yam, and therefore propose a hypothetical model of the genetic regulatory network associated with tuber expansion in yam, which can provide the basis for molecular regulatory mechanism to study the expansion development in Dioscorea opposite species.

Plant Material

Dioscorea spposite cultivar Guihuai 16 was planted in April at the farm of Guangxi University in Nanning, Guangxi Province, China (22°53＇06.7〞N 108°21＇36.6〞E). Its healthy tubers germination, planting patterns and tuber growth stages were consistent as described by Gong [2] in 2016. All tuber growth stages were proved to have stable agronomic characteristics on three years of field observation. Tuber samples were collected at tuber initiation stage after field planting approximately forty days (initiation stage); additional samples were collected at sixty days (expansion stage). Anatomical characteristics of Guihuai 16 tuber at different developmental stages, including initiation and expansion stages were studied by Luo [68], and its histological characteristics consist of epidermis, parenchyma, and scattered vascular bundle. Five plants were selected randomly from every repetition each time. Then, the distal ends (5mm long) of five fresh tubers from a repetition were washed with distilled water, cut down into pieces, and mixed as a biological repetition. All samples were immediately frozen in liquid nitrogen, and stored at -80°C until further use for the construction of various libraries. Two tuber developmental stages, including initiation and expansion stage, were used for the construction of various libraries. Three types of libraries were constructed, including one transcriptome de novo analysis libraries for mixed sample of initiation and expansion stage (named Total_1), six RNA-Seq analysis libraries(initiation stage named GH16_I, expansion stage named GH16_E), and small RNA analysis(initiation stage named C, and expansion stage named T in original data, in order to be consistent with the RNA-Seq data analyze, we renamed GH16_I, GH16_E for initiation and expansion stage, respectively). Each sample had three biological replicates with two technical replicates for each biological replicate.

Construction and sequencing of transcriptome de novo, RNA-Seq and small RNA libraries

Information about the sample pool can be found on the subsections of plant materials above. Total RNA was isolated from distal ends of tuber samples using MiniBEST reagent (TaKaRa, Dalian, China), and RNA integrity was assessed by an Agilent 2100 BioAnalyzer. Construction of the transcriptome de novo and RNA-Seq libraries were performed at the Beijing Genome Institute (BGI; Shenzhen, China). Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads, fragmented, and reverse transcribed into cDNA, and the full-length cDNA was used T4 DNA ligase (Invitrogen, USA), and products were enriched by PCR to generate the cDNA library. Finally, the cDNA library was examined using an Agilent 2100 Bioanalyzer before sequencing on BGISEQ-500 system platform at BGI-Shenzhen. Paired-end reads were generated with a length of 90 bp for each cDNA libraries.

For the small RNA libraries, RNA bands of around 18-30nt (small RNAs) in length were isolated. Then small RNAs were ligated to a 5’-adaptor and a 3’-adaptor, and small RNAs were transcribed into cDNA. 14 rounds of PCR amplification were performed to enrich the cDNA fragments. After target fragments (100~120 bp) were purified, it was used for cluster generation and sequenced using BGISEQ-500 system platform at BGI-Shenzhen. BGISEQ-500 is a next-generation sequencing (NGS) technology platform developed by the Beijing Genomics Institute (BGI) in 2016 [69]. It is based on combinatorial Probe-Anchor Synthesis and improved DNA Nanoballs technology to generate short reads at a large scale. This platform has already been used for genome sequencing, RNA-Seq, and small RNA analysis.

Raw reads from mRNAs were subjected to remove the reads with adaptors or containing more than 5% unknown (‘N’) bases. Then, the low-quality reads (defined as reads having more than 20% of bases with quality≤15) were trimmed. The filtered clean reads were saved in FASTQ file, which has been submitted to Sequence Read Archive (SRA) of NCBI. Next, the Trinity software [70] was used to perform de novo assembly of the high-quality cleaned reads, and Tgicl software [71] was assembled to remove the redundant Trinity generated contigs. Final assembled datasets were performed to annotate genes function and expression. Small RNA reads were, also screened from raw sequencing reads by removing adaptors, poly-A sequences, and low-quality bases. Sequences shorter than 18nt or longer than 32nt were removed after trimming. Final filtered clean reads were saved in the FASTQ file, which was used for sequential analysis.

Processing and mapping of mRNA-Seq and small RNA libraries sequencing data

Functional databases NT, NR, GO, KOG, KEGG, SwissProt, and InterPro were used to annotate genes function, while Blastn, Blastx, Diamond, Blast2GO, and InterProScan5 were used to align genes. RNA-Seq reads were mapped to the transcriptome de novo analysis data using bowtie2 results with default parameters after pre-processing mRNA-Seq data [72]. Gene expression levels were presented as FPKM (fragments per kilobase of exon per million fragments mapped) values [73], and

The small RNA clean reads were submitted to the sRNA databases by using Bowtie, Cmsearch, and Rfam software [74]. The clean reads aligned to the sRNA databases with more than 90% matched length and a maximum alignment score of 90 were used from further analysis.

Identification of known and novel miRNAs, and predicting the targets of miRNAs

To identify the isoforms of known miRNA, cleaned sRNA reads were performed using the miRProf tool, and conserved miRNAs were identified using known plant miRNAs registered in miRBase [75]. To identify novel miRNAs and their precursors, unique reads were submitted to the RIPmiR [76]. Potential miRNA targets were predicted using the psRoot and Target Finder [77].

Identification and functional annotation of differentially expressed genes and miRNAs

The expression levels of mRNA genes were measured as FPKM, and genes with expression levels >5 FPKM were retained for statistical analysis. miRNA read counts were normalized to reads per million (TRM). Differentially expressed genes (DEGs) and miRNAs (DEMs) were identified by DEGseq [78]. Therefore, DEGs and DEMs were identified by stringent thresholds (|the FDR < 0.01 , |log2 (fold change)|≥1 and P-value < 0.005) . All DEGs were subjected to gene ontology and KEGG pathway analysis. In addition, enrichment GO and KEGG pathway analyses were implemented using p ≤ 0.05 and q-value < 0.9 as the threshold. Enrichment analysis of DEGs was analyzed using GO and KEGG databases to obtain a detailed description of the DEGs. Each DEG in tuber was predicted by aligning the gene sequences against the Plant Transcriptional Factor Database. The DEGs were classified according to their TF families. A miRNA-target gene regulatory network was constructed using Cytoscape_v3. 2. 1 program.

Validation of the DEGs and DEMs data using RT-qPCR

To identify the accuracy and reliability of mRNA and miRNA data, RT-qPCR was used to measure the expressions of DEGs and DEMs. Total RNA used for RNA-Seq and small RNA analysis previously was reversely transcribed into cDNA with PrimeScript™ RT reagent kit and Mir-X™ miRNA first-strand synthesis kit (Takara, Dalian, China), separately, according to the manufacture methods. RT-qPCR was performed using primers (Additional file 1: Table S6, S7) on real-time PCR detection system (BIO-RAD). Quantification expression for DEGs were detected by using the TB Green^™ Premix Ex Taq^™ II (Tli RNaseH Plus) (Takara, Dalian, China), 10 μl reaction solution containing 1x Green TB Premix Ex Taq II, 10μM primer.

One-third dilution of the cDNA sample was used, and the reaction conditions were: 30s at 95°C followed by 40 cycles of 30s at 95°C and 30s at 60°C. Quantification expression for DEMs were detected by the Mir-X miRNA qRT-PCR SYBR kit using 25μl reaction solution, containing 2× SYBR advantage premix, 50×ROX dye, 10μM miRNA-specific Primer, 2μl cDNA samples were used and 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 was used as reference gene, the relative expression levels were collected by the comparative Ct protocol.

DEG: Differentially expressed gene; DEM: Differentially expressed microRNA; GO: Gene ontology; KEGG: Kyoto encyclopedia of genes and genomes; miRNA: MicroRNA; qRT-PCR: Real-time quantitative RT-PCR; RNA-Seq: RNA sequencing; TF: Transcription factor; GA2ox: gibberellin 2-beta-dioxygenase; GA3ox: gibberellin 2-beta-dioxygenase; SnRK2: serine/threonine-protein kinase SRK2;PP2C: protein phosphatase 2C; B-ARR: myb family transcription factor; CRE1: farnesyl diphosphate synthase; BSK: BR-signaling kinase; BKI1: BRI1 kinase inhibitor 1;CTR1: serine/threonine-protein kinase CTR1; MYC: transcription factor MYC; JAZ: jasmonate ZIM domain-containing protein;TGA: transcription factor TGA;NPR1: regulatory protein NPR1.

Acknowledgment

Not applicable

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 81860670), the Guangxi Innovation Team Project of Tubers of Modern Agricultural Industrial Technology System of China (nycytxgxcxtd-11-01), and the Guangxi Natural Science Foundation of China (Grant No. 2018GXNSFDA281029).

Availability of data and materials

The materials of this study were provided by the College of Agriculture at Guangxi University. Correspondence and requests for materials should be addressed to Long-Fei He ([email protected]). The raw sequencing data have submitted to the NCBI SRA database (PRJNA533985). The assembled genes data also were submitted to GenBank TSA database(GHUN00000000).

Authors’ contributions

DX and LFH have designed experiments. YYZ, Saba, and SZL have analyzed the sequencing data for transcriptome assembly. JZ and AQW have developed qPCR experiments. YYZ and LFH have written the manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

We declare that the research was conducted in the absence of any commercial or financial relationships that could be a potential conflict of interest.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Author Details

^aCollege of Agriculture, Guangxi University, Nanning, 530004, P.R. China.

^bGuangxi Key Laboratory for Agro-Environment and Agro-Product Safety, Nanning, 530004, P.R. China.

^cGuangxi Colleges and Universities Key Laboratory of Crop Cultivation and Tillage, Nanning, 530004, P.R. China.

Aksenova NP, Konstantinova TN, Golyanovskaya SA, Sergeeva LI, Romanov GA. Hormonal regulation of tuber formation in potato plants. Russ J Plant Physl. 2012;59(4):451-466.
Gong MX, Luo HL, Wang AQ, Zhou YY, Huang WJ, Zhu PL, He LF. Phytohormone Profiling during tuber development of Chinese yam by ultra-high performance liquid chromatography–triple quadrupole tandem mass spectrometry. J Plant Growth Regul. 2016;36(2):362-373.
Matsumoto R, Kikuno H, Shiwachi H, Toyohara H, Takebayashi Y, Jikumaru Y, Kamiya Y. Growth of vine cuttings and fuctuations of concentrations of endogenous plant hormones in water yam (Dioscorea alata L.). Trop Agric Dev. 2013;57(1):23-30.
Yoshida Y, Takahashi H, Kanda H, Kanahamaz K. Effect of seed tuber weights on the development of tubers and flowering spikes in Japanese yams (Dioscorea japonica) grown under different photoperiods and with plant growth regulators. J Jpn Soc Hortic Sci. 2007;76(3):230-236.
Chen SW, Shiwachi H, Sanada A , Toyohara H. Theobroxide and day-length effects on the growth of yam (Dioscorea spp.). J.ISSAAS. 2010;16(1):22-30.
Shiwachi H, Ayankanmi T , Asiedu R. Effect of day length on the development of tubers in yams (Dioscorea spp.). Tropical Science. 2002;42:162-170.
Ovono PO, Kevers C, Dommes J. Axillary proliferation and tuberisation of Dioscorea cayenensis-D-rotundata complex. Plant Cell Tiss Org. 2007; 91(2):107-114.
Chen FQ, Fu Y, Wang DL, Gao X, Wang L. The effect of plant growth regulators and sucrose on the micropropagation and microtuberization of Dioscorea nipponica Makino. J Plant Growth Regul. 2007;26(1):38-45.
Gong MX, Luo HL, Yuan HJ, Wei SQ, Yang XD, He LF. Effects of exogenous gibberellin and paclobutrazol on tuber expansion and bulbil formation of Chinese yam. Acta Horticulturae Sinica. 2015;42(6):1175-1184.
Behera KK, Pani D, Sahoo S, Maharana T, Sethi BK. Effect of GA₃ and urea treatments on improvement of microtuber production and productivity of different types of planting material in greater yam (Dioscorea alata L.). Not Bot Horti Agrobot. 2009;37(2):81-84.
Hamadina EI, Craufurd PQ, Battey NH, Asiedu R. In vitro micro-tuber initiation and dormancy in yam. Ann Appl Biol. 2010;157(2):203-212.
Bazabakana R, Wattiez R, Baucher M, Diallo B, Jaziri M. Effect of jasmonic acid on developmental morphology during in vitro tuberization of Dioscorea alata (L.). Plant Growth Regul. 2003;40(3):229-237.
Kim SK, Kim JT, Jang SW, Lee SC, Lee BH, Lee IJ. Exogenous effect of gibberellins and jasmonate on tuber enlargement of Dioscorea opposita. Agron Res. 2005; 3(1):39-44.
Yan H, Yang L, Li Y. Axillary shoot proliferation and tuberization of Dioscorea fordii Prain et Burk. Plant Cell Tiss Org. 2011;104(2):193-198.
Liang RF, Li CZ, Zhang J, He LF, Wei BH, Gan XQ, He HY. Changes of matter accumulation and relative enzymatic activity during yam tuber development. Acta Agronomica Sinica. 2011;37(5):903-910.
Loko YL, Bhattacharjee R, Agre AP, Dossou-Aminon I, Orobiyi A, Djedatin GL, Dansi A. Genetic diversity and relationship of Guinea yam (Dioscorea cayenensis Lam.–D.rotundata Poir.complex) germplasm in Benin (West Africa) using microsatellite markers. Genet Resour Crop Ev. 2017; 64(6):1205-1219.
Girma G, Gedil M, Spillane C. Morphological, SSR and ploidy analysis of water yam (Dioscorea alata L.) accessions for utilization of aerial tubers as planting materials. Genet Resour Crop Ev. 2017;64(2):291-305.
Tamiru M, Yamanaka S, Mitsuoka C, Babil P, Takagi H, Lopez-Montes A, Sartie A, Asiedu R, Terauchi R. Development of genomic simple sequence repeat markers for yam. Crop Sci. 2015; 55(5):2191-2200.
Akakpo R, Scarcelli N, Chair H, Dansi A, Djedatin G, Thuillet AC, Rhone B, Francois O, Alix K, Vigouroux Y. Molecular basis of African yam domestication: analyses of selection point to root development, starch biosynthesis, and photosynthesis related genes. BMC Genomics.2017; 18(1):782.
Zhao XT, Zhang XL, Guo XB, Li SJ, Han LL, Song ZH, Wang YY, Li JH, Li MJ. Identification and validation of reference genes for qRT-PCR studies of gene expression in Dioscorea opposita. Biomed Res Int. 2016;3089584.
Kondhare KR, Malankar NN, Devani RS, Banerjee AK. Genome-wide transcriptome analysis reveals small RNA profiles involved in early stages of stolon-to-tuber transitions in potato under photoperiodic conditions. BMC Plant Biology. 2018;18(1):284.
Van Harsselaar JK, Lorenz J, Senning M, Sonnewald U, Sonnewald S. Genome-wide analysis of starch metabolism genes in potato (Solanum tuberosum L.). BMC Genomics. 2017;18(1):37.
Gong L, Zhang H, Gan X, Zhang L, Chen Y, Nie F, Shi L, Li M, Guo Z, Zhang G et al. Transcriptome profiling of the potato (Solanum tuberosum L.) plant under drought stress and water-stimulus conditions. Plos One. 2015; 10(5):e0128041.
Shan J, Song W, Zhou J, Wang X, Xie C, Gao X, Xie T, Liu J. Transcriptome analysis reveals novel genes potentially involved in photoperiodic tuberization in potato. Genomics. 2013; 102(4):388-396.
Li J, Liang Q, Li CF, Liu MD, Zhang YS. Comparative transcriptome analysis identifies putative genes involved in dioscin biosynthesis in Dioscorea zingiberensis. Molecules. 2018;23(2):454
Wu ZG, Jiang W, Mantri N, Bao XQ, Chen SL, Tao ZM. Transciptome analysis reveals flavonoid biosynthesis regulation and simple sequence repeats in yam (Dioscorea alata L.) tubers. BMC Genomics. 2015;16:346.
Liu C, Liu XX, Xu WL, Fu WM, Wang FD, Gao JW, Li QY, Zhang ZG, Li JJ, Wang SF. Identification of miRNAs and their targets in regulating tuberous root development in radish using small RNA and degradome analyses. 3 Biotech. 2018;8(7):311.
Lakhotia N, Joshi G, Bhardwaj AR, Katiyar-Agarwal S, Agarwal M, Jagannath A, Goel S, Kumar A. Identification and characterization of miRNAome in root, stem, leaf and tuber developmental stages of potato (Solanum tuberosum L.) by high-throughput sequencing. BMC Plant Biology. 2014; 14:6.
Sun YY, Qiu Y, Zhang XH, Chen XH, Shen D, Wang HP, Li XX. Genome-wide identification of microRNAs associated with taproot development in radish (Raphanus sativus L.). Gene. 2015;569(1):118-126.
Mallory AC, Bartel DP, Bartel B. MicroRNA-directed regulation of Arabidopsis AUXIN RESPONSE FACTOR17 is essential for proper development and modulates expression of early auxin response genes. Plant Cell. 2005;17(5):1360-1375.
Sarkar D. Photoperiodic inhibition of potato tuberization: an update. Plant Growth Regul. 2010;62(2):117-125.
Martin A, Adam H, Diaz-Mendoza M, Zurczak M, Gonzalez-Schain ND, Suarez-Lopez P. Graft-transmissible induction of potato tuberization by the microRNA miR172. Development. 2009;136(17):2873-2881.
Tamiru M, Natsume S, Takagi H, White B, Yaegashi H, Shimizu M, Yoshida K, Uemura A, Oikawa K, Abe A. Genome sequencing of the staple food crop white Guinea yam enables the development of a molecular marker for sex determination. BMC Biology. 2017;15(1):86.
Hu RB, Yu CJ, Wang XY, Jia CL, Pei SQ, He K, He G, Kong YZ, Zhou GK. De novo transcriptome analysis of miscanthus iutarioriparius identifies candidate genes in rhizome development. Front Plant Sci. 2017;8:492.
Yu RG, Wang J, Xu L, Wang Y, Wang RH, Zhu XW, Sun XC, Luo XB, Xie Y, Everlyne M, et al: Transcriptome profiling of taproot reveals complex regulatory networks during taproot thickening in radish (Raphanus sativus L.). Front Plant Sci. 2016;7:1210.
Malhotra N, Sood H, Chauhan RS. Transcriptome-wide mining suggests conglomerate of genes associated with tuberous root growth and development in Aconitum heterophyllum Wall. 3 Biotech. 2016;6(2):152.
Yang M, Zhu LP, Pan C, Xu LM, Liu YL, Ke WD, Yang PF. Transcriptomic analysis of the regulation of rhizome formation in temperate and tropical lotus (Nelumbo nucifera). Sci Rep. 2015;5:13059.
Ursache R, Nieminen K, Helariutta Y. Genetic and hormonal regulation of cambial development. Physiologia plantarum. 2013;147(1):36-45.
Bhalerao RP, Fischer U. Environmental and hormonal control of cambial stem cell dynamics. J Exp Bot. 2017;68(1):79-87.
Nookaraju A, Pandey SK, Upadhyaya CP, Heung JJ, Kim HS, Chun SC, Kim DH, Park SW. Role of Ca2⁺-mediated signaling in potato tuberization: An overview. Bot Stud. 2012;53(2):177-189.
Chen T, Wu X, Chen Y, Li X, Huang M, Zheng M, Baluska F, Samaj J, Lin J. Combined proteomic and cytological analysis of Ca2⁺-calmodulin regulation in Picea meyeri pollen tube growth. Plant Physiology. 2009;149(2):1111-1126.
Feng JL, Yang ZJ, Chen SP, El-Kassaby YA, Chen H. Signaling pathway in development of Camellia oleifera nurse seedling grafting union. Trees-Struct Funct. 2017;31(5):1543-1558.
Gargantini PR, Giammaria V, Grandellis C, Feingold SE, Maldonado S, Ulloa RM. Genomic and functional characterization of StCDPK1. Plant Mol Biol. 2009;70(1-2):153-172.
Yu XC, Li MJ, Gao GF, Feng HZ, Geng XQ, Peng CC, Zhu SY, Wang XJ, Shen YY, Zhang DP. Abscisic acid stimulates a calcium-dependent protein kinase in grape berry. Plant physiology. 2006;140(2):558-579.
Ulloa RM, Raices M, MacIntosh GC, Maldonado S, Tellez-Inon MT. Jasmonic acid affects plant morphology and calcium-dependent protein kinase expression and activity in Solanum tuberosum. Physiologia Plantarum. 2002;115(3):417-427.
Kim MS, Kim HS, Kim YS, Baek KH, Oh HW, Hahn KW, Bae RN, Lee IJ, Joung H, Jeon JH. Superoxide anion regulates plant growth and tuber development of potato. Plant Cell Rep. 2007;26(10):1717-1725.
Kumar Meena M, Kumar Vishwakarma N, Tripathi V, Chattopadhyay D. CBL-interacting protein kinase 25 contributes to root meristem development. J Exp Bot. 2019;70(1):133-147.
Cui XY, Du YT, Fu JD, Yu TF, Wang CT, Chen M, Chen J, Ma YZ, Xu ZS. Wheat CBL-interacting protein kinase 23 positively regulates drought stress and ABA responses. BMC Plant Biology. 2018;18:93.
Hirt H. Connecting oxidative stress, auxin, and cell cycle regulation through a plant mitogen-activated protein kinase pathway. P Natl Acad Sci USA. 2000;97(6):2405-2407.
Polyn S, Willems A, De Veylder L. Cell cycle entry, maintenance, and exit during plant development. Current opinion in plant biology. 2015;23:1-7.
Guo M, Simmons CR. Cell number counts--the fw2.2 and CNR genes and implications for controlling plant fruit and organ size. Plant Sci. 2011;181(1):1-7.
Paque S, Mouille G, Grandont L, Alabadi D, Gaertner C, Goyallon A, Muller P, Primard-Brisset C, Sormani R, Blazquez MA et al. AUXIN BINDING PROTEIN1 links cell wall remodeling, auxin signaling, and cell expansion in Arabidopsis. Plant Cell. 2014;26(1):280-295.
Dong TT, Song WH, Tan CT, Zhou ZL, Yu JW, Han RP, Zhu MK, Li ZY. Molecular characterization of nine sweet potato (Ipomoea batatas Lam.) MADS-box transcription factors during storage root development and following abiotic stress. Plant Breeding. 2018;137(5):790-804.
Song WH, Tan CT, Zhu MK, Xu T, Dong TK, Li ZY. Cloning and expression Analysis of Stress Treatments of a MADS-box Transcription Factor in Sweet Potato. Bulletin of Botanical Research. 2017;37(4):587-595,627.
de Lucas M, Daviere JM, Rodriguez Falcon M, Pontin M, Iglesias Pedraz JM, Lorrain S, Fankhauser C, Blazquez MA, Titarenko E, Prat S. A molecular framework for light and gibberellin control of cell elongation. Nature. 2008; 451(7177):480-484.
Heo J, Chang KS, Kim IA, Lee M, Lee SA, Song S, Lee MM, Lim J. Funneling of gibberellin signaling by the GRAS transcription regulator SCARECROW-LIKE 3 in the Arabidopsis root. P Natl Acad Sci USA. 2011; 108(5):2166-2171.
Daviere JM, Achard P. A Pivotal Role of DELLAs in Regulating multiple hormone signals. Mol Plant. 2016;9(1):10-20.
Li WT, He M, Wang J, Wang YP. Zinc finger protein (ZFP) in plants-A review. Plant Omics. 2013;6(6):474-480.
Que F, Wang GL, Li T, Wang YH, Xu ZS, Xiong AS. Genome-wide identification, expansion, and evolution analysis of homeobox genes and their expression profiles during root development in carrot. Functional & Integrative Genomics. 2018;18(6):685-700.
Ponniah SK, Thimmapuram J, Bhide K, Kalavacharla V, Manoharan M. Comparative analysis of the root transcriptomes of cultivated sweetpotato (Ipomoea batatas L. Lam) and its wild ancestor (Ipomoea trifida Kunth G. Don). BMC Plant Biology. 2017;17:9.
Ursache R, Nieminen K, Helariutta Y. Genetic and hormonal regulation of cambial development. Physiol Plant. 2013;147(1):36-45.
Zhang Z, Ogawa M, Fleet CM, Zentella R, Hu J, Heo JO, Lim J, Kamiya Y, Yamaguchi S, Sun T. SCARECROW-LIKE 3 promotes gibberellin signaling by antagonizing master growth repressor DELLA in Arabidopsis. P Natl Acad Sci USA. 2011;108(5):2160-2165.
Aydinoglu F, Lucas SJ. Identification and expression profiles of putative leaf growth related microRNAs in maize (Zea mays L.) hybrid ADA313. Gene. 2019;690:57-67.
Yang SN, Wang J, Chen ZX, Tian XJ, Zhang JW, Long Y, Pei XW, Yuan QH. Clone and function analysis miR160f in common wild rice(Oryza rufipogon Griff.). Biotechnology Bulletin. 2014;11:114-118.
Bao ML, Bian HW, Zha YL, Li FY, Sun YZ, Bai B, Chen ZH, Wang JH, Zhu MY, Han N. miR396a-mediated masic helix-loop-helix transcription factor bHLH74 repression acts as a regulator for root growth in Arabidopsis seedlings. Plant and Cell Physiology. 2014;55(7):1343-1353.
Casadevall R, Rodriguez RE, Debernardi JM, Palatnik JF, Casati P. Repression of growth regulating factors by the microRNA396 inhibits cell proliferation by UV-B radiation in Arabidopsis leaves. Plant Cell. 2013;25(9):3570-3583.
Shi M, Hu X, Wei Y, Hou X, Yuan X, Liu J, Liu Y. Genome-wide profiling of small RNAs and degradome revealed conserved regulations of miRNAs on auxin-responsive genes during fruit enlargement in peaches. Int J Mol Sci. 2017;18(12).
Luo HL, Gong MX, Zhou YY, He LF: Analysis of starch accumulation and differentially expressed proteins during the development of Chinese yam tuber. Journal of South China Agricultural University. 2018;39(6):61-69.
Goodwin S, McPherson JD, McCombie WR: Coming of age: ten years of next-generation sequencing technologies. Nat Rev Genet. 2016;17(6):333-351.
70. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011; 29(7):644.
71. Pertea G, Huang X, Liang F, Antonescu V, Sultana R, Karamycheva S, et al.TIGR Gene Indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformatics. 2003;19(5):651.
72. Cock PJ, Fields CJ, Goto N, Heuer ML, Rice PM. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids R 2010;38(6):1767-1771.
73. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature M 2012; 9(4):357-359.
74. Nawrocki EP, Eddy SR. Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics. 2013;29(22):2933-2935.
75. Nawrocki EP, Burge SW, Bateman A, Daub J, Eberhardt RY, Eddy SR, Floden EW, Gardner PP, Jones TA, Tate J et al. Rfam 12.0: updates to the RNA families database. Nucleic acids research. 2015;43(Database issue):D130-137.
76. Breakfield NW, Corcoran DL, Petricka JJ, Shen J, Sae-Seaw J, Rubio-Somoza I, Weigel D, Ohler U, Benfey PN. High-resolution experimental and computational profiling of tissue-specific known and novel miRNAs in Arabidopsis. Genome research. 2012;22(1):163-176.
77. Wu HJ, Ma YK, Chen T, Wang M, Wang XJ. PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids R 2012;40(Web Server issue):W22-28.
78. Wang L, Feng Z, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. 2010;26(1):136-138.

Table 1 Statistic analysis of clean reads for mRNA in tuber initiation and expansion stages in yam

Sample	Total_1	GH16_E_r1	GH16_E_r2	GH16_E_r3	GH16_I_r1	GH16_I_r2	GH16_I_r3	Sum
Total Raw Reads(Mb)	74.71	66.43	66.43	66.35	66.43	66.42	66.52	－－
Total Clean Bases(Gb)	6.67	6.56	6.57	6.56	6.57	6.57	6.58	－－
Clean Reads Q20(%)	96.4	98.52	98.5	98.57	98.63	98.56	98.54	－－
Clean Reads Q30(%)	88.22	92.67	92.63	92.95	93.12	92.86	92.78	－－
Clean Reads Ratio(%)	89.34	98.7	98.87	98.9	98.91	98.92	98.96	－－
Total Mapped Reads(%)	－－	82.05	81.77	82.03	83.07	83.34	83.17	－－
Total Expressed Genes	－－	30012	29893	29811	29658	29905	29839	32026
Total trinity Transcripts	54781
Total Tgicl Genes	32207
GeneGC(%)	44.63
GeneN50	1508

Table 2 Statistic analysis of clean reads for small RNA sequencing in tuber initiation and expansion stages in yam

Sample name	Tota reads	Clean reads	Mapped reads	Known miRNA	Novel miRNA	Total miRNA
GH16_E_r1	28256295	26276062	4593044	20	49	69
GH16_E_r2	29343139	27085820	5032588	20	49	69
GH16_E_r3	28395798	26406723	4642869	22	50	72
GH16_I_r1	27687839	25438600	6388211	18	50	68
GH16_I_r2	28663585	26268556	5872589	19	47	66
GH16_I_r3	28610515	26482287	6086348	21	49	70
Sum	170957171	157958048

Table 3 GO enrichment analysis of DEGs
Type	ID	Term	Gene Number	Rich Ratio	P value
Molecular function	GO:0005198	Structural molecule activity	183	0.62	1.83867E-07
	GO:0003735	Structural constituent of ribosome	143	0.61	1.71036E-05
	GO:0003700	DNA binding transcription factor activity	116	0.59	0.001099451
Cellular component	GO:0005576	Extracellular region	91	0.70	2.22612E-07
	GO:0043228	Non-membrane-bounded organelle	285	0.58	1.68443E-06
	GO:0043232	Intracellular non-membrane-bounded organelle	285	0.58	1.68443E-06
	GO:0005840	Ribosome	157	0.61	1.90586E-05
	GO:0005618	Cell wall	49	0.72	4.36888E-05
	GO:0030312	External encapsulating structure	49	0.72	4.36888E-05
	GO:0048046	Apoplast	32	0.78	0.000075892
	GO:0022625	Cytosolic large ribosomal subunit	22	0.81	0.000355453
Biological process	GO:0044262	Cellular carbohydrate metabolic process	77	0.69	1.68372E-05
	GO:0005975	Carbohydrate metabolic process	171	0.60	0.00025216
	GO:0044264	Cellular polysaccharide metabolic process	53	0.70	0.000285538
	GO:0045786	Negative regulation of cell cycle	11	1.00	0.000452873
	GO:0006073	Cellular glucan metabolic process	50	0.69	0.000495823
	GO:0010383	Cell wall polysaccharide metabolic process	15	0.88	0.001069447
	GO:0010410	Hemicellulose metabolic process	15	0.88	0.001069447
	GO:0010411	Xyloglucan metabolic process	13	0.93	0.000841262

Table 4 Functional classification and pathway assignment of differentially expressed DEG by GO in expansion stage
Tissue development [GO:0009888]	Gene ID	log2(GH16_E/GH16_I)
Actin-related protein 3	CL1179.Contig2_Total_1	-2.00
Anaphase-promoting complex subunit 10	CL1997.Contig4_Total_1	1.82
Alpha-tubulin	CL28.Contig3_Total_1	-5.08
Tubulin alpha chain	CL3054.Contig1_Total_1	-1.81
Glutamine synthetase nodule isozyme	Unigene128_Total_1	-4.98
Phosphoenolpyruvate carboxylase 2	Unigene13453_Total_1	-1.23
Phosphoenolpyruvate carboxylase 2	Unigene13455_Total_1	1.32
Anaphase-promoting complex subunit 10	Unigene18039_Total_1	-1.15
Anaphase-promoting complex subunit 6	Unigene2600_Total_1	-1.06
Actin-related protein 2	Unigene4771_Total_1	-1.10
ATPase ASNA1 homolog	Unigene5972_Total_1	-1.32
Homeobox protein knotted-1-like 3	Unigene6778_Total_1	3.48
Root morphogenesi, root system development, root development [GO:0009888,GO:0022622,GO:0048364]
Cytoplasmic tRNA 2-thiolation protein 2	CL1237.Contig1_Total_1	-1.02
Cytoplasmic tRNA 2-thiolation protein 2	CL1237.Contig2_Total_1	-2.28
Mediator of RNA polymerase II transcription subunit 32	CL2787.Contig2_Total_1	2.02
Succinate dehydrogenase assembly factor 2	CL3034.Contig1_Total_1	-1.04
Guanine nucleotide-binding protein	Unigene18752_Total_1	-3.18
Enhanced ethylene response protein 5	Unigene2193_Total_1	-1.75
ATPase ASNA1	Unigene5972_Total_1	-1.32

Table 5 Differentially expressed miRNAs involved in expansion stage
miRNA id	Read count(GH16_I)	Read count(GH16_E)	Expression(GH16_I)	Expression(GH16_E)	log2Ratio (GH16_E/GH16_I)	P value
novel_mir43	54974	33041	3682.20	1731.72	-1.97	0.00000
novel_mir35	52437	32873	3515.56	1728.51	-1.91	0.00000
miR160a-5p	18032	121	1171.81	6.42	-8.45	0.00000
miR168a-5p	55653	59949	3708.90	3153.82	-1.13	0.00000
miR168	20208	106873	1334.10	5618.38	1.17	0.00000
novel_mir1	16275	7760	1085.59	407.58	-2.30	0.00000
miR396b	1097	21915	72.84	1176.08	3.09	0.00000
miR535a_1	24978	25221	1654.31	1315.07	-1.22	0.00000
miR396a-5p	12295	7317	816.34	374.44	-1.98	0.00000
novel_mir33	6408	4269	425.49	224.55	-1.82	0.00000
novel_mir15	9666	10433	642.81	549.12	-1.12	0.00000
novel_mir11	4852	3769	323.07	199.13	-1.60	0.00000
novel_mir25	5833	6318	386.43	333.34	-1.12	0.00000
novel_mir32	3196	2381	212.72	125.72	-1.66	0.00000
novel_mir5	2039	1094	135.31	58.03	-2.13	0.00000
novel_mir50	2777	2281	183.85	120.27	-1.52	0.00000
novel_mir10	3197	3009	213.02	158.49	-1.32	0.00000
novel_mir44	3381	3299	224.65	173.72	-1.27	0.00000
novel_mir36	1859	1246	124.13	65.29	-1.81	0.00000
novel_mir31	1632	1049	108.56	55.06	-1.87	0.00000
novel_mir18	1701	1305	114.05	68.55	-1.62	0.00000
novel_mir19	1491	1071	99.15	56.35	-1.71	0.00000
novel_mir12	1666	1308	111.28	69.05	-1.58	0.00000
novel_mir45	1478	1093	99.12	57.63	-1.67	0.00000
miR160	926	451	62.61	24.59	-2.27	0.00000
novel_mir48	1107	690	74.52	36.21	-1.92	0.00000
novel_mir26	1488	1201	99.79	63.01	-1.54	0.00000
miR160b_1	463	30	30.27	1.62	-5.18	0.00000
novel_mir46	1268	926	83.82	48.75	-1.69	0.00000
novel_mir42	694	275	46.22	14.50	-2.57	0.00000
novel_mir49	676	265	45.07	13.82	-2.59	0.00000
novel_mir24	1045	779	69.79	40.94	-1.66	0.00000
novel_mir23	377	109	25.16	5.68	-3.03	0.00000
miR168a-3p	1330	1489	87.83	78.66	-1.07	0.00000
novel_mir20	374	2206	24.76	116.77	1.33	0.00000
novel_mir30	964	999	63.39	52.50	-1.18	0.00000
novel_mir37	173	45	11.52	2.37	-3.18	0.00000
novel_mir47	90	10	6.03	0.50	-4.40	0.00000
novel_mir4	394	447	25.80	23.70	-1.05	0.00000
novel_mir34	10	200	0.68	10.70	3.09	0.00000
novel_mir8	197	207	13.12	10.72	-1.16	0.00000
miR168_1	30	5	1.98	0.26	-3.82	0.00000
novel_mir9	123	142	8.24	7.48	-1.03	0.00000
miR160h_1	7	0	0.46	0.00	-5.04	0.00041

Additional file 1:

Table S1. Differentially expressed genes (DEGs) in expansion stage compared with initiation stage. These genes were manually categorized according to their corresponding NR annotation, GO annotation, and KEGG annotation.

Table S2. The differentially expressed genes involved in signal transduction pathways in expansion stage (FPKM>5).

Table S3. The differentially expressed transcription factors involved in expansion stage.

Table S4. Analysis of clean reads for small RNA sequencing in initiation and expansion stages.

Table S5. Identification of miRNAs target genes in expansion stage.

Table S6. The primers of 15 DEGs validated by RT-qPCR analysis.

Table S7. The primers of 7 DEMs validated by RT-qPCR analysis.

Additional file 2:

Fig S1: Length distribution of small RNA sequences in small RNA libraries.

Download PDF

Journal Publication

published 03 Feb, 2020

Read the published version in BMC Genomics →

Editorial decision: Minor revision
22 Oct, 2019
Review #2 received at journal
16 Oct, 2019
Review #1 received at journal
16 Oct, 2019
Reviewer #2 agreed at journal
17 Sep, 2019
Reviewer #1 agreed at journal
13 Sep, 2019
Reviewers invited by journal
12 Sep, 2019
Editor assigned by journal
04 Sep, 2019
Submission checks completed at journal
03 Sep, 2019
Editor invited by journal
03 Sep, 2019

You are reading this older preprint version

Read the latest preprint version →

Integrated mRNA and miRNA transcriptome analysis reveals a regulatory network for tuber expansion in Chinese yam (Dioscorea opposite)

Status:

Journal Publication

Version 3

Abstract

Figures

Background

Results

Discussion

Conclusions

Materials and Methods

Abbreviations

Declarations

References

Tables

Additional Files

Supplementary Files

Status:

Journal Publication

Version 3