A comprehensive assessment of the transcriptome of cork oak (Quercus suber) through EST sequencing

Sequencing

Assembly

A stringent assembly pipeline was implemented and is summarized in Figure 1. The assembly methodology is described in the Materials and Methods section, consisting of two stages: first each library was assembled individually, and secondly all assembled libraries were further assembled (assembly of assemblies). The choice of this two-step protocol lied in the asynchronous nature of the libraries being sequenced, and the need to deal with future libraries that are expected to be generated for other conditions and stress types. The choice of parameters in our protocol maximized the number of contigs and their length (in MIRA -‐AL:egp = no:mrs = 85 reduces gap penalties and permits longer matches; -‐AS:mrpc = 1 allows for single read contings, thus increasing the number of contigs), was extensively validated, and is described in greater detail in a companion paper (in preparation). We opted for de novo assembly, as the lack of a closely related species with a completely sequenced genome resulted in poor assembly (not shown). The assembly statistics for each library are shown in Table 2. A total of 577,852 putative unigenes was achieved, including 501,257 contigs and 76,122 singlets. Each library produced from 8,442 up to 50,522 putative unigenes. These were all subjected to one additional assembly step (see Material and Methods section), which reduced the number of putative unigenes to approximately 159,298 unigenes. The final unigene length distribution is shown in Figure 2A. An average unigene length of 148.5 bp was found, which is smaller than those obtained in another oak using a combination the same sequencing platform with Sanger sequencing [15, 16] (see Table 3). A BlastP of all the unigenes the NR database finds Plant best hits in 97.3% of the cases, with the remaining being hits to other species that are likely contaminations not removed by our pipeline. A plot with the species distribution of these non-plant species is found on CorkOakDB.org.

Coverage and depth

The large number of libraries used, together with the choice of a two-step assembly, resulted in a high redundancy. Most of the nearly 5 million filtered ESTs were assembled into a large number of unigenes (~159 K). We obtained an average coverage depth of 3.9 (number of times each nucleotide was sequenced), with a maximum depth of 429 (25% percentile = 1; 75% percentile = 5). This is higher than other recent tree EST projects using the same sequencing platform (e.g. [22]), likely due to the extensive number of libraries sequenced in this project, prepared from multiple tissues, developmental stages and stress conditions. After the two rounds of assembly, 61,687 high quality reads remained unassembled and were treated as singletons. Thus, 65% of our unigenes derive from contigs, higher than other recent comparable projects (see Table nine in [15]).

In the absence of a complete genome sequence, it is impossible to know the true coverage of the cork oak gene space offered by this project. However, when we queried the proteomes of Arabidopsis thaliana and Populus trichocarpa using BLASTp to determine the potential number of unique genes detected, using a cut off of e < 10^-5, we found that 65% of cork oak unigenes hit 23,482 out of 27,379 predicted proteins in A. thaliana (85%), and 30,318 out of 45,555 in P. trichocarpa (67%) [23]. These numbers represent a rough estimate of the upper (85%) and lower (67%) boundaries one can expect from the Q. suber transcriptome coverage. This figure doesn’t change significantly if we use a more lenient cut off of e < 10^-2, where we hit 24,093 (79%) and 30,719 (67%), respectively. A high degree of redundancy in our unigenes is suggested, as multiple unigenes hit the same target genes in either species. The remaining 55,921 unigenes cannot find any hit in either A. thaliana or P. trichocarpa, representing about 35% of the cork oak transcriptome. These include small unigenes that would not achieve significance in BLASTp comparisons (see Figure 2A), as well as potential novel genes not present in these two genomes. This number could be eventually overestimated, if we consider some under-assembly in our libraries.

