Overall statistical comparison of zoarcid cDNA libraries

The setup comprises 454 sequencing data of two cDNA libraries, one prepared from liver and heart of P. brachycephalum, and one prepared from liver of Z. viviparus[23], generated with similar techniques and protocols (c.f. Material and Methods). The average read length in the library of P. brachycephalum was about 100 bases longer than that for Z. viviparus (Table1), due to the usage of Titanium-chemistry for sequencing. Subsequently, a higher yield of longer contigs with 487 bases in average for P. brachycephalum was found, compared to 342 bases in average per contig for Z. viviparus (Figure2A). Other basic parameters like number of reads per contig (Figure2B) as well as the average coverage were comparable (Figure2C), constituting a suitable basis for alignment analyses of contigs.

Database parameter	P. brachycephalum	Z. viviparus
Total Number of reads	481,802	~400,000
Number assembled	338,993	349,102
Average of read length	321	221
Number of all contigs	65,565	53,447
Number of non-redundant contigs	64,634	53,313
Percentage of redundancy	1.42	0.25

Overall statistical analyses of the sequencing runs of Z. viviparus data of the library generated by Kristiansson et al. 2009, in comparison to the library of the Antarctic eelpout P. brachycephalum. The resulting reads were set up to contigs with similar assembly techniques. Further analyses of the contig quality could be seen in Figure2.

Using the BLAST2GO tool[26] for evaluating annotations of contigs, statistics of expect values (e-values) show qualitatively comparable profiles, But quantitatively, sequences of P. brachycephalum show higher similarity to known sequences (Figure3A) accompanied by higher coverage per hit of high-scoring segment pairs (HSP) (Figure3B). For the purposes of this study we filtered by BLAST[27] for sequences of known functions, i.e., showing known translations, resulting in 19,460 contigs for P. brachycephalum and 16,315 for Z. viviparus.

Functional coverage compared to a phylogenetically related reference genome

The overall functional repertoires of the cDNA library sets were evaluated in terms of KOG/COG categories (clusters of orthologous groups) as obtained from the local alignment tool rpstBLASTn applied on the BLAST-annotated sequences against a database of metazoan-specific orthologies (meNOGs)[28, 29]. As a scale for the comparison, coding sequences from the closely related species G. aculeatus (Figure1A, B) were evaluated accordingly. From the reference genomic data (ENSEMBLE 61) 27,628 coding sequences were extracted to estimate the functional coverage of zoarcid libraries. To compare the appearance of functional terms quantitatively, the counts were pooled into superordinated categories of KOG/COGs (Figure4). The count data of both zoarcid transcriptomes were in a comparable range but covered only a fraction of the reference genome. Strongest inter-zoarcid differences were seen for categories F (nucleotide transport and metabolism) and V (defense mechanisms), showing higher counts in P. brachycephalum. The overall coverage compared to the scales of G. aculeatus genome is estimated to be 40% to 60%, with a pronounced lesser coverage in categories U, V, D and B. Furthermore, a notable deviation is found in category N (cell motility), where both eelpouts show more counts than are present in the reference genome. Contrary, in category W (extracellular structures) Z. viviparus have more counts compared to the reference set, whereas the number of counts for P. brachycephalum nearly equals it.

To identify species-specific functional patterns in gene expression, we searched for overrepresented meNOGs in each library set. The most exclusive terms in Z. viviparus are meNOG13752 coding for complement protein, meNOG11371 encoding for the cytochrome-c-oxidase (COI) and meNOG18150 for extracellular matrix constituent lubricant protein contributing to the structural extracellular matrix integrity. In P. brachycephalum the most specific term is a b-box-type zinc finger (meNOG24554) followed by meNOG15201 encoding for a tripartite motif protein and meNOG06310 which encodes a carboxylase function. Further terms in this species were ubiquitin (meNOG05057) and the heavy chain of myosin (meNOG 07153) (data not shown). The term for myosin (meNOG07153) is noticed as specific for P. brachycephalum consistent with the fact that we included heart RNA in the library construction.

Amino acid usage

