The mungbean variety “Zhonglu 1” was used in this study. Mungbean seeds were washed three times with distilled water, and then placed in 12-cm Petri dishes and imbibed in distilled water at 20°C in a dark biochemical incubator. Seeds from each imbibition time were dried silica gel at different dehydration times (5–50 h). Fifty seeds were placed in one plate and three biological replicates were prepared for each time-point sample. Then all the seeds were used for further seed germination and the germination rates and water contents were detected during the dehydration. For RNA-seq, after 3, 6, 18, and 24 h of imbibition, the mungbean seeds were dried for 24 h and collected for sequencing with three replicates for each time-point sample. In order to compare the differentially expressed genes (DEGs) during the imbibition, a sample of each stage without dehydration was used as a control (“ck”), while the seeds with dehydration treatment were used as experiment groups (“sy”). The seeds were collected and immediately frozen in liquid nitrogen, then stored at −80°C until further use.

Total RNA was extracted from the seeds using an RNeasy Plant Mini Kit (QIAGEN) following the manufacturer's protocol. The extracted RNA was treated with DNase I (New England BioLabs) to remove the contaminated DNA. The concentration and quality of each sample was detected by Nanodrop2000. Using Oligo (dT), mRNAs were enriched from total RNA and were randomly cleaved into short fragments. The cDNAs were synthesized using random hexamer primers, RNase H and DNA polymerase I. The double-stand cDNAs were purified and ligated to adaptors for paired-end sequencing. The quality and quantity of the libraries were detected using an Agilent 2100 Bioanalyzer and ABI real time RT-PCR equipment. The qualified cDNA libraries were sequenced by Biomarker Technologies (Beijing, China) using an Illumina HiSeq 2500 platform with PE125.

Raw reads were checked with FastQC package (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/), and adaptor sequences and low quality reads were removed. A stringent filtering criterion was used to minimize the effects of sequencing errors during assembly. The reads 70% of the bases in a read having high phred quality scores (≥20) were used for assembly and reads lengths shorter than 50 bp were discarded. The obtained clean reads of all twenty four samples were assembled by Trinity v2.1.0 (Released on Sep 29, 2015) with a paired-end model (Grabherr et al., 2011). To annotate the assembled unigenes, sequences were aligned by BLASTx (E-value <1e⁻⁵) to protein databases, including the non-redundant protein (NR) database, Clusters of eukaryotic Orthologous Groups of proteins (COG) database, Swiss-Prot and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. The transcript abundance was normalized by the Fragments Per Kilobase of transcript per Million mapped reads (FPKM) value.

The clean reads of each sample were mapped back to assembled contigs using bowtie2 (Langmead and Salzberg, 2012). The assembled contigs with more than 10 mapped reads were subjected to differential expression analysis. The expression difference of each transcript between different samples was calculated based on a MARS (MA-plot-based method with Random Sampling model) model using the DEGseq package. FDR (false discovery rate) values less than 0.01 and |log2(fold change)|≥2 were considered as significant differences at expression level. To understand the dynamics changes during the different imbibition stages, hierarchical cluster analysis of expression patterns were performed by MultiExpreriment Viewer (v4.1) (Eisen et al., 1998). Using the STRING protein-protein interaction database, we can construct the network with some of the DEGs and TFs. The interaction networks were then visualizd by the software Cytoscape (Shannon et al., 2003).

Functional annotation of the differential unigenes were performed to search against the NR, Swiss-Prot, GO (Gene Ontology), and KEGG databases. WEGO software was used for functional classification of the results of GO annotations. GO enrichment of DEGs was carried out by ArigGO (Du et al., 2010). GO functional enrichment and KEGG pathway enrichment analysis were also tested at a significance cutoff of p-value. All the p-values were adjusted with the criterion of Bonferroni correction. We selected the corrected p-value of 0.05 as the threshold to determine significant enrichment terms of the gene sets.

For the quantitative RT-PCR of the DEGs, 2 μg of total RNA was used to synthesize the cDNA using the RevertAid First Strand cDNA Synthesis Kit (Fermentas). Quantitative RT- PCR was performed using the FastStart Universal SYBR Green Master (Roche) according to the manufacturer's instructions on the StepOne plus Real time PCR Platform (Applied Biosystems). The qRT-PCRs were carried out with the following protocol: 95°C for 10 min, followed by 40 cycles of 95°C for 15 s, and at 60°C for 60 s. The VgActin was used as the internal control, which has been shown to be one of the reference genes in mungbean. After amplification, the melting curve was determined for each specific product. Three independent biological replicates for each sample and three technical replicates for each biological replicate were analyzed. All the primers used were listed in Table S1. Significant differences in the expression level between control (ck) and experimantal groups (sy) were evaluated using Student's t-test.

