Three homozygotes generated from anthers exhibiting different ploidy levels were used in comparison with seedlings cultured from seeds. For DGE sequencing analysis, the triploid line showing retarded growth was designated as Munro_3, the dodecaploid was designated as Munro_12, the hexaploid line was designated as Munro_6, and the seedlings used as the control were designated as Munro_10. Three biological repeats were performed per libraries.

The third true leaves were selected from the first branch of the bamboo stem, and the middle parts of the leaves were immersed into formaldehyde, acetate, and alcohol (FAA) fixative. Small pieces cut from the midrib and intercostal regions of bamboo leaves at median level were dehydrated through a graded ethanol series and embedded in Spurr's resin (Spurr, 1969). Transverse and longitudinal Sections Conclusion μm in thickness were cut with an ultramicrotome (LKB 2188) and stained with a 1% aqueous solution of toluidine blue. Morphology was observed under a light microscope (Olympus BX61). Images were recorded using the Charge-coupled device (CCD) attached to the microscope and analyzed with image processing software (CorelDraw X5). The mean value was averaged from the data from at least 20 individual cells per ploidy level. Stomatal size and density were measured by the silicone rubber impression technique (Weyers and Johansen, 1985). Twenty stomata per genotype lengthwise were measured for stomatal size. the number of stomata per square millimeter was calculated for stomatal density.

Growth properties include plant height, stem diameter, and shoot number were measured and calculated for three periods of shoot development (May to June, July to August, and September to October). Photosynthesis and transpiration rates were examined with an Li-6400 XT Portable Photosynthesis System with an Artificial light source (PAR, 800 μmol m⁻² s⁻¹; Flow, 500 μmol s⁻¹) in the greenhouse at 9:00–11:00 a.m. during the period of rapid growth of Ma bamboo.

Three biological repeats were performed per index. There are no data about the growth and physiological characteristics of the triploid lines due to their lack of vigor. Statistical analysis was carried out to quantitate the differences in characteristics at different ploidy levels using SPSS software version 15. The differences between the means were scored by Duncan's multiple range test (Duncan, 1955).

To compare the differential gene expression patterns among the four libraries, the tag distribution for gene expression level in each library was normalized to obtain an effective library size. Significant differentially expressed transcripts (DETs) with p < 0.05 and log2 fold-change > 1 were extracted using edgeR (Empirical analysis of Digital Gene Expression in R) (Robinson et al., 2010).

Transcript levels were quantified in reads per kilobase of exon model per million mapped reads (RPKM) (Mortazavi et al., 2008). The DETs were blast in the Swiss-prot database, NCBI Non-redundant (Nr) database and Nt database. The differentially expressed sequences were assigned to GO terms describing biological processes, molecular functions and cellular components by Blast2GO (Conesa et al., 2005; Conesa and Götz, 2008). The data represented a GO analysis at level 2, illustrating general functional categories.

Gene functions of the differentially expressed sequences were predicted and classified according to the COG database. KEGG pathways were also analyzed using the online KEGG Automatic Annotation Server (KAAS), http://www.genome.jp/kegg/kaas/. The results of KEGG analysis contained KO assignments and KEGG pathways. Heatmaps of gene expression were generated using Heatmap Illustrator, version 1.0.

The expression levels of 11 genes identified by DGE were validated by quantitative RT-qPCR performed on a 7300 Real Time PCR (Applied Biosystems, Carlsbad, CA, USA). The RNA samples were same as for the sequencing mentioned above. The cDNA was synthesized from 500 ng of total RNA using a PrimeScript 1st Strand cDNA Synthesis kit (TakaRa Bio, Dalian, China). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an endogenous control, as it was the housekeeping gene (Liu et al., 2012). The primers for RT-qPCR were designed using the Primer5 software (Table S1). RT-qPCR was performed using SYBR® Premix Ex Taq™ Kit (TakaRa Bio, Dalian, China). Each cDNA samples were analyzed in triplicate The 2^ΔΔCt method was used to calculate relative expression levels (Livak and Schmittgen, 2001; Moriya et al., 2007).

Network construction was performed to identify modules of highly correlated genes by weighted correlation network analysis (WGCNA) (Zhang and Horvath, 2005; Langfelder and Horvath, 2008). Expression values transformed to log2 values based on the DGE data. Quality control were analysised by clustering and principal component analysis according to the gene expression levels in the DGE data set. The correlation networks were performed by the WGCNA package (version 1.46) in R (version 3.1.3) (Venables et al., 2012). Eigengenes were calculated to visualize gene expression patterns for each gene coexpression module. The co-expressed genes with strong interconnections were defined as hub genes. Subsequently, differentially expressed genes were selected from the hub genes and their edges were analyzed. The results of gene co-expression analysis were visualized by Cytoscape software (version 3.2.1) (Saito et al., 2012).

