Genome-Wide Molecular Clock and Horizontal Gene Transfer in Bacterial Evolution

We describe a simple theoretical framework for identifying orthologoussets of genes that deviate from a clock-like model of evolution.The approach used is based on comparing the evolutionary distanceswithin a set of orthologs to a standard intergenomic distance,which was defined as the median of the distribution of the distancesbetween all one-to-one orthologs. Under the clock-like model,the points on a plot of intergenic distances versus intergenomicdistances are expected to fit a straight line. A statisticaltechnique to identify significant deviations from the clock-likebehavior is described. For several hundred analyzed orthologoussets representing three well-defined bacterial lineages, the ${alpha}$ -Proteobacteria, the ${gamma}$ -Proteobacteria, and the Bacillus-Clostridiumgroup, the clock-like null hypothesis could not be rejectedfor $~$ 70% of the sets, whereas the rest showed substantial anomalies.Subsequent detailed phylogenetic analysis of the genes withthe strongest deviations indicated that over one-half of thesegenes probably underwent a distinct form of horizontal genetransfer, xenologous gene displacement, in which a gene is displacedby an ortholog from a different lineage. The remaining deviationsfrom the clock-like model could be explained by lineage-specificacceleration of evolution. The results indicate that althoughxenologous gene displacement is a major force in bacterial evolution,a significant majority of orthologous gene sets in three majorbacterial lineages evolved in accordance with the clock-likemodel. The approach described here allows rapid detection ofdeviations from this mode of evolution on the genome scale.

INTRODUCTION

Top

Introduction

The classical molecular clock concept holds that the sequenceof a given gene (protein) evolves at a constant rate as longas its biological function remains unchanged (5, 26, 56, 57).The actual evolutionary rates for the entire set of genes ina genome show a broad distribution, with the rates for the slowestgenes differing from the rates for the fastest genes by at least2 orders of magnitude (18, 21, 26, 28). However, the molecularclock hypothesis holds that, within the same set of orthologs(i.e., homologous genes related via vertical inheritance) (16,44), the rate does not change, within some inevitable dispersion.The theoretical foundation of the molecular clock concept isthe neutral theory of molecular evolution, which holds that,at least when no functional changes occur, the great majorityof fixed substitutions in nucleotide and protein sequences areneutral, their rate being determined by the fixed, gene-specificfunctional constraints (13, 26, 39). The modern reincarnationof neutralism, the near-neutral theory, predicts a greater dispersionof the molecular clock because the probability of fixation ofslightly deleterious mutations critically depends on the effectivepopulation size, which is prone to major fluctuations (13, 39).Numerous tests of the molecular clock revealed both a generalpattern conforming to the theory and numerous violations. Ithas been repeatedly shown that the molecular clock is overdispersed;i.e., the variation of evolutionary rates within a set of orthologsis greater than that predicted by a null model based on thePoisson distribution (8, 18, 19).

Deviations from the molecular clock are thought to result fromlineage-specific acceleration of evolution, which could be dueeither to functional changes entailing relaxation of purifyingselection or positive selection or to increased mutational pressurecaused, at least in part, by effective population size effects(5). These phenomena cause overdispersion of the molecular clock,which is manifested in unequal lengths of tree branches comingout of the same node, under the assumption that the topologyof the phylogenetic tree for a given set of orthologs is known.Typically, the same species tree topology is assumed for allgenes. This approach is likely to be valid for the multicellulareukaryotes, which, historically, have been the objects of theanalyses that led to the molecular clock concept. However, recentcomparative genomic studies strongly suggest that in additionto the regular pattern of vertical inheritance, evolution ofprokaryotic genomes is dramatically affected by horizontal genetransfer (HGT) (6, 11, 12, 27, 30, 32, 38, 55). On many occasions,HGT seems to occur between evolutionarily distant organisms,although it has been argued that there could be a decreasinggradient of the HGT rate from closely related species to distantlyrelated species; the existence of such a gradient could be oneof the reasons why a species tree can be constructed at all,in spite of extensive HGT (20).

There seem to be certain connections between the amount of HGTand the biology of prokaryotic genes. In particular, the so-calledcomplexity hypothesis holds that HGT is much less common amonggenes that encode subunits of macromolecular complexes, suchas those involved in translation, transcription, and replication,than in genes coding for metabolic enzymes (25). While thisprediction might hold statistically, subsequent studies haveshown that there are very few, if any, genes that are completelyrefractory to HGT. In particular, evidence of HGT has been obtainedfor several ribosomal proteins, translation factors, and themajor RNA polymerase subunits (3, 4, 24, 34).

From a comparative genomic perspective, HGT events have beenclassified into three categories: (i) acquisition of genes thatare novel to a given phylogenetic lineage; (ii) acquisitionof paralogs of genes preexisting in the given lineage; and (iii)xenologous gene displacement (XGD), in which the original genefrom a given set of orthologs is displaced by a member of thesame set of orthologs from a different lineage (30).

