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%) and 1879 (7.31%) genes were aligned to Elaeis guineensis, Phoenix dactylifera, Ananas comosus, Musa acuminata subsp and Malaccensis (Fig. 2c). Similar species distributions were also observed for yam tuber in previous research , 8229 (16.2%) genes in D. zingiberensis had the most hits from Elaeis guineensis, followed by Phoenix dactylifera (6857, 13.5%), Musa acuminate (2692, 5.3%).
After filtering low-quality reads and adaptor sequences, 6.57 Gb, 6.57 Gb, 6.58 Gb, 6.56 Gb, 6.57Gb, and 6.56Gb 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 1: Fig S1). Subsequently, 6388211 (25.11%), 5872589 (22.36%), 6086348 (22.98%) reads in tuber initiation stage and 4593044 (17.48%), 5032588 (18.58%), 4642869 (17.58%) reads 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 2: 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%) (Fig. 3c). 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.
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 2: 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 3: 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 2: 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 3: 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 2: 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 3: 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 2: 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 3: Table S2). These results indicated that many functional genes were involved in expansion stage of yam tuber.
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 4: 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 5: Table S4) and 68, 66, 70 total miRNAs were detected in tuber initiation stage (GH16_I), and 69, 69, 72 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 5: 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 6: 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 7 and 8: 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.