Transcriptome Reveals the Specificity of Phyllostachys edulis ‘Pachyloen’ Shoots at Different Developmental Stages

Phyllostachys edulis ‘Pachyloen’ can have a stalk wall thickness of up to 2.5 cm at a height of 1.3 m, which is 1.8 times that of normal Moso bamboo (Phyllostachys edulis); this serves as an excellent cultivar, comprising both wood and bamboo shoots. We collected bamboo shoot samples of Phyllostachys edulis ‘Pachyloen’ and Moso bamboo on a monthly basis from September to April and used transcriptome sequencing to explore the differences in their development. The results showed that there were 666–1839 Phyllostachys edulis ‘Pachyloen’-specific genes at different developmental stages enriched in 20 biological processes, 15 cellular components, 12 molecular functions, and 137 metabolic pathways, 52 of which were significant. Among these, 27 metabolic pathways such as tyrosine metabolism and their uniquely expressed genes were found to play important roles in the thickening of Phyllostachys edulis ‘Pachyloen’. This study provides insights into the mechanisms underlying the thickening of the culm wall of Phyllostachys edulis ‘Pachyloen’.


Introduction
Bamboo plants are non-timber forest products of important economic, ecological, and cultural value [1][2][3][4][5]. Approximately 2.5 billion people worldwide rely economically on bamboo, and the annual trade generated by bamboo amounts to 2.5 billion dollars [6,7]. Bamboo plants are different from trees in that bamboo stems are hollow and there is no periodic radial growth process. Phyllostachys edulis 'Pachyloen' is a novel variety of Phyllostachys edulis (normal Moso bamboo) that exhibits several advantageous features of timber and bamboo shoot including bamboo wall thickening, subsolid culm, large biomass, good tasting shoot, and stable heredity [8]. It was first discovered in Jiangxi Province in southeastern China in 1995. The thickness of its culm wall is almost twice as that of the normal phenotype of the species, and its yields are comparatively higher [9]. The major biological characteristics of P. edulis 'Pachyloen' including the growth pattern of its bamboo shoots, photosynthetic physiology, transpiration dynamics, cold resistance, anatomical characteristics, and physical and mechanical properties have been studied systematically since its discovery [10][11][12].
The thickness of bamboo stalks is closely related to the development of bamboo shoots [13]. The walls of bamboo stalks are primarily composed of thin-walled basic tissue cells and vascular bundles, with thick-walled cells forming the main body. Vascular bundle-related genes regulate the growth and development of vascular bundles and improve plant resistance [14]. The expression of GbaVd1 and GbaVd2 genes has been associated with the increase in the degree of lignification of transgenic Arabidopsis lines [15]. Bamboo research has entered the era of functional genomics with the sequencing of the Moso bamboo genome [16,17]. In addition, the cellular mechanisms underlying the growth of the underground shoots of Moso bamboo have been systematically studied [18][19][20][21]. The molecular mechanisms underlying the regulation of bamboo flowering have also been elucidated using RNA-Seq [22][23][24][25]. However, the mechanism underlying the apparent thickening of the culm walls of P. edulis 'Pachyloen' remains unclear.
In the present study, the shoots of P. edulis 'Pachyloen' and Moso bamboo (control) in different growth stages were analyzed by transcriptome sequencing using the whole genome sequence of Moso bamboo as a reference. As a result, several uniquely expressed genes and their related pathways were discovered in P. edulis 'Pachyloen'. The results presented in this study have important scientific and theoretical implications for the characterization of the mechanism underlying the thickening of bamboo culm walls.

Plant Materials
Six to fifteen pieces of P. edulis 'Pachyloen' shoots were collected between September 2019 and April 2020 on a monthly basis from the experimental bamboo forest (113.1124063 • E, 28.2698183 • N, 44.9 m above sea level) in Lukou Town, Changsha County, Hunan Province, China, based on the speed of bamboo shoot growth. Moso bamboo shoots were also collected for use as a control at the same time and under the same environmental conditions. The samples were collected at 1 cm from the apex of the shoot. The volume and degree of development of the experimental shoots were consistent across the samples. The collected shoots were immediately washed with distilled water, wrapped in foil, placed in a liquid nitrogen tank, and stored at −80 • C until RNA extraction. The samples of P. edulis 'Pachyloen' were marked as H9, H10, H11, H12, H1, H2, H3, and H4; similarly, the Moso bamboo samples were marked as M9, M10, M11, M12, M1, M2, M3, and M4. The sample collection time and dates are outlined in Table 1. Three biological repeats, each with three technical repeats, were performed for all samples. Tender shoots lacking sheaths (50-100 mg) were ground in liquid nitrogen. Total RNA was extracted using the TRIzol reagent (Invitrogen Scientific, Inc., Carlsbad, CA, USA) in accordance with the manufacturer's protocol. The quality of the total RNA of the sample was detected using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA) with an RNA integrity number (RIN) value of ≥7 required for all samples. The libraries were constructed using a VAHTS mRNA-Seq V2 Library Prep Kit for Illumina ® (Vazyme Biotech Co., Ltd., Nanjing, Jiangsu, China) in accordance with the manufacturer's instructions. The libraries were then sequenced on an Illumina sequencing platform (HiSeq™ 2500). The Illumina raw sequencing data were submitted to the National Center for Biotechnology Information (NCBI) Short Reads Archive (SRA) database under accession number SRP169499.