Obviously, if an HGT event, particularly XGD, goes unnoticedin the course of phylogenetic analysis, an apparent gross violationof this molecular clock will be seen when evolutionary ratesare measured in the affected set of orthologous genes on thebasis of the assumed species tree topology. HGT is detectedthrough anomalies in the topology of phylogenetic trees of individualsets of orthologs or by so-called surrogate approaches, which,in the case of HGT between distant species, are based primarilyon the phyletic distribution of the homologs of a given gene(7, 30, 41). Simply put, unexpected phyletic patterns (e.g.,the presence of orthologs of a given gene in all or nearly allsequenced bacterial genomes but in only one archaeon) suggestthat there has been HGT (in this case from a bacterium to thearchaeon) (27, 30). These patterns can be expressed either interms of presence-absence only or, more quantitatively, by comparingthe significance levels of taxon-specific best hits. The generalvalidity of this approach seems to be supported by the biologicalplausibility of some of the trends in the apparent horizontalgene fluxes that were detected by phyletic pattern analysis.Thus, hyperthemophilic bacteria showed a clear preponderanceof genes of possible archael origin compared to mesophiles (2,29, 36), and probable gene transfer from eukaryotic hosts tosome bacterial pathogens also has been inferred (17, 40). Inprinciple, phylogenetic tree analysis is supposed to be a moreprecise indicator of probable HGT events than similarity-basedsurrogate methods because of inaccuracies in the latter resultingfrom the lack of exact correspondence between sequence similarityand phylogenetic affinity (31). A genome-wide phylogenetic analysis,aimed specifically at detection of horizontally transferredgenes, has been described (43). However, it is well known thatphylogenetic analysis is fraught with its own slew of artifacts,such as long branch attraction, particularly when fast methods,such as minimal evolution or neighbor joining, are employed(14). In addition, phylogenetic analysis can be prohibitivelyexpensive computationally when it is attempted on the genomescale, especially with powerful methods, such as complete maximum-likelihoodanalysis, and large sets of species. Therefore, surrogate, similarity-basedmethods have proved to be extremely useful, at least as a rapid,first-tier strategy that allows workers to delineate a set ofHGT candidates.

We were interested in investigating a surrogate approach togenome-wide study of prokaryotic evolution, which combines atest of the validity and an analysis of the rate distributionof the molecular clock with detection of potential HGT eventsand lineage-specific acceleration of evolution. Using the Clustersof Orthologous Groups (COGs) database for proteins (46, 47),we analyzed the molecular clock behavior of COGs from threemajor bacterial lineages, the ${alpha}$ -Proteobacteria, the ${gamma}$ -Proteobacteria,and the low-G+C-content gram-positive bacteria. We found thatclock-like evolution was dominant in all three groups, but wealso detected many anomalies, some of which are best explainedby XGD.

MATERIALS AND METHODS

Top

Materials and Methods

Sequences. Three unequivocally identified and well-characterized bacteriallineages, the ${gamma}$ -Proteobacteria, the ${alpha}$ -Proteobacteria, and theBacillus-Clostridium group of gram-positive bacteria, were selectedfor the present study. The ${gamma}$ -proteobacterial set included sixspecies: Escherichia coli K-12, Haemophilus influenzae, Pasteurellamultocida, Salmonella enterica serovar Typhimurium LT2, Vibriocholerae, and Yersinia pestis. The ${alpha}$ -proteobacterial set includedseven species: Agrobacterium tumefaciens C58 Cereon, Brucellamelitensis, Caulobacter crescentus CB15, Mesorhizobium loti,Rickettsia conorii, Rickettsia prowazekii, and Sinorhizobiummeliloti. The Bacillus-Clostridium set included eight species:Bacillus halodurans, Bacillus subtilis, Clostridium acetobutylicum,Listeria innocua, Lactococcus lactis, Staphylococcus aureusN315, Streptococcus pneumoniae TIGR4, and Streptococcus pyogenesM1 GAS. For each of these species sets, we identified a setof COGs (46, 48) in which each of the relevant species was representedby exactly one protein (i.e., all species present, no paralogs).Additionally, the constituent proteins were required to havea sufficient alignable length (at least 60 amino acids in conservedblocks [see below]). This search resulted in 563 COGs for the ${gamma}$ -proteobacterial set, 274 COGs for the ${alpha}$ -proteobacterial set,and 234 COGs for the Bacillus-Clostridium set; the overlap forthe three sets comprised 114 COGs.

Alignments. Multiple alignments of sequence families within the bacterialgroups were produced by using the MAP program (23). Sequencefamilies involving wider taxonomic sampling of proteins werealigned by using the T-Coffee program (37). Multiple alignmentswere aggressively filtered for potential incorrectly alignedpositions; only conserved blocks with no gaps containing 10or more positions were retained for further analysis (53).

