s2fieldsofstudy: Biology Medicine

text:
Membrane testosterone binding sites in prostate carcinoma as a potential new marker and therapeutic target: Study in paraffin tissue sections

Background Steroid action is mediated, in addition to classical intracellular receptors, by recently identified membrane sites, that generate rapid non-genomic effects. We have recently identified a membrane androgen receptor site on prostate carcinoma cells, mediating testosterone rapid effects on the cytoskeleton and secretion within minutes. Methods The aim of this study was to investigate whether membrane androgen receptors are differentially expressed in prostate carcinomas, and their relationship to the tumor grade. We examined the expression of membrane androgen receptors in archival material of 109 prostate carcinomas and 103 benign prostate hyperplasias, using fluorescein-labeled BSA-coupled testosterone. Results We report that membrane androgen receptors are preferentially expressed in prostate carcinomas, and they correlate to their grade using the Gleason's microscopic grading score system. Conclusion We conclude that membrane androgen receptors may represent an index of tumor aggressiveness and possibly specific targets for new therapeutic regimens.

Background
The biological activity of androgen occurs predominantly through binding to an intracellular androgen receptor (iAR) protein, a member of the nuclear receptor superfamily, functioning as a ligand-activated transcription factor [1]. However, in recent years, a number of reports indicate additional rapid androgen actions, including the rapid activation of kinase-signaling cascades, modification of actin cytoskeleton and modulation of intracellular calcium levels [see [2,3], for reviews]. These actions are considered to be non-genomic in nature because they occur in cells lacking functional iARs in the presence of inhibitors of transcription and translation, or are observed too rapidly to involve changes in gene transcription [4].
We have previously identified membrane androgen receptor (mAR) sites in prostate and breast cancer cell lines [5][6][7]. In a preliminary work performed on freshly prepared epithelial cells from prostate carcinoma, non-cancerous peritumoral tissue, and benign prostate hyperplasia (BPH), we have shown that mARs are expressed preferentially on carcinoma cells [8]. Activation of mARs induces actin cytoskeleton polymerization and redistribution [5,9], secretion of Prostate Specific Antigen (PSA) and apoptosis [6,7]. Membrane androgen receptor sites were responsible for the induction of testosterone-BSA-induced apoptosis in T47D breast cancer cells [7]. It is interesting to note that mARs may be different from iARs, as they are not recognized by antibodies against the latter [5], and are not inhibited by a series of commonly used antiandrogens in vitro [9] or in vivo [6].
In the present work we developed a method for the in situ detection of mARs in formalin-fixed and paraffin-embedded tissues. The aim of the present study was to: (1) Validate a method for the detection of mARs in formalin-fixed paraffin-embedded specimens of prostate tumors, and (2) explore the distribution of mARs in prostate carcinomas. More specifically, we have tried to reply to the following questions: (a) Are mARs equally distributed in cases of prostate carcinoma and BPH, and (b) are mARs expressed in prostate carcinoma related to the gravity of the disease, as expressed by the Gleason's score?

Methods
One hundred and nine (109) cases of prostate carcinoma (age 45-93 years, mean ± SD 71.17 ± 9.02, median 71 years, followed for 4 to 72 months) and one hundred and three (103) cases of BPH (age 37-92, mean ± SD 70.44 ± 8.38, median 70 years), retrieved from the archives of the Pathology Department of the University Hospital of Crete, were analyzed retrospectively. Prostate carcinoma specimens corresponded to 23 radical prostatectomies, 15 transurethral resections harboring incidentally encountered carcinoma, and 71 transrectal biopsies. Three cases of carcinoma diagnosed in transurethral resection specimens were rejected for technical reasons. Survival of patients ranged from 1-72 months (mean ± SD 31.4 ± 9.32, median 29.5 months). Benign prostate hyperplasia cases corresponded to 96 transurethral resection and 7 prostatectomy specimens. Hematoxylin and Eosin (H&E) stained slides were reviewed by two investigators independently and blindly to the patients' clinical data. The Gleason's sum (combined score, Gleason's score) [10,11] of carcinoma cases was reevaluated by two observers. The results obtained from the histological study of carcinomas using the Gleason's grading system are presented in Table  1. In case of discrepancy between the two observers, the final decision was reached by consensus.
A set of three serial, 3 µm thick, tissue sections, embedded in Paramat extra (BDH Lab Supplies, Poole, UK), were taken on negatively charged slides (SuperFrost Plus, Kindler O GmbH, Freiburg, Germany) from each representative paraffin block. Two of the slides were used for the needs of mAR detection, while the third was stained with H&E for morphological study, and for grading in the case of carcinomas.
In order to evaluate the presence of mAR in tissue sections, we had to (partially) regenerate membrane proteins. In this aim, in an initial set of experiments, we have assayed different combinations of melting temperatures and times of incubation. Our findings indicated that the optimal conditions for the identification of the highest specific binding on the target molecules (mARs) were achieved by mild melting of the embedding medium (Paramat extra) at 42.5°C for 20 min, followed by dewaxing in xylene (three times, 10 min. each) and rehydration in decreasing concentrations of ethanol. Then, specimens were washed in distilled water for 20 min. and incubated at 37°C in a citrate buffer 0.01 M (pH 6.2) for two hours. Finally, they were washed in Tris Buffered Saline (TBS, 10 mM, NaCl 150 mM, pH 7.4) and processed for the detection of mARs. As in the present study we used albuminconjugated steroids, and in order to minimize the nonspecific binding of BSA we pre-incubated slides with 3% BSA for 40 min. Then the slides were washed in TBS and one of them (test slide) was incubated for 1 hour with 10 -6 M Testosterone-BSA-FITC while the other (control slide) with 10 -6 M BSA-FITC (both reagents were obtained from Sigma, St Louis, MO) in TBS.
Another potential source of non-specific staining could be the association of the ligand (testosterone) with iARs, due to the fact that microscopic tissue sections contain broken (sectioned) cells. Since our previous experiments have shown that classical antiandrogens react with iARs but not with mARs [6,9], we routinely used cyproterone acetate (a specific androgen receptor antagonist, Sigma, St Louis, MO) at a concentration of 10 -4 M, diluted in the testosterone-BSA-FITC or BSA-FITC solution in order to block iARs (see Results). After one-hour incubation, slides were rinsed with TBS, coverslipped using Polyvinyl Alcohol Mounting Medium with Dabco Antifading (Fluka Biochemika) and examined under a fluorescence microscope.
An introductory comparative study using both a confocal laser scanning microscope (Leica TCS SP) and a fluorescence microscope (Nikon Microphot-FXA fluorescence microscope, mounted with a Nikon DX-DB2 photographic camera) showed no significant advantage of the confocal over the conventional fluorescence microscope. The latter was therefore used for the routine detection of mARs. Results were expressed as a percentage of mAR positive cells in reference to total number of cells examined in 10 sequential high power fields (X400) in a representative area of each test slide. The localization and pattern of staining were recorded.
Clinical data, namely age, previous treatment, follow-up and survival, were retrieved from the patients' files. Statistical analysis was performed by the Systat V 10.0 (SPSS, Chicago, IL) program. In the case of positive/negative cases of BPH and prostate carcinoma, we have used significance limits of the four-fold table test, using the χ 2 test, with 1 degree of freedom. The same test was used for the comparison of positive and negative cases in each group. For exploring the correlation between positive/negative cases of mAR expression and other parameters (age, treatment, survival, Gleason's sum) the Spearman rank correlation coefficient was calculated. The Mann-Whitney U test was used to compare the distribution of continuouslyscaled outcomes by mAR expression. All statistical evaluations were double-sided.
This study was approved by the Ethics Committee of the University Hospital of Heraklion, and the Ethics Committee of the University of Crete.

Validation of the method
Previous results indicated that mARs can be detected in freshly prepared prostate carcinoma tissue [8]. The method presented here for the in situ detection of mARs is based on the binding of a ligand (FITC-BSA-bound testosterone) on mARs. Membrane androgen receptors are membrane proteins [5,12], expected to be (partially) denatured in microscopic tissue sections by formalin fixation, paraffin embedding and high temperatures, conditions known to modify integral membrane proteins. Therefore, we attempted to renaturate these proteins by incubating the tissue sections at relatively low temperatures during deparaffinization and renaturation process (mild tissue dewaxing and protein renaturation). Different combinations of melting temperatures (from 40° to 60°C), incubation times, xylene treatment and acidification were tried. As indicated in the Methods section, the best results were obtained with a mild paraffin melting and protein renaturing temperature of 42.5°C, followed by mild acidification in citrate buffer to dissociate any ligands loosely-bound to mARs or iARs [13] and extensive washing with TBS pH 7.4.
When slides were incubated with testosterone-BSA-FITC without any treatment, they presented a high non-specific binding ( Figure 1, upper panel), due probably to the nonspecific absorption of BSA to membrane and/or intracellular structures or molecules. For this reason, slides were preincubated with a 3% BSA solution prior to ligand binding. As shown (Figure 2A and 2B) this treatment resulted in an attenuation of conjugated or unconjugated BSA-FITC to membranes.
Another matter of concern was that iARs, regularly present in prostate carcinoma cells, could bind the ligand since, in tissue sections, one expects an interaction of the conjugated ligand with membrane and (due to sectioning) intracellular structures and/or molecules, including iARs. Our previous results with prostate cancer cell lines have shown that interaction of testosterone-BSA with mARs is not influenced (at the level of binding or the resulting sig-naling cascade and secretion) by antiandrogens [6,9] Therefore, to eliminate a possible interaction of BSA con-jugated ligand with iARs, we used 10 -4 M cyproterone acetate in the binding experiments.
This treatment decreased nuclear, membrane and intracellular staining (Figure 1). The treatment of slides with testosterone-BSA-FITC in the presence of 3% BSA and 10 -4 M cyproterone acetate decreased non-specific staining in nearly background levels, making more discernible the interaction of testosterone-BSA-FITC with specific membrane-binding sites (compare Figure 1 with Figure 2A right panel).
According to the above results, for the detection of mARs in tissue sections, we applied BSA (3%) and cyproterone acetate (10 -4 M) throughout the present study. These conditions were found to accurately discriminate between positive (>75% of cells labeled by testosterone-BSA-FITC) and negative (<10% of cells labeled by testosterone-BSA-FITC) cases.

