Ruibo Hu1
Changjiang Yu1Xiaoyu Wang1Chunlin Jia2Shengqiang Pei1Kang He1Guo He1Yingzhen Kong3
Gongke Zhou1*- 1Key Laboratory of Biofuels, Qingdao Engineering Research Center of Biomass Resources and Environment, Qingdao Institute of Bioenergy and Bioprocess Technology, Chinese Academy of Sciences, Qingdao, China
- 2Shandong Institute of Agricultural Sustainable Development, Jinan, China
- 3Key Laboratory of Tobacco Genetic Improvement and Biotechnology, Tobacco Research Institute of Chinese Academy of Agricultural Sciences, Qingdao, China
HIGHLIGHT
• De novo transcriptome profiling of five tissues reveals candidate genes putatively involved in rhizome development in M. lutarioriparius.
Miscanthus lutarioriparius is a promising lignocellulosic feedstock for second-generation bioethanol production. However, the genomic resource for this species is relatively limited thus hampers our understanding of the molecular mechanisms underlying many important biological processes. In this study, we performed the first de novo transcriptome analysis of five tissues (leaf, stem, root, lateral bud and rhizome bud) of M. lutarioriparius with an emphasis to identify putative genes involved in rhizome development. Approximately 66 gigabase (GB) paired-end clean reads were obtained and assembled into 169,064 unigenes with an average length of 759 bp. Among these unigenes, 103,899 (61.5%) were annotated in seven public protein databases. Differential gene expression profiling analysis revealed that 4,609, 3,188, 1,679, 1,218, and 1,077 genes were predominantly expressed in root, leaf, stem, lateral bud, and rhizome bud, respectively. Their expression patterns were further classified into 12 distinct clusters. Pathway enrichment analysis revealed that genes predominantly expressed in rhizome bud were mainly involved in primary metabolism and hormone signaling and transduction pathways. Noteworthy, 19 transcription factors (TFs) and 16 hormone signaling pathway-related genes were identified to be predominantly expressed in rhizome bud compared with the other tissues, suggesting putative roles in rhizome formation and development. In addition, a predictive regulatory network was constructed between four TFs and six auxin and abscisic acid (ABA) -related genes. Furthermore, the expression of 24 rhizome-specific genes was further validated by quantitative real-time RT-PCR (qRT-PCR) analysis. Taken together, this study provide a global portrait of gene expression across five different tissues and reveal preliminary insights into rhizome growth and development. The data presented will contribute to our understanding of the molecular mechanisms underlying rhizome development in M. lutarioriparius and remarkably enrich the genomic resources of Miscanthus.
Introduction
Currently, global climate change and depleting fossil fuel reserves have spurred increased concerns in the pursuit of alternative renewable energy resources such as biofuels (Yuan et al., 2008). Miscanthus is a temperate perennial C4 grass that belongs to the Andropogoneae tribe, which consists of many economically important crops including maize, sorghum, and sugarcane. Miscanthus has been regarded as a promising bioenergy crop for lignocellulosic biofuel production due to its superior characteristics such as high biomass yield, perennial growth habitus, low water and fertilizer requirement, abiotic stress tolerance and broad adaptation to marginal land etc. (Lewandowski et al., 2000; Parveen et al., 2013; Robson et al., 2013; Lee and Kuan, 2015; Xing et al., 2016). The Miscanthus genus contains more than 14 species, most of which originate from southern and eastern Asia (Hodkinson et al., 2002). Among these species, Miscanthus lutarioriparius that is endemic to central China, is characterized by overwhelmingly higher biomass yield and outstanding restoration capability of marginal land compared with M. sacchariflorus and M. sinensis (Liu et al., 2012; Yan et al., 2012; Mi et al., 2014). Thus, M. lutarioriparius is currently the most widely cultivated species for biomass production in bioenergy and biorefinery in China (Mi et al., 2014).
In particular, one of the most striking characteristics of M. lutarioriparius beyond its high productivity, is its conspicuous rhizome compared to other Miscanthus species. Rhizome cutting is the primary means of propagation for M. lutarioriparius because of its self-incompatibility. Rhizomes are modified subterranean stems, which have both apical and lateral meristems capable of generating adventitious shoots and roots (Yoshida et al., 2016). Rhizomes play important roles in the persistence of many perennial grass species, serving as the primary reservoir for energy storage and propagation organ (Jernstedt and Bouton, 1985; Li and Beuselinck, 1996; Glover et al., 2010). To some extent, rhizomatous growth is a desirable attribute in the establishment and persistence to various abiotic stresses for plants cultivated in adverse growth conditions (e.g., marginal land) (Sacks et al., 2006; Su et al., 2008; Zhou et al., 2014). However, because of the advantages of vigorous propagation, rhizomatous growth is also adopted by some of the most noxious weeds to grow vigorously and rapidly by strong rhizomes, which renders them to invade a wide range of regions and has raised increasing ecological concerns worldwide (Paterson et al., 1995).
Despite of the important roles of rhizome in plant growth and development, the genetic and molecular mechanisms underlying rhizome initiation and growth remain largely unknown. Rhizome development is a very complicated process, which is synergistically regulated by both intrinsic cues and environmental stimuli. Recently, with the advent of next-generation sequencing (NGS) technology, high-throughput RNA sequencing (RNA-seq) has become a powerful and cost-efficient means to discover putative functional genes involved in diverse biological processes, especially for plant species without a reference genome. RNA-seq has also been employed to identify genes related to rhizome development in various rhizomatous species including Oryza longistaminata (Hu et al., 2011), Sorghum (Sorghum halepense and Sorghum propinquum) (Jang et al., 2006; Zhang et al., 2014), reed (Phragmites australis) (He et al., 2012), bamboo (Phyllostachys praecox) (Wang et al., 2010), lotus (Nelumbo nucifera) (Cheng et al., 2013; Yang et al., 2015), Cangzhu (Atractylodes lancea) (Huang et al., 2016), Ginger (Zingiber officinale) and turmeric (Curcuma longa) (Koo et al., 2013) etc. In addition, transcriptome profiling of genes related to the rejuvenation of spring rhizomes in M. × giganteus and the rhizome-specific genes in M. sacchariflorus has also been carried out (Barling et al., 2013; Kim et al., 2014). For instance, several TFs belonging to YABBY, TCP and homeobox families were identified to be specifically expressed in rhizome tips of Oryza longistaminata, suggesting putative roles in rhizome development (Hu et al., 2011). PpHB1, a homolog of Homeodomain-leucine Zipper (HD-ZIP) transcription factor (TF), was implicated to be responsible for procambial development and rhizome bud formation in bamboo (P. praecox) (Wang et al., 2010). The rhizome morphogenesis of lotus (N. nucifera) was revealed to be mediated by photoperiod as well as light spectra (Masuda et al., 2007). Twenty-two genes associated with pathways in photoperiod, starch metabolism and hormone signaling were preferentially expressed in lotus rhizome, indicating possible roles in rhizome development (Yang et al., 2015). These genes provided excellent candidates for further functional characterization to unravel their roles in rhizome differentiation, growth and development.
Although M. lutarioriparius has been considered as an important feedstock resource for biofuel and biorefinery, the genetic resources especially transcriptome data are currently very limited, which significantly restricts the genetic improvement of this species for desirable traits. Previously, transcriptome profiling has only been reported for genes related to water-use efficiency and high photosynthesis efficiency for M. lutarioriparius (Fan et al., 2015; Xing et al., 2016). Rhizomnousness is one of the most distinctive characteristics of M. lutarioriparius compared to the other Miscanthus species. Rhizome cutting is the primary means of propagation for M. lutarioriparius in biofuel production. Elucidating the molecular mechanisms underlying rhizome initiation and development is of fundamental importance, because it will not only contribute to our better understanding of this important biological process, but also serves as theoretical basis for rapid and efficient propagation of this species through rhizome cutting. However, the transcriptome profiling of genes involved in rhizome development has not been addressed thus far in M. lutarioriparius.
In this study, we performed deep transcriptome sequencing using an Illumina HiSeq 2500 platform to characterize gene expression in rhizomes and other four different tissues of M. lutarioriparius. To the best of our knowledge, this is the first report on the transcriptome profiling of candidate genes related to rhizome development in M. lutarioriparius. These findings could contribute to our better understanding of the molecular mechanisms underlying rhizome initiation and development in M. lutarioriparius. In addition, the data generated also present as valuable resources for further functional genomics research in Miscanthus species.
Materials and Methods
Plant Materials
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.
RNA Isolation
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.
RNA Sequencing and De novo Assembly
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).
Functional Annotation of Unigenes
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−5.
Analysis of Differential Expression Genes (DEGs)
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 Plotting of DEGs
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.
GO and KEGG Enrichment Analysis
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.
Identification of Transcription Factors (TFs)
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−5.
Regulatory Network Construction
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).
Quantitative Real-Time RT-PCR (qRT-PCR) Analysis
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.
Results
Sequencing and De novo Assembly
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).
Functional Annotation and Classification
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−5. 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−45, while 42.5% of the mapped hit showed moderate homology with E-values between 1.0E−45 and 1.0E−5 (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.
Figure 1. Distribution of unigenes annotated to the NCBI NR protein database. (A) E-value distribution of top BLAST hits against the NR database for each unigene. (B) Similarity distribution of top BLAST hits for each unigene. (C) Species distribution of top BLAST hits for matched unigene sequences.
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.”
Figure 2. Histogram presentation of GO classification. Bars represent the numbers of unigenes matched to each GO term of three categories: Biological process (red), Cellular component (blue), and Molecular function (green).
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.
Figure 3. Histogram presentation of KOG classification. Unigenes with significant matches in the KOG database were classified into 26 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.
Figure 4. Classification of unigenes in KEGG pathways. The metabolism pathways are assigned into five categories. (A) Cellular processes. (B) Environmental information processing. (C) Genetic information processing. (D) Metabolism, and (E) Organismal systems.
Gene Expression Profiling
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).
Figure 5. Hierarchical clustering of differentially expressed genes and profile plots of clusters. (A) Hierarchical clustering of differentially expressed genes across five different tissues based on Z-score normalized FPKM values. Blue color represents lower expressed genes, while red color represents higher expressed genes. (B) Expression profile plots of 12 significant clusters. L, leaf; S, stem; R, root; LB, lateral bud; RB, rhizome bud.
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.
Tissue-Enriched and Differential Expressed Genes (DEGs)
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).
Figure 6. Comparative analysis of tissue-enriched and differentially expressed genes. (A) Pie chart showing the distribution of tissue-enriched genes in each tissue. (B) Venn diagram showing the up-regulated genes in tissues compared to rhizome bud. (C) Venn diagram showing the down-regulated genes in tissues compared to rhizome bud. L, leaf; S, stem; R, root; LB, lateral bud; RB, rhizome bud.
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.
GO and KEGG Enrichment Analysis
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.
Figure 7. Bar plots showing KEGG pathway enrichment of genes preferentially expressed in rhizome bud. Top 20 KEGG pathways with the most significant enrichment for each tissue compared to rhizome bud are shown as bar plots. L, leaf; S, stem; R, root; LB, lateral bud; RB, rhizome bud.
Identification of Hormone Signaling-Related Genes in Rhizome 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).
Figure 8. Expression profile of hormone signaling-related genes across five different tissues. The expression of genes involved in hormone signaling and transduction-related pathways was Z-score normalized and hierarchically clustered in heatmap. A color scale is shown at the top. Blue color indicates lower expression, while red color indicates higher expression. L, leaf; S, stem; R, root; LB, lateral bud; RB, rhizome bud.
Identification of Transcription Factors in Rhizome Development
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).
Figure 9. The differential expression of TFs across five tissues. The differential expression of TFs were depicted in heatmap based on Z-score normalized FPKM values. Blue color indicates lower expressed genes, while red color indicates higher expressed genes. L, leaf; S, stem; R, root; LB, lateral bud; RB, rhizome bud.
Figure 10. Profiling and regulatory network of TFs and hormone signaling genes highly expressed in rhizome bud. (A) Expression profile of 16 hormone signaling-related genes highly expressed in rhizome bud. (B) Hierarchical clustering of 19 TFs highly expressed in rhizome bud. (C) Predictive regulatory network among TFs and hormone signaling-related genes highly expressed in rhizome bud. The model was constructed based on gene co-expression and visualized by Cytoscape. TFs are in orange color and hormone signaling genes are in green color. L, leaf; S, stem; R, root; LB, lateral bud; RB, rhizome bud.
Validation of the Expression of DEGs by qRT-PCR
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 (R2 = 0.9219, P < 0.05) (Figure S4). These results verified the accuracy and reliability of the RNA-Seq analysis.
Discussion
As a promising lignocellulosic bioenergy plant, one of the superior characteristics of M. lutarioriparius is its vigorous rhizome as compared with the other Miscanthus species (e.g., M. sacchariflorus and M. sinensis). Rhizomatousness is an ideal trait for M. lutarioriparius to rapidly set up canopy and exert superior stress tolerances especially when cultivated in marginal land. In addition, the M. lutarioriparius rhizome also serves as the primary propagation organ as well as the main storage organ for reservoirs of nutrients in fall season. Thus, understanding of the biological processes involved in rhizome development and growth is of fundamental importance. In this study, we provided the first de novo transcriptome data from five tissues of M. lutarioriparius with the main objective to identify genes putatively involved in rhizome bud development.
By transcriptome profiling analysis, 3,062 genes were identified to be preferentially expressed in M. lutarioriparius rhizome bud (Figure 6C). KEGG enrichment analysis revealed that a large number of genes were involved in the pathways related to “phenylpropanoid biosynthesis,” “starch and sucrose metabolism,” “amino sugar and nucleotide sugar metabolism,” and “hormone signaling and transduction” (Figure 7). These pathways are principally involved in cell wall biosynthesis, cell proliferation, nutrient accumulation, primary metabolism and hormone signaling. The enrichment of these pathways in rhizome is consistent with the physiological roles of M. lutarioriparius rhizome in asexual propagation and serves as the storage organ of carbohydrate metabolisms in fall season. Transcriptome analysis has already been carried out in several rhizomatous species (e.g., reed, lotus, bamboo and sorghum) including two Miscanthus species (M. × gigantues and M. sacchariflorus) to identify candidate genes involved in rhizome development (Jang et al., 2006; Wang et al., 2010; He et al., 2012; Barling et al., 2013; Cheng et al., 2013; Kim et al., 2014; Zhang et al., 2014; Yang et al., 2015). A large number of highly expressed genes in rhizomes have been reported to be predominantly associated with “primary metabolism,” and “hormone signaling and stress response” etc. (Jang et al., 2006; Wang et al., 2010; He et al., 2012; Cheng et al., 2013; Yang et al., 2015). However, we did not find significant enrichment of pathways involved in response to stresses in M. lutarioriparius rhizome in the current study. This may be partly attributed to the situation that the M. lutarioriparius plants we used for this study were maintained in optimized growth conditions and irrigated regularly. Similar results were also reported for M. × giganteus transcriptome, in which genes involved in “hormone signaling” pathways were significantly enriched in spring rhizomes, whereas genes associated with “amino acid metabolism and seed maturation” pathways were much higher represented in fall rhizomes (Barling et al., 2013). Moreover, the highly expressed genes in rhizomes of several herbal plants such as Cangzhu (A. lancea), Ginger (Z. officinale), and turmeric (C. longa) were associated with pathways in secondary metabolisms (Koo et al., 2013; Huang et al., 2016). This is not unexpected considering the fact that rhizomes are the main storage organs for abundant amounts of secondary metabolisms in these rhizomatous herbal plants.
Although substantial progresses have been gained in our understanding of biological processes in rhizome development, the molecular mechanisms underlying rhizome development still remain largely unknown. However, as rhizome and tiller are originally derived from the axillary meristems from the lowermost part of the shoot, it can be assumed that genes controlling tiller development should also play a role in rhizome development and growth. This assumption is supported by the comparative transcriptome analysis of lateral bud and rhizome bud in this study. Among the five different tissues analyzed, lateral bud and rhizome bud shared the highest overall correlation coefficient in expression (Figure 5 and Figure S3). Correspondingly, these two tissues also shared the least number of DGEs compared to the comparisons among the other tissues (Figures 6B,C). These results suggested that a certain percent of genes might play largely overlapped roles in governing the morphogenesis and development of rhizome bud and lateral bud in M. lutarioriparius. However, their definitive roles awaits further functional characterization by transgenic studies.
TFs play important roles in various plant developmental processes. In the last two decades, studies especially in model species Arabidopsis and rice have revealed that a substantial numbers of TFs play crucial role in the regulation of tiller initiation and plant architecture. For example, TEOSINTE BRANCHED 1 (TB1) encoding a TCP domain transcription factor was originally identified as a negative player in regulating axillary bud outgrowth in maize (Doebley et al., 1997). The homologs of TB1 were also identified in other plant species including Arabidopsis (Aguilar-Martínez et al., 2007), rice (Takeda et al., 2003), and sorghum (Kebrom et al., 2006). REVOLUTA (REV), a member of HD-ZIP III TF gene family, is required for the initiation of axillary meristem in Arabidopsis (Talbert et al., 1995; Otsuga et al., 2001). LATERAL SUPPRESSOR (LAS) encoding a GRAS family transcription factor specifically regulates the initiation of axillary meristem in Arabidopsis (Greb et al., 2003). MONOCULM1 (MOC1), an ortholog of LAS in rice, is a key player involved in axillary bud initiation and tiller outgrowth (Li et al., 2003). Although a substantial set of genes involved in tiller development have been characterized and the underlying molecular mechanisms largely elucidated, it still remains to be answered to what degree does the molecular mechanism controlling tiller development relevant to that of rhizome growth. Interestingly, our transcriptome profiling revealed that a subset of TFs belonging to NF-YB, MYB and GRAS families were predominantly expressed in rhizome bud (Figure 10B). However, no homologs to these above mentioned TFs except for one GRAS family member were discovered in our analysis. We further analyzed the TFs that were specifically expressed in lateral bud. The results showed that several TFs homologous to TB1, REV, and LAS were preferentially expressed in lateral bud (Figure 9). These results indicated that the rhizome development of M. lutarioriparius is substantially different from that of tiller development.
In addition, accumulating evidence indicate that hormones affect diverse aspects of plant development processes, especially in the determination of tiller outgrowth and the formation of storage organs. For example, exogenous application of GA can significantly promote rhizome elongation in two rhizomatous grass species (Ma and Huang, 2016). In addition, genes associated with GA signaling pathways are enriched in rhizomes in several transcriptome analysis (Hu et al., 2011; Cheng et al., 2013). CTK plays important roles in the regulation of axillary bud initiation and outgrowth (Choi and Hwang, 2007), and exogenous application of CTK can stimulate rhizome formation in tall fescue (Festuca arundinacea) (Ma and Huang, 2016). IAA and CTK are implicated to act antagonistically in the regulation of axillary buds growth. Apically derived auxin inhibits the lateral bud outgrowth, whereas CTK relieves the inhibition. It is the homeostatic balance of their levels that determinates the initiation and outgrowth of axillary buds (Shimizu-Sato et al., 2009). In this study, 16 genes associated with ABA, IAA, BR, and GA signaling and transduction were preferentially expressed in rhizome bud, implicating the putative roles of these hormones in rhizome development in M. lutarioriparius. Furthermore, these hormone signaling genes were predicted to constitute a complicated regulation network with four NF-YB TFs identified in this study (Figure 10C). These genes identified provided excellent candidates for further functional characterization toward unraveling the molecular mechanisms underlying rhizome development in Miscanthus species.
Conclusions
In this study, a de novo assembly of transcriptome data from five tissues of M. lutarioriparius was performed with the main objective to provide preliminary insights into rhizome development. The comparative transcriptome analysis revealed differential and tissue-enriched genes, which were enriched in KEGG pathways associated with primary metabolism and hormone signaling. Noteworthy, 19 TFs and 16 hormone signaling and transduction-related genes were identified to be preferentially expressed in rhizome bud, highlighting the involvement of TFs and hormones in the regulation of rhizome development. Taken together, the transcriptome dataset presented here identified a subset of candidate genes putatively associated with rhizome formation and growth in M. lutarioriparius, thus laid a foundation for further functional genomics studies on rhizome development in Miscanthus species.
Author Contributions
RH designed the experiment, performed data processing and drafted the manuscript. CY and GH helped in bioinformatics analysis and data interpretation. XW, SP, and KH prepared the materials and performed the experiments. CJ participated in the design of the study, helped in data processing, and revision of the manuscript. YK assisted in results interpretation and manuscript preparation. GZ conceived the study and revised the manuscript. All authors read and approved the final version of the manuscript.
Conflict of Interest Statement
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.
Acknowledgments
This work was supported by the National Key Technology Support Program of China (2013BAD22B01), the Joint Funds of the National Natural Science Foundation of China (U1432126), the Science and Technology Major Project of Shandong (2015ZDJS03001), the Youth Innovation Promotion Association of CAS (2014187), the Taishan Scholar Program of Shandong (to GZ), and the Youth Talent Plan of Chinese Academy of Agricultural Sciences (to YK).
Supplementary Material
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017.00492/full#supplementary-material
Figure S1. Length distribution of assembled transcripts and unigenes.
Figure S2. GO enrichment of genes differentially expressed in rhizome bud compared to each of the other four tissues.
Figure S3. Pearson's correlation coefficient among different tissues.
Figure S4. Validation of expression between transcriptome profiling and qRT-PCR analysis for 24 genes.
Table S1. Primer sequences used in qRT-PCR analysis.
Table S2. Sequencing statistics of the transcriptome from five tissues.
Table S3. Annotation statistics of unigenes in publicly available databases.
Table S4. FPKM values of TFs and hormone signaling-related genes in five tissues.
Availability of Supporting Data
All the sequencing data used in this research are available from the NCBI Sequence Read Archive (SRA) database with accession number SRP094875.
References
Aguilar-Martínez, J. A., Poza-Carrion, C., and Cubas, P. (2007). Arabidopsis BRANCHED1 acts as an integrator of branching signals within axillary buds. Plant Cell 19, 458–472. doi: 10.1105/tpc.106.048934
Barling, A., Swaminathan, K., Mitros, T., James, B. T., Morris, J., Ngamboma, O., et al. (2013). A detailed gene expression study of the Miscanthus genus reveals changes in the transcriptome associated with the rejuvenation of spring rhizomes. BMC Genomics 14:864. doi: 10.1186/1471-2164-14-864
Cheng, L., Li, S., Yin, J., Li, L., and Chen, X. (2013). Genome-wide analysis of differentially expressed genes relevant to rhizome formation in lotus root (Nelumbo nucifera Gaertn). PLoS ONE 8:e67116. doi: 10.1371/journal.pone.0067116
Choi, J., and Hwang, I. (2007). Cytokinin: perception, signal transduction, and role in plant growth and development. J. Plant Biol. 50, 98–108. doi: 10.1007/BF03030617
Doebley, J., Stec, A., and Hubbard, L. (1997). The evolution of apical dominance in maize. Nature 386, 485–488. doi: 10.1038/386485a0
Fan, Y., Wang, Q., Kang, L., Liu, W., Xu, Q., Xing, S., et al. (2015). Transcriptome-wide characterization of candidate genes for improving the water use efficiency of energy crops grown on semiarid land. J. Exp. Bot. 66, 6415–6429. doi: 10.1093/jxb/erv353
Glover, J. D., Reganold, J. P., Bell, L. W., Borevitz, J., Brummer, E. C., Buckler, E. S., et al. (2010). Agriculture. Increased food and ecosystem security via perennial grains. Science 328, 1638–1639. doi: 10.1126/science.1188761
Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883
Greb, T., Clarenz, O., Schafer, E., Muller, D., Herrero, R., Schmitz, G., et al. (2003). Molecular analysis of the LATERAL SUPPRESSOR gene in Arabidopsis reveals a conserved control mechanism for axillary meristem formation. Genes Dev. 17, 1175–1187. doi: 10.1101/gad.260703
He, R., Kim, M. J., Nelson, W., Balbuena, T. S., Kim, R., Kramer, R., et al. (2012). Next-generation sequencing-based transcriptomic and proteomic analysis of the common reed, Phragmites australis (Poaceae), reveals genes involved in invasiveness and rhizome specificity. Am. J. Bot. 99, 232–247. doi: 10.3732/ajb.1100429
Hodkinson, T. R., Chase, M. W., and Renvoize, S. A. (2002). Characterization of a genetic resource collection for Miscanthus (Saccharinae, Andropogoneae, Poaceae) using AFLP and ISSR PCR. Ann. Bot. 89, 627–636. doi: 10.1093/aob/mcf091
Hu, F., Wang, D., Zhao, X., Zhang, T., Sun, H., Zhu, L., et al. (2011). Identification of rhizome-specific genes by genome-wide differential expression analysis in Oryza longistaminata. BMC Plant Biol. 11:18. doi: 10.1186/1471-2229-11-18
Huang, Q., Huang, X., Deng, J., Liu, H., Liu, Y., Yu, K., et al. (2016). Differential gene expression between leaf and rhizome in Atractylodes lancea: a comparative transcriptome analysis. Front. Plant Sci. 7:348. doi: 10.3389/fpls.2016.00348
Jang, C. S., Kamps, T. L., Skinner, D. N., Schulze, S. R., Vencill, W. K., and Paterson, A. H. (2006). Functional classification, genomic organization, putatively cis-acting regulatory elements, and relationship to quantitative trait loci, of sorghum genes with rhizome-enriched expression. Plant Physiol. 142, 1148–1159. doi: 10.1104/pp.106.082891
Jernstedt, J., and Bouton, J. (1985). Anatomy, morphology, and growth of tall fescue rhizomes. Crop Sci. 25, 539–542. doi: 10.2135/cropsci1985.0011183X002500030026x
Kebrom, T. H., Burson, B. L., and Finlayson, S. A. (2006). Phytochrome B represses Teosinte Branched1 expression and induces sorghum axillary bud outgrowth in response to light signals. Plant Physiol. 140, 1109–1117. doi: 10.1104/pp.105.074856
Kim, C., Lee, T. H., Guo, H., Chung, S. J., Paterson, A. H., Kim, D. S., et al. (2014). Sequencing of transcriptomes from two Miscanthus species reveals functional specificity in rhizomes, and clarifies evolutionary relationships. BMC Plant Biol. 14:134. doi: 10.1186/1471-2229-14-134
Koo, H. J., McDowell, E. T., Ma, X., Greer, K. A., Kapteyn, J., Xie, Z., et al. (2013). Ginger and turmeric expressed sequence tags identify signature genes for rhizome identity and development and the biosynthesis of curcuminoids, gingerols and terpenoids. BMC Plant Biol. 13:27. doi: 10.1186/1471-2229-13-27
Lee, W. C., and Kuan, W. C. (2015). Miscanthus as cellulosic biomass for bioethanol production. Biotechnol. J. 10, 840–854. doi: 10.1002/biot.201400704
Lewandowski, I., Clifton-Brown, J. C., Scurlock, J. M. O., and Huisman, W. (2000). Miscanthus: european experience with a novel energy crop. Biomass Bioenergy 19, 209–227. doi: 10.1016/S0961-9534(00)00032-5
Li, B., and Beuselinck, P. (1996). Rhizomatous Lotus corniculatus L: II. morphology and anatomy of rhizomes. Crop Sci. 36, 407–411. doi: 10.2135/cropsci1996.0011183X003600020032x
Li, X., Qian, Q., Fu, Z., Wang, Y., Xiong, G., Zeng, D., et al. (2003). Control of tillering in rice. Nature 422, 618–621. doi: 10.1038/nature01518
Liu, W., Yan, J., Li, J. Q., and Sang, T. (2012). Yield potential of Miscanthus energy crops in the Loess Plateau of China. GCB Bioenergy 4, 545–554. doi: 10.1111/j.1757-1707.2011.01157.x
Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expressiondata using real-time quantitative PCR and the 2−ΔΔCT method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262
Ma, X., and Huang, B. (2016). Gibberellin-stimulation of rhizome elongation and differential GA-responsive proteomic changes in two grass species. Front. Plant Sci. 7:905. doi: 10.3389/fpls.2016.00905
Maere, S., Heymans, K., and Kuiper, M. (2005). BiNGO: a Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics 21, 3448–3449. doi: 10.1093/bioinformatics/bti551
Masuda, J., Ozaki, Y., and Okubo, H. (2007). Rhizome transition to storage organ is under phytochrome control in lotus (Nelumbo nucifera). Planta 226, 909–915. doi: 10.1007/s00425-007-0536-9
Mi, J., Liu, W., Yang, W., Yan, J., Li, J., and Sang, T. (2014). Carbon sequestration by Miscanthus energy crops plantations in a broad range semi-arid marginal land in China. Sci. Total Environ. 496, 373–380. doi: 10.1016/j.scitotenv.2014.07.047
Otsuga, D., deGuzman, B., Prigge, M. J., Drews, G. N., and Clark, S. E. (2001). REVOLUTA regulates meristem initiation at lateral positions. Plant J. 25, 223–236. doi: 10.1046/j.1365-313x.2001.00959.x
Parveen, I., Wilson, T., Donnison, I. S., Cookson, A. R., Hauck, B., and Threadgill, M. D. (2013). Potential sources of high value chemicals from leaves, stems and flowers of Miscanthus sinensis ‘Goliath’ and Miscanthus sacchariflorus. Phytochemistry 92, 160–167. doi: 10.1016/j.phytochem.2013.04.004
Paterson, A. H., Schertz, K. F., Lin, Y. R., Liu, S. C., and Chang, Y. L. (1995). The weediness of wild plants: molecular analysis of genes influencing dispersal and persistence of johnsongrass, Sorghum halepense (L.) Pers. Proc. Natl. Acad. Sci. U.S.A. 92, 6127–6131. doi: 10.1073/pnas.92.13.6127
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Robson, P., Jensen, E., Hawkins, S., White, S. R., Kenobi, K., Clifton-Brown, J., et al. (2013). Accelerating the domestication of a bioenergy crop: identifying and modelling morphological targets for sustainable yield increase in Miscanthus. J. Exp. Bot. 64, 4143–4155. doi: 10.1093/jxb/ert225
Sacks, E., Dhanapala, M., Tao, D., Cruz, M. S., and Sallan, R. (2006). Breeding for perennial growth and fertility in an Oryza sativa/O. longistaminata population. Field Crop. Res. 95, 39–48. doi: 10.1016/j.fcr.2005.01.021
Shimizu-Sato, S., Tanaka, M., and Mori, H. (2009). Auxin-cytokinin interactions in the control of shoot branching. Plant Mol. Biol. 69, 429–435. doi: 10.1007/s11103-008-9416-3
Su, K., Bremer, D., Fry, J., and Keeley, S. (2008). Rooting characteristics and canopy responses to drought of turfgrasses including hybrid bluegrass. Agron. J. 100, 949–956. doi: 10.2134/agronj2007.0292
Takeda, T., Suwa, Y., Suzuki, M., Kitano, H., Ueguchi-Tanaka, M., Ashikari, M., et al. (2003). The OsTB1 gene negatively regulates lateral branching in rice. Plant J. 33, 513–520. doi: 10.1046/j.1365-313X.2003.01648.x
Talbert, P. B., Adler, H. T., Parks, D. W., and Comai, L. (1995). The REVOLUTA gene is necessary for apical meristem development and for limiting cell divisions in the leaves and stems of Arabidopsis thaliana. Development 121, 2723–2735.
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621
Wang, K., Peng, H., Lin, E., Jin, Q., Hua, X., Yao, S., et al. (2010). Identification of genes related to the development of bamboo rhizome bud. J. Exp. Bot. 61, 551–561. doi: 10.1093/jxb/erp334
Wu, J., Mao, X., Cai, T., Luo, J., and Wei, L. (2006). KOBAS server: a web-based platform for automated annotation and pathway identification. Nucleic Acids Res. 34, W720–W724. doi: 10.1093/nar/gkl167
Xing, S., Kang, L., Xu, Q., Fan, Y., Liu, W., Zhu, C., et al. (2016). The coordination of gene expression within photosynthesis pathway for acclimation of C4 energy crop Miscanthus lutarioriparius. Front. Plant Sci. 7:109. doi: 10.3389/fpls.2016.00109
Yan, J., Chen, W. L., Luo, F., Ma, H. Z., Meng, A. P., Li, X. W., et al. (2012). Variability and adaptability of Miscanthus species evaluated for energy crop domestication. GCB Bioenergy 4, 49–60. doi: 10.1111/j.1757-1707.2011.01108.x
Yang, M., Zhu, L., Pan, C., Xu, L., Liu, Y., Ke, W., et al. (2015). Transcriptomic analysis of the regulation of rhizome formation in temperate and tropical lotus (Nelumbo nucifera). Sci. Rep. 5:3059. doi: 10.1038/srep13059
Yoshida, A., Terada, Y., Toriba, T., Kose, K., Ashikari, M., and Kyozuka, J. (2016). Analysis of rhizome development in Oryza longistaminata, a wild rice species. Plant Cell Physiol. 57, 2213–2220. doi: 10.1093/pcp/pcw138
Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11:R14. doi: 10.1186/gb-2010-11-2-r14
Yuan, J. S., Tiller, K. H., Al-Ahmad, H., Stewart, N. R., and Stewart, C. N. Jr. (2008). Plants to power: bioenergy to fuel the future. Trends Plant Sci. 13, 421–429. doi: 10.1016/j.tplants.2008.06.001
Zhang, T., Zhao, X., Wang, W., Huang, L., Liu, X., Zong, Y., et al. (2014). Deep transcriptome sequencing of rhizome and aerial-shoot in Sorghum propinquum. Plant Mol. Biol. 84, 315–327. doi: 10.1007/s11103-013-0135-z
Keywords: Miscanthus lutarioriparius, rhizome formation, RNA-Seq, lateral meristem, transcription factor
Citation: Hu R, Yu C, Wang X, Jia C, Pei S, He K, He G, Kong Y and Zhou G (2017) De novo Transcriptome Analysis of Miscanthus lutarioriparius Identifies Candidate Genes in Rhizome Development. Front. Plant Sci. 8:492. doi: 10.3389/fpls.2017.00492
Received: 04 November 2016; Accepted: 21 March 2017;
Published: 12 April 2017.
Edited by:
Michael Deyholos, University of British Columbia, CanadaReviewed by:
Xun Xu, BGI-Shenzhen, ChinaLiangcai Peng, Huazhong Agricultural University, China
Juan Yan, University of Chinese Academy of Sciences, China
Copyright © 2017 Hu, Yu, Wang, Jia, Pei, He, He, Kong and Zhou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Gongke Zhou, zhougk@qibebt.ac.cn