Transcriptome analyses of a Chinese hazelnut species Corylus mandshurica

Background Corylus was renowned for its production of hazelnut and taxol. To understand the local adaptation of Chinese species and speed up breeding efforts in China, we analyzed the leaf transcriptome of Corylus mandshurica, which had a high tolerance to fungal infections and cold. Results A total of 12,255,030 clean pair-end reads were generated and then assembled into 37,846 Expressed Sequence Tag (EST) sequences. During functional annotation, 26,565 ESTs were annotated with Gene Ontology (GO) terms using Blast2go and 11,056 ESTs were grouped into the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using KEGG Automatic Annotation Server (KAAS). We identified 45 ESTs that were homologous to enzymes and transcription factors responsible for taxol synthesis. The most differentiated orthologs between C. mandshurica and a European congener, C. avellana, were enriched in stress tolerance to fungal resistance and cold. Conclusions In this study, we detected a set of genes related to taxol synthesis in a taxol-producing angiosperm species for the first time and found a close relationship between most differentiated genes and different adaptations to fungal infection and cold in C. mandshurica and C. avellana. These findings provided tools to improve our understanding of local adaptation, genetic breeding and taxol production in hazelnut.


Background
Corylus is an important genus, both economically and ecologically, in China. There is currently more than 4 million acres of natural hazel groves in northeastern China alone [1]. Its nuts, rich in unsaturated fat and vitamins, are widely consumed. Its leaves are used by local farmers to feed domestic silkworm [2]. Its stocks have been successfully used for grafting Castanea henryi to improve chestnut production and quality [3] and have also been shown to be an ideal substitute for logs of Carpinus cordata in Ganoderma culture [4]. A part from its clear economic importance, Corylus plays an important role in soil and water conservation owing to its strong root system and contributes to the sustainability of forests in this region [2].
Corylus species are also important sources of taxol (also named as Paclitaxel), which is an effective yet relatively expensive medicine for treatment of breast, ovarian and lung cancer [5][6][7]. Taxol was originally isolated from the bark of Pacific yew [8] and then later found to be present in the yew genus Taxus [9]. It was initially believed to occur only in gymnosperms, but was recently identified in leaves and fruits of a hazelnut species (C. avellana) [10]. Further studies validated this finding by showing that in vitro hazel cell cultures produce taxol and taxanes, indicating that it is not exclusively produced by symbiotic fungi [11][12][13]. Taxol was recently discovered in another hazelnut species, C. mandshurica (synonym to C. sieboldiana) [14]. However, except for the Corylus species as well as a few other species like Magujreothamnus speciosus, Morinda citrifolia, Justicia gendarussa and Yunnanopilia longistaminata [15,16], few angiosperm species have been reported to contain taxol or its derivatives. Interests in taxol production from hazel trees, especially from its leaves, have grown rapidly with the aim of conserving endangered yew species [17].
C. mandshurica is widely distributed in northeastern China and its nuts are characterized by a thin husk and high shelling percentage [18]. The nuts from this species are of higher quality in flavor and taste, and therefore command a higher price than the nuts from C. avellana.
Moreover, C. mandshurica is highly resistant to Eastern Filbert Blight [19], a fungus that causes seriously damage to most commercially grown cultivars of C. avellana in the US [20], and has exceptional cold resistance; it is able to survive a frigid winter of -48°C [21,22]. All these traits make it a very desirable target for developing improved selections and breeding material [18,23]. Interspecific crossing and breeding experiments have been attempted between C. mandshurica and the commercial species C. avellana [1,[22][23][24][25][26]. Molecular breeding aided by microsatellite marking has also been reported [27,28].
Next generation sequencing is a quick and cost-effective method for surveying the complete coding sequence of a genome. Much progress has been made in obtaining longer sequence reads, and many tools and algorithms have been developed to allow assembly of short reads. Despite the ever-increasing sequencing data, the Expressed Sequence Tags (ESTs) from C. avellana have only recently been released [29] and remain the only available large-scale sequencing data for the Corylus genus. In this study, we sequenced the leaf transcriptome of C. mandshurica native to China. Our aims were (1) to explore how homologous genes of two hazel species have differentiated to give their contrasting adaptations, and (2) to identify the possible genetic basis of taxol production in the genus Corylus by transcriptome assembly and gene annotation.

Sequencing and assembly
After strict quality control, 12,255,030 clean pair-end reads were assembled into 37, 846 ESTs longer than 200 bp using Trinity [30]. The contig N50 was 799 bp and 8,328 ESTs had longer sequences. A total of 37,652 coding DNA sequences (CDSs) were predicted to have an average length of 431 bp using Orfpredictor [31]. Comparison of our assembly with the Jefferson transcriptome assembly on C. avellana was shown in Table 1. It could be seen that the contig N50 of our assembly was slightly lower, which was partly due to the large increase in the number of assembled sequences. It was apparent from EST length distribution ( Figure 1) that our assembly had more sequences at all length intervals beyond 160 bp. The assembled EST sequences of C. mandshurica in fasta format were available in Additional file 1.
Given the recently released genome of Betula nana [32], we used BLAT [33] to map our transcriptome assembly against this genome that currently consisted of 551,915 contigs. We found that 32,078 ESTs mapped to 32,849 contigs. In comparison, 25,073 ESTs out of the Jefferson transcriptome assembly mapped to 25,841 contigs, with 19,908 contigs shared between the two hazelnut transcriptome assemblies (Figure 2). The ESTs that mapped to unique contigs might represent different genes specifically expressed in each species or different fragments of the same genes due to the fragmentary nature of the current Betula genome and the limited sequencing depth of the transcriptomes. Thus, our transcriptome analysis revealed many novel EST sequences for the Corylus genus that could not be identified from the Jefferson transcriptome assembly and helped locate the genomic locus for each EST, which had important implications for the development of further breeding markers of the Corylus species.

Functional annotation
To functionally classify the assembled ESTs, a homologybased approach was adopted in transcriptome annotation. A total of 30,536 ESTs gave hits on performing BLASTX searches [34] in the NCBI non-redundant protein database using an E-value cutoff of 1e-5, accounting for 80.7% of all assembled sequences. When sorting the top blast hits by species, Vitis vinifera was ranked first with 10,321 top blast hits, followed by Populus trichocarpa and Ricinus communis with 5,537 and 5,155 top blast hits, respectively ( Figure 3). In addition, 26,565 ESTs were annotated with Gene Ontology (GO) terms using Blast2go [35] and 11,056 ESTs were annotated into Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways with KEGG Automatic Annotation Server (KAAS) [36] using the Single-directional Best Hit (SBH) method.
Identification of highly differentiated genes in the transcriptomes of C. mandshurica and C. avellana Since C. mandshurica and C. avellana were closely related species but with contrasting adaptations, our first goal was to identify which genes were highly differentiated. Using the available EST sequences for these two species, we performed a reciprocal blast to obtain best hit orthologs and compared both the sequence identities and presence of Insertion/Deletion (INDEL) according to BLASTN outputs. Since 98.7% of orthologs showed a sequence identity higher than 90% (Figure 4), we set a sequence identity of 90% as the low threshold in ortholog validation to exclude the presence of distantly related homologs. Because we were interested in orthologs with relatively great divergence between C. mandshurica and C. avellana, we took orthologs with the low sequence identity (less than 97%), which account for 10.4% of all orthologs, as the highly differentiated genes. Furthermore, INDEL might cause different reading frames in coding regions of the two sets of orthologous ESTs [37]. It might also cause mRNA secondary structure change [38] in the coding and noncoding regions with alternative roles in transcriptional polyadenylation site selection [39], pre-mRNA splicing [40], mRNA stability, translation efficiency and protein folding [41,42]. Therefore, INDEL was also used as an indicator for sequence divergence. Thus, orthologs with gaps in the BLAST alignment were taken as another set of highly differentiated genes. Next, we performed separate GO enrichment analyses on these two types of differentiated orthologs using the WEGO web server [43]. GO    The distribution of sequence identity in blast regions for all orthologs. 98.7% orthologs have sequence identity no less than 90%; 87.9% orthologs have sequence identity no less than 97%. To exclude distantly homologous sequences, orthologs with sequence identity less than 90% are discarded in GO enrichment analysis.  types of highly differentiated orthologs were available in Additional files 2, 3, 4). According to the GO enrichment analyses ( Figures 5  and 6), orthologs from most statistically significant GO terms were conserved in the respective sequences as they generally contained a low percentage of orthologs with a sequence identity lower than 97% or orthologs with INDEL. The conserved GO categories comprised GO terms in cellular process, developmental process and metabolic process (except hormone metabolic process) in the biological process domain for both species. All these processes were essential for plant survival. The divergent GO categories comprised immune, pollination and response to stress in the biological process domain for the GO enrichment of orthologs with sequence identities lower than 97%. The orthologs with INDEL were enriched in hormone-related or various stimuli-related GO terms in the biological process domain (including hormone metabolic process and responses to various stimuli, especially response to stress). These findings suggested that C. mandshurica and C. avellana had become genetically differentiated whilst adapting to their different habitats.
Stress response genes were more prone to both sequence substitution and insertion/deletion, with occurrences of 25.6% and 25.2% among all differentiated ESTs, respectively. A close examination of GO terms under response to stress (Figures 7 and 8) revealed that three main categories displayed increased sequence divergence, including genes participated in defenses to bacteria and fungi, genes involved in cold tolerance and genes related to salt/drought/water stress. As C. mandshurica has better adapted to fungal infection and cold stress than C. avellana, further study of the highly divergent genes in C. mandshurica could identify the key genes responsible for the resistance to fungal infection and cold pressure. However, because orthologs were not necessarily one to one match between two species when gene duplication occurred after speciation, the identified orthologs could be either true alleles or different copies of the same family in the genomes. Under the latter scenario, differential expressions of the genes at both time and space should be carefully examined, which might also represent one of the adaptation mechanisms in this species.