We performed a serial clustering at increasing levels of identity in order to evaluate the degree of redundancy in our assembly (Figure 2C). We found that at the protein level, there was a sharp decrease in the number of clusters at 95% identity, indicating that approximately 8000 predicted peptides show a high identity between each other, comparable to that found in other oak species [15]. This could indicate a recent event of polyploidization giving rise to many highly similar genes. Alternatively, and probably most likely, this could be accounted by the high genetic diversity among the multiple unrelated trees used to prepare the libraries [9]. Sequencing errors not fully resolved due to the relatively low coverage of many unigenes could also be responsible for this result. In the first scenario our decision to filter off redundancies at the cDNA level at 98% could have been excessive, leading to the underestimation of the predicted number of unigenes. In contrast, the second and third scenarios would suggest that 95% is insufficient and we are overestimating the number of unigenes that may be closer to 151,000. We do not have enough data to favour any of these scenarios, in particular because all three may co-exist. We have thus chosen the 98% cDNA clustering as a conservative parameter that we hope does not over-cluster paralogues. With future data accumulation, it will be easier to fuse unigenes than to resolve incorrectly clustered paralogues.

Functional annotation

We mapped the cork oak unigenes to the functional classes defined in Gene Ontology (GO) [24]. We had 73,766 sequences mapped to at least one GO term and the unigenes covered a total of 2,273 different GO terms. Each unigene mapped to 3.66 terms on average. The vast majority of terms is present at low frequency, with a few functional classes dominating. The Biological process “Metabolism” was the most frequent, with other metabolic categories in the top five categories - metabolism related categories cover 68% of the terms assigned (Figure 3). Consistently, enzyme functions dominate the Molecular Functions (“Catalytic activity”, “Transferase activity”, “Hydrolase activity”) (Figure 3). These are in contrast with the combined ESTs of two other oaks, Q. petraea and Q. robur, where the classes Transport (Biological Process) and Nucleotide Binding (Molecular Function) dominate [15]. Note, however, that this difference may simply lie in the fact that in that study non-normalized libraries were used, resulting in under-representation of lowly expressed genes. Furthermore, this difference may also lie in the fact that in that study, nuclear and organelle transcriptomes were, to the best of our knowledge, assembled together, while we removed both chloroplast and mitochondrial sequences from our assembly. This is supported by the observation that in the GO Cellular Component classification, the “Plastid” class is the most frequent in the Q. petraea/Q. robur ESTs, while in the cork oak, intracellular classes dominate (“Cell”, “Intracellular”, “Cytoplasm”, etc.) (Figure 3).

We used a simple and conservative scheme for gene naming of the cork oak unigenes. Besides its accession number (see below for details), we gave it an unigene name based on its similarity to proteins in A. thaliana and P. trichocarpa (Table 4). We observed that for nearly 40% of the unigenes we could not assign a clear annotation at cut off of e < 10^-5 (Figure 4), consistent with the number of unigenes that are not similar to any gene in other model plants. Conversely, we could identify conserved domains in 44% of the unigenes, and could establish clear homology relationships to an additional 16% of the unigenes, in a total of 60% unigenes with clear functional assignments in GO.

Method		Assignment
BDBH		Ortholog
BLASTp search
Alignment length	identity
> 85%	> 35%	High confidence
> 70%	> 25%	Homolog
< 70%	> 30%	Conserved domain
< 70%	< 30%	Low confidence

If a gene is bi-directional best hit (BDBH) of X in A. thaliana (or P. trichocarpa), we term it ortholog of X; if it is similar to X in A. thaliana (or P. trichocarpa) using BLASTp and it aligns in 85% of its length with more than 35% identity, we term it a High confidence X in Q. suber, etc.

We were able to map Interpro domains to 108,341 unigenes (68%). Nearly half of the domains were widespread in evolution, being present in both Eukaryota and Bacteria (Figure 5). The other half was dominated by general Eukaryotic domains and less than 10% of the domains were plant specific. These results are comparable to those reported for the complete genomes of A. thaliana, P. trichocarpa and P. persica genomes, as well as to those of the transcriptomes of the closely related Quercus robur and Castanea mollissima which are also depicted in Figure 5.

Evolution