Translations of coding sequence segments were generated based on fish-specific orthologies (fiNOGs) via rpstBLASTn for both cDNA libraries to exclude effects of sequencing artefacts and frame shifts on the subsequent analyses. From the identified orthologous sequences in each zoarcid fish, matching sequence segment pairs were obtained by harvesting best BLASTp hits per model orthology (i.e. fiNOG), providing a set of shared orthologous sequences in both fish. For a more stringent pairing, i.e., excluding probably false pairs, translated segment pairs were filtered for ≥ 80% amino acid sequence identity. Accordingly, frequencies of amino acids were computed from 4,155 translated and realigned segment pairs, including non-identical amino acid positions, but excluding indels (insertions/deletions). Alterations in amino acid profiles became visible by assessing the total changes in the aligned non-synonymous codon positions (Figure5). For P. brachycephalum net losses for glutamic acid and asparagine as well as a net gain of serine were most obviously.

The observed imbalances were furthermore reanalyzed taking the ERK–proxy into consideration[30] to determine whether species-specific shifts are detectable within the observed pattern of differentially used amino acids. A significant preference for using E+R+K in Z. viviparus was detected when compared to the orthologous associate in P. brachycephalum (Wilcoxon test, p = 0.02), analysing the 4,155 translated sequences and taking only exchanged positions into account.

Codon usage

The set of 4,155 accordingly re-aligned orthologous coding sequence pairs in the transcriptomic eelpout-libraries were analyzed by means of a Within Canonical Analysis (WCA) to distinguish synonymous codon usage between the two species[31]. Comparing the composition of the triplets coding for synonymous amino acids, a shift of the P. brachycephalum transcriptome sequence segments became visible in the direction of a preferred usage of A or T in the third position of codons (Figure6). The small shift between the blue covariance ellipse relative to the red one is exclusively due to the first principal component, reflecting the G/C-content on the third codon position (GC3). Furthermore, with codon usage differences from 4,155 sequence segment pairs (filtered for > 80% identical translated positions, depicted in Figure7), more positive shifts (23 in the left panel) of codons ending with A/T (AT3, marked “>”) were detected compared to the amount of negative shifts (Figure7 in the right panel). For this analysis, non-synonymous codon positions were not taken into account. A Fisher-test for the distribution of a contingency table of AT3 vs. GC3 codons in the left and right panel of this graph reveals a p-value of 5*10^-4, indicating a significantly non-random distribution of shifts.

The overall number of aligned amino acid residue pairs in this set was 414,978, constituting an average of 100 aligned codon pairs per segment with a mean of 96 conserved residues. Comparing 402,219 synonymous codon pairs with conserved nucleotides in positions 1 and 2, P. brachycephalum shows a net GC3-loss of 1,030, re-confirming an imbalance towards a preferential A/T in the third position of the codons in this cold-adapted fish. The analysed coding sequence segments showed a mean GC-fraction of 53.196% for Z. viviparus and 53.094% for P. brachycephalum mirroring the observed pattern on the third codon position.

Functional coverage compared to a reference genome

It is noteworthy that translated coding sequences of G. aculeatus are part of the eggNOG orthologies. Therefore, we expect that a comparison in a subset of these orthologies, i.e., metazoan (meNOGs), is useful to serve as a scale for inter-comparisons of the present cDNA libraries. Nevertheless, comparisons with G. aculeatus genomic coding sequences might show an overestimation of coverage due to the unresolved multiplicities by splice variants present in the assemblies of the transcriptomes. This fact might explain the observed higher coverage of functional terms in the KOG/COG categories N (cell motility) and W (extracellular structures) of the zoarcid fish compared to the reference genome.

Screening for occurrence of library-specific terms in principle allows for detection of signatures of the environmental conditions under which the transcriptomes were obtained, and not necessarily indicates presence or absence of genes in one of the species genomes. As a result transcripts related to COI were detected in elevated proportion in Z. viviparus, which is against the expectation that the normalization step in library preparation should reduce over-abundant RNA-variants in both libraries to a comparable degree. The same holds for the Antarctic species, for which the frequency of ubiquitin annotation appears in a higher rank. Discovering library specific functional terms, which would be expected to be included in similar ranks in both transcriptomes, reflects differences in the efficiency of normalization.