Genes responsible for taxol synthesis
According to KEGG annotation, 29 ESTs were found to be involved in the terpene synthesis pathway. These included genes involved in isopentenyl-PP (IPP) synthesis in both the mevalonate and MEP/DOXP pathways and genes responsible for geranyl-PP and geranyl-geranyl-PP (GGPP) synthesis. The committing step for taxol production was the conversion of GGPP to taxa-4(5)-11(12)-diene in the diterpenoid biosynthesis pathway; however, genes involved in this reaction, as well as the following processes, were absent from our KEGG annotation. This was also encountered in the KEGG annotation of C. avellana transcriptome. Nonetheless, 31 ESTs (Table 2) were found to be homologous to the prototype genes participating in taxol synthesis in yew species, with sequence identities ranging from 23.93% to 50.32%. This was similar to the sequence identities of 40%~44.1% reported in some taxol-producing fungi [44] and was close to the maximal sequence identity of around 40%~49.3% found between these genes and the available proteins from other plant species in the NCBI non redundant protein database (Table 3). In addition, 6 ESTs were found to be homologous to WRKY, and 8 ESTs homologous to JAMYC. These two transcriptional factors had been reported to induce taxol synthesis [45,46].
Overall, our study reported for the first time large-scale identification of genes involved in the terpenoid pathway in Corylus, which would facilitate understanding of taxol synthesis in angiosperms, although further experiments were required to clarify the roles of these genes in such processes. On the other hand, it should be noted that not all sequences of the genes related to taxol synthesis were revealed by the present transcriptome analyses because of difficulties in normalizing all cDNAs before sequencing when the level of leaf mRNA expression in the taxol synthesis pathway was very low. In addition, some taxadiene synthase genes might only be expressed in response to external stimuli, such as naturally occurring fungal infection or artificial chemical induction [47]. Such genes would not be detected by the present approach. Since the family of terpene synthases were highly diversified across plants [48], it would be interesting to investigate the reasons why taxol production was shared by these special gymnosperm and angiosperm plants. Horizontal gene transfer was a likely cause of such convergent evolutions via symbiotic organisms. For example, three genes from different taxol-producing fungi (two from Ozonium sp. BT2 and one from Cladosporium cladosporioides) isolated from the inner tree barks [49][50][51] had been shown high sequence   (Table 3). Undoubtedly, these unsolved questions merited further study, especially from genomic scanning and experimental tests.

