Variability in Leaf Color Induced by Chlorophyll Deficiency: Transcriptional Changes in Bamboo Leaves

1. Introduction

The diversity of leaf characteristics, including leaf color, is indeed a significant area of research in plant science [1,2]. The synthesis and functionality of chlorophyll are key factors determining leaf color [3]. During leaf development, the synthesis and degradation of chlorophyll are regulated by multiple genes [4]. In some plants, even though the normal chlorophyll synthesis genes may be present at the genetic level, abnormal gene expression regulation can lead to aberrant chlorophyll synthesis or functionality, ultimately affecting leaf color [5]. The synthesis and functionality of chlorophyll involve complex gene interactions and regulatory networks [6]. Disruptions in these interactions and networks by genetic factors can result in abnormal chlorophyll synthesis or functionality, consequently impacting plant leaf color [1,7].

In this study, we focus on Bambusa multiplex f. silverstripe (hereafter referred to as Silverstripe) a natural variant of Bambusa multiplex widely distributed in South China. Silverstripe displays a variety of leaf colors, including green, white, and green-white. Notably, phenotypic variation due to genetic information may lead to changes in gene expression caused by antagonistic interactions. Transcriptome studies can provide insights into the mechanisms underlying the chlorophyll deficiency and identify differentially expressed genes (DEGs) among different tissues [8]. Therefore, Silverstripe provides us with a unique opportunity to delve into the molecular mechanisms behind chlorophyll anomalies and genetic factors. By applying RNA-seq, we contribute new insights to unravel the mysteries behind plant leaf color diversity, laying the foundation for future research and applications. Exploring this field will enhance our understanding of plant development and leaf formation, ultimately providing valuable support for the development of plant breeding.

2. Methods

2.1. Plant Materials and Transcriptome Sequencing

The samples were collected from the Bamboo Botanical Garden of Fujian Agriculture and Forestry University, located in Cangshan District, Fuzhou City, Fujian Province, China (N26°05', E119°14'). For RNA extraction, fresh leaves were collected from five Silverstripe plants under the same normal maintenance conditions and free from any diseases, including individuals with green leaves (G), white leaves (W), and green-white leaves (GW). We ensured that the selected leaves from each plant exhibited similar growth levels, even though their sizes varied (Figure 1A). The collected leaves were promptly mixed and then flash-frozen in liquid nitrogen. Subsequently, they were then stored at -80°C. Total RNA from the leaves was extracted using the RNA prep Pure Plant Kit (Tiangen, Beijing, China). The quality of total RNA was assessed by 2% agarose gel electrophoresis. RNA concentration was measured using the NanoPhotometer^® spectrophotometer (IMPLEN, CA, USA) and the Qubit^® RNA Assay Kit with a Qubit^® 2.0 fluorometer (Life Technologies, California, USA). RNA integrity was evaluated using the RNA Nano 6000 Assay Kit on the Agilent^® Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Library construction was carried out using the NEB-Next^® Ultra^™ RNA Library Prep Kit for Illumina^® (NEB, USA). All samples were sequenced on the Illumina^® 6000 platform, generating 150 bp paired-end reads. RNA-seq for each type of leaf was performed with three technical replicates.

2.2. Assembly and Function Annotation

The initial processing of raw RNA-Seq reads was performed using fastp v0.23.2 [9] to remove sequences with adapters and those where the percentage of bases with a Qphred ≤ 20 exceeded 50% of the total length, as well as sequences with an N-base percentage > 15%. The remaining clean reads were utilized for reference-based transcriptome assembly. De novo transcriptome assembly was carried out using Trinity v2.5.1 with default parameters [10]. Subsequently, Corset v1.07 [11] was employed for clustering to obtain non-redundant transcript groups. The longest transcript from each transcript group was chosen as a unigene. Custom Python scripts were utilized to rename each unigene for subsequent analyses. The completeness of the transcriptome was assessed using the Benchmarking Universal Single Copy Orthologs (BUSCO) tool v5.1.2 [12].