Other terms, specific for the cDNA library of Z. viviparus are related to complement protein (meNOG13752) and extracellular matrix constituent protein (meNOG18150). These terms summarize functions for immune response and integrity of the tissue, comprising cell communication, growth factors and wound healing. The term for tripartite motif (TRIM) proteins, which is predominantly found in the P. brachycephalum library, summarizes sequences of a protein family associated to pathogen-recognition, regulation of the concomitant transcriptional responses, constituting an important part of the immune system[34]. These differences may indicate adaptive forces acting on the immune system or simply different life histories of the sampled specimens. Responses of the innate immune system in fish were apparently more robust and diverse than in higher vertebrates[35]. Similarly, e.g., Atlantic cod comprises a unique immune system with substantial gene losses in adaptive components without being exceptionally susceptible to disease under natural conditions[36]. Consequently specific adaptations in the immune systems in the two eelpout species due to the different habitat conditions seem possible.

Another overrepresented term in P. brachycephalum is a carboxylase, which belongs to the KOG/COG category I - lipid transport and metabolism, as listed in the meNOG tables. The lipid fraction in liver of P. brachycephalum kept at 0°C is threefold higher than in Z. viviparus at habitat temperature[37], in line with the general finding that Antarctic species have high capacities for catabolism of fatty acids[38]. Furthermore, a recent study demonstrated relationships between preferred energy fuels and transcript levels of respective genes in P. brachycephalum[39]. Although the overrepresentation of the carboxylase fits into the current picture, further genes could have been expected to occur substantiating this evidence. This mismatch may diminish when the entire genome sequence data will become available. In summary, some of the library specific functions identified in our study may point to transcriptomic and possibly genomic contrasts due to the different environmental conditions of the discrete habitats of the two zoarcid fish under study.

Amino acid usage

With the advent of more genomic data from species from various habitats analyses of temperature effects on amino acid usage profiles and the correlated GC content became possible[13, 40]. Gu and Hilser[41] analyzed intra-protein interaction energies on the level of primary and secondary structural sequence segments in thermophilic and psychrophilic proteomes, supporting the flexibility hypothesis of cold adaptation. Furthermore, local and global adaptation patterns were differentiated, depending on different protein families.

By comparing the amino acid composition of a comprehensive set of orthologous sequence pairs globally we discovered a distinct pattern of replacements, which may result in structurally important changes in the functional structure of proteins supporting cold adaptation of primary protein structures (Figure5, Table2). In the following section all observed changes of amino acid usage are discussed in respect of the Antarctic species when compared to the Z. viviparus transcriptome.

The most obvious shifts found for P. brachycephalum are the net loss of glutamic acid and asparagine and the gain of serine (compare values of Table2 and Figure5). The frequency of all acidic residues is reduced and contrasts with slightly higher frequencies for the basic amino acids lysine and arginine as well as histidine. Due to the eminent loss of glutamic acid, the total amount of charged residues for stabilizing salt bridges is reduced. This may contribute to weaker interaction involved in substrate binding and protein interaction. The net loss of glutamic acid and a preferential exchange with 85 aspartic acid residues preserves the acidic function at reduced flexibility of the side chain. These exchanges may cause increments in charge density at the surface favouring increased solvent interaction as postulated earlier[11, 42]. A reduction of salt links results from the exchange of the charged glutamic acid by polar glutamine (30 Q over E), in line with earlier assumptions drawn from analyses of individual proteins[11, 42]. Arginine as a bulky amino acid was identified to be less abundant in cold-adapted species[11, 42]. However, replacements of 30 serine into arginine and a slight increase of arginine on a global level contrasts with both former assumptions.

Based on genomic comparisons Wang and Lercher[30] developed a simple, reduced predictor based on replacement frequencies of the charged amino acids glutamine (E), arginine (R) and lysine (K), which were found most often exchanged in a thermal cline. Higher values for the ERK-proxy are found to be characteristic for hyperthermophiles compared to thermophiles and mesophiles. However, the ERK signal in our study is purely based on the large E change (Figure5), as a minor change of R+K in the reverse direction does not confirm the ERK-hypothesis perfectly. We interpret this partial contrasting outcome as a consequence of extrapolation, i.e., application of a thermophile-mesophile hypothesis below the mesophile scope seems to be appropriate only for the E signal. Furthermore, we are aware of the finding of Lobry and Necsulea[40] that thermophile-mesophile signatures of cold adaptation are more prominent[43] than mesophile-psychrophile signatures, at least in analyses of codon-usage trends in coding sequences of complete genomes.