Conclusions
In the present study, the transcriptome of C. mandshurica was de novo assembled with Trinity and functionally annotated with Blast2go and KAAS. We found that highly differentiated genes between C. mandshurica and C. avellana correlated with local adaptation of the two species. In addition, a set of genes that might contribute to taxol production were identified and genetic mechanisms for taxol synthesis in distantly related plants were discussed. Thus, our study broadened the available transcriptome resources for Corylus, and provided meaningful information for researchers interested in taxol synthesis and high tolerance of C. mandshurica to fungal infection and cold stress.

Sequencing and assembly
Total RNA was extracted from leaves of C. mandshurica according to the CTAB protocol. The integrity of RNA was detected on an Agilent 2100 Bioanalyzer. The initial 20 μg of total RNA was purified using polydT conjugated beads to extract polyA-tagged mRNA, which was subsequently cleaved into~200 bp fragments by treatment with divalent cations at 75°C. The first strand cDNA synthesis was carried out using reverse transcriptase (Invitrogen) with random hexamer primers, and the second strand using RNase H (Invitrogen)  and DNA polymerase I (New England BioLabs). Sequencing was performed on an Illumina Genome Analyzer II. After removal of adapter sequences, raw reads were filtered according to stringent criteria [52]. The clean reads generated were used for all subsequent analyses. Trinity was used to assemble the paired-end short reads into contigs.

