Transcriptome Sequencing and Analysis for Culm Elongation of the World’s Largest Bamboo (Dendrocalamus sinicus)

Dendrocalamus sinicus is the world’s largest bamboo species with strong woody culms, and known for its fast-growing culms. As an economic bamboo species, it was popularized for multi-functional applications including furniture, construction, and industrial paper pulp. To comprehensively elucidate the molecular processes involved in its culm elongation, Illumina paired-end sequencing was conducted. About 65.08 million high-quality reads were produced, and assembled into 81,744 unigenes with an average length of 723 bp. A total of 64,338 (79%) unigenes were annotated for their functions, of which, 56,587 were annotated in the NCBI non-redundant protein database and 35,262 were annotated in the Swiss-Prot database. Also, 42,508 and 21,009 annotated unigenes were allocated to gene ontology (GO) categories and clusters of orthologous groups (COG), respectively. By searching against the Kyoto Encyclopedia of Genes and Genomes Pathway database (KEGG), 33,920 unigenes were assigned to 128 KEGG pathways. Meanwhile, 8,553 simple sequence repeats (SSRs) and 81,534 single-nucleotide polymorphism (SNPs) were identified, respectively. Additionally, 388 transcripts encoding lignin biosynthesis were detected, among which, 27 transcripts encoding Shikimate O-hydroxycinnamoyltransferase (HCT) specifically expressed in D. sinicus when compared to other bamboo species and rice. The phylogenetic relationship between D. sinicus and other plants was analyzed, suggesting functional diversity of HCT unigenes in D. sinicus. We conjectured that HCT might lead to the high lignin content and giant culm. Given that the leaves are not yet formed and culm is covered with sheaths during culm elongation, the existence of photosynthesis of bamboo culm is usually neglected. Surprisedly, 109 transcripts encoding photosynthesis were identified, including photosystem I and II, cytochrome b6/f complex, photosynthetic electron transport and F-type ATPase, and 24 transcripts were characterized as antenna proteins that regarded as the main tool for capturing light of plants, implying stem photosynthesis plays a key role during culm elongation due to the unavailability of its leaf. By real-time quantitative PCR, the expression level of 6 unigenes was detected. The results showed the expression level of all genes accorded with the transcriptome data, which confirm the reliability of the transcriptome data. As we know, this is the first study underline the D. sinicus transcriptome, which will deepen the understanding of the molecular mechanisms of culm development. The results may help variety improvement and resource utilization of bamboos.


Introduction
Currently, due to depleting fossil reserves and increasing emission of greenhouse gases, it is obvious that utilization of renewable feedstock is one necessary step towards a sustainable development of our future [1]. Therefore, the efforts to exploit a variety of potential plant feedstocks for the energy, chemical, and food ingredients have been made particularly from agricultural and forestry biomass resources [2,3,4]. Dendrocalamus sinicus, belonging to Bambusoideae of Gramineae, is the world's largest bamboo species with strong woody stems (maximal diameter 30 cm and maximal height 33 m), which is mainly distributed in the southwest region of China [5]. It is called Da bamboo by local residents, and its culms yield per unit area is 5-8 times as that of Phyllostachys pubscens that is now popularized as a major economic bamboo species in China [6]. In its native area, D. sinicus is economically important as a raw material for furniture, construction, and industrial paper pulp [7]. Because of easy propagation, fast growth, and high productivity, it is considered as one of the most potential renewable nonwoody lignocellulosic feedstocks for bioenergy and biorefinery [8].
The culm growth of D. sinicus is extremely amazing: shooting initiate from early June, and height growth of culm end in late August. By the brief three months, new shoots can reach the height and diameter of adult mother bamboo-near to 30 cm of diameter and 30 m of height. To our knowledge, no any plant in this planet has faster growth rate compared with D. sinicus.
In terms of the bamboo growth, although some molecular researches were reported, such as the protein expression profiles [9], cloning of a set of cDNA sequences [10,11,12,13,14,15,16], Expressed Sequence Tags (EST) [17], RNA-seq [18,19], monoclonal antibody banks [20] and draft genome [21], the molecular mechanism related to culm development is still unclear. Especially, to date, there was almost no molecular data of D. sinicus except SSR marker development [22], which seriously restricted its variety improvement and resource development.
Transcriptomic analyses are extremely efficient methods for identifying differential expression genes at the whole-genome level. In contrast to the traditional fragment analysis techniques, RNAs-seq have some advantages: (1) high efficient and low cost; (2) cataloguing all kinds of transcripts including mRNAs, noncoding RNAs, and small RNAs; (3) investigating the transcriptional structure of genes, splicing patterns, and gene isoforms; (4) studying posttranscriptional modification and mutations; and (5) precisely quantifying gene expression in large-scale at the same time of sequencing. Moreover, RNA-seq is genome-independent and is especially useful for analyzing the transcriptome of a species without complete genome information [23,24].
In this study, the high throughput RNA-Seq was performed to elucidate the molecular processes involved in the rapid culm elongation of D. sinicus during an integrated growing season. Also, based on the transcriptome analysis, a series of results were displayed including novo assembly, gene characterization, gene classification, gene enrichment, SSR and SNP marker identification. Given the unique growth characteristics of D. sinicus, this study will lay foundation for deep understanding of bamboo growth.

