Annotation of expressed sequence tags for the East African cichlid fish Astatotilapia burtoni and evolutionary analyses of cichlid ORFs

The exceptionally diverse species flocks of cichlid fishes in the East African Great Lakes Tanganyika, Malawi and Victoria are prime examples for adaptive radiations and explosive speciation [1–3]. More than 2,000 cichlid species have evolved in the last few million years in the rivers and lakes of East Africa [1, 4–6]. Together with an additional ~1,000 species that are found in other parts of Africa, in South and Central America, in Madagascar, and in India, the family Cichlidae represents one of the most species-rich families of vertebrates. In addition to their unparalleled species-richness, cichlids are famous for their ecological, morphological and behavioral diversity [1, 2, 7], for their propensity for rapid speciation [5], for their capacity for sympatric speciation [8, 9], and for the formation of parallel characters in independently evolved species flocks [10–12]. For these reasons, the cichlid fishes are an excellent model system to study basic dynamics of evolution, adaptation and speciation. However, while the phylogenetic relationships between the main cichlid lineages are largely established and some of the cichlids' evolutionary innovations have been identified [1, 2, 4, 7, 13], little is known about the genomic and transcriptional basis of the evolutionary success of the cichlids.

The cichlid model system provides many advantages for evolutionary genomic research. The hundreds of closely related yet morphologically diverse species in East Africa's cichlid species flocks are even more powerful than a 'mutagenic screen' (to which these species assemblages have been compared [1, 12]) in that they represent combinations of alleles that confer a selective advantage under various ecological pressures. Because of the possibility to produce viable crosses between different cichlid species in the lab [14], these alleles can be tied to particular phenotypic traits by means of classical genetic experiments [15–18]. The close relatedness of the different species allows the design of primer sets for the amplification of particular genomic DNA regions such as candidate gene loci, microsatellites, or SNPs, which are applicable to a wide range of species [17, 19–21]. The same is true for expression profiling with cDNA microarrays that, once developed for one species, can be used for any East African cichlid species [22].

A variety of genomic resources have already been established for East African cichlid species. Genetic maps are available for the Nile tilapia Oreochromis niloticus [23, 24] and the Lake Malawi species Metriaclima zebra [17]. BAC libraries have been constructed for O. niloticus [25] and M. zebra (available at the Hubbard Center for Genome Studies), for the Lake Victoria haplochromine Paralabidochromis chilotes [26] and for Astatotilapia burtoni from Lake Tanganyika and surrounding rivers [27]. cDNA microarrays are available for A. burtoni [22] and for Lake Victoria haplochromines [28, 29]. Also, EST sequencing projects have been initiated [30], and a BLAST server for cichlid resources has been established [31]. Recently, the National Institute of Health (NIH) has committed to sequencing four cichlid genomes. A detailed description of genomic resources developed for cichlid fishes is available at [32].

Expressed sequence tags (ESTs) derived from the partial sequencing of cDNA clones provide an economical approach to identify large numbers of genes that can be used for comparative genomic and gene expression studies as well as for the detection of splice variants [33, 34]. Furthermore, EST projects facilitate genome annotation and are therefore often applied in addition to genome sequencing projects. Due to the large amount of data available in public databases, ESTs emerge as important resources for comparative genome-wide surveys both among closely and more distantly related taxa [35, 36]. A series of software applications have been developed to date to perform such EST-based analyses [37–39]. Since ESTs reflect the coding portions of a genome, they can also be used to test for different evolutionary rates in particular genes when comparing different lineages, and to detect genes that have undergone positive selection [35]. It is generally assumed that genes with a statistically significant increase in substitution rates have experienced relaxed functional constraints, while genes, which have not undergone accelerated substitution rates, have experienced purifying selection and, thus, could not accumulate substitutions at random. Positive Darwinian selection, on the other hand, is a phenomenon where selective pressure is favoring change. Natural selection is commonly thought of as a process of editing genetic change so that only a small number of mutational events are retained in natural populations. Under positive selection, the retention of mutations is much closer to the rate at which mutations occur.

