Identification of genes related to the phenotypic variations of a synthesized Paulownia (Paulownia tomentosa × Paulownia fortunei) autotetraploid

Gene 553 (2014) 75–83

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

Identification of genes related to the phenotypic variations of a synthesized Paulownia (Paulownia tomentosa × Paulownia fortunei) autotetraploid Yongsheng Li a,b, Guoqiang Fan a,b,⁎, Yanpeng Dong a,b, Zhenli Zhao a,b, Minjie Deng a,b, Xibing Cao a, Enkai Xu a, Suyan Niu b a b

Institute of Paulownia, Henan Agricultural University, Zhengzhou, Henan 450002, People's Republic of China College of Forestry, Henan Agricultural University, Zhengzhou, Henan 450002, People's Republic of China

a r t i c l e

i n f o

Article history: Received 22 February 2014 Received in revised form 22 September 2014 Accepted 29 September 2014 Available online 7 October 2014 Keywords: Paulownia Phenotype Sequencing Tetraploid Transcriptome

a b s t r a c t Paulownia is a fast-growing deciduous tree native to China. It has great economic importance for the pulp and paper industries, as well as ecological prominence in forest ecosystems. Paulownia is of much interest to plant breeder keen to explore new plant varieties by selecting on the basis of phenotype. A newly synthesized autotetraploid Paulownia exhibited advanced characteristics, such as greater yield, and higher resistance than the diploid tree. However, tissue-specific transcriptome and genomic data in public databases are not sufficient to understand the molecular mechanisms associated with genome duplication. To evaluate the effects of genome duplication on the phenotypic variations in Paulownia tomentosa × Paulownia fortunei, the transcriptomes of the autotetraploid and diploid Paulownia were compared. Using Illumina sequencing technology, a total of 82,934 All-unigenes with a mean length of 1109 bp were assembled. The data revealed numerous differences in gene expression between the two transcriptomes, including 718 up-regulated and 667 down-regulated differentially expressed genes between the two Paulownia trees. An analysis of the pathway and gene annotations revealed that genes involved in nucleotide sugar metabolism in plant cell walls were down-regulated, and genes involved in the light signal pathway and the biosynthesis of structural polymers were up-regulated in autotetraploid Paulownia. The differentially expressed genes may contribute to the observed phenotypic variations between diploid and autotetraploid Paulownia. These results provide a significant resource for understanding the variations in Paulownia polyploidization and will benefit future breeding work. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Polyploidization is the result of genome duplication events. It is widespread in angiosperms during their evolution and has significantly influenced the speciation and evolution of flowering plants (Wendel, 2000). Duplication of genomic content either of the same (autopolyploidy) or divergent (allopolyploidy) genomes is often associated with the emergence Abbreviations: A, adenosine; Bp, base pair(s); C, cytidine; CAB, chlorophyll a/b-binding; CCA1, circadian clock associated1;DEGs, differentially expressed genes; FDR, false discovery rate;FPKM, fragments per kbper million; G,guanosine; GAE, UDP-GlcA 4-epimerase;GalUA, D-galacturonate; LHY, late elongated hypocotyl; MYB, myeloblastosis; N50, median lengths; NGS, next-generation sequencing; Oligo, oligodeoxyribonucleotide; PCR, polymerase chain reaction; PaWB, Paulownia witches' broom; PTF4, autotetraploid Paulownia tomentosa × Paulownia fortunei; PTF2, diploid Paulownia tomentosa × Paulownia fortunei; PE, UDP-pectinesterase; PG, polygalacturonase; PHY, phytochrome; PIF3, phytochrome-interacting factor 3; Q30, proportion of nucleotides with quality value larger than 30; T, thymidine; TOC1, timing of CAB expression 1; UDP-GalUA, UDP-D-galacturonate; UDP-GlcUA, UDP-D-glucuronate. ⁎ Corresponding author at: College of Forestry, Henan Agricultural University, Henan 450002, People's Republic of China. E-mail address: [email protected] (G. Fan).

http://dx.doi.org/10.1016/j.gene.2014.09.057 0378-1119/© 2014 Elsevier B.V. All rights reserved.

of novel phenotypes not present in the diploid progenitors (Levy and Feldman, 2002). Among the known polyploid plants, some woody plant species, such as Betula, Populus, and Robinia, have shown morphological modifications and elevated resistance to abiotic and biotic stresses (Mu et al., 2012). The special traits have made the use of these woody polyploid species in forestry both possible and practical. Paulownia belongs to the monogeneric family Paulowniaceae and encompasses several species with similar characteristics (CastellanosHernández et al., 2009). These species are fast-growing deciduous trees with great economic importance. It was reported that some Paulownia species can grow over 18 ft as sprouts in one year, and can produce 1 m3 of wood at the age of 5–7 years (Ates et al., 2008). Because of its high biomass production, this species is an excellent source for the pulp and paper industries, and it has been adapted to grow in many countries in East Asia, Europe, and America (Ipekci and Gozukirmizi, 2003). Besides being a valuable timber species, Paulownia tomentosa and Paulownia fortunei and their hybrids were cultivated for afforestation and mine site reclamation because of their fast growth rate and deep root system. Moreover, the strong tolerance of Paulownia to drought, salt, and heavy metal stress has been of great interest for a long time (Jiang et al., 2012).

76