Plant materials
Culm samples were collected from natural population of D. sinicus located at Ninger county in Yunnan province of China (23°07 0 *23°09 0 N, 101°04 0 *101°08 0 E), where the D. sinicus distributed typically. The sampling site belongs to public forest, and this field sampling was permitted by Ninger County Government. Based on the depiction by Banik (2015) [25], the internodal elongation begins at the basal portion of the culm and then gradually proceeds to the top, that is, elongation is mainly due to the intercalary meristem present at the node. Development, maturation, and aging gradually completed from basal to top internodes. Thus, the basal internode represents the higher lignified degree compared with middle and top internodes. In intercalary growth, the immature axis increases in length by the elongation of cells in zones of secondary meristems each located just above the node. Thus, different internode represents different developmental condition. Therefore, basal culms from different internode originated from a same bamboo were harvested. A total of nine internodes were selected and triplicate replicate was mixed. Samples were washed with deionized water, and wiped with filter paper, then, immediately frozen in liquid nitrogen and stored at −80°C until analysis.

RNA extraction, library construction and RNA-seq
Total RNA of each sample was isolated separately using Trizol Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. The purified RNA concentration was quantified by a spectrophotometer (UV-Vis Spectrophotometer, Quawell Q5000, San Jose, CA, USA), and the purity and degradation of total RNA were checked on 1% agarose gels before proceeding. For maximizing the diversity of transcriptional units, RNA from each sample was mixed into a single pool. The mRNA-seq library was constructed using Illumina's TruSeq RNA Sample Preparation Kit (Illumina Inc, San Diego, CA, USA). Shortly, mRNA was purified by oligo (dT) magnetic beads. After purification, the mRNA was fragmented into small pieces using divalent cations under elevated temperature and the cleaved RNA fragments were used for first strand cDNA synthesis using reverse transcriptase and random primers. This was followed by second strand cDNA synthesis using DNA polymerase I and RNaseH. Then, cDNA fragments were perfored an end repair process and ligation of adapters. The product was purified and enriched with PCR to create the final cDNA library [26]. Finally, the cDNA library was sequenced by the Illumina HiSeq™ 2000 sequencing platform.

