A genomic survey of the fish parasite Spironucleus salmonicida indicates genomic plasticity among diplomonads and significant lateral gene transfer in eukaryote genome evolution

The first sequence data from Spironucleus on the genomic level

Variable G+C content along the genome

On the whole, the S. salmonicida genome was found to be G+C poor; the average G+C contents were 35.9%, and 39.3% within the 5070 GSSs and the eleven finished sequences, respectively. However, the eleven finished contigs showed large variation in G+C content, with average G+C contents ranging from 24.9% to 58.3% (named Sp1-11, Figure 1). This probably cannot be explained by the presence within each contig of genes with the higher values; the coding content of Sp2 is 96% and its G+C content is only 31%, while the coding content of Sp8 is 75% and its G+C content is 58% (Additional file 1). In fact, the G+C content sometimes varies drastically within each contig. For example, the average G+C content in the region between 1–5 kbp from the 5' end of Sp1 is around 30%, while between 6–11 kbp its G+C content is almost 70% (Figure 1). Similar sharp shifts in G+C content along the genome are observed in some other contigs (e.g. Sp7, Sp10, and Sp11: Figure 1). These shifts in G+C content do not seem to be correlated with the presence and absence of genes: both G+C rich and poor regions are fairly dense with genes, with the exception of the G+C poor contig Sp4 and the 5' end of Sp9. However, these observations are consistent with the existence of different mutational biases in different regions of the S. salmonicida genome; the G+C contents of both the third synonymous position of protein coding genes and the non-coding regions varies as expected if mutational bias shapes the G+C content [35], with the GC3 values of genes higher than the average G+C content for genes in G+C rich regions and lower in G+C poor regions (Figure 1). A similar pattern is observed for non-coding regions (Additional file 1), and the codon usage analysis indicates large variations of the G+C contents among genes (see below).

S. salmonicida has a compact genome

Analysis of the eleven completely sequenced contigs showed that the coding density varies significantly between contigs (Figure 1 and Additional file 1). Some parts of the genome appear to be very gene-dense, with short intergenic regions. In fact, six cases of putative overlapping open reading frames (ORFs) were identified, similar to earlier findings in G. lamblia [36]. The lengths of the overlapping regions range from 1 to 101 bp; in four cases they are encoded in the same direction, while in the other two cases the 5' regions overlap. Experimental studies are needed to determine which initiation codons actually are used in these putative genes to determine if they are truly overlapping transcribed ORFs that encode functional proteins. Nevertheless, the whole genome does not seem to be packed with genes; no protein coding sequence could be detected in the almost 5-kbp contig Sp4, nor in a long region upstream of the 18S ribosomal RNA gene in Sp9 (Figure 1). In G. lamblia, the genes for ribosomal RNA have been found to be associated with telomeres [37]. However, we were unable to detect any sequence similarity to genes that are associated with telomeres in G. lamblia, or any telomere repeats, in the long region upstream of the rRNA genes in S. salmonicida Sp9.

The transcription of protein-encoding genes in Giardia is atypical for eukaryotes. The promoters are short; less than 70 bp is needed for efficient expression of most genes, even if they are stage-specific [38]. The 5'-UTRs of the transcripts are usually only 1–10 nt, and the genes lack TATA boxes or other cis-acting promoter elements characteristic of typical eukaryotic promoters; AT-rich sequences around the G. lamblia ATG initiation codons determine the transcriptional start site and are essential and sufficient for promoter activity [39]. The S. salmonicida promoters apparently lack TATA-boxes and other typical eukaryotic promoter motifs (Figure 2A). However, AT-rich stretches are found close to the translational start sites, or in the first 50 bp upstream (Figure 2A) and the intergenic regions are very short (Figure 1). Thus, the S. salmonicida transcriptional machinery appears to resemble the 'atypical' machinery previously characterized in G. lamblia [39].

The position of the polyA tail in relation to the termination codon could be mapped in 134 EST clones. The 3' untranslated regions (UTR) of S. salmonicida genes appear to be short; the average distance between the termination codon and the first A in the polyA tail is 13.2 bases and in 122 of the 134 cases the distance is between 9 and 14 nucleotides (Figure 2B). Such very short 3' UTRs have previously been observed in G. lamblia, but this is in stark contrast to most other studied eukaryotes where 3' UTRs have been proposed to play important roles in regulating messenger RNA (mRNA) stability and translation [40]. In G. lamblia, AGTRAAY has been identified as a consensus polyadenylation signal [41]. Our analyses do not identify any strong consensus sequence outside the conserved unique termination codon TGA, although we observe a preference for A in the position two bases upstream of the termination codon, and a T immediately downstream (Figure 2B). Interestingly, the three most conserved nucleotides in the G. lamblia polyadenylation signals are the central TRA nucleotides, which is similar to the situation in S. salmonicida. Furthermore, a preference for a U-rich region close to the polyA tail is also apparent, a pattern conserved in many eukaryotes but not Giardia [41, 42]. It has been reported that a fraction (~20%) of complementary DNAs (cDNAs) from G. lamblia represent polyadenylated sterile anti-sense transcripts that may result from a loose control of transcription [43]. The rarity of such sterile transcripts in our cDNA library (~1% among 2341 EST) indicates that this phenomenon is not universal for diplomonads. Indeed, our data indicate considerable variation in the specificity of the transcription and polyadenylation machineries within diplomonads.

The availability of both EST and genomic sequences was used to identify putative introns. An intron should be detected as a sequence gap in an alignment between the EST and the corresponding genomic sequence because the EST sequences typically represent the spliced version of the gene sequence. This approach was used to identify potential introns using similarity searches with the EST sequences as probes against the GSS and finished contig genomic sequences. 290 out of 885 ESTs matched with an e < e^-10 (i.e., almost all were identical sequences) and manual inspection failed to detect any introns. The S. salmonicida genome appears to be very intron-poor. Given that only a few introns have been detected within the diplomonad G. lamblia [20, 21], a paucity of introns could represent the ancestral state of the diplomonads. However, two introns were detected in one of the first genes sequenced from Carpediemonas [44], an excavate taxon and possible sister lineage to diplomonads. Thus, introns may have been more frequent in a common ancestor of diplomonads and Carpediemonas, and that the ancestral diplomonad lineage experienced genome-wide intron loss, an evolutionary phenomenon that is not uncommon [45]. A basic spliceosomal organisation is present in Giardia suggesting that the spliceosome is ancestral to extant eukaryotes [46]. Although our analyses did not identify any proteins likely to be involved in splicing in S. salmonicida, further sequencing and more detailed analyses are needed to determine whether introns and a functional spliceosomal apparatus are present in S. salmonicida.

Coding capacity of S. salmonicida