Bamboo leaves commonly consist of the epidermis, chlorophyll, and vascular tissue. As confirmed by microscopic analysis, the epidermis in Dendrocalamus latiflorus Munro is composed of both long and short cells; the former contain sinuous walls, while the latter often occur in pairs or in isolation, with non-sinuous walls (Soderstrom and Ellis, 1988). Here we found that plants of all ploidies examined showed similar microstructure (Figure 1). Bulliform cells, which are an intrinsic part of the epidermis that allow the leaf to bend or fold inward to avoid excess water loss, were regularly distributed (e.g., Figure 1B2). The vascular system in the bamboo leaf was oriented in the longitudinal direction, with collateral vascular bundles of various sizes enclosed in a double sheath. The outer sheath consisted of parenchymatous cells, and the inner sheath comprised sclerenchymatous pericycle formed by cells with small dimensions and thickened walls (e.g., Figure 1B1). Sheath extension of a sclerenchymatous nature could be observed in the vascular bundles of bamboo leaves (e.g., Figure 1B2). The midrib of the bamboo leaf was distinct and contained several vascular bundles. The remaining space between the vascular tissue and epidermis was mainly filled with two types of thin-walled chlorophyll cells, i.e., fusoid cells and chlorophyllous parenchyma (e.g., Figure 1B2). Fusoid cells were found on both sides of the vascular tissue in perpendicular orientation. These cells likely function in the storage and transport of water in the mesophyll, considering their location next to vascular bundles (Vieira et al., 2002). Chlorophyllous parenchyma with vertically orientated invaginations filled the spaces between the epidermis and fusoid cells.

There were significant variations in the anatomical features of bamboo leaves from plants with different ploidy levels (Figure 2). The average thickness of the dodecaploid bamboo leaves was significantly greater than that of the other leaves, whereas the triploid had the thinnest leaves (Figure 2A). The dodecaploid bamboo leaves had the thickest upper epidermis, whereas that of triploid leaves and seedlings was the thinnest (Figure 2B). No obvious difference in the thickness of the lower epidermis was observed in this study (Figure 2C). The dodecaploid bamboo also showed the biggest fusoid cells (Figures 2D,E). The stomatal size increased from triploids to hexaploids to dodecaploids (Figure 2F), whereas there was no obvious difference in stomatal density between lines.

We measured plant growth properties including plant height, stem diameter, and shoot number continuously from May to October 2013, finding that there were significant differences between seedlings and plants generated from anthers. Although the hexaploid and dodecaploid plants generated from anthers were shorter and had smaller stem diameters than seedlings, they tended to have more shoots. In particular, the number of shoots in dodecaploid bamboo was more than double that of seedlings (Figure 3).

Chlorophyll a, chlorophyll b, and total chlorophyll levels were previously shown to be significantly higher in dodecaploid lines than in hexaploid lines (Qiao et al., 2013). We found that the photosynthesis and transpiration rates exhibited the same trend, as these values frequently increased with ploidy level in D. latiflorus (Figure 4).

To characterize the digital gene expression (DGE) profiles in plants of three different ploidy levels generated from anthers and seedlings grown from seeds, we constructed four DGE libraries and sequenced them using Illumina deep sequencing technology. We calculated the gene expression levels, which are represented by RPKM values, as well as the number of genes belonging to different RPKM intervals (Figure S1). In the triploid (Munro_3) library, the number of genes with RPKM>100 accounted for 40% of expressed genes but only 6% of the total number of genes. Genes with RPKM values between 0 and 10 represented 10% of expressed genes but only approximately 50% of the total number of genes (Figures S1A,B). Genes in hexaploid (Munro_6), dodecaploid (Munro_12), and seedlings (Munro_10) showed similar tendencies. These results indicate that few genes were expressed at high levels, while most were expressed at relatively low levels. Our data reflect the non-uniformity and redundant characteristics of mRNA expression. Figure S2 shows the sequencing saturation levels of the four libraries, revealing that our DGE sequencing covered most of the expressed genes. We tested the randomness of the cDNA fragments by calculating the number of reads matched to different areas of the reference genes (Figure S3). The reads fully covered different loci and were distributed equally, implying that the randomness of the cDNA conformed to the requirements for sequencing and analysis. All of the data were uploaded to an ftp site (ftp.biomarker.com.cn) under the category /zhuory/Munro_DGE (Please contact R. Zhuo for ftp access).

We obtained more than 3.0 million raw tags per library (Table S2). The number of clean tags in each library ranged from 3.03 to 3.34 million. The number of clean tags mapped to transcripts in each library ranged from 2.31 to 2.55 million, accounting for 76.44–76.86% of the total. The number of clean reads that perfectly mapped to transcripts ranged from 1.08 to 1.21 million in different samples, accounting for 45.82–47.36% of the total. Among the four libraries, the number of indel reads ranged from 52,329 to 57,731, whereas the number of indels and mismatched reads ranged from 139,182 to 151,839.

We performed pair-wise comparisons of the four libraries, implementing a total of six comparisons. We found that 981–3,015 transcripts had significant differences in expression, and a total of 8,396 transcripts were differentially expressed (Table 1). The largest differences in expression occurred between Munro_10 and Munro_6. The smallest difference was observed between Munro_10 and Munro_3, in which only 981 DETs were identified.