Evolutionary distances between genes and genomes and phylogenetic trees. Maximum-likelihood distances between individual protein sequenceswere computed for each of the COGs analyzed by using the PAMLpackage, with the JTT substitution model corrected for observedamino acid frequencies and the ${alpha}$ parameter of ${gamma}$ -distribution ofintraprotein evolutionary rate variability set to 1.0; thisestimate includes a correction for possible multiple substitutionsin the same site (54). Additionally, multiple alignments consistingof 21 sequences each were produced for the 114 COGs in whichall three groups were represented, and pairwise distances weresimilarly computed from this alignment. Phylogenetic trees forindividual COGs were constructed by using the maximum-likelihoodmethod implemented in the Tree-Puzzle package, with the expectedlikelihood weight determined for the tree topology involvingsplit versus topologies in which the respective lineage remainedmonophyletic (42).

To calculate intergenomic evolutionary distances, all intergenicdistances obtained from the same pair of genomes in the setof 114 conserved COGs were pooled, and the median of the distributionof these distances was taken to represent the intergenomic distance(21, 52). Neighbor-joining and least-squares trees were reconstructedfrom the pairwise genome distance matrices by using the programsNEIGHBOR and FITCH of the PHYLIP package, respectively (15).

Limited tree analysis. A comprehensive phylogenetic analysis of a protein family, evenif it is limited to a set of orthologs from completely sequencedgenomes, is rarely feasible. There are several reasons thatusually preclude this type of analysis, as follows: in manycases, the sequences of orthologs from the most distant lineages,such as bacteria and archaea, are only weakly similar, whichcomplicates the construction of an alignment suitable for buildinga phylogenetic tree; the large number of sequences hampers theuse of advanced methods for reconstruction of the phylogeny;and the presence of in-paralogs of various ages hinders theinterpretation of results. Therefore, we elected to performa limited tree analysis for the cases where the split distanceanalysis indicated a likely split of a set of orthologs froma particular bacterial lineage into two subsets with distinctevolutionary histories; these groups are referred to below asleft ({L}) and right ({R}) sets. For each sequence from thegiven COG, global pairwise alignment scores (calculated by usingthe ALIGN program with default parameters) (35) for alignmentswith the other COG members were obtained (sequences from closelyrelated species were excluded). Members of the COG not belongingto either {L} or {R} (i.e., the sequences of orthologs outsidethe lineage analyzed) were arranged into two ordered lists,<H^L> and <H^R>, according to their similarity scoresagainst {L} and {R}, respectively. Multiple alignments thatincluded {L}, {R}, and the top five sequences from <H^L>and <H^R> were constructed by using the T-Coffee program(37). Maximum-likelihood trees were constructed by using theProtML program of the MOLPHY package (1). For all internal branchesseparating the {L} and {R} subsets (see Fig. 4), the RELL bootstrapvalues were determined. The highest bootstrap value observedon such a branch indicated the level of support for the separationhypothesis.

View larger version (9K):
[in this window]
[in a new window]
FIG. 4. Limited tree procedure. Formally, let us consider an internal branch B, which partitions all nodes into four subsets, {S₁ } to {S₄}. If the partition satisfies the criteria (i) {S₁} contains all of {L} and none of {R}, (ii) {S₂} contains all of {R} and none of {L}, and (iii) {S₃} and {S₄} contain neither {L} nor {R}, the branch provides the evolutionary separation between {L} and {R} regardless of the true position of the tree root.

RESULTS AND DISCUSSION

Top

Results and Discussion

Theory. (i) Molecular clock and deviations from clock-like evolution. The approach developed here is based on comparing the evolutionarydistances between genes (proteins) in a given COG to the distancesbetween the corresponding genomes. There is no single correctway to calculate the intergenomic distances. In principle, areference gene believed to be a good clock, such as rRNA, ora set of genes believed to evolve under the same model, suchas genes for ribosomal proteins (3, 22, 51), can be employedfor this calculation. We chose one of the genome-wide approachesto evolutionary distance determination, in which the medianof the distribution of the distances between orthologs froma given pair of genomes is used as a proxy for the intergenomicdistance (51, 52) (see Materials and Methods).

Under a perfect molecular clock and strict vertical inheritance,the evolutionary rates of all genes differ from each other onlyby proportionality constants, so the distance between any pairof genes (d_AB) is proportional to the distance between genomesA and B (D_AB):

(1)
where v is the relativerate of evolution for the given gene (COG) (v = 1 for COGs thatevolve at the median rate). Obviously, in this case, a scatterplot of d_AB versus D_AB for all pairs of genomes is a straightline with the intercept at 0 and the slope equal to v (Fig.1A).