Y. Li et al. / Gene 553 (2014) 75–83

Genetic improvement of woody species plays a critical role in the modern timber industry; thus, many newly engineered cultivars have been produced because of their outstanding traits. In 2010, an autotetraploid P. tomentosa × P. fortunei was synthesized from its diploid using colchicine treatment. Compared with its diploid, the autotetraploid exhibited some morphological changes including increased thickness and size of leaves, larger stomata size, and lower stomata density per leave (Fan et al., 2010). Significant progress has been made in the study of polyploidy, and novel gene expression and regulation are considered to be the primary factors that lead to the phenotypic variations (Adams and Wendel, 2005). Thus, understanding the molecular mechanisms that trigger biosynthesis in cells and tissues is a pressing need. Over the past decades, molecular markers such as restriction fragment length polymorphisms, and DNA markers generated by polymerase chain reaction (PCR) have been used widely in genetic studies of Paulownia. These traditional methods have relatively low sensitivity for the less abundant transcripts, and have low reproducibility (Ozsolak and Milos, 2010). With the development of the high-throughput next-generation sequencing (NGS) technology, complete plant genomes or transcriptomes can be sequenced with great depth, providing information about biological processes in species with, as yet, unknown genome information (Strickler et al., 2012). Moreover, NGS technology which offers a critical opportunity to generate molecular resources is more efficient and economical than traditional technologies. It is believed to be a powerful tool for advanced research in de novo transcriptome sequencing for polyploid plants (Mardis, 2008). So far, very few studies have applied NGS technology on Paulownia species except the one which focused on the key gene expression variations induced by Paulownia witches' broom (PaWB) phytoplasma infection (Liu et al., 2013; Mou et al., 2013). Variations of gene expression patterns related to cell wall biosynthesis and photosynthesis were believed to be the main factor behind the PaWB symptoms (Mou et al., 2013). However, the information concerning the process and the underlying mechanisms in Paulownia polyploidization need to be explored. In the present study, comparisons between the transcriptomes of autotetraploid P. tomentosa × P. fortunei and its diploid were performed using an Illumina/Solexa Genome Analyzer. The obtained sequences will serve as a valuable resource for: (1) understanding novel forms of gene expression through analysis of the metabolic pathways associated with the differentially expressed genes (DEGs) in Paulownia polyploidization, and (2) discovering novel genes that could be applied in plant breeding programs to develop superior varieties. 2. Materials and methods 2.1. Plant materials The experimental materials were obtained from the Institute of Paulownia (Henan Agricultural University, Zhengzhou, Henan, China). Samples for transcriptomic analysis were collected from tissue-cultured plants. The apical leaves of growing seedlings of autotetraploid P. tomentosa × P. fortunei (PTF4) and its diploid (PTF2), and their duplicate samples (PTF4-1, and PTF2-1) were collected and then sub-cultured in 1/2 MS medium supplemented with 25 g·L−1 sucrose under a photoperiod of 16:8 h/light:dark, at 20 °C. After 30 days, top leaves were clipped from three individual Paulownia species and frozen immediately in liquid nitrogen, and stored at −80 °C for RNA extraction. 2.2. RNA extraction Total RNA was extracted separately from each Paulownia sample (8 mg) using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), and purified using a RNeasy MiniElute Cleanup Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. RNA quantity was assessed using a NanoVue UV–vis spectrophotometer (GE Healthcare Bio-Science, Uppsala, Sweden) at an absorbance of 260 nm. RNA quality

