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-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-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-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.

Results and discussion

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.

Table 1

Comparison of transcriptome assembly and coding sequence prediction for and

	EST		CDS
	C. mandshurica	C. avellana	C. mandshurica	C. avellana
Average Length	580	532	431	377
Length Range	201 ~ 6821	80 ~ 5490	30 ~ 4890	42 ~ 4143
Numbers	37846	28255	37652	28167
N50 Length	799	961	594	651
Sequences (longer than N50)	8328	4991	8028	4945

Figure 1

The length distribution of assembled ESTs from transcriptomes of and . ESTs are counted at an interval of 80 bp, with ESTs longer than 2000 bp counted as 2000 bp. It is clear that more ESTs are assembled in the transcriptome of C. mandshurica than C. avellana at length intervals after 160 bp. (Trinity has a minimum length of 200 bp in transcriptome assembly).

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.

Figure 2

The contigs from genome mapped by transcriptomes of and . The contigs mapped uniquely by ESTs from C. mandshurica may represent novel genes specifically expressed in it or different fragments of the same genes unidentified in the transcriptome of C. avellana.

Functional annotation

To functionally classify the assembled ESTs, a homology-based 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.

Figure 3

The percentage of top blast hits by species.

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 terms for all orthologs and the two types of highly differentiated orthologs were available in Additional files 2, 3, 4).

Figure 4

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.

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.

Figure 5

GO enrichment on orthologs with sequence identity lower than 97%. Direct parent GO terms in biological process domain are displayed for simplicity. The GO terms in the biological process domain that show increase in ESTs with lower sequence identity are related to immune response and response to stress.

Figure 6

GO enrichment on orthologs with INDEL. Direct parent GO terms in biological process domain are displayed for simplicity. The GO terms in the biological process domain that show increase in ESTs with INDEL are related to hormone metabolic process and response to various stimuli (including response to stress).

Figure 7

GO terms under response to stress in GO enrichment on orthologs with low sequence identity. GO terms related to defense against bacteria and fungi (including innate immune response), tolerance to cold and heat, and salt/drought/water stress show increase in ESTs with lower sequence identity.

Figure 8

GO terms under response to stress in GO enrichment on orthologs with INDEL. GO terms related to defense against bacteria and fungi (including innate immune response), tolerance to cold and heat, and salt/drought/water stress show increase in ESTs with INDEL.

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].

Table 2

ESTs homologous to genes involved in taxol synthesis in

EST ID	Hit Protein GI	Identity (%)	Length	Description
comp37211_c1_seq1	386304248	36.51	189	10-deacetylbaccatin III-10-O-acetyl transferase, partial
comp42386_c0_seq1	28558088	28.85	104	3′-N-debenzoyl-2′-deoxytaxol N-benzoyltransferase
comp70594_c0_seq1	339521621	23.93	422	C-13 phenylpropanoid side chain CoA acyltransferase
comp70669_c0_seq1	28380187	34.95	432	taxa-4(20),11(12)-dien-5alpha-ol-O-acetyltransferase
comp68118_c0_seq1	28380187	28.44	450	taxa-4(20),11(12)-dien-5alpha-ol-O-acetyltransferase
comp53308_c0_seq1	386304662	38.04	163	taxadienol acetyl transferase, partial
comp119236_c0_seq1	53690152	45.98	87	taxadien-5-alpha-ol-O-acetyltransferase
comp68580_c0_seq1	53690152	30.06	173	taxadien-5-alpha-ol-O-acetyltransferase
comp37211_c0_seq1	53690152	32.39	142	taxadien-5-alpha-ol-O-acetyltransferase
comp83331_c0_seq1	53690152	28.84	215	taxadien-5-alpha-ol-O-acetyltransferase
comp64789_c0_seq1	53759170	42.6	446	taxadiene 5-alpha hydroxylase
comp57975_c0_seq1	386304485	50.32	155	taxadiene 5nalpha hydroxylase, partial
comp172528_c0_seq1	38201489	36.47	85	taxa-4(5),11(12)-diene synthase
comp193967_c0_seq1	15080743	46.03	63	taxadiene synthase
comp63152_c1_seq1	386304920	29.69	128	taxadiene synthase, partial
comp40035_c0_seq1	24266823	47.89	71	5-alpha-taxadienol-10-beta-hydroxylase
comp53405_c1_seq1	24266823	49.47	95	5-alpha-taxadienol-10-beta-hydroxylase
comp133851_c0_seq1	44903417	32.47	77	5-alpha-taxadienol-10-beta-hydroxylase
comp110423_c0_seq1	60459952	41.38	87	taxane 13-alpha-hydroxylase
comp69534_c0_seq1	60459952	33.79	441	taxane 13-alpha-hydroxylase
comp38773_c1_seq1	60459952	45.83	96	taxane 13-alpha-hydroxylase
comp36415_c0_seq1	60459952	34.19	427	taxane 13-alpha-hydroxylase
comp57975_c1_seq1	60459952	44.93	69	taxane 13-alpha-hydroxylase
comp104139_c0_seq1	60459952	42.17	83	taxane 13-alpha-hydroxylase
comp93979_c0_seq1	75297723	38.03	71	Taxane 14b-hydroxylase
comp61533_c1_seq1	380039801	33.57	143	taxane 14b-hydroxylase
comp143394_c0_seq1	380039801	29.17	120	taxane 14b-hydroxylase
comp74596_c0_seq1	380039801	29.51	122	taxane 14b-hydroxylase
comp84707_c0_seq1	380039801	35.22	230	taxane 14b-hydroxylase
comp36896_c0_seq1	67633430	30.84	467	taxoid 2-alpha-hydroxylase
comp192945_c0_seq1	238915468	43.75	64	taxoid 7-beta-hydroxylase
comp36946_c0_seq1	365776087	55	60	transcription factor WRKY
comp64330_c0_seq1	365776087	54.24	59	transcription factor WRKY
comp78449_c0_seq1	365776087	66.67	54	transcription factor WRKY
comp67132_c0_seq1	365776087	56.34	71	transcription factor WRKY
comp59687_c0_seq1	365776087	41.22	131	transcription factor WRKY
comp68275_c0_seq1	365776087	56.72	67	transcription factor WRKY
comp104123_c0_seq1	222355764	29.87	154	JAMYC
comp69212_c0_seq1	222355764	41.69	710	JAMYC
comp69212_c0_seq2	222355764	46.72	259	JAMYC
comp69212_c0_seq3	222355764	41.69	710	JAMYC
comp124731_c0_seq1	222355764	100	27	JAMYC
comp69971_c3_seq1	222355764	29.9	204	JAMYC
comp38183_c0_seq2	222355764	30	160	JAMYC
comp83061_c0_seq1	222355764	35.71	112	JAMYC