View larger version (16K):
[in this window]
[in a new window]
FIG. 1. Clock-like model of evolution: predictions and detection of deviations. (A) Clock-like evolution with vertical inheritance. (B) Clock-like evolution with one HGT from an outside source. The red branches in the trees indicate the genes from species B and C that are inferred to have been transferred from the unknown, distant species, X, and the red points on the plot at the bottom right indicate the evolutionary distances that are taken as evidence of this HGT event.

Various events that could occur during evolution could distortthis idealized picture. Accelerated evolution on one terminalbranch of the tree would increase all distances involving thecorresponding species; as a result, the points correspondingto this species on the scatter plot are located above the mainline. Similar acceleration on an internal branch of the treewould result in an upward shift of the points correspondingto two or more species. We were particularly interested in casesof HGT of a gene encoding a protein of a given COG from an outsidesource, resulting in displacement of the original gene characteristicof the given lineage (i.e., XGD). In these cases, the evolutionarydistance between the transferred and nontransferred sequencesis determined by the distance between the transfer source andthe recipient lineage (D*), rather than by the intragroup distances.The scatter plot is expected to consist of two distinct setsof points: (i) points corresponding to nontransferred genesthat fit a line with the intercept at 0 and the slope v, likethose obtained with the clock model, and (ii) points correspondingto the transferred gene (i.e., the distances between this geneand all other members of the COG; these points fit a horizontalline with the intercept corresponding to the distance to thetransfer source) (Fig. 1B).

(ii) Statistical model. The accumulation of substitutions in protein sequences was treatedas a Poisson stochastic process. The number of substitutionsin a given gene, accumulated over a given time interval, isa Poisson-distributed variable with the variance proportionalto the expected value:

(2a)
If thefact that the observed distance between proteins encoded bygenomes A and B has a random error, E_AB (which includes Poisson-distributedsampling error, error in deriving evolutionary distances fromobserved differences, and branch-to-branch variation of evolutionaryrates), is taken into account, equation 1 becomes:

(2b)
By using equation 2a, equation 2b can be rewrittenas

(3a)
where e_AB is a random variablewith the expectation of 0 and the variance independent of D_ABand v. Similarly, if one of the two orthologs was acquired viaHGT from a distant source (Fig. 1B):

(3b)

(iii) Statistical analysis: vertical inheritance. Let us consider a set of N genomes {G}, each with a single orthologin a given COG. One can measure the distance between orthologsin the given COG (d_IJ) and the distance between the genomes(D_IJ) for all N' = N(N – 1)/2 pairs ([I,J]) from {G}.Minimizing the square error over all such pairs,

weget the optimal value of v:

(4a)
andthe fit error (s²) and residual variance (u²) of e_IJ:

(4b)

(iv) Statistical analysis: HGT. If genes in one of the lineages in a COG were acquired via HGTfrom a distant source, all intergenic distances in this COGfall into two groups: those that reflect vertical inheritance(D_AD and D_BC in Fig. 1B) and those that reflect the transferdistance (D_AB, D_AC, D_BD, and D_CD in Fig. 1B). Considering allpairs that correspond to vertical inheritance ([I,J]) and allpairs that correspond to HGT ([K,L]), the square error of theapproximation is:

Minimization of E²over v and D* yields:

(5a)

(5b)
and

(5c)

(v) Statistical analysis: baseline noise. If the pattern of a gene's inheritance is completely disjointedfrom the pattern of intergenomic relationships, neither equation3a nor equation 3b adequately describes the relationships betweenthe intergenomic and intergenic distances. In the absence ofa clear dependence between these variables, the scatter plotfor d_AB versus D_AB represents random scatter of points, andthe following simple equation applies:

(6a)
where is simply the mean intergenic distanceover all pairs. The baseline variance of e_IJ is:

(6b)

(vi) Statistical analysis of COG evolution. Each COG was analyzed with the three models described above.

(a) Noisy data, no clock-like evolution (equation 6a). The baseline variance of intergenic evolutionary distances (u_N²)was calculated by using equation 6b.

(b) Simple molecular clock (equation 3a). The residual variance of the straight line fit (u_C²) and thefit error (s_C²) were calculated by using equation 4b. The relativeevolutionary rate (v) was calculated by using equation 4a.

(c) Single significant deviation from molecular clock (equations 3a and 3b). The genomes were partitioned into two sets by breaking eachbranch of the species tree. For each of the possible splits,the residual variance (u_T²) and the fit error (s_T²) were calculatedby using equation 5c. The split with the minimal s_T² was accepted.The relative evolutionary rate (v) was calculated by using equation5a, and the distance to the transfer source was calculated byusing equation 5b. Additionally, to detect the transfers originatingfrom outside the group, the relative transfer distance (D_T)was calculated as D*/max(D_KL), with the maximum taken over allcross-group pairs.