Sequence Data Analysis and Assembly
The original image data were transformed into sequence data by base calling, which is defined as raw reads. Reads shorter than 50 bp, with a predicted error rate of >20%, or with full passes lower than 0 were removed [26]. The clean reads were aligned with the Moso bamboo reference genome. The data sources of the Moso bamboo reference genome are listed in Supplementary File S1. Coding sequences (CDS) of Moso bamboo with a mismatch of up to five bases were allowed to obtain the gene expression profile data for 'Pachyloen' at the different developmental periods.

Differential Expression Analysis
The differentially expressed genes (DEGs) (old change ≥ 2, FDR ≤ 0.001) detected at different developmental stages of 'Pachyloen' were subjected to Gene Ontology (GO) functional annotation and classification as well as to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway significance enrichment analysis. Hierarchical cluster analysis was used to annotate and analyze the GO results and explore patterns of gene expression in WEGO software [27]. R software (v3.6.1) was used to analyze the gene expression patterns of specific metabolic pathways.

Quantitative Real-Time PCR (qRT-PCR) Validation of Differential Expression
Ten putative genes (PH01000081G111, PH01000104G133, PH01000152G124, PH01000159G013, PH01000931G0650, PH01001141G0100, PH01001292G0400, PH01001299G0140, PH01001317G0160, and PH01001859G0010) were randomly selected for validation by quantitative real-time PCR (qRT-PCR) to validate the accuracy of the RNA-Seq results. The transcription shoots were analyzed by qRT-PCR in an ABI 7300 Real-Time PCR system (Applied Biosystems, Foster City, CA, USA) using SYBR Green Realtime PCR Master Mix. The tonoplast intrinsic protein (TIP41) gene was amplified using the forward primer 5 -AAAATCATTGTAGGCCATTGTCG-3 and the reverse primer 5 -ACTAAATTAAGCCAGCGGGAGTG-3 as an internal control [28]. One microgram of extracted RNA was reverse-transcribed using a Superscript II First Strand Synthesis Kit (Invitrogen) and random hexamer primers according to the manufacturer's instructions. For each sample, three replicate experiments were performed at a final volume of 10 µL (3 µL of ddH 2 O, 1 µL of cDNA, 5 µL of SYBR Mix, 0.5 µL of primers, and 0.5 µL of ROX Reference Dye). The cycling conditions were as follows: 95 • C for 10 min, 95 • C for 10 s, and 62 • C for 10 s for 50 cycles. The qRT-PCR primers used are listed in Supplementary File S2. The amplification results were analyzed using the comparative cycle threshold (Ct) method with the 2 −∆∆CT formula. The qRT-PCR results were calculated based on the mean values from three replicate experiments. Significant differences between treatments were evaluated based on their standard deviation. The error bars represent the standard error of the three biological replicates.

Read Mapping to the Reference Genome Dataset
The clean reads of H9-H4 and M9-M4 were mapped to the genes expressed in the Moso bamboo reference genome (Table 2). Approximately 27,000 to 28,000 genes were annotated to the reference gene in each sample, accounting for 88% of the genes in the Moso bamboo reference genome, which proved that the sequencing library constructed by cDNA fragments was comprised of high quality and reliable sequencing data. During the eight months of growth and development of P. edulis 'Pachyloen', the number of expressed genes was the highest in February at 28,911, and the lowest in April at 27,265. Similarly, for the Moso bamboo samples, the number of expressed genes in March was the highest at 28,448, and the lowest in April at 26,958. The gene expression patterns of P. edulis 'Pachyloen' and Moso bamboo based on the gene coverage and randomness of reads mapped to the reference genome and genes are outlined in Supplementary File S3.