Here we report the collection and annotation of more than 12,000 ESTs generated from two different cDNA libraries obtained from the East African cichlid species Astatotilapia burtoni, as well as a smaller cDNA library from the Lake Malawi species Metriaclima zebra. Astatotilapia burtoni has long been used as a model system to study cichlid spawning behavior [7, 40, 41], social interactions [41–44], neural and behavioral plasticity [45, 46], endocrinology [47], the visual system [48], as well as cichlid development and embryogenesis [49]. In addition, the phylogenetic position of A. burtoni makes this species an ideal model system for comparative genomic research [27]. Astatotilapia burtoni, which belongs to the most species-rich lineage of cichlids, the haplochromines, was shown to be a sister group to both the Lake Victoria region superflock (~600 species) and the species flock of Lake Malawi (~1,000 species) [4, 5, 50, 51]. Three highly specialized haplochromine species from two species assemblages, Paralabidochromis chilotes and Ptyochromis sp. "redtail sheller" from Lake Victoria and Metriaclima zebra from Lake Malawi, have already been established as genomic models [16, 26, 28, 30]. Important insight into cichlid (genome) evolution will be afforded by the comparison of their genomes to that of A. burtoni, which has a more generalist life style and is likely to resemble the ancestral lineage that seeded the cichlid adaptive radiations in these two lakes [4, 7].

For EST sequencing, we utilized a cDNA library from A. burtoni brain tissue ('brain') that was used for the construction of a cDNA microarray [22] and a newly generated normalized cDNA library constructed from different A. burtoni tissues at different developmental stages ('pinky'). We annotated the ESTs on the basis of similarity searches with BLAST and using the structured vocabulary provided by the Gene Ontology Consortium [52], based on molecular studies of gene function in various model organisms [53]. For evolutionary analyses, we combined our newly generated ESTs with all available sequence data for haplochromine cichlids [30] and a previously constructed library from skin tissue of the Lake Malawi species Metriaclima zebra (W. Salzburger, H. A. Hofmann & A. Meyer; unpublished data), which resulted in a total of more than 45,000 ESTs. We then compared the haplochromine ESTs to sequence data from two pufferfish species (Takifugu rubripes and Tetraodon nigroviridis), trout, and zebrafish, and identified those ESTs with cichlid specific differences in evolutionary rates with EverEST [37].

A total of 8,636 A. burtoni sequences assembled into EST contigs have an open reading frame (ORF) of at least 400 bp. Of these, 1,219 (14%) had matches in the Takifugu database and 7,417 (86%) had no matches when an expected value threshold (e-value) of < 1 × 10^-50 was used. 2,902 (34%) had matches in the Takifugu database with an expected value threshold of < 1 × 10^-15 and 3,460 (40%) had matches with an expected value of < 1 × 10^-5. Similar proportions were retrieved with other databases (Fig. 1).

Among the 8,363 A. burtoni assembled sequences, 2,977 could be annotated according to Gene Ontology (GO) terms. Additional files 1, 2 and 3 use the generic GO slim subset of terms ([54]; Generic GO slim; Mundodi and Ireland; downloaded 04/06/2007) that have been developed to provide a useful summary of GO annotation for comparison of genomes, microarrays, or cDNA collections when a broad overview of the ontology content is required. 2,692 ESTs could be assigned to genes listed in the molecular function ontology, 2,532 to genes listed in the biological process ontology, and 2,293 to genes listed in the cellular components ontology, when using an e-value of < 1 × 10^-12. Additional files 4, 5, and 6 provide more detail of the specific fine-grained terms. Because a single A. burtoni assembled sequence may be annotated in all three ontologies and according to multiple ontology terms, a total of 27,451 annotations have been applied (10,926 among biological process, 9,414 among molecular function, and 7,111 among cellular component).

For the comparative evolutionary analyses, we combined our newly generated ESTs with previously published data from Paralabidochromis chilotes and P. sp. "redtail sheller" [30] and about 1,000 sequences obtained from a Metriaclima zebra skin cDNA library (W. Salzburger, H. A. Hofmann & A. Meyer; unpublished data). When using this set of haplochromine cichlid ESTs as reference, we identified 759 open reading frames that are present in all six databases used for comparative analyses (haplochromine cichlids, Danio rerio, Homo sapiens, Oncorhynchus mykiss, Takifugu rubripes, and Tetraodon nigroviridis).