Two statistical tests were performed for each COG. The firsttest aimed at discriminating between H₀ (the data do not followeither of the two clock-like models) and H₁ (the data fit eitherthe simple-clock or single-transfer model). The ratio F_C = u²_N/min(u²_C,u²_T) was subjected to Fisher's test with (N' – 1, N' –3) degrees of freedom. If the value of F_C exceeded the criticallevel at the 0.05 level of significance (1.94 to 2.64, dependingon the group of species analyzed), H₀ was rejected, and thedata were considered to conform to the molecular clock model.

The second test discriminated between H₀ (the data fit the simple-clockmodel) and H₁ (the data fit the single-transfer model). Theratio F_T = s_C²/s_T² was subjected to Fisher's test with (N' –2, N' – 3) degrees of freedom. If the value of F_T exceededthe critical level at the 0.05 level of significance (1.95 to2.66 depending on the group of species analyzed), H₀ was rejected,and the data were considered to conform significantly betterto the single-transfer model. In this case, the value of v calculatedby using equation 5a (rather than that calculated by using equation4a) was used to describe the relative evolutionary rate forthe COG in question; a D_T value of >1 was considered to bean indication of HGT from an outside source.

Empirical results. (i) Molecular clock and deviations from clock-like evolution in bacteria. We applied the theory described above to an analysis of theevolution of three major bacterial lineages, the ${alpha}$ -Proteobacteria,the ${gamma}$ -Proteobacteria, and low-G+C-content gram-positive bacteria.For each of these groups, the COGs that contained a single representativefrom each species were selected for analysis (Table 1). Examplesof relationships between the intergenomic and intergenic distancesare shown in Fig. 2. Figure 2B shows clear evidence of the evolutionaryheterogeneity of the COG in question. Table 1 shows a breakdownof the COGs analyzed according to their fit to one of the threemodels of evolution described above. For the great majorityof the COGs (90 to 99%), the data could be reconciled with eitherthe clock-like model or the model that involved one significantdeviation from the clock (Table 1). For $~$ 70% of the COGs in eachlineage, the clock-like model was found to be compatible withthe data, whereas the evolution of the rest of the COGs wasbetter explained when the group was split into two subsets withdifferent histories. For 13 to 22% of the COGs in the threelineages, the apparent acceleration of evolution in one of theclades (D_T > 1) was significant enough to suspect HGT froman outside source.

View this table:
[in this window]
[in a new window]
TABLE 1. Alternative evolutionary models for COGs

View larger version (14K):
[in this window]
[in a new window]
FIG. 2. Three representative examples of the relationships between intergenic and intergenomic distances. (A) Clock-like evolution. (B) Clock-like evolution with one strong deviation probably due to HGT from an outside source. (C) Random scatter of points: uncertain evolutionary scenario. The magenta lines in panels B and C show the horizontal trend lines with the best fit to the data; the evolutionary distances thought to reflect XGD in panel B are indicated by magenta diamonds.

(ii) Relative evolutionary rates in the three bacterial lineages. All three bacterial lineages analyzed displayed a wide rangeof relative evolutionary rates (i.e., evolutionary rates normalizedby using the intergenomic rate). The difference between theCOGs that evolved fastest and the COGs that evolved slowestreached almost 2 orders of magnitude for the ${gamma}$ -Proteobacteriaand almost 1 order of magnitude within the set of 108 conservedCOGs that fit one of the two evolutionary models (Table 2).The distributions of relative evolution rates were close tolog normal in all three lineages (Fig. 3). Because the medianintergenomic distance was calculated only for the 114 COGs thatare conserved in all three bacterial groups (see Materials andMethods for details), the medians of the rate distributionswithin each of the groups were not equal to one. Of the threegroups, the ${gamma}$ -Proteobacteria showed the fastest median rate andthe broadest variation in the relative evolutionary rates (Fig.3A). This is probably due to the fact that the selected setof ${gamma}$ -Proteobacteria species forms a tighter (more recently diverged)cluster than the sets of ${alpha}$ -Proteobacteria and Bacillus-Clostridiumgroup bacteria analyzed, as indicated by the comparison of thecalculated intergenomic distances (data not shown). Accordingly,the set of ${gamma}$ -proteobacterial COGs is more functionally diverseand includes faster-evolving genes than the COGs from the othertwo lineages. This conclusion is supported by the comparisonof the rate distributions for the full set of 555 ${gamma}$ -proteobacterialCOGs and the 108-COG subset that is conserved in all three lineages;not unexpectedly, the distribution in the latter set is considerablynarrower and is shifted toward lower rates (Fig. 3B). Withinthe conserved set of 108 conserved COGs, the relative evolutionaryrates are strongly correlated among the three lineages (Table3). However, statistically significant differences between thelineages were detected; in the same COG, the rate among ${gamma}$ -Proteobacteriatended to be the lowest, and the rate among the ${alpha}$ -Proteobacteriatended to be the highest (Table 3).

View this table:
[in this window]
[in a new window]
TABLE 2. Relative evolution rates

