Comparative Genomic Analysis of Transgenic Poplar Dwarf Mutant Reveals Numerous Differentially Expressed Genes Involved in Energy Flow

In our previous research, the Tamarix androssowii LEA gene (Tamarix androssowii late embryogenesis abundant protein Mrna, GenBank ID: DQ663481) was transferred into Populus simonii × Populus nigra. Among the eleven transgenic lines, one exhibited a dwarf phenotype compared to the wild type and other transgenic lines, named dwf1. To uncover the mechanisms underlying this phenotype, digital gene expression libraries were produced from dwf1, wild-type, and other normal transgenic lines, XL-5 and XL-6. Gene expression profile analysis indicated that dwf1 had a unique gene expression pattern in comparison to the other two transgenic lines. Finally, a total of 1246 dwf1-unique differentially expressed genes were identified. These genes were further subjected to gene ontology and pathway analysis. Results indicated that photosynthesis and carbohydrate metabolism related genes were significantly affected. In addition, many transcription factors genes were also differentially expressed in dwf1. These various differentially expressed genes may be critical for dwarf mutant formation; thus, the findings presented here might provide insight for our understanding of the mechanisms of tree growth and development.


Introduction
Photosynthesis is a process used by plants and other organisms to convert light energy into chemical energy and store it in the bonds of sugars [1,2]. The process of photosynthesis takes place in the chloroplasts and employs the green pigment chlorophyll [3][4][5]. There are two parts to photosynthesis, namely, the light reaction and the dark reaction [6][7][8]. The light reaction occurs in the thylakoid membrane and converts light energy into chemical energy. This chemical reaction must, therefore, take place in the light. The light-harvesting complex is a complex of subunit proteins that make up part of a larger supercomplex, the photosystem, which is the functional unit of photosynthesis [9][10][11]. These complexes consist of proteins and photosynthetic pigments and surround a photosynthetic reaction center to focus energy, attained from photons absorbed by the pigment, toward the reaction center using Förster resonance energy transfer. Chlorophylls and carotenoids are important components of light-harvesting complexes in plants. The energy harvested via the light reaction is stored by forming ATP, a compound used by cells for energy storage.
The dark reaction, which takes place in the stroma within the chloroplast, converts CO2 to sugar. This reaction does not directly require light, but it does require the products of the light reaction (ATP and NADPH) [12,13]. The dark reaction involves the Calvin cycle, in which CO2 and energy from ATP are used to form sugar [14][15][16].
Monosaccharides are used as energy reserves in plants. Although a plethora of substrates such as proteins and lipids can be oxidized in plants, respiration tends to be dominated by carbohydrate oxidation through the glycolytic pathway and the tricarboxylic acid (TCA) or citric acid cycle. Carbohydrate is converted into pyruvate and malate through two major pathways, i.e., glycolysis and the oxidative pentose phosphate pathway. The primary roles of the pentose phosphate pathway are to generate NADPH for use in biosynthetic reactions, and to provide ribose 5-phosphate for nucleotide synthesis and erythrose 4-phosphate for the synthesis of shikimic acid derivatives [17].
In our previous study, the Tamarix androssowii LEA gene (lea IV; DQ663481; TaLEA) was introduced into Populus simonii × Populus nigra to improve salt tolerance [18,19]. During these studies, it was noted that one of these mutants displayed a dwarf phenotype compared with the other transgenic lines. Next generation (NG) sequencing methods have emerged as a cost-effective high-throughput approach to sequence a very large number of expressed genes, even in small experiments [20][21][22][23]. To uncover the mechanism of dwarf formation, digital gene expression (DGE) libraries were constructed and sequenced from the following lines: wild type, transgenic lines 5 and 6 (XL-5 and XL-6, i.e., transgenic lines with normal phenotypes), and dwf1. In addition, we used gene ontology (GO) analysis to elucidate the functions of the genes that were found to be differentially expressed in the dwarf mutant.