Trend of Alterations in Gene Expression Quantity
January and February are typically the coldest months of the year. Information on the climate of the collection sites during sampling is provided in Supplementary File S4. In winter, bamboo shoots enter a dormant stage, known as the overwintering period, to resist the harsh conditions in the outside world. Bamboo shoots then enter the germination stage in March, during which the number of genes being expressed increases, reflecting an increase in fertility. The gene expression number of P. edulis 'Pachyloen' and Moso bamboo (control) increased at first, and then decreased with changes in phenology ( Figure 1). The number of genes expressed in P. edulis 'Pachyloen' each month was higher than Moso bamboo, except for March, which may be related to the physiological and metabolic activities of P. edulis 'Pachyloen'.
The highest number of shared genes was 27,549, in the samples collected in January, and the lowest was 25,735 in the samples collected in April. Several P. edulis 'Pachyloen'-specific genes were detected upon comparison with the genes from Moso bamboo as control genes. Information on the P. edulis 'Pachyloen'-specific genes and Moso-specific genes is provided in Supplementary File S5. P. edulis 'Pachyloen' had more unique genes than Moso bamboo throughout its growth and development, except in March. The identification of P. edulis 'Pachyloen'-specific genes and its functional annotation and metabolic pathway enrichment analysis play an important role in explaining the thickening mechanism of P. edulis 'Pachyloen'.  The unigenes were mainly enriched in 22 terms of biological processes (BP), 16 cellular components (CC), and 13 molecular functions (MF), according to the GO enrichment analysis of the DEGs (Figure 3). "Cellular process", "metabolic process", "cell", "cell part", "membrane", "binding", and "catalytic activity" were rich in unigenes. The seven GO terms with the most enrichment of P. edulis 'Pachyloen'-specific genes at the different developmental stages were counted for the comparison of the developmental differences between P. edulis 'Pachyloen' and Moso bamboo (Table  3).  (Figure 3). "Cellular process", "metabolic process", "cell", "cell part", "membrane", "binding", and "catalytic activity" were rich in unigenes. The seven GO terms with the most enrichment of P. edulis 'Pachyloen'-specific genes at the different developmental stages were counted for the comparison of the developmental differences between P. edulis 'Pachyloen' and Moso bamboo (Table 3).    Table 4). The   11,033 (H3/M3), and 2832 (H4/M4) DEGs were subjected to pathway enrichment analysis using the KEGG pathway database. The co-expressed genes and Moso bamboo-specific genes were enriched in 68 and 35 significant metabolic pathways (p ≤ 0.05). Flavone and flavonol biosynthesis, phenylpropanoid biosynthesis, and phenylalanine metabolism were the top three of the number of DEGs annotated to each metabolic pathway as a percentage of the total number of genes in the metabolic pathway. The P. edulis 'Pachyloen'-specific genes were enriched into 52 significantly enriched metabolic pathways throughout the development process (H9-H4; p ≤ 0.05) ( Table 4). The metabolic pathways of the co-expressed genes, the P. edulis 'Pachyloen'-specific genes, and the Moso bamboo-specific genes were analyzed to determine the key metabolic pathways involved in the development of bamboo thick-walled germplasm. As a result, 27 significant metabolic pathways were considered to be related to the thickening of P. edulis 'Pachyloen' (Table 5), 15 metabolic pathways were only involved in the expression of P. edulis 'Pachyloen'-specific genes, and 12 metabolic pathways were involved in DEGs and P. edulis 'Pachyloen'-specific genes.      The free amino acid content and biomass of P. edulis 'Pachyloen' were significantly greater than that of the control (Moso bamboo). Therefore, we analyzed the pathways of the tyrosine metabolism, sulfur metabolism, and biosynthesis of unsaturated fatty acids from 27 significant metabolic pathways. Fifteen P. edulis 'Pachyloen'-specific genes were found to be involved in the regulation of wall thickening ( Table 6). Our findings on the other 24 metabolic pathways are provided in the Supplementary Materials. Table 6. P. edulis 'Pachyloen'-specific genes in three significant pathways.

Quantitative Real-Time PCR (qRT-PCR) Validation of Differential Expression
The expression profiles of 10 randomly chosen genes are provided in Supplementary File S6. The relative gene expression was calculated using 2 − Ct and amplified using a LightCycler ® 96 fluorescence quantitative PCR instrument for three biological replicates. The qRT-PCR verification of the 10 selected genes was consistent with the sequencing results.