Within the group of polar amino acids the picture is heterogeneous due to the large gain of serine and the net loss of asparagine. Comparing the total counts of gains and losses of polar with unpolar residues a gain of 102 polar residues remains. This indicates that despite this pooling unpolar to polar shift -supporting increased solvent interaction- cannot be resolved as a strong signal within our study.

By comparing a limited set of prokaryotic proteome sequences amino acid frequencies canonically discriminating psychropilic, mesophilic and thermophilic species were determined after modelling of candidate genes on existing three-dimensional structures[43]. In this way the analysis of particular orthologuous sequences were extended to part of the proteomes, and amino acid frequencies were assessed separately in structural categories ‘buried’ and ‘surface’. In general, a trend towards polar residues (in particular serine) resemble the finding of a (solely significant) preferential usage of serine and a net gain of polar residues for a cold-adapted species. Our study confirms these findings as a cold induced positive serine shift is one of the dominating signals in our study. Therefore an increased protein surface-solvent interaction in cold-adapted proteins can be hypothesized, at least on the basis of the net serine signal. However, for a confirmation, further sequence analyses are required to gain structural discrimination of solvent exposure. A diverse pattern is detectable for threonine, accounting to a small net gain similar to serine: the surplus of 82 exchanges over alanine significantly increases the polarity, possibly contributing to less hydrophobic interactions within and improved solvent interactions at the surface of proteins. Similarly, the net gain of threonine over asparagine contributes to the overall reduction of the latter and may support better solvent interaction through increased polarity. A surplus of 30 serine exchanges for threonine was detected. At preserved polarity this exchange is likely to increase the accessibility of the hydroxyl group. In contrast, methionine is preferred over threonine (46 M over T), decreasing the polarity. In summary, the multitude and diversity of exchanges within the polar category of amino acids together with the strong signal found in an earlier study[43] point to their importance for thermal adaptation of proteins.

Within the group of unpolar residues no prominent net gains or losses were detected (Figure5). It should be noted that a moderate gain in the usage of proline for P. brachycephalum is detectable, contrasting conclusions of existing views[11, 42] as this amino acid causes a large negative impact on structural flexibility. Several amino acid exchanges with conserved functions became apparent in the hydrophobic domain: net gains of 31 alanine in exchange for leucine reflect the preference of shorter residues, i.e., reducing entropic and enthalpic net contributions of side chains to structural stability of protein cores[8, 10, 11]. The finding of 51 valine over alanine and 65 isoleucine over valine contrast with the global pattern of reduced hydrophobic interaction. However, it is likely that an enlargement in unpolar residue length may enhance the protein–lipid interaction, e.g., within membranes, and must not be in contradiction with the flexibility hypothesis.

As only minor changes in sequence seem to be necessary to adapt kinetic properties, adaptive changes may be difficult to be identified unless proteins of closely related species from different thermal habitats are compared[5, 8]. Furthermore, studies based on structural attribution are restricted to proteins with available 3D structures, which is possibly biased towards soluble proteins with enzymatic function as these proteins are overrepresented in common databases. The approach of our study considers all types of expressed proteins within the transcriptome including membrane proteins, structural proteins, etc. without further assumptions.

We claim that our study is less distorted by a probable bias from a specific selection of species, in that we restrict our analyses on model orthologous sequences from two comparable transcriptomes. Furthermore, we focussed on aligned segments, excluding insertions, i.e., we hypothesize that adaptation can be analysed on a purely local basis at single amino acid positions.

In summary, the observed amino acid substitutions represent a global net pattern of molecular shifts in the local, i.e., position-specific amino acid composition for two species of the same family, inhabiting thermally very distinct sites.

Codon usage

A holistic approach for studying temperature-dependent evolutionary profiles on genetic up to proteomic levels was initially subjected to archaea, bacteria and only some eukaryotes on a large thermal scale[14, 44], uncovering trends of a higher GC usage for species with a higher OGT. As the GC content reflects the degree of hydrogen bonding in nucleic acid chain molecules, higher GC contents lead to an increased thermal stability in (deoxy-) ribonuclein acid chains[45, 46]. Thermal adaptation was analyzed in structural ribosomal RNAs for prokaryotic 16SrRNA finding rising GC contents in species with higher OGTs[47]. In vertebrates similar observations were detected in the 18SrRNA[48] with a trend for higher GC content in endotherm animals compared to poikilotherms showing a correlation to the environmental habitat temperature of the species under study. Similarly, the coding parts of genomes of cold blooded vertebrates and mammals are proposed to be separated by a “major compositional transition” in the GC content, resulting in nearly 100% GC3 levels in mammals[15]. In contrast, genomes of poikilotherm vertebrates are characterized throughout by lower GC3 levels.

