RNA-Seq analysis of differential gene expression in Betula luminifera xylem during the early stages of tension wood formation

Plant materials

B. luminifera clone (1V25-2) plants used in this study were grown in an environmentally controlled greenhouse of Zhejiang A&F University, Hangzhou, China. The greenhouse was set to a 16-h light/8-h dark cycle, and the air temperature was held between 20 and 28 °C. For reference transcriptome sequencing, different samples were collected from one-year-old trees at 3 pm between April and May 2014. Among these samples, male and female inflorescences and mature leaves were sampled from the upper crown on 4th April and 15th April 2014, respectively. Roots from hydroponic plants and the lignified stems from the upright current-year branches of the upper crown were collected on 12th May 2014. For the RNA-Seq analysis of TW formation, TW induction was performed through the bending method as described by Huang et al. (2014) and Paux et al. (2005). Four time-points, including 0 h, 6 h, 48 h and 7 d, were set up for bending treatment. At each time point, five straight one-year-old plants were bended at the 15th internode (counted from the top) to a 45° angle and fixed with ropes to induce TW formation. The 7 d bending treatment was first started at 3 pm on 12th May 2014, the 48 h bending treatment was then started at 3 pm on 17th May, and the 6 h bending treatment was started 9 am on 19th May. All of the bending treatments ended at 3 pm on 19th May, and the xylem samples were collected immediately after the completion of the bending treatment. The xylem samples from five straight plants without bending treatment (0 h) were also collected at the same time and taken as control (NW). The xylem tissues (approximately 2 mm in thickness) were scraped from the upper side of the bent stems, as well as from the equivalent position of the control after the removal of the bark. To identify the genes preferentially expressed in developing xylem, leaf blades (Leaf) were also prepared from mature leaves by removal of the central vein. All samples used for RNA isolation were immediately frozen in liquid nitrogen and stored at −80 °C.

Anatomical analysis

Xylem samples were collected after the bending treatment for six months. To determine the G-layers for the presence of TW, 20 µm thick sections were examined using scanning electron microscopy (Foston et al., 2011), and 10 µm sections in thickness were double-stained with Safranin O and Fast Green FCF. The lignified walls are stained red with Safranin O, whereas Fast Green stains the pure cellulosic G-layer blue–green. Furthermore, fluorescence microscopy was used to directly observe the autofluorescence of lignin in unstained sections (Selig et al., 2012). The wall thickness of fiber was measured using a modification of Franklin’s method (McLaughlin & Schuck, 1991). At least five different plants were used, and 50 fibers from each plant were measured using a light microscope, and the data were analyzed for significant differences using SPSS 19.0 (IBM) with Student’s t test.

Measurements of sucrose, fructose and IAA

The Sucrose/Fructose/D-Glucose assay kit (Megazyme, Wicklow, Ireland) was used to determine sucrose and fructose contents by following the manufacturer’s instructions. The contents of indole-3-acetic acid (IAA) were determined as described by Yang et al. (2001). Each treatment had five biological replicates. The data were analyzed for significant differences using SPSS 19.0 (IBM) with Student’s t test.

Total RNA isolation

Total RNA from the different tissues was extracted with the PureLink™ Plant RNA Reagent (Invitrogen, CA, USA) according to the manufacturer’s protocol. To remove DNA contaminants, total RNA samples were then treated with RQ1 RNase-free DNase (Promega, WI, USA). RNA integrity was determined using a 2100 Bioanalyzer RNA Nanochip (Agilent, CA, USA) with an RNA Integrity Number (RIN) value of more than 7.0. The RNA concentration was determined using a NanoDrop ND-1000 spectrophotometer (Nano-Drop, DE, USA).

Construction of a cDNA library and reference transcriptome sequencing

A mixture of equal amounts of RNA from mature leaves, stems, roots, inflorescences and TWs was used for reference transcriptome sequencing. The cDNA library with fragments between 200 and 700 bp was constructed using the Illumina kit according to the manufacturer’s instructions (Illumina, CA, USA). Transcriptome sequencing was performed using Illumina HiSeq 2000 (100-bp paired-end reads) at Beijing Genomics Institute (BGI, Shenzhen, China).

Analysis of reference transcriptome sequence data