For functional annotation, the transcripts were annotated using eggNOG v2.1.12 [13], based on diamond alignment (e-value ≤ 0.001), and InterProScan v5.63 with default parameters. Transcripts that remained unannotated were subjected to a BLAST v2.10.0 search (e-value ≤ 1e-5) [14] against the Arabidopsis thaliana (TAIR10) coding sequences (CDS) database. Additionally, transcripts were mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database using KofamKOALA [15] and the KAAS [16] server with default parameters.

2.3. Alignment, Quantification, and Identification of DEGs

In order to assess the utility of the reference transcriptome for gene expression quantification, we utilized the align_and_estimate_abundance.pl script from Trinity. This script utilizes bowtie2 v2.4.2 [17] for transcriptome alignment and RSEM v1.3.3 [18] for transcript abundance estimation. TBtools v2.012 [19] calculates the transcripts per kilobase million (TPM) normalized gene expression estimations, TPM values were used for the Principal Component Analysis (PCA) analysis, while log2 transformed TPM values were used for gene expression graphs.

For statistical differential expression tests, count data were further analyzed using the R package DESeq2 v1.41.12 [20], |log2 (Fold Change)| ≥1 and FDR (False Discovery Rate) < 0.01 were used as screening criteria. Subsequently, common DEGs among different types were identified using the Venn module within TBtools.

2.4. Expression Clustering and Functional Enrichment

For the clustering of common DEGs, TPM values were utilized, and the analysis was performed using Mfuzz v2.60.0. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was conducted using ClusterProfile v4.10.0 [21]. To create gene expression graphs, the heatmap module within TBtools was employed.

3. Results

3.1. Summary of Sequencing Data and Assembly

After trimming and quality assessment of the raw reads, we obtained a total of 76.17 GB of clean reads (Table S1). Following the assembly process, we identified 336,879 non-redundant unigenes. Subsequently, the unigenes were renamed from Bm000001 to Bm336879. Among them, 100951 (29.97%) of the unigenes ranged from 201 to 500 bp in length, 111452 (33.08%) of the unigenes ranged from 501 to 1000 bp in length, 77117 (22.89%) of the unigenes ranged from 1001 to 2000 bp in length, and 47359 (14.06%) unigenes were > 2000 bp in length (Table S2). The N50 value of the assembly was 1,586 bp, and the Benchmarking Universal Single-Copy Orthologs (BUSCO) value was 89.0% (Figure S1).

3.2. Expression Quantification and Classification

To elucidate the transcriptional mechanisms underlying the heterogeneity of leaf color in bamboo, we conducted an RNA-seq experiment covering three leaf colors (Figure 1A). By aligning and quantifying the sequencing data, we obtained count data and calculated TPM values (Table S3; Table S4). To further investigate the differences between the chlorophyll-deficient type and the normal type, we performed PCA analysis on the expression data of the nine transcriptome samples and generated a correlation heatmap to study their variation and similarity. In the scatterplot between PC1 and PC2, it is evident that the three replicates of each leaf color type cluster closely together, indicating similarities in gene expression levels. Additionally, there is a notable distinction between G and GW or W, whereas the GW and W exhibit closer distribution and smaller differences (Figure 1B).

Furthermore, the correlation heatmap demonstrates high within-group correlation for samples of the same leaf color type, while the inter-group correlation is comparatively lower (Figure 1C). In line with the PCA analysis results, the correlation between W and GW is significantly higher than the correlation between G and either of them.

3.3. Identification of DEGs

Significant differences in transcript abundance were identified by DESeq2 based on the count data (Figure S2). Among the different combinations, the G vs. W combination had the highest number of DEGs. There were a total of 22,523 upregulated genes and 28,978 downregulated genes in this combination. In contrast, the GW vs. W combination had the lowest number of DEGs, with 2,846 genes upregulated and 1,326 genes downregulated. Additionally, in the G vs. GW combination, there were 12,220 upregulated genes and 23,069 downregulated genes (Figure S1). These results suggest that as the differences in leaf color become more pronounced in Silverstripe, the number of DEGs increases.