In this study we have identified more than 1300 novel genes, distributed among all functional categories, representing a twenty-fold increase of the number of genes with assigned functions identified in the S. salmonicida genome (Table 1). Among these were genes involved in translation: 74 genes encoding ribosomal proteins, many translation factors, and fifteen different tRNA synthetases were identified (Additional file 2). Twenty tRNA genes were identified, which cover 24 of the 63 sense codons allowing for the normal wobble rules for codon-anticodon pairing. No tRNA was identified that decodes UAG or UAA codons, which encode glutamine in Spironucleus, although such tRNA genes have been identified previously [25]. Interestingly, all 64 codons appear to have the potential to code for incorporation of amino acids into proteins in S. salmonicida; both tRNA identification programs used in our analyses identified a putative tRNA with the anticodon UCA, annotated as a selenocysteine (Sec) tRNA. The identification of a putative Sec tRNA, which is a central component of selenoprotein biosynthesis [47], in our data set suggests that Spironucleus is able to use the single stop codon (UGA) to incorporate this rare amino acid into selenoproteins. Indeed, the usage of selenocysteine seems to be widespread feature in protists; Sec tRNAs have recently been identified in Dictyostelium, Tetrahymena, and Plasmodium [48, 49].

The identified proteins within the cellular, environmental information, and genetic information process categories in our data set clearly reflect a coding potential similar to a typical eukaryote, including eukaryotic translation and transcription machineries, many proteins involved in eukaryotic signal transduction pathways, a large family of dynein proteins, and five genes encoding members of the Rab small GTPase family (Additional file 2). This is indeed expected from the current view of diplomonad phylogeny and cell biology (see Background section). Interestingly, only eight Rab genes were found in G. lamblia, while a large number of Rab GTPases were found in Trichomonas vaginalis, most likely related to the apparently more-complex endomembrane system in this parabasalid [50]. Two of the five Rab proteins identified in our survey lack identifiable orthologs in the G. lamblia genome, but branch with T. vaginalis Rab sequences in phylogenetic trees (data not shown), suggesting that S. salmonicida has retained Rab proteins that have been lost in the lineage leading to G. lamblia.

Absences of genes are very difficult to infer from partial genome data, but some general trends may be observed. For example, relatively few enzymes involved in amino acid metabolism were detected (Additional file 2). Aminoacyl-tRNA synthetases are classified into this category, but they are atypical since they are essential for protein synthesis. Furthermore, only a single protein (malate dehydrogenase) associated with the tricarboxylic acid (TCA) cycle was found. However, malate dehydrogenase actually functions in a pyruvate synthesis pathway in G. lamblia [2], rather than in the TCA cycle, suggesting a similar role in S. salmonicida. In contrast, several glycolytic proteins are present in our data set (Additional file 2). This pattern of metabolic proteins is expected from a fermentative phagotrophic heterotroph which has access to organic compounds from its host, and indeed is very similar to the pattern found in Entamoeba histolytica [51]. These similarities between S. salmonicida and E. histoytica are almost certainly due to independent adaptations to such an environment in the two lineages. Indeed, a considerable fraction of the metabolic proteins was found to be more closely related to prokaryotic rather than eukaryotic homologs in the phylogenomic analyses, suggesting acquisition of these genes by gene transfer (see further discussion below).

Does S. salmonicida possess mitosomes?

A new organelle, the mitosome, has recently been identified in G. lamblia, that is probably a remnant of a mitochondrion [19]. Currently, the only known function of mitosomes in Giardia is iron-sulfur (Fe-S) cluster synthesis [52, 53]. The phylogenetic relationship between G. lamblia and S. salmonicida, the likely mitochondrial origin of the organelle, together with mitochondria and their derived organelles currently being thought to be universally present in all extent eukaryotes [10] strongly suggest that the common ancestor of the two species contained an organelle with mitochondrial ancestry. Here we identified two strong cases of candidate mitosomal proteins in S. salmonicida, a chaperon GroEL (or Hsp60 or Cpn60) a key protein for mediating protein folding in mitochondria [28] and a cysteine desulfurase (called Nifsp in Saccharomyces cerevisiae), a key enzyme of the Fe-S cluster synthesis pathway [54] (Additional file 2). Orthologous proteins of Hsp60 and Nifsp are known to localize and function in S. cerevisiae mitochondria and G. lamblia mitosomes [52, 53]. Localization studies of these two proteins should indeed be very useful to investigate whether S. salmonicida contains mitosomes.

A full-length ORF encoding a dynamin-like protein that is found on contig Sp3 (Figure 1) is another protein potentially linked with the mitosome. The G. lamblia genome encodes a single dynamin-like protein and phylogenetic analyses recover the two diplomonad sequences as monophyletic (data not shown). The function of the G. lamblia single dynamin-like protein is currently unknown. Yet, a single dynamin-like protein encoding gene is also found in three kinetoplastid genomes (two Trypanosoma and one Leishmania) [55] and the microsporidium Encephalitozoon cuniculi [56], whereas most eukaryotic genomes encode several dynamin-like paralogues that function either in membrane trafficking or organelle division (mitochondria and plastids) [57]. Interestingly, the single dynamin-like protein from Trypanosoma brucei was shown to be involved in mitochondrial division [58]. Hence the single dynamin-like proteins in Spironucleus, Giardia and Encephalitozoon could all be involved in mitosome division.

A large family of cysteine-rich proteins

The predicted amino acid sequences of 149 genes were found to contain more than 10% cysteine and were classified as cysteine-rich proteins. Most of these were annotated as conserved hypothetical proteins, although some could be assigned to a functional category (Additional file 2). Most of the cysteine residues were found as CXXC motifs (Figure 3 and data not shown). Such an arrangement is similar to the large protein family of variant-specific surface proteins (VSP) found in G. lamblia [59], which is estimated to comprise 2.4% of the genome [60]. However, less than 10% of the cysteine-rich proteins in S. salmonicida showed highest sequence similarities to G. lamblia cysteine-rich proteins, and the two conserved motifs of G. lamblia VSP proteins, CRGKA and GGCY [59], were not found within any of the amino acid sequence of S. salmonicida cysteine-rich proteins. Indeed, in similarity searches the majority were most similar to putative proteins discovered in the ciliate Tetrahymena thermophila genome project (data not shown) [61]. These observations indicate that diplomonads vary greatly in their cysteine-rich proteins; the gene families indeed appear to have expanded independently in the Giardia and Spironucleus lineages. The biological role of VSP proteins is not well understood, but VSP proteins are immuno-dominant in G. lamblia and expressed on the cell surface [59]. Among the 149 cysteine-rich proteins coding sequences two are likely to represent full-length ORFs. These are strong candidate surface proteins since they possess a transmembrane (TM) domain and a cysteine-rich domain made of furin-like and/or epidermal growth factor (EGF)-like domains typical of surface or secreted proteins, a feature shared with G. lamblia VSPs (Figure 3). A total of 14 sequences encoding partial proteins with cysteine-rich domains possess similar structural organization as the one shown in Figure 3, i.e. they possess a TM domain with the cysteine-rich domains likely to face the extra-cellular milieu. No TM domains could be found in the remaining cysteine-rich proteins within our dataset. However, in several cases these partial cysteine-rich proteins show high similarity to cysteine-rich proteins containing TM domains, suggesting that they could represent partial surface protein sequences (derived from EST or GSS). This observation is consistent with the genome of S. salmonicida encoding a large gene family of cysteine-rich variant surface proteins, as established for the VSP protein family in G. lamblia where only one of the 150 VSP genes is expressed per cell. In G. lamblia the VSP gene family displays antigenic variation and new VSP genes are induced at a relatively high frequency, which is important for escaping the host's adaptive immune system [2, 59, 62]. Further sequence data and analyses, as well as expression and cellular localization studies are needed to understand the function of this large protein family in S. salmonicida.