Prostate carcinoma specimens preferentially express mARs
On examination of the specimens, positive cases of both carcinomas and BPHs were found to present peripheral, membrane type, fluorescence. In contrast, morphologically normal-looking prostate epithelial tissue showed no specific fluorescence. In Figure 2, representative cases of prostate carcinoma are presented, showing a positive (Figure 2A) or negative ( Figure 2B) staining for mARs. Positive cases showed peripheral, continuous and/or stippled, membrane type, fluorescence. This extended over more than 75% of carcinoma cells of the tissue section. Negative cases, on the other hand, showed no more than 10% reacting cells. Intensity of staining in positive cells was not considerably different among cases. Thus, reading and translating mAR positivity seems to be a matter of nearly a "black and white" effect and membrane pattern recognition.
The tissue stroma sometimes showed a diffuse positive reaction, probably due to connective tissue non-specific binding. This non-specific staining did not interfere, though, with the estimation of mAR expression in epithelial cells of the tissue section. Ninety percent (90%) of BPHs were negative for mAR, while about 38% of carcinomas were positive ( Figure 3) (χ 2 = 23.24, p < 0.0001). No correlation was found between mAR expression and age of the patient, treatment before or after surgery and survival.

Membrane Androgen Receptor sites are related to the level of differentiation of prostate carcinoma
As mentioned above, 38% of prostate carcinomas were found to express mARs (χ 2 6.68, p < 0.1, among all cases of carcinoma, n = 109). Thus, mAR expression cannot be used for the discrimination between malignant and benign prostate cases. To evaluate the possible role of mARs in the evolution of prostate carcinoma, we correlated the presence or absence of mARs to the differentiation of the prostate carcinoma, estimating the Gleason's sum of each case on the corresponding H&E stained slide. Using correlation analysis, our results showed that mAR expression is positively related to higher Gleason's sum (Mann-Whitney U = 991, p < 0.04, Spearman's rho = 0.236, p < 0.01, Figure 4). It is worthy to note that 50% of poorly differentiated carcinomas (Gleason's sum 8 and 9, in our series) expressed mARs while in the Gleason's sum 9 carcinomas, mARs expression predominated. As seen in Figure 4, there is an increased expression of mARs with increasing Gleason's sum (15%, 40%, 44%, 36%, 71%, for Gleason's sums from 5 to 9 respectively). No difference was found between Gleason's sum 4+3 and Gleason's sum 3+4 regarding mAR expression (11 out of 21 cases of Gleason's sum 3+4, and 9 out of 20 cases of Gleason's sum 4+3 expressed mAR).
Different percentages of mAR positivity were revealed when prostate carcinoma cases were stratified according to the type of tissue sample. The number and percentages of positive carcinoma cases in each type of specimen is presented in Table 2. Radical prostatectomies showed 13.04% positivity and transurethral resection specimens showed 26.22%, while almost half of the transrectal biopsies were found to express mARs (47.80%). The number of transrectal biopsies (68 studied) was higher correlated to the number of cases of the other two groups (radical prostatectomies 23 and transurethral resections 15). Statistical analysis showed that there is a significant differ-ence in mAR expression between the three types of tissue samples (binomial distribution comparison, p < 0.01).

Discussion
The classical mode of action of steroid hormones includes their entrance into hormone-dependent cell and binding to their intracellular cognitive receptors, which subsequently dimerize and translocate to the nucleus, acting as specific transcription factors [14]. This sequence of events requires time for their completion (>2 hours). In recent years, however, steroids were found to exert a number of effects in seconds or minutes, a phenomenon not explained by their intracellular action. Ion mobilization, activation of specific signaling cascades or secretion are included among these rapid actions and are considered non-genomic, independent of transcription or translation [ [15], see [16], for a review]. A number of possible explanations for the non-genomic steroid actions include the presence of membrane binding sites, different from intracellular receptors [17], the anchorage of intracellular receptors to the plasma membrane through post-translational modification of receptor molecules [18], the interaction of intracellular steroid receptors with other membrane-bound receptors, or the activation of signaling molecules through intracellular receptor-bound steroids.
Recently, a progesterone binding membrane site has been isolated and found to belong to the seven-loop G-protein coupled receptor superfamily [19]. Data on other steroids binding membrane sites have also been reported [2,4,[20][21][22], although their molecular nature has not been elucidated until now. More recently, we have reported the identification of mARs in prostate cancer cell lines [5,6] and in a small series of human prostate carcinomas [8].
The above data indicate that the detection of mARs could be of value in the diagnosis; and possibly mARs themselves could be exploitable in a new therapeutic approach of the mAR positive prostate carcinomas, as prostate carcinoma xenografts in nude mice are significantly reduced in size by testosterone-BSA treatment of animals [6]. In the present study we investigated whether mARs could also be detected in formalin-fixed and paraffin-embedded specimens of prostate carcinoma. The main disadvantage in mARs detection is that, until the elucidation of their molecular structure, no antibodies interacting with them are available. In addition, antibodies to iAR do not crossreact with mAR [7]. Therefore, for the time being, the only appropriate method for mAR detection is that based on Distribution of positive and negative mARs in BPH and pros-tate carcinoma Figure 3 Distribution of positive and negative mARs in BPH and prostate carcinoma. Ninety % of benign prostate hyperplasia cases (n = 103) were negative for mAR receptor expression. On the other hand, prostate carcinoma cases (n = 106) showed mAR expression in a significant percentage (38% of these cases were positive) (χ 2 = 16.7, p < 0.0001).
the use of labeled non permeable ligands, such as testosterone-BSA-FITC. Previous results of experiments of our team performed in formalin-or PFA-fixed cell lines (LNCaP, DU-145, and PC12) and circulating leukemic lymphocytes have shown that mAR could be effectively detected on the cell surface. In view of these data, in the present study we used mild deparaffinization, in combination with citrate buffer pretreatment, known to induce dissociation of ligands or iAR loosely-bound on membrane sites [7,13]. Our results indicate that, under such conditions, mARs could be detected in prostate carcinomas on formalin-fixed and paraffin embedded tissue sections.
A problem encountered in our attempt to reveal in situ mARs was that BSA-FITC bound non-specifically to/or it was absorbed by stromal structures or cell and/or intracellular components. In order to decrease this non-specific interaction, we pre-incubated tissue sections on slides with high concentrations of BSA, thus eliminating or efficiently decreasing this interaction. Intracellular ARs (iARs) do not interfere with the detection of mARs by the method presented here. Indeed, as found in our previous studies, these latter sites do not react with antiandrogens, which can block intracellular receptors [6,9]. We have therefore exploited this differential antiandrogen selectivity in order to block selectively iAR binding, by the use of high concentrations of the antiandrogen cyproterone acetate. In addition, we report that, although confocal microscopy could be considered more precise and accurate than the -cheaper in hardware-conventional fluorescence microscopy, once the researcher becomes accustomed to mAR expression pattern, conventional fluorescence microscopy can be a valuable alternative in revealing the presence of mAR in tissue sections.
Our findings indicate that intensity is not a factor to take into account in the evaluation of mAR expression, since, in individual positive cases, a rather uniform staining intensity was observed. In the case of advanced prostate carcinomas there is a heterogeneity of cell populations expressing iAR [24,25], suggesting iAR counting in hot spots. Such heterogeneity was not observed in our study of mARs, and thus positivity was estimated in 10 sequential high power fields in representative tissue section regions (X400).
The most appropriate samples for the study of mAR expression seem to be the transrectal biopsies (50% positivity), which constituted most of our specimens (71 out of 106 finally studied), followed by the TUR specimens (26.66% positivity), and the radical prostatectomy specimens (13.04% positivity). This signifies that the quality of fixation (as depended on time and room temperature) is of critical importance in the study of mARs.
We report that only 38% of prostate carcinoma cases are positive for mAR staining. In this respect, mAR positivity could not be a discriminant factor for prostate cancer diagnosis. Stratification of cases, according to the Gleason's sum revealed that mAR-positive staining correlates with less differentiated tumors. No significant difference in mAR expression between Gleason's sum 4+3 (9 out of 20 cases) and Gleason's sum 3+4 (11 out of 21 cases) prostate carcinomas was found. Nevertheless, no solid conclusions can be drawn, due to the relatively small number of cases studied.
Correlation of mAR positivity with Gleason's sum in prostate carcinomas Figure 4 Correlation of mAR positivity with Gleason's sum in prostate carcinomas. Membrane androgen receptors are preferentially expressed in prostate carcinomas and their presence is correlated to Gleason's sum.  Membrane androgen receptors are not only expressed in the prostate. They have been found in osteoblasts, lymphoid cells and Sertoli cells [reviewed in [12]] independently of the presence or not of iARs. Preliminary data from our group studies indicate that mARs are expressed in higher grade carcinomas of the breast, B-lymphomas, pheochromocytomas and in AC133 + cord blood stem cells (unpublished observations). In addition, in estrogen receptor negative (ER-negative) carcinomas of the breast we have found a strong staining for mAR [7]. In this respect, mAR might be considered rather as a marker of "immature" and behaviorally aggressive cells than a tissue specific marker. If this comes true, it could explain the finding that tumors with high Gleason's sum preferentially express mARs.

Conclusion
Our results indicate that mARs detection cannot discriminate between prostate carcinoma and BPH, not being exploitable in the differentiation between low-grade carcinomas and BPHs. Nevertheless mARs could be used for a biologic grading of prostate carcinomas, since the more aggressive ones-of higher histopathological grade-preferentially express mARs. Results of recent studies of our team indicate that agonists of mARs-being expressed in cases of both iAR-positive and iAR-negative prostate cancer cell lines-might act as specific therapeutics directing carcinoma cells towards apoptosis [6]. In addition, the activation of mAR-signaling cascade leads to apoptosis [6], while activation of iARs is antiapoptotic [26]. This, in combination with the finding that prostate carcinomas of low Gleason score have a significantly higher iAR content than those with high Gleason score [27], signifies a different biological role for mARs in relation to iARs. Hence, the possibility that combined examination of iARs and mARs in prostate carcinomas, in the context of a test clinically applicable and easy to perform, could lead to a more accurate biologic grading of predictive and/or prognostic importance. To test the prognostic power of mARs, further studies in prostate carcinomas are needed. Moreover, mARs' study in large series of prostate carcinomas is imperative in order to fully elucidate the pathway of their action and reveal whether they represent a possible target for new therapies of these carcinomas, especially in view of our results on regression of prostate tumor xenografts by testosterone-BSA [6].


extfieldsofstudy: Medicine Biology

text:
An emerging phylogenetic core of Archaea: phylogenies of transcription and translation machineries converge following addition of new genome sequences

Background The concept of a genomic core, defined as the set of genes ubiquitous in all genomes of a monophyletic group, has become crucial in comparative and evolutionary genomics. However, it is still a matter of debate whether lateral gene transfers (LGT) may affect the components of genomic cores, preventing their use to retrace species evolution. We have recently reconstructed the phylogeny of Archaea by using two large concatenated datasets of core proteins involved in translation and transcription, respectively. The resulting trees were largely congruent, showing that informational gene components of the archaeal genomic core belonging to two distinct molecular systems contain a coherent signal for archaeal phylogeny. However, some incongruence remained between the two phylogenies. This may be due either to undetected LGT and/or to a lack of sufficient phylogenetic signal in the datasets. Results We present evidence strongly favoring of the latter hypothesis. In fact, we have updated our transcription and translation datasets with five new archaeal genomes for a total of 6384 and 2928 amino acid positions, respectively, and 25 taxa. This increase in taxonomic sampling led to the nearly complete convergence of the transcription-based and translation-based trees on a single phylogenetic pattern for archaeal evolution. In fact, only a single incongruence persisted between the two phylogenies. This concerned Methanopyrus kandleri, whose placement remained strongly biased in the transcription tree due to its above average evolutionary rates, and could not be counterbalanced due to the lack of availability of closely related and/or slower-evolving relatives. Conclusion To our knowledge, this is the first report of evidence that the phylogenetic signal harbored by components of the archaeal translation apparatus is confirmed by additional markers belonging to a second molecular system (i.e. transcription). This rules out the risk of circularity when inferring species evolution by small subunit ribosomal RNA and ribosomal protein sequences, since it has been suggested that concerted LGT may affect these markers. Our results strongly support the existence of a core of proteins that has evolved mainly through vertical inheritance in Archaea, and carries a bona fide phylogenetic signal that can be used to retrace the evolutionary history of this domain. The identification and analysis of additional molecular markers not affected by LGT should continue defining the emerging picture of a genuine phylogenetic core for the third domain of life.

Background
The discovery that Lateral Gene Transfers (LGT) play a major role in the evolution of prokaryotic organisms has raised concerns about the possibility of reconstructing species phylogenies [1]. Some biologists even believe that LGT have obscured the phylogenetic record to such an extent that the task may be in fact hopeless [1,2]. However, others have argued that the careful selection of a 'core' of genes that have been refractory to transfer may help solving (at least partly) this conundrum [3][4][5][6]. The genomic core concept, i.e. the set of homologous genes present in all -or most-genomes of a phylogenetically coherent group, has become crucial in comparative and evolutionary genomics [7]. Indeed, the identification of 'genomic cores' can provide crucial information on the composition of ancestral genomes [8,9], as well as on organisms evolution at various phylogenetic depths [6,10,11]. However, homology-based analyses to define core genes cannot discriminate between vertically transmitted components and horizontally exchanged ones (i.e. "cryptic orthologous replacements" [7]). Thus, it is still a matter of debate whether Lateral Gene Transfers (LGT) may affect the components of genomic cores, preventing their use to retrace species evolution. Nevertheless, the extent and nature of the horizontal component of genomic cores can be identified by molecular phylogeny. The use of conserved gene cores to retrace species evolution has mainly focused on the translation apparatus, since the ribosome appears to be one of the best conserved macromolecular machines in the living world. The concatenation of either bacterial and archaeal ribosomal protein sequences has produced global phylogenies that are roughly similar to those obtained with both small and large ribosomal subunit rRNA genes (16S and 23S rRNA) [4,5,11]. Moreover, careful individual analyses have indicated that ribosomal proteins have been apparently never exchanged between the three Life domains, and rarely between different lineages within domains [4,5,12,13]. However, it may be argued that concerted LGT involving rRNA and ribosomal protein genes, since they belong to the same macromolecular machinery, could escape detection in such analyses. Nevertheless, this hypothesis could be discarded if phylogenies based on additional sets of genomic core proteins belonging to other molecular machineries are congruent with those of the components of the translation apparatus.
Accordingly, we have recently performed an in-depth analysis of proteins involved in transcription and translation from Archaea [11]. Individual phylogenies of these proteins confirmed that the components of these informational molecular machineries are little affected by LGT in the archaeal domain, and permitted the assembly of two large concatenated datasets of likely vertically-transmitted genes to reconstruct the phylogeny of the third domain of life [11]. The trees based on the 'translation' dataset (53 ribosomal proteins, Figure 1A) and the 'transcription' dataset (11 RNA polymerase subunits and 3 transcription factors, Figure 1B) were globally congruent, suggesting that the two informational systems contain a coherent phylogenetic signal for the archaeal phylogeny [11]. However, a number of incongruent nodes remained between the two trees ( Figure 1A and 1B). First, the hyperthermophilic methanogen Methanopyrus kandleri was close to other methanogens in the translation tree ( Figure  1A), whereas it emerged with a strong statistical confidence at the base of the euryarchaeal phylum in the transcription tree (Bootstrap Value BV = 90%, Figure 1B). A second incongruence concerned the position of the euryarchaeon Archaeoglobus fulgidus, since this archaeon was grouped, albeit with weak support (BV = 41%), with Thermoplasmatales in the translation tree ( Figure 1A), whereas in the transcription tree it was strongly placed as sister group to the clade composed of Methanosarcinales and Halobacteriales (BV = 100%) ( Figure 1B). Finally, although in both phylogenies Methanobacteriales and Methanococcales were located in-between Thermococcales and a large clade comprising Thermoplasmatales, Archaeoglobus, Methanosarcinales and Halobacteriales, they were paraphyletic in the translation tree ( Figure 1A) whereas they were monophyletic in the transcription tree ( Figure 1B).
In the case of M. kandleri, we suggested that the discrepancy between the translation and transcription trees was likely due to the very fast evolutionary rate of its RNA polymerase subunits (reflected by the very long branch of M. kandleri in the transcription tree, Figure 1B). Such an accelerated rate of evolution may be due to the lack in this archaeon of the critical transcription factor TFS [11,14]. A Long Branch Attraction (LBA) artefact [15] between the very long branch of M. kandleri and the outgroup (i.e. Crenarchaeota) may thus be responsible for the basal position of this methanogen in the transcription tree. In contrast, the incongruence between the two trees in the position of A. fulgidus, and in those of Methanobacteriales and Methanococcales may be either due to undetected LGT, and/or result from an insufficient phylogenetic signal in the two protein datasets. Only in the latter case should an increased taxonomic sampling help resolving this incongruence, whereas if LGTs are responsible, the addition of more taxa should not increase resolution and will possibly add more confusion. The recent sequencing of several new genomes from Euryarchaeota now permits tackling these two alternatives.

Results and discussion
We have updated our previous datasets of the components of the translation and transcription machineries [11] to include a total of 25 Archaea. In particular, we included the psychrophilic methanogen Methanogenium frigidum [16] and the mesophilic methanogen Methanococcoides burtonii [16]-two lineages belonging to Methanomicrobiales [17] and Methanosarcinales [18], respectively. We also included the halophile Haloferax volcanii [19], the Thermococcale Thermococcus gammatolerans (Yvan Zivanovic and Fabrice Confalonieri, personal communication), and Nanoarchaeum equitans, a highly divergent archaeon that has been suggested as the representative of a new archaeal phylum, the Nanoarchaeota [20][21][22]. As in our previous studies [5,11] we did not include any eukaryotic outgroup in order to limit biases due to LBA.
As described previously [11], separated phylogenetic analyses were performed on each of these new datasets in order to identify and remove potential lateral gene transfer (LGT) events (data not shown). Despite the fact that  most relationships were largely unresolved in several trees due to the small size of most proteins, we checked for any possible strongly supported departure from undisputed nodes in the archaeal phylogeny, such as the clades of Thermoplasmatales, Halobacteriales, Sulfolobales, Thermococcales, Methanosarcinales and Methanococcales. Following the addition of novel taxa, no new clear-cut case of LGT could be observed with respect to these nodes, confirming that transfers are indeed very rare for these markers [5,11]. A few proteins gave an instable placement for Nanoarchaeum equitans. However, since the position of this taxon in the archaeal phylogeny has not yet been firmly tested, we did not judge these proteins as clear-cut cases of LGT. The 53 ribosomal proteins and the 14 proteins involved in transcription were thus concatenated into two large 'translation' and 'transcription' datasets, whose sizes were 6384 and 2928 amino acid positions, respectively.
Exhaustive Maximum Likelihood (ML) searches were performed on the two updated translation and transcription fusion datasets, with a few constraints given to undisputable nodes (i.e. supported by BV = 100% in preliminary Neighbor Joining and ML heuristic analyses (not shown)). The best ML topologies for the translation and transcription datasets are presented in Figure 2A and 2B, respectively. Topologies not significantly less likely than the ones presented in Figure 2 differed by minor rearrangements on nodes that are feebly supported by bootstrap, such as the branching order within halobacteriales in the transcription tree, or the grouping of methanopyrus with methanococcales/methanobacteriales in the translation tree (not shown).
Interestingly, the addition of five new archaeal taxa led to convergence of the transcription and translation trees on a coherent phylogenetic pattern, to the exclusion of the position of M. kandleri-still emerging after Thermococcales in the translation tree (Figure 2A), but recovered as a very long branch at the base of the euryarchaea in the transcription tree ( Figure 2B)-(the same trees were obtained when removing M. kandleri from the datasets, data not shown). This only incoherence between the two phylogenies is most likely due to the fact that the LBA artefact affecting the position of M. kandleri in the transcription tree persisted even after increase in taxonomic sampling, due to unavailability of closely related and slower evolving species. Nanoarchaeum equitans emerged as a separate branch distinct from those leading to Crenarchaeota and Euryarchaeota domains, in both translation and transcription trees (Figure 2A and 2B), supported by strong bootstrap values (BV = 100%). This position is congruent with previous results based on ribosomal proteins concatenation [22]. T. gammatolerans branched off of at the base of Thermococcales, that were confirmed as the first emerging euryarchaeal phylum, as in our previous studies [5,11]. Interesting, Methanobacteriales and Methanococcales formed now a monophyletic group in both translation and transcription trees (BV = 55% and 82%, Figure 2A and 2B, respectively). This suggests that the paraphyly of these groups observed in our previous translation tree ( Figure  1A) was likely incorrect due to a lack of phylogenetic signal rather than to a LGT bias. The Methanobacteriales/ Methanococcales monophyletic group is sister to a large cluster including both methanogenic and non-methanogenic species: A. fulgidus, the three Thermoplasmatales, the three Halobacteriales and the five Methanomicrobia (Methanomicrobiales and Methanosarcinales) (BV = 97% and BV = 77%, Figure 2A and 2B, respectively). This supports the hypothesis of an ancient origin of methanogenesis in Archaea followed by subsequent loss in some lineages (A. fulgidus, Thermoplasmatales and Halobacteriales). Moreover, the position of A. fulgidus, while left uncertain in our previous analyses (i.e. either sister-group of Thermoplasmatales in the translation tree, or of Halobacteriales/Methanomicrobia in the transcription tree, Figure 1A and 1B, respectively) was now robustly indicated as sister to Methanomicrobiales, Methanosarcinales and Halobacteriales in both translation and transcription trees (BP = 85% and 96%, respectively, Figure 1A and 1B). The strong placement of A. fulgidus in our updated translation tree is likely due the stabilisation of the node following addition of new taxa. This result further supports the hypothesis of a late and independent emergence of aerobic respiration in Euryarchaeotes (Halobacteriales), possibly via the recruitment of bacterial genes. Finally, both translation and transcription trees confidently grouped M. burtonii and M. frigidum with the three Methanosarcina (BV = 100% and BV = 100%, Figure 2A and BV = 100% and BV = 79%, Figure 2B, respectively) within the Methanomicrobia group. The very close relationship between M. burtonii and the three Methanosarcinales constitutes a novel phylogenetic argument justifying its inclusion in the order Methanosarcinales, at present based only on 16S rRNA phylogeny [18].

Conclusion
The congruence we obtained between the archaeal phylogenies based on the components of the translation and transcription machineries strongly supports the existence of a core of genes that evolved mainly through vertical inheritance in Archaea, and carry a bona fide phylogenetic signal that can be used to infer the phylogeny of this domain. Our results confirm also that the addition of new taxa strongly improves phylogenetic inference, and support the idea that evolutionary considerations should be included in the choice of new genomes to be sequenced. However, our conclusions should not be considered as the "last word" on the subject. For example, the misplacement of M. kandleri at the base of Euryarchaea in the tran-   scription tree was not cured by the increase in taxonomic sampling. The inclusion of sequences from slower evolving and close relatives, when they will be available, may help resolving this bias. Similarly, the very long branch displayed by N. equitans (Figure 2A and 2B) suggests that its placement as a separate branch distinct from that leading to Euryarchaeota and Crenarchaeota, although congruent between the transcription and translation trees, should be taken with caution due to the risk of an LBA artefact. The analysis of the components of additional molecular systems and the inclusion of more taxa may eventually lead to a confident placement for these two interesting species in the archaeal phylogeny.
Finally, our results make us confident that the construction of a phylogeny that retraces the vertical history of the archaeal domain is a feasible task. The identification and analysis of additional molecular markers not affected by LGT on large phylogenetic scales and their phylogenetic analysis by approaches that minimise reconstruction artefacts should continue defining the emerging picture of a genuine phylogenetic core for the Archaea. The application of a similar strategy to the bacterial and eukaryal domains could also lead to a bona fide reconstruction of their respective evolutionary histories.

Methods
In order to update the datasets of our previous analysis [11], we included the two methanogens Methanogenium frigidum [16] and Methanococcoides burtonii [16], the halophile Haloferax volcanii [19], the Thermococcale  [23]. For each dataset, novel sequences were manually added to previous alignments by using the ED program of the MUST package [24]. Regions were the alignment was ambiguous were removed from the each dataset.
Trees were computed by a number of different approaches. Neighbor-Joining (NJ) trees were calculated by the NEIGHBOR program of the PHYLIP package [25], using Maximum Likelihood (ML) distance matrices (JTT model including a Γ-correction) computed by TREE-PUZ-ZLE 5.1 [26]. Heuristic ML trees were computed using PHYML with the JTT model including a Γ-correction [27]. Exhaustive tree topology searches with limited constraints were performed using PROTML of the MOLPHY package [28]. The likelihoods and branch lengths of ML topologies were performed by TREE-PUZZLE (JTT model including a Γ-correction). For exhaustive ML searches, constraints (asterisks in Figures 1 and 2) were given to undisputable nodes (supported by BV = 100%), based on preliminary NJ and ML heuristic analyses (not shown).
The SEQBOOT program of the PHYLIP package [25] was used for the generation of bootstrapped datasets, and PUZZLEBOOT [29] and CONSENSE in the PHYLIP package [25] were used for bootstrap value calculations on 1000 replications and consensus tree reconstructions, respectively.
Datasets and their corresponding phylogenies are available on request from CB.


text:
The rarity of gene shuffling in conserved genes

The incidence of gene shuffling is estimated in conserved genes in 10 organisms from the three domains of life. Successful gene shuffling is found to be very rare among such conserved genes. This suggests that gene shuffling may not be a major force in reshaping the core genomes of eukaryotes.

Background
How do genes with new functions originate? This remains one of the most intriguing open questions in evolutionary genetics. Three principal mechanisms can create genes of novel function: point mutations and small insertions or deletions in existing genes; duplication of entire genes or domains within genes, in combination with mutations that cause functional divergence of the duplicates [1][2][3]; and recombination between dissimilar genes to create new recombinant genes (see, for example [4,5]). We here choose to call only this kind of recombination gene shuffling, excluding, for example, duplication of domains within a gene. In such a gene shuffling event, the parental genes may be either destroyed or preserved [6]. Gene shuffling is clearly the most potent of the three causes of functional innovation because it can generate new genes with a structure drastically different from that of either parental gene. Laboratory evolution studies show that gene shuffling allows new gene functions to arise at rates of orders of magnitudes higher than point mutations [7,8].
Much is known about rates of point mutations [9] and of gene duplications [10,11]. In contrast, the rate at which gene shuffling occurs is relatively unexplored, despite the importance of shuffling for functional innovation. To be sure, anecdotal evidence suggests that successful gene shuffling occurs and that it creates genes with new functions [4]. In particular, proteins are often mosaics of domains that are characterized by sequence and structural similarity [12][13][14][15][16][17][18][19]. Many domains occur in multiple proteins of different functions, suggesting that new proteins can arise through the combination of domains of other proteins, a process requiring recombination. In addition, many studies have systematically identified one subclass of gene-recombination events -gene fusions [20][21][22][23][24]. These studies count gene fusion events in a genome of interest relative to multiple, often very distantly related, species. Because fused genes often have similar functions, identification of fusion events can aid in inferring gene functions. Here we address a question that goes beyond the above studies: how frequent is gene shuffling in comparison with other forces of genome change, such as gene duplication? This problem is difficult because of the many possible outcomes of recombination events. These outcomes fall into three principal categories, gene fusions, domain deletions, and domain insertions (Figure 1a). To identify these outcomes systematically on a genomic scale is computationally intensive, which has limited our analyses to a modest number of genomes (Table 1).
One can identify gene-shuffling events either from protein sequence information or from information about protein structure. Structure-based approaches [12][13][14][15] have the advantage of being able to detect recombination events where sequence similarity between a recombination product and its parents has eroded beyond recognition. However, because two very distantly related structural domains can also have arisen through convergent evolution [25,26], identifying common ancestry of two domains based on structure alone can be problematic. As a further limitation, structure-based approaches can only identify recombination events that respect the boundaries of protein domains, whereas some successful recombination events may occur within domains [27][28][29]. In addition, structural information is not available for all genes. For example, the Pfam database of protein domains [30] contains no structural information for more than 40% of proteins in budding yeast (Saccharomyces cerevisiae). Structure-based approaches may thus miss many shuffled genes. Because of these issues we chose a sequencebased approach which allows us to search for shuffling events without making restrictive assumptions regarding their nature. Essentially, our search imposes no restrictions on shuffling except that it must merge in a single gene two protein-coding sequences that were previously a part of two different genes. We thus avoid assuming that shuffling occurs only at domain boundaries or with certain recombination mechanisms without precluding either possibility. Our analysis can also account for gene-duplication events in either parental or recombined genes.
We here identify gene-shuffling events that have occurred in a 'test' species T since its divergence from a reference species R1. A gene in the test genome whose parts match more than one gene in the reference genome is a candidate for a geneshuffling event that has occurred since the common ancestor of the two genomes. Our analysis also uses a third genome (reference genome R2) to prevent gene fission or gene loss in the reference genome R1 from resulting in spurious identification of gene shuffling events. Because R2 is an outgroup relative to T and R1, it allows us to detect such events in R1 (see Figure 1b). Like any comparative sequence-based approach, our analysis depends on detectable sequence similarity among genes. In other words, our analysis excludes rapidly evolving genes.

Little gene shuffling in closely related genomes
Mosaic proteins are not rare in most genomes, which suggests that successful gene shuffling might be frequent on an evolutionary timescale. We thus searched four closely related genomes for shuffled genes. These genomes fit three essential criteria for this analysis: close taxonomic spacing; availability of complete genome sequence; and, most important, reliable gene identification. (Gene identification is notoriously unreliable in higher organisms because of their complex gene structure.) These species were the four yeasts Saccharomyces cerevisiae, S. paradoxus, S. bayanus and S. mikatae [31], which diverged from their common ancestor between 5 and 20 million years ago (Mya) [31]. We found multiple candidate genes for shuffling in different T-R1 pairs of these four species. However, almost all of these candidates proved spurious for a variety of reasons: First, some of them occurred in two or more species in a manner inconsistent with these species' phylogeny, or they matched more closely a single reference species gene than their two putative parents. Both observations make recent recombination an unparsimonious explanation for a gene's origin (Figure 1b). Second, the putative shuffled domains in some candidate genes had a synonymous, or silent, nucleotide divergence from their parental domains that differed by a factor of two or more. However, the recombined parts of a shuffled gene should show equal sequence divergence to their respective parental genes, because they have identical divergence times (namely the time since T and R1 shared a common ancestor). We used silent nucleotide substitutions as an indicator of sequence divergence because such substitutions are under little or no selection and thus accumulate in an approximately clock-like fashion [32]. Use of amino-acid changing (nonsynonymous) substitutions (K a ) as an indicator led to similar conclusions. After exclusion of all such spurious genes, only two potential shuffled genes remained in our analysis, which indicates a low incidence of gene shuffling.

Shuffled genes in distantly related genomes
Because our analysis of yeast genomes suggests that gene shuffling may be rarer than one might expect, the need arises to study more distantly related genomes. This raises two principal problems. First, such an analysis will miss events where either parental or shuffled genes have diverged beyond sequence recognition since two genomes shared a common ancestor. We thus emphasize that our analysis applies only to 'core' genomes: genes so well conserved that their homology even among distantly related species is beyond doubt. The incidence of shuffling among more rapidly evolving genes may be different and cannot be estimated with this approach. In this regard, we also note that our analysis cannot simply use multiple outgroups for a given test genome [20][21][22][23][24] to solve this problem, because doing so has the potential to misestimate shuffling rates by making wrong assumptions about the most parsimonious placement of such events (especially among prokaryotes, where horizontal transfer of shuffled genes may occur). For the remainder of our analysis, we chose ten distantly related genomes ( Table 1) that best met the joint requirements of well known phylogenetic relationships and reliably annotated genome sequences (which is often problematic for the higher eukaryotes).
In addition to raising problems, the comparison of distantly related genomes also has one advantage: such genomes are more likely to be annotated independently from each other than are closely related genomes. In a group of closely related genomes, the first sequenced genome may often be used as a guidepost to annotate the other genomes, which may lead to errors (for instance, by misidentifying a shuffled region as an intron).
The number of shuffled genes we found is modest even for anciently diverged species pairs. For example, only 82 Identifying gene shuffling  gene-shuffling events among the 5,800 genes considered ( Table 2) may have been preserved in Caenorhabditis elegans since its common ancestor with Drosophila melanogaster, which lived around 600 Mya [33]. Similarly, only four surviving recombination events (out of 2,300 genes) may have occurred in the budding yeast Saccharomyces cerevisiae since its split from the fission yeast Schizosaccharomyces pombe more than 300 Mya [34]. We emphasize that all these numbers refer to shuffled genes that have 'survived': extant genome sequences alone are insufficient for estimating the frequency of the recombination events themselves, since the products of these events often will not become fixed in populations.
One further observation indicates the rarity of gene shuffling: most shuffled genes contain at least one domain of low sequence similarity to a parental gene. The above analysis is based on identifying sequence domains as homologous in a parental and recombined gene if they show more than 35% amino-acid sequence identity. Increasing this identity threshold to 40% can reduce the number of candidate shuffled genes dramatically (see Table 2). For instance, it removes 28 of 48 shuffled fruit fly genes and half of the shuffled fission yeast genes. This observation underscores that shuffling is rare among highly conserved genes: otherwise we would see higher sequence similarities among parental/recombined domain pairs. Figure 2 shows representative examples of shuffled genes, illustrating some of the types of recombination diagrammed in Figure 1a. For example, Figure 2b shows the budding yeast his4 gene, which is involved in histidine biosynthesis. This (apparent fusion) gene appears to combine the functions of the two fission yeast genes his7 (a phosphoribosyl-AMP cyclohydrolase) and his2 (a histidinol dehydrogenase) [35]. Figure 2c shows the fruit fly gene Aats-tyr, a tyrosyl-tRNA synthetase [36]. This gene is a likely recombination product of a predicted worm methionyl-tRNA synthetase gene mrs-1 [37] and a second worm gene Y105E8A.19 of unknown function. A list of all shuffled genes identified in these ten genomes is available in Additional data file 1.

Gene shuffling and structural domains
Because our approach is based on sequence domains, we wished to find out whether the recombined regions of shuffled genes match structural protein domains. If so, this would indicate that successful recombination events -events preserved in the evolutionary record -occur mostly at structural domain boundaries. To address this question, we used the Pfam database [30,38] of protein domains to identify domains in our shuffled genes that were significant at E ≤ 10 -5 . These Pfam domains were compared to the sequence alignments that we used to identify shuffled genes in the first place. As Figure 3 shows, the boundaries of recombined sequence domains and Pfam structural domains tend to coincide (P < 0.001 using a domain randomization approach, see Materials and methods). However, Figure 3 also suggests that not all successful recombination events occur at structural domain boundaries. Experimental and computational work on individual proteins [27] supports the notion that successful recombination occurs preferentially, but not exclusively, at structural domain boundaries. This was done because K s is a pairwise distance, meaning that it gives the sum of the divergences from the common ancestor to T and from the common ancestor to R1. The same was done for the K a analysis.

Gene shuffling and exon-intron structure
The exon-shuffling/introns-early hypothesis [39][40][41] predicts that exon-intron boundaries delimit functional domains and hence that recombination events that preserve exons would be more likely to yield functional recombinant proteins. Long introns also increase the probability of a DNA-level recombination event preserving exons (since in this case the number of possible DNA-level recombination events leading to the same recombinant protein may be quite large), a further reason to expect an association of shuffling boundaries and exon boundaries. The two multicellular eukaryotes (Drosophila and C. elegans) have a sufficient number of introns to allow us to test this prediction by comparing the boundaries of recombined sequence domains to the positions of introns in the sequences in question. However, contrary to these expectations, we found no tendency for our shuffling boundaries to associate with exon-intron boundaries (P > 0.1, domain randomization test; see Materials and methods).

The incidence of gene shuffling
We cannot estimate the incidence of gene shuffling in absolute (geological) time, because divergence dates for most of our test species are unknown or highly uncertain. In addition, the rarity of gene-shuffling events further complicates such estimates. However, we can obtain order-of-magnitude estimates of the incidence of gene shuffling relative to the Representative examples of shuffled genes identified A fusion of the fission yeast genes his7 (a phosphoribosyl-AMP cyclohydrolase) and his2 (a histidinol dehydrogenase) to produce the budding yeast his4 gene, which is involved in histidine biosynthesis. The budding yeast gene appears to combine the functions of the two fission yeast genes [35]. (c) The fruit fly gene Aats-tyr is a tyrosyl-tRNA synthetase (Flybase annotation) [36]. It is a probable recombination product of a predicted worm methionyl-tRNA synthetase gene mrs-1 (WormBase annotation) [37] and a second worm gene Y105E8A.19 of unknown function. (d) C. elegans gene ceh-20, which encodes a homeodomain protein. This gene appears to be the result of a domain exchange between the Drosophila genes exd (extradenticle, also a homeodomain protein) and Pkg21D (cGMPdependant protein kinase). (e) E. coli b4343, a hypothetical protein apparently formed via a domain exchange between Salmonella genes STY4850 (annotated as a DEAD-box helicase-related protein) and STY4851 (hypothetical protein). The numbers in the recombinant gene box are amino-acid positions in the protein product, indicating the portion of the protein derived from each of its 'parental' proteins. incidence of other mutational events important in genome evolution. One such event is gene duplication, whose incidence has been estimated previously [10,11].
To compare the incidence of gene duplication to that of gene shuffling, we cannot rely on the silent nucleotide divergence among duplicate genes to estimate the rate of duplication, as is commonly done [10,11], because several of our study genomes are very distantly related. We thus estimated the rate at which gene duplications occurred in a test species T since its common ancestor with R1 using the following approach. We identified, for each test species T, all genes that had only a single homolog in the reference species R1. We denote the number of these reference species genes as g. Second, for each of these genes i we determined the number n i of test species genes homologous to gene i. If this number n i is greater than 1, then the test species homolog of gene i underwent one or more duplications since the common ancestor of T and R1. We estimated the (minimal) number of duplication events necessary to establish a gene family of size n i as . The total estimated number d of gene duplications for the g reference species genes then calculates as the sum . Values for g and d are shown for each reference species in Table 2. One can view the ratio d/g as the per-gene incidence of gene duplication.
We then used this ratio to estimate the ratio of gene-shuffling events per gene duplication event (Table 1). To do so, we first had to estimate the number of gene-shuffling events per gene, which we obtained by dividing the number s of gene-shuffling events in a test species T (Table 1) by the average number h of genes in T or R1 with detectable sequence similarity to genes in the other genome (Table 2). This approach of estimating the number of gene-shuffling events per gene compensates for the reduced ability to recognize gene homology in distantly related genomes. The ratio of shuffling events per duplication can then be calculated as (s/h)/(d/g). Figure 4a compares this ratio for the organisms we studied. The bacteria analyzed share with the archaeans a high incidence of gene shuffling relative to duplication, while the eukaryotes show a much lower incidence. The Bacillus species (B. anthracis and B. cereus) have a much higher relative incidence of gene shuffling than any other species pair we studied.
Other mutations useful to calibrate the incidence of gene shuffling are nonsynonymous (amino-acid replacement) and synonymous (silent) mutations on DNA. Synonymous substitutions are an indicator of divergence time between two genes or species because they are subject to few evolutionary constraints and thus may change at an approximately constant (neutral) rate [32]. We estimated the incidence of gene shuffling relative to synonymous substitutions by first determining the average fraction, K s , of synonymous nucleotide changes per synonymous nucleotide site for 100 orthologous genes in a T-R1 species pair. We then simply divided the number of gene-shuffling events per gene (s/h) by this average K s (Figure 4b). The evolutionary distance of two of our species pairs (E. coli vs Salmonella and B. anthracis vs B. cereus was sufficiently low to allow us to directly calculate the average synonymous divergence for 100 pairs of randomly selected single-copy orthologs in the test and reference species (see Materials and methods). For the other species pairs, most synonymous sites are saturated with substitutions [32]. In these cases, we thus extrapolated the value of K s between R1 and T from that observed between either of these species and a third, closely related species (see Materials and methods for details). We emphasize that this procedure would be unsuitable to make evolutionary inferences for any one gene, because it introduces considerable uncertainty into our estimates. It is, however, adequate to identify the approximate, genome-wide patterns we are concerned with.
Finally, we also estimated, completely analogously, the number of gene-shuffling events per unit amino-acid divergence (K a = 1). These results are summarized in Table 1 and Figure 4c. The incidence of gene shuffling relative to silent and amino-acid divergence varies less systematically among the domains of life than that of gene shuffling relative to duplication. However, it is again apparent from these analyses that successful gene shuffling is very rare for conserved genes. For some species, crude estimates of the absolute geological time needed for two sequences to accumulate a Association between recombined sequence domains and Pfam structural domains pairwise divergence of one silent nucleotide substitution per silent site are available. In the fruit fly this amount of time is approximately equal to 64 million years [32]. During this period of time, we would expect only 5,864 × 2.8 × 10 -3 = 16 gene shuffling events to occur (Tables 1 and 2; 5,864 is the average number of fruit fly and worm genes in our core gene set). By way of comparison, even using our very conservative method of counting duplicate genes, we would expect 146 gene duplications in this period. Similarly, in the yeast S. cerevisiae, where K s = 1 synonymous substitutions accumulate every 100 million years [42], one would expect three shuffling events during this period of time (2,365 × 1.3 × 10 -3 ), as compared to 200 gene duplications. We emphasize that these are order-of-magnitude estimates that mainly serve to underscore the rarity of successful gene shuffling.
A multidomain protein may include both distinct and repeated structural domains [13]. Multidomain proteins with repeated domains raise a special problem for identifying gene-shuffling events: a shuffling event followed by domain duplication might lead us to miss a shuffled gene because our local alignments of that gene to its parental genes would only include one copy of the duplicated domain and hence might reveal less than the 50% of alignable amino-acid residues we require (see Materials and methods). To assess whether this problem would substantially bias our results, we examined candidate shuffled genes that had been excluded by our criterion (that is, those having between 10% and 50% of their residues alignable). We asked whether a failure to account for domain duplication was responsible for their exclusion. After adding potentially duplicated domains to the aligned regions of these genes, we found that only a handful of them (two genes in Drosophila and three in C. elegans) met our 50% alignability threshold. Failure to account for domain duplications internal to a gene is thus not the reason for our low estimates of the incidence of gene shuffling.
Several lines of evidence show that successful gene shuffling is very rare for genes conserved between the distantly related genomes we studied. For the single-celled yeasts -currently the only group of very closely related eukaryotes with sufficiently reliable genome annotation -shuffling appears rare in the genome as a whole. In most of the genomes we analyzed, gene shuffling is much rarer than other important kinds of mutations affecting gene structure, such as gene duplication. For example, in the time that it takes to accumulate K s = 0.01 synonymous substitutions per synonymous site, other research indicates that ten fruit fly genes and 164 worm genes undergo duplication [10]. In contrast, each lineage has only a 50% chance of undergoing a successful gene shuffling event in the same amount of time (if one assumes our estimates can be applied to the entire genome). We note that our estimates of duplication rates are more conservative than those of others [10], partly because we limit ourselves to single-gene duplications. The fact that we still see a lower Incidence of gene shuffling relative to various other mutational events  incidence of shuffling than duplication (Figure 4) is thus all the more remarkable.

Discussion
The rarity of gene shuffling relative to gene duplication has a simple potential explanation. A gene duplication creates a copy of a gene while preserving an original that is able to exercise its function. In contrast, unless a recombination event is accompanied by gene duplication, the original (parental) genes disappear in the event. An organism may survive a recombination event only if neither parental gene was essential to its survival and reproduction or if the recombinant gene(s) can carry out the function of both parental genes. This rarity of successful gene shuffling stands in stark contrast to the frequency of DNA recombination itself, which is a ubiquitous process accompanying DNA replication and repair. This suggests that the vast majority of recombined genes have deleterious effects on the organism, which may be particularly true for the highly conserved genes examined here.
We emphasize that the rarity of gene shuffling we find is not in contradiction with earlier studies that have identified multiple gene fusions -a special, simple case of gene shuffling -in fully sequenced genomes [20][21][22][23][24]. These studies identified prokaryotic gene fusion events in one test genome relative to multiple, often very distantly related, reference genomes. Any such approach may find many fusion events even if such events are rare. Our data also do not rule out the possibility that shuffling played an important role in forming the conserved eukaryotic core genome, because the pertinent geneshuffling events would have occurred before the divergence of the eukaryotic species pairs we examined. (The identification of such ancient shuffling events may require an approach based on protein structure.) Furthermore, our results are not in contradiction with anecdotal evidence for the abundance of gene-shuffling events in some functional categories of genes [32]. The reason is that our results pertain to the average incidence of gene shuffling among conserved genes. Some genes may be shuffled at a much greater rate. Indeed, structural studies of multidomain proteins tend to find a few domains which co-occur with a wide variety of other domains (indicating the common shuffling of such promiscuous domains), whereas many other domains co-occur with only one or a few other domains (rare shuffling) [43]. Similarly, a lack of reliable genome annotation made it impossible to reliably identify gene-shuffling events in vertebrate genomes, where gene shuffling may be more frequent overall [44].
Perhaps the central caveat to our results regards sources of ascertainment bias. The comparison of distantly related genomes alone introduces a powerful source of ascertainment bias: we can only analyze gene-shuffling events for genes that have been sufficiently conserved to be recognizable in both genomes. However, shuffling might be more common among rapidly evolving genes. An additional possible source of bias is that after a successful gene-shuffling event the rate of amino-acid substitutions may be elevated as a result of directional selection on the newly created gene. Such a bias would cause us to underestimate shuffling frequencies in distantly related species even for conserved genes. Nonetheless, our results from the four closely related genomes argue against such a bias, because shuffling also appears rare in these genomes.
Another caveat is that our ability to identify successful geneshuffling events depends on the continued presence of both parental genes in the reference genome. Genomes, however, occasionally lose genes. For instance, recent work has suggested that S. cerevisiae has lost roughly 10% of its genes since its last common ancestor with S. pombe [45]. If gene loss in other organisms occurs at comparable rates, our approach may slightly underestimate the number of recombination events in a lineage. However, note that gene loss affects our estimates of gene shuffling and gene duplication in similar ways, thus compensating for any such bias.
We used a second reference genome R2 to be able to exclude gene-fission events in reference genome R1. Such events can lead to misidentification of recombination in the test species and have been documented in several organisms [20]. Unfortunately, this approach fails if the same recombination event occurred twice, once in the lineage leading to reference species R2 and once in the lineage leading to test species T. Such a case of parallel evolution or homoplasy would lead us to misidentify a recombination event in T as a gene-fission event in R1. However, because successful gene shuffling is very rare in general, and because the required recombination event would have to occur at exactly the same position twice, this possibility is probably not a major confounding factor in our analysis.
A fourth caveat lies in the possibility that some of our recombinant genes may result from two independent recombination events. Our algorithm can identify such genes, but given the high sequence divergence of recombinant domains it may often be impossible to resolve the order of the individual recombination events. The generally small number of recombinant proteins implies that genes produced by two or more recombination events would be extremely rare. Indeed, among 203 identified recombinant genes, a mere 16 show matches to more than two parental genes, making these the only cases with indications that more than one recombination process was involved in their creation.
Finally, our approach to estimating the rate of gene duplication identifies only duplications of single-copy genes in the reference species. Multicopy genes may undergo duplication more frequently. We may thus have underestimated the number of gene duplications. As a result, the incidence of gene shuffling relative to gene duplication may be even lower than indicated by our estimates.

Recombination and introns
Our findings speak to a long-standing debate in molecular evolution, a debate that revolves around the origin of introns.
Introns are stretches of DNA that do not code for proteins and that separate exons, the protein-coding regions of genes.
According to one point of view, introns originated early in the evolution of life, perhaps as early as the common ancestor of prokaryotes and eukaryotes [39][40][41]. According to this perspective, introns may have acted as spacers between exons and thus greatly facilitated recombination among exons to create new proteins. The opposite point of view is that introns arose late in life's evolution, perhaps as late as eukaryotes themselves [28,29] and thus had no role in gene shuffling earlier in life's history. Genes in two of our test genomes have a sufficient number of introns to test the hypothesis that introns facilitate gene shuffling. Neither of these genomes showed an association between gene-shuffling boundaries and exon position. In addition, neither of these genomes showed an elevated incidence of gene shuffling. Although based on a small number of genomes, this finding casts doubt on the importance of introns for gene shuffling, and it suggests that other aspects of genome architecture may be more important. One potential example is the organization of functionally related prokaryotic genes into operons. The close proximity of such genes may facilitate their reorganization and the generation of new functions, whether through simple fusion or fission or through more radical change.

Natural selection or drift?
A nonhomologous recombination event that gives rise to a shuffled gene occurs in only one individual of a potentially large population. Does a shuffled gene typically rise to high frequency and become fixed through natural selection or genetic drift? To answer this question, one could in principle study the relationship between the rate at which fixed shuffled genes arise and population size (taking account of differences in nonhomologous recombination rates among species). Three possibilities exist in principle. First, there may be no relation between population size and the rate at which fixed shuffled genes arise. This would be the case if most gene-shuffling events are strictly neutral [46] or if they have very large beneficial effects. Second, there may be a positive relation between the rate at which fixed shuffled genes arise and population size. This would be the case if most shuffling events are mildly beneficial. The reason is that in this case selection favoring the fixation of a shuffled gene has to overcome the effects of genetic drift, which are weakest in large populations. Finally and perhaps most likely, there may be a negative association between population size and the rate at which fixed shuffled genes arise. This would be the case if most shuffling events are mildly deleterious. Such a negative association has been observed for several indicators of genome structure such as genome size, transposable element load, and rates of preservation of duplicated material [47].
Unfortunately, insufficient data are available to distinguish rigorously between these possibilities. First, estimates of effective population sizes N e , based on estimates of N e µ [47] and the mutation rate µ [9], exist only for three of our five pairs of genomes (E. coli-Salmonella, S. cerevisiae-S. pombe and D. melanogaster-C. elegans). Second, we have insufficient information on recombination rates (whose variation among genomes needs to be taken into account). Specifically, although estimates of homologous recombination rates are available for a few of our organisms [48][49][50][51], gene shuffling occurs strictly by nonhomologous recombination, whose rate need not have a simple relationship with the homologous recombination rate. A third difficulty is that recombination rates and mutation rates are conventionally measured per cycle of DNA replication, whereas we would require per-year estimates as well as estimates of absolute divergence times between our taxa of interest to make appropriate comparisons.
Despite such insufficient data, we can make the qualitative observation that the observed incidence of shuffling does not follow a simple pattern: For example, S. cerevisiae has a relatively high effective population size (N e µ is approximately half of that for E. coli [47] while µ is actually higher than that of E. coli [9]) and a high homologous recombination rate compared to C. elegans or E. coli [48][49][50][51], and yet it shows the lowest incident of gene shuffling of any of our taxa. In the slightly deleterious scenario above, we would instead expect yeast to show an incidence of shuffling greater than that of E. coli, while in the slightly beneficial scenario we would expect it to show an incidence greater than that of the multicellular eukaryotes.
A second qualitative observation is that the incidence of gene shuffling is not elevated in higher organisms relative to the rate of nucleotide substitutions. (The higher incidence of gene shuffling relative to gene duplication in prokaryotes from Figure 4a may be a consequence of the lower rate of gene duplication in these taxa.) This is consistent with the hypothesis that the fate of most shuffled genes is driven by natural selection rather than genetic drift. In other words, most shuffling events may not be neutral. This is again plausible if one considers that most gene-shuffling events change a gene's structure drastically. A corollary of this hypothesis is that preserved shuffled genes have been preserved for a reason -the benefit they confer to an organism. While rare in number, shuffled genes may thus be of great importance in organismal evolution.
Our analysis of gene shuffling has left many open questions, most notably about the association between the rate of sequence evolution and the rate of gene shuffling. To arrive at firm answers for this and other questions, we must be able to study shuffling rates not only for conserved proteins but also for rapidly evolving proteins. Such studies will require closely related genome sequences with reliable gene identification derived independently for each genome.

Identifying shuffled genes
Our method identifies shuffled genes in a test genome (T) relative to a reference genome (R1). Table 1 shows the ten test genomes -two archaeal, six prokaryotic, and four eukaryotic genomes -we used in this analysis. Every pair of genomes R1-T occurs twice in Table 1, because one of two genomes can be used either as the test or the reference genome. To exclude spurious recombination events that reflect gene loss or fission in R1, the method also employs a second reference genome, R2. The two archaeans in our analysis were Pyrococcus horikoshii [52] and Methanocaldococcus jannaschii [53].
To identify sequence homology between all genes in these genomes we used the Washington University implementation of gapped BLASTP [66,67], followed by exact pairwise local alignment using the Smith-Waterman algorithm [68] with a gap-opening penalty of 10 and a gap-extension penalty of 2, and the BLOSUM 62 scoring matrix [69]. We excluded from further analysis all gene pairs with BLAST E-values greater than 10 -6 , fewer than 50 aligned amino acids, amino-acid identity in the alignment of less than 35%, or alignments consisting of more than 50% low-complexity sequences as determined by the SEG program [70,71].
The requirement of 35% sequence identity may appear to bias our estimates of shuffling incidence between distantly related taxa. However, because we calculate these values relative to the total number of genes with the same (35%) degree of sequence identity between the test and reference genome (h), this bias is most likely to be small.
The result of this procedure is a list of partially or fully matching genes in the two species T and R1. We used this list to identify shuffled genes in the test genome T. Specifically, for each gene in the test genome T we searched for pairs of genes in the reference genome R1 that matched the test species gene, but in nonoverlapping or minimally overlapping regions. (To account for edge errors in local alignments, we allowed regions to overlap by a maximum of 20 residues). After having identified any such gene, we verified that it did not also have full-length homologs in the reference genome, because otherwise gene shuffling would not be the most parsimonious explanation of the gene's origin. We developed a special-purpose algorithm for this search [72], which identifies, for any one gene, the combination of local alignments to genes in the reference genome that covers the maximum number of residues in the shuffled gene. This algorithm can identify shuffled genes (genes to which two or more reference species genes contributed), but it will also return only a single alignment if this alignment is longer than any combination of non-overlapping alignments.

Three criteria for validating shuffling events
We used three additional criteria to validate candidates for shuffled genes. First, we computed the proportion of a shuffled gene's amino-acid residues that could be aligned to its (parental) reference species genes. If this proportion is small, a gene may be too highly diverged for us to confidently ascertain that it is a recombination product. We excluded genes where this proportion was smaller than 50%. This requirement may appear restrictive, but additional analyses show that our conclusions hold even if it is completely eliminated. For example, eliminating this criterion increases the number of shuffled genes by a factor ranging from 1 (no increase, E. coli) to 4.2 (C. elegans), but the eukaryotes surveyed still show an incidence of shuffling smaller than the duplication rate, while the prokaryotes show similar frequencies of shuffling and duplication. We have maintained the 50% requirement throughout our main analysis to err on the side of caution: Putative shuffled genes with very short alignable regions to a parental gene are more likely to be false positives. They also do not belong in the set of genes conserved between T and R1, which is our focus here.
To motivate our second validation criterion, we note that in the eukaryotic test genomes some shuffled genes had undergone duplication. We identified gene duplicates as gene pairs with amino-acid divergences K a < 1 using a previously described and publicly available tool [73]. We counted each gene family of shuffled genes only once to avoid doublecounting duplicates of shuffled genes.
A third indicator of true recombination is the divergence of different sequence domains within a putative shuffled gene. The recombined parts of a shuffled gene should show equal sequence divergence from their respective parental genes, if these parts have diverged in a clock-like fashion. The two principal indicators of DNA sequence divergence are the number of silent nucleotide substitutions at synonymous sites (K s ) and the number of non-synonymous substitutions at amino-acid replacement sites (K a ) [32]. We used the methods of Muse and Gaut and Goldman and Yang [74,75] to estimate these divergence indicators for our putative shuffled genes. Silent substitutions are subject to much weaker selection pressures than amino-acid replacement substitutions, and their rate of accumulation is thus more clock-like [32]. For our analysis of four closely related yeast species, we calculated the distance (K a or K s ) between each sequence domain of a shuffled gene and its counterpart parental gene. We excluded candidate shuffled genes from further analysis if these distances (K a or K s ) differed by more than a factor two between the different domains. For our highly divergent species (Table 1) we could not use synonymous divergence K s because most domain pairs had become saturated with silent substitutions. The (unavoidable) disadvantage of using only amino-acid divergence K a is that natural selection on aminoacid changes can cause non-clock-like evolution of domains. However, even with this caveat, our approach identifies true recombination events. Otherwise, our shuffled genes would not show the highly significant statistical association in amino-acid divergence K a we observe between sequence domain pairs ( Figure 5; Pearson's r 0.47, Spearman's s 0.45, P < 0.0001 for both). Note that we did not exclude shuffled genes from our analysis of distantly related genomes based on unequal K a estimates.

Detecting domain duplication in putative shuffled genes
Internal domain duplications could cause a shuffled gene to fail the first of the three test criteria above -50% alignable residues with its parental genes -which would cause us to miss such genes. To assess how serious a problem such missed genes might be, we relaxed the first of the above test criteria, requiring only 10% alignable residues between a test gene and two or more potential parental genes. For test species genes that met this criterion, we then excised the portion of the test species gene that aligned to a reference gene, and recomputed a local alignment between this trimmed gene and the reference gene. If the test-species gene contained internal duplications of the reference gene, the trimmed sequence should still align to the reference gene. We added any new alignments found in this way to the original alignment combination and assessed whether the resulting combination of alignments met the required threshold of 50% alignable residues. For alignments that met our other prescreening criteria for candidate shuffled genes (50 amino acids in length and > 35% amino-acid identity) we iterated this excision procedure to determine the total number of repeated domains. This analysis identified only five additional shuffled genes (see above) and led us to conclude that internal duplications are not a major confounding factor in our analysis.

Estimating relative frequencies of gene shuffling
We estimated the incidence of gene shuffling relative to two other frequent kinds of evolutionary events -gene duplications and nucleotide substitutions. To do so, we first needed to account for differences in genome size among our study species. We thus estimated the number of gene shuffling events per gene. This estimate poses a problem that stems from the different (and highly uncertain) times since common ancestry of our different T-R1 species pairs. The longer the time since common ancestry, the fewer genes (shuffled or not) with recognizable sequence similarity two species will share. Thus, when simply dividing the number of shuffled genes s by the total number of genes in a test genome, one may wrongly estimate the number of gene-shuffling events per gene. To account for this problem, we divided the total number of gene-shuffling events s by the number of recognizable homologs h shared between species T and R1 to obtain the number of gene shuffling events per gene, s/h. To obtain h itself, we determined the total number of genes in T with at least one homolog -using the criteria outlined earlier -in the genome of R1, the total number of genes in R1 with at least one homolog in T, and averaged these two numbers.
We then related this number of gene shuffling events per gene to the number of gene duplications in T since its common ancestor with R1. Because our test and reference genomes are only distantly related, we could not rely on the silent nucleotide divergence among duplicate genes to estimate the rate of duplication, as is commonly done [10,11]. We thus used the following, alternative, method. First, for each test species T, Similarity in sequence divergence between regions of shuffled genes Figure 5 Similarity in sequence divergence between regions of shuffled genes. The amino-acid divergences (K a ) of recombined domains to their respective parental counterparts are correlated. One outlying observation (K a1 = 3.26 and K a2 = 0.36) is not shown in this plot but was included in the calculation of correlation coefficients. Excluding this observation increases the Pearson correlation coefficient to 0.61 and leaves the Spearman correlation coefficient unchanged (P < 0.0001 for both). we identified all genes which match only a single gene in R1 at a BLAST threshold of 10 -6 and had 70% or more of their sequences aligned to that gene. The number of unique reference species genes matching one or more test species genes gives a baseline number g of genes before duplication. Second, for each such unique reference species gene i we identified duplicate pairs of test species genes that showed a pairwise amino-acid divergence K a which was less than either gene's amino-acid divergence from the putative ortholog i. For many genes i, more than one pair of test species genes fulfilled these criteria. Such genes represent families of n i duplicates of gene i. We estimated the (minimal) number of duplication events necessary to establish such a gene family of size n i as . The total estimated number d of gene duplications for the g reference species genes then calculates as the sum . Genes in very large families are more likely to undergo duplication than genes in smaller families [76]. For this reason, we excluded reference species genes i with more than ten duplicate genes from this analysis. Including such genes would tend to increase our estimated rate of gene duplication, meaning that the results shown in Figure 4a are conservative estimates of the excess of duplicates relative to shuffled genes. To estimate the rate of gene shuffling events relative to gene duplication events, we then divided the number of gene shuffling events per gene (s/h, obtained above) by the number d/g of gene duplication events per gene.
We also estimated the number of gene-shuffling events per unit of silent substitutions K s that accumulate in a gene. Two of our species pairs (B. anthracis-B. cereus and E. coli-S. enterica) allowed us to estimate synonymous divergence K s directly. For these two species pairs, we first identified 100 pairs of single-copy genes in each genome that are unambiguous orthologs [32,77]. We then divided the number of shuffled genes per gene by the average synonymous divergence K s of the orthologs with unsaturated synonymous divergence (> 97 for both species pairs) to obtain an estimate of the number of gene-shuffling events per unit change in K s .
For the three other genome pairs, this approach was not feasible because of their mostly saturated synonymous divergence. We thus had to estimate the average synonymous divergence K s between T and R1 by extrapolating from the synonymous divergence between T and a more closely related species. This approach relies on previous work [10] which indicates that the genome-wide average ratio K a /K s of aminoacid divergence to synonymous divergence approaches an asymptotic value for large numbers of distantly related genes. For each of the three T-R1 species pairs, our approach to estimating K s in this way consists of two steps. We first estimated the average amino-acid divergence K a between 100 randomly chosen unique single-copy orthologs of a T-R1 species pair. For the second step, we first identified an organism C with fully sequenced genome that is closely related to either T or R1. The genome of C should be sufficiently closely related to estimate K s reliably, but sufficiently distantly related to reliably estimate the asymptotic ratio of amino-acid to silent divergence. This organism was Saccharomyces paradoxus for S. cerevisiae; Pyrococcus furiosus [78] for P. horikoshii; and Caenorhabditis briggsae [79] for C. elegans. For each of these closely related genome pairs, we chose at random 100 single-copy gene pairs that were unambiguous orthologs. We calculated the average ratio K ac /K sc for these orthologs. We then used this value to extrapolate the average fraction K s of synonymous substitutions between genes in T and R1 as K s = K a /(K ac /K sc ). This is the estimated average synonymous substitution rate between a T-R1 species pair. We then related the rate of gene shuffling to this extrapolation of K s . Specifically, we estimated the number of shuffling events per gene per one K s as (s/h)/(K s /2). (The reason for dividing the average K s by 2 is that our approach estimates the number of gene-shuffling events only for one of the two species of a T-R1 species pair.) We are well aware of the shortcomings of this approach, which averages heterogeneous substitution rates and assumes that the ratio of amino acid to silent divergence is constant within the taxonomic group considered. However, we emphasize that we use the approach here only to arrive at order-of-magnitude estimates of the incidence of gene shuffling. We also note that although we have not explicitly taken codon usage bias into account, the use of codon position-specific nucleotide frequencies (which partially account for such a bias, at a cost of larger estimate variances) increased all of our estimated average K s values without changing the patterns seen in Figure 4. Thus, the values in Figure 4 are conservative in the sense that the actual incidence of shuffling relative to K s may be lower than shown.
Third, and finally, we also estimated the number of geneshuffling events per unit of amino-acid replacement substitutions K a that accumulate in a gene. To do so, we divided s/h as above by one-half the average amino-acid divergence (K a /2) of 100 unambiguous orthologs in a T-R1 species pair.

Comparison of identified domains to Pfam database
Because our analysis was sequence and not structure-based, we used the Pfam database of structural protein domains [30,38] to evaluate how well the recombined sequence domains we identified matched structural domains. To do so, we queried all shuffled genes against the Pfam database and retained identified Pfam domains with E ≤ 10 -5 . We then compared the location of these structural domains to the location of the sequence domains in the shuffled genes. Specifically, we calculated the distance between each alignment domain and its closest Pfam domain. a large extent. Parametric tests are not appropriate to evaluate the statistical significance of this association. We thus applied a gene-randomization procedure that created new alignment domains within each gene, domains whose starting and ending positions are uniformly distributed but cover the same proportion of the gene as did the original domains. Each randomized gene possessed the same number of simulated domains as observed domains, but with different positions and lengths. We then calculated D for the simulated alignment domains and compared its distribution to the empirically observed values of D. We applied an analogous approach to test whether exon/intron boundaries and shuffling boundaries are associated in our two multicellular eukaryotes (C. elegans and Drosophila). This approach substituted the closest exon for the closest Pfam domain and applied the same randomization procedure.

Additional data files
Additional data file 1, available with the online version of this paper, contains a table listing all shuffled genes included in our analysis of the ten distantly related genomes.
Additional File 1 A table listing all shuffled genes included in our analysis of the ten distantly related genomes A table listing all shuffled genes included in our analysis of the ten distantly related genomes. A table listing all shuffled genes included in our analysis of the ten distantly related genomes. Click here for file