Comprehensive EST analysis of the symbiotic sea anemone, Anemonia viridis

ESTs generation and analysis

A cDNA library was made from both symbiotic (gastroderm) and non-symbiotic (epiderm) tissues of the sea anemone Anemonia viridis (Additional file 1). In order to maximize the diversity of genes expressed under both normal and stress conditions within the symbiotic association, sea anemones were subjected to different environmental conditions before RNA extraction and cDNA library construction (light/dark sampling, thermal stress, hyperoxia conditions, see Materials & Methods). Out of 50,304 sequenced clones, a total of 41,247 readable sequences were produced, corresponding to a sequencing success of 82%, and 39,939 high quality ESTs were generated (Figure 1 and Additional file 3). A. viridis repetitive sequences were identified and masked using RepeatMasker before assembly. Trimmed masked ESTs were further subjected to cluster analysis using the TIGR-TGICL pipeline with default parameters [32]. A total of 30,087 ESTs were assembled into 4,652 contigs (Additional file 4) while 9,852 sequences remained as singletons. The mean trimmed EST sequence length was 626 bp, but sequences could be assembled into contigs of up to 4 kb (CL4Contig4, described as ribosomal protein L8, 3,887 bp). Most contigs (2,959) were composed of three ESTs or less, with an average length of 895 bp, suggesting that they covered only parts of the transcript sequences. Apart from the functional genomics resources available for the reef-building corals Acropora palmata and Montastrea faveolata [24], most cnidarian genomic studies have been performed on non-symbiotic tissues (aposymbiotic developmental stages or adults). To our knowledge, our study provides the largest EST collection from a symbiotic cnidarian.

Repeated elements

Based on reverse transcriptase domain search, 46 transposable elements (retrotransposons/retroposons) were identified (data not shown). Among these, 4 were of prokaryotic origin (without any similarity to N. vectensis sequences) and 42 were of metazoan origin. While 2 of them were almost identical to N. vectensis transposable elements, 21 only had slight similarity to N. vectensis sequences (BlastX or BlastN with E-value of 1.10^-25 to 1.10^-10), and 19 showed no similarity to N. vectensis sequences but had been previously identified in other Metazoa.

Self BlastN analysis on our EST dataset also identified another abundant set of repeated elements. Over 1,500 repetitive sequences were identified, with 412 repeated more than 10 times in our UniSeqs dataset. A typical sequence would be around 200 bp length, harbouring 40–50 bp inverted terminal repeats at its extremities. These repeats often formed very long and stable hairpins, depending on their orientation (Figure 2). They were found in 19.8% of UniSeqs and were overrepresented within contigs (62% of contigs contained at least one repeat). When open reading frame could be identified, these A. viridis repeated elements were always mapped in the 3'UTRs of the protein coding genes. Using a BlastN analysis, none of these sequences were found either in the genome of N. vectensis or in the large EST datasets of Hydra magnipapillata, and search in Repbase did not show any homologs. However, members of these A. viridis repeated elements could be found in the limited EST datasets available for the symbiotic coral Acropora millepora and the sea anemone Anthopleura elegantissima. Further structural and functional characterisations of these repeated elements are currently being made.

Comparative analysis to related taxa

