Estimating the extent of horizontal gene transfer in metagenomic sequences

Horizontal gene transfer (HGT) is believed to be a very important phenomenon in prokaryotic evolution, as it enables the acquisition of new genes or sets of genes that can accelerate evolution and adaptation to new environments or changing conditions.

In recent years, the availability of complete genome sequences from a wide range of prokaryotic species has made it possible to find many instances of potentially transferred genes in different genomes [1]. Nevertheless, the actual extent and importance of these events, and therefore the relevance of HGT in the evolution of prokaryotes, is still a matter of debate [1–6].

HGT has been extensively studied for sequenced genomes, but its real extent in natural populations in different environments remains unknown. Many of these communities are interesting from the ecological, biotechnological and/or clinical points of view, and therefore it is very important to understand how their members interact, evolve and innovate. HGT can be a very relevant factor since it may influence the dynamics of a community, and a recent report emphasises its important contribution to niche specialisation and adaptation [7]. However, the real impact of HGT in natural populations is unknown, and to date we do not know whether these events preferentially influence communities in particular environments. The advent of metagenomic sequencing [8, 9] provides enough environmental DNA sequences to address the study of global patterns of prokaryotic evolution in natural communities, including HGT.

Two types of methods have been used to detect HGT events in genomic sequences: Phylogenetic methods, based on the examination of the phylogenies of individual genes or proteins; and compositional methods, based on the analysis of DNA composition, which is assumed to contain some evolutionary information in the form of species-specific signatures.

In theory, phylogenetic analysis is a powerful tool used to recognize HGT events. On building the phylogenetic tree of a family of sequences from different species, putative HGT can be identified when the position of a specie in the tree does not match that of a reference phylogeny, which can be obtained in different ways: from universally conserved genes such as 16S rRNA or some tRNA synthetases, or by supertrees of conserved genes merged together [10]. Phylogenetic methods are therefore usable with individual genes, as long as a set of orthologues is known. In the practice, phylogenies are often noisy and/or inconclusive for several reasons (lack of homologues for specific genes, phylogenetic artefacts such as long-branch attraction, inclusion of paralogues, uneven mutation rates, etc). Therefore, the results of a phylogenetic analysis must be carefully examined before proposing HGT. This task is more difficult for metagenomic sequences because the origin of the sequence is almost always unknown (since the metagenome is a mixture of sequences from different species) and, consequently, the construction of a reference tree is not straightforward. Nevertheless, at least one attempt has been made to study environmental sequences by means of the comparison of individual phylogenies [11], which must be taken as a valuable proof of concept.

Several authors have used these phylogenetic approaches to detect and quantify HGT in genomes, orthologous groups of genes and protein families [12–16]. These methods estimate a rather low extent of HGT in all cases, usually below 2%. A possible explanation for this is that, since the barriers to HGT seem much easier to overcome when closely related species are involved [17], it is likely that most HGT events are produced between species belonging to the same taxon. These events are much more difficult to detect using phylogenetic analysis as the resolution of the trees decreases for such a small range. Consequently, it is likely that specific HGT events cannot be conclusively determined using phylogenetic approaches, leading to an underestimation in the frequency of HGT.

Compositional methods are based on the assumption that individual genomes have characteristic features regarding their DNA composition, which has been termed DNA or genomic signatures. Each genome has its own pattern of GC content and codon usage due to a combination of environmental and genetic factors [18–20]. Recently transferred regions can be detected since they possess a different DNA signature to that of the host genome which reflect their foreign origin [21–23]. Therefore, HGT events in complete genomes could be detected by recognising zones of different or unusual compositional patterns. In these cases, the genomic signature will be closer to those of the donor genomes, which in some cases can reveal the phylogenetic origin. As time passes, the foreign sequence is shaped by the biases in the replication and reparation machinery, which causes its signature to progressively resemble that of the host. This process is known as amelioration, and implies that the original sequence can no longer be recognised as that from a foreign origin on the basis of nucleotide composition alone [24]. Therefore, these methods perform better to detect recent HGT events, where genomic signatures differ significantly between the donor and recipient genomes.

Many compositional methods based on several approaches and algorithms exist [21–23, 25–29]. The percentage of HGT estimated by compositional methods is significantly higher than by methods based on phylogenies. In addition, different methods have been shown to predict very different sets of genes as having been horizontally transferred [28, 30–32]. This apparently discouraging result has been recently interpreted to mean that each method is tuned to detect gene transfers of different ages [32]. Normally, compositional methods have been used to analyse large genomic sequences since a significant amount of sequence is needed to extract the global features. This is a serious drawback for its use with metagenomic sequences.