In order to identify sequences that evolve significantly more rapidly or more slowly in the haplochromine cichlid, we applied the triangle method implemented in EverEST [37] to calculate the p-distance for each of these 759 ORFs in all fish species relative to the human ortholog. There were 22 cases in which more than one haplochromine sequence was found. In these cases, we used the longest sequence for further analyses. The relative p-distances for three fish species were then mapped in ternary diagrams. An example of such a ternary diagram is shown in Fig. 2a, in this case showing the relative p-distances of cichlid, Takifugu rubripe s, and Danio rerio amino acid sequences with respect to the homologous Homo sapiens genes. Figure 2b depicts a diagram with Oncorhynchus mykiss amino acid sequence divergence instead of haplochromine cichlid. The ternary diagrams show that in all combinations most genes are clustered around the center of the respective triangle, which indicates that, in general, the p-distances relative to the human outgroup are similar in all fish species.

When compared to the green-spotted pufferfish (Tetraodon nigroviridis) and fugu (Takifugu rubripes) (always with human as outgroup), 49 gene fragments appeared to have a significantly faster rate of evolution in haplochromine cichlids, and 213 had a slower rate. In the comparison including zebrafish and fugu, 52 genes were found to have evolved faster and 185 genes slower in cichlids. When trout and zebrafish were used, 69 genes were faster and 139 genes evolved slower. In a comparison including trout and fugu, 68 genes appeared to have a faster rate in haplochromines, and 132 had a slower rate. In total 69 genes were found to have evolved faster, and 213 genes appeared to have evolved with a significantly slower mutation rate in haplochromines compared to other fish species. Altogether, about 22% of the surveyed ESTs were found to have haplochromine specific rate differences in at least one of the comparisons suggesting that these genes might play a role in lineage specific features of haplochromine cichlids. A set of 170 cichlid genes appeared in all comparisons. Forty-eight cichlid genes were found to have a higher rate of amino-acid substitution compared to the other fish species included in this study, while 122 cichlid genes were found to have a slower rate. Cichlid sequences that match Danio rerio, Takifugu rubripes, Tetraodon nigroviridis, and Oncorhynchus mykiss genes and have a significantly higher or lower p-distance compared to the other fish genes relative to the human outgroup are listed in Additional files 7 and 8, respectively.

A histogram of the abundance of amino acid sequence divergences of all five fish species with respect to homologous human genes is depicted in Fig. 3. The p-distances appear normally distributed. With 0.211, cichlids show the lowest average distance followed by Oncorhynchus mykiss (0.216), Danio rerio (0.239), Takifugu rubripes (0.242), and Tetraodon nigroviridis (0.258). The average distance of all five fish species to Homo sapiens is 0.233. We also used the 482 redundant sequences that were found in all three large haplochromine cichlid EST datasets (P. chilotes and P. sp. "redtail sheller" [30]; Astatotilapia burtoni, this study) to calculate mean pairwise p-distances. Within these three cichlid species, we found a mean p-distance of 0.14 between A. burtoni and P. chilotes, 0.17 between A. burtoni and P. sp. "redtail sheller", and 0.08 between the two Lake Victoria species P. chilotes and P. sp. "redtail sheller".

We then calculated K_a/K_s ratios for all genes with a higher or slower rate of base substitution in cichlids. K_a/K_s ratios greater than one, which are indicative of positive selection in that gene, were found in four genes that evolve more slowly in cichlids compared to the other fish species. The highest K_a/K_s ratio (3.77) was found in the neuroendocrine convertase subtilisin/kexin type 1 that is responsible for processing large precursor proteins into mature peptide hormones [55, 56]. In claudin 3, a member of the claudin family involved in the formation of tight junctions in various tissues [57], the K_a/K_s ratio was 1.55. A K_a/K_s ratio of 1.30 was observed in the catalyzing enzyme glutathione peroxidase 3, and a ratio of 1.19 was found in ménage a trois 1 (MNAT1), which is a member of the CDK7-cyclin H complex that functions in cell cycle progression [58], basal transcription, and DNA repair.

