# Genome-Wide Identification of Genes Probably Relevant to the Uniqueness of Tea Plant (Camellia sinensis) and Its Cultivars

Authors: Yan Wei; Wang Jing; Zhou Youxiang; Zhao Mingming; Gong Yan; Ding Hua; Peng Lijun; Hu Dingjin
Journal: International Journal of Genomics
 (2015)

---

## Abstract

Tea (Camellia sinensis) is a popular beverage all over the world and a number of studies have focused on the genetic uniqueness of tea and its cultivars. However, molecular mechanisms underlying these phenomena are largely undefined. In this report, based on expression data available from public databases, we performed a series of analyses to identify genes probably relevant to the uniqueness of C. sinensis and two of its cultivars (LJ43 and ZH2). Evolutionary analyses showed that the evolutionary rates of genes involved in the pathways were not significantly different among C. sinensis, C. oleifera, and C. azalea. Interestingly, a number of gene families, including genes involved in the pathways synthesizing iconic secondary metabolites of tea plant, were significantly upregulated, expressed in C. sinensis (LJ43) when compared to C. azalea, and this may partially explain its higher content of flavonoid, theanine, and caffeine. Further investigation showed that nonsynonymous mutations may partially contribute to the differences between the two cultivars of C. sinensis, such as the chlorina and higher contents of amino acids in ZH2. Genes identified as candidates are probably relevant to the uniqueness of C. sinensis and its cultivars should be good candidates for subsequent functional analyses and marker-assisted breeding.


---

## Body

## 1. Introduction

Tea (Camellia sinensis) is one of the most popular beverages in the world. It belongs to the genusCamellia (Theaceae: Ericales) and originated from East Asia [1]. The cultivation ofC. sinensis probably started more than 2,000 years ago [1]. The extensive secondary metabolites in tea leaves, including polyphenols, theanine, and volatile oils, are good for people’s health [2]. Nowadays, many cultivars ofC. sinensis, such as Longjing 43 (LJ43), Zhonghuang 1 (ZH1), and Zhonghuang 2 (ZH2), are cultivated extensively in China. In addition toC. sinensis, the genusCamellia includes many other species of great value, such asC. oleifera andC. azalea.C. oleifera has been cultivated for thousands of years in China. It is a kind of small tree with multiple trunks and branches. Its seeds can be pressed to yield edible tea oil, more than 80% of which is monounsaturated fat [3]. Unlike the aforementioned two species,C. azalea is very precious and was discovered in Guangdong Province decades ago. This is a very beautiful plant with great ornamental value.A transcriptome is all the transcripts expressed in one or a population of cells at a certain time. With the advent of next-generation sequencing technology, a great number of transcriptomes, especially those from nonmodel species, have been reported. Given their important economic and ornamental values, transcriptomes of the above threeCamellia species have been reported and some of their unique characteristics were identified. Specifically, Shi et al. described the transcriptome ofC. sinensis and identified candidate genes related to natural product pathways that are important to tea quality, such as genes involved in flavonoid, theanine, and caffeine biosynthesis pathways [4]; Wang et al. compared the biochemical and transcriptomic differences between LJ43 and ZH2 to uncover mechanisms underlying their phenotypic differences [5]; and Xia et al. characterized the transcriptome from tender shoots, young leaves, flower buds, and flowers ofC. oleiferaand detected many genes potentially related to lipid metabolism [6]. Despite this body of work, the molecular events underlying the uniqueness ofC. sinensis and its cultivars remain largely undefined.Identifying coding regions harboring mutations probably relevant to the uniqueness of a species or a cultivar has attracted the interests of many biologists for decades. Generally, nonsynonymous substitutions are harmful for their carriers and will be eliminated rapidly. However, a small number of nonsynonymous substitutions will benefit their carriers and genes harboring those mutations are termed positively selected genes. When a population extends its range or is moved by human activity to a new environment with environmental factors that are different to the original one, a series of genes may be subject to positive selection and may ultimately result in a new species. Comparing the numbers of nonsynonymous (d
N) and synonymous (d
S) substitutions per site is often used for diagnosing the extent and direction of selection on sequence evolution, with the ratio d
N
/
d
S >1, =1, and <1 denoting positive evolution, neutral evolution, and purifying evolution, respectively [7, 8]. Church et al. demonstrated that using likelihood-based variable selection models is feasible for comparing sequence pairs [9]. In fact, many studies have reported analyses identifying positively selected genes for several different species [10–12]. It seems highly likely that a number of genes may be subject to positive selection during the speciation ofC. sinensis. Moreover, during the cultivation ofC. sinensis, many important agronomical genes may be subject to artificial selection and may therefore result in a new cultivar.In this report, based on the transcriptomes reported for the three species of the genusCamellia, we identified genes potentially relevant to the uniqueness ofC. sinensis. Furthermore, candidate genes probably relevant to the divergence of LJ43 and ZH2 were also detected. Our results should be important for understanding the uniqueness ofC. sinensis and its cultivars and provide hints for subsequent breeding.