Both the phylogenetic and compositional approaches have their advantages and drawbacks, even though the two are complementary since they use independent sets of information. They make predictions based on different features and can, therefore, be combined to provide a more comprehensive insight into the relevance of HGT processes.

Here, we propose a double approach to measure the incidence of HGT in short sequences of unknown origin as these are obtained by metagenomic sequencing of natural prokaryotic communities. Both a compositional and phylogenetic method were derived and tested in order to apply them to complex mixtures of sequences. The methods were then used to study the sequences of several available metagenomic projects.

We have developed two different methods to detect possible HGT events, both of which aim to find transitions in the DNA sequences. The first is a phylogenetic method that relies on the taxonomic assignment of the ORFs found in the contig to identify the instances of ORFs assigned to disparate taxa (related approaches can be found in [39–41]). The second is a compositional method that scans DNA, and looks for zones of differential oligonucleotide usage. We now go on to describe both methods in detail.

Phylogenetic assignment method

For the homologues found for a given putative ORF, a distance score is calculated as the inverse of the bit score of their alignment. The taxonomy for each homologue is obtained from the GenBank database and chopped to the desired taxonomic rank (in decreasing order: kingdom (first rank), phylum (second), class (third), order (forth) and family (fifth)). When more than one homologue is found for the same species, only the best hit is kept, thus effectively excluding paralogues. The ORF is assigned to a given taxon if the lowest distance scores are those that correspond to the homologues from it. Only two homologues from other taxa are allowed to have a lower score. If this is not the case, we often find that the ORF has many close homologues from a single taxon, in addition to several others from different taxa. This could be the consequence, for instance, of a high degree of transfers in that particular protein family. In such cases, the ORF is assigned to the taxon whose average distance score is at least 15% lower than that of any other. Examples of this procedure are depicted in Additional file 4.

In any case, an ORF requires at least three homologues with the taxon to be assigned to it. To avoid assignments when only distant homologues are found, a minimum identity is required between the ORF and the homologues (a default value is 30% of identity), and a minimum distance score is also needed. All homologues below these levels are not considered.

If the ORF could not be assigned by this procedure, a list of candidate taxa is kept for it, containing the average distance scores for all taxa for which at least one homologue has been found. This list will be subsequently used to improve the assignments.

Usually, between 10–40% of the ORFs of real metagenomes can be accurately assigned in this way. It is possible to take advantage of the fact that the contigs often contain more than one ORF to increase the number of assignments while maintaining a very high accuracy. For instance, let us suppose that a contig contains three ORFs, and we can clearly assign a common taxon for two of them. If the third ORF could not be assigned, we could assume that it belongs to the same taxon, provided it has close homologues from it (Figure 1). Therefore, we have introduced an additional stage that takes into account the assignments for the neighbouring ORFs. The unassigned ORF is scored in the following way: the first five taxa in the list of candidate taxa are scored with 1/r, where r is their rank in that list. The scores of the neighbouring ORFs are then added to this. If the neighbouring ORF has already been assigned to a given taxon, a score of one is added only to that taxon. An ORF is assigned to the taxon with the highest score, given that this is higher than 1.5. For instance, a frequent scenario is that of an unassigned ORF having both its closest homologue and a neighbouring ORF belonging to a given taxon. In this case, the final score for that taxon will be equal to two, and the ORF will be assigned to it. An ORF cannot be assigned to a taxon that it is not in its candidate list, which implies that the ORF must have at least one homolog belonging to the taxon.

This procedure increases the number of assignments in 10–12% for random metagenomes, and precision remains unaltered. For real metagenomes, the increase in the number of assignments depends greatly on the average number of ORFs in the contigs (see Additional file 1), which ranges between 3% for soil and 20% for Sargasso sea metagenomes.

Once the ORFs have been assigned in this way, a possible HGT event will be predicted when a contig contains sequences from different taxa (Figure 2A). Unassigned ORFs are not considered.

In addition, a contig is assigned to a particular taxon if 75% of its ORFs, or more, belong to it, which can be used to study the taxonomic composition of metagenomes.