Further analysis revealed that there were 1,430 common DEGs shared among the three combinations (Figure 1D). Moreover, cluster analysis showed that these DEGs could be categorized into three clusters (Figure 1E). Furthermore, there were 405 genes and 609 genes that exhibited the highest expression levels in the W and G groups, respectively, and their expression patterns consistently upregulated or downregulated with the degree of leaf whitening (Figure 1E; Figure 1G). This implies that these genes may be associated with the proportion of leaf white area and are considered as key genes related to chlorophyll deficiency.

3.4. Function Annotation and Enrichment

To gain a deeper understanding of the molecular mechanisms behind physiological changes and to explore the biological significance of key genes in more detail, we initiated functional annotation for all unigenes. A total of 170,574 unigenes were successfully associated with GO terms (Table S5), accounting for 50.63% of all unigenes. However, only 53,046 unigenes, representing 15.75% of the total unigenes, could be effectively annotated to corresponding KEGG pathways. These pathways were distributed across 121 different categories (Table S6).

To gain a more in-depth understanding of the functional roles of the key genes in leaf color differentiation, we conducted KEGG pathway and GO enrichment analyses. In W, the relative upregulated key genes were subjected to GO functional enrichment, revealing their enrichment in 145 GO terms, which could be categorized into five main classes: amino acid and organic acid transport (e.g., amino acid transport, amino acid transmembrane transport), cell wall metabolism (e.g., cell wall macromolecule catabolic process, cell wall macromolecule metabolic process, chitinase activity), wall-related enzyme activity (e.g., chitinase activity and galactosyltransferase activity), metabolic pathways (L-phenylalanine biosynthetic process, chitin metabolic process), and the regulation of biological processes (e.g., RNA-directed RNA polymerase complex, regulation of DNA-templated transcription) (Figure 2A; Table S7). Additionally, KEGG analysis only enriched eight pathways. This might be due to the relatively smaller number of unigenes included in the KEGG annotation background. However, similar to the GO enrichment results, these pathways can be simplified into categories such as transcription factors (TFs) regulatory pathways (e.g., Transcription factors and Plant hormone signal transduction), organic acid metabolism (e.g., Phenylalanine, tyrosine, and tryptophan biosynthesis and alpha-Linolenic acid metabolism), and Starch and sucrose metabolism (Figure 2A; Table S7).

Among the downregulated key genes, GO functional enrichment analysis revealed their enrichment in 175 GO terms, primarily falling into four categories: photosynthesis-related pathways (e.g., photosystem, photosynthetic membrane), chloroplast-related processes (e.g., chloroplast thylakoid, plastid stroma), redox reactions (e.g., oxidoreductase complex, photosystem II oxygen evolving complex, oxidoreductase activity), and metabolic processes (e.g., glycine catabolic process, serine family amino acid catabolic process) (Figure 2B; Table S8). Correspondingly, KEGG pathway analysis enriched 17 pathways, largely related to photosynthesis (e.g., Photosynthesis - antenna proteins, Photosynthesis proteins), amino acid metabolism (e.g., Alanine, aspartate, and glutamate metabolism, Arginine biosynthesis), and carbon and energy metabolism (e.g., Glyoxylate and dicarboxylate metabolism, One carbon pool by folate) (Figure 2B; Table S8).

The result of GO and KEGG functional enrichment indicates that as chlorophyll content decreases in the leaves, genes involved in TFs regulatory and organic acid metabolism pathways are significantly upregulated, while genes related to photosynthesis and carbohydrate metabolism pathways may be downregulated as the extent of chlorophyll deficiency increases.

4. Discussion

In previous studies on variegated plants, the rich diversity in leaf colors is primarily attributed to variations in the types, quantities, and distribution of pigments within leaf cells [22]. However, Silverstripe as a variegated plant, it characterized by a fact that white leaves resulting from chlorophyll deficiency, which represents a significant departure from color changes induced by pigment accumulation. Therefore, studying the functional genes associated with chlorophyll deficiency phenotype provides novel insights into plant leaf color variation.