Insights into the molecular basis of S. salmonicida pathogenicity

Spironucleus species have been described as a cause of disease in hosts that include birds, fish and mice [1]; S. salmonicida, for example, is a problem for the fish industry [3, 4]. However, essentially nothing is known about virulence factors in Spironucleus. In healthy fish, the parasite Spironucleus vortens is commonly found in the flagellated stage (trophozoites) in the lumen of the upper intestine, where it remains attached to the intestinal mucosa, controlled by the mucosal immune system of the host. This stage of infection can cause diarrhea, but villus atrophy, ulceration or inflammation are thought to be absent [63], which is similar to observations during G. lamblia infections [62]. Antigenic variation among the identified cysteine-rich proteins in S. salmonicida could potentially protect the parasite against secretory antibodies at the mucosal surfaces. Nitric oxide (NO) and reactive oxygen species (ROS) are important factors in the host's protection against mucosal pathogens [62]. S. salmonicida encodes at least two rubrerythrins, two A-type flavoproteins and arginine deiminase (Additional file 2). These proteins, which are coded by genes putatively acquired from prokaryotes via gene transfer (see below), protect other microbes against ROS and NO [62, 64, 65], and similar roles are likely in S. salmonicida.

In the laboratory, Spironucleus muris is transmitted as cysts in fecal material like Giardia [66]. Cysts of the genus Spironucleus share many morphological features with those of Giardia including the presence of two to four nuclei, flagellar axonemes, and a distinct cyst wall. Cysts of S. muris even display the same immunostaining as Giardia cysts when labeled with antibodies specific for Giardia cyst wall [66]. We identified one homolog of G. lamblia cyst wall proteins, and two enzymes involved in the production of the sugar components of the cyst wall (glucosamine-6-phosphate isomerase and UDP-N-acetylglucosamine pyrophosphorylase – Additional file 2). The presence of these three genes indicate that S. salmonicida may also be able to form cysts, an interesting possibility since nothing yet is known of how S. salmonicida is spread and how it initiates infection in fish. These genes are stage-specifically expressed in the G. lamblia life cycle, suggesting that studies of their expression might identify regulatory regions of stage-specific genes in Spironucleus. Characteristic features of Giardia cysts include fragments of the adhesive disk but this structure, involved in attachment of Giardia, is conspicuously absent from Spironucleus cells. None of the genes encoding disc-specific proteins in G. lamblia were identified in this genomic survey, consistent with the absence of this structure in S. salmonicida.

In contrast to G. lamblia, trophozoites of S. salmonicida can enter the blood causing systemic and organ infections [67]. Cysteine proteases are essential for pathogenicity and invasion of the intestine of E. histolytica [68]. Therefore the several cysteine proteases that were found in our genomic survey are potential virulence factors possibly involved in invasion/tissue destruction by S. salmonicida. Secretion of cysteine proteases was indeed recently shown to occur from trophozoites of G. lamblia upon interaction with epithelial cells [69]. Cysteine protease activities are also important for excystation and encystation of Giardia [70, 71], suggesting yet another potential role for these enzymes in S. salmonicida.

Distinct codon usage patterns in highly versus weakly expressed genes

Selection on codon usage in highly expressed S. salmonicida genes