Expressed sequence tags are important genomic resources and their numbers in public databases such as GenBank are rapidly increasing. Full-length cDNA and EST sequencing projects typically accompany genome sequencing projects, as these data are essential for the recognition and annotation of genes, the characterization of the transcriptome, the identification of intron-exon boundaries and the detection of splice variants in eukaryotes, etc.[33, 34, 59–61]. In addition, the standardized procedure of cDNA library construction and normalization, and the comparably low costs of large-scale DNA sequencing facilitate EST projects in organisms for which the whole genome sequencing has not (yet) been completed. Thus, EST sequencing projects outnumber genome-sequencing projects – particularly in groups with larger genome sizes such as plants and vertebrates – leading to a large body of sequence data available for comparative analyses. Large-scale EST analyses have been used in many other contexts, such as primary gene expression assays [62, 63], the estimation of the total number of genes in an organism [64], cDNA microarray annotations [65], or the construction of genetic linkage maps [66–68]. Expressed sequence tags can further be used for phylogenomics [36, 69], and for the identification of microRNAs [70].

Despite their many advantages, there are also some problems associated with ESTs. For example, EST sequences typically cover only parts of a gene, so that two sequences of the same gene might not necessarily overlap. That only fragments of a gene are available also leads to problems with homology-based analyses such as BLAST. Then, EST sequences often contain the untranslated regions (UTRs) that are present in mRNAs but do not translate into amino acids. Finally, it is often difficult to figure out the proper reading frame, particularly in shorter ESTs, which impedes certain analyses. A combination of multiple EST projects (as we have done here) helps to alleviate some of the shortcomings inherent in EST data.

We have sequenced, annotated and conducted evolutionary analysis of ESTs of haplochromine cichlids for several reasons. First, this large set of sequence data for cichlid ORFs provides insight into the genome of a representative of haplochromine cichlids, which are a main model system for the study of adaptive evolution and explosive speciation [1–3]. Second, we wanted to extend the existing genomic resources for Astatotilapia burtoni such as a genomic BAC library [27] by establishing cDNA libraries from different tissues. Furthermore, these cDNA libraries provide the basis for annotated cDNA microarrays that are being used for expression analyses in a variety of cichlid species [22, 28, 71]. Finally, we were interested in identifying genes with a different evolutionary rate in the rapidly radiating cichlid lineage compared to other fish species, as well as in identifying genes that show the signature of adaptive evolution in cichlids.

Of the two A. burtoni cDNA libraries that were used for EST sequencing, the normalized mixed tissue library ('pinky') was of better quality. Not only were there much fewer redundant sequences as compared to the brain library, which was mainly due to the normalization step, but also the average insert size was larger and the average read length was longer. Altogether, about 85% of the sequenced cDNA clones led to high-quality ESTs of a length of >200 bp (86% in pinky, and 85% in brain). In the BLAST searches against Takifugu rubripes, Tetraodon nigroviridis, and Danio rerio, between 14% (when compared to T. rubripes; e-value ≤ 10^-50) and 43% (when compared to D. rerio; e-value ≤ 10^-5) of the A. burtoni ESTs led to hits (Fig. 1). This lies well within the range of other EST sequencing projects [63, 65, 72].

About 8,600 A. burtoni ORFs (or 75% of the high quality ESTs) were longer than 400 bp, and about 3,000 sequences could unambiguously be annotated and classified following the vocabulary provided by the Gene Ontology Consortium [Additional files 1, 2, 3, 4, 5, 6]. According to the Gene Ontology classification, it appears that a broad range of genes involved in functions, processes and compartments are represented in our EST set. This cichlid specific GO slim offers several advantages. First, it offers a rapid visual interpretation of gene subsets. Second, because the cichlid specific slim is built from those sequences used to build a cDNA microarray, it offers maximal power when testing for over- or under-representation of gene lists while reducing the need for correction for multiple hypothesis testing. Finally, it allows for a less experimenter-biased interpretation of microarray results, or other genomics analyses in a manner that can be easily compared between experiments.

One of our main goals was to characterize genes in haplochromine cichlids that show a faster or slower rate of base substitutions in cichlids compared to other fish species, as this is indicative of a relaxed or reinforced selection regime, respectively [35]. To this end, we combined our newly generated ESTs with previously published sequences for Lake Victoria haplochromine cichlids [30] and about 1,000 sequences obtained from a Metriaclima zebra skin cDNA library, which resulted in a total of about 45,000 ORFs. By means of homology searches against human, the two pufferfishes, trout, and zebrafish using local BLAST, we identified a set of 759 ORFs that are present in all species and that show a sufficient degree of homology (e-value ≤ 10^-50) for further analyses with EverEST [37]. The number of genes with a cichlid-specific faster or slower rate of molecular evolution (always with human as outgroup) varied when different fish taxa were used in addition to the cichlid ORFs. However, we found a set of 170 genes (48 "faster" and 122 "slower"; Additional files 7, 8) that appeared in all comparisons and are, thus, good candidates for playing an important role in the evolution of (haplochromine) cichlid fishes.