View larger version (26K):
[in this window]
[in a new window]
FIG. 3. Distributions of the relative evolutionary rates in the three bacterial lineages analyzed and in the conserved set of COGs represented in all lineages. (A) Distributions of the relative evolutionary rates in the three lineages. gamma, ${gamma}$ -Proteobacteria; alpha, ${alpha}$ -Proteobacteria; bacil, Bacillus-Clostridium group. (B) Distributions of the relative evolutionary rates in the full set of ${gamma}$ -Proteobacteria (gamma-555) and in the subset of COGs represented in all lineages (gamma-108). The distributions for the three bacterial lineages and the corresponding log-normal approximations (dashed lines) are color coded. The scale for the horizontal axis is logarithmic.

View this table:
[in this window]
[in a new window]
TABLE 3. Evolution of the same COG in different groups

(iii) Apparent violations of molecular clock and anomalies in gene phylogenies. All COGs that fit the two-mode model better than the simpleclock-like model were analyzed further by using the limitedtree procedure (see Materials and Methods for details).

Briefly, the purpose of this procedure was to investigate whetherthe two sets of species separated by the statistical analysisdescribed above were also separated by a strongly supportedinternal branch in the phylogenetic tree for the given COG (Fig.4). Not surprisingly, the COGs that had the longest split distance(D_T) also showed a consistent tendency to have anomalies inthe reconstructed phylogenies (Table 4).

View this table:
[in this window]
[in a new window]
TABLE 4. COGs with most pronounced deviations from the clock-like model

Table 4 lists the top 10 COGs with the most pronounced deviationfrom the clock-like model (i.e., the greatest D_T value) foreach of the bacterial lineages analyzed. These COGS were furthersubjected to complete phylogenetic analysis, including statisticaltests (see Materials and Methods), in order to determine thenature of the anomalies detected. Upon detailed inspection ofthe trees produced with the three methods, the cases were classifiedinto those that most likely reflected HGT resulting in XGD,those that were best explained by accelerated evolution in aparticular lineage, and those for which the evolutionary scenarioremained uncertain. Of the 30 cases analyzed, only 5 did notappear to be resolved well enough to reach conclusions on theevolutionary mechanism; one case involved acceleration of evolution,and the remaining cases seemed to be best explained by invokingHGT, mostly resulting in XGD. The statistical tests comparingthe likelihoods of tree topologies with different positionsof the corresponding branches validated the anomalous clusteringand thus supported the case for HGT (Table 4). In functionalterms, the majority of the inferred cases of HGT involved metabolicenzymes, in agreement with the complexity hypothesis, whichpostulated that genes coding for proteins that are not involvedin tight interactions with other proteins as subunits of macromolecularcomplexes are more readily subject to HGT (25). It is also noteworthythat 5 of the top 30 cases included amianoacyl-tRNA synthetases,a category of enzymes which appears to be particularly proneto HGT (49, 50).

In COG0060 (isoleucyl-tRNA synthetase, IleS), both Rickettsiaspecies and C. acetobutylicum are separated from their sistergroups ( ${alpha}$ -Proteobacteria and gram-positive bacteria, respectively)and partition into the archaeon-eukaryote part of the tree (Fig.5A). No alternative topology placing Rickettsia with other ${alpha}$ -Proteobacteriaand/or Clostridium with bacilli showed a detectable degree ofsupport in statistical tests. The apparent origin of many bacterialIleRS proteins via HGT from eukaryotes is well documented (10,49, 50); the resulting inconsistency of gene-to-gene distancesis readily captured by the analysis presented here.