Because the cDNA library was made from symbiotic tissue, we expected to find ESTs related both to the cnidarian host and to its dinoflagellate symbionts. The analysis pipeline used in this study is presented in Figure 1. Sequences were compared with SwissProt (2008.03) and Uniprot KB (TrEMBL+Swissprot) (2008.03) databases using BlastX with a cutoff E-value of < 1.10^-10 to retrieve functional annotations. Of the assembled dataset, 6,238 UniSeqs had putative similarities while 8,266 had no similarity to any sequences in the chosen databases. Sequences were also compared with NCBI-indexed prokaryotic nucleic sequences (2007.08 release), using BlastN with a stringent E-value of < 1.10^-15 to assess the proportion of prokaryotic sequences. Finally, a specific search of the UniSeqs for virus proteins returned 7 hits. Relative contribution is shown in Figure 3. First of all, a relatively high proportion of sequences (57%) remained that showed no significant similarities to previously described genes and were therefore considered as 'unknown'. This is somewhat comparable with results obtained from other cnidaria, Acropora palmata, Acropora millepora, Montastrea faveolata, and Nematostella vectensis [19, 24]. Most of the UniSeqs identified in Symbiodinium sp were also of unknown origin [28]. Metazoan hits were found for 32.3% of UniSeqs (75% of annotated sequences). Among these, most of the annotated UniSeqs (4,266 out of 4685) matched with Nematostella vectensis predicted proteins. However, we also identified sequences that were clearly from the host (first Blast hits of metazoan origin with cutoff E-value < 1.10^-50), but these had no significant similarity to predicted proteins of N. vectensis. Three of these (a glycoprotein, a ferroxidase and an amine-oxidase) were studied in more detail (Figure 4). PCR and sequencing were first performed on genomic DNA from both A. viridis epiderm and in vitro cultured Symbiodinium, which confirmed the animal origin of these sequences (Figure 4A). The first sequence studied is related to the ependymin glycoprotein family (more specifically to the Mammalian Ependymin-Related Proteins group or MERP1), which has been well described in vertebrates due to its involvement in the regeneration processes. Ependymins are secretory proteins that can bind calcium and that were found predominantly in the cerebrospinal fluid of teleost fish. A bound form has been described, associated with the extracellular matrix. Recent data demonstrated that these proteins are also present in non-vertebrate deuterostomes and protostomes [35], and that positive selection may have shaped their evolution. Figure 4B illustrates an amino acid alignment of A. viridis MERP sequences to homologous proteins found in publicly available databases (Bayesian analysis of the MERP sequence confirmed the phylogenetic relationship, not shown). The presence of ependymin proteins in basal organisms clearly indicates for the first time that this protein family is far older than previously thought (first described as chordate-specific, then deuterostome-specific, and finally also found in protostomes).

In addition, there are at least 3 distinct MERP homologues in A. viridis, all of animal origin. As this is never the case in other phyla (only one homologue was found in any other species), we could hypothesize that MERP are important for symbiosis in cnidarians. The second sequence is related to a copper-dependent ferroxidase family protein (hephaestin) involved in copper detoxification, which has been described both in eukaryotes (unicellular eukaryotes and vertebrates) and prokaryotes. Nevertheless, no significant similarity was found in D. melanogaster or in C. elegans genomes, nor in any available protostome database, suggesting that this has been lost during protostome evolution. The third sequence is related to a prokaryotic amine-oxidase (permease) family protein, which has also been identified in the genome of C. elegans and other Caenorhabditidae but not in any other animal sequence database. Such species-specific occurrence could be explained by the presence of a similar prokaryote in both nematoda and A. viridis associated flora (sequence contamination) or a parallel gene transfer event during their evolution. However, molecular phylogenetic analyses (neighbor joining, maximum parsimony and maximum likelihood) showed that cnidarian and rhabditidae nematode sequences cluster to outside other prokaryote representative gene sequences (Additional file 5). The latter rather suggest that amine-oxidase genes were maintained in these taxa but lost in other metazoa.

Quite a large proportion of UniSeqs (9%) were uniquely shared with prokaryotes, of which Proteobacteria was the most prominent bacterial group (Pseudomonas, Bordetella, and Burkholderia). The holobiont has already been described as a dynamic assemblage, made up of the animal host, zooxanthellae, endolithic algae and fungi, Bacteria, and Archea [7]. Such "prokaryotic sequences" could therefore be assigned to the sea anemone-associated flora. In addition, a small but significant number of A. viridis sequences, recognized as being of prokaryotic origin based on Blast analysis, had already been identified in the N. vectensis genome. Such genes are clearly similar in sequence to prokaryotic homologs, although they contain introns. In N. vectensis and A. millepora they were proposed as "ancient genes", conserved in cnidarians but lost in other animal genomes [19]. Although our results are in line with this interpretation, comparative genomic studies among cnidarians, such as A. viridis, could help to identify the most probable evolutionary scenario between maintenance of "non-metazoan" genes in cnidarians or lateral gene transfer events, followed by rapid intron acquisition.