Based on the identified segment pairs on the coding sequence level, we focused on a subset of transcript sequences of equal size in both libraries to uncover temperature-related patterns in the codon usage. This filter is supposed to allow for a more specific statistical comparison in that it is not distorted by effects of insertions. We furthermore filtered for highly similar and aligned sequences (80% minimum translated sequence identity) to reduce false pairing. Admittedly, total codon frequencies are on comparable levels for the two species (Figure7). However, by aggregating all pairwise comparisons in a single Fisher’s test we were able to detect a significant signal (p = 5*10^-4) of increased AT3 for the Antarctic eelpout (Figure7). A set of biochemically meaningful hypotheses may be proposed to explain the latter findings. Firstly, under energetically constrained conditions such as in polar environments, proof reading of DNA might be too cost-inefficient and an increment of A/T in the last codon position encoding for synonymous amino acids might be advantageous without having any impact on protein function. A subsequently lowered GC content is therefore discussed for the cold adapted diatom Fragilariopsis cylindrus (Thomas Mock, personal communication). Secondly, translation itself can be repressed kinetically through increased GC in codons, over-stabilizing enzymatic transition states within the chain of reaction steps of translation in the cold. Thirdly, to avoid over-stabilized secondary structure elements in messenger RNA in the cold lowering GC content could be an evolutionary resort. For example, the molecular function of cold shock proteins (CSPs) has been described to prevent mRNA from frozen, i.e., over-stabilized hairpin conformations in the cold[49]. An adaptation to the cold by means of reducing GC in the mRNA would imply a reduced pressure to express CSPs permanently at high levels, with implications for metabolic cost.

Phylogenetic analyses of fish

Animal collection and sample preparation of Pachycara brachycephalum

Specimens of P. brachycephalum were caught with baited traps at the position 62°10.9' S 58°20.8' W during expedition ANTXV/3 with the RV “Polarstern” in 1998. The animals were brought to the Alfred Wegener Institute in Bremerhaven and kept at 0°C in re-circulated seawater at 34 PSU (practical salinity units). The fish were fed ad libitum with Crangon crangon once a week; feeding was terminated exactly one week before sampling.

For sampling, the fish were anaesthetized by exposure to MS222 (0.2 mg*l^-1) before being killed. The fish (n=9) had a mean body length of 26 ± 0.83 cm (± SEM) and a mean body weight of 69.77 ± 5.78 g. Tissue samples were quickly excised, frozen instantaneously in liquid nitrogen and stored at −80°C until further processing.