Isolation and Morphological Characterization of the dwf1 Phenotype
Among the 15 independent TaLEA transgenic lines, one grows extremely slowly compared to the wild type and the other transgenic lines, named dwf1, (Figure 1A-F). In addition to slow-growth, dwf1 also exhibits dark green leaves. Chlorophyll content analysis indicated that dwf1 had higher chlorophyll content ( Figure 1G). These characteristics indicate that dwf1 had multiple morphological defects ( Figure 1).

Figure 1.
Phenotype of dwf1. (A) Seedlings of XL-5, XL-6, dwf1 and wild-type. Under the same growth condition, dwf1 exhibited a dwarf phenotype compared with the XL-5, XL-6 and wild-type; (B) Leaves of dwf1 and wild-type. The leaf shape of dwf1 was slightly affected; (C-F) Seedlings of dwf1 and wild-type grown in vitro. (F) Dwf1 exhibited reduced root biomass and fewer lateral roots compared with (E) wild-type; and (G) Chlorophyll content of XL-5, XL-6, dwf1 and wild-type. Error bar indicates standard error (SD). Single factor ANOVA, followed by Tukey Kramer post-hoc test was employed in this study. Sixty leaves from six individuals of each transgenic line and wild-type were measured. Letters (a and b) on the graph denote significant differences (p < 0.05).

Identification and Validation of Differentially Expressed Genes
To identify genes showing significant changes in expression, we compared the levels of gene expression in XL-5, XL-6 and dwf1 with the wild type. We identified 1804 (1489 upregulated and 315 downregulated), 2183 (1874 upregulated and 309 downregulated), and 2163 (1318 upregulated 315 downregulated), 2183 (1874 upregulated and 309 downregulated), and 2163 (1318 upregulated and 845 downregulated) differentially expressed genes in XL-5, XL-6 and dwf1, respectively. We further analyzed the intersection of differentially expressed genes. As shown in Figure 2A, most (947) of the differentially expressed genes (DEGs) in XL-5 and XL-6, shared similar expression patterns. In contrast, only 210 genes had similar expression patterns among these three transgenic lines. Instead, 1246 DEGs were dwf1 unique. Of these dwf1-unique genes, 158 had opposite expression profiles between dwf1 and the other two transgenic lines, XL-5 and XL-6. 1088 DEGs (491 upregulated, 597 downregulated) were only identified in dwf1. Detailed information is presented in Figure 2. , and dwf1 indicate differently expressed genes in each sample in compared to wild-type (WT); "XL-5, XL-6" indicates common differently expressed genes in XL-5 and XL-6; "XL-5, XL-6, and dwf1" indicates common differently expressed genes in XL-5, XL-6 and dwf1; "dwf1 special" indicates differentially expressed genes that were only identified in dwf1; (B) Expression patterns of dwf1 special up-regulated genes in XL-5 and XL-6; (C) Expression patterns of dwf1 special down-regulated genes in XL-5 and XL-6.
To validation the digital gene expression (DGE), 18 differentially expressed genes were randomly selected for qPCR analysis. Fourteen of them showed similar expression patterns to those detected by DGE data (Figure 3), thus the results were reliable for identification of differentially expressed genes.