Surprisingly, only a small fraction of our dataset could be assigned to unicellular eukaryote sequences (putative Symbiodinium sp, 3.6% of annotated sequences). Two well accepted explanations have been proposed: i) poor representation of dinoflagellate sequences in databases, leading to wrong assignment after Blast analysis; and ii) technical bias due to the Symbiodinium cell wall impairing complete RNA extraction with standard methods, thus leading to an under-representative number of cDNAs in our library.

GC content was calculated for all UniSeqs to better assign a source to the transcripts (species of origin). Figure 5 represents the distribution of GC content in the three major taxa (metazoans, unicellular eukaryotes and associated prokaryotic flora), as well as for sequences without BlastX hits (unknown). Two distinct distribution curves were observed, with the metazoan ESTs averaging 40% GC, while bacterial ESTs averaged 63%. Despite our small number of Symbiodinium sequences, their GC content (58%) agreed with that already published for dinoflagellates, which have a relatively high GC content compared with other eukaryotes: Triplett et al [36] calculated a GC content of 64% for Heterocapsa pygmaea and an EST analysis of Alexandrium tamarense estimated that coding region GC-content was 60.8%, whereas GC-content in the 3'-UTR was slightly less, at 57.6% [37]. GC-content analysis of our unknown sequences revealed that ESTs were clustered around two peaks corresponding to 35% and 58% GC content for the larger and smaller peaks respectively. These data strongly suggest that most of the unknown sequences are likely to be of animal origin. As the GC content of 3' UTRs is usually lower than that of 5' UTRs or CDSs, these results are also consistent with our hypothesis that these sequences are mainly 3' UTR (due to cloning strategy), which also explains their low level of annotation. On the other hand, the smaller peak (58% GC) could be attributed to dinoflagellate sequences, although no dinoflagellate genome is currently available to test this, and the most closely related available genomes are those of apicomplexans (Plasmodium falciparum and Toxoplasma gondii), which are also compact genomes of parasitic protozoans.

To gain some insight into genes potentially involved in symbiosis among cnidaria, we performed comparative Blast analyses on the subpart of A. viridis sequences that were of animal origin (12,448 sequences). They were subjected to tBlastX (cutoff E-value of < 1.10^-10) against cnidarian EST sequences available in NCBI-dbEST (2009.05). A parallel BlastX (cutoff E-value of < 1.10^-10) was performed against the N. vectensis predicted proteins dataset. Table 2 gives the number of positive Blast hits (presence vs absence) obtained from A. viridis sequences against ESTs from selected cnidarian species. This comparison highlighted two species with a high number of homologous sequences, the sea anemone M. senile and the symbiotic stony corals A. millepora and A. palmata. In the Figure 6 Venn diagram we compared A. viridis with two non symbiotic species (N. vectensis and M. senile), and the symbiotic Acropora species: 1,535 sequences were common to the four species, while only 129 sequences (of which 88 have no BlastX hit) were specific to the symbiotic species, which probably contain candidate genes required for the symbiosis.

Functional classification of ESTs