Upon the identification of DEGs among different color types, it was observed that the number of downregulated genes in the chlorophyll-deficient W and GW phenotypes, as compared to the green G phenotype, exceeded the number of upregulated genes (Figure S1). Chlorophyll is a fundamental component of photosynthesis, and its deficiency can impact energy conversion and biomass accumulation [23]. Consequently, chlorophyll deficiency likely mediates gene expression changes at the RNA level, resulting in disruptions in various biological functions, including photosynthesis. The functional enrichment analysis of key genes that exhibit consistent upregulation in the chlorophyll-deficient phenotype reveals their association with transcription factor regulation and organic acid metabolism (Figure 2A). Plants can adapt to chlorophyll deficiency by adjusting metabolic pathways such as carbon and organic acid metabolism [24]. Some intermediates of the citric acid cycle are known to regulate certain enzymes in the chlorophyll synthesis pathway [25]. Notably, one of the precursors required for chlorophyll synthesis is 5-aminolevulinic acid (ALA) [26], and the citric acid cycle, also known as the tricarboxylic acid cycle, is one of the pathways for ALA production [25]. This suggests that the regulation of organic acid metabolism may influence the rate of chlorophyll synthesis. In the following discussion, we will explore the ALA biosynthesis pathway, the transcriptional changes related to chlorophyll deficiency, and their impact on photosynthesis and sugar metabolism.

4.1. Porphyrin Biosynthesis in Chlorophyll Deficiency in Bamboo Leaves

Chlorophyll is the primary pigment in plant chloroplasts, actively participating in photosynthesis to convert light energy into chemical energy. Chlorophyll synthesis is a complex process involving multiple steps and precursor molecules, such as ALA, porphyrins, magnesium protoporphyrin, protopheophorbide, and protochlorophyllide [27]. Additionally, pyrrole is a nitrogen-containing organic compound and a crucial precursor for pigments like chlorophyll and hemoglobin [28]. The pyrrole metabolic pathway refers to a series of biochemical reactions from the precursor substance ALA to chlorophyll or hemoglobin [29].

Unexpectedly, many structural genes related to the pyrrole biosynthetic pathway are instead increased expression in tissues with a higher degree of chlorophyll deficiency (Figure 2C). We identified 28 key genes in this pathway, with the majority of genes upregulated in the white phenotype, and only six genes downregulated. One possible inference is that these genes may be involved in negatively regulating the synthesis of various enzymes in the pyrrole biosynthetic pathway or accelerating the degradation of some chlorophyll precursor molecules. Additionally, when chloroplasts are deficient or damaged, plants may increase the synthesis of other pigments to compensate for the loss of photosynthesis [30]. Furthermore, the pyrrole metabolic pathway also provides essential functional molecules for plants, such as hemoglobin, vitamin B12, and carotenoids, which play crucial roles in plant respiration, metabolism, antioxidant defense, and photoprotection [28]. Therefore, the upregulation of these genes in expression may contribute to maintaining the supply of chlorophyll precursors, thus sustaining a certain level of chloroplast function and physiological stability. However, in the ALA pathway, a small subset of genes also exhibits higher expression in the green phenotype (Figure 2C). Since light is essential for the expression of hundreds of chloroplast- and photosynthesis-related genes, chlorophyll deficiency may result in some light-responsive genes not being activated normally [31]. This could potentially lead to the negative regulation of a few ALA synthesis genes by light-responsive genes.

4.2. Chlorophyll Deficiency Induces Alterations in the Photosynthetic Process

The expression of pigment synthesis protein genes is followed by the expression of genes related to the light-dependent reactions of photosynthesis and the final carbon metabolism [32]. Chlorophyll is a core component of photosynthesis, and its loss affects energy conversion and biomass accumulation [23]. The results indicate that the downregulated DEGs in white leaves are involved in the light reactions of photosynthesis (Figure 2D), as expected due to the lack of chlorophyll in white leaves, hindering the proper functioning of photosynthesis.

The genes predominantly downregulated in white leaves are those that play a role in the light reactions of photosystem I and the reduction phase of the Calvin cycle (Figure 2D; Figure 2E), such as NADP+ reductase and glyceraldehyde-3-phosphate dehydrogenase. These genes are responsible for converting light energy into chemical energy and reducing fixed carbon dioxide to glyceraldehyde-3-phosphate. The downregulation of these genes implies a reduction in the efficiency of light energy utilization and carbon fixation in white leaves, leading to a decrease in carbohydrate levels. This aligns with reported effects of chlorophyll deficiency on plant physiological metabolism [33,34].