Analysis of Illumina sequencing results
The raw reads from mRNA-seq were filtered by discarding the reads with adaptor contamination, masking low-quality reads with ambiguous 'N' bases and removing the reads in which more than 10% bases had a Q-value <20 [27]. The clean reads were assembled into contigs using the Trinity program [28]. In the Trinity method, an optimized k-mer length of 25 was used for de novo assembly. According to the paired-end information of the sequences, the contigs were linked into transcripts. Subsequently, the transcripts were clustered based on nucleotide sequence identity, and the longest transcripts in the cluster units were regarded as unigenes to eliminate redundant sequences.
ESTScan was used to predict its coding regions and sequence direction if a unigene was not aligned in any of the above databases [29]. To reflect the molar concentration of a transcript by normalizing for RNA length and for the total read number, the expression abundance of the unigenes was represented using the number of fragments per kilobase of transcript per million fragments mapped (FPKM) [30]. Functional annotation by gene ontology terms (GO, http:// www.geneontology.org) was analyzed by Blast2GO [31] and WEGO [32] software. The KEGG pathways annotation was performed using Blastall software against the KEGG database. KEGG pathways were retrieved from KEGG web server (http://www.genome.jp/kegg/).

Detection of SSR and SNP markers
To dig the potential SSR markers, the assembled sequences were searched using MISA software (http://pgrc.ipk-gatersleben.de/misa/). The parameters were designed for the identification of perfect dinucleotide motifs with a minimum of six repeats, and tri-, tetra-, penta-, and hexanucleotide motifs with a minimum of five repeats [33,34]. Also, SOAPsnp was used to detect the single-nucleotide polymorphism (SNP) in the transcript. The program can assemble consensus sequence for the genome of a newly sequenced individual based on the alignment of the raw sequencing reads on the Unigenes. The SNPs can then be identified on the consensus sequence through the comparison with the Unigenes [35].

Gene validation and expression analysis
To validate the transcript sequence and analyze the temporal expression, six genes were mined from the D. sinicus transcriptome database including PsbA, PsaB, PetB, PetF, Lhca2 and F5H. Of which, the first five genes were involved in photosynthesis and the last one gene was associated with lignin biosynthesis. The expression level of genes in different internodes was detected by Real-time qPCR, and UBQ (Ubiquitin) gene used to as the internal control. The mRNA was reverse transcribed into complementary DNA (cDNA) according to the instruction of the Superscript III First-Strand Synthesis System (Invitrogen, Carlsbad, USA). The primers used for qPCR were designed using primer premier 3.0 (http://bioinfo.ut.ee/primer3-0.4.0/) and synthesized from Sangon Biotech Company (Shanghai, China). The corresponding gene names, sequences and primers used for RT-qPCR analysis were displayed in File S1. The qPCR reaction was implemented according to the protocol of SYBR 1 Premix DimerEraser™ (TaKaRa, Dalian, China) by 7300 Real Time PCR System (Applied Biosystems, CA, USA). The cycling conditions of PCR reaction were recommended by the manufacturer (30 s at 95°C, 40 cycles of 95°C for 5 s, and 60°C for 31 s). Three biological replicates of each sample and triplicates of each reaction were performed.

Transcriptome sequencing and de novo assembly
For getting a broad gene library related to culm development, RNA of nine developmental stages of culms was pooled. A total of 65.08 million raw reads and 4.59 gigabase pairs (Gbp) with an average GC content of 53.21% were obtained by a stringent quality check. The reads with Q 20 and no ambiguous "N" were defined as high-quality reads. Using SOAPdenovo [36], 54,807,622 high-quality reads were assembled into 240,630 contigs after removal of adaptor sequences and exclusion of contaminated or short reads. By the Trinity de novo assembly program, short-read sequences were assembled into 81,744 unigenes with a mean length of 723 bp, of which, 63,130 unigenes (77.23%) with length 200-1000 bp, 14,161 unigenes (17.32%) with length 1000-2000 bp, and 4453 unigenes (5.45%) with length > 2000 bp ( Table 1). The N50 values of contigs and unigenes were 310 bp and 1095 bp, respectively. The number of reads based on the relative position in gene (5 0 -3 0 ) presented normal distribution (Fig 1). There was a positive relationship between the length of a given unigene and the number of reads (Fig  2), indicating a randomly fragmented transcriptome in this study. What's more, the raw paired-end sequence data with FASTQ format was deposited in the NCBI Sequence Read Archive (SRA) database with accession number SRA302259, facilitating the access and use of the D. sinicus transcriptome sequencing data.

Functional annotation and classification
A threshold of 10 −5 was adopted by performing a BLASTX search against diverse protein databases, including the NCBI nonredundant protein (Nr) database, NCBI non-redundant nucleotide sequence (Nt) database, UniProt/Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), Cluster of Orthologous Groups of proteins (COG) and Gene Ontology (GO), and 78.71% of unigenes (64,338) were annotated. According to the BLASTX results, 56,587 (69.22%) unigenes showed significant similarity to Nr protein database, and 61,256 (74.94%)  Table 2). Total of 22,569 (40%) unigenes displayed significant homology with sequences of Oryza, and 18% and 14% of the mapped sequences have a high similarity with sequences of Brachypodium distachyon and Setaria italica, respectively. Interestingly, only 509 unigenes (1%) have a high similarity with sequences of P. edulis (Fig 3). Overall, GO was used to classify the functions of the assembled transcripts and describe gene products in terms of their associated biological processes, cellular components, and molecular functions [37]. According to GO classification, 42,508 unigenes were divided into the three categories including 52.57% in biological processes, 28.81% in cellular components, and 18.62% in molecular functions (Fig 4). The four largest biological processes were cellular processes, metabolic processes, response to stimulus, and biological regulation. Under the cellular components, the major classifications were cell, cell part, organelle, and membrane. The most of the genes were classified into the molecular functions of binding, catalytic activity, transporter activity, and nucleic acid binding transcription factor activity. These results indicated that culm development was involved in many fundamentally biological regulation and metabolism.

SSR and SNP marker detection
A total of 7,353 sequences containing 8,553 SSRs were identified from 81,744 unigenes, with 997 unigene sequences containing more than one SSR, and 487 unigene sequences containing SSR present in compound formation. Among them, the number of trinucleotide motifs and dinucleotide motifs were the maximum, accounting for 53.83 and 31.12%, respectively (Table 3). In addition, the most abundant repeat type was AG/CT (1,941), followed by CCG/ CGG (1,669), AGG/CCT (918), AGC/CTG (748), GA (611), respectively. However, in the all unigene sequences, the number of quad-, penta-and hexa-nucleotide motifs was less than 5% (Fig 6).
By detection of 81,744 Unigenes (total length 59,137,193 bp) sequence information, a total of 81,534 SNPs were identified with the frequency of 1/725, including transition 52,462 and transversion 29,072. A/G and C/T are the main types, accounting for 32.43% and 31.91% of all  (Table 4).

Functional genes involved in lignin biosynthesis
As is well known, the lignin content of bamboo is higher than most herbaceous plants [38], which may be due to the differences in the number or level of expression of key enzymes involved in lignin biosynthesis [18]. In the present study, 388 unigenes encoding 16 key enzymes, from the 81,744 unigenes, were identified involved in lignin biosynthesis (Table 5 and Fig 7). The numbers of putative unigenes related to lignin biosynthesis were compared based on Da bamboo of this study, Ma bamboo transcriptome data, Moso bamboo cDNA Database, and rice genes identified from the genome sequences, of which, genes encoding peroxidase remarkably had the most abundance. From the rice genome database [39], we found abundant of transcripts including 26     phylogenetic tree was constructed by alignments of the protein sequences containing condensation domain, which can examine the phylogenetic relationship between the condensation domain proteins in D. sinicus and other plants (Fig 8). Based on the phylogenetic tree, twentyseven sequences identified from D. sinicus were divided into two clades (A and B). In clade A, seventeen sequences were mainly classed to dicotyledon, however, in clade B, ten sequences were generally grouped with monocotyledon including Elaeis guineensis and Phoenix dactylifera, implying the structural or functional diversity of HCT unigenes in D. sinicus.
Generally, the transcript number involved in lignin biosynthesis in Moso bamboo is smaller than that in Da bamboo and Ma bamboo. The significant difference among Da bamboo, Ma bamboo and Moso bamboo is not surprising because the 10,608 Moso FL-cDNAs actually represent only one third to one fourth of the estimated total of Moso bamboo genes. Also, the Moso FL-cDNA libraries were constructed from shoots, leaves and roots from germinating seeds which were not the most representative tissues for high lignin content, while Da bamboo and Ma bamboo materials used for the transcriptome sequencing covered as many tissues as possible including culms of different developing periods [18]. The above results mean that the unique features of gene expression involved in lignin biosynthesis of Da bamboo.

Functional genes involved in transcription factor and photosynthesis
According to the previous research, many transcription factor families play key roles in plant growth, development and immunity [40,41,42]. After BLASTn analysis, many unigenes were putatively identified as transcription factors, including ERF, Myb, Zinc finger, WRKY, Homeobox, bZIP, bHLH, NAC, MADS, etc (Table 6). Remarkably, Zinc finger had the most abundance in Ma bamboo and Da bamboo, while ERF was the most abundant in Moso bamboo. In the Da bamboo, the order of the number of transcription factor was as follows: Zinc finger, bHLH, Myb, WRKY, Homeobox, NAC, ERF, MADS, bZIP and CBF/NF-Y/archaeal histone, which was completely different from that of Moso bamboo, and was slightly different when compared with Ma bamboo. Taking into account that Da bamboo and Ma bamboo belong to sympodial bamboo and Moso bamboo was classified as scattered bamboo, the  distinctions of abundance and order of the number of transcription factor might reflect the species differences.
A total of 109 transcripts encoding photosynthesis were identified (Fig 9A), including photosystem I and II, cytochrome b6/f complex, photosynthetic electron transport and F-type ATPase, and 24 transcripts were characterized as antenna proteins that regarded as the main tool for capturing light of plants (Fig 9B), which imply that culms involve in the photosynthesis during its height development, and that culm photosynthesis plays a certain role in the biomass production. In addition, to confirm the reliability of transcriptome data, the expression level of 14 unigenes, including PsbA, PsaB, PetB (chloroplast-encoded genes involved in photosynthesis), PetF, Pgk, Lhca2, Lhca3, Lhca6 (nuclear-encoded genes involved in photosynthesis), AtpC (nuclear-encoded genes involved in Calvin cycle), and F5H, CAD, COMT, 4CL, HCT (associated with lignin biosynthesis) were detected using real-time quantitative PCR (Fig 10). The results showed the expression level of 14 genes accorded with the transcriptome data.

Discussion
The significance of the BLAST comparison depends in part on the length of the query sequence, and short reads obtained from sequencing would rarely be matched to known genes [43]. In the present study, 78.71% of unigenes were annotated, which is relative high in contrast to other uncharacterized plant such as Litsea cubeba (56.00%) [37], Siraitia grosvenorii (59.90%) [26], P. heterocycla (69.75%) [19] and D. latiflorus (78.90%). The unmapped unigenes can be ascribed to the short sequence reads generated by the sequencing technology and the relatively short sequences of the resulting unigenes, most of which probably lack the conserved functional   domains [44]. There are other possible reasons that some of these unigenes might be non-coding RNAs [18], and that the insufficient sequences of bamboo in public databases also influence the annotation efficiency [44]. Given the existence of species-specific genes, even among the species with high genomic synteny, many genes had big difference in their syntenic regions [45]. Thus, although the draft genome of moso bamboo [21], transcriptome data of other bamboo species [18,19] were reported, the present study is necessary to elucidate the clum development of D. sinicus. Thiel et al. (2003) [46] indicated that approximately 3-7% of expressed genes contain putative SSR motifs, mainly within the un-translated regions of the mRNA. SSRs within gene sequences may have different putative functions, for example, SSR variations in 59-untranslated regions (UTRs) could regulate gene expression by affecting transcription and translation; SSR expansions in the 39-UTRs cause transcription slippage and produce expanded mRNA; intronic SSRs can affect gene transcription, mRNA splicing, or export to cytoplasm; SSRs within genes should be subjected to stronger selective pressure than other genomic regions [47]. In the present study, approximately 9% of the 7,353 unigenes contain single sequence repeat (SSR), which is similar to that of Ma bamboo (12.8%), and is much smaller than that of Moso bamboo (24%) and rice (44%). This result might be related to species difference and sampling strategy.
The synthesis of monolignols, lignin precursor, is implicated in many enzymatic reactions, which involves the general phenylpropanoid pathway starting with the deamination of phenylalanine and leading to the production of hydroxycinnamoyl CoA esters [48]. Although lignin is the most abundant phenylpropanoid derived from the hydroxycinnamoyl-CoA esters, the latter are also the precursors of a wide range of end products including flavonoids, anthocyanins and condensed tannins, which vary according to species, cell type and environmental signals [49]. In order to produce monolignols, hydroxycinnamoyl-CoA esters undergo successive hydroxylation and O-methylation of their aromatic rings involving the following enzymatic activities: shikimate O-hydroxycinnamoyltransferase (HCT); caffeoyl shikimate esterase (CSE); p-coumarate 3-hydroxylase (C3 0 H); caffeoyl CoA 3-O-methyltransferase (CCoAOMT); ferulate 5-hydroxylase (F5H) and caffeate/5-hydroxyferulate O-methyltransferase (COMT) [50,51].
HCT uses p-coumaroyl-CoA and caffeoyl-CoA as preferential substrates to transfer an acyl group to the acceptor compound shikimic acid, yielding p-coumaroyl shikimate. Usually, the differential expression of HCT gene was involved in defense against abiotic and biotic stress, for example, down regulation of HCT gene of Leymus chinensis's leaves resulted from mechanical wounding [52], and down regulation of HCT gene was found when compared the susceptible with resisitant genotypes to leaf miner (Leucoptera coffella) [53]. Hoffmann et al. (2004) [54] pointed out that silencing of the HCT gene induced a dwarf phenotype and changes in lignin composition, suggesting that HCT catalyzes the reactions both immediately preceding and following the insertion of the 3-hydroxyl group by C3H into monolignol precursors. In this study, due to specific expression of HCT compared to other bamboos and rice, we conjecture that it might involve in the high lignin content and giant culm in D. sinicus.
In comparison to most other eukaryotes, plant genomes contain a higher proportion of recently duplicated genes [55]. These duplicates are mostly derived from segmental, whole genome duplication, and tandem duplication events [56,57]. Except segmental duplication and whole genome duplication, tandem duplication, which produces duplicates that are located in close proximity, has contributed significantly to the expansion of plant gene families [58,59]. Compared with the whole genome duplication, tandem duplications have more frequent occurrence and are responsible for much of the gene copy number and allelic variation within a population [60,61]. Although each tandem duplication event only affects a small number of genes, tandemly duplicated genes constitute approximately 14% of all duplicates in Arabidopsis [58]. In particular, genes involved in stress responses have an elevated probability of retention in a single-lineage fashion following tandem duplication, suggesting that these tandem duplicates are likely important for adaptive evolution to rapidly changing environments [62]. For example, in the Eucalyptus grandis, of the tandem duplicated genes EgrHCT1-5, only EgrHCT4 and 5 are likely to be associated with developmental lignification; the other three have only residual expression in xylem and are inducible by abiotic stresses to different levels, suggesting neo-or subfunctionalization following duplication from a common ancestor [62]. In the present study, HCT gene families were divided into obviously different and regular clades, which is likely consistent with gene tandem duplication. This needs further study by quantitative analysis of gene expression. Like as arbor species with stem, bamboos have culm that originated from the mother rhizome. The culm is covered by numerous overlapping sheaths that will shed off with the culm development. The growth process of culm is very strange. It initially hid under the ground surface, however, after growth starting, reached full height within a single growth season at an amazing elongation rate, e.g. for some bamboos, of over 100 cm per day [63]. It was also observed that culm did not show any diameter increment during or after the elongation period. The diameter with which it emerges remains unchanged throughout its life [25]. Because it is leafless and enwrapped by the rigid and nontransparent sheath during the culm development, the action of reserve nutrients was considered as the main source for the growth of bamboo culm and rhizome [25]. However, the transcripts encoding photosynthesis identified from culms indicated that culms photosynthesis associated with its height development.
Woody tissue photosynthesis has been hypothesized to have enabled the maintenance of biodiversity and terrestrial plant productivity throughout the Pleistocene, when atmospheric CO 2 concentrations were only c. 200 ppm [64]. Even today, with atmospheric CO 2 concentrations c. 400 ppm, recycling of CO 2 by woody tissue photosynthesis can instantaneously offset between 7% and 98% of carbon loss in stems and branches of different species in summer, with a median of 72% [65,66]. The contribution of the woody tissue photosynthesis to the overall tree carbon balance largely differs between and even within tree species, depending on size, age and environment [67]. Irradiance is considered to be one of the major limiting factors because photosynthesis requires adequate photosynthetically active radiation to pass through the epidermal, peridermal, and/or rhytidomal layers to reach the light harvesting complexes of the chloroplasts [65,66,68]. In general, the productive months of culm emergence have days of longer photoperiod. Thus, it appears that the period of culm emergence also depend on climatic condition of locality [25]. Heavy shade of upper tree canopy inhibits bamboo regeneration and growth. It has been observed that clumps growing in the open sites produce culms of much better quality and quantity than the clumps growing under heavy shade [69]. Given that the leafless culm enwrapped by sheath and the culm usually grown under crowns, it seems that culms hardly capture enough light for assimilating carbon. According to views reported by Ávila et al. (2014) [65], the stem photosynthesis can be distinguished to two types: stem net photosynthesis (SNP), which includes net CO 2 fixation by stems with stomata in the epidermis and net corticular CO 2 fixation in suberized stems, and stem recycling photosynthesis (SRP), which defines CO 2 recycling in suberized stems. It is widely accepted that the respiration of stem inner tissues is the CO 2 source for SRP; however, labeled carbon transported upward by a xylem sap from roots to leaves is fixed in branches of Platanus occidentalis [70] and in detached branches of Populus deltoids [71], suggesting strongly that CO 2 recycled by SNP or SRP may come also from the roots. Because development, maturation, and aging gradually performed from basal to top internodes, the expression pattern of functional gene among internodes diversified, which involved in upregulated type, down-regulated type, fluctuant type, and discontinuous type. By real-time quantitative PCR analysis (Fig 10), discontinuous expression were displayed in some unigenes, including PsbA (chloroplast-encoded gene related to photosynthesis), AtpC (nuclear-encoded genes involved in Calvin cycle), and CAD, COMT, HCT (associated with lignin biosynthesis), indicating that the comlex expression pattern of genes resulted in development of D. sinicus. As has been noted, both SNP and SRP might involve in bamboo development. What's more, the contribution rate of stem photosynthesis accounting for culm development and the carbon source of stem photosynthesis await further investigation.
Supporting Information S1