Several complementary approaches were utilized to annotate the DGE tags. Among the 8396 DETs, 5411, 6497, and 3379 were annotated in the Nr database, Nt database and Swiss-Prot database, respectively. Based on Nr annotation, we carried out Gene Ontology (GO) (Ashburner et al., 2000) analysis. 1118 unigenes were assigned one or more GO terms, including 59.1% in the category molecular functions, 30.1% in biological processes and 10.8% in cellular components (Figure S4). For molecular functions category, Enzyme activity was the most highly represented GO term, followed by binding activity. Genes involved in physiological processes and cellular processes were highly represented in the biological processes category. Cell was the most highly represented category in the cellular components category.

In addition, we subjected all DETs to searches against the Cluster of Orthologous Groups (COG) of proteins database for functional prediction and classification. In total, 1,020 sequences were classified into 25 protein functional families involved in molecular processing, biochemistry metabolism, signal transduction, cellular structure, and so on (Figure 5). The category general function (24.6%) represented the largest group, followed by carbohydrate transport and metabolism (11.7%), replication, recombination and repair (10.9%) and transcription (10.8%).

Annotations of differentially expressed Ma bamboo genes were also subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway analysis. A total of 1,103 unigenes were predicted to 128 pathways involved in crucial physiological processes such as photosynthesis, TCA cycle, phenylpropanoid biosynthesis, as well as genetic information processing such as DNA replication, transcription, translation and so on.

Based on the above annotation results, we are concerned with the DETs involved in regulation of cell proliferation and cell expansion, including transcription regulation, protein biosynthesis and modification, hormonal regulation and cell-wall loosening (Table S3). Among them, we found that most genes of the expansin family were expressed at the highest level in the dodecaploid lines (Figure 6).

To validate the results obtained by digital gene expression tag profiling, we chose representative genes differentially expressed from the four libraries and performed quantitative reverse-transcription PCR (qRT-PCR) to assess their characteristics. Three genes, including the genes encoding trehalose-6-phosphate synthase (Eastmond et al., 2002), neutral trehalase (Liu et al., 2011), and endoglucanase (Westereng et al., 2011), which participate in carbohydrate transport and metabolism, were predicted to be dramatically differentially expressed based on DGE analysis. Our qRT-PCR results showed a similar tendency, although the relative ratios of gene expression differed (Figure 7). The expression of uridine kinase, which is involved in nucleotide transport and metabolism (Chu et al., 2012), was in accordance with the DGE results. Moreover, the differential expression of oxidoreductase genes involved in energy production and conversion, and laccase-like genes associated with lignin biosynthesis, was also validated by qRT-PCR. Overall, there was good agreement between the DGE data and the qRT-PCR results.

To identify coexpression of differentially expressed genes, a network was constructed based on the DGE dataset. Eight gene coexpression modules containing 8,016 unigenes were identified after excluding non- and weakly-expressed genes (Figure S5). These gene modules range in size from 1,969 (turquoise module) to 495 (pink module) genes. Based on this network analysis combined the phenotypic data, we focused on two network hub genes, expansin B3 (EXPB3, BMK12426) and TCP (BMK57883), which play prominent roles in cell growth and differentiation. EXPB3 was expressed at the highest levels in the dodecaploid lines and at the lowest levels in the triploids. The expression of TCP exhibited the opposite trend. In this network, the hub genes EXPB3 and TCP were coexpressed and directly connected with 184 edges and 53 edges, respectively.

We also constructed a subnetwork containing 146 genes by identifying 11 categories related to cell growth, transcription factor, defense response, photosynthesis, and signal transduction, among others (Figure 8). Eleven genes involved in cell growth, such as homologs of β-EXPANSIN (BMK15074, BMK2905, BMK47902), UBC (UBIQUITIN-CONJUGATING ENZYME, BMK18623), and GIDRP88 (GROWTH INHIBITION AND DIFFERENTIATION RELATED PROTEIN 88, BMK24015), may be related to cell wall expansion, cell cycle control, and cell differentiation. Ten photosynthesis-related genes were found in the network, including genes encoding two photosystem II proteins (BMK23946 and BMK31724), chlorophyll a-b binding protein (BMK66421), and glutamyl-tRNA reductase (BMK29030). In addition, 13 and 22 genes belonged to the categories transcription factor and other binding activities, respectively. These genes, such as myeloblastosis (MYB) superfamily genes (BMK38028 and BMK7065), Dof (BMK56897), and basic-leucine zipper (bZIP) gene HY5, regulate many pathways involved in plant growth. We identified 32 genes related to defense responses in the network, which may regulate cellular status to help plants adapt to the environment. Finally, we identified numerous genes with unknown functions in the network (not shown in the figure), perhaps due to the relative lack of gene annotations for bamboo species; these genes represent good candidates for further investigation. Importantly, we can utilize this network to identify novel genes that interact with EXPB3 and TCP in future studies.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.