Additionally, this may reflect the inhibitory response of leaves to genes related to photosynthesis under conditions of chlorophyl deficient. This regulatory mechanism may involve complex signaling pathways, including light signal transduction and plant hormone regulation [6]. These findings remind us that the differential gene regulation in white leaves is associated with maintaining plant energy balance and growth, providing deeper insights into understanding the regulatory network of photosynthesis.

4.3. Differential Glycolysis in Chlorophyll-Deficient Phenotypes

An intriguing finding is the downregulation of key genes in white leaves associated with glycolysis. This suggests that, despite defects in photosynthesis, white leaves may sustain survival by adjusting glycolysis. A previous study emphasized that plants, when facing disruptions in photosynthesis, regulate the activity of sugar metabolism pathways to maintain energy balance [35]. Glycolysis is crucial for plant growth and development as it not only provides energy but also participates in various biochemical processes [36]. A part of the glycolytic pathway, known as the Pentose Phosphate Pathway (PPP), shares some common intermediates with the Calvin cycle (Figure 2F). These intermediates can be interconverted in different pathways to meet the cell's energy and synthesis requirements [37,38].

On the other hand, deficiencies in protochlorophyllide or any of these other components lead to absent or aberrant prolamellar bodies [39]. Photosynthesis in chloroplasts is a highly regulated process, and the downregulation of sugar metabolism-related genes in white leaves may be due to the potential generation of oxygen radicals when chlorophyll is absent, leading to oxidative stress [6]. This oxidative stress results in lipid peroxidation and cell membrane damage, thereby reducing sugar metabolism to minimize reactive oxygen species production. We observed that the downregulated genes in white leaves are primarily involved in the glycolytic pathway, such as phosphoglyceraldehyde kinase and phosphoglycerate mutase, which both participate in the transfer of phosphate groups in sugar metabolism. The downregulation of these genes implies lower glycolytic activity and pentose phosphate pathway activity in white leaves, possibly as a response to cell damage and energy consumption caused by oxidative stress. Additionally, the regulation of aging by carbon signaling and the accumulation of sugars in aging leaves serve as signals for the initiation or acceleration of plant senescence [35,40]. The downregulation of sugar metabolism genes may contribute to increased sugar accumulation to some extent, ultimately leading to premature senescence in chloroplast-deficient tissues.

5. Conclusion

The diversity of leaf characteristics, including leaf color, is a prominent subject of study in plant science. Leaf color is predominantly determined by the synthesis and functionality of chlorophyll, a key component of photosynthesis. The regulation of chlorophyll synthesis and degradation involves complex gene interactions, and disruptions in these processes can lead to abnormal chlorophyll synthesis and impact leaf color. This study focuses on Bambusa multiplex f. silverstripe, a natural variant with various leaf colors, including green, white, and green-white leaves. The variations in leaf color are attributed to genetic factors and their influence on gene expression. By employing RNA-seq, we investigate the molecular mechanisms behind chlorophyll anomalies and genetic factors in Silverstripe.

Overall, the white phenotype exhibits lower levels of photosynthetic and sugar metabolism activity compared to the green phenotype, with potential implications for cellular structural damage and reduced chlorophyll synthesis efficiency. However, these mutant tissues demonstrate slower growth rates and shorter lifespans, potentially compromising their survival advantages [1,34]. Fortunately, under relatively minor environmental pressures, these mutant variants can maintain stable populations. As horticultural industries evolve, there is a growing appreciation for the ornamental value of these color mutants, leading to their cultivation and selective breeding. This may contribute to positive selection for genes associated with chlorophyll synthesis defects. Given the limited research on the regulatory relationships between mutation sites and transcriptional control in color mutant individuals, further studies from a whole-genome perspective are warranted to explore mutation sites related to chlorophyll deficiency in bamboo.