All the sequencing data have been deposited in the NCBI short read archive (SRA) database under the accession number: SRP077637.

In many plants, seeds provide a convenient system for DT, even at the early stages of germination. With increasing imbibition time, the DT was gradually lost. In order to investigate the gene expression associated with DT, different imbibition times were characterized (Figure 1). For mungbean seeds, we observed that the germination rate after 3 h of imbibition was similar to that of seeds without imbibition (Figure 1A), showing 100% during all dehydration times. Germination rate of seeds after 6 h imbibition gradually decreased with increasing dehyration time, while the germination rates were sharply reduced for seeds after 18 and 24 h of imbibition (Figure 1A). In contrast, the moisture contents were high in seeds of the 18 and 24 h of imbibition (Figure 1B). Comparison of the radicles and germs, revealed that the lengths of 0, 3, and 6 h of imbibition had little difference, with significantly decreased in the seeds of 18 and 24 h imbibition (Figures 1C,D). These results suggested that the 3–18 h of imbibition were the key period for DT in mungbean. Therefore, studies on the expression patterns of these imbibition times will reveal the molecular mechanism of DT in mungbean.

To obtain a comprehensive transcriptome profile of DT, 24 libraries of mungbean seeds were constructed with or without dehydration from different time points. Totally, 128.7 Gb clean reads with an average of 4.03 Gb for each sample (Table S2). Using the de novo assembly program Trinity, the contig sequences were created. As a result, a total of 2,544,450 contigs from the 24 samples were assembled (Table 1). And then a total of 420, 228 transcripts and 92,155 unigenes were found. The assembled unigenes had a length distribution from 201 to 39,104 with an average length of 956.02 bp (Figure S1 and Table 1). For each time-point sample, the three biological replicates, showed a high correlation with each other (Figure S2). These assembled unigenes in this study will provide useful resources for further studies in mungbean.

For annotation, a similarity search of all unigenes were performed against protein sequences available at NCBI-NR, Swiss-Prot, Pfam and COG databases using BLASTx algorithms with an E-value threshold of 1e⁻⁵ (Table 2). As a result, a total of 62, 811 sequences were annotated, of them, 55,222 (83.90%), 32,041 (51.01%), 45,140 (71.87%), and 30,763 (48.98%) unigenes were aligned against the three protein databases, respectively. GO and KEGG analysis were also performed, generating 32,855 and 18,924 unigenes, respectively (Table 2). The length of these annotated unigenes was mainly ranged from 300 to 1000 bp (65.51%) (Table 2).

To investigate the gene expression patterns of seeds during the different stages of imbibition, FPKM values were used to normalize the reads from RNA-seq. Thus, the comparison of gene expression patterns between the control and experimental groups or among the experiment groups were performed (Figure 2A). Therefore, the numbers of shared DEGs in mungbean seeds between the adjacent stages after imbibition were 40,217 and 39 for 3 vs. 6 h, 6 vs. 18 h, and 18 vs. 24 h, respectively (Figure 2A). Since 3–18 h imbibition was the key peirod for mungbean seeds DT, the DEGs among these stages were further for in-depth analysis. A total of 121 DEGs were shared at the 3 vs. 6 h and 6 vs. 18 h of imbibition (Figure 2B). A heatmap of the shared DEGs were also showed the difference of their expression patterns, with significant divergence for 3, 6, 18, and 24 h of imbibition (Figure 2C). Clustering analysis of the total DEGs among 3, 6, and 18 h (3310 DEGs) revealed 12 different clusters with differentially expression patterns (Figure S3 and Table S3). Clustes 6, 11, and 13 showed highly expressed DEGs in the experimental groups. Particularly, most DEGs of the cluster 13 were induced at 3 and 6 h of imbibition (Figure 3 and Figure S3). While, other DEGs were particularly down-regulated at 3 and 6 h of imbibition (cluster 4 and cluster 8), suggesting that these DEGs might involved in DT of mungbean seeds (Figure 3).

Previous studies have reported that accumulation of protective molecules, including LEA proteins, maturation proteins and antioxidant systems contribute to DT (Maia et al., 2011; Huang et al., 2012). In this study, most of the LEA genes (c145025.graph_c0 and c166481.graph_c0) were up-regulated after 3 or 6 h of imbibition (sy3h and sy6h) (Figure 4, Table 4, and Table S4). And some of the peroxidases and multicopper oxidases showed differentially expression during the dehydration (Figure 4 and Table S4). Another important component of DT may be the accumulation of a high level of soluble sugars, which were implicated by correlation as adaptive agents for DT during seed development and germination (Blackman et al., 1992). In the present study, sugar transporters (c163377.graph_c0 and c176714.graph_c0) were up-regulated at 6 h of imbibition (Figure 4 and Table S4). Interestingly, some of the DEGs encoding methyltransferases were showed differentially expression at 3 or 6 h of imbibition (Table 4, Figure S4 and Table S4). And the expression levels of these DEGs were reduced at the stages (sy3h and sy6h), indicating that methylation invovled in the regulation of DT in mungbean (Figure S4).