Protein GIs, instead of their accession numbers, are provided here for convenience in table layout. These can be queried at NCBI protein databases.

Table 3

Proteins most homologous to genes involved in taxol synthesis in species outside

Query Protein GI	Hit Protein GI	Identity (%)	Taxonomy	Hit Protein GI	Identity (%)	Taxonomy
15080743	-	-	-	62511183	48.53	Abies grandis
38201489	-	-	-	62511183	48.66	Abies grandis
386304920	-	-	-	62511183	49.3	Abies grandis
24266823	56609042	97.99	Ozonium sp. BT2	75319884	43.97	Picea sitchensis
44903417	56609042	99.2	Ozonium sp. BT2	75319884	44.17	Picea sitchensis
53759170	56609042	66.38	Ozonium sp. BT2	75319884	44.44	Picea sitchensis
60459952	56609042	63.9	Ozonium sp. BT2	75319884	45.32	Picea sitchensis
67633430	56609042	56.34	Ozonium sp. BT2	75319884	40.89	Picea sitchensis
75297723	56609042	60.72	Ozonium sp. BT2	75319884	42.5	Picea sitchensis
238915468	56609042	56.2	Ozonium sp. BT2	75319884	40	Picea sitchensis
380039801	56609042	61.12	Ozonium sp. BT2	75319884	42.92	Picea sitchensis
386304485	56609042	66.22	Ozonium sp. BT2	75319884	47.11	Picea sitchensis
28380187	62461771	98.39	fungal sp. BT2*	148906373	43.98	Picea sitchensis
28558088	62461771	60.23	fungal sp. BT2	148906373	45.62	Picea sitchensis
53690152	62461771	60.83	fungal sp. BT2	148906373	44.25	Picea sitchensis
339521621	62461771	60	fungal sp. BT2	148906373	39.83	Picea sitchensis
386304662	62461771	98.23	fungal sp. BT2	148906373	44.07	Picea sitchensis
386304248	169135276	98.45	Cladosporium cladosporioides	148906373	44.64	Picea sitchensis
365776087	-	-	-	167859869	43.18	Picea abies
222355764	-	-	-	148906957	46.99	Picea sitchensis

* Ozonium sp. BT2 and fungal sp. BT2 are the same fungus species. [49,51].

Columns 2-4 show top hit protein information from fungi; columns 5-7 show top hit protein information from plants. Methodically, protein queries are blasted against NCBI nonredundant protein database and protein hits from the two designated sources with top sequence identity are recorded. Protein GIs, instead of their accession numbers, are provided here for convenience in table layout. These can be queried at NCBI protein databases. The reason for only three identified hit proteins from fungi is possibly due to the absence of genome data for taxol-producing fungi.

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-51] had been shown high sequence identities (98.39%, 98.45% and 99.2%) to the corresponding taxol genes in yew species (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.

Methods

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.

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).

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

HM and JL analyzed the data and wrote the manuscript. ZL acquired the leaf sample. BL prepared the mRNA for sequencing. QQ provided helpful suggestion in data analysis. All authors read and approved the final manuscript.

Additional file 1

The assembled EST sequences of in fasta format.

Click here for file

Additional file 2

GO terms for orthologs with sequence identity lower than 97.

Click here for file

Additional file 3

GO terms for orthologs with INDEL.

Click here for file

Additional file 4

GO terms for all orthologs.

Click here for file