Multiple Genome Comparison within a Bacterial Species Reveals a Unit of Evolution Spanning Two Adjacent Genes in a Tandem Paralog Cluster

It has been assumed that an open reading frame (ORF) represents a unit of gene evolution as well as a unit of gene expression and function. In the present work, we report a case in which a unit comprising the 3′ region of an ORF linked to a downstream intergenic region that is in turn linked to the 5′ region of a downstream ORF has been conserved, and has served as the unit of gene evolution. The genes are tandem paralogous genes from the bacterium Staphylococcus aureus, for which more than ten entire genomes have been sequenced. We compared these multiple genome sequences at a locus for the lpl (lipoprotein-like) cluster (encoding lipoprotein homologs presumably related to their host interaction) in the genomic island termed νSaα. A highly conserved nucleotide sequence found within every lpl ORF is likely to provide a site for homologous recombination. Comparison of phylogenies of the 5′-variable region and the 3′-variable region within the same ORF revealed significant incongruence. In contrast, pairs of the 3′-variable region of an ORF and the 5′-variable region of the next downstream ORF gave more congruent phylogenies, with distinct groups of conserved pairs. The intergenic region seemed to have coevolved with the flanking variable regions. Multiple recombination events at the central conserved region appear to have caused various types of rearrangements among strains, shuffling the two variable regions in one ORF, but maintaining a conserved unit comprising the 3′-variable region, the intergenic region, and the 5′-variable region spanning adjacent ORFs. This result has strong impact on our understanding of gene evolution because most gene lineages underwent tandem duplication and then diversified. This work also illustrates the use of multiple genome sequences for high-resolution evolutionary analysis within the same species.