Pathway Analysis
Pathway analysis was performed using KEGG Automatic Annotation Server (KAAS) (http://www.genome.jp/tools/kaas/) [25]. The 1246 dwf1-specific differentially expressed genes were mapped to more than 100 pathways. As shown in Table 2, the results were coincident with the GO enrichment analysis. Numerous genes were mapped to photosynthesis, photosynthesis-antenna proteins, carbon fixation in photosynthetic organisms, and porphyrin and chlorophyll metabolism. Furthermore, carbohydrate metabolism was affected. The most representative pathways were glycolysis/gluconeogenesis, the citrate cycle (TCA cycle), the pentose phosphate pathway, starch and sucrose metabolism, and pyruvate metabolism (Table S1).

Transcription Factors
Transcription factors are essential for the regulation of gene expression. They are vital for many important cellular processes. In this study, among the 1246 dwf1-unique DEGs, we identified 90 transcription factors, including 55 upregulated and 35 downregulated genes (Table 3). Among these 90 genes, 18 were MYB or MYB-related transcription factors, 12 were bHLH family proteins, eight were bZIPs, seven were NACs, five were C2H2 transcription factors, four were WRKYs, and the remaining genes included AP2, B3, RAV and others.

Dwf1 Showed a Unique Gene Expression Profile Compared to the Other Transgenic Lines
In our previous study, a TaLEA gene (GenBank: DQ663481.1) was introduced into poplar (Populus simonii × P. nigra) by Agrobacterium tumefaciens-mediated transformation [18,19]. All of the 11 transgenic lines grew normally except dwf1, so we analyzed the transcriptional differences among XL-5, XL-6, dwf1 and wild type. As expected, XL-5 and XL-6 shared more DEGs (Figure 2), and XL-5, XL-6 and dwf1 shared fewer. Instead, a larger number of dwf1-unique DEGs were identified (1246 genes). More than 150 genes showed a contrasting expression profile in dwf1 compared to XL-5 and XL-6. The unique expression pattern of dwf1 indicated that the molecular events occurring in dwf1 were different from the other transgenic lines. Indeed, further analysis of the dwf1-unique DEGs may provide useful insights for further investigation of the molecular mechanisms underlying this mutant.

T-DNA-Affected Genes Might Play a Role in Dwarfism Formation
As only dwf1 exhibited a dwarf phenotype, the T-DNA insertion sites of dwf1 may play a role in the formation of this mutant. The T-DNA-affected genes or the interaction between TaLEA and T-DNA-affected genes may have affected the dwf1 gene expression pattern. In our previous work, two T-DNA insertion sites and a T-DNA-affected gene (RAV) were isolated, and qRT-PCR analysis showed that RAV was specifically upregulated in dwf1 [18,19]. Genes associated with starch, sucrose, galactose and glycerolipid metabolism were different expressed in dwf1 compared to wild-type. These results were consistent with this study, however only dwf1 and wild-type were analyzed in previous studies and it was impossible to identity dwf1-unique DEGs. In order to overcome this difficulty, we sequenced other transgenic lines (XL-5 and XL-6), allowing more useful dwf1-unique DEGs relevant for uncovering the mechanism of dwarfism.

Genes Involved in Photosynthesis and Photosynthesis-Antenna Proteins
Many genes involved in photosynthesis and photosynthesis-antenna proteins were specifically differentially expressed in dwf1 ( Table 2). Photosystem II PsbW protein (POPTR_0005s24030.1), photosystem II Psb27 protein (POPTR_0002s05720.1), photosystem II 10 kDa protein (POPTR_0001s42970.1) and photosystem II 10 kDa protein (POPTR_0011s14540.1) were differentially expressed in dwf1. PsbW, Psb27, and photosystem II 10 kDa protein are the main components of photosystem II, and photosystem II is the first protein complex in the light-dependent reaction. The process of photosynthesis requires a variety of energy transducing protein complexes that transport electrons and pump protons, leading to the formation of ATP and NADPH. Photosystem II (PS II) is one of these protein complexes, and it has been the focus of many studies. PS II transfers electrons from water to plastoquinone, and this process generates a pH gradient. Plastoquinone (PQ) carries electrons from PS II to the cytochrome bf complex. This complex traps light and uses it to reduce the QB plastoquinone and to synthesize molecular oxygen from water [26][27][28][29]. Photosystem I (PS I) is the second photosystem in plants, and it contains many subunits. PS I subunit XI (POPTR_0002s24070.1) and PS I subunit PsaN (POPTR_0007s04160.1) were found to be specifically upregulated in dwf1 (Table 2). ATPases are a class of enzymes that catalyze the decomposition of ATP into ADP and a free phosphate ion. This dephosphorylation reaction releases energy [30,31]. In this study, an F-type H + -transporting ATPase subunit b (POPTR_0018s01830.1) was found to be specifically upregulated in dwf1.
The antennae of plants consist of a large number of protein-bound chlorophyll molecules that absorb photons and transfer their energy to the reaction center. The antennae of plants consist of an inner part and an outer part. The outer part, formed by the light harvesting complexes (LHCs), collects light. The inner part of the antenna, consisting of core complexes, is an integral constituent of the reaction centers. LHCs are formed by polypeptides, which bind chl-a, chl-b, xanthophylls and carotenes [32,33]. In this study, seven LHC chlorophyll a/b binding proteins were found to be specifically upregulated in dwf1, including six LHC II chlorophyll a/b binding proteins and one LHC I chlorophyll a/b binding protein ( Table 2).
The altered expression profiles, especially the upregulated expression profile of these photosynthesis-related genes, indicate that dwf1 might have an increased photosynthetic efficiency. We therefore measured the photosynthetic rates and photosynthetic fluorescence parameters of wild type, XL-5, XL-6 and dwf1. Unexpectedly, the altered gene expression in dwf1 did not result in any detectable differences in photosynthetic efficiency or photosynthetic fluorescence parameters (Figure 4).

Genes Involved in Carbohydrate Metabolism
According to the DGE data, the expression of numerous genes involved in carbohydrate metabolism was specifically altered in dwf1, including many key enzymes involved in glycolysis/gluconeogenesis, the citrate cycle (TCA cycle), pentose phosphate pathway, pentose and glucuronate interconversions, fructose and mannose metabolism, galactose metabolism, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, and pyruvate metabolism (Table S1). Glucose-6-phosphate isomerase is an enzyme that catalyzes the conversion of glucose-6-phosphate into fructose 6-phosphate in the second step of glycolysis. Glucose-6-phosphate isomerase is also involved in the pentose phosphate pathway and starch and sucrose metabolism [34][35][36]. In this study, a gene that encodes a glucose-6-phosphate isomerase was found to be specifically downregulated in dwf1. Phosphoglucomutase is an enzyme that transfers a phosphate group on the α-D-glucose monomer from the 1' to the 6' position in the forward direction or the 6' to the 1' position in the reverse direction. This reaction is involved in glycolysis/gluconeogenesis, the pentose phosphate pathway, galactose metabolism, and starch and sucrose metabolism [37][38][39]. In this study, two genes that encode phosphoglucomutase were found to be specifically downregulated in dwf1. Glucose is a major product of photosynthesis, and is stored mainly in the form of starch granules in plastids such as chloroplasts, and especially in amyloplasts. Starch biosynthesis in higher plants is catalyzed by four classes of enzymes, including ADP-Glc pyrophosphorylase (AGPase), starch synthase (SS), starch branching enzyme (BE), and starch debranching enzyme. SS elongates glucans by adding Glc residues from ADP-Glc to the glucan nonreducing ends through α-1,4 linkages. In this study, a starch synthase was found to be specifically upregulated in dwf1.

Transcription Factors
Transcription factors are shaped so precisely that they are able to bind preferentially to the regulatory regions of only a few genes out of the thousands present in a cell's DNA. Transcription factors are involved in almost all biological processes, such as development, responses to intercellular signals, responses to the environment, and pathogenesis. In this study, 90 transcription factors were found to be specifically differentially expressed in dwf1. Among these genes, 17 displayed an expression profile in dwf1 that was the opposite of that observed in XL-5 and XL-6, and 73 were found to be differentially expressed in dwf1, while the expression of these genes was not altered in XL-5 or XL-6 ( Table 3).
MYB transcription factors are known to play an important role in the control of phenylpropanoid metabolism [40,41]. Plant MYB proteins also play a role in hormonal responses during seed development and germination. In this study, more than ten MYB and MYB-related transcription factors were found to be differentially expressed. We also found that WRKY, NAC, bZIP, AP2, RAV and other transcription factors were differentially expressed in dwf1 [42][43][44][45]. All of these differentially expressed genes may play a vital role in dwarf mutant formation.

Plant Growth Conditions and Transgenic Plants
The transgenic poplar (Populus simonii × P. nigra) lines were produced as previously described [46]. The seedlings were planted in a mixture of turfy peat and sand (2:1 v/v) and transferred to a greenhouse at 75% relative humidity and a constant temperature of 24 °C. All wild-type and transgenic plants used for analysis were three months old. The third and fourth young leaves, counted from top, were harvested from six plants of each transgenic line and wild-type. Leaves harvested from two plants of each transgenic line and wild-type were pooled together as one replicate. The leaves were immediately immersed in liquid nitrogen and stored at −70 °C for subsequent RNA extraction.

Chlorophyll Content
A commercial chlorophyll meter (SPAD-502, Spectrum Technologies, Chicago, IL, USA) was used to estimate the chlorophyll contents of leaves. Six independent plants from each transgenic line and the wild type were randomly selected. And then ten leaves of each plant were measured.

RNA Isolation and RNA-seq Library Preparation
Total RNA was isolated from each sample using the CTAB method [47]. All RNA samples were quantified and examined spectrophotometrically for protein contamination (A260/A280 ratios) and reagent contamination (A260/A230 ratios). Extracted RNA samples were selected based on the 28S/18S rRNA band intensity (2:1), spectroscopic A260/A280 readings between 1.8 and 2.0, and A260/A230 readings greater than 1.5. According to the protocol, we generate three libraries for each sample. Subsequently, sequencing by synthesis (SBS) was performed using four types of nucleotides identified by labeling with four different colors.

Data Analysis
Raw reads were filtered through the Illumina pipeline. After that, clean tags were mapped to poplar reference sequences (http://www.phytozome.net/poplar) [48] using SOAP (Short Oligonucleotide Analysis Package) [49] A library was generated containing all possible fragments CATG + 17 bases in length among the reference gene sequences. All clean tags were mapped to the reference sequences, allowing for only 1 bp mismatch. Clean tags that mapped to multiple genes were discarded. The number of unambiguous clean tags was calculated for each gene and normalized to the number of transcripts per million clean tags (TPM) [50,51].
Differential expressed genes (DEGs) were obtained by pairwise comparison of transgenic and WT libraries using edgeR/Bioconductor [24]. The significance of DEGs is based on a false discovery rate (FDR) of 0.05 and a log2 fold change of 1.
GO analysis was performed for functional categorization of differentially expressed genes using agriGO software (http://bioinfo.cau.edu.cn/agriGO/index.php) [52], and the p-values were corrected by applying the FDR correction to control falsely rejected hypotheses during GO analysis. The hypergeometric test was used as the statistical test, and an FDR corrected p-value ≤0.05 was used as the cutoff value [25].
Pathway analysis was performed using KEGG Automatic Annotation Server (KAAS) (http://www.genome.jp/tools/kaas/) [25]. Differentially expressed genes were mapped to terms in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Real-Time RT-PCR
cDNA was fist synthesized from purified RNA using Reverse Transcription Kit (Qiagen, Hilden, Germany), and then 1:10 diluted using Rnase-free water. Quantitative PCR was performed using 1 μL of diluted cDNA as a template in a 20 μL of SYBR Green (Toyobo, Osaka, Japan) mixed reaction. The relative expression levels of the selected genes were calculated using the relative 2 −ΔΔCt method [53]. α-tubulin and actin served as the internal reference genes. Results represent mean standard deviation of the three experimental replicates.

Conclusions
This study has demonstrated the usefulness of the digital gene expression (DGE) approach to identify differentially expression genes between dwf1, XL-5, XL-6, and WT. The dwarf mutant was found to specifically exhibit differential expression of numerous genes involved in photosynthesis and photosynthesis-antenna proteins, genes involved in carbohydrate metabolism, and genes encoding transcription factors. These findings might provide a strong basis for future research on dwarf mutant formation.