M. lutarioriparius was propagated asexually via rhizomes from one plant collected from Changsha, Hu'nan province, China. Plants were grown in a growth chamber under 16 h light/ 8 h darkness photoperiod at 25°C. Leaves (fully expanded, upper apex), roots, stems (2nd top internode), rhizome buds (~3 cm) and shoot lateral buds (4th basal internode) were separately collected from 1-year-old potted plants. Each tissue has two replicates from six individual plants. All samples were immediately frozen in liquid nitrogen and stored at −80°C until use.

Total RNA was isolated using Trizol reagent (Invitrogen) according to the manufacturer's instructions. The extracted RNA was treated with DNase I (Promega) to remove the contaminated DNA. The quality of RNA was initially evaluated by electrophoresis in 1.5% agar gel, then quantified using the NanoDrop2000 spectrophotometer (ThermoFisher Scientific) and Agilent 2100 Bioanalyzer (Agilent Technologies). Only samples with the RNA integrity number (RIN) values higher than 8.0 were subjected to further analysis. Two biological replicates were used for RNA extraction and further transcriptome sequencing.

Ten cDNA libraries (two for each tissue) were constructed from purified mRNA using a NEBNext Ultra RNA Library Prep Kit following the manufacturer's protocols. The cDNA libraries were sequenced using the Illumina Hiseq 2500 platform by Novogene Company (Beijing, China). Raw reads obtained after sequencing were filtered by removing adapter containing reads, empty reads, and low quality reads (Q < 20). The remaining clean reads were de novo assembled into non-redundant unigenes using the Trinity program (Grabherr et al., 2011).

Assembled unigenes were annotated using BLAST alignment against public databases, including the NCBI non-redundant protein database (NR), NCBI non-redundant nucleotide database (Nt), Swiss-Prot, Protein Family (Pfam), Gene Ontology (GO), eukaryotic Orthologous Groups (KOG) database, the Kyoto Encyclopedia of Genes and Genomes (KEGG), and KEGG Ortholog database (KO) with E-value threshold of 10⁻⁵.

The expression level of each unigene was normalized and calculated as the value of fragments per kilobase of exon per million fragments mapped (FPKM) using the edgeR package (Robinson et al., 2010; Trapnell et al., 2010). The differentially expressed genes (DEGs) were determined by the criterion of fold change ≥2.0 and False Discovery Rate (FDR) ≤ 0.05. The tissue-enriched genes were defined as FPKM value is larger than 10 and the fold change of relative expression is higher than 3.0 compared to other tissues. Venn graph was drawn using an online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Heatmap was generated with Z-score normalized FPKM values of the DEGs using the online OmicShare tools (www.omicshare.com/tools). The K-means clustering was conducted based on Pearson correlation of gene expression profiles.

The GO enrichment analysis of DEGs was implemented by the GOseq R packages based on Wallenius non-central hypergeometric distribution (Young et al., 2010). The KEGG pathway enrichment analysis of DEGs was performed using KOBAS software (Wu et al., 2006). Both analysis was tested at a significance cutoff of FDR ≤ 0.05.