We compared the gene content of the cork oak, as estimated by our EST sequencing project, with that of 31 completely sequenced plant genomes. We used BLASTp at e < 10^-5 and also at the permissive cut off of e < 10^-2 to determine how many predicted proteins in those species are similar to at least one cork oak unigene. The results of this analysis are shown in Figure 6, indicating a broad concordance with the generic taxonomic/evolutionary distance of the species. This result does not change when we use a more permissive cut off of e < 10^-2 (not shown).

We compared the unigenes derived from the cork oak with those of the red oak (Q. rubra), the pedunculate Oak (Q. robur - also known as English or French oak) and the Chinese chestnut (Castanea mollissima). For this comparison, the data from the Fagaceae Genome Web was used, for Q. rubra and C. mollissima which include multiple tissues also sequenced using the 454 pyrosequencing platform (http://www.fagaceae.org/node/87455 and http://www.fagaceae.org/node/181796/, respectively), and the data for Q. robur, which included 454 and Illumina generated sequences, and was obtained from http://www.ufz.de/trophinoak/index.php?de=31205[16, 26]. We used our own assembly pipeline on these sequences to ensure that no additional differences were introduced on methodological grounds. The comparison is shown in Figure 7. The total number of distinct unigenes is higher in the cork oak project, probably reflecting the higher number of tissues and conditions sampled in our libraries, as well as incomplete assembly due to library biases and genetic heterogeneity of the samples. We verified that between 77% and 82% of the unigenes from those species are similar to at least one unigene in the cork oak, as expected from evolutionarily close species. The remaining 18% - 23% of the unigenes of the red and english oaks and chestnut tree are likely species-specific, but may also be partially accounted by an incomplete coverage of the Q. suber. The large number of cork oak unigenes that does not find a hit in the other transcriptomes (30% - 44% at e < 10^-5) does however suggest that, most likely, this is not a major factor. This cork-oak-specific set represents a mixture of small reads that fail to attain statistical significance (e.g. from incomplete assembly), as well as a putative set of cork oak-specific genes. Note that when we compare Q. suber with a completely sequenced genome of the Prunus persica, 94% of the P. persica genes find a hit in Q. suber, further suggesting that incomplete coverage of the gene space was probably not a major problem of our project.

Database and interface

To support the assembly and annotation pipeline we have a data warehouse system that records the data and metadata associated with each step of the pipeline. This is described in a companion paper (in preparation). From this warehouse we generated a public portal as a community resource for cork oak genomics, which is found at http://www.corkoakdb.org. The assembled genes, the proteins they encode, and the functional annotations are made accessible through a web interface, partially shown in Figure 8. The gene view features sequence data, cDNA and protein, as well as plots of base-by-base coverage information for the unigene. Users are shown pre-computed phylogenetic profiles against other plants according to two distinct methods, the bi-directional best BLAST hit and the inparanoid, two standard methods to identify orthologs and paralogues [27]. The gene view further includes functional annotations, namely GO annotations, Interpro domain assignments, KEGG pathways and best BLAST hits against general and plant-specific databases. Genes of interest can be discovered by searching specific fields or by running a nucleotide or protein BLAST search against the Cork Oak database.

Samples, collection and preparation

Within this initiative, in order to guarantee high transcript coverage and to increase gene diversity, total RNA was isolated from Quercus suber biological samples obtained from different organs and tissues at varying developmental stages (roots, leaves, buds, flowers, fruits, phellogen, vascular tissue, good and bad quality cork), as well as from plants that had been exposed to infection with Phytophthora cinnamomi, symbiosis with Pisolithus tinctorius mycorrhizal fungus and different abiotic stresses (cold, heat, drought, salinity and oxidative stress). Furthermore, total RNA was also isolated, at two distinct dates (May and September), from annual shoots of 30 years old Quercus suber x cerris hybrid trees that either produce or don’t produce cork, in order to cover different developmental stages of the phellogen meristem. No approval or licenses were required for sample collection. In each library, plant material from half-siblings (e.g. abiotic and biotic stress libraries) or from several unrelated trees was used. All the plant material used was from Portuguese trees except for those trees used to detect polymorphism, which were from different Mediterranean countries [28]. The detailed conditions applied in each situation are described in http://www.corkoakdb.org/libraries. The full set of libraries is described in Table 1.

cDNA preparation, library normalization and pyrosequencing

Total RNA from each tissue/condition was used as the source of starting material for cDNA synthesis and production of normalized cDNA libraries intended for 454 sequencing. Briefly, the total RNA quality was verified on Agilent 2100 Bioanalyzer with the RNA 6000 Pico kit (Agilent Technologies, Waldbronn, Germany) and the quantity assessed by fluorimetry with the Quant-iT RiboGreen RNA kit (Invitrogen, CA, USA). A fraction of 1–2 μg of total RNA was used for cDNA synthesis with the MINT cDNA synthesis kit (Evrogen, Moscow, Russia), a strategy based on the SMART double-stranded cDNA synthesis methodology using a modified template-switching approach that allows the introduction of known adapter sequences to both ends of the first-strand cDNA. Amplified cDNA was then normalized with TRIMMER cDNA Normalization kit (Evrogen, Moscow, Russia) using the Duplex-Specific Nuclease-technology [20, 29].

Normalized cDNA was quantified by fluorescence and sequenced in 454 GS FLX Titanium according to the standard manufacturer’s instructions (Roche-454 Life Sciences, Brandford, CT, USA) at Biocant (Cantanhede, Portugal).

Sequence processing and assembly

The implemented sequence analysis strategy included an initial pre-processing stage, performed on each library, where contaminant, low quality, redundant and repeat-full sequences were removed and each library assembled. This was followed by a multilibrary assembly (described below, and summarized in Figure 1). Initially, each read, respective quality scores and ancillary information, were extracted from the sequencing machine output (.sff), using open source software sff_extract (http://bioinf.comav.upv.es/sff_extract/). Reads of each sample were selected using a Python pipeline that screens the reads for primer sequences, classifying them by sample origin and allocating them in different files. For each sample we generated a file with the sequences (.fasta) and the corresponding file with the quality scores (.qual). At this stage we removed adaptors and reads smaller than 40 bp. Thereafter, artificial duplicates associated with pyrosequencing were removed using cd-hit-454 [30] at a threshold of 98%, and Seq-trim [31] was used to remove small sequences (length < 100 bp) or sequences with low quality (QV > 20, quality window = 10), as well as poly-A or poly-T tails, and adaptors.

In the following step, contaminant sequences were removed. For this, a database of possible types of contaminants was prepared (ContaminantsDB - see supplementary material for details) and queried with the Q. suber reads using BLASTn (5, -E 3 -e 1e-09 -q -5 -b 1 -G 3). Reads that found a match in this database, were subsequently blasted against a database of plant proteins (PlantDB - see supplementary material for details) using the same parameters as before. If the hit (match) e-value in ContaminantsDB was smaller than hit (match) e-value in Plant DB, the read was considered as a contaminant and removed from the pipeline. The remaining reads continued in the pipeline to be screened for repetitive elements, using the program RepeatMasker 3.2.9 (http://www.repeatmasker.org) against PlantRepeatsDB [32]. Whenever sequences were masked in more than 90% of their length they were discarded.

The final step of the preprocessing stage was the classification of all the trimmed reads into potential mitochondrial, chloroplastidial or nuclear sequences. For this, a BLASTn (-e = 0.001) was first performed against a database containing coding region sequences from complete plant mitochondrial genomes (from Arabidopsis thaliana, Medicago truncatula and Populus tricocharpa). The sequences that presented a hit were considered potential mitochondrial sequences and were kept in a FASTA file reserved for this organelle sequences. A similar process was then applied against a database of coding region sequences of plant complete plastidial genomes (same organisms).

Assembly

We chose MIRA 3.2.0 [33] to assemble the resulting sequences, as this has been shown to have higher coverage than other assemblers [34]. For each library, we obtained contigs and singletons with the following parameters: --job = denovo, est, accurate, 454; --GE:not = 20; --SK:not = 20; 454_SETTINGS -LR:mxti = no, -CL:qc = no:cpat = no:mbc = yes, --AL:egp = no:mrs = 85, -OUT:sssip = yes, -AS:mrpc = 1. Following this step, all the contigs and singlets resulting from the assembly of each library were then clustered to remove redundancy using CD-HiT [35], and the resulting non-redundant sequence collection was re-assembled using the same parameters as before. The resulting sequences were considered to be Unigenes, and at this point they were given an unigene accession number. Libraries L20 and L21 were not used in the analysis presented in this manuscript, but are available in the full assembly on the CorkOakDB.

Protein prediction

In order to be able to translate the nucleotide sequences to protein sequences, the pipeline first performs a Blast search (blastx) against a RNA database [36], to remove non-protein coding unigenes. It then queries all Viridiplantae protein sequences existing in the Uniprot database [37]. The program Prot4EST [38] then takes the outputs of these BLAST searches and translates the sequences into putative peptide sequences. Those unigenes without significant hits are translated using the program ESTscan [39], and for the remaining untranslated sequences, the longest ORF of the 6 frames is selected.

Sequence naming

In order to assign names to the genes/proteins found, putative peptides were used to query, using BLASTp at a cut off of e < 10^-5, a database of Uniprot sequences from A. thaliana and P. tricocharpa. Whenever a putative peptide does not have a hit, it is considered “Predicted hypothetical protein”. If a similar hit is detected, then the protein name is assigned to the putative peptide in Q. suber together with a label that describes the level of confidence of the annotation (see Table 4).

Functional annotation

In order to obtain domains and functional sites of putative peptides, an Interpro search was executed [40]. The Interpro database [41] integrates different classification methods based on amino-acid patterns and profiles, protein family fingerprints, protein sequences and structural domains, as well as functional information. The Interpro database 28.0 was downloaded and searches were run locally. Afterwards, a BLAST (BLASTp) search against non-redundant protein database was executed and results entered the program Blast2GO [42]. We used the pipeline version of the B2G called B2g4pipe, obtaining GO-terms and E.C. Numbers. The same pipeline was used to assign Interpro domains for the transcriptomes analysed in Figure 5.

Database implementation

A MySQL relational database was deployed, using the InnoDB engine to allow rollback of transactions in case of failure. This was essential, given the progressive nature of the data loading. Every EST sequence was stored in the database, and as each step of the pipeline was ran, the results were added to the corresponding tables, up to the functional annotation of assembled unigenes, as well as metadata related to the EST libraries. Some intermediate output data, such as large FASTA and XML files, were kept on the file system. The web interface is powered by a Python application built on Django (an open source web framework), HTML/CSS and Javascript. KEGG data is displayed using the KEGG SOAP API.

Accession numbers and unigene naming

Accession numbers on the corkoakDB have the following format QS_000000, for unigenes, and QS_P_000000 for putative peptides. Whenever the sequences are putative mitochondrial or potential chloroplast sequences they start with QSm or QSc, respectively.

Evolutionary analysis

Comparisons to other organisms were made using predicted proteomes obtained from the superfamily database [43] release 1.75. We used BLASTp for the comparisons, always filtering for low complexity regions and using the cut offs indicated in the text. We used the standard NCBI’s taxonomic tree as a reference for Figure 6. Red oak libraries were obtained from the Fagaceae genomics web (http://www.fagaceae.org/node/87455) and processed using our own pipeline, resulting in 38,346 predicted unigenes. We then used BLASTp with a cut off at e = 0.01 to determine how many unigenes from the cork oak were similar to at least one unigene in the red oak.

Availability of supporting data

All sequenced ESTs were submitted to the sequence read archive (http://www.ncbi.nlm.nih.gov/sra) with the accession number ERP001762, and accession name “Cork Oak”.