Transcription factors (TFs) play crucial roles in the regulation of gene expression and plant stress responses. In the present study, we also found that many TFs, including HSP, ARF, AP2, MYB, GRAS, WRKY, and NAC, were differentially expressed during the dehydration (Figure 4, Table 4, and Table S4). Totally, 14 of the 18 HSP TFs were up-regulated after 3 h of imbibition (Figure 4 and Table S4). Many of the AP2 or ARF TFs, which have important functions in the transcriptional regulation of growth and development as well as stress response, were also differentially expressed during the dehydration process (Figure 4).

To infer the biological processes and functions of genes associated with DT, we conducted GO analysis of the DEGs at different stages. In total, 32,855 transcripts were assigned to 3458 functional terms and the GO terms of DEGs from the adjacent stages of the key period of DT were enriched (Figures 5, 6 and Table S5). The assigned functions of expressed transcripts covered a broad range of GO categories. Particularly, “extracellular region part” and “transcription regulator” were over-represented in DEGs of 3 vs. 6 h, while “cell killing,” “translation regulator,” and “growth” were the predominant groups for the DEGs at 6 vs. 18 h (Figure 5). GO enrichment analysis also showed that more DEGs were involved in “developmental process” and “response to stimulus” processes during early imbibition (3 vs. 6 h) (Figure 6A). Moreover, for 6 vs. 18 h after imbibition, “oxidation reduction,” “methylation,” and “cell wall polysaccharide metabolism” were the predominant categories (Figure 6B). Thus, the DEGs along with the time of dehydration treatment can be divided into two phases: response to stimulus, and oxidation reduction or metabolic process (Figure 6 and Table S5). The antioxidant systems and membrane structural stabilizing were essential to DT in resurrection plants (Moore et al., 2009).

In order to mine the genes involved in DT, 18,924 non-redundant potential protein sequences were searched against the KEGG database (Table S6). A total of 324 KEGG pathways were identified for the DEGs from the adjacent stages. The top 10 pathways enriched for the DEGs of 3 vs. 6 h and 6 vs. 18 h after imbibition were shown in Table 3. Most pathways were involved in metabolic processes, such as carbon metabolism (ko01200), starch and sucrose metabolism (ko00500) and amino sugar and nucleotide sugar metabolism (ko00520). However, there were some different KEGG pathways between 3 vs. 6 h and 6 vs. 18 h after imbibition. For the DEGs between 3 and 6 h after imbibition, pathways of fatty acid metabolism, pyruvate metabolism, and carbon fixation in photosynthetic organisms were enriched. Whereas, more DEGs were involved in the pathways of 6 vs. 18 h, and phenylpropanoid biosynthesis, phenylalanine metabolism, cysteine and methionine metabolism and ribosome pathways were overrepresented (Table 3). Sucrose was known membrane protectants, stabilizing the cellular processes, while some amino acid can be compatible solutes (Moore et al., 2009). Therefore, the results suggested different pathways were occurred during the loss of DT. Furthermore, COG analysis also showed that DEGs related to “transcription,” “cell wall/membrane,” and “lipid metabolism” were differentially expressed at 3–18 h of imbibition (Figure S5).

To assess the reliability of the RNA-seq data, quantitative RT-PCR analysis was performed using gene-specific primers of nine randomly selected genes. As expected, LEA protein (c166481.graph_c0) and seed maturation protein (c160456.graph_c0), which are involved in DT, showed highly expression levels at all the experiment groups. Sugar transporter (SUT, c163377.graph_c0) and copper/zinc superoxide dismutase (SODC, c166497.graph_c0) were induced at sy6h (Figure 7). Interestingly, photosystem II CP47 reaction center protein (c177934.graph_c1) was significantly highly expressed at 3 and 6 h of imbibition. Two enzymes, fatty acid hydroxylase (c168761.graph_c0) and aldehyde dehydrogenase (ADH, c169901.graph_c0) were also validated to have differentially expression patterns during the different stages (Figure 7). In addition, two TFs (HSP20, c87928.graph_c0 and AP2, c172061.graph_c0) were also up-regulated at 6 h of imbibition (Figure 7). The expression profiles of these nine detected genes showed the same trend and consistent results between the qRT-PCR and RNA-seq. Network analysis were also revealed that some of the DEGs were involved in diverse signaling pathways mediated by hormones (auxin and ABA) (Figure 8A). The TFs were closely connected with LEA and SUT, which might suggest a contribution to the regulation of DT in mungbean.

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