View larger version (20K):
[in this window]
[in a new window]
FIG. 5. Maximum-likelihood trees for detected cases of probable XGD. Colors indicate species belonging to the bacterial lineage(s) for which a significant deviation from the clock-like model was detected, namely, ${alpha}$ -Proteobacteria and gram-positive bacteria (A), gram-positive bacteria (B), ${alpha}$ -Proteobacteria (C), and ${gamma}$ -Proteobacteria (D). The triangles represent collapsed clades (the numbers of proteins are indicated in brackets). (A) COG0060 (isoleucyl-tRNA synthetase, IleS). The following proteins and species are included: RP617 of R. prowazekii, RC0953 of R. conorii, CC0701of C. crescentus, BMEII1043 of B. melitensis, AGc1230 of A. tumefaciens, SMc00908 of S. meliloti, mlr8250 of M. loti, CAC3038 of C. acetobutylicum, SA1036 of S. aureus, BS_ileS of B. subtilis, BH2545 of B. halodurans, lin2127 of L. innocua, L0350 of L. lactis, SPy1513 of S. pyogenes, SP1659 of S. pneumoniae, FN0067 of F. nucleatum, TM1361 of Thermotoga maritima, aq_305 of Aquifex aeolicus, DR1335 of Deinococcus radiodurans, BB0833 of Borrelia burgdorferi, and TP0452 of Treponema pallidum. (B) COG1217 (predicted membrane GTPase involved in stress response, TypA). The following proteins and species are included: CAC1684 of C. acetobutylicum, SA0959 of S. aureus, BS_ylaG of B. subtilis, BH2632 of B. halodurans, lin1055 of L. innocua, L0370 of L. lactis, SPy1527 of S. pyogenes, SP0681 of S. pneumoniae, FN0634 of F. nucleatum, and DR1198 of D. radiodurans. (C) COG1282 (NAD/NADP transhydrogenase beta subunit, PntB). The following proteins and species are included: RP074 of R. prowazekii, RC0104 of R. conorii, CC3303 of C. crescentus, BMEII0325 of B. melitensis, AGc4531 of A. tumefaciens, SMc03938 of S. meliloti, mlr5184 of M. loti, slr1434 of Synechocystis sp., and all3408 of Nostoc sp. (D) COG1526 (uncharacterized protein required for formate dehydrogenase activity, FdhD). The following proteins and species are included: FdhD of E. coli K-12, ZfdhD of E. coli O157:H7, ECs4821 of E. coli O157:H7 EDL933, STM4038 of S. enterica serovar Typhimurium, YPO4060 of Y. pestis, HI0005 of H. influenzae, PM0410 of P. multocida, VC1519 of V. cholerae, PA5180 of Pseudomonas aeruginosa, RSp1050 and RSc2367 of Ralstonia solanacearum, CC1242 of C. crescentus, and Cj1508c of Campylobacter jejuni.

COG1217 (predicted membrane GTPase involved in stress response,TypA) is another case of the C. acetobutylicum sequence beingseparated from the sequences of the rest of the Bacillus-Clostridiumgroup. The C. acetobutylicum protein clusters with the proteinfrom Fusobacterium nucleatum (an arrangement which is observedin phylogenetic trees of many genes [Wolf, unpublished data])and is separated from its cognate group by the cyanobacterial,actinobacterial, and Deinococcus branches (Fig. 5B). Althoughthe support for this topology was not overwhelming (expectedlikelihood weight, 0.67 to 0.89 for the topology shown in Fig.5B versus 0.10 to 0.29 for the alternative topology placingboth Clostridium and Fusobacterium at the base of the gram-positivebacteria), HGT from the Fusobacterium lineage to Clostridiumresulting in XGD seems to be the best interpretation of thistree.

In COG1282 (NAD/NADP transhydrogenase beta subunit, PntB) theAgrobacterium sequence clusters with ${gamma}$ -Proteobacteria and Neisseria,whereas the rest of the ${alpha}$ -Proteobacteria cluster with Ralstonia(Fig. 5C). No alternative topology is viable according to thetests performed. Thus, in Agrobacterium, the original ${alpha}$ -proteobacterialgene apparently has been displaced with the ${gamma}$ -proteobacterialortholog.

COG1526 (uncharacterized protein required for formate dehydrogenaseactivity, FdhD) exemplifies the frequently observed separationof Vibrio from the rest of the ${gamma}$ -Proteobacteria (Fig. 5D). Inthis case, the Vibrio protein is closely related to the Ralstoniaortholog, and the Vibrio-Ralstonia clade joins the gram-positivebranch instead of the proteobacterial branch. No topology joiningVibrio with other ${gamma}$ -Proteobacteria is supported. In this case,at least two HGT events have to be inferred: displacement ofthe original proteobacterial gene in either the Vibrio or theRalstonia lineage by the ortholog from gram-positive bacteriaand a secondary HGT between Vibrio and Ralstonia.

We also examined whether the presence a molecular clock violationin a COG is correlated for the three bacterial lineages studied.In the set of 108 conserved COGs, the Pearson correlation coefficientswere –0.02 between ${gamma}$ - and ${alpha}$ -Proteobacteria, –0.11between ${gamma}$ -Proteobacteria and gram-positive bacteria, and 0.17between ${alpha}$ -Proteobacteria and gram-positive bacteria. The distributionof the number of anomalies detected in each COG (range, 0 to3) nearly perfectly agreed with the expectation under the independencehypothesis [P( ${chi}$ ²) > 0.5]. Thus, somewhat counterintuitively,HGT appeared to occur independently in different lineages, andthere was no obvious HGT propensity that could be considereda characteristic of evolution of an entire COG. Nor did we detectany significant connection between the violations of the molecularclock detected and the relative evolution rates. The distributionsof v values were very similar for split and nonsplit sets, andthe difference between the means was insignificant for all threegroups according to the t test (data not shown).

