Different segments of stem bark exhibit distinct morphological features
Ramie fibers continuously develop along the stem during plant’s growth, while the internodes of the stem show obvious elongation only until the plant is fully elongated. To analyze the developmental stages of fiber formation and the gene expression profiles, three parts of ramie’s shoot were harvested, including top bud (TB), internode elongating region (ER) of stem and internode fully elongated region (FER) of stem (Fig.1A). The top buds and the barks peeled off from both ER and FER regions were used for histological analysis and RNA extraction. The scheme for RNA-seq data analysis is illustrated in Fig. 1B. The cross and longitudinal sections of TB, ER and FER were analyzed (Fig.1C and D). In the TB sample, ramie has amphicribral vascular bundle, which is different from flax or hemp plants but is similar to woody plant with continuous cambia within and outside the vascular bundles (Fig. 1C). The vascular structure in TB is characteristic of multiple layers of primary phloem without obvious boundary between vascular bundles. In ER and FER, clear differences were observed between these two regions (Fig. 1D). Firstly, FER has thicker bark than ER, and FER barks consist of more enlarged cells and more layers of phloem tissues. Secondly, fiber cells show thicker cell wall in FER without an increase in cell size; the thickness of the fiber cell wall is about 5.38 µm in FER vs. 1.87 µm in ER (Supplemental table 1). Thirdly, the cell wall of the fibers from FER phloem contains more lignin than that from ER, which is indicated by stronger red color of safranin dye staining (Fig. 1D-b and 1D-d). The differences among these three samples indicate different developing stages of phloem fiber cells. Therefore, we used these samples for gene expression profiling attempting to identify genes important for fiber development in ramie.
Assembly of de novo transcriptome and identification of unigenes
Thirty-three RNA samples were collected and subjected to the next generation sequencing (NGS), and the RNA-seq data, including the 9 submitted SRA files (SRR9112644-SRR9112651), were analyzed. More than 5G sequences with clean bases from each sample was obtained, and thus the total analyzed clean bases were about 1.7E+11. The genome size of Zhongzhu No. 1 is approximately 340 Mb [3, 4]. Therefore, the depth of the RNA-seq data used in this study is sufficient for a high quality de novo assembly of transcriptome for the expressed genes from the top bud and stem bark tissues. The 10 species with the most matching reads to our RNA-seq data were shown in Fig. 2A. Among all the reads generated, 3048 reads match with those in Boehmeria nivea, and the highest matching ratio (28%) was found to be with Morus notabils. Overall, there were 59486 unigenes assembled with the length longer than 300 bp, 47016 unigenes longer than 500 bp, and 31395 unigenes longer than 1000 bp. The GC content distribution of all unigenes was shown in Fig. 2B, and two peaks appeared between the range of 30% and 45%. The correlation analysis showed that the three replicates of each sample were closely correlated (Fig. 2C). The detailed size distribution of all unigenes was illustrated in Fig. 2D and 2E. The sequence of each unigene was subsequently processed by blast to NR, SWISSPROT and KOG databases, respectively, and the annotations were obtained according to the most similar protein or gene with e<1e-5.
Identification of differentially expressed genes (DEGs) and expression patterns among TB, ER and FER
DEGs among the three tissues were identified following the scheme shown in Fig. 1B. When compared with TB, there were 4138 unigenes up-regulated and 6638 unigenes down-regulated in the ER, and 3853 unigenes up-regulated and 5075 unigenes down-regulated in the FER (Fig. 3A). The VENN diagram showed that the DEGs among these 3 samples were grouped in 6 distinct clusters (Fig. 3B). The heatmaps of the expression of these clustered genes were shown in Fig. 4A, and the schematic map of the expression patterns and GO analysis were illustrated in Fig. 4B.
The cluster 1 and 2 contain the most DEGs with 4354 up- and 2046 down-regulated unigenes only in TB (Figure 4A and 4B). The cluster 1 DEGs consist of the unigenes with higher expression level in TB but lower expression level in both bark regions. GO analysis showed that these DEGs are involved in meiotic chromosome segregation and cell division and stomatal or leaf development (Fig.4B). GO analysis of the DEGs in cluster 2 showed that up-regulated transcription factors or transcription processes and the plant-type secondary cell wall biogenesis are among the top categories. Fifty-five unigene contigs for cell wall components or cell wall biogenesis and modification related factors were identified in the cluster 2 (Table S2). These factors include Cellulose Synthase A Catalytic Subunit 3 and 8 (CesA 3 and 8), Fasciclin-like Arabinogalactan Protein (FLA), beta-galactosidase (BGAL), several pectinesterase/pectinesterase inhibitors (PMEs/PMEIs) and the enzymes for the synthesis of other cell wall components, such as glucuronoxylan glucuronosyltransferase, galacturonosyltransferase, endochitinase, callose synthase, xyloglucan glycosyltransferase, XTHs, etc. (Tab. S2).
The cluster 3 and 4 show the unigenes up- or down-regulated only in FER (Fig. 4A and 4B). In these two clusters, there were 93 unigenes down-regulated and 476 unigenes up-regulated only in the FER. In cluster 3, a small amount of unigenes for membrane construction were down-regulated in FER. In cluster 4, relatively more unigenes were up-regulated in FER comparing with the down-regulated unigenes in cluster 3. Among these up-regulated unigenes, ethylene signaling pathway genes were the most enriched unigenes. There were totally 39 transcription factors up-regulated only in FER, and 18 out of the 39 were ethylene activating unigenes (Tab. S3). The DEGs only in the ER were clustered in cluster 5 and 6. Interestingly, some phloem development related unigenes were found to be down-regulated only in ER when compared with those in both TB and CER (Fig.4B).
In addition to the expression patterns analyzed among TB, ER and FER, DEGs between TB and ER or FER were also analyzed and GO analyses were performed. The top 10 items of three GO terms were shown in Fig. S1 and S2. When compared with TB sample, the barks of FER showed gene expression patterns distinct from the barks of ER.
GO analysis of DEGs between ER and FER
There were 1628 up-regulated unigenes and 757 down-regulated unigenes identified in ramie’s bark of FER when compared with ER (Fig. 5). GO analysis shown in Fig. 6 revealed the top 10 up-regulated biological processes including phloem development, response to chitin, ethylene-activated signaling pathway, DNA replication, salicylic acid mediated signaling, defense response, protein transmembrane transport and vasculature development, and the top 10 down-regulated biological processes including cytoplasmic translation, tricarboxylic acid cycle, indole glucosinolate metabolic process, plant-type secondary cell wall biogenesis, etc.. The up-regulated genes in the activation of ethylene signaling pathway in FER is listed in Tab. S4. Overall 21 unigenes or contigs of 14 Ethylene Respond Factors (ERFs) were up-regulated in FER, which include ERF1, ERF1A, ERF1B, ERF2, ERF3, ERF5, ERF17, ERF22, ERF53, ERF61, ERF71, ERF109, PAR2-13 and PAP2-4 (Tab. S5).
KEGG analysis of DEGs between ER and FER
The KEGG analysis of total DEGs from FER vs. ER revealed additional information to the GO analysis. The KEGG analysis indicated that these DEGs are involved in the pathways of starch and sucrose metabolism, citrate cycle, nitrogen metabolism, cysteine and methionine metabolism, ribosome, diterpenoid biosynthesis, phenylpropanoid biosynthesis, DNA replication, cell cycle, etc. (Fig. 7 and Tab. S7).
From the KEGG analysis, we found that the expression of 23 unigenes encoding 11 enzymes in the starch and sucrose metabolisms differed between ER and FER. These enzymes include sucrose synthase (EC184.108.40.206), sucrose-phosphate synthase (EC220.127.116.11), bata-amylase EC18.104.22.168, endoglucanase (EC22.214.171.124), bata-glucosidase (EC126.96.36.199), glucan endo-1, 3-beta-glucosidase (EC188.8.131.52), glucose-6-phosphate isomerase (EC184.108.40.206), phosphoglucomutase (EC220.127.116.11), UTP-glucose-1-phosphate uridylyltransferase (EC18.104.22.168), trehalose phosphatase (EC22.214.171.124) and trehalase (EC126.96.36.199) (Fig. 8). Most of these enzyme-encoding unigenes were up-regulated in FER, which suggests that multiple pathways for free D-glucose production might be enhanced in FER. In addition, other sugar producing processes such as sucrose-6P, maltose and dextrin might also be enhanced in FER. The increase in these sugar precursors could be important in providing building materials for the secondary cell wall biogenesis in ramie.
More lignin accumulation in FER was observed by the staining with safranin dye (Fig. 1D-d), and KEGG analysis identified the up-regulation of enzymes responsible for phenylpropanoid biosynthesis in this region (Fig.7 and 9). In the pheylpropanoid biosynthesis pathway (KO00940), 16 up-regulated unigenes encode enzymes including peroxidase (EC188.8.131.52), trans-cinnamate 4-monooxygenase (EC184.108.40.206), anthranilate N-methyltransferase (EC220.127.116.11), flavonoid 3',5'-methyltransferase (EC18.104.22.168), beta-glucosidase (EC22.214.171.124), caffeylshikimate esterase and vinorine synthase, some of which were among the DEGs contributing to secondary cell wall synthesis (Fig. 9 and Tab. S7). Most enzyme-encoding genes involved in the biosynthesis of lignin were up-regulated in FER, and up-regulation of several PMEs in FER was also identified (Tab. S7).
Interestingly, in the diterpenoid biosynthesis pathway, unigenes encoding enzymes such as Ent-kaurenoic acid oxidase 2 (KAO2) (EC126.96.36.199) and gibberellins 20 oxidase (GA20ox) (EC188.8.131.52) for converting the precursors to active GA isoforms were up-regulated, while the transcript level of the enzyme gibberellins 2-beta-dioxygenese 8 (GA2ox8) (EC184.108.40.206) for inactivation of GAs was decreased in FER (Fig.10A). These results suggest that a higher accumulation of active GAs in FER than in ER. We subsequently determined 10 types of GA molecules, including GA 1, 3, 4, 7, 9, 15, 19, 20, 24 and 53, in TB, ER and FER samples by LC-MS-MS. All the GAs but GA4 were detected in ramie samples. The active form GA7 had a very low concentration, while GA1 was the most abundant active GA. GA9 was the most abundant precursor. Our results showed that FER samples had the most abundant GA precursors and active GA molecules, while the TB samples had lowest contents of GAs (Fig.10B).