Finally, this set of UniSeqs was annotated using similarity searches in both nucleotide and protein databases, as well as domain searches. Gene ontology terms were then assigned automatically using customized scripts based on InterProScan search results. We also submitted N. vectensis predicted proteins to the same InterProScan analysis for a comparative approach. To homogenize annotation level we only kept the root domain, using the hierarchical domain organization available from EBI. InterPro protein domains were identified in 40% of proteins. Based on similarity searches (Figure 3), we assigned putative taxa (A. viridis, Symbiodinium or associated flora) to each unique sequence. Figure 7 represents the distribution of assigned gene ontology terms. Comparison of domain occurrence in A. viridis with that in N. vectensis demonstrated that our library is highly representative of the Actinaria transcriptome. We also identified increased abundance of a number of molecular functions in A. viridis compared with N. vectensis: translation regulator activity (0.46 vs 0.23%), oxidoreductase activity (6.86 vs 4.98%), antioxidant activity (0.07 vs 0.03%), and structural molecule activity (3.24 vs 1.28%). It is noteworthy that some of these molecular functions are involved in symbiotic interactions with zooxanthellae, such as protein binding [24], nutrient transport and antioxidant activity [8]. For example, A. viridis counters the oxidative stress produced by the photosynthetic activity of its symbionts by using a great diversity of antioxidant defences [38]. Some biological processes were also increased in the A. viridis dataset compared with N. vectensis: secretion (0.81 vs 0.36%) and catabolism (0.65 vs 0.26%). The cytoplasmic components were also found to be more highly represented in A. viridis compared with N. vectensis (20.42 vs 9.59%).

All these differences between the present A. viridis dataset and the N. vectensis genome may reflect the crucial role of trophic exchanges between the sea anemone and its dinoflagellate symbionts, as well as specific host adaptations.

Sampling and cDNA library construction

To maximize the diversity of genes expressed in the symbiotic association under both normal and stress conditions, sea anemones were subjected to different controlled stress conditions before RNA extraction and cDNA library construction.

Specimens of the Mediterranean sea anemone, Anemonia viridis (Forskål, 1775), were collected close to Villefranche-sur-mer (France). During an initial acclimatation period of at least 4 weeks, animals were maintained in tanks of running seawater at 17 ± 1°C with a light intensity of 250 μmol quanta m^-2s^-1 (overhead metal halide lamps Philips HQI TS 400W), on a 12 h light/12 h dark cycle (starting at 8 am). Animals were fed once a week. After the acclimatation period, animals were treated and sampled (five tentacles from each specimen, four specimens from each experiment except for the hyperoxia stress) as follows:

The mRNAs were then extracted using Trizol Reagent (Invitrogen), as described in Richier et al. [39], and equal amounts of mRNA from the different conditions described above were combined before cDNA library construction. A. viridis cDNA library was generated at the RZPD (Deutsches Ressourcenzentrum für Genomforschung GmbH, Berlin, Germany) by oligodT, random priming and directionally cloning into pSPORT1. From this library, 50,304 clones were picked, replicated and subsequently sequenced at the Genoscope (French national sequencing centre). High-throughput sequencing of 5'end was performed using a Big-Dye terminator cycle sequencing kit and M13 reverse primer on an ABI-3730 Genetic Analyzer (Applied Biosystems) following the manufacturer's protocol (Genoscope, Evry, France). 41,247 chromatograms were thus generated and further analyzed.

EST processing and analysis

EST sequences were processed using SURF analysis pipeline tools (SURF: SeqUence Repository and Feature detection, developed by the SIGENAE team, Dehais Patrice and Eddie Iannucelli, INRA, Toulouse). Basically, SURF provided an integrated solution, from chromatogram data storage to cloned insert detection, by integrating several dedicated bioinformatic software programs (sequence base calling, vector detection, etc.) in order to produce relevant nucleotide sequences according to base quality and feature detection. The chromatogram files were exported to PHRED for base calling [40, 41]. Cloned insert detection was made according to different detected features (vector, adaptator, poly(A) or poly(T) tails and repeat) and their respectively positions, using third party programs (Crossmatch and RepeatMasker). Only inserts with more than 100 bp, with a Phred score >20, and not belonging to a low complexity area were exported into a fasta format with its corresponding quality file. Additional extremity trimming was made using the "trimseq" command (EMBOSS package). Low complexity regions and repeats were masked using the RepeatMasker program [42]. For this purpose, two different libraries were used: the RepeatMaskerLib (RepBase Update of 2007.09.24, [43]) and a custom library of A. viridis. This custom library was made both by using CENSOR [44] to retrieve publicly available repeats, and by running a BlastN of all ESTs against themselves to identify the most abundant repeat regions.

