Reads mapping and analysis
The bamboo shoot development can be divided into six stages: dormancy, germination, development stage I, II and III, and shoot stage (33, 34) (Fig. 1). To characterize and generate the expression profiles for rhizome bud development, cDNA samples from rhizome lateral bud of Moso Bamboo in germination stage or early shoot stage were prepared and sequenced.
High throughout RNA-seq reads were mapped onto the published Moso bamboo draft genomic sequences (22) (Table 1). With a pre-treatment to the reads generated, there were 11,863,411 merged single reads and 3,885,776 paired-end reads in germination stage sample, and 20,854,206 merged single reads and 6,910,356 paired-end reads in shoot stage sample. Among them, 14,769,297 (94.9%) reads in germination stage sample and 25,909,031 (93.3%) reads in shoot stage sample can be mapped onto the genome sequences. Only 32.6% of mapped reads in germination stage sample and 31.9% of mapped reads in shoot stage sample were mapped onto the previously annotated genome regions.
We thus used Cufflinks to reannotate novel genes by integrating information on reads mapping patterns and previous annotation and obtained a total list of 55,869 genes (35). Accordingly, 13,939,159 reads (94.4%) in germination stage sample and 24,674,009 reads (95.2%) in shoot stage sample were mapped to these genetic regions, indicating that the reannotation covered the vast majority of reads. In gene saturation line (Fig. 2A), the increase of newly covered genes reached plateau as more reads were sequenced, which indicated that our sequencing depth was high enough to cover most genes expressing in our samples. The average reads covered per site of gene relative positions showed that although a part of sequences located to both terminates of genes were lost, our data had no obvious 5 prime or 3 prime bias (Fig. 2B and 2C).
Identification of novel genes in RNA-seq
We identified 22,107 novel genes (25,808 transcripts or 60,454 exons) within inter-genetic regions with an RPKM cutoff of 1.0 and compared the distributions of exon and intron length between the novel genes and previously annotated genes from published bamboo genome (Fig. 3). Both novel and annotated genes shared a similar distribution of exon length with two peaks at about 200bp and 700bp, consistent with recent data. However, there were larger gene density at 700bp in our newly identified genes due to a proportion of single long exon genes (Fig 3A). The distribution of intron length in novel genes and annotated genes were also similar, with a high sharp peak at about 100bp and a low peak at about 500bp (Fig. 3B). Such similarity in both exon and intron length between novel genes and old genes confirmed the validity of our newly identified genes.
Novel genes were grouped according to the number of exons, with 5,770 multiple exon genes (MEG) and 16,337 single exon genes. In MEG set, Trinity predicted 5,986 non-redundant coding sequences (36), 5,149 of which can be aligned to and annotated by Oryza sativa japonica genomic data with an e-value less than 1e-5. We clustered 3,153 newly annotated MEG and 23,213 previously annotated genes in Mapman. In each category, the number of newly annotated MEG were proportional to the number of previously annotated genes (~10%, Fig 3C). The majority of newly annotated MEG were involved in basal metabolism, protein and RNA synthesis and cell signaling. In addition, novel genes were identified in categories that were closely related to meristem development such as cell wall (26), cell division and organization (77) and hormone metabolism (44).
In the single exon genes set, 7,962 protein sequences were predicted by Trinity (36) and 1,737 non-coding RNA sequences were found by blast against Rfam (37), including 1,055 rRNA, 58 tRNA, 12 micro-RNA and 37 snoRNA (Table 2). The results above indicated that the sequences identified in this study captured the abundance of novel genes (both functional proteins and non-coding RNA) expressed in Moso bamboo rhizome lateral buds and early shoots. For most plants, only less than ~37% of the genes are intronless (38, 39). The MEG set was considered more reliable. Thus, we continued with MEG set for further functional analysis.
Identification of alternative splicing (AS) events
We also used our data to identify AS events in bamboo. A total of 28,217 AS events were identified by Astalavista-3.2 in the two samples (40) (Table 3), a number similar to previous study on AS events in rhizome bud tissue using PacBio technology (58). The intron retention events (IR), alternative 3’ splice sites events (Alt 3’ ss), exon skipping events (ES) and alternative 5’ splice sites events (Alt 5’ ss) constituted the most frequent AS events types, which also consistent with previous findings (56) (59). About 37.2.% and 32.2% of the AS junctions were previously annotated splicing sites in germination stage sample and shoot stage sample, respectively, while 38.0% and 40.6% were assigned into unknown splicing sites in annotated genes. The rest 24.9% and 27.2% junctions were located in novel genes. Together, we identified a total of 10,588 AS related splicing variants. 9263 variants originated from ~ 40% of the annotated genes from previously published bamboo genome (22), whereas the rest 1325 were from the novel genes.
AS creates protein isoforms and leads to increased protein complexity and functional diversity. A large portion of AS events in winter bamboo shoot development lead to domain lost in hormone associated genes, potentially regulating hormone signaling (57). Wang et al. divided alternative splicing events into two classes according to their impacts on domain architectures (41), altering the length of domain (Class I) and retaining/deleting the domain (Class II). To investigate if any protein families ‘prefer’ AS during RNA transcription in Moso Bamboo, we performed clustering analysis to our data. Proteins in PKinase and PKinase_Tyr, DEAD and PP2C families had enriched Class I AS, whereas WD40, TPR_11, IQ and zf-CCCH domain families were prone to Class II AS (Fig.4A).
Identification of single nucleotide polymorphisms (SNPs) and base insertions and deletions (Indels)
We next examined single-nucleotide polymorphisms (SNPs), and insertions and deletions (indels) presented in our data. There were 18602 homozygous SNPs and 25979 heterozygous SNPs compared to the draft genome. 44.2% homozygous SNPs and 67.8% heterozygous SNPs located at the coding regions of previously annotated genes. We observed more nonsynonymous SNPs (5019) than synonymous SNPs (3208) in homozygous SNPs, while the opposite was seen in heterozygous SNPs (10415 synonymous SNPs and 7206 nonsynonymous SNPs). The nonsynonymous SNPs were more frequently seen in proteins with PKinase domains, WD40 domains, nucleic acid-binding domains and Zinc finger FYVE domains (Fig. 4B).
We have also identified 6544 homozygous indels and 15347 heterozygous indels. Contrary to SNPs, only 5.34% (350) homozygous indels and 11.8% (1811) heterozygous indels located at annotated coding regions. It was also of notice that most of these indels (268 homozygous and 1430 heterozygous) did not cause frame shifts. After function clustering, we found that indels were enriched in protein kinases with ARM-like or RRM-dom domains (Fig. 4C). Premature termination codons (PTCs) often originate from nonsynonymous SNPs or frameshift. It is thus not surprising that protein families enriched in nonsynonymous SNPs and indels also possess PTC more frequently (Fig. 4C).
Genes with two or more nonsynonymous alleles are highly variant genes. Only less than 10% (1518) genes in our RNA-seq data were highly variant, suggesting a close relationship between the sequenced transcriptome and draft genome. These genes are enriched in processes of transportation, hormone metabolism, nucleotide metabolism, minor carbohydrate metabolism and photosynthesis, especially in ABC transport and multidrug resistance system and Jasmonate signaling pathway, indicating a faster gene evolution in these categories.
Functional analysis on differentially expressed genes
We next investigated how gene expression was differentiated between germination stage and shoot stage. With a RPKM cutoff of 1.0 (both RPKM values of genes in the two samples should be larger than the cutoff), we looked at how major pathways in plant growth change between the two stages (Fig. S1).
Transcription factors
Transcription factors regulate gene expression through binding to cis elements in response to environmental and developmental changes. It was estimated by Plant Transcription Factor Database that Moso Bamboo has about 1,768 transcription factors from 54 families. We found 138 transcription factors in 35 families displaying expression level changes between germination stage and early shoot stage. As high as three quarters of the transcription factors were highly expressed in germination stages including ARF, WRKY, NAC and DOF families, and the rest quarter transcription factors were highly expressed in early shoot stage, including NF-YA, HSF and MIKC type MADS families (Fig. 5A).
Rhizome meristem development
The underground development of the rhizome buds determines the size and yield of mature bamboo. The apical meristem maintains cell proliferation and differentiation for stem elongation and the auxiliary meristem differentiates and develops to branches, leaves and sheaths. We found several negative regulators of auxiliary meristem development LAX1, D14 and LOB were highly expressed during germination stage whereas genes promoting auxiliary meristem development (IPA1 and EP2) were upregulated during shoot stages (Fig. 5B). Our data suggested that auxiliary meristem development was inhibited during germination and promoted during shooting stage.
Carbohydrate metabolism
The rhizome buds are heterotrophs during the germination and early shoot stages, consuming the energy from metabolizing carbohydrates synthesized in above-ground tissues. We compared the gene expressions of the three main carbohydrate metabolism pathways in the two rhizome-bud developmental stages, glycolysis, tricarboxylic acid (TCA) cycle and pentose phosphate pathway (PPP) (Fig. 6).
Glycolysis is the major pathway for glucose metabolism in rhizome bud cells. The phosphorylation of glucose by hexokinase is the first rate-limiting and irreversible step, which initiates glycolysis. The resulting glucose-6-phosphate (G6P) is transformed to fructose-6-phosphate before going through another phosphorylation by ATP-dependent phosphofructkinase (PFK). Finally, a third rate-limiting kinase, pyruvate kinase, transfers the phosphate group from phosphoenolpyruvate to ADP and produces ATP and pyruvate, a TCA cycle precursor. We observed all six putative hexokinase related sequences were greatly expressed at germination stage and three of them were more than 2 G-fold over early shoot stage (Fig. 6A). Other two rate-limiting kinases also saw strongly (PFK, XLOC_038519, -2.23123) and moderately (pyruvate kinase, XLOC_034008, -1.395) increased activity during germination stage.
The pyruvate produced from glycolysis can form acetyl-CoA and go through a series of oxidization in TCA cycle, the most effective energy-producing pathway. We noted that a majority of putative genes in all steps of TCA cycle were upregulated during germination stage (Fig. 6B). A similar trend was also observed in PPP pathways, as the gene expression were elevated during germination stage, including the rate-limiting 6-phosphogluconate hydrogenase (Fig. 6C). Thus, the carbohydrate metabolism appeared more active during bud germination compared to early shoot development.
Hormone
Plant hormones play vital roles in the growth of bamboo shoot. We found putative genes annotated to indole-3-acettic acid (IAA), Brassinosteroid (BR) and gibberellin acid (GA) synthesis and signal transduction pathway displayed significant expression level changes during Moso Bamboo rhizome bud underground development.
IAA, the main auxin produced in plants, participates in cell proliferation, cell cycle regulation, cell elongation and cell wall development of the apical meristem. Using Trptophan as a precursor, IAA can be synthesized via multiple routes (Fig. 7). We found decreased expression level in putative gene annotated to major IAA synthesis pathways (CYP79B2 and CYP79B3 in IAOx route, IAM1 in IAM route and TAA1 in IPA route). The decline was further confirmed using reverse transcription qPCR to measure the expression of those putative genes throughout the five underground development stages of Moso bamboo rhizome bud (Fig. 7). It is of interest that the relevant gene expressions started to drop even at very early stages of underground development, prompting the question how IAA synthesis is regulated and how IAA level is maintained in rhizome buds.
BRs are a group of steroidal phytohormones that promote cell growth and stress responses, participating in all stages of plant development. Both our transcriptome analysis data and qPCR quantification showed that the putative genes involved in BR signaling pathway displayed systematic increase in expression in early shoot stages (Fig. 8). Thus, a strong role of BR in early shoot development is indicated.
GA stimulates seed germination and meristem development, and in later stages, triggers transition to flowering and fruiting. The transcriptional expression of GA3ox2 and GA20ox3, the key genes responsible for GA synthesis, were greatly elevated at early shoot stage, which is consistent with the trend seen in qPCR results (Fig. 9). The GA signaling activity leads to the decline of DELLA, a protein that inhibits cell growth and promote branching. Plants defective in GA synthesis caused accumulation of DELLA and displayed dwarf and branchy phenotypes. In our study, the DELLA gene expression decreased from germination to early shoot, suggesting an active GA signaling pathway during various underground stages of rhizome bud development (Fig. 9).
GA Stimulates Basal Internode Elongation in Moso bamboo
To confirm our findings from transcriptome analysis, we investigated the effects of GA signaling on fast growth in Moso bamboo. Moso bamboo seedlings were collected from Anji Bamboo Garden and grown in laboratory greenhouse till the establishment of five to seven true leaves before sprayed with GA, the growth inhibitor PAC and water, consecutively for 6 months (Fig. 10A).
The GA treatment significantly boosted the bamboo seedling vertical growth compared to the water treatment control (Fig. 10A, right), consistent with previous study (52). The PAC control, on the contrary, displayed darker leaves, shorter internode stem length and shorter distance from base to the first true leaf, indicating a slower stem elongation (Fig. 10A, middle). Neither treatment affected bamboo leaf width (LW), leave length (LL) and sheath length (SL) (Fig. S2, A and B).
We further dissected and compared the distribution of internode growth in different treatments (Fig. 10, B and C). The water treated wildtype had increased stem length for the first three basal internodes and the stem length sharply decreased as it came closer to the apical side. Treating GA did not disturb the basal stem growth pattern but resulted in longer and more-elongating basal internodes. PAC treatment completely abolished the basal stem growth in all internodes.
DELLA is a critical growth inhibitor during plant development. Same as in other monocot, we only detected one DELLA sequence in Phyllostachys edulis, SLR1 (Figure S3). In wildtype A. thaliana or DELLA defective mutant rga, the DELLA protein level remains undetectable as it is often rapidly degraded. However, an antibody against A. thaliana DELLA is able to detect the accumulation of DELLA in ga3, a mutant defective in GA synthesis, in Arabidopsis (Fig. S2C). Consistent with our previous findings, the accumulation of DELLA was only observed using the same antibody in PAC treated bamboo stem protein extract where the inhibition of growth was seen (Fig. 10D). Thus, GA and its signaling pathway is actively involved in Moso bamboo fast growth.