The task of detecting HGT in metagenomic sequences is difficult. The available methods are designed to work with complete genomes (or at least with a considerable amount of sequence), and it is assumed that one knows the species from which the sequences are taken: phylogenetic methods attempt to identify deviations with regard to the expected phylogeny, while compositional methods detect differences in the compositional profile of each ORF as compared to that of the complete genome. These methods are not applicable to metagenomic sequences. Metagenomes are usually composed of a complex mixture of short sequences belonging to many different species.

Not knowing the source species for each sequence in the metagenome is the most serious hurdle for all metagenomic studies. One of the most relevant advances in metagenome analysis will be to solve the problem of binning, that is, the assignment of metagenomic sequences to source species or taxa. Although some attempts have been made by compositional methods to address this task [47–49], the problem remains unresolved. As it stands, more complete genomes would be required for training, as would longer sequences than those currently supplied by metagenomic projects [50]. Nevertheless, it is likely that the length of the sequences obtained in these projects will not change very much in the near future since metagenomics is tightly linked to the use of fast and cheap pyrosequencing machinery [51], which currently obtains short readings and contigs (although considerably efforts are being made to improve read lengths). Therefore, either more powerful methods or different approaches will be needed to address the task of assigning sequences to species or taxa. We have presented herein the usage of a highly precise method, based on sequence homology, which can provide an interesting alternative for complementing compositional approaches. The advantages of methods based on sequence homology and phylogenies are that they are unaffected by the short length of the contigs, they can work with partial sequences, and can even classify relatively unknown taxa with high accuracy. It is also possible to take advantage of the presence of several ORFs in the contig to support the assignments, thus increasing precision and, especially, sensitivity.

It is very likely that both methods presented herein significantly underestimate the number of HGT events. Compositional methods can only detect recent transferences before amelioration erases the distinctive signatures of the sequences. When dealing with metagenomic sequences, the situation is even more challenging since the analysis must only be done for a very short segment of the sequence, and we lack the contextual information provided by the rest of the genome. On the other hand, phylogenetic methods usually rely on distinguishing incongruent phylogenies, identifying the branches of a tree that are in atypical positions with respect to a reference. Once again, this is a difficulty for metagenomic sequences because the species source for the sequence is unknown. Therefore, it is not possible to assess its position in the tree. We rely on assigning the most probable taxon for each ORF to distinguish instances of neighbouring sequences of different likely origins. The main disadvantage of this procedure is that it has a limited sensitivity. It is very difficult and risky to perform assignments at deep taxonomic ranks, and it is probably impossible to do so at the species level. Therefore, we will only detect events involving distant taxa, and will overlook those that happen between closely related species or taxa, which are probably the most abundant.

The results presented in this work must be taken as an indication of possible trends for HGT events in prokaryotic samples from diverse environments. However, it must be taken into account that many different factors, some yet unknown, can influence the prevalence of HGT in a community. For instance, we do not know how the composition of the community shapes the likelihood of HGT events. HGT could be expected to be widespread in a very diverse community in which many different genes and subsystems are present, and transferences could greatly accelerate evolution and adaptation to environmental challenges [52]. On the other hand, HGT might be easier in a community composed of closely related species for which the barriers to the transference can be more easily surmounted [17]. Furthermore, some species seem to be more prone to HGT events than others [53–55]. The prevalence of HGT events is certainly influenced by both environmental characteristics and the composition of the community, and as yet, we do not know the balance between both factors. Furthermore, we must realize that we are observing putatively transferred genes in the genomes of the species present in these communities. Obviously, these species are present in other communities and environments, and the putative HGT events may have occurred in any of these. For instance, we have found transferences between different archaea in the human gut metagenome (Additional file 2). The only archaea that has been found in measurable amounts in the human gut is Methanobrevibacter smithii. This indicates that these putative transfers could have happened in some other environment in which M. smithii can co-exist with other archaea, for example, sewage waters. A similar circumstance has been found in another study of bacteria from the human intestine [7].

This study attempts to provide the first methodological step for the study of the dynamics of HGT in prokaryotic communities. Further studies should deal with the influence that community composition and environmental factors have on HGT rates by addressing questions such as: Is HGT more widespread in more diverse communities? Or is it favoured in cases where the community is dominated by a few closely related taxa? To what extent do environmental challenges affect the rate of gene transfer? Is there a correlation between the amount of transferences and the presence of bacteriophages?

The answer to these questions will shed light on the true role of HGT in adaptation and innovation, and will greatly improve our understanding of the dynamics of natural prokaryotic communities.