Introduction
Gene duplication has long been recognized as an important mechanism in the evolution of genes. Since the essential role of gene duplication in the emergence of novel genes was proposed (Ohno 1970), numerous works have studied the mechanisms by which the duplicated genes diversified (Li 1997). Recent progress in genome sequencing has confirmed the significance of gene duplication by revealing a wide prevalence of multiple homologous genes within a genome, often referred to simply as paralogs (Lynch and Conery 2000;Friedman and Hughes 2001;Gevers et al. 2004). Comparison of multiple genome sequences has provided ample evidence that the present paralogs have undergone various patterns of diversification (Prince and Pickett 2002;Taylor and Raes 2004).
In bacteria, many paralogous gene groups have been shown to be involved in genome rearrangements that help the bacteria adapt to ever-changing environments. These processes are referred to as phase variation or antigenic variation, and the dynamics of paralog evolution in bacteria have been mostly examined with reference to these mechanisms (van der Woude and Baumler 2004; Villemur and Deziel 2005). Studies have revealed that these genes have become rearranged through simple processes such as inversion, deletion, or gene conversion, via various molecular mechanisms such as site-specific recombination, homologous recombination, or slipped-strand mispairing during DNA replication (Hughes and Norstrom 2005;Villemur and Deziel 2005). Many paralogous genes that have undergone these rearrangements have been discovered through whole-genome sequencing. Comparison of closely related prokaryotic genomes can help in the elucidation of the molecular mechanisms underlying the rearrangements of paralogous genes.
Paralogous genes are often present as a tandem cluster in prokaryotic genomes (Gevers et al. 2004;Reams and Neidle 2004). They often encode surface proteins (Kihara and Kanehisa 2000;Gevers et al. 2004), which may undergo phase variation or antigenic variation (Hughes and Norstrom 2005;Villemur and Deziel 2005). Because their genes were likely to have originated from duplication, tandem paralogs can be suitable targets for the study of paralog diversification. Some of these paralogous gene clusters are on genomic islands, which are likely to have been horizontally acquired, to be highly polymorphic among strains, and to confer strain-specific adaptive properties such as drug resistance or pathogenicity (Dobrindt et al. 2004).
Considering the above point, tandem paralogs found in the genomic islands of Staphylococcus aureus are suitable targets for the study of paralog evolution. Staphylococcus aureus is a Gram-positive bacterium with low GC content. It is a major human pathogen and is notable for its expression of a vast variety of toxins (Kuroda et al. 2001). Whole-genome sequences have been determined for more than 10 S. aureus strains (Kuroda et al. 2001;Baba et al. 2002Baba et al. , 2008Holden et al. 2004;Ohta et al. 2004;Gill et al. 2005;Diep et al. 2006;Highlander et al. 2007; http://genome.jgi-psf.org/finished_microbes/). Two genomic islands, mSaa and mSab, common and unique to this bacterial species, encode three tandem paralog clusters: exotoxins (ssl) and lipoproteins (lipoprotein-like [lpl]) encoded by mSaa and proteases (spl) encoded by mSab. All these paralogs encode secreted proteins that are inferred to be pathogenicity related (Williams et al. 2000;Kuroda et al. 2001;Reed et al. 2001;Chavakis et al. 2007;Kulig et al. 2007). The copy number and sequence composition varies among strains. In our previous study (Tsuru et al. 2006), we compared these three clusters from seven strains available at the time and inferred an involvement of homologous recombination in their evolution. Among the three, the lpl cluster was unique in that highly conserved nucleotide sequences, representing a possible recombination site, were discovered at the middle of the proteincoding region of every paralogous open reading frame (ORF) and in that rearrangements there seemed to be far more extensive than in the other two clusters.
In the present study, we made a detailed comparison of the lpl tandem paralog clusters of nine sequenced S. aureus genomes to examine their evolutionary processes. Contrary to the general belief that a unit of evolution of a gene is an ORF itself, we came to the conclusion that a unit composed of the ''3# half of an ORF and the 5# half of a downstream ORF'' serves as a unit of evolution in this cluster.

Materials and Methods
Homology Search for lpl Genes The annotated lpl amino acid sequences of S. aureus N315 were used as the query sequences to search for homologs of lpl genes using BlastP and TBlastN. The search was carried out against the complete genome sequences of 520 bacterial and 46 archaeal species (31 July 2007 data) from the National Center for Biotechnology Information (NCBI) Genome database (http://www.ncbi.nlm.nih.gov/ Genomes/). All hits with an e value , 10 À5 were collected. Pseudogenes of lpl were identified through these analyses, some of which had been referred as such in the original annotations. Against the lpl homologs, manual refinement at the nucleotide sequence level was carried out by macroscopic pairwise genome comparison using CGAT software ) and by multiple sequence alignment using ClustalW version 1.83 (http://www.ebi.ac.uk/Tools/ clustalw/). The initiation and termination positions for ORFs (including pseudogenes) were reassigned consistently to resolve discrepancy between the original annotators. The resulting coordinates for all the lpl ORFs are listed in supplementary table S1 (Supplementary Material online).

Nomenclature of Genes
Names for ORFs and pseudogenes used throughout this study were set by modifying the corresponding locus_tag in RefSeq database (http://www.ncbi.nlm.nih.gov/RefSeq/): a three-letter genome name used in the Kyoto Encyclopedia of Genes and Genomes database (http://www.genome.jp/ kegg/) was followed by the locus_tag number, for example, SAA0410 for SAUSA300_0410 in RefSeq. The strains in which the relevant homologs were identified and their threeletter genome names are as follows. In S. aureus, SAU is for strain N315 (Kuroda et al. 2001), SAV for Mu50 (Kuroda et al. 2001;Ohta et al. 2004), SAM for MW2 (Baba et al. 2002), SAR for MRSA252 (Holden et al. 2004), SAS for MSSA476 (Holden et al. 2004), SAC for COL (Gill et al. 2005), SAO for NCTC8325 (http://microgen.ouhsc.edu/ s_aureus/s_aureus_home.htm), SAA for USA300 (Diep et al. 2006), and SAB for RF122 (Herron-Olson et al. 2007). In Staphylococcus epidermidis, SEP is for strain ATCC12228 (Zhang et al. 2003) and SER for RP62A (Gill et al. 2005). In Staphylococcus haemolyticus, SHA is for strain JCSC1435 (Takeuchi et al. 2005). A pseudogene was indicated by adding ''p'' at the end of the name, for example, SAA0412p.

Sequence Comparison
An initial multiple sequence alignment of the nucleotide sequences for all the lpl ORFs was constructed using Clus-talW with the default parameters. Then a Neighbor-Joining (NJ) phylogeny for the alignment was constructed using MEGA 4.0 (Tamura et al. 2007; http://www.megasoftware. net/) with pairwise deletion mode for gaps and with a maximum composite likelihood model for substitutions ( More detailed analyses were carried out with lpl ORFs on genomic island mSaa of S. aureus as detailed below. For the nucleotide sequences of the ORFs, multiple sequence alignment was once again constructed using ClustalW. A multiple sequence alignment for the predicted amino acid sequences was also constructed, omitting those for the pseudogenes. The nucleotide sequence alignment was then manually refined considering their encoding amino acid sequences, using a function implemented in MEGA 4.0. The resulting nucleotide and amino acid sequence alignments are shown in supplementary figure S3A and B (Supplementary Material online), respectively. The nucleotide sequences of the relevant regions including the ORFs and the intergenic regions are identical between N315 and Mu50 and between MW2 and MSSA476; therefore, the sequences of Mu50 and MSSA476 were omitted here and in the following analyses.

Definitions Related to the Structure of lpl ORFs on mSaa
The presence of the central conserved region was visualized by similarity plots of aligned nucleotide/amino acid sequences constructed using PLOTCON (http://emboss. sourceforge.net/), in which a similarity score, with window size of 5, was calculated with EDNAFULL score file for nucleotide sequences and with EBLOSUM62 for amino acid sequences ( fig. 2). A region conserved both at the nucleotide sequence level and at the amino acid sequence level was determined in the alignments by visual inspection and defined as the central conserved region (supplementary fig. S3A and B, Supplementary Material online). Divergent regions to its 5# side and to its 3# side were defined as the 5#-variable region and the 3#-variable region, respectively.

Phylogenetic Comparison and Grouping of 5#-Variable Region and 3#-Variable Region
Nucleotide sequences of the 5#-variable region and the 3#-variable region were once again aligned using ClustalW, and NJ phylogenies for them were constructed using MEGA 4.0 with the complete deletion mode for gaps and with the maximum composite likelihood model for substitutions. The phylogenies were compared with each other by connecting operational taxonomic units (OTUs) in a pair of the 3#-variable region of an ORF and the 5#-variable region of its downstream ORF ( fig. 4A) and in a pair within an ORF ( fig. 4B). Shapes of the trees were modified by flipping and rerooting to reduce the number of crosses of the connecting lines in each comparison with aid of TreeMap 2.0 (http://www.it.usyd.edu.au/~mcharles/software/treemap/ treemap.html). Grouping was made for each variable region based on the comparison of phylogenetic trees in figure 4A. Pairwise identities with respect to the above grouping and with all the sequences were calculated using MEGA 4.0 (table 1).
For statistically testing whether our tree-based grouping reflects significant linkage between a 3#-variable region and its downstream 5#-variable region, we performed Fisher's exact test implemented in R packages version 2.7.0 (http://www.R-projects.org/), which is based on the FOR-TRAN program FEXACT (Mehta and Patel 1986;Clarkson et al. 1993).

Comparison of Intergenic Regions
Sequences of the intergenic regions were compared using ClustalW for each group. The presence of intergroup similarities was detected with aid of NJ trees using MEGA4.0 and multiple sequence alignments using ClustalW ( fig. 5). Pairwise identities were calculated using MEGA4.0 (table 1). Ribosome-binding sites were predicted in the multiple alignments, referring to the characterized consensus sequences (Novick 1991).

Intergenomic Comparison
Intergenomic comparison was carried out using a map of ORFs shown in figure 6A. For indels, boundaries of homology and nonhomology were compared by ClustalW as described previously (Tsuru et al. 2006) to identify recombination sites ( fig. 6B and C).

Analyses of Another Tandem Gene Cluster for Hypothetical Proteins in S. aureus
The above definitions related to gene structure and phylogenetic comparison were repeated for another tandem gene cluster for hypothetical proteins, homologs of SA1317 (figs. 8-10). Grouping was carried out for each of the two variable regions separately so as the mutual evolutionary distance remains equal to or shorter than 0.15 within a group ( fig. 9).

Distribution of lpl Homologues in Staphylococcus Genomes
In a previous study, we examined diversity of the lpl gene cluster on the genomic island mSaa of several S. aureus strains with a sequenced genome (Tsuru et al. 2006). It had been reported that the homologs of these lpl genes were found in other loci than the genomic island and that they compose the largest group of paralogous genes in S. aureus strain N315 (Kuroda    Baba et al. 2004). Therefore, we first carried out a search for their homologs against the microbial genome sequence database (NCBI Genome database; 31 July 2007 data, see Materials and Methods). This revealed that significantly related lpl homologs were found only in the genomes of the genus Staphylococcus. All the nine sequenced strains of S. aureus examined (N315, Mu50, MW2, MRSA252, MSSA476, COL, NCTC8325, USA300, and RF122) and two sequenced strains of S. epidermidis (ATCC12228 and RP62A) carry multiple lpl homologs within their genome, and the sequenced strain of S. haemolyticus (JCSC1435) carries one lpl homolog (supplementary table S1, Supplementary Material online). In the only remaining Staphylococcus genome that has been sequenced, Staphylococcus saprophyticus strain ATCC15305 (Kuroda et al. 2005), and in the other bacterial or archaeal genomes, no homolog could be found.
In all the S. aureus genomes, paralogous lpl genes were found in four loci ( fig. 1A). First, in the genomic island mSaa, which we analyzed previously (Tsuru et al. 2006), the tandem cluster with three to ten lpl homologs is present in all nine strains (Locus 0). Second, one to five tandemly repeated lpl genes are found in a locus corresponding to SAU0092-SAU0096 of strain N315 in all the strains (Locus I). Third, only one lpl gene is present in a locus cor-responding to SAU0203 of N315 in all the strains except for RF122 and MRSA252 (Locus II). Fourth, one to four tandemly repeated lpl genes, sometimes with intervening ORFs, can be found in a locus corresponding to SAU2269-SAU2275 of N315 in all the strains (Locus III). Note that all the nine S. aureus strains show synteny along the entire genome and carry the lpl genes at the same loci in the same orientation.
In order to examine orthologous/paralogous relationships among the lpl homologs, all the nucleotide sequences were compared in a multiple alignment to construct an NJ phylogeny. The condensed version of the resulting tree is displayed in figure 1B and the uncondensed version in supplementary figure S1A (Supplementary Material online). The tree revealed the presence of three monophyletic groups completely corresponding to Locus I, Locus II, and Locus III of S. aureus. Within each group, sequences are relatively similar to each other. Homologs from S. epidermidis and S. haemolyticus were rather divergent and were apparently related to the family at Locus III of S. aureus. The genes on mSaa of S. aureus (Locus 0) did not form a closely related monophyletic group, rather a mixture of several phylogenetic groups. However, these are clearly distinct from the groups at the other three loci of S. aureus. These observations suggested that the lpl genes of S. aureus have evolved, to the first approximation, separately at each locus (see also Discussion). Therefore, in the following sections, we will focus on sequence comparison among lpl homologs from the mSaa genomic island of S. aureus.

Central Conserved Region in lpl Genes on mSaa
In our previous study (Tsuru et al. 2006), we identified a conserved sequence within the lpl genes on mSaa through a multiple dot-plot analysis using the genome sequences of seven strains available at the time. We here confirmed the presence of this conserved sequence in two additional strains, USA300 and RF122, through the same analysis (supplementary fig. S2, Supplementary Material online). The presence of the conserved sequence is seen as dots at the crossing points of a horizontal red line and a vertical red line in comparison within the same genome, as noted in our previous work (Tsuru et al. 2006).
To verify this conserved sequence, multiple alignments for nucleotide sequences and for predicted amino acid sequences were constructed for all the relevant lpl ORFs (supplementary fig. S3A and B, Supplementary Material online). These alignments revealed the presence of a region conserved at both the nucleotide sequence level and the amino acid sequence level, which we named the central conserved region ( fig. 2). This region is 132 nt long or 44 amino acids long without any gap in the alignments. The regions 5# and 3# of the central conserved region are less conserved and have variable lengths, and, accordingly, they were defined as the 5#-variable region and the 3#variable region, respectively.
Calculation of pairwise nucleotide sequence identities using all the relevant sequences (table 1, All) revealed that those for the central conserved region are high (minimum of 79% and an average of 88%). Out of the 1,128 pairwise The central conserved region is highlighted by gray shading. A predicted signal peptide region is indicated together with a conserved cysteine residue at the C-terminus.
relationships, 427 showed identity as high as 90%. Meanwhile, those for the 5#-variable region and the 3#-variable region show minimal values of 41% and 52%, respectively, whereas the averages are 65% and 66%, respectively.
Involvement of this conserved sequence in the rearrangement of the lpl cluster is visible in an intergenomic dot-plot comparison (supplementary fig. S2, Supplementary Material online) as the termination of long black lines at a crossover point of a horizontal red line and a vertical red line as noted earlier (Tsuru et al. 2006). Generally, homologous recombination requires two homologous sequences long enough and similar enough to each other. In Bacillus subtilis, the closest bacterium to S. aureus in which this process has been studied in detail, the minimal length is 70 bp (Khasanov et al. 1992), a size comparable to those reported for other prokaryotic systems (Shen and Huang 1986;Fujitani et al. 1995). Frequency of homologous recombination is very sensitive to homology length around this length: it was found to be proportional to the third power of the homology length (Fujitani et al. 1995). The frequency of homologous recombination decreases very rapidly as the two sequences diverge (Vulic et al. 1997;Majewski and Cohan 1998;Fujitani and Kobayashi 1999). Thus, the homologous recombination between the central conserved regions is likely to occur.

Phylogeny Comparison and Grouping of 5#-Variable Region and 3#-Variable Region
If the central conserved region served as a recombination site during diversification of this region, the crossingover events there should have changed the combination of the 5#-variable region and the 3#-variable region of an ORF ( fig. 3). These events would result in incongruent phylogenies between them within an ORF. On the other hand, these crossing-over events will not disturb the linkage between 3#-variable region of an ORF and 5#-variable region of the next downstream ORF ( fig. 3). They would result in a congruent phylogeny between these combinations. To test these predictions, a phylogeny for 5#-variable region and that for 3#-variable region were constructed and compared ( fig. 4).
In figure 4B for a pair of variable regions within an ORF, the two phylogenies were found to be significantly incongruent. In figure 4A for a pair encompassing adjacent ORFs, the phylogenies appeared to be more congruent to each other. We observed that the connecting lines in figure 4A could be divided into groups of parallel lines with only a few exceptions. By contrast, in figure 4B, such clear groupings of the parallel connecting lines were difficult to identify (see Materials and Methods).
We also observed that the grouping of the parallel lines in figure 4A coincided with monophyletic or paraphyletic grouping of each phylogeny. Based upon this observation, we were able to assign groups for the 5#-variable region and the 3#-variable region, respectively. The 3#-variable regions were grouped into 13 distinct monophyletic groups named ''A'' through ''M,'' whereas the 5#-variable regions were also grouped into 13 distinct monophyletic groups named ''a'' through ''m.'' The group ''B'' was further divided into a paraphyletic group ''B1'' and a monophyletic group ''B2'' because these two are paired with distinct groups of the 5#-variable region in figure 4A: group B1 is paired with group ''e'' and group B2 is paired with ''h,'' respectively. The internal branch lengths for the resulting groups were relatively short, which indicated sequences are similar to each other within each group. Indeed, pairwise nucleotide identities calculated for each of the resulting groups were calculated to be !80%, which is in contrast to those calculated using all the sequences (table 1).
In all, 12 distinct groups of the pairs could be assigned in the combination of the 3#-variable region of an ORF and the 5#-variable region of its downstream ORF, which are displayed in different colors of connecting lines in figure 4A. The corresponding monophyletic or paraphyletic groups in each phylogeny are displayed in the same, but faint, colors as the connecting lines. One pair, ''L'' and ''l,'' forms a group of only one member. The other groups have multiple members. The sequences of two variable regions are conserved within each group but are diverged between groups. There are only two exceptional pairs in this grouping, which are indicated by gray connecting lines in figure 4A. These are discussed in detail later (see Discussion). Our grouping is based on comparison of the two phylogenetic trees in figure 4A. In order to examine whether it represents a statistically significant relationship, we performed Fisher's exact test (Agresti 1992; see Materials and Methods). The probability calculated for our data is 2.2 Â 10 À16 , which indicates that there is a highly significant linkage between a 3#-variable region and its downstream 5#-variable region.
This presence/absence of the conserved pairs in figure  4A and B supports the hypothesis that the two variable regions of one ORF have been shuffled during the diversification processes via crossing-over events at the central conserved region, whereas the two variable regions encompassing adjacent ORFs have been conserved during these processes.

Comparison of Intergenic Regions
If the linkage between the 3#-variable region of an ORF and the 5#-variable region of its downstream ORF FIG. 3.-An elementary process of diversification through homologous recombination between the central conserved regions. A crossingover will change combinations of the 5#-variable region and the 3#variable region of an ORF, but linkage of 3#-variable region of an ORF, its downstream intergenic region, and 5#-variable region of its downstream ORF will be maintained. has been maintained, their intervening sequence should also have been conserved ( fig. 3). To verify this, the intergenic regions were compared by multiple alignments (fig. 5). They turned out to be highly conserved within each group (see also table 1, intergenic region). Here ''A-a'' represents, for example, an intergenic region sandwiched by 3#variable region of ''A'' and 5#-variable region of ''a''.
In most of the cases, the sequences from different groups differed in length and composition, showing almost no similarity. A putative ribosome-binding site (Novick 1991) could be identified in all their sequences except for ''K-d,'' and their sequences from the different groups are likewise distinctive. There are, however, several cases in which the intergenic sequences from different groups are similar to each other: among ''B1-e,'' '' B2-h,'' and ''D-f''; between ''F-b'' and ''G-k''; and between ''H-j'' and ''L-l'' (fig. 5). Additionally, ''C-k'' and ''C-g,'' the exceptional pairs of the conservation, also showed similarity to each other and to ''A-a''.
Therefore, the conserved pair of the two variable regions can be extended to the conserved unit comprising a 3#-variable region, a downstream intergenic region, and a downstream 5#-variable region, which spans two adjacent ORFs. The units from different groups are substantially distinctive from each other, though some sequence families of intergenic regions are found to be common to a few different groups.

Intergenomic Comparison and Reconstruction of Genome Rearrangement Events
The result of the grouping based on the phylogeny comparison in figure 4A is displayed in a schematic map of ORFs in figure 6A. Gene orders and gene compositions are highly variable among strains, which indicates the occurrence of multiple rearrangement events in this region in the past. Combinations of the 5#-variable region and the 3#variable region within an ORF appear to have been extensively shuffled. For instance, ''f'' is paired with ''G'' in SAU0397/SAV0434 of N315/Mu50, but it is paired with ''C'' in SAA0411 of USA300. At the level of ORF, 32 distinct patterns were generated in the combination of each 13 groups of 3# and 5# regions ( fig. 6A). In contrast, the presence of the conserved pairs of the 3#-variable region of one ORF and the 5#-variable region of its downstream ORF, as displayed by the same coloring, can be clearly observed here. These observations also support the hypothesis that the linkage encompassing two adjacent genes has been well maintained in spite of the extensive rearrangements.
Among these rearrangements, two indels, one found between USA300 and COL and the other between COL and NCTC8325, are remarkable in that the variation between the strains can be explained by one deletion event. The close relationship between these three strains has been elucidated by phylogenetic analysis using the sequences of housekeeping genes (Baba et al. 2008).
For the above two indels, alignment of the sequences involved allowed determination of the recombination site ( fig. 6B and C). In the upstream region, a progeny sequence in the second line of the three perfectly aligned with one of the parental sequences in the first line but not with the other parental sequence in the third line. Then in the downstream region, the progeny sequence in the second line aligns perfectly with the parental sequence in the third line but not with the parental sequence in the first line. The transient region, around which recombination is likely to have taken place, coincides with the central conserved regions in both the cases. This indicates that these indels have been generated by a homologous recombination event at the central conserved region.

A Model of Paralog Cluster Diversification through Homologous Recombination between Central Conserved Regions
Using the observations and arguments presented above, we propose a model for the diversification processes in this tandem paralog cluster ( fig. 3 and fig. 7). During the sequence diversification process, the central conserved region somehow maintained its sequence being sandwiched by two variable regions. The central conserved region has provided a site for mutual homologous recombination. A recombination event was able to change combination of the 5#-variable region and the 3#-variable region within one ORF but unable to disturb the linkage of the 3#-variable region of one ORF, the intergenic region, and the 5#variable region of its downstream ORF, resulting in maintenance of this conserved unit ( fig. 3). These processes make it possible to generate an ORF with a novel combination of the two variable regions.
Multiple rounds of such crossing-over can cause various types of rearrangements including deletions, conversions, translocations, and substitutions, which are observed among the genomes studied ( fig. 6A). The occurrence of translocation can be explained by circle formation and ensuring reintegration, which is proposed as a mechanism to cause gene amplification or phase variation (Mahan and Roth 1989; Howell-Adams and Seifert 2000; Barten and Meyer 2001), though multiple rounds of unequal crossingover events can also result in translocation. An apparent gene conversion can be explained by repeated rounds of crossingover (Yamamoto et al. 1988(Yamamoto et al. , 1992. If horizontal gene transfer takes place between closely related strains, intermolecular recombination might also occur (Hacker and Kaper 2000;Ochman et al. 2000) and could contribute to these rearrangements, as illustrated in figure 7B. Genetic exchange via horizontal gene transfer has been suggested by sequence comparison of the hsdS gene of the Type I restriction-modification Paralog Evolution 2463 FIG. 5.-Multiple alignments of the lpl intergenic regions. ''A-a'' represents, for example, an intergenic region sandwiched by the 3#-variable region of ''A'' group and the 5#-variable region of ''a'' group. A putative ribosome-binding site (Novick 1991), which could be found in all except for ''K-d'', is indicated by asterisks. system, which is linked to the lpl cluster (Tsuru et al. 2006). This genomic island, however, has been considered to be no longer mobile because it harbors only the remnants of an integrase homologue (Baba et al. 2004). Another type of genomic island of S. aureus, called SaPIs, has intercellular mobility with the aid of a specific helper phage (Ruzin et al. 2001;Novick and Subedi 2007), but the helper phage for the mSaa island has not yet been identified. Natural transformation has not been reported in S. aureus. Occurrence of conjugation was proposed in order to explain large-scale chromosome replacement, suggested from multilocus sequence typing (Robinson and Enright 2004).
It is likely that the diversity of this region was generated through accumulation of these intragenomic and intergenomic processes. This model can explain both formation of highly distorted gene orders and shuffling of the two variable regions of an ORF.
Search for Other Tandem Paralogous Gene Clusters in S. aureus with the Same Diversification Pattern as the lpl Cluster One might expect that the present model for diversification is general to paralog clusters: any highly conserved sequence within tandem paralogous genes would create a site for recombination and could cause similar diversification to the lpl cluster. In order to explore this possibility, we searched signs of this kind of rearrangement in other tandem gene clusters.
Two other tandem clusters on the genomic islands of S. aureus, the ssl cluster on mSaa and the spl cluster on mSab, were examined in our previous study (Tsuru et al. 2006), which revealed that the possible recombination sites are located at the upstream region through the 5#-terminal region of the genes in both cases. The mechanism and evolutionary consequence of their diversification are thus different from those for the lpl genes.
We searched for other tandem paralogous genes in strain N315 of S. aureus by performing a self-self genome comparison in CGAT , which is based on all-against-all BlastN analysis and a local program to detect linked repetitive sequences. This screening uncovered 15 tandem clusters (table 2), which include the ssl, lpl, and spl clusters. We then carried out multiple sequence alignment within member genes of each cluster to examine whether they carry a highly conserved sequence sandwiched by variable sequences. One tandem cluster, composed of SAU1317, SAU1318, SAU1319, and SAU1321 in N315 (encoding hypothetical proteins) showed such a structure.
We carried out screening of their homologs by BlastN search (e value cutoff of 10 À5 ) against prokaryote complete genome sequences in the NCBI Genome database (31 July 2007 data). This screening revealed multiple (two or four) genes in all the nine S. aureus strains analyzed and a single homologous gene in each of the sequenced S. epidermidis strains and in the S. haemolyticus strain. Some of these genes were annotated as lipoproteins (Holden et al. 2004), although not all the genes encode a cysteine residue around the signal peptide-like sequence, a hallmark of bacterial lipoproteins (Sibbald et al. 2006) (data not shown). These genes are termed SA1317 homologs in this study. The multiple homologous genes in the S. aureus strains are present in tandem, sometimes with intervening ORFs and/or an intervening prophage. Similarity plots using all the sequences of the S. aureus genomes confirmed the presence of the central conserved region that is sandwiched by two variable regions ( fig. 8).
A phylogenetic comparison was carried out to detect a linkage encompassing two adjacent paralogous genes ( fig. 9). The tree for the 5#-variable region is not as well resolved as that for the 3#-variable region, presumably because of its short and variable length. The two phylogenetic trees seem to be more congruent for pairs encompassing adjacent genes ( fig. 9A) than for pairs within the same gene ( fig. 9B). We could not apply clear-cut grouping based on the congruence of the two trees for this cluster. Therefore, we grouped 5#-variable region and 3#-variable region separately based on each tree. The 3#-variable region was grouped into seven groups named ''T'' through ''Z,'' whereas the six 5#-variable region groups were named ''u'' through ''z,'' such that the evolutionary distance remains equal to or shorter than 0.15 within a group ( fig. 9). We then tried to connect these 3# groups and 5# groups. This grouping is also indicated in the schematic map in figure 10. This revealed the presence of sequence

(A) An unequal crossing-over between sister chromosomes ([b]) will cause deletion ([a] to [c]) and duplication ([a] to [d]). An additional unequal crossing-over ([e]) will result in apparent conversion ([a] to [f]). An intra-chromosomal, unequal crossing-over ([g]) can form a circle ([h]) and a deletion ([a] to [i]), and ensuing re-integration of the circle ([j]) will results in apparent translocation ([a] to [k]). The two routes of deletion formation ([a] to [c]; [a] to [i]) can result in apparent substitution ([c] and [i]). (B) Inter-molecular recombination involving horizontal gene transfer also explains formation of various types of rearrangements. Double cross-overs between a donor [l] and a recipient [m] will cause deletion ([m] to [n]; [m] to [o]), resulting in apparent substitution ([n] and [o]). Additional inter-molecular recombination between [n] and [o] can result in an apparent conversion ([m] to [p]). Another round of recombination between [p] and [n] can result in an apparent translocation ([m] to [q]).
conservation in the pairs of two variable regions (figs. 9 and 10). Fisher's exact test supported the significance of the linkage between these pairs; the P value is 0.000104.
Intergenomic comparison showed that the sequences of this cluster are similar to each other between N315 and Mu50; between USA300, COL, and NCTC8325; and between MW2 and MRSA476 ( fig. 10). This relatedness among the strains is also indicated in the phylogenetic trees shown in figure 9. Comparison between USA300, COL, and NCTC8325 showed that the linkage between ''W'' and ''u'' has been conserved regardless of the apparent insertion of prophages in USA300 and NCTC8325. The presence of intervening nonhomologous ORFs and/or an intervening prophage between SA1317 homologs does not seem to distort the linkage spanning two paralogs between the above closely related strains. On the other hand, comparison between the diverged strains did not provide direct evidence of their rearrangements.
Taken together, a significant relationship between the two variable regions encompassing genes was also identified also in the SA1317 homolog tandem cluster. However, we could not judge whether this linkage is formed by the diversification processes proposed in this study (figs. 3 and 7) because the current data sets provided no evidence of rearrangements.

Discussion
The phylogenetic trees for the 3#-variable region of the lpl ORF and the 5#-variable region of its downstream ORF ( fig. 4A) appeared congruent, but the congruence is not complete. This may be because the level of divergence among the groups is too extensive to reconstruct their evolutionary relationships in either or both of the trees (Nei and Kumar 2000). Another (not mutually exclusive) possibility is that the 5#-variable region may be too short to construct a reliable tree (Nei and Kumar 2000).
The conserved units comprising the 3#-variable region, the intergenic region, and the 5#-variable region are maintained within a group in almost all the cases (figs. 4A and 6). There are, however, some exceptions in their one-to-one correspondence. One example is present in the 3#-variable region of ''B''; there are two types, ''B1'' and ''B2'', which are paired with ''e'' and ''h'' groups of the 5#-variable region, respectively. A similar situation is found in the 5#-variable region of ''k''. ''k'' is paired with ''G'' group of the 3#-variable region in most  cases. However, in SAU0400/SAV0440-SAU0401/ SAV0441 of N315/Mu50, ''k'' is paired with ''C''. Such imperfect conservation of the unit suggests the involvement of mechanisms other than the homologous recombination at the central conserved region. Illegitimate recombination that requires only short or no homology may have taken place at their intergenic region. The putative recombination points lie close to the 5# end in both the cases.
The 3#-variable region of ''C'' group is involved in another exceptional pair: ''C'' is paired with ''k'' in SAU0400/ SAV0440-SAU0401/SAV0441 of N315/Mu50 as mentioned above but paired with ''g'' in SAA0411-SAA0412p in USA300 ( fig. 4A, gray connecting lines; fig. 6A). The prototype of the 5#-variable region of SAA0412p, a pseudogene, could be ''k''. Fast accumulation of mutations in the 5#variable region of the ancestor of SAA0412p may have led to its erroneous grouping into ''g''. Sequence similarity in their intergenic regions ( fig. 5) and the common occurrence of ''H-j'' pair in their downstream ( fig. 6) support this hypothesis.
One prediction from our diversification model is that the 5#-variable region at the upstream end of this cluster and the 3#-variable region at the downstream end should be conserved among strains. In the case of the lpl cluster, conservation of ''k'' at the upstream end is consistent with this prediction (fig. 6A). However, the 3#-variable region at the downstream end is not conserved ( fig. 6A), which indicates the involvement of a mechanism other than our model. Illegitimate recombination could be involved in the diversification of the downstream end. Interestingly, such a recombination event that involves a long homology at one joint and a short homology in the other joint has been reported in various bacteria (Kusano et al. 1997;de Vries and Wackernagel 2002;Prudhomme et al. 2002).
The present model attempts to explain the generation of a novel paralog by changing the combination of the two variable regions. It does not address the issue of the origin of the sequence diversity among these groups. The intergroup sequence diversity contrasts with the intragroup sequence conservation. Thus, it seems reasonable to assume that these represent two distinctive processes. The genomic island carrying this lpl cluster is likely to have been acquired by an ancestor of S. aureus because this island is common to this species but has not been found in any other Staphylococcus species (Baba et al. 2004). The diverse repertoire of the lpl genes may have been formed before the island was acquired; yet, we do not know the precise molecular mechanisms generating the diverse repertoire.
Bacterial surface proteins are often more variable in their amino acid sequence than proteins encoded by housekeeping genes, presumably due to diversifying selection, exerted by the host immune system (Caporale 2003). In FIG. 9.-Phylogenetic comparison of the 5# regions and 3# regions of SA1317 homologs. A nucleotide NJ phylogeny for 5#-variable regions and one for 3#-variable regions were compared with each other by connecting OTUs in a pair of 3#-variable region of an ORF and 5#-variable region of its downstream ORF (A) and in a pair within an ORF (B). The bootstrap values (%) were obtained from 1,000 resamplings. Groups of each phylogeny were assigned so as the mutual evolutionary distances remain equal to or shorter than 0.15 within a group.
S. aureus, IgG-binding protein A gene (spa) and seven S. aureus surface protein genes (sas genes) are known to be so and have been used for strain typing (Shopsin et al. 1999;Mazmanian et al. 2001). In order to examine the evolutionary rate of the lpl genes, we chose SAU0404 (in N315), SAA0418 (in USA300), and SAR0444 (in MRSA252), as a possible orthologous gene set in the lpl cluster (figs. 4 and 6), and compared their evolutionary distances with those of the seven concatenated housekeeping genes (arc, aroE, glpF, gmk, pta, tpi, and yqiL) used in multilocus sequence typing (Enright et al. 2000). The distances calculated for the lpl genes are much longer than those calculated for the concatenated sequences of the seven housekeeping genes: the average Kimura's two-parameter distance (Kimura 1980) for the lpl genes is 0.066 compared with 0.0077 for the housekeeping genes. Phylogenies based on these distances are presented in supplementary figure S1B (Supplementary Material online). This result indicates that the lpl genes have a relatively fast diversification rate. This fast sequence diversification, together with shuffling via recombination as proposed in our model, has very likely contributed to variability in this cluster.
The sequence divergence/conservation was detected using the whole sequences of the lpl region in this study. One might expect that a small portion might undergo gene conversion between different groups. To further explore this possibility, we performed detailed pairwise nucleotide sequence comparisons using Blast2 (Tatusova and Madden 1999). These results detected several examples indicating gene conversion (Tsuru T, Kobayashi I, unpublished data). However, contribution of such gene conversion seems so small that it did not affect the clear grouping in figure 4A. Likewise, the presence of significantly related lpl homologs in the loci of S. aureus other than mSaa and in some other Staphylococcus species leaves the possibility of interloci recombination and interspecies recombination. The FIG. 10.-Maps of SA1317 homolog clusters in various Staphylococcus aureus strains. The SA1317 homologs are drawn in bold lines. Naming of their 5#-variable region and 3#-variable region is after the tree-based grouping in figure 9. The larger intervening ORF, SAU1320 and its homologous genes, is observed in all the strains except for MRSA252, whereas the shorter intervening ORF, SAA1377 and its homologs, is observed in strains USA300, COL, and NCTC8325. SAB1350 and SAB1349 in RF122 are truncated genes homologous to SAU1320. Insertion of a prophage into the larger intervening ORF observed in USA300, NCTC8325, and MW2 is indicated by a black triangle. same detailed Blast2 analysis, however, did not detect any sequence homology indicating such interactions; the maximal length of the sequence homology equal or more than 95% identity was only 48 bp. Therefore, it is likely that the lpl genes of S. aureus have evolved separately at each locus, as was suggested by phylogeny analysis using full-length sequences in supplementary figure S1 (Supplementary Material online).
The sequences of the central conserved region defined by us are not always completely identical to each other along the entire length (table 1 and supplementary fig. S3 [Supplementary Material online]). This suggests that the precise recombination point varies in each event within the central conserved region because the different recombination points will explain the generation of sequence diversity of the central conserved region. In this study, we were able to estimate the recombination site for the two indels found among USA300, COL, and NCTC8325 ( fig. 6). In both these cases, the detected recombination points resided within the central conserved region, but their positions were different from each other. If we examined more sequences and could find recombination events between them, we might be able to test the above hypothesis. In addition, comparison of the sequence of very closely related strains will enable us to determine the position of recombination more precisely. Recently, entire genome sequences have become available for 4 more S. aureus strains: JH1 and JH9 (Mwangi et al. 2007), USA300-TCH156 (Highlander et al. 2007), and Newman (Baba et al. 2008). Additionally, ten genomes of USA300 derivatives were examined by comparative genome sequence analysis (Kennedy et al. 2008). These are closely related to each other or to the other sequenced strains. Although we studied the sequences of these strains to obtain more examples of recombination events, we could not obtain any additional insight into the lpl gene cluster: JH1 and JH9 are almost identical to N315/Mu50 with only a few point mutations throughout the relevant region, and USA300-TCH156 and Newman are almost identical to USA300. From this analysis, we concluded that testing the above hypothesis would require many more genome sequences of appropriate divergence.
The frequent recombination at the central conserved region of the lpl genes on mSaa suggests, on the other hand, that an intrastrain homogenization process might occur in this region. However, we could not detect such tendency in the phylogenetic analysis using the sequences of the central conserved region (data not shown).
Why has the central region been conserved? One possibility is that this region has been maintained because of its functional importance. At present, neither a functional role for the central conserved region nor a physiological role of these lpl paralogs has been discovered. Generally, bacterial lipoproteins are involved in cellular processes such as antibiotic resistance, adhesion to host cells or other bacterial cells, transport of substances, or intercellular communication (Sibbald et al. 2006). Several lipoproteins are recognized as an antigen by a host's immune system, and they thus affect bacterial survival and their pathogenicity (Henderson et al. 1996). Antigenic variation systems involving lipoproteins have been identified and characterized in Borrelia (Haake 2000;Norris 2006), Porphyromonas (Hall et al. 2005), and Mycoplasma (Lysnyansky et al. 2001;Glew et al. 2002;Ron et al. 2002). Tandemly encoded lipoproteins such as vsp lipoproteins in Borrelia hermsii (Restrepo et al. 1994; Barbour and Restrepo 2000) or p35 family lipoproteins in Mycoplasma penetrans (Neyrolles, Chambaud, et al. 1999;Neyrolles, Eliane, et al. 1999) cause antigenic variation. We do not know whether the interstrain variation in the lpl gene arrangement represents a type of antigenic variation, but we think it is likely because this explains why many copies of these paralogs have been maintained (Jordan et al. 2001;Hooper and Berg 2003). The relatively high evolutionary rate observed in this cluster also supports this view. If this is the case, the central region may be conserved for its function as a recombination site for programed diversification similarly to DNA repeats in some other antigenic variation systems. Shuffling of two variable regions may contribute to alteration of sets of Lpl proteins on the cell surface, partly because replacement of the signal peptides in the 5#-variable region will affect efficiency of secretion of each lipoprotein (Sibbald et al. 2006). The wide diversity in the sequence of the intergenic region including a putative ribosome-binding site ( fig. 5) might contribute to diversity in expression efficiency of these genes. Recent transcriptome analyses using the sequenced strains have reported expression of some of their lpl genes (http://www.bioinformatics.org/sammd/; Sobral et al. 2007), although this has not been confirmed by proteome analyses (Nandakumar et al. 2005;Gatlin et al. 2006). How their expression is regulated has not been elucidated. Precise analysis to identify their differential expression within a strain and between strains is necessary to understand their biological role.
A prominent feature of our model of the diversification process is the presence of the conserved unit encompassing adjacent ORFs. Mosaic gene formation by intergenomic and intragenomic recombination has been reported in various bacteria for antibiotic resistance genes (Maiden 1998;Normark BH and Normark S 2002) and antigen genes (Bessen and Hollingshead 1994;Hobbs et al. 1998;Lachenauer et al. 2000). For allelic variation in the por gene for porin protein in Neisseria gonorrhoeae, homologous recombination of a partial conserved sequence within the gene is proposed to be responsible for diversity (Cooke et al. 1998;Fudyk et al. 1999;Hamilton and Dillard 2006). However, a linkage of segments beyond adjacent ORFs has not been analyzed so far. The present characteristic mode of tandem paralog diversification, maintaining the 3# part of a gene, the intergenic region, and the 5# part of its downstream gene as a unit of evolution, is, to our knowledge, novel among studies of paralog rearrangements.
How general is the presented model of diversification among tandemly repeated genes? A further search for other tandem paralog clusters on S. aureus genome revealed another tandem gene cluster, SA1317 homolog, with a gene structure in which a highly conserved sequence is sandwiched between variable sequences ( fig. 8). A significant linkage between the variable regions encompassing the genes was also identified in this cluster (figs. 9 and 10), although whether this linkage is formed by the present diversification mechanism is not clear because no evidence of rearrangement was found there. We also examined the other well-characterized tandem paralog genes in other bacteria. In the p35 lipoprotein genes in M. penetrans, a conserved region is located at their 5# end encoding a signal peptide (Sasaki et al. 2002). In the vsp genes in Mycoplasma bovis, a conserved sequence is found upstream of the ORF (Lysnyansky et al. 2001). In both these cases, a different diversification process from our model must be operating.

Conclusion
We tend to regard an ORF as a unit of gene evolution as well as a unit of gene expression and its function. The present work with a tandem paralog cluster identified a unit consisting of the 3# half of a gene, a downstream intergenic region, and the 5# half of a downstream gene as the unit of evolution. This is because the central region of the gene provides a site to recombine the 5# half and the 3# half of a gene to generate variation.

Supplementary Material
Supplementary table S1, figures S1-S3, and legends for the supplementary figures are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).