High quality ESTs (39,939) were then assembled into contigs using the TIGR-TGICL tool [32].

Repeated elements analysis

Putative transposable elements (46 sequences) were first identified based on homology search after BlastX analysis against UniProt KB (E-value of < 1.10^-20). An additional local BlastN search was performed using the EST dataset both as query sequence file and target database. Repeated motif sequences, i.e. repeats occurring more than twice from non overlapping ESTs, were selected as a first screen repeats dataset. These were used, together with the Repbase repeats library, to mask our EST sequences before clustering and assembling (TGICL). The same first screen dataset was then BlastN compared with the assembled database, and repeated motif sequences occurring on more than two different UniSeqs (E-value of < 1.10^-20) were considered as A. viridis repeat sequences.

Gene Ontology annotation

Gene functions were automatically assigned to 39% of the predicted proteins (5,652 UniSeqs). This assignment was based on the identification of InterPro (IPR) domains [45] using InterproScan [46] and the following command line: iprscan -cli -i unisequences.fa -o unisequences.ipr.raw -seqtype n -goterms -iprlookup -format raw. For comparative analysis of IPR domains found in the A. viridis dataset, we also ran the program InterproScan on the predicted N. vectensis proteome dataset. To homogenize the granularity level of annotation between organisms for each non-overlapping set of domains found, we only kept the root domain. We used the hierarchical organization of domains proposed in the "Parent-Child" description available on the EBI public ftp server ftp://ftp.ebi.ac.uk/pub/databases/interpro/ParentChildTreeFile.txt. For example, all CYP proteins which have a P450 domain ([InterPro:IPR002949], [InterPro: IPR002397], [InterPro: IPR008070]) were counted at their root domain (IPR001128).

Comparative genomic analysis

All contig and singleton sequences were compared with several databases, using Blast: the Nematostella vectensis draft genome (Predicted proteins, http://www.jgi.doe.gov/), SwissProt (2008.03), TrEMBL (2008.03) and other ESTs from symbiotic cnidarian species (Acropora millepora, Acropora palmata, Aiptasia pallida, Montastrea faveolata) or from non symbiotic species (Metridium senile).

To confirm the origin of some selected genes, amplifications were performed on 10 ng of genomic DNA from A. viridis epidermal cells (non-symbiotic cells, animal origin [39]), cultured Symbiodinium cells extracted from A. viridis tentacles (non symbiotic cells, symbiont origin), and whole tentacle extracts (symbiotic cells, both animal and symbionts). Primers designed in the experiment are presented in Additional file 6 and were used in 40-cycle PCR reactions. Elongation factor 1 alpha and Elongation factor 2 were used as positive controls for nuclear-encoded genes from A. viridis and Symbiodinium spp, respectively, while psbA (photosystem II protein D1) was used as a positive control for chloroplast-encoded Symbiodinium spp genes. Sequence alignment was done using Multalin [47]. Signal peptide prediction was performed using SignalP [48]. Phylogenetic analyses were done using both MEGA 4.0 [49] and PHYML [50] software.

Accession numbers

The ESTs generated in this study were submitted to dbEST ([GenBank:FK719875–FK759813]). Accession numbers of sequences used in MERP alignment: Am_MERP, Acropora millepora, [GenBank:EZ013381.1]; Mm_MERP, Mus musculus, [REFSEQ:NP_598826]; Hs_MERP, Homo sapiens, [REFSEQ:NP_060019]; Ca_MERP1, Carassius auratus, [GenBank: X14134.1]; Ca_MERP2, Carassius auratus, [GenBank: J04986.1].