BACKGROUND
comparison of completely sequenced genomes belonging to a wide variety of prokaryotic species has revealed how fluid these genomes are .
for instance, synteny is rapidly lost  <cit>  as a consequence of incessant gene flows inside a genome and between species.
moreover, there is frequent gain and loss of genes, even between various strains of the same species or species belonging to the same genus .
however, a small proportion of orthologous genes, hereafter called positional orthologous genes , conserve their local order even between distantly related species  <cit> .
such observations gave birth to the concept of genomic context  <cit> .
accordingly, this rare conservation of gene order has been interpreted as the signature of functional relationships between the products of these stably associated genes  <cit> .
detecting these stable islands in a sea of moving genes is a complex and difficult problem in the field of comparative genomics.
the potential pitfalls increase with the deluge of newly published genome sequences.
several groups have described these difficulties, the primary one being reliably determining the orthologs among the whole set of homologs detected by exhaustive comparisons between each pair of completely sequenced genomes.
various strategies have been applied in the past to detect orthologs.
the most popular one has been the so-called bidirectional best hit method.
it is based on blast reciprocal-best-hits  between either a pair  <cit>  or a triplet of genomes  <cit> .
however, it has been shown that blast searches often return as the highest scoring hit a protein that is not the nearest phylogenetic neighbour of the query sequence  <cit> .
as a consequence, it has been repeatedly emphasized that actual orthologs could be missed, including when using rbh approaches .
according to fulton et al.
<cit>  roughly  <dig> in  <dig> rbh-predicted rat-mouse orthologs are false positives  and about  <dig> in  <dig> rbh-predicted orthologs for two pseudomonas species are probably similarly incorrect.
note however that there is presently no direct experimental validation of the effects of the koski-golding demonstration on the prediction accuracy of blast reciprocal-best-hits  approaches.
to alleviate these problems in detection of orthologs, it seems sounder to use more theoretically consistent approaches.
wall and coworkers  <cit>  have proposed a new method called rsd  that is based on estimates of evolutionary distance rather than blast scores.
other approaches based on phylogeny criteria have been proposed for eukaryotic genomes  <cit> , and more recently  for prokaryotic data  <cit> .
this work used two alternative and complementary approaches to reliably determine orthologs from the set of homologous genes detected after an exhaustive genomic comparison.
we first adapted the rsd method  <cit>  to our pipeline of programs for defining synteny blocks.
we then introduced a new automatic approach for differentiating between orthologs and paralogs based upon automated generation of phylogenetic trees for families of homologous genes in prokaryotes.
we compare the advantages and disadvantages of both the rsd and phylogeny approaches and show that, while neither method alone is always reliable, their union is giving reasonably complementary data.
this allows to recognize all bona fide pogs and to further identify all syntenic blocks whatever the phylogenetic distance separating the studied species.
recognizing such blocks of adjacent genes is essential to study them as ancient evolutionary units.
indeed, we further showed using statistical comparisons that some selective pressure is keeping pogs associated in a large majority of genomes.
we further investigated whether pogs afford a particular mode of evolution that would be different from that of other orthologous genes which are free to move independently.
RESULTS
an outline of our experimental approach
we designed a three-step experimental scheme to find the conserved syntenic blocks in a large assortment of  <dig> microorganisms  located at various distances from one another in the taxonomy space.
in step one, we compared the proteome  of each organism with all other proteomes to identify the full set of homologous proteins they share.
in step two, we developed new approaches for finding valid orthologs among all the homologs found in a comparison of each genome with all the others.
step three searches for a minimal conserved unit, i.e.
a pair of pogs that are adjacent in at least two genomes.
after extending these minimal units as far as possible, it becomes feasible to assess the relative number and size of synteny blocks in both nearby and distant species.
collecting homologous genes in  <dig> microorganisms
in this first step, we compared each organism's proteome with those of all the others to identify the full set of homologous proteins they share.
accordingly, we made  <dig> comparisons between every pairing of the whole set of  <dig> prokaryotic species  listed in table  <dig>  to ensure that the search for homologs was both exhaustive and exact, regardless of the taxonomic distances separating the compared species, we used the darwin allall program  <cit>  to compare each pair of proteomes.
as previously shown  <cit> , the maximum likelihood approach used by the darwin programs efficiently detects in one step all segments of homology, including those between evolutionarily-distant proteins  <cit> .
these features are crucial in work such as this that compares a large number of genomes separated by a wide range of phylogenetic distances.
an additional advantage of the darwin programs is that the evolutionary distance between similar proteins in species with a common ancestor uses the well-established pam units, a pam unit being defined as the number of accepted point mutations per  <dig> residues separating two sequences  <cit> .
the frequency with which any particular pair of  amino acids occurs at a given position in two properly-aligned homologous proteins can be used as a pam score to evaluate the evolutionary distance separating the two proteins.
accordingly, as previously described  <cit>  and detailed in methods, we applied the rules proposed by altschul  <cit>  to define how significant sequence similarities are.
we required that the pam distance separating a pair of homologs be less than  <dig> units and that each segment of homology have a minimal size of  <dig> residues.
moreover, in this search for homologs belonging to syntenic blocks, we further required that the pair of aligned amino acid sequences extend at least 80% of the length of the shorter matching protein.
the allall homolog search compared each of the  <dig>  proteins to the whole set encoded by the  <dig> prokaryotic genomes .
we found that the vast majority  of these proteins have at least one homolog in another microorganism, probably because many compared species have at least one  close taxonomic relative in the data set .
reliable determination of orthologous genes
in this second step, we tried to recognize a reliable set of orthologs among the full set of  <dig>  homologs the first step had found.
as explained in the introduction, this is not a trivial task; so we compared the results of two alternative strategies to avoid the pitfalls of applying each method independently.
identifying orthologs using a reciprocal smallest distance  approach
we used a reciprocal evolutionary distance to avoid the risks of errors inherent in similar previous approaches explained in the introduction.
we took advantage of specific properties of the darwin allall program  <cit>  that set it apart from blast algorithm.
also, the allall program calculates the pam evolutionary distances separating homologous genes, retaining the shortest distance found for each analyzed pair of proteins.
these features made it trivial to determine rsd ortholog pairs in a comparison of two proteomes, using a strategy similar to that of wall and coworkers  <cit> .
the rsd ortholog pair between protein a encoded by genome ga and protein b encoded by genome gb, was the one for which the pam distance separating a from b was smaller than that separating a from any other protein encoded by gb and b from any other protein encoded by ga.
this rsd approach produced a list of  <dig>  orthologs amounting to 69% of the proteins encoded by the  <dig> genomes studied in this paper.
these  <dig>  orthologs form a total of  <dig> , <dig> pairs .
however, we found that this one-to-one approach would miss a significant number of orthologous relationships when gene duplications arose after a recent speciation event.
therefore, we developed a second approach based on phylogenetic analyses that detects these so-called in-paralogs  <cit> .
identifying orthologs by a phylogenetic approach
this second approach analyses the evolutionary features of the detected homologous proteins after grouping them in families.
a transitive approach based on a graph algorithm for determining connected components linked the  <dig>  homologous proteins into  <dig>  families.
unfortunately, we observed that a few of these families were enormous, heterogeneous groupings that, in a significant number of cases, consisted of two or more groups of highly-connected proteins linked by a few probably-spurious edges.
although many of these bridges are indicating a biologically meaningful paralogy, they appear detrimental to our search for consistent and valid families required to identify genuine orthologous relationships.
to avoid any disadvantageous groupings of paralogous families and to reduce the size of too-large and heterogeneous families, we further applied a graph algorithm for bridge detection, and broke many unwanted bridges .
applying this strategy to the  <dig>  transitive families redistributed all the homologous genes into  <dig>  better defined families.
after this step, one extremely large protein "family" remained.
this heterogeneous cluster contained  <dig>  members that are mainly hydrophobic proteins such as transporters and other membrane proteins, including many  proteins with unknown function.
such a huge gathering of disparate proteins is biologically meaningless.
moreover, it was clearly worthless to analyze it by a tree approach due to its size and complexity .
therefore, we applied the mcl algorithm  <cit>  to break up this huge and heterogeneous cluster.
this left  <dig>  homologous proteins that are members of  <dig>  families:  <dig>  pairs of proteins  and all the homologous proteins belonging to each of the remaining  <dig>  families of size  <dig> or larger – including the largest one that contains  <dig> proteins.
after reconstructing a phylogenetic tree for each family, we applied an ad hoc algorithm  to identify the orthologous relationships within each phylogenetic tree, as shown in fig.
<dig> and detailed in methods.
as it has been repeteadly shown, there is no prokaryotic species tree that is recognized as a consensus.
accordingly, we used an approach that does not require such a model tree.
briefly, the treeortho algorithm analyses the leaves of a midpoint-rooted tree to determining for each node – whichever its deepness – whether its descendants arose from a gene duplication  or a speciation event .
for instance, the deepest node of the simplified tree in fig.
<dig> is viewed as a duplication event  because the two branches stemming from this node display leaves that correspond to the same set of species.
the descendents of the immediate next nodes being in separate species, these nodes identify speciation events .
probing all the family trees that we could analyze by this method recovered a total of  <dig>  proteins that have an ortholog in at least one another species.
in total, this phylogeny approach yielded a total of  <dig> , <dig> orthologous pairs .
comparing the rsd and phylogenetic approaches
fig.
<dig> shows a tree which illustrates the pros and cons of both approaches.
family  <dig> groups chemotaxis proteins chec in different euryarchaeota, proteobacteria, firmicutes, one actinobacteria and thermotoga maritima.
consider first the search for the orthologs of the b. subtilis chec protein .
the rsd approach finds only  <dig> orthologs , discarding the paralogs present in thermococcus kodakarensis, pyrococcus abyssi and the halophilic archaea haloarcula and halobacterium .
the phylogeny approach detects all these in-paralogs, going up to  <dig> retained proteins .
thus, in this case the phylogeny approach appears to be far more efficient, as expected.
let's take now the search for the orthologs of the chec protein in archaeoglobus fulgidus .
here the rsd approach finds  <dig> orthologs , and the phylogeny approach  <dig>  but only  <dig>  are found by both approaches.
by definition, seven paralogs present in the different archaeal species  were excluded by the rsd approach but included in the harvest of the phylogeny approach after identification of the nodes corresponding to a duplication event  and those corresponding to a speciation event.
however, the presence in the dataset of two different sequences  for the same organism  led the treeortho algorithm to identify as a gene duplication node  the separation between the two subtrees containing the respective m. acetivorans sequences.
as a consequence, all the nine descendents of the ancestral copy that led to the m. acetivorans chec <dig> sequence were excluded from the phylogenetic tally  as corresponding to a paralogy event.
these descendants are found by the rsd approach because there are only one remaining copy in arcfu  and one remaining copy  in each of the seven different bacteria branching in the chec2_metac subtree.
note that the ancestral chec duplication appears rather ancient as shown in a tree containing all archaeal chec sequences and rooted with a few bacterial homologs .
the two copies have been separated by at least three speciation events, the chec <dig> copy having highly diverged  and appearing now more similar to a group of bacterial chec.
this large divergence is further underlined by the respective gene context of both copies of chec.
although, there was an ancestral duplication of the whole operon, we found that the adjacency ched <dig> – chec <dig> is conserved in the various phylogenetically close bacteria .
on the contrary, gene chec <dig> is found adjacent to chea <dig> in m. acetivorans as it is in the other euryarchaeota.
such correct algorithmic interpretations that lead to reject rsd orthologs as being actual paralogs probably explain many of the cases of sequences that are discarded by phylogeny but recovered by rsd.
such instances could not occur in the  <dig> low-complexity trees that are made uniquely of orthologs, without paralogs, but may happen when analysing the  <dig> trees of medium complexity containing only orthologs and in-paralogs and the  <dig> trees of high complexity that contain in- and out-paralogs, as well as orthologs.
combining data obtained with the rsd and phylogenetic approaches
both approaches, thus, yield incomplete data.
the union of both result sets provides a total of  <dig> , <dig> pairs of orthologs .
because of the potential for confusion from the false positives returned by both methods, it must be noted that their intersection left  <dig> , <dig> pairs of orthologs that are endowed with a high degree of confidence that they are genuine.
however, such a stringent filtering removes all in-paralogs.
identifying pairs of adjacent orthologous genes
given the lists of bona fide orthologs, our third step is to determine how many form reciprocal pairs of strictly adjacent genes in at least two genomes.
we found that the most efficient and rapid way to collect all these pairs was the following strategy.
we stored all the data obtained from the comparison of the  <dig> pairs of genomes in a relational database syntebase  and we designed a specific sql request  to query syntebase to detect the quadruplets of proteins containing two pairs of orthologs that are encoded by genes which are adjacent in both genomes.
table  <dig> shows that we detected  <dig>  synteny blocks of size  <dig> with the rsd method and  <dig>  ones with the phylogenetic analysis.
the slight excess of synteny blocks found by the phylogeny approach is probably due to the detection of remote homologs in families that could not be found in rsd pairs due to our initial alignment threshold.
this is most probably the case, for instance, of the chec sequences encoded by shewanella oneidensis  and vibrio cholerae  which are not detected by rsd but found by tree analysis as orthologs of b. subtilis chec .
table  <dig> further shows that we get more pairs of orthologs with rsd but fewer blocks of synteny, underlining that some rsd-orthologs are most probably ancient paralogs .
summing all data obtained with both methods provides a total of  <dig>  blocks of two adjacent orthologs .
finally, we attempted to extend the detected neighbourhood relationships to larger syntenic blocks.
the quadruplets, stored in syntebase, were analyzed with synblock, an algorithm we designed to map the synteny blocks of size  <dig> and to merge those that share a common pair of orthologous genes.
note that, in this work where we aim to find out the conserved blocks of genes that reflect the ancestral gene order, we required strict gene adjacency, forbidding any insertion in a synteny block of a gene which would be unique to one of the pairs of compared genomes.
we further used these data entered into syntebase to make some statistical studies on the respective level of conservation of neighbourhood and other evolutionary features linked to gene order.
assessing the evolution rates of positional orthologous genes
measuring the distribution of the pam distances in the taxonomic space
the pam distance that separates two orthologs in various pairs of genomes was used to compare the rate of evolution of different classes of orthologs with different respects.
first, fig.
3a shows the distribution of the pam distances separating orthologous pairs when comparing e. coli with closely related enterobacteriaceae s. enterica, the firmicute b. subtilis, the member of the bacteroidetes/chlorobi group bacteroides thetaiotaomicron and the archaeon m. acetivorans.
as expected, this distribution is very dependent of the taxonomic proximity between species.
this is confirmed in fig.
3b and table  <dig> that display the taxonomic distributions of the mean of the pam distances separating pairs of orthologs and of the total number of orthologous pairs in the different genomic comparisons.
<dig> proteome size of e. coli k12: 4279
measuring the evolution rates of the different classes of orthologous genes
to check if the observed difference between the pam distances separating orthologous genes located inside and outside synteny blocks  is statistically significant for any pair of species, we further tested the whole set of compared genomes, using a bootstrap sampling approach as detailed in methods.
the null hypothesis  assumes that proteins encoded by pogs evolve at the same evolutionary rate as those encoded by orthologous genes located outside the synteny blocks.
the alternative hypothesis  assumes that the proteins encoded by orthologous genes located inside synteny blocks evolve more slowly than the proteins encoded by orthologous genes located outside synteny blocks.
to validate this study, we first used the phylogenetic approach, taking into account orthologs contained in trees of various complexities defined as follows.
we analyzed trees made up uniquely of orthologs , and trees containing only orthologs and in-paralogs .
we then did similar computations for the  <dig> genomes with the entire dataset obtained using the rsd ortholog detection method.
fig.
<dig> shows that in both rsd and phylogeny approaches, the overwhelming majority of the statistical tests reject the null hypothesis h <dig> in favour of h <dig>  we found that the remaining, untested cases  correspond to a few comparisons that are too small to be safely used in this bootstrap test.
this is the case, for example, of all comparisons  involving nanoarchaea equitans, an archaeon with a very small genome that is also very distant from all other species, including the other archaea.
we conclude that there is a universal trend in the evolution of prokaryotic genomes: genes located in regions with preserved gene order  evolve more slowly than genes found outside.
discussion
as soon as several prokaryotic genome sequences were completed, it became clear that very few of their genes had retained their local ordering through evolution.
this observation was reinforced as more and more genomes belonging to a wide spectrum of organisms – from very close to very distantly-related – were published.
this led to the concept of gene  context  <cit> .
as we prepare to be inundated by thousands of completely sequenced genomes in the next few years  <cit> , it becomes urgent to design a strategy for quickly evaluating the conservation of gene order in newly-published genomes, comparing them to all those already analyzed.
this paper describes a suite of programs  that automate all the essential steps for evaluating gene order conservation in prokaryotes, as well as presenting some results from their initial application.
after identifying in step  <dig>  the whole set of homologs shared by  <dig> prokaryotic species , we designed two independent strategies  to extract from this set the collection of orthologs that are shared by at least two genomes.
in a third step , we located the positional orthologous genes  in this collection, and we studied the specific properties that distinguish them from other orthologs.
correctly identifying positional orthologs
this is a crucial step, because mixing paralogs and orthologs can lead to errors in predicting ortholog pairs that seriously compromise studies of the mode of evolution of genes and genomes.
this is especially true when studying prokaryotes.
we already showed  <cit>  that many gene families stem from ancestral gene duplication events that occurred well before the divergence of present-day species.
therefore, orthologs are often more similar between closely related species, allowing them to be distinguished from paralogs.
however, when we compare genomes belonging to distant species, the difference in percent identity between orthologs and paralogs falls rapidly , and the risk of confusing errors quickly increases.
although there is a long history of research on identifying orthologs, this is an extremely active research area at the moment .
numerous papers appeared within months, each one claiming to present the "best method" .
a consensus does seem to be emerging about the superiority of approaches based on evolutionary studies over best bidirectional blast hits .
anyhow, the fact that so many "best methods" were published in a single year clearly underlines that the field is immature.
because no single approach is superior in all cases, our work applies two alternative and complementary techniques, reciprocal smallest distance  and phylogenetic approaches.
both approaches are based on well-established pam evolutionary distances, computed by a maximum likelihood approach  <cit> .
using these pam distances allows directly determining the ortholog pairs separated by rsd, comparing each homolog with all others by means similar to the method already proposed by wall et al.
<cit> .
in our  <dig> genomes, this method identified  <dig>  orthologs  forming  <dig> , <dig> pairs.
we expected that the phylogeny approach would systematically increase these figures; since it should collect not only all pairs of orthologs found by rsd but also so-called in-paralogs  <cit>  the rsd approach misses.
in fact, this does happen in many cases.
for example, as fig.
<dig> shows for orthologs of the b.subtilis chemotaxis chec protein, we increase the number from  <dig> found using the rsd approach to  <dig> with the phylogeny approach.
unfortunately, although this phylogeny strategy was found to be quite effective, it can create unexpected difficulties when comparing many protein sequences that form some heterogeneous and very large "families".
in fact, these huge agglomerates contained about the half of all found homologs.
we managed to decompose many of these into more valid and much smaller families.
however, this decomposition step caused the loss of a significant number of genuine pairs of orthologs.
on the other hand, phylogenetic analysis is advantageous when it rejects a significant number of homologous pairs that have been found as orthologs by the rsd approach, as being actual paralogs.
this will occur for instance when a deep node is correctly viewed as due to an ancestral gene duplication event .
moreover, the apparent excess of pairs of orthologs found by the rsd method  looks suspicious; since a significant number of them are discarded when looking at the number of pogs found by this approach.
therefore, it seems that phylogeny approach is more accurate to detect bona fide orthologs, excluding pseudoorthologs  and including in-paralogs.
however, it appears presently hazardous to choose a single "best method" to identify orthologs.
our present strategy has been to combine the results of our two alternative methods.
uses of the data fall between two possible options extremes: using the union and using the intersection of their results.
where the goal is to collect a maximum of putative synteny cases in the maximum number of genomes, it is appropriate to merge the orthologs determined by both methods.
alternatively, if the goal is to work with pairs of orthologs that are highly likely to be authentic, it would be best to use only the results that are common to both approaches, but such an intersection will unfortunately remove all in-paralogs.
evolutionary constraints on pogs
based upon these identified synteny blocks and the computed  evolutionary distances separating the proteins encoded by these adjacent genes , it was possible to compare the rates of evolution of the orthologous genes located inside and outside synteny blocks, respectively.
as shown in table  <dig> and figs.
<dig> and  <dig>  proteins encoded by genes located within synteny blocks are found to evolve at a lower rate than those encoded by genes located outside these blocks.
this observed difference appears to be consistent with the concept of genomic context in which it is assumed that neighboring genes are under a strong evolutionary pressure to maintain their adjacency.
such a pressure would exist, for example, if the proteins they encode have to interact in order to participate to the same biological process.
however, in most cases it is not yet clear what the actual cause of low evolutionary rate is.
our demonstration that pogs are more evolutionarily constrained than other orthologs that do not reside in synteny blocks confirms and extends the previous observation made by dandekar et al.
<cit>  that "the degree of sequence conservation is on average substantially higher than that in genes that do not exist as conserved genes pairs".
moreover, our demonstration agrees with both the variables and universals  that have previously been associated with differences in protein evolutionary rates .
we expect that persistently adjacent genes result from several of these universals:  they are probably essential since they are conserved in large parts  of the taxonomic space.
they will be co-transcribed  at the same level of expression.
in many cases, their products will interact in the cell, requiring a significant level of connectivity between them and strict conservation of residues to permit a sufficient density of structural and functional contacts.
we performed a preliminary evaluation of these hypotheses using experimental data obtained from the model organisms e. coli and b. subtilis.
of the  <dig> pogs shared by these two species,  <dig>  are essential e. coli genes,  <dig> non-essential and  <dig> unknown.
moreover,  <dig> of the shared e. coli and b. subtilis pogs  belong to the e. coli interactome and  <dig> of these are essential genes .
in total, around 23% of the  <dig> pogs e. coli shares with b. subtilis are both essential and part of its interactome, making them what have been called hub proteins  <cit> .
the different data obtained in this work suggest that blocks of adjacent pogs which display conserved co-occurence in various present-day genomes can be interpreted as ancient evolutionary units that have resisted over the years to the imperious fluidity forces.
thus, these evolutionary units could be viewed as the remnants of the skeleton of the gene order found in a common ancestor to all compared present-day microorganisms.
we further propose that pogs have a low evolutionary rate because they are constituents of functional modules as previously defined by hartwell et al.
<cit>  in a systems biology perspective .
being a part of such modules would require a coordinated evolution between their sequence and expression level to participate in a collective cell function.
CONCLUSIONS
maintaining the order of genes whose products interact in the same cell function is clearly beneficial to organisms, for instance to allow co-expression and co-regulation.
this requires the creation of clusters that have low recombination rates  and that the adjacency of functionally related genes is maintained by epistatic selection  <cit> .
however, as we compare species that are distant in the taxonomy space, there is an inverse relation between taxonomic proximity and the size of the synteny blocks.
this evidence agrees with poyatos and hurst  <cit>  contention that gene order is forever in flux.
therefore, the exceptional instances where a stable gene order is conserved, even for a few gene pairs, implies strong selective pressure that should be reflected in high levels of co-expression of essential products such as one sees, for instance, in the exemplary case of ribosomal proteins.

 comparison of completely sequenced genomes belonging to a wide variety of prokaryotic species has revealed how fluid these genomes are. for instance, synteny is rapidly lost as a consequence of incessant gene flows inside a genome and between species. moreover, there is frequent gain and loss of genes, even between various strains of the same species or species belonging to the same genus. however, a small proportion of orthologous genes, hereafter called positional orthologous genes, conserve their local order even between distantly related species. such observations gave birth to the concept of genomic context. accordingly, this rare conservation of gene order has been interpreted as the signature of functional relationships between the products of these stably associated genes. detecting these stable islands in a sea of moving genes is a complex and difficult problem in the field of comparative genomics. the potential pitfalls increase with the deluge of newly published genome sequences. several groups have described these difficulties, the primary one being reliably determining the orthologs among the whole set of homologs detected by exhaustive comparisons between each pair of completely sequenced genomes. various strategies have been applied in the past to detect orthologs. the most popular one has been the so-called bidirectional best hit method. it is based on blast reciprocal-best-hits between either a pair or a triplet of genomes. however, it has been shown that blast searches often return as the highest scoring hit a protein that is not the nearest phylogenetic neighbour of the query sequence. as a consequence, it has been repeatedly emphasized that actual orthologs could be missed, including when using rbh approaches. according to fulton et al. roughly in rbh-predicted rat-mouse orthologs are false positives and about in rbh-predicted orthologs for two pseudomonas species are probably similarly incorrect. note however that there is presently no direct experimental validation of the effects of the koski-golding demonstration on the prediction accuracy of blast reciprocal-best-hits approaches. to alleviate these problems in detection of orthologs, it seems sounder to use more theoretically consistent approaches. wall and coworkers have proposed a new method called rsd that is based on estimates of evolutionary distance rather than blast scores. other approaches based on phylogeny criteria have been proposed for eukaryotic genomes, and more recently for prokaryotic data. this work used two alternative and complementary approaches to reliably determine orthologs from the set of homologous genes detected after an exhaustive genomic comparison. we first adapted the rsd method to our pipeline of programs for defining synteny blocks. we then introduced a new automatic approach for differentiating between orthologs and paralogs based upon automated generation of phylogenetic trees for families of homologous genes in prokaryotes. we compare the advantages and disadvantages of both the rsd and phylogeny approaches and show that, while neither method alone is always reliable, their union is giving reasonably complementary data. this allows to recognize all bona fide pogs and to further identify all syntenic blocks whatever the phylogenetic distance separating the studied species. recognizing such blocks of adjacent genes is essential to study them as ancient evolutionary units. indeed, we further showed using statistical comparisons that some selective pressure is keeping pogs associated in a large majority of genomes. we further investigated whether pogs afford a particular mode of evolution that would be different from that of other orthologous genes which are free to move independently.  an outline of our experimental approach we designed a three-step experimental scheme to find the conserved syntenic blocks in a large assortment of microorganisms located at various distances from one another in the taxonomy space. in step one, we compared the proteome of each organism with all other proteomes to identify the full set of homologous proteins they share. in step two, we developed new approaches for finding valid orthologs among all the homologs found in a comparison of each genome with all the others. step three searches for a minimal conserved unit, i.e. a pair of pogs that are adjacent in at least two genomes. after extending these minimal units as far as possible, it becomes feasible to assess the relative number and size of synteny blocks in both nearby and distant species. collecting homologous genes in microorganisms in this first step, we compared each organism's proteome with those of all the others to identify the full set of homologous proteins they share. accordingly, we made comparisons between every pairing of the whole set of prokaryotic species listed in table to ensure that the search for homologs was both exhaustive and exact, regardless of the taxonomic distances separating the compared species, we used the darwin allall program to compare each pair of proteomes. as previously shown, the maximum likelihood approach used by the darwin programs efficiently detects in one step all segments of homology, including those between evolutionarily-distant proteins. these features are crucial in work such as this that compares a large number of genomes separated by a wide range of phylogenetic distances. an additional advantage of the darwin programs is that the evolutionary distance between similar proteins in species with a common ancestor uses the well-established pam units, a pam unit being defined as the number of accepted point mutations per residues separating two sequences. the frequency with which any particular pair of amino acids occurs at a given position in two properly-aligned homologous proteins can be used as a pam score to evaluate the evolutionary distance separating the two proteins. accordingly, as previously described and detailed in methods, we applied the rules proposed by altschul to define how significant sequence similarities are. we required that the pam distance separating a pair of homologs be less than units and that each segment of homology have a minimal size of residues. moreover, in this search for homologs belonging to syntenic blocks, we further required that the pair of aligned amino acid sequences extend at least 80% of the length of the shorter matching protein. the allall homolog search compared each of the proteins to the whole set encoded by the prokaryotic genomes. we found that the vast majority of these proteins have at least one homolog in another microorganism, probably because many compared species have at least one close taxonomic relative in the data set. reliable determination of orthologous genes in this second step, we tried to recognize a reliable set of orthologs among the full set of homologs the first step had found. as explained in the introduction, this is not a trivial task; so we compared the results of two alternative strategies to avoid the pitfalls of applying each method independently. identifying orthologs using a reciprocal smallest distance approach we used a reciprocal evolutionary distance to avoid the risks of errors inherent in similar previous approaches explained in the introduction. we took advantage of specific properties of the darwin allall program that set it apart from blast algorithm. also, the allall program calculates the pam evolutionary distances separating homologous genes, retaining the shortest distance found for each analyzed pair of proteins. these features made it trivial to determine rsd ortholog pairs in a comparison of two proteomes, using a strategy similar to that of wall and coworkers. the rsd ortholog pair between protein a encoded by genome ga and protein b encoded by genome gb, was the one for which the pam distance separating a from b was smaller than that separating a from any other protein encoded by gb and b from any other protein encoded by ga. this rsd approach produced a list of orthologs amounting to 69% of the proteins encoded by the genomes studied in this paper. these orthologs form a total of, pairs. however, we found that this one-to-one approach would miss a significant number of orthologous relationships when gene duplications arose after a recent speciation event. therefore, we developed a second approach based on phylogenetic analyses that detects these so-called in-paralogs. identifying orthologs by a phylogenetic approach this second approach analyses the evolutionary features of the detected homologous proteins after grouping them in families. a transitive approach based on a graph algorithm for determining connected components linked the homologous proteins into families. unfortunately, we observed that a few of these families were enormous, heterogeneous groupings that, in a significant number of cases, consisted of two or more groups of highly-connected proteins linked by a few probably-spurious edges. although many of these bridges are indicating a biologically meaningful paralogy, they appear detrimental to our search for consistent and valid families required to identify genuine orthologous relationships. to avoid any disadvantageous groupings of paralogous families and to reduce the size of too-large and heterogeneous families, we further applied a graph algorithm for bridge detection, and broke many unwanted bridges. applying this strategy to the transitive families redistributed all the homologous genes into better defined families. after this step, one extremely large protein "family" remained. this heterogeneous cluster contained members that are mainly hydrophobic proteins such as transporters and other membrane proteins, including many proteins with unknown function. such a huge gathering of disparate proteins is biologically meaningless. moreover, it was clearly worthless to analyze it by a tree approach due to its size and complexity. therefore, we applied the mcl algorithm to break up this huge and heterogeneous cluster. this left homologous proteins that are members of families: pairs of proteins and all the homologous proteins belonging to each of the remaining families of size or larger – including the largest one that contains proteins. after reconstructing a phylogenetic tree for each family, we applied an ad hoc algorithm to identify the orthologous relationships within each phylogenetic tree, as shown in fig. and detailed in methods. as it has been repeteadly shown, there is no prokaryotic species tree that is recognized as a consensus. accordingly, we used an approach that does not require such a model tree. briefly, the treeortho algorithm analyses the leaves of a midpoint-rooted tree to determining for each node – whichever its deepness – whether its descendants arose from a gene duplication or a speciation event. for instance, the deepest node of the simplified tree in fig. is viewed as a duplication event because the two branches stemming from this node display leaves that correspond to the same set of species. the descendents of the immediate next nodes being in separate species, these nodes identify speciation events. probing all the family trees that we could analyze by this method recovered a total of proteins that have an ortholog in at least one another species. in total, this phylogeny approach yielded a total of, orthologous pairs. comparing the rsd and phylogenetic approaches fig. shows a tree which illustrates the pros and cons of both approaches. family groups chemotaxis proteins chec in different euryarchaeota, proteobacteria, firmicutes, one actinobacteria and thermotoga maritima. consider first the search for the orthologs of the b. subtilis chec protein. the rsd approach finds only orthologs, discarding the paralogs present in thermococcus kodakarensis, pyrococcus abyssi and the halophilic archaea haloarcula and halobacterium. the phylogeny approach detects all these in-paralogs, going up to retained proteins. thus, in this case the phylogeny approach appears to be far more efficient, as expected. let's take now the search for the orthologs of the chec protein in archaeoglobus fulgidus. here the rsd approach finds orthologs, and the phylogeny approach but only are found by both approaches. by definition, seven paralogs present in the different archaeal species were excluded by the rsd approach but included in the harvest of the phylogeny approach after identification of the nodes corresponding to a duplication event and those corresponding to a speciation event. however, the presence in the dataset of two different sequences for the same organism led the treeortho algorithm to identify as a gene duplication node the separation between the two subtrees containing the respective m. acetivorans sequences. as a consequence, all the nine descendents of the ancestral copy that led to the m. acetivorans chec sequence were excluded from the phylogenetic tally as corresponding to a paralogy event. these descendants are found by the rsd approach because there are only one remaining copy in arcfu and one remaining copy in each of the seven different bacteria branching in the chec2_metac subtree. note that the ancestral chec duplication appears rather ancient as shown in a tree containing all archaeal chec sequences and rooted with a few bacterial homologs. the two copies have been separated by at least three speciation events, the chec copy having highly diverged and appearing now more similar to a group of bacterial chec. this large divergence is further underlined by the respective gene context of both copies of chec. although, there was an ancestral duplication of the whole operon, we found that the adjacency ched – chec is conserved in the various phylogenetically close bacteria. on the contrary, gene chec is found adjacent to chea in m. acetivorans as it is in the other euryarchaeota. such correct algorithmic interpretations that lead to reject rsd orthologs as being actual paralogs probably explain many of the cases of sequences that are discarded by phylogeny but recovered by rsd. such instances could not occur in the low-complexity trees that are made uniquely of orthologs, without paralogs, but may happen when analysing the trees of medium complexity containing only orthologs and in-paralogs and the trees of high complexity that contain in- and out-paralogs, as well as orthologs. combining data obtained with the rsd and phylogenetic approaches both approaches, thus, yield incomplete data. the union of both result sets provides a total of, pairs of orthologs. because of the potential for confusion from the false positives returned by both methods, it must be noted that their intersection left, pairs of orthologs that are endowed with a high degree of confidence that they are genuine. however, such a stringent filtering removes all in-paralogs. identifying pairs of adjacent orthologous genes given the lists of bona fide orthologs, our third step is to determine how many form reciprocal pairs of strictly adjacent genes in at least two genomes. we found that the most efficient and rapid way to collect all these pairs was the following strategy. we stored all the data obtained from the comparison of the pairs of genomes in a relational database syntebase and we designed a specific sql request to query syntebase to detect the quadruplets of proteins containing two pairs of orthologs that are encoded by genes which are adjacent in both genomes. table shows that we detected synteny blocks of size with the rsd method and ones with the phylogenetic analysis. the slight excess of synteny blocks found by the phylogeny approach is probably due to the detection of remote homologs in families that could not be found in rsd pairs due to our initial alignment threshold. this is most probably the case, for instance, of the chec sequences encoded by shewanella oneidensis and vibrio cholerae which are not detected by rsd but found by tree analysis as orthologs of b. subtilis chec. table further shows that we get more pairs of orthologs with rsd but fewer blocks of synteny, underlining that some rsd-orthologs are most probably ancient paralogs. summing all data obtained with both methods provides a total of blocks of two adjacent orthologs. finally, we attempted to extend the detected neighbourhood relationships to larger syntenic blocks. the quadruplets, stored in syntebase, were analyzed with synblock, an algorithm we designed to map the synteny blocks of size and to merge those that share a common pair of orthologous genes. note that, in this work where we aim to find out the conserved blocks of genes that reflect the ancestral gene order, we required strict gene adjacency, forbidding any insertion in a synteny block of a gene which would be unique to one of the pairs of compared genomes. we further used these data entered into syntebase to make some statistical studies on the respective level of conservation of neighbourhood and other evolutionary features linked to gene order. assessing the evolution rates of positional orthologous genes measuring the distribution of the pam distances in the taxonomic space the pam distance that separates two orthologs in various pairs of genomes was used to compare the rate of evolution of different classes of orthologs with different respects. first, fig. 3a shows the distribution of the pam distances separating orthologous pairs when comparing e. coli with closely related enterobacteriaceae s. enterica, the firmicute b. subtilis, the member of the bacteroidetes/chlorobi group bacteroides thetaiotaomicron and the archaeon m. acetivorans. as expected, this distribution is very dependent of the taxonomic proximity between species. this is confirmed in fig. 3b and table that display the taxonomic distributions of the mean of the pam distances separating pairs of orthologs and of the total number of orthologous pairs in the different genomic comparisons. proteome size of e. coli k12: 4279 measuring the evolution rates of the different classes of orthologous genes to check if the observed difference between the pam distances separating orthologous genes located inside and outside synteny blocks is statistically significant for any pair of species, we further tested the whole set of compared genomes, using a bootstrap sampling approach as detailed in methods. the null hypothesis assumes that proteins encoded by pogs evolve at the same evolutionary rate as those encoded by orthologous genes located outside the synteny blocks. the alternative hypothesis assumes that the proteins encoded by orthologous genes located inside synteny blocks evolve more slowly than the proteins encoded by orthologous genes located outside synteny blocks. to validate this study, we first used the phylogenetic approach, taking into account orthologs contained in trees of various complexities defined as follows. we analyzed trees made up uniquely of orthologs, and trees containing only orthologs and in-paralogs. we then did similar computations for the genomes with the entire dataset obtained using the rsd ortholog detection method. fig. shows that in both rsd and phylogeny approaches, the overwhelming majority of the statistical tests reject the null hypothesis h in favour of h we found that the remaining, untested cases correspond to a few comparisons that are too small to be safely used in this bootstrap test. this is the case, for example, of all comparisons involving nanoarchaea equitans, an archaeon with a very small genome that is also very distant from all other species, including the other archaea. we conclude that there is a universal trend in the evolution of prokaryotic genomes: genes located in regions with preserved gene order evolve more slowly than genes found outside. discussion as soon as several prokaryotic genome sequences were completed, it became clear that very few of their genes had retained their local ordering through evolution. this observation was reinforced as more and more genomes belonging to a wide spectrum of organisms – from very close to very distantly-related – were published. this led to the concept of gene context. as we prepare to be inundated by thousands of completely sequenced genomes in the next few years, it becomes urgent to design a strategy for quickly evaluating the conservation of gene order in newly-published genomes, comparing them to all those already analyzed. this paper describes a suite of programs that automate all the essential steps for evaluating gene order conservation in prokaryotes, as well as presenting some results from their initial application. after identifying in step the whole set of homologs shared by prokaryotic species, we designed two independent strategies to extract from this set the collection of orthologs that are shared by at least two genomes. in a third step, we located the positional orthologous genes in this collection, and we studied the specific properties that distinguish them from other orthologs. correctly identifying positional orthologs this is a crucial step, because mixing paralogs and orthologs can lead to errors in predicting ortholog pairs that seriously compromise studies of the mode of evolution of genes and genomes. this is especially true when studying prokaryotes. we already showed that many gene families stem from ancestral gene duplication events that occurred well before the divergence of present-day species. therefore, orthologs are often more similar between closely related species, allowing them to be distinguished from paralogs. however, when we compare genomes belonging to distant species, the difference in percent identity between orthologs and paralogs falls rapidly, and the risk of confusing errors quickly increases. although there is a long history of research on identifying orthologs, this is an extremely active research area at the moment. numerous papers appeared within months, each one claiming to present the "best method". a consensus does seem to be emerging about the superiority of approaches based on evolutionary studies over best bidirectional blast hits. anyhow, the fact that so many "best methods" were published in a single year clearly underlines that the field is immature. because no single approach is superior in all cases, our work applies two alternative and complementary techniques, reciprocal smallest distance and phylogenetic approaches. both approaches are based on well-established pam evolutionary distances, computed by a maximum likelihood approach. using these pam distances allows directly determining the ortholog pairs separated by rsd, comparing each homolog with all others by means similar to the method already proposed by wall et al.. in our genomes, this method identified orthologs forming, pairs. we expected that the phylogeny approach would systematically increase these figures; since it should collect not only all pairs of orthologs found by rsd but also so-called in-paralogs the rsd approach misses. in fact, this does happen in many cases. for example, as fig. shows for orthologs of the b.subtilis chemotaxis chec protein, we increase the number from found using the rsd approach to with the phylogeny approach. unfortunately, although this phylogeny strategy was found to be quite effective, it can create unexpected difficulties when comparing many protein sequences that form some heterogeneous and very large "families". in fact, these huge agglomerates contained about the half of all found homologs. we managed to decompose many of these into more valid and much smaller families. however, this decomposition step caused the loss of a significant number of genuine pairs of orthologs. on the other hand, phylogenetic analysis is advantageous when it rejects a significant number of homologous pairs that have been found as orthologs by the rsd approach, as being actual paralogs. this will occur for instance when a deep node is correctly viewed as due to an ancestral gene duplication event. moreover, the apparent excess of pairs of orthologs found by the rsd method looks suspicious; since a significant number of them are discarded when looking at the number of pogs found by this approach. therefore, it seems that phylogeny approach is more accurate to detect bona fide orthologs, excluding pseudoorthologs and including in-paralogs. however, it appears presently hazardous to choose a single "best method" to identify orthologs. our present strategy has been to combine the results of our two alternative methods. uses of the data fall between two possible options extremes: using the union and using the intersection of their results. where the goal is to collect a maximum of putative synteny cases in the maximum number of genomes, it is appropriate to merge the orthologs determined by both methods. alternatively, if the goal is to work with pairs of orthologs that are highly likely to be authentic, it would be best to use only the results that are common to both approaches, but such an intersection will unfortunately remove all in-paralogs. evolutionary constraints on pogs based upon these identified synteny blocks and the computed evolutionary distances separating the proteins encoded by these adjacent genes, it was possible to compare the rates of evolution of the orthologous genes located inside and outside synteny blocks, respectively. as shown in table and figs. and proteins encoded by genes located within synteny blocks are found to evolve at a lower rate than those encoded by genes located outside these blocks. this observed difference appears to be consistent with the concept of genomic context in which it is assumed that neighboring genes are under a strong evolutionary pressure to maintain their adjacency. such a pressure would exist, for example, if the proteins they encode have to interact in order to participate to the same biological process. however, in most cases it is not yet clear what the actual cause of low evolutionary rate is. our demonstration that pogs are more evolutionarily constrained than other orthologs that do not reside in synteny blocks confirms and extends the previous observation made by dandekar et al. that "the degree of sequence conservation is on average substantially higher than that in genes that do not exist as conserved genes pairs". moreover, our demonstration agrees with both the variables and universals that have previously been associated with differences in protein evolutionary rates. we expect that persistently adjacent genes result from several of these universals: they are probably essential since they are conserved in large parts of the taxonomic space. they will be co-transcribed at the same level of expression. in many cases, their products will interact in the cell, requiring a significant level of connectivity between them and strict conservation of residues to permit a sufficient density of structural and functional contacts. we performed a preliminary evaluation of these hypotheses using experimental data obtained from the model organisms e. coli and b. subtilis. of the pogs shared by these two species, are essential e. coli genes, non-essential and unknown. moreover, of the shared e. coli and b. subtilis pogs belong to the e. coli interactome and of these are essential genes. in total, around 23% of the pogs e. coli shares with b. subtilis are both essential and part of its interactome, making them what have been called hub proteins. the different data obtained in this work suggest that blocks of adjacent pogs which display conserved co-occurence in various present-day genomes can be interpreted as ancient evolutionary units that have resisted over the years to the imperious fluidity forces. thus, these evolutionary units could be viewed as the remnants of the skeleton of the gene order found in a common ancestor to all compared present-day microorganisms. we further propose that pogs have a low evolutionary rate because they are constituents of functional modules as previously defined by hartwell et al. in a systems biology perspective. being a part of such modules would require a coordinated evolution between their sequence and expression level to participate in a collective cell function.  maintaining the order of genes whose products interact in the same cell function is clearly beneficial to organisms, for instance to allow co-expression and co-regulation. this requires the creation of clusters that have low recombination rates and that the adjacency of functionally related genes is maintained by epistatic selection. however, as we compare species that are distant in the taxonomy space, there is an inverse relation between taxonomic proximity and the size of the synteny blocks. this evidence agrees with poyatos and hurst contention that gene order is forever in flux. therefore, the exceptional instances where a stable gene order is conserved, even for a few gene pairs, implies strong selective pressure that should be reflected in high levels of co-expression of essential products such as one sees, for instance, in the exemplary case of ribosomal proteins. 