## 2. Materials and Methods

### 2.1. Data Acquisition and Filtering

Paired-end Illumina short reads generated for the floral bud transcriptome ofC. azalea (PRJNA257896) [13] and 454 reads for transcriptomes of tender shoots, leaves, flowers, and flower buds ofC. oleifera[6] generated by Xia et al. (PRJNA239933), along with 454 reads of the new shoot transcriptome generated by Wang et al. for fourC. sinensis strains (PRJNA223181) [14], were downloaded from the Sequence Read Archive at the National Center for Biotechnology Information (SRA, http://www.ncbi.nlm.nih.gov/sra/). In addition, paired-end Illumina short reads of the transcriptome of the two cultivars ofC. sinensis, LJ43 (PRJNA261659, PRJNA240661, and PRJNA79643) and ZH2 (PRJNA261659), were downloaded from SRA.Quality controls were implemented using NGS QC Toolkit_v2.3.3 with default settings [15]. Low quality bases that reside in short reads generated using Illumina and Roche 454 platforms were filtered by scripts included in this package. Single-end 454 reads were further processed using Seqclean (http://sourceforge.net/projects/seqclean/files/latest/download) to trim the vector sequences included in UniVec (ftp://ftp.ncbi.nih.gov/pub/UniVec/).


### 2.2.De Novo Assembly

Ade novo assembly based on paired-end Illumina short reads was performed using Trinity [16] with default settings forC. azalea. For the assembly of 454 reads of the other two species, the iAssembler package was used, which employs MIRA (http://sourceforge.net/projects/mira-assembler) and CAP3 [17] and can assemble large-scale ESTs into consensus sequences with significantly higher accuracy [18]. For each species, contigs shorter than 200 bp were discarded in the ensuing analyses. For thede novo transcriptome assembly for each species, TransDecoder (http://sourceforge.net/projects/transdecoder/) was used for predicting the probable open reading frames.


### 2.3. Identification of Orthologous Genes and Alignment

Tentative orthologs among the three species were predicted using a transitive Reciprocal Best Hits (RBH) approach implemented in the Ortholuge pipeline [19] with default settings except the e-value for blastn being set to 1
e
-
9. For each ortholog group, we compared each nucleotide sequence with the corresponding protein sequence predicted forC. oleifera using Genewise [20] and used a customized Perl script to extract the matched coding regions and generate the proper alignment format for the subsequent PAML [8] analyses. We excluded alignments with premature stop codons or those shorter than 30 codons after deleting the gaps.


### 2.4. Evolutionary Analyses

The evolutionary rates (K
a, K
s, and K
a
/
K
s) for each ortholog group and each separate species in the genusCamellia were calculated using CODEML under the branch-free model (model = b). To test for selection acting onC. sinensis, we used CODEML’s branch-site models (model = 2, NSsites = 2), which allowed d
N
/
d
S to vary among codon sites and across branches of the phylogeny. By setting the branch leading toC. sinensis as the foreground branch, we compared a selection model that allowed a class of codons on that branch to have d
N
/
d
S > 1 (Model A2, fix_omega = 0, omega = 1.5) with a neutral model that constrained this additional class of sites to have d
N
/
d
S = 1 (Model A1, fix_omega = 1, omega = 1). A likelihood ratio test (LRT) with χ
2 approximation was used to determine the significance of difference between the nested models. All these analyses were implemented twice with different starting omega values to check for convergence. Alignments of positive results were checked and adjusted manually and subject to analyses once again to reduce the possibility of false positives.


### 2.5. Mutations between LJ43 and ZH2 and Experimental Verification

All Illumina paired-end reads of LJ43 and ZH2 were used for identifying mutation sites potentially relevant to their divergence. For these two cultivars, all reads were mapped back to the nonredundant transcriptome using BWA (-n 0.005 -k 5) [21], and all duplicate reads were removed using the MarkDuplicates program from the software package Picard (http://broadinstitute.github.io/picard/). Reads with minimum mapping quality < 20 were removed using SAMtools [22], and then a synchronized file was generated using the program mpileup2sync in the software package PoPoolation2 [23]. The synchronized file listed allele frequencies for every population at every base in the nonredundant transcriptome in a concise format. Sites with base frequencies more than five in LJ43 but absent from ZH2 were identified as unique mutations and were compared with the predicted transcript structure to decide if they were located in the coding regions and changed the encoding amino acids.Due to the fact that RNA-Seq may produce false positives, mutation sites residing in some genes and probably related to the phenotype differentiation were selected and subject to experimental verification. Briefly, total RNAs from the leaves of LJ43 and ZH2 were extracted with Trizol Reagent Kit (Invitrogen, Madison, USA) and were reverse transcribed into cDNA using the M-MLV RTase cDNA synthesis kit (Takara, Dalian, China); sequences involved in the pathway “Porphyrin and chlorophyll metabolism” harboring nonsynonymous mutations were selected and amplified using gene-specific primers (Table1). Sequences of the amplified productions were bidirectionally sequenced using an ABI 3730 Model DNA sequencer (Shanghai Sangon, China) and submitted to Genbank.Table 1
Profiles of the validation of genes harboring nonsynonymous mutations in ZH2 compared with LJ43.
Serial number
Genbank accession number (LJ43)
Genbank accession number (ZH2)
Description
Mutation site
LJ43
ZH2
Confirmed or not
Forward primer
Reverse primer
UN010458
KT427374
KT427366
Uroporphyrinogen decarboxylase
1301
A
G
No
CGCAGATTTGAGACTGGTTG
TGGCCAATCAAGTGAAGATG
UN010458
KT427374
KT427366
Uroporphyrinogen decarboxylase
1310
A
G
Yes
CGCAGATTTGAGACTGGTTG
TGGCCAATCAAGTGAAGATG
UN010458
KT427374
KT427366
Uroporphyrinogen decarboxylase
1421
A
G
Yes
CGCAGATTTGAGACTGGTTG
TGGCCAATCAAGTGAAGATG
UN014682
KT427373
KT427365
Ferrochelatase 1
1745
T
C
No
ACTTTCTCCAAGGCTGCTC
CCAACACGGGTACTAACG
UN016015
KT427372
KT427364
Glutamyl-tRNA reductase 1
674
G
C
Yes
AAGGCCAACTCAACAGAAGC
AGTTGGGCAAGGAGTCACTG
UN020635
KT427371
KT427363
Cytochrome c oxidase assembly protein COX15
446
G
C
Yes
GATCCCTCCATCGTCATCAT
CATGCAGCAGAAGCAAAAAG
UN028022
KT427367
KT427359
Chlorophyllase 1
346
T
C
Yes
TTACAGGAGCAATAGTAGGTT
CCATAGTAGAGGTGGAAAGA
UN031492
KT427368
KT427360
Chlorophyllase 2
176
A
G
Yes
GAATGTAAACCGCCAAGT
CATCCAAACAAGCCCTTA
UN043046
KT427369
KT427361
Chlorophyll(ide) b reductase NOL
183
A
T
No
TCTTCTGTTGTCTGGCGCTT
TCGATTCCTCCTAGCAACCAA
UN047188
KT427370
KT427362
Glutamyl-tRNA synthetase
128
T
A
No
CAAAAGCAAAAGCACCCAACC
GAACCCACCACCATACTCGC
UN047188
KT427370
KT427362
Glutamyl-tRNA synthetase
93
G
C
Yes
CAAAAGCAAAAGCACCCAACC
GAACCCACCACCATACTCGC


### 2.6. Pathway Analysis

KEGG pathways were assigned to all the ortholog groups using the KOBAS software [24]. After that, the K
a values for each pathway with more than fifteen genes assigned were compared betweenC. sinensis and the other two species separately using the binomial test. Contigs with mutation sites that encoded different amino acids between LJ43 and ZH2 or were identified as positively selected genes inC. sinensis were individually subject to pathway enrichment analysis using the KOBAS software.


### 2.7. Differentially Expressed Gene Families betweenC. sinensis andC. azalea

The identification of differentially expressed genes between species is blocked by the accurate assignment of ortholog relationships, especially for nonmodel species for which whole genome sequences are unavailable. To understand the pivotal genes probably relevant to the uniqueness ofC. sinensis such as the high content of flavonoid, theanine, and caffeine, we employed a strategy used for the identification of differentially expressed gene families cross species [25]. Briefly, the proteome ofArabidopsis thaliana was downloaded from Ensembl Plants (http://plants.ensembl.org/) and clustered as gene families using CD-HIT [26]. The sequence identity threshold was set to 0.6 and a representative sequence for each gene family was selected. Then, short Illumina reads for LJ43 andC. azalea were mapped back to thede novo assembled transcriptomes ofC. sinensis andC. azaleaseparately, using bowtie [27] with default settings. After that, each sequence of the transcriptomes of the two species was uniquely mapped to the representative of each gene family, and the number of mapped short reads was accumulated if two or more sequences were mapped to the same gene family for each species. Finally, edgeR [28] was used to identify the differentially expressed gene families between the two species, and the method of Benjamini and Hochberg was used to adjust p values for multiple comparisons [29]. Changes in expression patterns of gene families including genes involved in the pathways synthesizing “flavonoid,” “theanine,” and “caffeine” were scrutinized [4].


## 2.1. Data Acquisition and Filtering

Paired-end Illumina short reads generated for the floral bud transcriptome ofC. azalea (PRJNA257896) [13] and 454 reads for transcriptomes of tender shoots, leaves, flowers, and flower buds ofC. oleifera[6] generated by Xia et al. (PRJNA239933), along with 454 reads of the new shoot transcriptome generated by Wang et al. for fourC. sinensis strains (PRJNA223181) [14], were downloaded from the Sequence Read Archive at the National Center for Biotechnology Information (SRA, http://www.ncbi.nlm.nih.gov/sra/). In addition, paired-end Illumina short reads of the transcriptome of the two cultivars ofC. sinensis, LJ43 (PRJNA261659, PRJNA240661, and PRJNA79643) and ZH2 (PRJNA261659), were downloaded from SRA.Quality controls were implemented using NGS QC Toolkit_v2.3.3 with default settings [15]. Low quality bases that reside in short reads generated using Illumina and Roche 454 platforms were filtered by scripts included in this package. Single-end 454 reads were further processed using Seqclean (http://sourceforge.net/projects/seqclean/files/latest/download) to trim the vector sequences included in UniVec (ftp://ftp.ncbi.nih.gov/pub/UniVec/).


## 2.2.De Novo Assembly

Ade novo assembly based on paired-end Illumina short reads was performed using Trinity [16] with default settings forC. azalea. For the assembly of 454 reads of the other two species, the iAssembler package was used, which employs MIRA (http://sourceforge.net/projects/mira-assembler) and CAP3 [17] and can assemble large-scale ESTs into consensus sequences with significantly higher accuracy [18]. For each species, contigs shorter than 200 bp were discarded in the ensuing analyses. For thede novo transcriptome assembly for each species, TransDecoder (http://sourceforge.net/projects/transdecoder/) was used for predicting the probable open reading frames.


## 2.3. Identification of Orthologous Genes and Alignment

Tentative orthologs among the three species were predicted using a transitive Reciprocal Best Hits (RBH) approach implemented in the Ortholuge pipeline [19] with default settings except the e-value for blastn being set to 1
e
-
9. For each ortholog group, we compared each nucleotide sequence with the corresponding protein sequence predicted forC. oleifera using Genewise [20] and used a customized Perl script to extract the matched coding regions and generate the proper alignment format for the subsequent PAML [8] analyses. We excluded alignments with premature stop codons or those shorter than 30 codons after deleting the gaps.


## 2.4. Evolutionary Analyses

The evolutionary rates (K
a, K
s, and K
a
/
K
s) for each ortholog group and each separate species in the genusCamellia were calculated using CODEML under the branch-free model (model = b). To test for selection acting onC. sinensis, we used CODEML’s branch-site models (model = 2, NSsites = 2), which allowed d
N
/
d
S to vary among codon sites and across branches of the phylogeny. By setting the branch leading toC. sinensis as the foreground branch, we compared a selection model that allowed a class of codons on that branch to have d
N
/
d
S > 1 (Model A2, fix_omega = 0, omega = 1.5) with a neutral model that constrained this additional class of sites to have d
N
/
d
S = 1 (Model A1, fix_omega = 1, omega = 1). A likelihood ratio test (LRT) with χ
2 approximation was used to determine the significance of difference between the nested models. All these analyses were implemented twice with different starting omega values to check for convergence. Alignments of positive results were checked and adjusted manually and subject to analyses once again to reduce the possibility of false positives.


## 2.5. Mutations between LJ43 and ZH2 and Experimental Verification

All Illumina paired-end reads of LJ43 and ZH2 were used for identifying mutation sites potentially relevant to their divergence. For these two cultivars, all reads were mapped back to the nonredundant transcriptome using BWA (-n 0.005 -k 5) [21], and all duplicate reads were removed using the MarkDuplicates program from the software package Picard (http://broadinstitute.github.io/picard/). Reads with minimum mapping quality < 20 were removed using SAMtools [22], and then a synchronized file was generated using the program mpileup2sync in the software package PoPoolation2 [23]. The synchronized file listed allele frequencies for every population at every base in the nonredundant transcriptome in a concise format. Sites with base frequencies more than five in LJ43 but absent from ZH2 were identified as unique mutations and were compared with the predicted transcript structure to decide if they were located in the coding regions and changed the encoding amino acids.Due to the fact that RNA-Seq may produce false positives, mutation sites residing in some genes and probably related to the phenotype differentiation were selected and subject to experimental verification. Briefly, total RNAs from the leaves of LJ43 and ZH2 were extracted with Trizol Reagent Kit (Invitrogen, Madison, USA) and were reverse transcribed into cDNA using the M-MLV RTase cDNA synthesis kit (Takara, Dalian, China); sequences involved in the pathway “Porphyrin and chlorophyll metabolism” harboring nonsynonymous mutations were selected and amplified using gene-specific primers (Table1). Sequences of the amplified productions were bidirectionally sequenced using an ABI 3730 Model DNA sequencer (Shanghai Sangon, China) and submitted to Genbank.Table 1
Profiles of the validation of genes harboring nonsynonymous mutations in ZH2 compared with LJ43.
Serial number
Genbank accession number (LJ43)
Genbank accession number (ZH2)
Description
Mutation site
LJ43
ZH2
Confirmed or not
Forward primer
Reverse primer
UN010458
KT427374
KT427366
Uroporphyrinogen decarboxylase
1301
A
G
No
CGCAGATTTGAGACTGGTTG
TGGCCAATCAAGTGAAGATG
UN010458
KT427374
KT427366
Uroporphyrinogen decarboxylase
1310
A
G
Yes
CGCAGATTTGAGACTGGTTG
TGGCCAATCAAGTGAAGATG
UN010458
KT427374
KT427366
Uroporphyrinogen decarboxylase
1421
A
G
Yes
CGCAGATTTGAGACTGGTTG
TGGCCAATCAAGTGAAGATG
UN014682
KT427373
KT427365
Ferrochelatase 1
1745
T
C
No
ACTTTCTCCAAGGCTGCTC
CCAACACGGGTACTAACG
UN016015
KT427372
KT427364
Glutamyl-tRNA reductase 1
674
G
C
Yes
AAGGCCAACTCAACAGAAGC
AGTTGGGCAAGGAGTCACTG
UN020635
KT427371
KT427363
Cytochrome c oxidase assembly protein COX15
446
G
C
Yes
GATCCCTCCATCGTCATCAT
CATGCAGCAGAAGCAAAAAG
UN028022
KT427367
KT427359
Chlorophyllase 1
346
T
C
Yes
TTACAGGAGCAATAGTAGGTT
CCATAGTAGAGGTGGAAAGA
UN031492
KT427368
KT427360
Chlorophyllase 2
176
A
G
Yes
GAATGTAAACCGCCAAGT
CATCCAAACAAGCCCTTA
UN043046
KT427369
KT427361
Chlorophyll(ide) b reductase NOL
183
A
T
No
TCTTCTGTTGTCTGGCGCTT
TCGATTCCTCCTAGCAACCAA
UN047188
KT427370
KT427362
Glutamyl-tRNA synthetase
128
T
A
No
CAAAAGCAAAAGCACCCAACC
GAACCCACCACCATACTCGC
UN047188
KT427370
KT427362
Glutamyl-tRNA synthetase
93
G
C
Yes
CAAAAGCAAAAGCACCCAACC
GAACCCACCACCATACTCGC


## 2.6. Pathway Analysis

KEGG pathways were assigned to all the ortholog groups using the KOBAS software [24]. After that, the K
a values for each pathway with more than fifteen genes assigned were compared betweenC. sinensis and the other two species separately using the binomial test. Contigs with mutation sites that encoded different amino acids between LJ43 and ZH2 or were identified as positively selected genes inC. sinensis were individually subject to pathway enrichment analysis using the KOBAS software.


## 2.7. Differentially Expressed Gene Families betweenC. sinensis andC. azalea

The identification of differentially expressed genes between species is blocked by the accurate assignment of ortholog relationships, especially for nonmodel species for which whole genome sequences are unavailable. To understand the pivotal genes probably relevant to the uniqueness ofC. sinensis such as the high content of flavonoid, theanine, and caffeine, we employed a strategy used for the identification of differentially expressed gene families cross species [25]. Briefly, the proteome ofArabidopsis thaliana was downloaded from Ensembl Plants (http://plants.ensembl.org/) and clustered as gene families using CD-HIT [26]. The sequence identity threshold was set to 0.6 and a representative sequence for each gene family was selected. Then, short Illumina reads for LJ43 andC. azalea were mapped back to thede novo assembled transcriptomes ofC. sinensis andC. azaleaseparately, using bowtie [27] with default settings. After that, each sequence of the transcriptomes of the two species was uniquely mapped to the representative of each gene family, and the number of mapped short reads was accumulated if two or more sequences were mapped to the same gene family for each species. Finally, edgeR [28] was used to identify the differentially expressed gene families between the two species, and the method of Benjamini and Hochberg was used to adjust p values for multiple comparisons [29]. Changes in expression patterns of gene families including genes involved in the pathways synthesizing “flavonoid,” “theanine,” and “caffeine” were scrutinized [4].


## 3. Results and Discussion

Withde novo assembly methods, we obtained 246,972, 103,002, and 141,099 sequences forC. azalea,C. sinensis, andC. oleifera, respectively. The number of sequences obtained forC. azalea was similar to that reported previously [13]. However, using the Newbler software, other researchers obtained 60,479 and 120,425 sequences forC. sinensis andC. oleifera, respectively [6, 14]. The significantly greater number of sequences we obtained for these two species may be because Newbler performs best for restoring full-length transcripts [30, 31], but iAssembler can identify incorrectly assembled contigs and should be more conservative [18]. Moreover, the greater number of sequences forC. azalea was much more than the other two species and may result from the different sequencing platforms and assemblers, since short reads forC. azalea were generated by using Illumina sequencing technology and the long reads for the other two species were generated by 454 pyrosequencing.Using the transitive RBH method, 8,787 one-one-one ortholog groups were identified and 4,617 groups were selected for subsequent evolutionary analyses. Distribution of alignment length for each ortholog group reveals that most groups are shorter than 600 bp. This may be because sequences used in the current analyses were generated by RNA-Seq and have a bias for shorter sequence reads (Figure1).Figure 1
Length distribution of sequence alignments for the ortholog groups amongC. azalea,C. sinensis, andC. oleifera.Using a likelihood method and the binomial test, we found that the evolutionary rates of genes involved in metabolic pathways were not significantly different among the three species (data not shown). The likelihood method identified a total of 97 sequences as positively selected genes inC. sinensis, and pathways that these sequences may participate in are shown in Supplemental Table S1 (in the Supplementary Material available online at http://dx.doi.org/10.1155/2015/527054). The pathways are mainly related to the metabolism of carbohydrates, lipids, amino acids, and some secondary metabolites, such as “glycolysis/gluconeogenesis,” “steroid biosynthesis,” “fatty acid degradation,” “arginine and proline metabolism,” “histidine metabolism,” and “butanoate metabolism.” However, genes participating in pathways synthesizing important secondary metabolites of the species, such as “flavonoid,” “theanine,” and “caffeine,” were not identified in the current analyses. These results may be because positively selected genes involved in these pathways were not included in the original dataset and the conservatism of the method in the condition of few sequences [32]. Further investigations employing the transcriptomes of more species may address this issue.Changes of expression patterns of some pivotal genes may also contribute to the uniqueness ofC. sinensis. To test this hypothesis, differentially expressed gene families betweenC. sinensis (LJ43) andC. azalea were identified. Specifically, profiles of the genes that participate in the synthesis of the three important secondary metabolites inC. sinensisare shown in Table 2. We found that most gene families including the genes that participate in the pathways synthesizing “flavonoid,” “theanine,” and “caffeine” were significantly upregulated and expressed inC. sinensis. In contrast, the significance of gene families including those genes that were downregulated inC. sinensis was not so notable, except the family including the gene “leucoanthocyanidin reductase.” Thus, we suggest that the uniqueness ofC. sinensis may result from the upregulation of some pivotal genes. However, the reliability of the results needs further investigation since these results are based on the comparison of expression patterns of gene families and employing the combined expression data from different tissues.Table 2
Profiles of the expression patterns of gene families (C. sinensis (LJ43) versus C. azalea) including genes involved in the pathways synthesizing flavonoid, theanine, and caffeine.
Secondary metabolic pathway
Gene name
Number of members of the family
Fold change (log
⁡
2)
p value
p adj.
Flavonoid biosynthesis
Phenylalanine ammonia lyase
4
1.73
0.01
0.04
Cinnamate 4-hydroxylase
1
0.82
0.21
0.40
4-Coumarate CoA ligase
3
-
0.0
2
0.98
0.99
Flavanone 3-hydroxylase
1
-
0.06
0.92
0.97
Flavonol synthase
3
1.98
0.003
0.02
Leucoanthocyanidin reductase
2
-
2.35
0.0006
0.004
Glutamine synthetase
6
1.15
0.08
0.20
Theanine biosynthesis
Gamma-glutamyl transpeptidase
1
-
0
.
20
0.76
0.87
Alanine aminotransferase
1
2.06
0.002
0.01
S-Adenosylmethionine decarboxylase
3
-
0.10
0.88
0.94
Caffeine biosynthesis
GMP synthase
1
0.65
0.32
0.53
5
′-Nucleotidase
1
1.56
0.02
0.07
Adenylosuccinate synthase
1
2.29
0.0008
0.005
AMP deaminase
1
0.002
1
1
S-Adenosylmethionine synthase
4
0.76
0.25
0.44Mutation analyses found polymorphisms present between LJ43 and ZH2 for more than 10,000 sites. Further investigations showed that 3,655 mutations were located in the coding regions and 2,021 of them were nonsynonymous mutations. Pathway analysis showed that genes harboring these nonsynonymous mutations were involved in a number of pathways. In particular, the pathways “alanine, aspartate, and glutamate metabolism,” “Porphyrin and chlorophyll metabolism,” “glycine, serine, and threonine metabolism,” “valine, leucine, and isoleucine degradation,” and “flavonoid biosynthesis” were well represented (Supplemental Table S2). Gene ontology analysis showed that genes with those nonsynonymous mutations were significantly enriched in GO terms “chloroplast part” and “chloroplast stroma” (Supplemental Table S2). To validate the mutations identified using high-throughput analyses, the sequences involved in the pathway “Porphyrin and chlorophyll metabolism” were selected and subjected to experimental verification. Bidirectional sequencing confirmed seven of the eleven nonsynonymous sites residing in these genes (Table1). Biochemical analyses revealed that contents of free amino acids and flavonoid are different in the yellow-leaf tea cultivar ZH2 and normal cultivar LJ43 [5]. Our study suggests that nonsynonymous mutations residing in the coding regions of some genes may also take part in the formation of differences between the two cultivars, in addition to the differential expression of some other genes [5].



---

*Source: 527054-2015-10-12.xml*