The TF families were identified by BLASTX against known plant TFs identified in PlnTFDB database (http://plntfdb.bio.uni-potsdam.de/v3.0) with E-value threshold ≤ 10⁻⁵.

TFs and hormone signaling genes preferentially expressed in rhizome buds were selected for the generation of regulatory network. The Pearson coefficient was calculated between any two genes based on their FPKM values from five different tissues. A threshold of coefficient higher than 0.85 was adopted to discriminate significant co-expression between gene pairs. The network was displayed using Cytoscape software (Maere et al., 2005).

Twenty four unigenes were selected for the expression verification by qRT-PCR analysis. All reactions were carried out in 96-well plates in the LightCycler 480 detection system (Roche) using the SYBR Premix Ex Taq II (TaKaRa) kit. Amplification procedure is 95°C for 30 s, followed by 40 cycles of 95°C for 5 s, 60°C for 15 s, and 72°C for 10 s. The expression was normalized using reference gene ACTIN11 and determined by the 2-deltadelta Ct method (Livak and Schmittgen, 2001). All reactions were performed with three replicates. All primers used are shown in Table S1.

To gain a comprehensive overview of M. lutarioriparius transcriptome, 10 cDNA libraries from leaf, stem, root, lateral bud and rhizome bud were constructed and sequenced using Illumina HiSeq 2500 platform. After filtering out adaptor sequences, ambiguous and low-quality reads, ~6.5 gigabase (GB) clean reads were acquired for each sample. The total clean reads summed to 66 GB for all samples (Table S2). Using the Trinity program, all clean reads were de novo assembled into 308,155 transcripts with an average length of 1,080 bp and an N50 length of 1,843 bp. These transcripts were further assembled into 169,064 unigenes. The size of the unigenes ranged from 201 to 17,395 bp, with a mean length of 759 bp and an N50 value of 1,339 bp. Among these unigenes, 35,599 (21.1%) were longer than 1,000 bp, and 100,775 (59.6%) were shorter than 500 bp (Figure S1).

To perform annotation of the assembled unigenes, BLAST search was carried out against seven public protein/nucleotide databases with an E-value threshold of 10⁻⁵. In total, 103,899 (61.5%) unigenes were successfully annotated in at least one of these seven databases, while ~38.5% unigenes remain unmapped. The detailed gene annotation info was provided as Table S3.

According to the E-value distribution of significant hits against the NR database, 57.5% of the matched sequences showed significant homology with E-values less than 1.0E⁻⁴⁵, while 42.5% of the mapped hit showed moderate homology with E-values between 1.0E⁻⁴⁵ and 1.0E⁻⁵ (Figure 1A). In terms of similarity distribution, 63.1% of the matched sequences had a similarity higher than 80%, while 26.9% of the hits showed a similarity ranging from 18 to 80% (Figure 1B). Further analysis of homologies among different plant species revealed that the annotated unigenes had the highest homology with sequences from Sorghum biocolor (29.5%), followed by Zea mays (15.4%) and Setaria italic (6.1%) (Figure 1C). The top two species with the highest hits match are from the Poaceae family, indicating that the assembly and annotation of M. lutarioriparius transcriptome is proper and reliable.

To gain insight into the functional categorization of the assembled unigenes, GO classification was performed based on the NR annotation. A total of 69,973 unigenes could be assigned to GO terms, which were classified into 49 functional groups under three principle categories, i.e., Biological process, Molecular function and Cellular components (Figure 2). Under the Biological process category, “cellular process,” “metabolic process,” and “single-organism process” were predominantly represented. Within the Cellular components category, “cell,” “cell part,” and “organelle” were the most highly represented categories. For the Molecular function category, the most abundant of genes were associated with “binding” and “catalytic activity.”

To evaluate the effectiveness of the annotation and the completeness of the transcriptome library, all the assembled unigenes were subjected to search against the Eukaryotic Orthologous Groups (KOG) database. Based on sequence homology, unigenes were assigned into 26 KOG categories (Figure 3). Clusters “general function predicted only” and “posttranslational modification, protein turnover and chaperones” represented two of the largest ones, followed by “translation, ribosomal structure and biogenesis” and “signal transduction mechanisms.” In contrast, the clusters of “cell mobility” and “unnamed proteins” represented the smallest categories.

KEGG pathway analysis was performed to identify the biochemical pathways in M. lutarioriparius. Based on sequence homology, a total of 34,668 (20.5%) unigenes were matched in the KEGG database and assigned to 32 pathways, covering five major KEGG categories (Figure 4). The three most represented pathways were “translation,” “signal transduction,” and “carbohydrate metabolism,” followed by “carbohydrate metabolism” and “overview,” whereas “signaling molecules and interaction” and “sensory system” pathways represented the smallest categories.

To gain a global view of gene expression pattern across different tissues, a heatmap was plotted using the Z-score normalized FPKM values (Figure 5A). All biological replicates were clustered together with almost identical expression profiles, indicating the reliability of sample collection and analytical procedure. In addition, the overall expression of rhizome buds and lateral buds samples correlated well compared to the correlation between the other samples, which suggest that these two tissues are much closely related (Figure S2).

To further investigate gene expression profiles across different tissues, K-means clustering was carried out to classify differentially expressed genes (Figure 5B). Totally, 12 clusters with distinctive expression patterns were identified. The most abundant group is cluster 12 comprised of 2,222 genes, whose expression was highest in lateral bud, stem and rhizome bud. The second largest group is cluster 10, which contained 1,889 genes, whose expression is much lower in root. In contrast, the collection of genes in cluster 7 and 9 mostly showed relatively higher expression in root. The collection of genes in cluster 1, 2, and 3 was predominantly expressed in leaf. Cluster 11 genes exhibited higher expression in both lateral bud and rhizome bud. These different gene expression profiles indicate that each tissue was associated with specific gene clusters.

Knowledge of genes that are specifically or preferentially expressed in tissues can provide insights into the specialized developmental processes in these tissues. From this perspective, we sought to investigate the Differentially Expressed Genes (DEGs) in each tissue examined. All unigenes that have FPKM values higher than 10 and the fold change of relative expression is higher than 3.0 compared to other tissues were deemed as tissue-enriched genes. Totally, 4,609, 3,188, 1,679, 1,077, and 1,218 genes were identified to be specifically expressed in root, leaf, stem, rhizome bud and lateral bud, respectively (Figure 6A).

In order to identify genes specifically involved in the process of rhizome initiation and development, we further focused on the identification of DEGs in rhizome bud. As shown in the venn diagram of the distribution of DEGs in different tissues compared to rhizome bud (Figures 6B,C), root shared the largest number of differential expressed genes compared to rhizome bud with 42,443 genes up-regulated and 13,979 genes down-regulated, respectively. In contrast, lateral bud shared the smallest number of DEGs compared to rhizome bud, with 6,838 up-regulated and 9,658 down-regulated genes, respectively. In total, 1,004 common genes were up-regulated whereas 3,062 genes were down-regulated in the other four tissues compared to rhizome bud. These DGEs in rhizome buds may have important functional implications in rhizome development.

To gain insights into the functional categorization and metabolic pathways involved in rhizome development, DEGs in rhizome bud were subjected to enrichment analysis based on GO and KEGG pathways. Significantly enriched GO terms were found between rhizome bud and the other four tissues. Comparisons of these enriched GO terms in rhizome bud compared to leaf, lateral bud, stem and root revealed that the numbers and types of GO terms varied substantially (Figure S3). For example, “response to stimulus” and “DNA binding” were significantly enriched in rhizome bud compared to leaf, whereas “carbohydrate metabolic process” and “catalytic activity” were significantly enriched in a comparison between rhizome bud and stem. Nevertheless, several commonly enriched GO terms in rhizome bud compared to the other tissues, such as “cell,” “cell part,” “intracellular part” in Biological process category, “cellular process” and “metabolic process” in Cellular components category, and “binding,” and “catalytic activity” in Molecular functions category.

We further mapped the DEGs in rhizome bud to the KEGG database and analyzed the enrichment of metabolic pathways. The top 20 enriched pathways in rhizome bud compared to the other tissues were listed in Figure 7. Among these pathways, a large number of genes were involved in the pathways related to “phenylpropanoid biosynthesis,” “starch and sucrose metabolism,” and “amino sugar and nucleotide sugar metabolism.” Remarkably, a number of genes were enriched in “plant hormone signal transduction” pathways in rhizome bud, suggesting that hormones may play crucial roles in the regulation of rhizome bud development.

Hormones have been implicated to paly crucial roles in diverse aspects of plant development processes. To gain insights into the functional roles of hormones during M. lutarioriparius rhizome development, we mapped the DGEs to hormone signaling and transduction pathways and analyzed their expression in different tissues (Figure 8). A total of 175 genes were identified to be associated with the biosynthesis, metabolism and signaling of eight hormones, including abscisic acid (ABA), auxin (IAA), ethylene (ETH), brassinosteroid (BR), jasmonic acid (JA), cytokinin (CTK), gibberellin acid (GA), and salicylic acid (SA). The genes associated with IAA biosynthesis and metabolism represented as the largest group with 43 members, followed by ABA signaling-related genes with 40 members. In contrast, the genes associated with GA signaling was the smallest group with only five members. All these hormone signaling and transduction-related genes exhibited differential expression across the different tissues examined. Noteworthy, a subset of 16 genes showed preferential expression in rhizome bud, suggesting a putative role of hormones signaling in the regulation of rhizome growth and development. Among these genes, seven genes were associated with IAA biosynthesis and metabolism, followed by six genes related to ABA, and two genes were involved in BR signaling (Figure 10A and Table S4).

Accumulating evidence indicate that transcription factors (TFs) play critical roles in various plant development processes. To provide insights into the regulatory network underlying rhizome development, we examined the expression of TFs in different tissues especially their dynamic expression in rhizome bud. Totally, 690 TFs belonging to 44 different families were found to be differentially expressed in the tissues examined (Figure 9). The TFs displayed distinct expression patterns across the five tissues. A majority of TFs showed relatively broad expression patterns in all the tissues, while some exhibited distinctive tissue-specific patterns. Noticeably, 19 TFs were preferentially expressed in rhizome bud, implying that they might play important roles in the regulation of rhizome bud initiation and growth (Figure 10B and Table S4). These rhizome-specific TFs belong to different families with the largest members from NUCLEAR FACTOR Y (NF-Y) family. Furthermore, we constructed a putative regulation network based on gene co-expression patterns for the TFs and hormone signaling-related genes preferentially expressed in rhizome bud. The model indicated that four NF-YB TFs may putatively regulate the expression of six hormone signaling genes related to ABA and IAA. (Figure 10C).

To validate the reliability of the expression profiling obtained by RNA-seq, 24 genes with differential expression in rhizome bud were selected for quantitative real-time RT-PCR (qRT-PCR) analysis. The genes chosen for qRT-PCR analysis included 12 TFs and 12 hormone signaling-related genes. For all these genes, the results of qRT-PCR exhibited almost identical expression patterns as compared to the transcriptome profiling results. In addition, good agreement between qRT-PCR and RNA-seq results was also supported by a significant positive correlation between them as revealed by Pearson correlation analysis (R² = 0.9219, P < 0.05) (Figure S4). These results verified the accuracy and reliability of the RNA-Seq analysis.

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.