was assessed spectrophotometrically at an absorbance ratio of OD260/280 and OD260/230, and the RNA integrity was confirmed by 1% agarose gel electrophoresis. 2.3. cDNA library generation and RNA sequencing Briefly, 8 μg of total RNA was sequenced using Illumina RNA sequencing (RNA-seq) technology. The cDNA library that was used in the sequencing was constructed as follows. (1) After the total RNA was collected, Oligo (dT) attached magnetic beads were used to isolate poly(A) mRNA. (2) The mRNAs were digested to produce short fragments by adding fragmentation buffer. (3) The cleaved RNA fragments to be used as templates were synthesized into first-strand cDNA using a random hexamer primer, followed by synthesis of the second-strand cDNA. (4) Short fragments were purified with a QiaQuick PCR extraction kit (Qiagen, Valencia, CA, USA) and resolved in EB buffer for end reparation and the addition of a poly(A) tail. (5) Suitable fragments with lengths from 300 ± 50 bp were isolated by agarose gel electrophoresis and selected as sequencing templates for PCR amplification. (6) The four cDNA libraries were sequenced on an Illumina HiSeq™ 2000 platform (Beijing Genomics Institute (BGI), Shenzhen, China). 2.4. De novo transcriptome assembly and functional annotation All the RNA-seq raw reads were cleaned by removing the reads with adaptors, low quality sequences, and reads with average length less than 20 bases. Then, the clean reads from each library were de novo assembled using the Trinity software package (Grabherr et al., 2011) into contigs, scaffolds, and finally unigenes. Next, the unigenes were taken to the process of sequence splicing without redundancy, and assembled until they could no longer be extended by the TIGR Gene Indices clustering tools (TGICL). Finally, the unigenes from the two different libraries were combined to generate the All-unigene dataset. The expression of the All-unigenes was calculated using the fragments per kb per million (FPKM) fragments method as described previously (Mortazavi et al., 2008). For the functional annotations, the All-unigene sequences were compared with the sequences in the NCBI Nr (http://www.ncbi. nlm.nih.gov/) databases, and Swiss-Prot (http://www.expasy.ch/ sprot), clusters of orthologous groups (COG) (http://www.ncbi. nlm.nih.gov/COG/), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.ad.jp/kegg/) using BLAST with an E-value cutoff of less than 1.0E− 5. The best alignment results were used to annotate the transcriptome sequences and to determine the orientation of the sequences. When the results from the different databases conflicted, the orientations of the All-unigene sequences were determined using the priority order of Nr, Swiss-Prot, KEGG, and COG. For the All-unigene sequences that were not aligned to any of these databases, ESTScan software (Iseli et al., 1999) was used to determine their coding regions and sequence orientation. The All-unigene sequences were assigned Gene Ontology (GO) terms (http://www. geneontology.org) using the Blast2GO software (Conesa and Götz, 2008), and metabolic pathways were compared with the KEGG database (Kanehisa et al., 2008). 2.5. Analysis of DEGs The comparison of gene expression levels between the autotetraploid and diploid of Paulownia transcriptome sequences was determined using the FPKM method as follows: 6

FPKM ¼

10 C : NL 103

Y. Li et al. / Gene 553 (2014) 75–83

The NOISeq method was used to select DEGs by presenting pair-wise comparison of the samples (Tarazona et al., 2011). A threshold q ≥ 0.8 and an absolute value of log2Ratio N 2 were used as the criteria to judge the significant DEGs.

77

were longer than 500 bp, and 34,264 (41.32%) were longer than 1000 bp (Table 2). The correlation coefficients of the expression of both the duplicate diploid and tetraploid were shown in Fig. 1. Both duplicates showed linear correlations with their corresponding ones. The Pearson r values were 0.90 and 0.81, respectively.

2.6. Gene ontology and pathway enrichment analysis of DEGs 3.2. Annotation of non-redundant unigenes To determine the significance of the GO functional classification and the metabolic pathways in the DEGs compared with the whole transcriptome background, a hypergeometric test of GO and pathway enrichment was performed as follows: m−1 X

P ¼ 1−

i¼0

  M i

N−M n−i

N 

 :

n

In this formula, N is the number of all the genes with GO/KEGG functional annotation, n is the number of DEGs in N, M is the number of all the genes with a specific GO/KEGG annotation, and m is the number of DEGs in M. 2.7. Quantitative real-time (q)RT-PCR analysis of DEGs To confirm the results of the Illumina sequencing analysis, the expression levels of selected candidate genes were measured by qRTPCR. Three developmental stages of Paulownia plants (first stage, 30-day-old in vitro plantlets; second stage, 6-month-old seedlings; third stage, 12-month-old seedlings) were used to isolate total RNA, and each stage of the plant has two biological replications. RNA from the leaf samples of the two Paulownia genotypes (PTF4 and PTF2) was extracted using TRIzol reagent (Sango, Shanghai, China) according to the manufacturer's protocol. Then, the purified RNA was reverse transcribed to first-strand cDNA using a PrimeScript RT reagent Kit (Takara Biotechnology Dalian Co., Ltd., Dalian, China). Next, the resultant cDNA was amplified in a Bio-Rad CFX96TM Real-Time System (Bio-Rad, Hercules, USA) under the following conditions: 50 °C for 2 min, 95 °C for 30 s, followed by 40 cycles of 94 °C for 15 s and 60 °C for 1 min. Three biology duplications of each sample (each sample contained leaves from at least three individual plants) and triplicates of each reaction were analyzed for each of the candidate genes. The relative gene expression was analyzed using the ΔΔCT method and the average threshold cycle (Ct) was normalized with 18 rRNA. 2.8. Data deposition The Paulownia Illumina sequence data have been deposited in the NIH Short Read Archive database (http://www.ncbi.nlm.nih.gov/sra) under accession number SRP034738. 3. Results 3.1. Transcriptome sequencing and assembly The Illumina/Solexa sequencing generated 38,272,198 raw reads comprising 3,827,219,800 nucleotides for PTF2 library, 52,653,774 raw reads comprising 4,738,839,660 nucleotides for PTF2-1 library, 36,440,822 raw reads comprising 3,644,082,200 nucleotides for PTF4 library and 55,267,938 raw reads comprising 4,974,114,420 nucleotides for PTF4-1 library. The Q30s were 93.65%, 98.43%, 91.55%, and 98.52% for PTF2, PTF2-1, PTF4 and PTF4-1, respectively. All the high-quality reads in the four libraries were assembled into unigenes (88,043 for PTF2, 62,266 for PTF2-1, 77,958 for PTF4, and 63,600 for PTF4-1) with average lengths of 932 bp, 597 bp, 863 bp, and 656 bp, respectively (Table 1). Together the assembled sequences produced a total of 82,934 non-redundant unigenes (All-unigenes) with a mean length of 1109 bp, and N50 of 1832 bp. Among these unigenes, 49,900 (60.17%)

Because genomic information for Paulownia is limited, a BLASTX search with a threshold of 10−5 was performed to predict the functions of the All-unigenes against six public databases (Nr, Nt, Swiss-Prot, KEGG, COG and GO). The results indicated that 57,350 unigenes (70% of All-unigenes) shared similarities with known sequences. Among the mapped unigenes, 74.9% showed strong homology (E-value b 1.0E−30), and 26.1% showed some homology (1.0E−30 b E-value b 1.0E−5) with sequences in Nr (Fig. 2A). The similarity distribution of the aligned sequences showed that 25.7% of the unigenes had similarities higher than 80%, and 74.3% had similarities of 17–80% (Fig. 2B). The species distribution showed that 33.8% of the unigenes had highest matches to sequences from Lycopersicon esculentum, followed by Vitis vinifera (28.1%), Amygdalus persica (6.4%), and Ricinus communis (5.4%) (Fig. 2C). 3.3. Functional classification of DEGs In total, 1385 DEGs were identified between the PTF4 and PTF2 transcriptomes (Supplementary file 1). Among these DEGs, 718 were up-regulated and 667 were down-regulated in the autotetraploids compared with the diploids. To understand the functions of the DEGs in the autotetraploids and diploids, GO and KEGG terms and pathways were used to predict the functions of the DEGs. A total of 872 DEGs (63.0%) were assigned terms under the three main GO categories, biological process, cellular component, and molecular function. Under the biological process category, metabolic process (506, 58.06%) and cellular process (557, 63.83%) were the most highly represented terms. Under cellular component, cell (612, 70.24%) and organelles (516, 59.12%) were significantly represented. Catalytic activity (452, 51.81%) and binding activity (404, 46.30%) were the most highly represented terms in the molecular function category (Fig. 3). In addition, oxidoreductase activity, response to karrikin, cellular carbohydrate catabolic process, and polyol metabolic process showed the most significant differences (p-value b 1.0E−11) between PTF4 and PTF2. A total of 653 DEGs were assigned to 119 KEGG pathways (Supplementary file 2). Of these, biosynthesis of secondary metabolites (ko01110), circadian rhythm-plant (ko04712), photosynthesisantenna proteins (ko00196), and galactose metabolism (ko00052) showed the most significant differences (q-value b 1.0E−11) between PTF4 and PTF2. The DEGs related to morphological differences between PTF2 and PTF4 were categorized into three groups. The first group of DEGs was associated with the nucleotide sugars required for cell wall biosynthesis and included two genes that code for pectinesterase (EC:3.1.1.11) and polygalacturonase (EC:3.2.1.15) (Fig. 4). The second group of DEGs was associated with circadian rhythm and included four genes that code for phytochrome-interacting factor 3 (PIF3), pseudo-response regulator 1 (TOC1), MYB-related transcription factor LHY (LHY), and circadian clock associated 1 (CCA1) (Fig. 5). The third group of DEGs was associated with lignin biosynthesis and included seven genes that code for four enzymes, 4-coumarate–CoA ligase (EC: 6.2.1.12), cinnamoyl-CoA reductase (EC:1.2.1.44), cinnamyl-alcohol dehydrogenase (EC:1.1.1.195), and peroxidase (EC:1.11.1.7) (Fig. 6). 3.4. Verification of RNA-seq by qRT-PCR To confirm the results of the transcriptome sequencing analysis, qRT-PCR was performed on nine randomly selected candidate genes,

78

Y. Li et al. / Gene 553 (2014) 75–83

Table 1 Overview of the sequencing and assembly. PTF2-1

PTF4-1

PTF2

PTF4

52,653,774 4,738,839,660 98.43 0.00 46.18

55,267,938 4,974,114,420 98.52 0.00 45.57

38,272,198 3,827,219,800 93.65 0.01 46.39

36,440,822 3,644,082,200 91.55 0.01 47.35

Contigs Number of contigs Total nucleotides (nt) in contigs Average length of contigs (nt) Length of N50 (nt)

106,721 34,863,584 327 501

110,689 37,055,004 335 542

149,151 52,042,494 349 702

137,177 47,805,523 348 711

Unigenes Number of unigenes Total nucleotides (nt) in unigenes Average length of unigenes (nt) Length of N50 (nt)

62,266 37,191,802 597 1006

63,600 41,698,691 656 1134

88,043 82,012,558 932 1727

77,958 67,293,279 863 1616

All-unigenes Number of All-unigenes Total nucleotides (nt) in All-unigenes Average length of All-unigenes (nt) Length of N50 (nt)

82,934 91,942,755 1109 1832

Statistics of data production Number of clean reads Total nucleotides (nt) Q20 percentage (%) N percentage GC percentage (%)

which were used to examine the expression pattern during Paulownia genome duplication. The real-time PCR primers are shown in Table 3. Of the selected genes, two were down-regulated (CL751.Contig2 and CL4377.Contig2) and associated with starch and sucrose metabolism pathway, three (CL9674.Contig8, CL4544.Contig2, and CL1717) were up-regulated and associated with the circadian rhythm pathway, and four (Unigene35753, Unigene33630, CL9867.Contig2_All, and CL541.Contig1) were up-regulated and associated with the phenylpropanoid biosynthesis pathway. The expression patterns revealed by qRT-PCR analysis were consistent with those predicted by analysis of the DEG sequences (Fig. 7). This result validated the gene expression variations obtained by analysis of the unigene sequences.

(Soltis and Soltis, 1999). In the present study, the transcriptome of leaves from tissue-cultured Paulownia was sequenced and assembled using Illumina/Solexa sequencing technology. A total of 1385 DEGs between the diploid and autotetraploid P. tomentosa × P. fortunei plants were obtained. Some of these genes were predicted to be involved in nucleotide sugar metabolism in the plant cell wall, light signal pathway, and biosynthesis of structural polymers, which may be associated with some of the phenotypic variations of autotetraploid Paulownia. To better identify biological relevant DEGs, two libraries from diploid and tetraploid P. tomentosa were created, respectively. The expression's logarithmic values of both the duplicated samples of diploid and autotetraploid Paulownia were significantly correlated. This showed that our experiments were biologically replicable.

4. Discussion Polyploidization is an important mechanism of speciation that results from genome doubling. It is a persistent force in plant evolution, and many existing plant species have undergone a polyploidization event (Gaut and Doebley, 1997). Colchicine is a chemical that can cause chromosome doubling. It has been used to create new plant polyploids, thereby enriching species resources (Chen, 2007). Although the process of genome doubling may induce silencing and loss of duplicated genes, plant polyploids often show novel phenotypic effects compared with their diploid contributors (Adams and Wendel, 2005). A phenotypic comparison between diploid and autotetraploid Paulownia revealed that the higher ploidy exhibited increased leaf size, fiber length, and stomata sizes (Zhang et al., 2013). These observations were consistent with those in other plant species, which suggested that polyploidization may affect diverse organisms in similar ways

4.1. DEGs involved in plant cell wall biosynthesis Leaf size is regulated by cell division and cell expansion. Therefore, larger leaf size might be the result of either a higher cell number because of increased cell division, or larger cells as the result of increased cell growth (Gonzalez et al., 2010; Nahirñak et al., 2012). The comparative significance of cell number versus cell growth in determining leaf size variation has been reported previously. Leaf size enlargement of Arabidopsis thaliana was reported to be determined chiefly by an increase in cell number (Gonzalez et al., 2010). However, according to our earlier study of tetraploid Paulownia, increased cell size was mainly responsible for the enlarged leaf size, as well as for greater thickness of the upper and lower epidermis, and palisade tissue compared with diploid Paulownia (Zhang et al., 2013).

Table 2 Length distribution of unigenes. Unigene length

100–500 bp 500–1000 bp 1000–1500 bp 1500–2000 bp ≥2000 bp

PTF2-1

PTF4-1

PTF2

PTF4

All

Total number/ percentage

Total number/ percentage

Total number/ percentage

Total number/ percentage

Total number/ percentage

40,154 11,050 5461 3019 2582

64.49% 17.75% 8.77% 4.85% 4.15%

38,640 11,637 6214 3729 3380

60.75% 18.30% 9.77% 5.86% 5.31%

44,162 14,250 10,114 7921 11,596

50.16% 16.19% 11.49% 9.00% 13.17%

41,248 12,688 8887 6313 8822

52.91% 16.28% 11.40% 8.10% 11.32%

33,034 15,636 11,362 8995 13,907

39.83% 18.85% 13.70% 10.85% 16.77%

Y. Li et al. / Gene 553 (2014) 75–83

79

Fig. 1. Correlation coefficients of the gene expression of duplicate samples. Y-axis represents the logarithmic value of diploid (A) or tetraploid (B) expression, while X-axis represents the logarithmic value of the corresponding duplicate samples.

Fig. 2. Characteristics of homology search of Illumina sequences against the nr database. (A) E-value distribution of BLAST hits for each unique sequence with a cut-off E-value of 1.0E−5. (B) Similarity distribution of the top BLAST hits for each sequence. (C) Species distribution is shown as a percentage of the total homologous sequences with an E-value of at least 1.0E−5.

80

Y. Li et al. / Gene 553 (2014) 75–83

Fig. 3. Functional classification of PTF unigenes based on three main Gene Ontology (GO) categories: biological process, molecular function and cellular component. The right Y-axis indicates the number of genes in each category. The left Y-axis indicates the percentage of a specific category of genes in the corresponding GO category.

The plant cell wall is composed entirely of polysaccharides (cellulose, hemicellulose, and pectin), structural proteins, and various enzymes; its materials vary with cell type, environmental conditions, and developmental stages (Cosgrove, 2005). Because of the cellulose cell wall, plant cell enlargement is strictly restrained; thus, plant organ size is largely associated with cell wall development. In Paulownia, photosynthesis incorporates fixed carbon into the cell wall carbohydrate, which comprises various nucleotide sugars. Nucleotide sugar pathways involve a series of sophisticated enzymatic reactions. For example, D-galacturonate (GalUA), one of the main nucleotide sugars, is biosynthesized by sequential conversion from UDP-D-glucose. Twelve of the down-regulated DEGs in the PTF4 transcriptome were associated with this process, including genes that code for pectin cell wall degradation enzymes, such as UDPpectinesterase (PE), and polygalacturonase (PG) (Micheli, 2001). One of the precursors necessary for the synthesis of pectin is UDP-Dgalacturonate (UDP-GalUA), the activated form of GalUA. UDP-GalUA is synthesized from UDP-D-glucuronate (UDP-GlcUA) by UDP-GlcA 4-epimerase (GAE). It was reported previously that plant GAEs could be type II membrane proteins targeted to the Golgi, suggesting that GAEs might provide UDP-GalUA to galacturonosyltransferases in pectin synthesis (Goldberg et al., 1996). Pectins are complex polysaccharides generally present in the primary cell wall and in the spaces between cell walls (middle lamellae) where they bind adjacent cell walls together. Pectins form ~35% of the dry weight of dicot cell walls, and when bound to cellulose, pectins confer rigidity to cell walls (Hayashi, 1989). Three genes (CL751.Contig2,

Fig. 4. Pectin biosynthetic pathway in Paulownia. Down-regulated expressed genes in PTF4 are in green boxes.

CL5499.Contig1, and CL9437.Contig2) that code for PE and PG were down-regulated by 2.955-, 2.359-, and 2.331-fold, respectively, in PTF4, which implies that less cell wall degradation occurs in PTF4 compared with in PTF2. Modulation of the degrading enzymes for polysaccharides might be an alternative factor that regulates cell wall extension. This process may explain the enlargement of leaf size observed in autotetraploid Paulownia. 4.2. DEGs involved in light signaling pathway Light is one of the most important environmental factors. Through light signaling pathways, light provides a variety of signals for growth and development in higher plants. In the present study, five genes involved in light signal transduction were up-regulated in PTF4, including genes that code for phytochrome (PHY), phytochrome-interacting factor 3 (PIF3), timing of CAB expression 1 (TOC1), circadian clock associated 1 (CCA1), and late elongated hypocotyl (LHY). Paulownia has evolved sophisticated systems of numerous photoreceptors that regulate various light responses through their interaction with different signaling proteins, such as phytochromes, cryptochromes, and phototropins. Among the identified signaling proteins that interact with the phytochromes, PIF3, a basic helix-loop-helix protein, is the most widely characterized (Ni et al., 1998). By binding to both PHYA and PHYB, PIF3 acts as a light-dependent regulator by inducing light-regulated genes. Previous hypotheses have proposed that PIF3 could have dual functions as both a positive transcription factor for chloroplast development and a negative regulator for hypocotyl elongation (Stephenson et al., 2009). Moreover, a study on A. thaliana seedlings under diurnal light/dark conditions showed that PIF3 could act as a growth promoter by targeting growth-related genes (Soy et al., 2012). Among the DEGs in the present study, a gene (CL9674.Contig8) that codes for PIF3 was found to be up-regulated in PTF4. This finding may suggest that enhanced light reception in autotetraploid Paulownia could contribute to the development of certain organs and increase their sizes. Circadian rhythm generates expressed genes that allow the plant to sense environmental light changes during the day/night cycle. Molecular genetic studies have revealed an autoregulatory loop comprising a pseudo response regulator, TOC1, and two MYB transcription factors, CCA1 and LHY, which together connect the morning- and evening-phase circuits (Levy and Dean, 1998). It was suggested that the expression of

Y. Li et al. / Gene 553 (2014) 75–83

81

Fig. 5. Light signal pathway in Paulownia. Up-regulated expressed genes in PTF4 are in red boxes.

LHY and CCA1 may rise in the morning and inhibit the expression of the evening-expressed TOC1. Accordingly, when LHY and CCA1 levels fall, the expression of TOC1 rises (Alabadı́ et al., 2001). In the present study, two genes (CL4544.Contig2 and CL1717.Contig4) that code for TOC1 and LHY/CCA1, respectively, were found to be up-regulated in PTF4 and were predicted to be involved in the circadian rhythm pathway. Several metabolic developmental processes in plants show a circadian rhythm in leaf movement and stomata opening which is correlated with the CO2 assimilation rate (Wang and Tobin, 1998). The upregulated genes related to circadian rhythm in PTF4 might enhance

the photosynthesis rate in PTF4 and, thus, increase the growth rate. This finding coincides with the results of a previous study that the photosynthesis rate of several autotetraploid Paulownia species was higher than that of diploid Paulownia (Zhang et al., 2013). 4.3. DEGs involved in lignin biosynthesis In vascular plants, the phenylpropanoid pathway is one of the most important secondary metabolism pathways through which fixed carbon from primary metabolism is used to biosynthesize lignin and other

Fig. 6. Lignin biosynthesis pathway occupied with transcripts coding for corresponding enzymes. Up-regulated expressed genes in PTF4 are in red boxes.

82

Y. Li et al. / Gene 553 (2014) 75–83

Table 3 Primers used for real-time quantitative polymerase chain reaction verification. Annotation

Primer

CL751.Contig2

Sequence definition

Size 894

Pectinesterase

CL4377.Contig2

495

Polygalacturonase

CL9674.Contig8

615

Phytochrome-interacting factor 3

CL4544.Contig2

978

Zinc finger protein CONSTANS

CL1717.Contig4

657

Circadian clock associated 1

Unigene35753

933

Peroxidase

Unigene33630

762

Cinnamyl-alcohol dehydrogenase

CL9867.Contig2

924

Cinnamoyl-CoA reductase

CL541.Contig1

1632

CGCTGTCGGGCTGGTTTC CGCTCCTCGGCAAGTTCC AGGTTCTTGTCAATGTGGATGG ACTAAATAACGACGCCCTTCAG CACACTTCCCTTCTCCAATG TTCCTCAAAGCACTTAATCCC GGTGGAGTGGTGGATGATTAC TTCGCTTGTGATACTGGATGG CTGTATGGTCGTGCTTGG GGGTAAGGATGTAGAGGTTTC AATCTCAACGCCCTCATTTCTG TGTTTCCGCCGCTTCTGG ATCTTGGTGCTCCTCTCC TTGCTTCCTTCTTCTTACTAATAG TGCCGTGTTCTATTCCTTC GTCAATGGTGTCTGTCTGG CGTCGTGCTGTGCGTTCTG CGGTGCCTCTGTATCAACTCC

4-Coumarate–CoA ligase

metabolites. Lignin deposited in the spaces between the cellulose, hemicellulose, and pectin components of cell walls, confers the mechanical strength needed for plant growth and resistance to pathogen invasion (Boudet, 2007). During lignin biosynthesis, phenylalanine is converted to a series of acids, comprising cinnamic acid, caffeic acid, ferulic acid, 5-hydroxyferulic acid, and sinapic acid. These acids are then catalyzed

by 4-coumarate–CoA ligase with CoA to form the corresponding CoA acid derivatives and are reduced later by cinnamoyl-CoA reductase and cinnamyl alcohol dehydrogenase to form alcoholic compounds. The last step of this process is the biosynthesis of lignin monomers by peroxidase. Lignin deposition in the cell wall is known to play an important role in plant development, especially during elongation. Suppression of either

Fig. 7. Relative expression levels of nine selected DEGs between PTF2 (blue) and PTF4 (red) were determined by qRT-PCR. The following genes were validated: (A) CL751.Contig2, and (B) CL4377.Contig2 were up-regulated; (C) CL9674.Contig8, (D) CL4544.Contig2, (E) CL1717.Contig4 (F) Unigene35753, (G) Unigene33630, (H) CL9867.Contig2_All, and (I) CL541.Contig1 were down-regulated. Error bars represent the standard deviation of the qRT-PCR signals from each of the three replicates. *Represents statistically significant differences between the diploid and the tetraploid(p-value was less than 0.05).

Y. Li et al. / Gene 553 (2014) 75–83

of the reducing enzymes would inhibit lignin biosynthesis, and likely affect plant growth. For example, the repression of 4-coumarate-CoA in Pinus radiata resulted in retarded growth, and affected the phenotype (Wagner et al., 2009). In the present study, five up-regulated genes (CL541.Contig1, CL9867.Contig2, Unigene 33630, CL3589.Contig 3, and Unigene35753) that coded for 4-coumarate–CoA ligase, cinnamoyl-CoA reductase, cinnamyl alcohol dehydrogenase, and peroxidase, respectively, were predicted to be involved in lignin biosynthesis. These results may provide an explanation for organ size enlargement in autotetraploid Paulownia. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2014.09.057.

Acknowledgments This work was supported by the Joint Funds of the National Natural Science Foundation of China (NSFC) (Grant No. U1204309), the Fund of the Transformation Project of the National Agricultural Scientific and Technological Achievement of China (Grant No. 2012GB2D000271), and the Central Financial Forestry Science Promotion Project (Grant No. GTH [2012] 01).

References Adams, K.L., Wendel, J.F., 2005. Polyploidy and genome evolution in plants. Curr. Opin. Plant Biol. 8, 135–141. Alabadı́, D., Oyama, T., Yanovsky, M.J., Harmon, F.G., Mas, P., Kay, S.A., 2001. Reciprocal regulation between TOC1 and LHY/CCA1 within the Arabidopsis circadian clock. Science 293, 880–883. Ates, S., Ni, Y., Akgul, M., Tozluoglu, A., 2008. Characterization and evaluation of Paulownia elongata as a raw material for paper production. Afr. J. Biotechnol. 7, 4153–4158. Boudet, A.M., 2007. Evolution and current status of research in phenolic compounds. Phytochemistry 68, 2722–2735. Castellanos-Hernández, O.A., Rodríguez-Sahagún, A., Acevedo-Hernández, G.J., RodríguezGaray, B., Cabrera-Ponce, J.L., Herrera-Estrella, L.R., 2009. Transgenic Paulownia elongata SY Hu plants using biolistic-mediated transformation. Plant Cell Tissue Organ Cult. 99, 175–181. Chen, Z.J., 2007. Genetic and epigenetic mechanisms for gene expression and phenotypic variation in plant polyploids. Annu. Rev. Plant Biol. 58, 377. Conesa, A., Götz, S., 2008. Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int. J. Plant Genomics 2008, 1–12. Cosgrove, D.J., 2005. Growth of the plant cell wall. Nat. Rev. Mol. Cell Biol. 6, 850–861. Fan, G.Q., Zhai, X.Q., Wei, Z.Z., Yang, Z.Q., 2010. Induction of autotetraploid from the somatic cell of Paulownia tomentosa × Paulownia fortunei and its in vitro plantlet regeneration. J. Northeast For. Univ. 38, 22–26. Gaut, B.S., Doebley, J.F., 1997. DNA sequence evidence for the segmental allotetraploid origin of maize. Proc. Natl. Acad. Sci. 94, 6809–6814. Goldberg, R., Morvan, C., Jauneau, A., Jarvis, M., 1996. Methyl-esterification, deesterification and gelation of pectins in the primary cell wall. Prog. Biotechnol. 14, 151–172. Gonzalez, N., De Bodt, S., Sulpice, R., Jikumaru, Y., Chae, E., Dhondt, S., Van Daele, T., De Milde, L., Weigel, D., Kamiya, Y., 2010. Increased leaf size: different means to an end. Plant Physiol. 153, 1261–1279.

83

Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., 2011. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat. Biotechnol. 29, 644–652. Hayashi, T., 1989. Xyloglucans in the primary cell wall. Annu. Rev. Plant Biol. 40, 139–168. Ipekci, Z., Gozukirmizi, N., 2003. Direct somatic embryogenesis and synthetic seed production from Paulownia elongata. Plant Cell Rep. 22, 16–24. Iseli, C., Jongeneel, C.V., Bucher, P., 1999. ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. ISMB, pp. 138–148. Jiang, Z.F., Huang, S.Z., Han, Y.L., Zhao, J.Z., Fu, J.J., 2012. Physiological response of Cu and Cu mine tailing remediation of Paulownia fortunei (Seem) Hemsl. Ecotoxicology 21, 759–767. Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., Katayama, T., Kawashima, S., Okuda, S., Tokimatsu, T., 2008. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, D480–D484. Levy, Y.Y., Dean, C., 1998. Control of flowering time. Curr. Opin. Plant Biol. 1, 49–54. Levy, A.A., Feldman, M., 2002. The impact of polyploidy on grass genome evolution. Plant Physiol. 130, 1587–1593. Liu, R.N., Dong, Y.P., Fan, G.Q., Zhao, Z.L., Deng, M.J., Cao, X.B., Niu, S.Y., 2013. Discovery of genes related to witches broom disease in Paulownia tomentosa × Paulownia fortunei by a de novo assembled transcriptome. PLoS One 8, e80238. Mardis, E.R., 2008. The impact of next-generation sequencing technology on genetics. Trends Genet. 24, 133–141. Micheli, F., 2001. Pectin methylesterases: cell wall enzymes with important roles in plant physiology. Trends Plant Sci. 6, 414–419. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-seq. Nat. Methods 5, 621–628. Mou, H.Q., Lu, J., Zhu, S.F., Lin, C.L., Tian, G.Z., Xu, X., Zhao, W.J., 2013. Transcriptomic analysis of Paulownia infected by Paulownia witches'-broom phytoplasma. PLoS One 8, e77217. Mu, H.Z., Liu, Z.J., Lin, L., Li, H.Y., Jiang, J., Liu, G.F., 2012. Transcriptomic analysis of phenotypic changes in birch (Betula platyphylla) autotetraploids. Int. J. Mol. Sci. 13, 13012–13029. Nahirñak, V., Almasia, N.I., Fernandez, P.V., Hopp, H.E., Estevez, J.M., Carrari, F., Vazquez-Rovere, C., 2012. Potato snakin-1 gene silencing affects cell division, primary metabolism, and cell wall composition. Plant Physiol. 158, 252–263. Ni, M., Tepperman, J.M., Quail, P.H., 1998. PIF3, a phytochrome-interacting factor necessary for normal photoinduced signal transduction, is a novel basic helix-loop-helix protein. Cell 95, 657–667. Ozsolak, F., Milos, P.M., 2010. RNA sequencing: advances, challenges and opportunities. Nat. Rev. Genet. 12, 87–98. Soltis, D.E., Soltis, P.S., 1999. Polyploidy: recurrent formation and genome evolution. Trends Ecol. Evol. 14, 348–352. Soy, J., Leivar, P., González‐Schain, N., Sentandreu, M., Prat, S., Quail, P.H., Monte, E., 2012. Phytochrome-imposed oscillations in PIF3 protein abundance regulate hypocotyl growth under diurnal light/dark conditions in Arabidopsis. Plant J. 71, 390–401. Stephenson, P.G., Fankhauser, C., Terry, M.J., 2009. PIF3 is a repressor of chloroplast development. Proc. Natl. Acad. Sci. 106, 7654–7659. Strickler, S.R., Bombarely, A., Mueller, L.A., 2012. Designing a transcriptome nextgeneration sequencing project for a nonmodel plant species1. Am. J. Bot. 99, 257–266. Tarazona, S., Garcia-Alcalde, F., Dopazo, J., Ferrer, A., Conesa, A., 2011. Differential expression in RNA-seq: a matter of depth. Genome Res. 21, 2213–2223. Wagner, A., Donaldson, L., Kim, H., Phillips, L., Flint, H., Steward, D., Torr, K., Koch, G., Schmitt, U., Ralph, J., 2009. Suppression of 4-coumarate–CoA ligase in the coniferous gymnosperm Pinus radiata. Plant Physiol. 149, 370–383. Wang, Z., Tobin, E.M., 1998. Constitutive expression of the CIRCADIAN CLOCK ASSOCIATED 1 (CCA1) gene disrupts circadian rhythms and suppresses its own expression. Cell 93, 1207–1217. Wendel, J.F., 2000. Genome evolution in polyploids. Plant Molecular EvolutionSpringer, pp. 225–249. Zhang, X.S., Zhai, X.Q., Fan, G.Q., Deng, M.J., Zhao, Z.L., 2013. Observation on microstructure of leaves and stress tolerance analysis of different tetraploid Paulownia. J. Henan Agric. Univ. 46, 646–650.