accordingly, we used the data obtained with either methods or their intersection to number the orthologs that are adjacent in for each pair of genomes, the positional orthologous genes, and to further study their properties. this is a complex and difficult problem in the field of comparative genomics but will help to better understand the way prokaryotic genomes are evolving. once all these synteny blocks have been detected, we showed that pogs are subject to more evolutionary constraints than orthologs outside synteny groups, whichever the taxonomic distance separating the compared organisms. we have developed a suite of programs that automate three essential steps to study conservation of gene order, and validated them with a set of bacteria and archaea that cover the majority of the prokaryotic taxonomic space. the suite of programs described in this paper allows a reliable detection of orthologs and is useful for evaluating gene order conservation in prokaryotes whichever their taxonomic distance. the first computes the reciprocal smallest distance using the pam distances separating pairs of homologs. detecting synteny blocks requires reliable methods to determining the orthologs among the whole set of homologs detected by exhaustive comparisons between each pair of completely sequenced genomes. we applied two strategies to extract from the set of homologs a collection of valid orthologs shared by at least two genomes. the second method groups homologs in families and reconstructs each family's evolutionary tree, distinguishing bona fide orthologs as well as paralogs created after the last speciation event. we identified the whole set of shared homologs between two or more species and computed the evolutionary distance separating each pair of homologs.
 comparison of completely sequenced microbial genomes has revealed how fluid these genomes are. detecting synteny blocks requires reliable methods to determining the orthologs among the whole set of homologs detected by exhaustive comparisons between each pair of completely sequenced genomes. this is a complex and difficult problem in the field of comparative genomics but will help to better understand the way prokaryotic genomes are evolving.  we have developed a suite of programs that automate three essential steps to study conservation of gene order, and validated them with a set of bacteria and archaea that cover the majority of the prokaryotic taxonomic space. we identified the whole set of shared homologs between two or more species and computed the evolutionary distance separating each pair of homologs. we applied two strategies to extract from the set of homologs a collection of valid orthologs shared by at least two genomes. the first computes the reciprocal smallest distance using the pam distances separating pairs of homologs. the second method groups homologs in families and reconstructs each family's evolutionary tree, distinguishing bona fide orthologs as well as paralogs created after the last speciation event. although the phylogenetic tree method often succeeds where rsd fails, the reverse could occasionally be true. accordingly, we used the data obtained with either methods or their intersection to number the orthologs that are adjacent in for each pair of genomes, the positional orthologous genes, and to further study their properties. once all these synteny blocks have been detected, we showed that pogs are subject to more evolutionary constraints than orthologs outside synteny groups, whichever the taxonomic distance separating the compared organisms.  the suite of programs described in this paper allows a reliable detection of orthologs and is useful for evaluating gene order conservation in prokaryotes whichever their taxonomic distance. thus, our approach will make easy the rapid identification of pogs in the next few years as we are expecting to be inundated with thousands of completely sequenced microbial genomes. 