Raw reads were filtered to obtain high-quality clean reads by removing adaptor sequences, reads containing more than 5% ambiguous bases and low-quality reads containing more than 20% bases with a Q-value <10. The de novo assembly of clean reads was carried out with the short-read assembly program Trinity with default parameters (Grabherr et al., 2011). Trinity first combines reads with a certain length of overlap to form longer fragments, which are called contigs. Then, the reads are mapped back to contigs; with paired-end reads it is able to detect contigs from the same transcript and the distances between these contigs. Finally, Trinity connects the contigs and obtains sequences that cannot be extended on either end. Such sequences were defined as Unigenes and were used for a blastx search against four public protein databases with an E-value ≤1e−5. These protein databases included the NCBI non-redundant proteins (NR) database (http://www.ncbi.nlm.nih.gov), the Swiss-Prot protein database (http://www.expasy.ch/sprot), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg), and the Clusters of Orthologous Groups of proteins (COG) database (http://www.ncbi.nlm.nih.gov/COG). Unigenes with highest sequence similarity were retrieved for subsequent analysis. Based on NR annotation, Gene Ontology (GO) functional classification was obtained by WEGO (http://wego.genomics.org.cn/) (Kotake et al., 2011) via GO terms assigned by Blast2GO (version 2.3.5) with an E-value threshold of 1e−5 (Karpinska et al., 2004). The possible functions of Unigenes were also predicted and classified by aligning to the COG database with an E-value threshold of 1e−5. Pathway assignments were carried out according to the KEGG database using blastx with an E-value threshold of 1e−5.

Profiling differentially expressed genes

The RNA was isolated from the xylem sample of each plant, and then the RNA samples from the five individual plants of each time point were equally mixed to construct four TW cDNA libraries. Another cDNA library using leaf RNA samples derived from leaf blades without the central vein was also constructed to identify the xylem-specific genes. A total of five cDNA libraries were prepared according to the manufacturer’s instructions (Illumina, San Diego, CA, USA). These cDNA libraries were sequenced via Illumina HiSeq2000 (1 ×50-bp read length) by BGI in Shenzhen, China. Raw reads were first filtered by removing adaptors, reads containing more than 10% ambiguous bases, and low-quality reads containing more than 50% bases with a Q-value ≤5. All clean reads were mapped to the reference transcriptome sequences described above using SOAPaligner 2.0 with the parameters “-m 0 –x 1,000 –s 40 –l 35 –v 3 –r 2” (Li et al., 2009). Then, the number of mapped clean reads for each Unigene was calculated and normalized to a Reads Per Kilobase per Million clean reads (RPKM) value (Mortazavi et al., 2008). Based on the RPKM values differentially expressed Unigenes during the early stages of TW formation were identified according to the report of Audic and Claverie (Audic & Claverie, 1997). The false discovery rate (FDR) was calculated via SAS 8.0 and was used to identify the p-value threshold in multiple tests and analyses (Benjamini & Yekutieli, 2001). Compared to NW, the differentially expressed genes (DEGs) during the early stages of TW formation were screened with the threshold of FDR <0.001 and the absolute value of log₂(ratio) ≥1. Overall expression patterns of Unigenes, excluding those with nonsignificant expression changes or with changes over 32-fold, were visualized in 3D with MATLAB (R2012). Among the DEGs, those that showed up- or down-regulated expression at all tension time-points were taken as up- or down-regulated Unigenes, respectively. Then, the GO functional classifications of the up-regulated and down-regulated Unigenes were performed using WEGO.

In addition, B. luminifera Unigenes, which are potentially involved in the cellulose and lignin biosynthetic pathways and auxin regulation, were identified according to the method described by Shi et al. (2011). To reflect dynamic changes and absolute expression patterns at the early stages of TW formation, eight different colors were used to represent different RPKM values of these Unigenes. The RPKM values were converted to logarithm base 2, and clustering was performed using MultiExperiment Viewer (MeV) (version 4.9.0) with the algorithm “K Means Clustering”. Cellulose and lignin biosynthetic pathways were produced manually as described by Huang et al. (2012).

To verify the quality of assembled Unigenes, 20 genes involved in cellulose and lignin biosynthesis (Table S1) were isolated using the SMARTer™ RACE cDNA amplification kit (Clontech, CA, USA) according to the manufacturer’s instructions. Based on the sequences of RACE cloning, full-length cDNAs of these genes were cloned by RT-PCRs using the PrimeScript^®RT-PCR kit (TaKaRa, Dalian, China). The primers used are listed in Table S2. The open read frames (ORFs) of the putative full-length cDNA sequences were predicted using the online ORF finder program (https://www.ncbi.nlm.nih.gov/orffinder/), and the amino acids were deduced using DNASTAR 7.0. Phylogenetic analysis of genes was performed using MEGA4.0 (Tamura et al., 2007). The GenBank accession numbers of other gene sequences are listed in Table S3 . To evaluate the quality of assembled Unigenes, full-length sequences of these cloned genes were aligned with related Unigenes using the program AlignX (Invitrogen, CA, USA).

Co-expression analysis

Co-expressed gene networks were constructed by means of Mutual Rank (MR). First, Unigenes supported by less than 10 reads in each subsample were filtered, and the sample redundancy was calculated with RPKM. Then, PCCs (Pearson’s correlation coefficients) between gene expression patterns were calculated following known formulas (http://atted.jp/help/coex_cal.shtml) and were converted into MRs with the method described by (Obayashi et al., 2007). For each guide gene, the top 20 related genes were selected to produce the co-expressed gene relationships. Finally, co-expression networks were visualized using Cytoscape 2.8.3 (Smoot et al., 2011).

Quantitative RT-PCR validation and expression analysis

Gene expression patterns in different tissues were analyzed by qRT-PCR using the CFX96™ Real-Time PCR Detection System (Bio-Rad). Twenty genes involved in cellulose and lignin biosynthesis and the IAA signaling pathway were validated. All PCR reactions were performed with five replicates. The actin gene (GenBank accession number: FJ410442), which is stably expressed in different tissues/organs, was chosen as an internal reference for normalization. The primers used for qRT-PCR are listed in Table S2.

Reference transcriptome analysis of paired-end sequencing data

First, a reference transcriptome of B. luminifera was generated by paired-end sequencing. A total of 70,358,340 clean reads with an average length of 90 bp were obtained and assembled into 45,700 Unigenes with an average length of 775 bp and a combined length of 35.4 Mb. Of all the Unigenes, 23,149 (50.7%) were over 500 bp long, and 6,025 were longer than 1,500 bp (Fig. S1A). The mean sequencing depth of the Unigenes was 34.02, while 87.9% of the Unigenes were mapped by more than 10 reads, and 39.1% of them were mapped by more than 100 reads (Fig. S1B). To verify the quality of the assembled Unigenes, 20 genes involved in cellulose and lignin biosynthesis were cloned and then compared with related Unigenes by sequence alignment as shown in Table S1, 27 Unigenes were found to cover whole or partial regions of the cloned genes with sequence identities over 98%. For example, the Unigene 4901.6 covered 99.3% of the BlCesA2 with a sequence identity up to 99.6%. This indicated that the assembled Unigenes of the B. luminifera transcriptome were high in quality. Therefore, they could be used as a reference transcriptome for transcriptomic profiling with RNA-Seq.

Functional annotation and classification of the reference transcriptome

To annotate the assembled Unigenes, sequence similarity searches were performed against four public databases. A total of 35,135 (76.9%) sequences showed significant similarity to known proteins in at least one of the searched databases, and 9,600 (21.0%) had significant matches in all four databases (Fig. 1A). According to the NR annotations, 81.6% of the mapped sequences showed strong homology (E-value <1e−20), and 54.4% showed very strong homology (E-value <1e−50) to available plant protein sequences (Fig. 1B). As shown in Fig. 1C, over 87% of the Unigenes could be annotated with proteins from four top-hit species, i.e., Vitis vinifera, Ricinus communis, Populus trichocarpa, and Glycine max. Based on the GO annotations, 27,531 Unigenes were assigned to 59 subcategories (Fig. S2). In the three main GO categories (cellular component, biological process and molecular function), ‘cell’, ‘metabolic process’ and ‘catalytic activity’ were the most abundant groups, respectively. A total of 18,125 Unigenes belonging to the ‘metabolic process’ group could be used to identify the genes involved in the secondary growth of the cell wall during wood formation.

Moreover, the assembled Unigenes were compared against the COG database for in-depth analysis of phylogenetically widespread domain families. Of all the 45,700 Unigenes, 12,224 had a COG classification. Among the 25 COG categories, ‘general function prediction only’ (4,075 Unigenes, 33.3%) and ‘nuclear structure’ (three Unigenes, 0.02%) represented the largest and smallest groups respectively (Fig. S3). It is noteworthy that 1,125 Unigenes were assigned to the category of cell wall/membrane/envelope biogenesis.

Pathway-based analysis, performed by searching against the KEGG database, helps to further understand the biological functions and interactions of genes. There were 19,474 Unigenes mapped to 128 KEGG pathways. The pathways with the highest Unigene representation were metabolic pathways (Ko01100, 4,653 Unigenes, 23.89%), followed by secondary metabolites (Ko01110, 2,199 Unigenes, 11.29%) and plant-pathogen interactions (Ko04626, 1,176 Unigenes, 6.04%) (Table S4). These annotations would be a valuable resource for further research on specific processes, functions and pathways of wood formation. As an example, phenylpropanoid biosynthesis, an important pathway associated with lignification, was examined. Overall, 384 Unigenes were mapped to this pathway, which contains more than 50% of the known enzymes (Fig. S4).

Anatomical characteristics of bending-induced tension wood

Xylem samples of B. luminifera were taken after six months of bending and were transversely sectioned. The lumen diameters of fiber cells in TW were obviously smaller than those in NW, while the fiber cell walls in TW were much thicker (Figs. 2A, 2B, 2D and 2E). Further measurement indicated that no significant difference was detected in the fiber width between TW and NW, but the mean wall thickness of TW was 1.5-fold that of NW (Fig. 2G). This means that cell wall thickening mainly occurred on one side of the fiber lumen. Transverse sections were also double-stained with Safranin O and Fast Green FCF, by which the lignin of the cell wall is stained red with Safranin O, and celluloses are stained blue–green with Fast Green. As shown in Fig. 2C, the secondary cell walls (SCWs) in TW were almost all stained bluish green, whereas those in NW were almost all red (Fig. 2F). This indicated the presence of a cellulosic G-layer in TW. These findings were consistent with previous reports regarding other species, such as Populus euramericana, Eucalyptus globulus, Eucalyptus gunnii and Fraxinus mandshurica, where the G-fibers were formed under similar treatments (Jiang et al., 2009; Jourez, Riboux & Leclercq, 2001; Paux et al., 2005; Toghraie et al., 2006). In general, typical TW can be clearly generated when B. luminifera is under bending treatment.

Changes in sucrose, fructose and IAA contents during the early stage of TW formation

The hydrolysis of sucrose catalyzed by sucrose synthase (SUS) generates UDP-glucose to provide the substrate for cellulose biosynthesis. As shown in Fig. 3A, the sucrose contents at three time-points presented nonsignificant differences compared with NW. Fructose, another product of sucrose decomposition, was obviously accumulated in TW after 48 h bending, and the corresponding differences between TW and NW reached 1.34 and 1.39 mg g⁻¹DW, respectively (Fig. 3B). Early studies also demonstrated that indole-3-acetic acid (IAA) was tightly involved in TW formation of hardwood species (Kennedy & Farrar, 1965; Leach & Wareing, 1967; Robnett & Morey, 1973). In comparison with NW, a significant reduction in the IAA concentration in TW was observed after 7 d of bending, but no obvious changes occurred after 6 h or 48 h of bending (Fig. 3C).

Global transcriptome analysis during TW formation

Identification and annotation of differentially expressed genes (DEGs)

Having generated a reference transcriptome for B. luminifera, our next goal was to perform quantitative measurements of the transcriptome in the developing xylem at the early stages of TW formation. To achieve this, the global gene expression was determined by RNA-Seq in different xylem samples, including 0 h (NW, normal wood) and three different tension woods (TW6 h, bent for 6 h; TW48 h, bent for 48 h; and TW7d, bent for 7 d). According to the RPKM and FDR methods (Audic & Claverie, 1997; Mortazavi et al., 2008), 21,621 Unigenes were expressed during TW formation, with 4,348, 4,983 and 5,977 Unigenes showing differential expression between NW and TW6 h, NW and TW48 h, and NW and TW7d, respectively (Fig. 4A). Taken together, 13,273 DEGs were found at the early stages of TW formation. Among these DEGs, 35.5% were also preferentially expressed in NW compared with the Leaf sample (Fig. 4A). The expression of DEGs exhibited multiple patterns, which are shown in Fig. 4B. During this process, 3,650 and 5,749 DEGs showed significantly up-regulated and down-regulated expression, respectively. It is worth noting that the lowly expressed Unigenes could have changes over 1,000-fold. The up-regulated and down-regulated Unigenes could be categorized into 38 GO functional groups (Fig. 4C). The GO terms of metabolic process and catalytic activity were two major responsive groups in which the numbers of up-regulated Unigenes were significantly higher than those of the down-regulated Unigenes, suggesting that metabolic processes were enhanced by the bending treatment.

TW is characterized by a relatively higher cellulose content and a lower degree of lignification (Foston et al., 2011; Pilate et al., 2004; Wada et al., 1995); therefore expression analysis of genes involved in cellulose and lignin biosynthesis could provide a basis for elucidating the molecular mechanism of TW formation. In total, 125 Unigenes encoding 21 enzymes related to cellulose and lignin biosynthesis were expressed at the early stages of TW formation, and 77 of them were differentially expressed, including 42 up-regulated genes, 29 down-regulated genes and six genes without clear expression patterns (Tables 1 and 2 and Table S5). Among these DEGs, 43 (55.8%) Unigenes were preferentially expressed in NW (Tables 1 and 2). In the biosynthetic pathway of cellulose, 15 Unigenes encoding the related enzymes presented significantly up-regulated expression (Fig. 5, Table 1 and Table S5), which could be attributed to the fact that TW produces the cellulose-rich G-layer. Specifically, five Unigenes (4901.4, 4901.5, 4901.7, 824.3 and 824.4) encoding CesA and SUS exhibited dramatic up-regulation (—log₂(TW/NW)— ≥4) and were strongly expressed in the developing xylem. In the monolignol biosynthesis pathway, 66 related Unigenes were identified, and six Unigenes annotated as CCR, CAD, and SAD were apparently down-regulated (Fig. 6, Table 2 and Table S5). In contrast, 27 genes were significantly up-regulated during the early stages of TW formation (Fig. 6, Table 2 and Table S5). For example, the transcripts of eight DEGs annotated as PAL, C4H, C3H and HCT were increased significantly, seven of which were preferentially expressed in the NW compared with Leaf (Fig. 6).

Table 1:

Candidate genes involved in cellulose biosynthesis related metabolism.

Gene	Enzyme	KO id (EC-No.)	No. Alla	No. Changeb	No. Upc	No. Downd	No.Xyleme
CesA	Cellulose synthase	K10999 (EC:2.4.1.12)	21	13	7	5	5
SUS	Sucrose synthase	K00695 (EC:2.4.1.13)	7	2	2	0	2
GTF	Glucosyltransferase	K05841 (EC:2.4.1.173)	4	0	0	0	2
UGP	UDP-Glucose Pyrophosphorylase	K00963 (EC:2.7.7.9)	3	3	2	1	1
PGM	Phosphoglucomutase	K01840 (EC:5.4.2.8)	5	3	1	2	0
HK	Hexokinase	K00844 (EC:2.7.1.1)	5	2	1	1	0
FRK	Fructokinase	K00847 (EC:2.7.1.4)	9	6	1	5	3
KOR	endo-1,4-beta-glucanase	K01179 (EC:3.2.1.4)	5	5	1	1	2

DOI: 10.7717/peerj.5427/table-1

Notes:

aNo. All indicates the total number of Unigenes analyzed.

bNo. Change indicates the number of Unigenes with expression significantly changed during tension wood formation.

cNo. Up indicates the number of Unigenes with expression significantly up-regulated during tension wood formation.

dNo. Down indicates the number of Unigenes with expression significantly down-regulated during tension wood formation.

eNo. Xylem indicates the number of DEGs preferentially expressed in the xylem (${log}_2\frac{\text{Leaf}}{\mathrm{NW}}<-1$).

Table 2:

Candidate genes involved in lignin biosynthesis related metabolism.

Gene	Enzyme	KO id (EC-No.)	No. Alla	No. Changeb	No. Upc	No. Downd	No. Xyleme
PAL	Phenylalanine ammonia lyase	K10775 (EC:4.3.1.24)	3	3	2	0	2
C4H	Cinnamate 4-hydroxylase	K00487 (EC:1.14.13.11)	2	2	2	0	1
4CL	4-coumarate CoA ligase	K01904 (EC:6.2.1.12)	11	7	3	4	3
C3H	p-coumarate 3-hydroxylase	K05280 (EC:1.14.13.21)	2	2	2	0	2
HCT	p-hydroxycinnamoyl CoA shikimate /quinate p-hydroxycinnamoyl transferase	K13065 (EC:2.3.1.133)	2	2	2	0	2
CCoAOMT	Caffeoyl CoA O-methyltransferase	K00588 (EC:2.1.1.104)	5	4	3	0	1
F5H	Ferulate 5-hydroxylase	K09755 (EC:1.14.-.-)	2	2	2	0	2
COMT	Caffeic acid O-mothyltransferase	K05279 (EC:2.1.1.76)	7	3	1	2	3
CCR	Cinnamoyl CoA reductase	K09753 (EC:1.2.1.44)	14	6	2	4	4
CAD	Cinnamyl alcohol dehydrogenase	K00083 (EC:1.1.1.195)	5	2	1	1	2
SAD	Sinapyl alcohol dehydrogenase	K00083 (EC:1.1.1.195)	1	1	0	1	1
SAMS	S-adenosylmethionine synthetase	K00789 (EC:2.5.1.6)	5	5	5	0	3
POD	Peroxidase	K00430 (EC:1.11.1.7)	7	4	2	2	2

DOI: 10.7717/peerj.5427/table-2

Notes:

aNo. All indicates the total number of Unigenes analyzed.

bNo. Change indicates the number of Unigenes with expression significantly changed during tension wood formation.

cNo. Up indicates the number of Unigenes with expression significantly up-regulated during tension wood formation.

dNo. Down indicates the number of Unigenes with expression significantly down-regulated during tension wood formation.

eNo. Xylem indicates the number of DEGs expressed preferentially in the xylem (${log}_2\frac{\text{Leaf}}{\mathrm{NW}}<-1$).

Indole acetic acid (auxin) is a key regulator of wood development (Nilsson et al., 2008). The expression patterns of auxin-regulated genes in this study are shown in Fig. 7. In total, 42 Unigenes annotated as auxin-responsive factors (ARFs) were expressed during the bending treatment (Table S6), 25 of these ARFs were preferentially expressed in NW when compared to Leaf (Fig. 7). Compared with NW, 26 (61.9%) DEGs were found, including 14 down-regulated Unigenes (Fig. 7). It is worth noting that the transcript levels of 13 Unigenes and 18 Unigenes were changed significantly after 6 h and 48 h of bending treatment, respectively. Among the 34 Unigenes annotated as auxin-induced proteins (Fig. 7 and Table S6), 22 Unigenes presented the significant expression changes, including seven up-regulated and 11 down-regulated Unigenes (Fig. 7 and Table S6). Furthermore, when 6 h and 48 h of bending were applied, 12 and 16 DEGs in TW were detected, respectively. In addition, 10 Unigenes encoding auxin-transport carrier proteins were expressed over the period of TW formation, and only three Unigenes annotated as efflux carrier proteins presented obvious expression changes including one up-regulated gene (Fig. 7, Table S6).

Identification of genes co-expressed with cellulose synthesis genes

To derive useful biological information, 12 up-regulated DEGs coding for CesA, SUS, KOR, UGP and FRK were selected as “guide genes” to identify co-expressed genes using RPKM data from five libraries. According to the MR method, a co-expression network was built based on the 20 most strongly correlated genes for each guide gene, and 215 Unigenes with similar expression trends were examined (Fig. S5). Of 198 Unigenes with annotations, 18 Unigenes coded for transcription factors (TFs), such as MYB TFs, Knotted1-likehomebox (KNOX) TFs and OVATE family proteins (OFPs) (Table S7). For instance, Unigene 3087 encoding a MYB TF was mainly expressed in the xylem and significantly increased during the early stages of TW formation (Table S7). The expression of Unigene 3087 was highly similar to that of BlCesA1 (Unigene 4901.5). Furthermore, Unigenes 13166 and 19576 were both annotated as OFPs, and their expression patterns were similar to those of Unigene 4487 (BlKOR1) and Unigene 6504.2 coding for UGP, respectively (Table S7).

Validation of the gene expression profiles by qRT-PCR

To validate the gene expression profiles, 20 genes were selected for the transcript analysis by qRT-PCR, and the results were compared to the RNA-Seq data (Fig. 8). As a result, the expression of 15 genes (BlCesA1, BlCesA3, BlCesA4, BlCesA8, BlSUS2, BlKOR1, Unigene14665, Unigene431.16, Unigene13207, Unigene5935, Unigene5662.1, Unigene9691, Unigene7840, Unigene611.1 and Unigene13162) generally agreed with the results of RNA-Seq results (Fig. 8). While, the expression levels of BlCesA2, BlCesA5, BlCesA6, BlCesA7 and Unigene6481.2 exhibited inconsistency in some samples, and the fold differences were all less than one (Fig. 8). These results showed similar expression patterns for most genes between RNA-Seq and qRT-PCR, and indicated the reliability of the RNA-Seq data.