To investigate the codon usage in more detail, we plotted the effective number of codons calculated with a method that accounts for variations in the nucleotide composition (N_c') [75] against the G+C content in the third synonymous position (GC3_s), and performed correspondence analyses on the relative synonymous codon usage, for both the Spironuclues and Giardia datasets (Figure 4). The highly expressed S. salmonicida genes are centered in the plot with GC3_s values around 0.4, and show lower N_c' values than the weakly expressed genes with similar GC3_s values, indicative of non-random synonymous codon usage (Figure 4A). For G. lamblia, the putatively highly expressed genes exhibit slightly higher GC3_s and lower N_c' values than the majority of the genes (Figure 4B). However, a separation of the genes into different categories can be observed: the vast majority of the weakly expressed genes have GC3_s values between 0.4 and 0.7, while the highly expressed genes have slightly higher values (Figure 4B). Correspondence analyses showed that the G+C content in the third synonymous position is the main determinant of the codon usage in both S. salmonicida and G. lamblia; GC3_s showed a strong correlation with the first axis (R² = 0.969 and 0.944, respectively: data not shown). Expression levels strongly influence the second axis in the correspondence analysis for S. salmonicida; the highly expressed genes are located at the top separated from the majority of the weakly expressed genes, when axis 2 is plotted against axis 1 (Figure 4C). No such trend is observed for G. lamblia, for which the putative expression levels seem mainly to be correlated with axis 1 (GC3_s) (Figure 4D). These analyses clearly indicate that there is selection on codon usage in S. salmonicida genes classified as putatively highly expressed, and suggest that their G. lamblia homologs also are under selection (Figure 4 and Table 2). Translation efficiency has previously been suggested as the cause of the selection in highly expressed genes in G. lamblia, based on a much smaller data set [74], even though only estimates of expression levels were available for a few of the genes. Our study provides, for a much larger dataset, a connection between the subset of genes with a strong codon bias, and an indirect indication of expression levels (cDNA abundance), corroborating the earlier conclusion and extending it to the distantly related diplomonad S. salmonicida.

A shift in codon usage within diplomonads

Interestingly, the favored codons differed between S. salmonicida and G. lamblia for ten amino acids, while the same codon is most abundant for eight (two amino acids are encoded by only a single codon) (Table 2). In nine of the cases where there is a difference between the two species, S. salmonicida preferentially uses a codon ending with U, while G. lamblia uses a C-ending codon in eight cases and a G-ending codon once. The tenth case is a two-fold degenerate amino acid where S. salmonicida and G. lamblia prefer A- and C-ending codons, respectively (Table 2). Thus, in general S. salmonicida favors A+U rich codons, while G. lamblia prefers codons with G or C in the third position in putatively highly expressed genes. There is a distinction between the base compositional biases in G. lamblia and S. salmonicida; the average G+C content of the S. salmonicida sequences in our dataset is 36%, and the majority of the weakly expressed genes have GC3_s values around 20%, while the weakly expressed G. lamblia genes have GC3_s around 50% close to the genomic G+C content of G. lamblia [60] (Figure 4A, B). This is an expected pattern if mainly mutational processes have shaped the codon usage [72, 73]. However, the genomic G+C content differences seem to have influenced the codon preferences in the two diplomonad lineages in the highly expressed genes as well. In the absence of information about codon usage in a closely related outgroup, it is difficult to determine in which of the two lineages reassignment of the optimal codons has occurred. At any rate, the usage of an alternative genetic code in Spironucleus but not Giardia [26] indicates that codon usage have been remodeled in the lineage leading to Spironucleus, at least for the stop and glutamine codons. The canonical stop codon UAA is indeed the most common codon for glutamine in our dataset (Table 2), as expected in a G+C poor genome if the codon usage is close to mutational equilibrium [72, 73].

Unusual codon usage in a subset of S. salmonicida genes

We observed large variation among GC3_s values for weakly expressed genes in S. salmonicida. Unexpectedly, the N_c' values for genes with GC3_s above 0.5 are negatively correlated with GC3_s (Figure 4A); if only the background nucleotide composition shaped the codon usage for these genes, the N_c' values should be close to the maximum. This pattern indicates that there is selection on the codon usage in the group of weakly expressed genes with high GC3_s values. However, selection for the set of optimal codons for highly expressed genes does not explain this pattern, as they include both G+C rich and G+C poor codons (Table 2). Also, the large spread of the genes to the left on axis 1 (which is strongly correlated with GC3_s – Figure 4C) argues against selection for an alternative set of optimal codons in these genes. Furthermore, the deviating GC3_s values in these genes are probably not due to recent introductions into the S. salmonicida genome from heterogeneous sources by lateral gene transfer (LGT). Although some genes with high GC3_s values do show indications of gene transfer in the phylogenomic analyses, the majority do not (Figure 4A, C). In fact, many of these genes have homologs in Giardia, indicating that they were present in the common ancestor of the two lineages, strongly suggesting that the deviated GC3_s and N_c' values arose during the evolution of diplomonads on the branch leading to S. salmonicida.

Analyses of the finished contigs also yielded unexpected results; the base compositional bias in segments of the S. salmonicida genome seems to be markedly towards G+C, while the largest part of the genome shows a compositional bias towards A+T (Figure 1). This pattern is in contrast with the A+T rich genomes of D. discoideum, E. histolytica, and Plasmodium falciparum, which show a relatively uniform nucleotide composition across all chromosomal regions [76, 77]. The genes with deviating G+C content are correlated with these genomic regions of high G+C content (Figure 1). Furthermore, genes in G+C rich regions have higher GC3_s values than the average G+C content for the genes (their bars are located above the sliding window analysis in Figure 1). Conversely, as expected, genes in G+C poor regions have low GC3_s values. Thus, genes with unexpected codon usage (i.e. high GC3_s values) are located in specific regions of the genome. However, neither mutational or selection scenarios can easily explain their non-random codon usage, nor the fluctuations in G+C content along the genome. Heterogeneity in overall G+C content (and/or GC3s) across the genome has been also observed in Saccharomyces cerevisiae [78] and in vertebrate genomes (e.g. isochores) [79]. For S. cerevisiae (a unicellular eukaryote with a compact genome, that is more directly comparable to S. salmonicida than vertebrates) a number of possible explanations have been advanced such as regional variation in propensities for mutation or recombination and partitioning of the genome into "distinct replicational and transcriptional domains" during specific cell cycle phases that may experience different chemical environments [see [78] and references therein]. Longer genomic fragments, as well as comparative data from more-closely related diplomonads, will be needed to characterize the pattern further in S. salmonicida and to tease apart possible causes for this peculiar phenomenon.

Phylogenomic analysis reveals frequent lateral gene transfer

We and others have previously studied the occurrence of LGT in diplomonads on a single-gene basis [16, 17, 29–31], [80–82]. Here we present a systematic phylogenomic analysis of the complete EST and GSS data sets with the goal of identifying S. salmonicida genes that have been involved in LGT. We used the PhyloGenie package [83] to assemble aligned amino acid data sets automatically from databases including both published sequences and data from ongoing eukaryotic genome projects (see Methods). Among the 711 genes with three or more identified homologs for which more than 100 amino acid positions could be aligned, 84 were retained as genes putatively involved in LGT affecting the S. salmonicida lineage, on the basis of unexpected phylogenetic positions (Figure 5, Additional files 3, 4, 5, 6). Cases where two S. salmonicida genes branched as a monophyletic group in the phylogenetic tree to the exclusion of other sequence were interpreted as gene duplication event after the gene transfer. Taking this into account, the 84 genes putatively involved in LGT corresponded to 68 unique gene transfer events at most; multiple genes may indeed have been introduced in a single event in some cases.

In principle, contamination in our libraries could falsely indicate a pattern of frequent gene transfer affecting the S. salmonicida genome. However, we strongly believe this not to be case. 57 of the 84 putatively transferred genes formed a monophyletic group with G. lamblia homologs in the phylogenetic analyses indicating that these genes were present in the common ancestor of the two species (Additional Files 3, 4, 5, 6). Furthermore, on average five in-frame UAG and UAA sense codons (which are termination codons in all characterized prokaryotes [84]) were detected in 21 of the remaining 27 putatively transferred genes, strongly suggesting a true S. salmonicida origin also for these genes. Four of the six genes without UAG or UAA codons formed a monophyletic group with a S. salmonicida gene which does utilize such codons, indicating a likely origin by gene duplication. We cannot formally exclude contamination on the basis of presence of alternative codons for the remaining two sequences (gTor1213bT7 and gZap260bT7, trees #59 and #64 in Additional file 6), but notice that the absence of in-frame UAG and UAA codons in these genes could likely be due to a recent introduction from a donor utilizing the canonical genetic code. At any rate, together these observations indicate that the phylogenetic pattern we interpret as putative LGT events are unlikely to result from substantial contamination of our libraries.

The strength and type of support for the inferred gene transfer events varied between genes in the phylogenomic analysis. Some Spironucleus genes were strongly associated with genes from distantly related organisms. For example, the S. salmonicida carotenoid isomerases were found to be nested within proteobacterial sequences (Figure 5D). Such anomalous positions of the diplomonads genes strongly suggest gene transfer events. Many datasets contain sequences from only one or a few eukaryotic taxa other than diplomonads (Figure 5A, C, E), and in some, diplomonads are the only eukaryotes present (Figure 5B–D). Although the separation of diplomonad from the other eukaryotic sequences is sometimes weak (i.e. Figure 5E), we suggest that LGT is a more likely explanation for their origin than differential gene losses and/or phylogenetic artefacts, as we have argued previously for genes with a patchy phylogenetic distribution among eukaryotes [80].

There are indeed examples of inferences of LGT based on unexpected patterns of phyletic distribution that have previously been shown to be wrong. Re-analyses of the putative bacteria-to-vertebrate gene transfer events in the human genome using phylogenetic analyses instead of similarity searches showed that more parsimonious alternative explanations existed for the vast majority, if not all, of the reported cases [85]. Phylogenetic analyses may also lead to false positive interpretations of LGT. For example, phylogenetic analysis indicated an LGT origin of the T. vaginalis hydrogenosomal NuoF protein distinct from the mitochondrial homologs [86], a conclusion shown likely to be incorrect upon re-analysis using more realistic models of sequence evolution [87]. On the other hand, re-analyses using extended taxon sampling and more advanced models of sequence evolution of four genes putatively affected by prokaryote-to-protist LGT resulted in support for a larger number of inter-domain gene transfer events, rather than a fewer number as expected if differential gene loss, and/or phylogenetic artefacts caused the initial indications of gene transfer [30, 80]. Obviously, inferences about LGT from the phylogenetic analyses depend on the accuracy of the phylogenetic method, as well as on the breadth of organismal representation in sequence databases. Therefore, it is inevitable that the list of putative LGTs will include both false positives and false negatives (Additional files 3, 4, 5, 6). Nevertheless, our phylogenomic analysis indicates that vertical inheritance combined with gene duplication and gene loss is rather unlikely to have produced the observed phylogenetic relationships and patchy taxonomic distributions and that gene transfer is the simplest hypothesis currently available to explain these trees. Hence our analyses indicate that LGT has been an important mechanism in the evolution of the S. salmonicida genome.

Acquisitions of prokaryotic metabolic genes in the evolution of protist genomes

The vast majority of genes inferred to be involved in gene-transfer events are metabolic genes (74%) or encode conserved hypothetical proteins (20%) (Additional file 3). This bias is not unexpected, as a higher rate of transfer for metabolic genes compared to informational genes has been observed for prokaryotes [88, 89] and, on a smaller scale, in diplomonads [30] and E. histolytica [51]. In the latter study it was suggested that this protist parasite had expanded its metabolic repertoire via gene acquisitions from prokaryotes present in the environment, and the same may apply for S. salmonicida. For example, amino acids have been identified to be a source of energy in diplomonads via the arginine dihydrolase pathway [2] especially under limited oxygen conditions [90]. Among eukaryotes, this pathway has been detected only in Trichomonas and diplomonads. The genes for two enzymes in the pathway, arginine deiminase and ornithine transcarbamoylase, appear to have been acquired from prokaryotes in a common ancestor of diplomonads and parabasalids, whereas carbamate kinase may have distinct LGT origins in the two groups (Figure 5A and trees #10 and #30 in Additional file 4). Similarly, S. salmonicida encodes rubrerythrins and A-type flavoproteins, proteins putatively involved in the protection against the host's defense (see above) [62, 64, 65] that are also found in anaerobic prokaryotes and eukaryotes including Entamoeba (Figure 5E and Additional files 3, 4, 5, 6) [51, 80]. The phylogenetic distributions of these genes suggest a role related to an anaerobic lifestyle, consistent with the hypothesis that these acquisitions provided metabolic advantages to the recipient protist lineages. Indeed, distantly related lineages that live in the same environment appear to frequently transfer genes [89]. Increased taxon sampling and more detailed studies of the ecology and metabolism/biochemistry of the organisms represented by sequences in our trees will likely recognize shared environments between the donor and recipient lineages for additional acquired metabolic genes (Additional file 3).

Continuous exchange of genes with all domains of life

In this study we have identified 84 genes of putatively lateral origin, and in previous phylogenomic studies 96 E. histoytica genes, almost 50 kinetoplastid genes, 24 in Cryptosporidium parvum and 148 genes in anaerobic ciliates were found to represent candidate LGT genes [51, 55, 91, 92]. Thus, numerous inter-domain gene transfer events have been identified in four divergent protist lineages – Spironucleus, Entamoeba, kinetoplastids (Trypanosoma and Leishmania) and anaerobic ciliates – supporting the idea that gene acquisition from prokaryotic organisms is a common evolutionary mechanism in unicellular eukaryotes [22, 93–95]. In all these studies, metabolic genes have been found to be the most common functional category among the genes implicated in LGT genes. This observation, which is consistent with the complexity hypothesis [88], could indicate that metabolic adaptation is a selective force for inter-domain gene transfer events from prokaryotes to protists [96, 97].

G. lamblia homologs that form clades with the S. salmonicida sequences were found in 49 of the 68 putative LGT events, while no closely related G. lamblia homologs were found in the remaining 19 cases (Figure 5 and Additional file 3). Thus, at least 72% of the putative LGT events happened before the divergence of the Spironucleus and Giardia lineages. The remaining 19 events could have taken place in the S. salmonicida lineage after the divergence from the lineage leading to G. lamblia, although gene losses in the latter could also explain the observed pattern. If these represented relatively recent gene transfers, they should show a codon-usage pattern distinct from that of other S. salmonicida genes. Some in fact do show deviant codon usage indices (Figure 5A, B) suggesting a more recent acquisition compared to the other cases. However, there are many other S. salmonicida genes with deviant codon usage that do not show any indication of LGT (Figure 4A, B), suggesting that codon usage pattern alone is a poor indicator of LGT in S. salmonicida, as shown in other systems [98]. In 12 cases where both Spironucleus and Giardia sequences are present, T. vaginalis sequences are found as a sister group to the diplomonad cluster. This pattern has previously been observed for candidate laterally transferred genes [16, 17, 80] and almost certainly reflects a common ancestry of diplomonads and parabasalids to the exclusion of other sampled eukaryotic lineages, although the possibility of gene transfer between the two lineages cannot be formally excluded. However, a diplomonad-parabasalia relationship has recently been shown to be robust in several phylogenetic analyses of concatenated protein sequences [13–15]. As some acquired genes are unique to S. salmonicida, others shared with G. lamblia, and some are also shared with T. vaginalis, this strongly suggests that genes have been acquired by LGT in these excavate lineages continuously throughout their evolutionary history. Further characterization of the phylogenetic distribution of these genes may provide information about organismal relationships within Excavata [99].

In the four previous phylogenomic studies of LGT in protists [51, 55, 91, 92] only prokaryote-to-eukaryote transfer events were examined. Our selection procedure allowed us to detect additional evolutionary scenarios where a prokaryotic gene is transferred to a eukaryote, and then possibly further transferred to a second eukaryotic lineage. For example, the arginine deaminase gene was likely acquired from a prokaryote by an ancestor of diplomonads and parabasalids, and then transferred to the lineage leading to Dictyostelium (Figure 5A). Seventeen genes were putatively affected by intra-domain LGT, corresponding to 13 unique events (Additional file 3). Twelve of these cases are exchanges between the lineage leading to S. salmonicida and amoebozoan lineages; five cases involve gene exchange with the Dictyostelium lineage, six cases with the Entamoeba lineage, and a single case with the Mastigamoeba lineage (Figure 5A, Additional files 3, 4, 5, 6). In sharp contrast, only a single intra-domain gene transfer event representing other eukaryotes was detected, although opisthokonts are comparatively well sampled with genome sequences, and several diverse and complete, or nearly complete, chromalveolate sequences were present in our data sets (Additional files 3, 4, 5, 6). It is noteworthy that only genes that have been acquired from prokaryotes and have subsequently been transferred within the eukaryotic domain will be detected using our selection procedure. Additional eukaryote-to-eukaryote transfers may actually have occurred among the genes that did not show any indication of prokaryote-to-eukaryote transfer. A few examples of gene transfers between unicellular eukaryotic organisms have indeed been published recently [22, 80, 100, 101]. The 13 putative intra-domain transfers identified here further suggest that LGT between protists is an evolutionary mechanism that should not be neglected.

Sources of RNA and DNA, library construction, and sequencing

S. salmonicida (ATCC 50380) [6] (previously known as S. barkhanus [5, 105]) was grown in axenic culture following the ATCC protocol. Messenger RNA from approximately 10⁸ cells of S. salmonicida was isolated using the Dynabead mRNA system (Dynal), and cDNA was synthesized and cloned into the lambda Uni-Zap XR vector (Stratagene), according to the manufacturer's instructions. This procedure requires that the polyA tail is present in the mRNA to be cloned. Infection of SolR cells was carried out and plated on selective LB media containing carbenicillin, isopropyl β-D-1-thiogalactopyranoside (IPTG) and 5-bromo-4-chloro-3-indolyl-beta-D-galactopyranoside (X-gal) for blue/white selection. Positive colonies were picked manually and plasmid DNA was purified using Perfectprep Plasmid Isolation Kit with the Perfectprep Plasmid 96 Vac system (Eppendorf). The DNA was quantified by electrophoresis on a 1% agarose gel. Sequencing was carried out on an ABI 377 (Applied Biosystems) using the ABI Big Dye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems) and the SK primer.

Genomic DNA was purified using standard protocols. A genomic DNA library was constructed for GSS sequencing. The DNA was partially digested using Sau 3AI and ligated to lambda Zap Express vector pre-digested with BamH I (Stratagene) according to the instructions. Mass in vivo excision of the pBK-CMV phagemid from the Zap Express vector was performed using XL-1 Blue MRF' cells and ExAssist phage with the XLOLR strain, and the cells were plated on LB with kanamycin, IPTG and X-gal for blue/white selection. Cells were picked and the plasmids were purified and sequenced as described above using the T3 and T7 primers. Roughly half of the sequencing was done from this library. However, since it was revealed that the insert sizes of this library were too short for efficient random sequencing, a second GSS library was constructed using 20 μg genomic DNA. The DNA was physically sheared using a nebulizer, blunt ended, ligated to the insert DNA and cloned into plasmid vectors using the TOPO Shotgun Subcloning Kit (Invitrogen). One Shot chemically competent cells (Invitrogen) were transformed with an aliquot of the ligation mix and spread on LB plates containing kanamycin and X-gal for blue/white selection. Positive clones were manually picked and plasmids were purified using 96-well Plasmid Preparation Kit (Millipore). The plasmids were quantified and sequenced as described above using T3 and T7 primers.

Finally, S. salmonicida DNA was partially digested with Bam HI and cloned into the lambda DASHII vector (Stratagene) according to the manufacturer's instructions, to make a library with larger inserts (8–28 kbp). The library was grown in MRA cells and spread on eight 150-mm plates each containing approximately 70,000 plaques. Three plaques were purified and the inserts were amplified by PCR, nebulized and the fragments were shotgun cloned into puc18 Sma I/BPA plasmids (Amersham Pharmacia). The ligations were transformed into XL2-Blue MRF' ultracompetent cells. Plasmids were purified and sequenced as described above using M13f and M13r primers.

Sequence analysis

Sequencing reactions were initially screened manually for exclusion of unsuccessful reactions. Sequences from 2408, 3204 and 3684 clones from the Uni-Zap XR EST, the Zap Express GSS and TOPO Shotgun Subcloning GSS libraries, respectively, passed this initial screen. The average read lengths of these were 540, 399 and 366 bp, respectively.

The 2408 EST sequences were screened against the vector sequences and the Escherichia coli genome sequence using the Phred software [106] to check for contamination of DNA from the library host, and clustered using the Phrap assembler [106] with the default settings. The contigs were quality trimmed manually, while the singletons were trimmed using the Phred software [106] with the default quality cut-off. Hybrid clones were identified as containing stretches of nine or more As within the insert, and only the 5' end of the hybrid was retained. 27 clones with inserts in the opposite direction were identified as starting with a stretch of Ts (15), or as having significant matches in the wrong direction to known genes (12); these were reversed. Finally, nine EST sequences were trimmed according to putative frame-shifts detected in similarity searches. After removal of contigs shorter than 100 bases 884 sequences (298 contigs and 586 singletons) corresponding to 502589 unique bp and based on 2341 clones remained for further analysis. 20 and 11 of these overlapped with published genes or genes from the finished contigs (see below), respectively, and were excluded, leaving 853 ESTs for further analyses.

300 ESTs ended with a stretch of eight or more As, indicative of mRNA polyadenylation, were identified. Of these, 134 were annotated as coding based on sequence similarity to genes in the databases (see below). The region around the 3' end of the gene was aligned based on the position of the termination codons (UGA), and a sequence logo was created using WebLogo [107].

The 6888 GSS sequences were trimmed using the Phred software [106] with the default quality cut-off, and the vector sequences were masked and removed. All sequences were screened against the E. coli genome sequence. The 5275 GSS chromatograms that passed this screen (no host contamination was detected), and represented sequences longer than 100 bp, were included in the assembly. The average read lengths of the clones from the Zap Express GSS and TOPO Shotgun Subcloning GSS libraries that passed this screening were 543 and 472 bp, respectively. The Phrap assembler was used with the default settings, yielding 1008 contigs consisting of two or more sequences and 2208 singletons. Some of the longer contigs were selected for complete sequencing. Regions that required additional sequencing were identified within these contigs, and amplified from genomic DNA using PCR with primers designed based on the GSS assembled sequence. The PCR was performed under standard conditions using the following parameters: denaturing at 95°C for 5 min, followed by forty cycles of denaturing at 94°C for 30 s, annealing at 48°C for 1 min, and extension at 72°C for 2 min, and a final extension at 72°C for 10 min. To extend the contigs, sequence gaps covered by clones were identified and amplified using PCR with specific primers based on the obtained shotgun sequences. The PCR products were purified using the Qiaquick PCR Purification Kit (Qiagen) and sequenced using the PCR primers as described above. This procedure resulted in eight continuous stretches of genomic sequence between 2365 and 9005 bp in length. Together with the three sequenced lambda clones these yielded 80504 bp of finished genomic sequence in eleven contigs which were covered by at least one-fold high quality sequence in each direction (Additional file 1). 205 of the 5275 GSS sequences overlapped with the completely sequence contigs and were excluded from further analysis. The remaining 5070 sequences were treated as single reads, and covered 2539160 bp of non-unique sequence after quality trimming.

To estimate the amount of unique high-quality sequence within our survey project, we performed assemblies with only the quality trimmed sequence using the Staden package [108]. The 5275 GSS clones yielded 2070953 bp in 3566 contigs with one or more sequences, and the addition of the EST clones gave 2549915 bp of unique high-quality sequence. Lander and Waterman proposed that the number of "islands" (E) in a genomic project with perfectly representative libraries with equal length of clones (sequences) is E = Ne^-cσ [109], where N is the number of clones sampled, c is the redundancy of coverage = LN/G, σ = 1 - T/L, L is the clone length (500.8 bp in our case), G is the haploid genome size in bp, and T is the amount of overlap in bp needed to detect overlap (20 bp). G is the only unknown parameter, and the equation can be solved: G = N(L-T)/ln(N/E). Using the information from our assembly (N = 5275; L = 500.8 bp; T = 20 bp; E = 3566), the genome of S. salmonicida was estimated to be 6.5 Mbp in size.

Databases

All databases used in the analyses were downloaded in January 2005. In addition to the non-redundant protein database downloaded from the National Center for Biotechnology Information (NCBI) [110], protein databases from various nearly completed genome sequencing projects of diverse eukaryotes were downloaded: T. vaginalis (parabasalid), Tetrahymena thermophila (ciliate), Trypanosoma brucei, and Trypanosoma cruzi (kinetoplastids) from the Institute for Genomic Research (TIGR) [111]; Dictyostelium discoideum (mycetozoan) from dictyBase [112]; Chlamydomonas reinhardtii (green alga), Phytophthora sojae (oomycete), and Thalassiosira pseudonana (diatom) from the DOE Joint Genome Institute (JGI) [104]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were downloaded from the Kyoto University Bioinformatics Center [113] for functional annotation purposes.

Gene identification

Various approaches were used to identify coding genes among the sequences. Similarity searches using BLASTx, version 2.2.13 [114], with the default settings with the nucleotide sequences from the single-read GSS and the EST contigs against the amino acid sequence databases were performed. 1179 single-read GSS sequences and 502 EST contigs yielded matches with E values better than e^-5, which were considered a significant indication of a coding sequence. Using the indication of the frame and direction from the alignment between the query and database sequence in the BLAST result files, the putative coding DNA and amino acid sequences were extracted from these sequences. The 3891 single-read GSS sequences and 380 EST contigs that failed to yield a match better than e^-5 were not analyzed further.

To identify coding sequences in the finished contigs, similarity searches against the amino acid sequence databases and the EST sequences were performed using the same cutoff as above. Similarity searches on the nucleotide level against the non-redundant nucleotide databases were performed at NCBI to identify non-protein coding sequences. Finally, open reading frames longer than 450 bp with the expected GC pattern in the three coding positions were annotated as putatively coding and included in subsequent analyses. This procedure identified 38 protein-coding genes in the eleven finished contigs.

tRNA genes were detected using two programs: tRNAscanSE version 1.23 [115] with the maximum sensitivity (-C option), and ARAGORN version 1.1 [116] with the default parameters. 20 tRNA genes were identified with both methods, all within GSS clones. For 19 of these, both programs assigned the same anticodon, while one differed in anticodon identity (see Results). Two tRNAs were found to encode putative introns. No tmRNA was found. Finally, similarity searches at the nucleotide level identified two genes encoding ribosomal RNA within one of the finished contigs.

Creation of a non-redundant dataset of S. salmonicida protein-coding genes

In order to explore the coding potential of S. salmonicida, we created a non-redundant dataset of S. salmonicida genes that contained only the longest amino acid sequence from each gene among the GSS, EST, finished contigs and previously published genes. Three of the 30 previously published genes were represented by long ESTs and were excluded. This procedure resulted in538 unique protein-coding genes. 175 of the 1179 genes identified within the GSSs were excluded because they were represented by longer ESTs or published sequences, while 30 were excluded because the extracted ORF was shorter than 100 bp. Furthermore, 176 of the remaining 974 genes identified within GSS sequences were excluded because they overlapped with a longer identical gene within the GSSs, leaving 798 unique protein coding genes identified within the GSS data. Thus, in total 1335 unique protein-coding genes were identified and included in further analyses.

Codon-usage analysis

1153 of the 1335 unique protein-coding sequences were 300 bp or longer and were included in the codon-usage analysis. Codon-usage indices were calculated for each gene. N_c is a measure of the effective number of codons used in a gene [117]. A modified version, N_c', has been developed which also account for the background nucleotide composition of the gene; this is advantageous in situations where the composition varies among the genes analyzed [75]. The N_c' values were calculated using the software INCA [118]. The GC3_s values – the frequency of G+C in synonymously third codon positions (i.e. Met, Trp and termination codons are excluded) – were calculated using the program CodonW [119]. The variation of codon usage among genes was explored using the correspondence analysis tool within CodonW. To avoid identification of trends in codon usage due to biased amino acid usage among the genes, the correspondence analysis was performed on the relative synonymous codon usage (RSCU) values for each gene. The RSCU value for a codon is the observed frequency divided by the frequency expected if all synonyms codons for that amino acid were used equally. This should remove the effects of amino acid composition on codon usage. RSCU values close to 1 indicate a lack of bias, while much higher and much lower values indicate preference and avoidance of that particular codon, respectively. The correspondence analysis plots genes according to their RSCU values in a 61-dimensional space, and then identifies the major trends as the axes through this multidimensional hyperspace which account for the largest fractions of variation among genes.

To examine the variation of codon usage within diplomonads, similarity searches of the 1153 S. salmonicida amino acid sequences were performed against the G. lamblia protein sequences. 533 S. salmonicida genes gave matches with E values < e^-20, corresponding to 438 unique genes. These were considered homologs (not necessarily orthologs) and the nucleotide sequences were obtained for codon-usage analysis as described above.

The number of times each gene was detected in the EST library was used as an indication of expression level in the cell. 35 genes were represented by more than ten individual EST sequences and were considered highly expressed, while 682 genes were not detected in the EST data and were considered weakly expressed. Codon-usage tables were calculated for each of these groups, both for S. salmonicida and the G. lamblia homologs. G. lamblia genes represented by several S. salmonicida genes were assigned the highest number of ESTs among the genes.

Gene annotation

The Kyoto Encyclopedia of Genes and Genomes (KEGG) project [113], a bioinformatic knowledge base for systematic analysis of gene functions [120], was used for gene annotation. The KEGG system links genomic information with higher-order functional information using a system with KEGG Orthology (KO) numbers. Gene sequences in the GENES database, which includes collection of gene catalogs from completely sequenced genomes, are assigned KO numbers which are used to extract higher-order information, such as functional categories and pathways, from the KO database [120].

Similarity searches using BLASTx, version 2.2.13 [114], with the default settings were performed for the 1335 unique S. salmonicida protein-coding genes against the GENES database, using the E value cutoff e^-5. 249 returned a best match to which a KO -number had been assigned, and an additional 301 had non-best hits with assigned KO numbers. Functional annotations were extracted from the KO database for these 550 genes. However, 27 genes were annotated in the functional category "Human Disease"; five of these could be manually assigned to functional categories, while 22 were classified as conserved hypothetical proteins. Thus, 528 genes could be assigned a function using the KEGG databases. Similarity searches against all databases were performed for the remaining 785 genes, using the E value cutoff e^-20. 219 of these had matches, 63 of which could be used for functional annotation. No meaningful functional information could be found for remaining 156; these were annotated as conserved hypothetical proteins. In addition, the 553 genes which had best hits with E values between e^-20 and e^-5 were annotated as conserved hypothetical proteins, and the remaining 13 genes that did not have any match better than e^-5 were annotated as hypothetical proteins. To summarize: 591 protein-coding genes could be assigned a putative function, 731 were conserved hypothetical proteins, and 13 hypothetical proteins. The annotations for all proteins mentioned in the text have been refined manually.

The fraction of cysteine within the conceptually translated amino acid sequences was found to vary considerably between genes. While the majority of genes would encode a protein with cysteine content below 5%, 149 of the 1335 unique genes would encode a protein with more than 10% cysteine. These were annotated as cysteine-rich proteins. 122 are conserved hypothetical proteins, while 27 were assigned to functional categories.

A data set of 1044 yeast proteins localized in the mitochondrion [121] was used to search against the identified S. salmonicida proteins using BLASTp, version 2.2.13 [114], with the default settings. The matches with e values below e^-5 were searched against the total yeast proteome with the same cutoff, which identified 71 reciprocal best matches between the datasets. Based on the quality of the BLASTp hit (or protein length, to short to give any confidence in the identified HSP), phylogenetic trees and functional annotation of BLAST hits, only two entries were considered as strong candidate mitochondrial/mitosomal proteins and one as possible candidate proteins (see main text).

Phylogenomic and phylogenetic analyses

One goal of this project was to identify genes in the S. salmonicida genome that have been affected by LGT. We used phylogenetic analyses to identify gene transfers, as similarity searches are only poor indicators of such events [122]. We used the PhyloGenie package [83] to perform efficient phylogenetic analyses for a large number of genes. It takes amino acid sequences in standard (fasta) format, and performs similarity searches (BLAST) against a protein sequence database consisting of publicly available databases and/or sequence data released from ongoing genome projects. The program uses HMMER [123] to build a Hidden Markov Model (HMM) profile from the alignment in the BLAST result file, and uses it to search the full-length BLAST hits; sequences scoring better than a selected threshold are included in the dataset. Finally, an alignment of these sequences is made using the profile. No attempt was made to exclude putative paralogs from the datasets.

Cysteine-rich proteins were excluded from the phylogenomic analysis, since they were found to be problematic to align probably due to their frequent possession of highly repetitive sequences. PhyloGenie was run on the 1174 annotated non-cysteine-rich genes with homologs in the databases, using default settings except that coverage of the query sequence compared to the database sequence was not used as a selection criteria (coverage = -1) since coverage is meaningless with partial gene sequences such as EST and GSSs, and the maximum number of sequences in the dataset was set, for practical reasons, to 200 (seqs = 200). Where more than 200 sequences pass the E-value cutoff, the 200 with the lowest E-values in the HMM-based search will be retained. A custom-implemented function (T. Frickey, personal communication) was used that allows for exclusion of similar sequences from the dataset based on the taxonomic level (i.e. only one sequence is retained as a representative of all Bacillus sequences above a set identity). To reduce the size of the dataset as much as possible without losing too much information, an identity cutoff of 80% was used regardless of the taxonomic description in the sequence (taxlevel = 0, maxsim = 0.8).

These settings yielded 932 datasets with four sequences or more, among which 711 included 100 or more aligned amino acid positions assessed on the S. salmonicida input sequence. Only these were analyzed further, to avoid misleading results based on too-short alignments. A Perl script provided in the PhyloGenie package [83] was modified to use the fast maximum likelihood program PHYML version 2.4.4 [124]. Trees were obtained for each of the 711 datasets using this script and the default PHYML settings, except that the Whelan and Goldman (WAG) substitution model [124] was used together with a mixed four-category discrete-gamma model of among-site rate variation plus invariable sites (WAG + Γ + Inv).

The trees were automatically rooted using the approach implemented in Phylogenie; taxonomic information for all sequences in each tree is used to put the root at the most-basal node that is most-distant from the sequence that was used to select the dataset [see Figure 4 in [83] for details]. Phatg, the tree-browsing program within the PhyloGenie package [83], was used to identify 136 putative cases of LGT among the optimal trees, using the criterion that a diplomonad sequence was required to group with homologs from up to four different non-diplomonad eukaryotes within a prokaryotic clade. Manual inspection of the positive trees excluded from further analysis 31 trees for which the input S. salmonicida sequence failed to show convincing indications of LGT, or the G. lamblia homolog branched with a prokaryotic sequence. Bootstrap analyses were performed for the remaining 105 datasets, using 100 replicates and the methods and settings described above. Visual inspection of the consensus bootstrap trees yielded 84 that showed indications of LGT affecting S. salmonicida. Among these, two appear to be clones derived from different parts of the same S. salmonicida gene, and 14 are part of diplomonad gene clusters which may have formed via gene duplication events after a putative transfer event. Thus, our phylogenomic analysis identified 68 putative LGT events (Additional files 3, 4, 5, 6).

Nucleotide sequence accession numbers

The sequences reported here were deposited in GenBank at the National Center for Biotechnology Information [Genbank:EC585128-EC586011, GenBank:DX913336-DX918405, GenBank:DQ812518-DQ812528].