General discussion and conclusions. Comparative genomics allows researchers to restate old questionsof evolutionary biology on a grander scale and, perhaps, ina more biologically meaningful way. Thus, rather than analyzingthe molecular clock for one or a few selected protein families,it is possible to address the issue at the level of completegenomes and to ask what fraction of the genes follow the clock-likemodel of evolution and which genes demonstrably deviate fromit. Here we describe a simple theoretical framework that allowedus to classify orthologous sets of bacterial genes (COGs) intothese two categories. A similar relative rate test was describedby Syvanen in his analysis of the acceleration of evolutionof rRNA in eukaryotes (45).

We found that for several hundred COGs analyzed representingthree well-defined bacterial lineages, the ${alpha}$ -Proteobacteria,the ${gamma}$ -Proteobacteria, and the Bacillus-Clostridium group, theclock-like null hypothesis could not be rejected for $~$ 70% ofthe COGs, whereas the rest showed substantial anomalies. Itshould be noted that the null hypothesis employed here is, infact, a soft clock, which, strictly speaking, does not requireconstancy of evolutionary rates for the genes analyzed. Allthat is required is that the rate distribution remains constant(i.e., all genes are allowed to accelerate or decelerate synchronously)(21). Notably, we also found that, within the set of 108 COGsthat were represented by a single ortholog in all species analyzedfrom the three lineages, the relative evolutionary rates werestrongly correlated among the lineages. This observation emphasizesthe general validity of the soft genomic clock.

We also analyzed the nature of the observed anomalies and foundthat, for the most conspicuous anomalies, the majority weremost readily explained by HGT from phylogenetically distantlineages. Importantly, this is a conservative estimate becausewe analyzed only a special set of well-behaved COGs, which containedexactly one ortholog from each of the species included. HGTevents have been classified into three broad categories: (i)acquisition of genes new to the recipient lineage, (ii) acquisitionof paralogs of resident genes, and (iii) XGD, in which a residentgene is displaced by an ortholog from a different lineage (30).The present analysis was designed to identify only cases ofXGD (although in some exceptional situations we obtained indicationsof paralog acquisition where the COG analyzed seemed to containhidden paralogs). The results suggest that XGD occurs duringevolution of $~$ 10 to 15% of the bacterial genes. This relativelylow fraction of HGT among single-ortholog bacterial genes iscompatible with the notion that, at least within well-definedclades, such as the ${gamma}$ -Proteobacteria, the ${alpha}$ -Proteobacteria, andthe Bacillus-Clostridium group, these genes may be combinedto produce organismal phylogenies, preferably after exclusionof the genes with detected HGT (9, 33, 51). However, the fractionof likely HGT detected here is considerably greater than thatrecently reported for single-ortholog gene sets from ${gamma}$ -Proteobacteria(33). It seems likely that the difference is a consequence ofthe criterion used for selection of orthologous gene sets inthe latter study, in which only genes that are highly conservedwithin ${gamma}$ -Proteobacteria were examined. This criterion probablyresulted in exclusion of orthologous sets deviating from theclock-like behavior.

An unexpected finding of this study is the lack of significantcorrelation among the three bacterial lineages analyzed withrespect to the deviations from the clock model and the probableoccurrence of XGD. An observation of such a deviation in anyone of the lineages was a poor predictor of deviations in otherlineages. This result seems to be poorly compatible with thecomplexity hypothesis (25) and similar notions concerning thedependence of HGT on biological function and is more in linewith the ideas on random scatter and lineage-specific trendsof HGT events (55). It should be noted that we analyzed a limitedgene set and only one form of HGT (XGD). Functional correlatesare likely to emerge in larger-scale studies, but our presentresults indicate that these connections are far from being absolute.

The approach described here allows rapid identification of orthologousgene sets whose evolution significantly deviates from the softclock model. Interpretation of these deviations as lineage-specificacceleration of evolution, XGD, or a combination of the tworequires detailed phylogenetic analysis. Nevertheless, we believethat this methodology has its own advantages and could be usefulin the study of genome-wide evolutionary trends. In particular,this approach allows workers to detect significant deviationsfrom the clock-like model of evolution in a particular lineagewithout using any information on species outside that lineage,such as the (often unknown) source of the potential HGT. Morepractically, the procedure described here could be suitablefor removing anomalous COGs from multigene sets employed forconstruction of organismal phylogenies.

ACKNOWLEDGMENTS

We thank Eva Czabarka and Georgy Karev (National Center forBiotechnology Information) for helpful discussions on the statisticalanalysis of the data.

P.S.N., M.S.G., and A.A.M. were partially supported by grantsfrom the Howard Hughes Medical Institute (grant 55000309), theprograms "Molecular and Cellular Biology" and "Origin and Evolutionof the Biosphere" of the Russian Academy of Sciences, and theFund for Support of Russian Science (MSG).

FOOTNOTES

* Corresponding author. Mailing address: National Center for Biotechnology Information, National Library of Medicine, National Institutes of Health, Bethesda, MD 20894. Phone: (301) 435-5913. Fax: (301) 435-7794. E-mail: koonin{at}ncbi.nlm.nih.gov.

REFERENCES

Top