Functional annotation
The EST sequences were searched against the NCBI non redundant protein database using BLASTX with an E-value cutoff of 1e-5. The blast output in XML format was then annotated by Blast2go using default parameters. Kyoto Encyclopedia of Genes and Genomes (KEGG) was a universally acknowledged database for delineating networks of macromolecular interaction within cells. Pathway annotation was conducted using the KEGG Automatic Annotation Server (KAAS) web server and Single-directional Best Hit (SBH) method against representative sets for eukaryotes. GO enrichment was analyzed with WEGO.

Ortholog identification and comparison
Bi-directional BLASTN searches were performed for the transcriptomes of C. mandshurica and C. avellana. The reciprocal best blast hits were considered as orthologs. Orthologs with sequence identity lower than 90% were discarded in further GO analyses in order to exclude distant homologs due to the incomplete and fragmentary nature of transcriptomes. Two types of sequence variations were studied in GO enrichment analyses. One type focuses on orthologs with relatively low sequence identity, which includes 10.4% of all orthologs with sequence identity less than 97%. The other focuses the presence of gaps in local alignments of orthologs as shown in BLASTN outputs. GO terms of all orthologs and these two types of orthologs were extracted from Blast2go outputs. GO enrichment analyses were carried out on WEGO server. GO terms with p-value of Pearson Chi-square test below 0.05 was considered statistically significant.

The identification of genes involved in taxol synthesis
Genes related to taxol syntheses were identified by extensively parsing gene descriptions in the XML-formatted BLASTX output using key words of all the corresponding enzymes. The potential genes were further manually verified.
In order to compare these sequences with homologous genes in other species, we used the prototype genes responsible for taxol synthesis in yew as query sequences to search against the NCBI non redundant protein database using BLASTP. The top protein hits from fungi and plants were extracted.

DATA Availability
Reads are deposited at NCBI SRA (SRR857924).

Additional files
Additional file 1: The assembled EST sequences of C. mandshurica in fasta format.
Additional file 2: GO terms for orthologs with sequence identity lower than 97.