RNA-seq and phytohormone analysis reveals the culm color variation of Bambusa oldhamii Munro

Plant materials

The middle and lower culm internode epidermis samples that were removed from B. oldhamii (LZ) were labeled LZ_1, LZ_2, and LZ_3, and those from B. oldhamii f. revoluta W.T. Lin & J. Y. Lin (HLZ) were labeled HLZ_1, HLZ_2, and HLZ_3, representing three biological replicates of each type of bamboo. The culm skin samples were frozen in liquid nitrogen immediately for further phytohormone detection, RNA-seq, and relative gene expression. Total RNA was isolated using plant RNA isolation kits (Tiangen Biotech, Beijing, China).

Library construction and sequencing

Library construction and sequencing steps were performed based on the Illumina HiSeq platform for RNA-seq protocols (https://www.illumina.com) at Wuhan Metware Biotechnology Co., Ltd. (Wuhan, China). The output data contained raw reads in fastq format, that were then processed for quality control, including filtering and trimming of low confidence bases, biased nucleotide composition, adapters, duplicates and low-quality reads to acquire clean reads. Trinity (2.6.6) (Grabherr et al., 2011) and Corset (1.07) (Davidson & Oshlack, 2014) were used to assemble the clean data and process the transcript cluster analysis, and the longest transcripts in each cluster were filtered out as unigenes. All of the following analyses were performed using unigene sequences.

Gene annotation

Nr (NCBI nonredundant protein sequences), Pfam (Protein family), KOG (Protein family), Swiss-Prot, Trembl, KEGG (Kyoto Encyclopedia of Genes and Genomes), and GO (Gene Ontology) were used for gene annotation. The Nr, KOG, Swiss-Prot, Trembl, KEGG, and GO annotations were performed using BLAST (v2.7.1) with an e-value = 1e−5, and the Pfam annotation was performed using the hmmscan command of the HMMER 3.2 package with e-value = 0.01. Transcription factor annotation was performed with iTAK (1.7a) (Zheng et al., 2016) with default parameters.

Gene expression analysis

We used the assembled transcriptome of Trinity as a reference and then mapped the clean reads of each sample to the reference with RSEM software. FPKM (Fragments Per Kilobase of transcript per Million fragments mapped) values were calculated to estimate gene expression and abundance after normalization of the mapped reads and transcript lengths. The R package Pheatmap was used to draw a heatmap with normalized log2(FPKM+1) data and clusters of expression patterns with kmeans_k = 10. The color from red to blue indicates gene expression from high to low.

Differential expression analysis

After acquiring the abundance information and performing normalization, gene expression between the groups was compared. DESeq2 (1.22.2) was used to calculate the differentially expressed genes between the LZ and HLZ groups, which were corrected with FDR (False Discovery Rate) by Benjamini–Hochberg methods. Differentially expressed genes were filtered with the condition of —log2(Fold Change)— ≥ 1 and FDR < 0.05.

Validation of RNA-seq analysis via qRT-PCR

The RNA-seq results were validated for selected genes using qRT-PCR assays. cDNA was synthesized with HiScript^® II Q RT SuperMix for qPCR (Vazyme, China). Quantitative Real-Time PCR (qRT-PCR) was performed on a LightCycler^® 480 II Real-Time PCR system (Roche International Diagnostics system, Switzerland) using the Unique Aptamer™ qPCR SYBR^® Green Master Mix. The components of the qRT-PCR were as follows: SYBR Premix Ex Taq (2x) (10 µl), forward primers (0.5 µM), reverse primers (0.5 µM), cDNA template (2 µl), and ddH₂O to 20 µl. Then, qRT-PCR was performed as follows: initial denaturation at 95 °C for 5 min; 40 cycles of denaturation at 95 °C for 10 s and annealing at 72 °C; and finally, steps for melt-curve analysis (95 °C for 15 s, 60 °C for 60 s, 95 °C for 15 s). Actin was used as the internal control (Zeng et al., 2015; Table S1), and relative expression was calculated with the 2^−ΔΔCT method (Livak & Schmittgen, 2001).

Detection of phytohormone contents

The quantification of endogenous abscisic acid (ABA), auxin (indole-3-acetic acid, IAA), jasmonic acid (JA), and gibberellic acid (GA₁, GA₃, GA₄, and GA₇) was performed by Genepioneer Biotechnologies Co., Ltd. (Nanjing, China) using an HPLC–MS/MS platform.

Statistical analysis

The enrichment of up- and downregulated genes was determined using GOseq and KOBAS (Mao et al., 2005). The cor.test function was used to calculate the correlation between phytohormone contents and gene expression within the corresponding periods. Bar chart data of the phytohormone contents are reported as the mean ± SEM (n = 9) with a significant difference (p < 0.05) according to unpaired t-tests.

Plant materials

B. oldhamii (LZ) is a species of clumping bamboo (Fig. 1) and it has entirely green culms. The B. oldhamii variety referred to as B. oldhamii f. reboluta W.T. Lin & J. Y. Lin (HLZ) has green culms with yellow stripes of random widths. The culm skin was removed from LZ and HLZ to research the correlation of culm color variation on the phytohormone contents and gene expression levels.

Transcriptome sequences and data output

Total RNA was isolated from the culm epidermis samples of B. oldhamii with three biological replicates marked LZ_1 - LZ_3 and from B. oldhamii f. revoluta with three biological replicates marked HLZ_1 - HLZ_3. After the cDNA library was constructed and sequenced, approximately 282 million raw sequence reads were obtained from the RNA-seq experiment, and 267 million clean sequence reads remained after filtering with a Q20 above 98% after quality control was performed. The error correction and GC content are shown as follows (Table 1).

De novo assembly

The clean data were used for de novo assembly with Trinity (Grabherr et al., 2011), and overall, 604,900 transcripts and 383,278 unigenes were generated (Table 2). The average length of the transcripts was 670 bp, with an N50 of 1,033 bp and an N90 of 264 bp. Most (80%) of the transcripts were between 200 and 1,000 bp, while the remaining 20% of the transcripts had a length longer than 1,000 bp. The average length of the unigenes was 906 bp, with an N50 of 1,221 bp and an N90 of 433 bp. About 70% of the unigenes were shorter than 1,000 bp, while the remaining 30% of the unigenes were longer than 1,000 bp.

Gene annotation and functional classification

All of the unigenes were annotated using seven databases (Table 3) (https://doi.org/10.6084/m9.figshare.16912324). The annotation results produced 38.64%, 70.18%, 44.92%, 69.23%, 40.03%, 57.24%, and 47.26% unigenes annotated in the KEGG, NR, SwissProt, Trembl, KOG, GO, and Pfam databases, respectively. Approximately 274,681 unigenes were annotated in at least one database.

The unigenes annotated with GO functions were assigned to three main ontologies: molecular function (MF), cellular component (CC), and biological process (BP). The terms cellular process (GO:0009987), metabolic process (GO:0008152), biological regulation (GO:0065007), and response to stimulus (GO:0050896) were the most common BP ontologies; the terms cell (GO:0005623), cell part (GO:0044464), and organelle (GO:0043226) were the most common CC ontologies; and binding (GO:0005488), catalytic activity (GO: 0003824), transporter activity (GO:0005215), and transcription regulator activity (GO:0140110) were the most common MF ontologies (Fig. 2A; Table S2). The unigenes with KOG annotation were categorized as posttranslational modification, protein turnover, and chaperones; signal transduction mechanisms; translation, ribosomal structure and biogenesis; energy production and conversion; transcription; intracellular trafficking, secretion, and vesicular transport (Fig. 2B; Table S3). The majority of the annotated genes were characterized as ribosome pathway (ko03010) (1,841), glyoxylate and dicarboxylate metabolism (ko00630) (1,778), metabolic pathways (ko01100) (1,774), biosynthesis of secondary metabolites (ko01110) (1,768), oxidative phosphorylation (ko00190) (1,764), phenylpropanoid biosynthesis (ko00940) (1,755), carbon metabolism (ko01200) (1,742), glycolysis/gluconeogenesis (ko00010) (1,723), MAPK signaling pathway-plant (ko04016) (1,720), biosynthesis of amino acids (ko01230) (1,714), and pyruvate metabolism (ko00620) (1,710) (Table S4).

Expression patterns of differentially expressed genes (DEGs)

Differential expression analysis between the LZ and HLZ groups revealed 449 differentially expressed genes; approximately 190 DEGs were downregulated in the HLZ group, and 259 DEGs were upregulated in HLZ group. The downregulated genes in HLZ were classified into 10 clusters (Fig. 3A; Table S5). The results showed that cluster 2 included 1 gene (unigene-21666.123476) annotated with metallothionein that was highly expressed in all six samples, and cluster 1 included 8 genes that were more highly expressed in LZ samples. unigene-21666.69211 (mitogen-activated protein kinase kinase kinase ANP1), unigene-21666.134631 (SAUR family protein), unigene-21666.128797 (EREBF-like factor), and others without a specific annotation were included. The upregulated genes in HLZ were classified into 10 clusters (Fig. 3B; Table S6), where cluster 5 included 2 genes, and cluster 1 contained 9 genes that might play a vital role, including unigene-21666.113648 (phenylalanine ammonia-lyase), unigene-21666.90795 (phenylalanine ammonia-lyase); unigene-21666.169778 (serine/threonine-protein kinase PBS1), unigene-21666.167280 (granule-bound starch synthase), unigene-21666.149531 (acyl-[acyl-carrier-protein] desaturase), unigene-21666.125526 (anthranilate O-methyltransferase), unigene-21666.107961 (serine/threonine-protein kinase PBS1), and others without specific annotations (Fig. 3). Random DEGs were chosen for qRT-PCR analysis to validate the accuracy of the RNA-Seq data results. The relative expression results showed a strong correlation between the RNA-Seq and qRT-PCR data (Fig. 4; Table S7).

Functional classification of all DEGs

All 449 DEGs were processed for GO, KOG, and KEGG classification. The cellular process (GO:0009987), metabolic process (GO:0008152), response to stimulus (GO:0050896), and biological regulation (GO:0065007) terms were the most common BP ontologies; cell (GO:0005623), cell part (GO:0044464), organelle (GO:0043226), and membrane (GO:0016020) were the most common CC ontologies; and binding (GO:0005488) and catalytic activity (GO:0003824) were the most common MF ontologies (Fig. 5A; Table S8). For the KOG classification, the majority of the DEGs were involved in signal transduction mechanisms, followed by posttranslational modification, protein turnover, chaperones, transcription, translation, ribosomal structure, and biogenesis (Fig. 5B; Table S9). Based on the KEGG classification results, the DEGs were mostly involved in metabolic pathways (ko01100), followed by biosynthesis of secondary metabolites (ko01110) and plant-pathogen interactions (ko04626) (Fig. 5C).

Functional enrichment of up- and downregulated DEGs

Differential expression analysis between the LZ and HLZ groups showed that 190 and 259 DEGs were down- and upregulated in the HLZ culm skin samples. The GO functional enrichment revealed that the downregulated DEGs were more enriched in monocarboxylic acid transport (GO:0008028), transporter complex (GO:1990351), transmembrane transporter complex (GO:1902495), replication fork (GO:0005657), organelle envelope lumen (GO:0031970), lipid transporter activity (GO:0005319), and monocarboxylic acid transmembrane transporter activity (GO:0008028) (Fig. S1; Table S10).

The upregulated DEGs were more enriched in biological processes and molecular functions. In particular, the top three terms were L-phenylalanine metabolic process (GO:0006558), L-phenylalanine catabolic process (GO:0006559), and entrainment of the circadian clock (GO:0009649), followed by photoreceptor activity (GO:0009881), phenylalanine ammonia-lyase activity (GO:0045548), NAD(P)H dehydrogenase (quinone) activity (GO:0003955), and carbon-nitrogen lyase activity (GO:0016840) (Fig. S2; Table S11).

The downregulated DEGs were enriched in the alpha-linolenic acid metabolism pathway (ko00592) (Fig. S3), and the upregulated DEGs were enriched in the unsaturated fatty acids biosynthesis (ko01040), linoleic acid metabolism (ko00591), biosynthesis of secondary metabolites (ko01110), phenylpropanoid biosynthesis (ko00940), and phenylalanine metabolism (ko00360) pathways (Fig. S4).

Phytohormones control culm color variation

To reveal how phytohormones control culm color variation in B. oldhamii, the content of endogenous abscisic acid (ABA), auxin (indole-3-acetic acid, IAA), jasmonic acid (JA), and gibberellic acid (GA₁, GA₃, GA₄, and GA₇) was detected (Fig. 6A). The majority of phytohormones were highly accumulated in LZ, especially GA₁ and GA₇, which were significantly accumulated in LZ; the contents of ABA were relatively highly accumulated in HLZ (Fig. 6A; Table S12). The phytohormones of GA₁, GA₇ and ABA were choosen for further study, and the relationship between phytohormones (ABA, GA₁, and GA₇) and gene expression was analyzed with corresponding samples. There were 18 genes whose expression was significantly related to ABA accumulation patterns (Table S13), while 35 and 91 genes were significantly related to GA₁ (Table S14) and GA₇ (Table S15) accumulation patterns respectively. The genes significantly related to ABA, GA₁ and GA₇ were more enriched in the GO terms of chloride channel complex (GO:0034707), ion channel complex (GO:0034702), protein kinase complex (GO:1902911), transporter complex (GO:1990351), anion channel activity (GO:0005253), chloride channel activity (GO:0005254), phosphatidylinositol bisphosphate binding (GO:1902936), and phosphatidylinositol-3,5-bisphosphate binding (GO:0080025) (Fig. 6B; Table S16). All the genes significantly related to ABA, GA₁ and GA₇ were also enriched in the circadian rhythm (plant) pathway (ko04712) (Fig. 6C; Table S17). In the circadian rhythm (plant) pathway, the bZIP transcription factor HY5 (unigene-21666.126192 and unigene-21666.239597) and protein FLOWERING LOCUS T (FT) (unigene-54902.0) were up-regulated. Moreover, the expression of genes significantly related to ABA, GA₁ and GA₇ showed that unigene-21666.90795 (phenylalanine/tyrosine ammonia-lyase, PTAL) was relatively highly expressed in HLZ (Fig. 6D). The gene of PTAL is a crucial gene in phenylalanine metabolism and phenylpropanoid biosynthesis pathway which might conduct bamboo culm color variation.

Differentially expressed transcription factors

Differentially expressed transcription factors were filtered out of the annotated unigenes. A total of 21 differentially expressed transcription factors were identified (Fig. 7; Table S18). The genes unigene-21666.57044 (bHLH), unigene-21666.218036 (MYB-related), and unigene-21666.41135 (SBP) were not detected in the HLZ samples and were only detected in the LZ samples. bHLH transcription factors can interact with MYB transcription factors to regulate anthocyanin biosynthesis (Hichri et al., 2010). However, unigene-21666.194872 (FAR1), unigene-21666.98709 (MADS-M-type), unigene-21666.139384 (MYB), and unigene-21666.212035 (TRAF) were not detected in the LZ samples and were only detected in the HLZ samples. MADS-box transcription factors are involved in flower promotion and development, and simultaneous death usually follows after mass production of bamboo flowers (Abe, Miguchi & Nakashizuka, 2001). Four genes belonged to the bZIP transcription factor family, among which three genes (unigene-21666.126192, unigene-21666.239597, and unigene-43769.0) were annotated as the transcription factor HY5. The target genes of HY5 participate in many biological signaling processes, such as light signaling, circadian clock, anthocyanin biosynthesis, and chlorophyll biosynthesis (Gangappa & Botto, 2016). HY5 could induce the expression of the structural genes CHS (chalcone synthase), CHI (chalcone isomerase), and FLS (flavonol synthase) to regulate anthocyanin biosynthesis (Gangappa & Botto, 2016). There were three genes (unigene-21666.139348, unigene-21666.15994, and unigene-61591.0) that were members of the MYB transcription factor family. MYB transcription factors could regulate anthocyanin biosynthesis (Wang et al., 2019) and combine with other transcription factors to form MYB-bHLH-WDR complexes to regulate flavonoid biosynthesis (Xu, Dubos & Lepiniec, 2015).