Molecular Characteristics of Development in P. edulis 'Pachyloen'
In terms of morphology, the meristem at the apex of the underground stem of P. edulis 'Pachyloen' does not differentiate into medullary tissue; rather, it differentiates into vascular bundle tissues. The apical meristems of P. edulis 'Pachyloen' shoots are not similar to the low and wide ridge-shaped meristems of Moso bamboos; instead, they are tall and narrow finger-shaped meristems. Moreover, the meristems of thick-walled bamboo shoots are larger, with a greater number of cells. Larger apical meristems can produce more auxin, which acts as a promoter of the differentiation and proliferation of vascular tissue [29,30]. The same phenomenon occurs in sunflower and Arabidopsis mutants [31,32]. Certain studies have reported that wheat stalk thickness is controlled by a single gene, complementary genes, or multiple genes. The rice cultivar "Kasalath" has chromosomal fragments that increase stem resistance to enhance lodging resistance [33]. Yano et al. detected a quantitative trait locus (QTL) named SCM3, which regulates the thickness of rice stem and stem wall [34]. However, the molecular mechanism underlying bamboo culm wall thickening is yet to be completely elucidated. In this study, we detected 27 key metabolic pathways associated with P. edulis 'Pachyloen' culm wall thickening, including tyrosine metabolism, sulfur metabolism, and unsaturated fatty acid biosynthesis. Plant endogenous hormones such as cyclin A, cyclin B, auxin, and gibberellin appeared to strongly influence cell division, cell elongation, and cell cycle. The P. edulis 'Pachyloen' genes identified play important regulatory roles in different periods.

Gene Expression during the Development of P. edulis 'Pachyloen'
The increase in the cell number of the culm wall plays an important role in promoting wall thickening in P. edulis 'Pachyloen' [35][36][37]. In this study, we investigated the changes in the morphology of thick-walled P. edulis 'Pachyloen' from shoots to bamboo using transcriptomic analysis to identify the key metabolic pathways involved in the formation of its thick culm walls. In general, the different growth and development states of plants are known to be directly reflected in their DEGs. In this study, the shoot of P. edulis 'Pachyloen' and Moso bamboo showed evident phenological characteristics and their underground growth process could be divided into differentiation periods: the overwintering period and germination period [38,39]. At the end of summer and early autumn, the side buds on the bamboo whip began to germinate and differentiate into shoots. In winter, the bamboo shoots entered the overwintering period and lay dormant. By the second year, the bamboo shoots continued to grow, and increased in size with rising temperatures. During this process, extremely complex physiological and biochemical reactions occurred in bamboo shoots [40,41]. We found that there was a lack of regular changes in the number of upregulated and downregulated DEGs related to the growth and development of P. edulis 'Pachyloen' between developmental stages.
From the winter (November to January) of the following year, the bamboo shoot buds were dormant, with a small number of genes being actively expressed due to low temperatures. The highest number of DEGs were expressed from March to April. Once the shoots had completed the underground development stage, cell proliferation increased and shoots emerged from the ground, marking the beginning of a rapid increase in cell volume. As bamboo shoots undergo the greatest number of changes when undergoing growth above ground, the number of downregulated genes was much higher during March and April than in other months.

Genetic Regulation of Important Moso Bamboo Traits
At present, there are no studies on the omics of bamboo culm wall thickening; however, studies have been conducted on the omics of rapid growth and flowering. Bamboo plants have the unique reproductive biological characteristics of a long flowering cycle and plant death after blooming. It has been reported in certain studies that different key genes play a central role in the regulation of different growth periods of flowers; 290 genes related to the transition from vegetative to reproductive growth and 47 genes related to the growth and development of flowers have been identified. MADS, GID1, GID2, LHY, and MYB were upregulated during flower development. The expression of Dof genes was observed to be significantly inhibited during anthesis and embryo formation, which suggests that Dof genes play a more significant role in the early stages of flower development [42]. Additionally, 6076 upregulated genes and 4613 downregulated genes were identified during the rapid growth stage of Moso bamboo. Concurrently, candidate genes related to transcription factors, plant hormones, cell cycle regulation, cell wall metabolism, and cell morphogenesis have also been identified [13]. There are 31,987 genes in the Moso bamboo reference genome [43]. It is obvious that a large number of single genes regulate multiple traits in Moso bamboo. The gene regulation network of traits should be continuously evaluated until it is completely characterized.

Conclusions
Transcriptome sequencing analysis was performed on the shoots of P. edulis 'Pachyloen' and Moso bamboo from September to April of the following year to determine the developmental differences at the omics level. In all, 666-1839 genes were uniquely expressed in P. edulis 'Pachyloen', while 1079-1356 genes were uniquely expressed in Moso bamboo. The P. edulis 'Pachyloen'-specific genes were enriched in 20 BPs, 15 CCs, 12 MFs, and 137 metabolic pathways including 52 significant metabolic pathways. In total, 27 metabolic pathways such as tyrosine metabolism were considered important in the thickening of P. edulis 'Pachyloen'.
This study provides a solid foundation for the further study of the thickening characteristics of P. edulis 'Pachyloen' as well as research on bamboo plants more generally.