Handling and killing of the fish were conducted in line with the recommendations of the American Veterinary Medical Association (AVMA). The work was approved by a competent German authority (Freie Hansestadt Bremen, reference number 522-27-11/02-00(93).

cDNA synthesis, normalization and pyrosequencing of Pachycara brachycephalum

Total RNA was extracted from 20–40 mg tissue with the Qiagen RNeasy kit with a modified protocol (proteinase K digestion after homogenisation) according to the manufacturer′s instruction (Qiagen, Hilden, Germany). Quantity and purity of the RNA were determined using the NanoDrop ND 1000 (Peqlab Biotechnologie, Erlangen, Germany). Integrity of the RNA was analysed by capillary electrophoresis (bioanalyser: Agilent, Waldbronn, Germany). Liver- and heart RNA samples of fish were pooled in equal amounts and used for the synthesis of cDNA.

Preparation of a random-primed and normalized cDNA for pyrosequencing was done with 100 μg of pooled total RNA from liver and heart (Vertis, Freising, Germany). Poly(A)+ RNA was prepared and first strand cDNA was amplified with a N6 randomized Primer. 454 sequencing adapter A (5^′-GCCTCCCTCGCGCCATCAG-3^′) and B (5^′- CTGAGCGGGCTGGCAAGGC-3^′) were ligated to the 5' and 3' ends of the cDNA, followed by an amplification for 21 cycles with a proof-reading enzyme. Normalization was carried out by one cycle of denaturation and reassociation of the cDNA, resulting in N1-cDNA. Reassociated ds-cDNA was subtracted from the remaining ss-cDNA by separation with a hydroxylapatite column. The purified (and normalized) ss-cDNA was amplified with 8 PCR cycles.

The resulting cDNA product was loaded to a preparative agarose gel (1.5%) and fragments in a size range of 450 – 650 bp were eluted for sequencing. The application was performed at Eurofins-MWG, Germany according to[50] using the GS FLX Titanium-technique (MWG-Biotech AG, Ebersberg, Germany). High-density pico reactions were performed on a half of a sequencing run.

Quality assessment and assembly of cDNA of P. brachycephalum

The sequencing of the normalized cDNA library of Pachycara brachycephalum resulted in 481,802 reads with an average read length of 321 bases (Table1). The reads were preprocessed first with a base calling quality control, followed by a polyA-clipping and a screening for remaining sequencing primer adapters in the 5′as well as in the 3′ends that were truncated. The assembly was done with standard settings for stringency and homology by the Mira Assembler Version 2.9.43 (MWG-Biotech AG, Ebersberg, Germany) according to the instructions by MIRA[51, 52]. The high quality reads (338,993) were assembled into 65,565 contigs while 123,038 singlets could not be matched against any other reads. The contigs have a mean length of 487 bases and 22,651 sequences were larger than 500 bases (cf. Figure2A). The larger contig length compared to Zoarces viviparus may be due to different sequencing protocol as we also used the FLX sequencer with a protocol including Titanium chemistry.

Functional annotation and comparisons of zoarcid cDNA libraries

The assembled contigs of Pachycara brachycephalum were loaded into the free software tool BLAST2GO[26] and BLASTed against the NCBI non-redundant database (BLASTx)[53] with an e-value cut-off of 1.0^-3 and a HSP cut-off length of 33 bases. In total 47,584 Sequences had a BLAST hit, while the e-value distribution (Figure3A) and the HSP/hit distribution (Figure3B) reflects the good quality of the P. brachycephalum library. We processed the library of contigs of Zoarces viviparus in the same way[23] and found in the 53,459 sequences 35,133 with a BLAST result. After the GO-mapping step in BLAST2GO we processed both libraries for annotation with no cut-off for the HSP/hit -coverage, resulting in 19,460 sequences for P. brachycephalum and 16,315 sequences for Z. viviparus.

The latter sets of sequences were used for evaluating amino acid usage differences and for synonymous/non-synonymous codon usage comparisons. For this, we generated partial sequence translations by rpstBLASTn against a collection of fish-specific orthologies (protein sequences from fiNOG from eggNOG version 3)[54] with e-value below 10^-9, which yields 52,729 (P. brachycephalum and 34,733 (Z. viviparus) hits in a total set of 7,374 fiNOG-orthologies. Pairwise alignments between sequences of the two libraries in translated form were then obtained from BLASTP of the obtained fiNOG-compatible translations. 5,352 pairwise alignments (best filtered HSP per fiNOG) with percent-identity above 80% and e-value below 10^-9 were obtained. Coding sequences for the aligned translated segments were aligned accordingly to the protein alignments. Removal of sequence pairs containing stop-codons and/or ambiguity nucleotide-codes yielded a set of 4,155 segments for further analyses. GC3 comparisons are based on aligned synonymous codon positions with unchanged nucleotides in the first and in the second codon positions.

Amino acid usage has been analysed from these alignments as well as synonymous codon usage (Within Canonical Analysis), both in R[55] with the packages seqinr and ade4 (after[31]). The re-aligned sequence segments contain pairwise information that can be compared in a position-specific way, allowing to quantify preferred amino acid usage, e.g., by counting amino acid specific replacement frequencies for amino acid i in P. brachycephalum to amino acid j in Z. viviparus at the same position in the alignment. A non-synonymous replacement imbalance is computed from the matrix.

For a functional characterisation of the libraries, rpstBLASTn was used against the database of metazoan orthologies computed from eggNOG-alignments (version 2,[29]). For comparison with a reference genome of a related fish, we analysed the coding sequences from the genome of G. aculeatus (BROADS1, ENSEMBLE Version 1.63 at EMBL) in the same way. Best hits to meNOG orthologies with e-value below 10^-20 were kept, pooling functional annotations into COG/KOG categories.

Access to data

The transcriptome of P. brachycephalum is published in the sequencing-read archive (SRA) at NCBI under Accession SRA049761.