When characterizing these genes further, by means of calculating K_a/K_s ratios, we found that four genes (or 2.35% of all deviating genes) showed the signature of adaptive evolution in the haplochromine lineage. The highest K_a/K_s ratio (3.77) was found in the neuroendocrine convertase subtilisin/kexin type 1, followed by claudin 3, (1.55), glutathione peroxidase 3 (1.50), and ménage a trois 1 (1.19). All gene fragments that show a K_a/K_s > 1 are found among the more slowly evolving genes. These genes are now candidate genes for further investigations. The gene with the highest K_a/K_s ratio appears particularly interesting. It is known that neuroendocrine factors, such as gonadotropin releasing hormone (GnRH), are involved in regulation of reproduction and behavior in A. burtoni [56, 73].

There was a relatively even distribution of rapidly evolving genes across all GO categories. Only three terms, "response to stimulus", "cellular process" and "actin cytoskeleton" deviated significantly from the distribution expected by chance alone. The most significant disproportionate under-representation for the rapidly evolving genes was the category of response to stimulus for which only 1 of the 86 possible annotated ORFs was included on the list.

The distribution across GO categories was highly non-uniform for the slowly evolving genes. Many categories from each ontology were represented by significantly more or fewer ORFs than would be expected by chance. Among those terms over-represented we found several relating to cellular processes such as protein kinase cascade, mitosis, and cell signaling as well as growth and cell adhesion, while metabolic process was under-represented along with the secretory pathway category.

The GO analysis highlights the possible categories of genes that may play an important role in the evolution of the haplochromine cichlid fishes. This analysis presents hypotheses to be tested through focused experimental or sequence analysis. An interesting contrast in GO analysis results was observed between the rapidly evolving genes that showed little tendency to derive from a particular class and slowly evolving genes that were more structured in their distribution. The lack of structure to the distribution of rapidly evolving genes may reflect the possibility that specialization among cichlids occurs along diverse biological pathways rather than a repeated divergence of a given biological process or molecular function. The GO categories that are over-represented among slowly evolving genes could represent genes whose functions are important for phenotypic plasticity or other traits linked to the successful adaptive radiation, while those categories that are under-represented by slowly evolving genes represent categories that are not as tightly constrained.

Our p-distance comparisons between the five fish species and human (as outgroup) also revealed that cichlids show the lowest average p-distance compared to Homo sapiens (Fig. 3). This might be an artifact that is due to the use of the haplochromine cichlid sequence as query for all BLAST searches. Alternatively, as we also found 122 slowly evolving genes in haplochromine cichlids, there might be a tendency in haplochromines to retain ancestral forms and functions. The pairwise average p-distance comparisons between the three cichlid species Paralabidochromis chilotes, Ptyochromis sp. "redtail sheller", and Astatotilapia burtoni revealed that the coalescence time between the two Lake Victoria species (0.08) is about half compared to their coalescence time with A. burtoni (0.14 and 0.17, respectively), which is in concordance to the phylogenetic relationships between these three taxa [4].

Here we report the sequencing and annotation of more than 11,000 ESTs from the East African haplochromine cichlid Astatotilapia burtoni. Our EST set comprises a broad range of genes involved in functions, processes and compartments. By combining the A. burtoni ESTs with publicly available ORFs from two Lake Victoria haplochromines and subsequent comparisons to other fish model systems, we identify a set of 170 genes with haplochromine-specific differences in evolutionary rates. These genes appear as good candidates for playing an important role in the evolution of the exceptional diversity found in (haplochromine) cichlids. Interestingly, genes that were more slowly evolving in the cichlid lineage were not evenly distributed across Gene Ontology categories; classes that are over-represented could represent genes whose functions are important for successful adaptive radiation. We also identify four genes with a K_a/K_s ratio greater than one, which are, hence, likely to have undergone positive selection in haplochromines. The A. burtoni ESTs provide novel insights into the genome of haplochromine cichlids and will serve as valuable resource for researchers working in the field of (cichlid) evolutionary genomics, particularly in the light of the forthcoming sequencing of four cichlid genomes.