Meta-analyses of studies of the human microbiota

Our body habitat-associated microbial communities are of intense research interest because of their influence on human health. Because many studies of the microbiota are based on the same bacterial 16S ribosomal RNA (rRNA) gene target, they can, in principle, be compared to determine the relative importance of different disease/physiologic/developmental states. However, differences in experimental protocols used may produce variation that outweighs biological differences. By comparing 16S rRNA gene sequences generated from diverse studies of the human microbiota using the QIIME database, we found that variation in composition of the microbiota across different body sites was consistently larger than technical variability across studies. However, samples from different studies of the Western adult fecal microbiota generally clustered by study, and the 16S rRNA target region, DNA extraction technique, and sequencing platform produced systematic biases in observed diversity that could obscure biologically meaningful compositional differences. In contrast, systematic compositional differences in the fecal microbiota that occurred with age and between Western and more agrarian cultures were great enough to outweigh technical variation. Furthermore, individuals with ileal Crohn''s disease and in their third trimester of pregnancy often resembled infants from different studies more than controls from the same study, indicating parallel compositional attributes of these distinct developmental/physiological/disease states. Together, these results show that cross-study comparisons of human microbiota are valuable when the studied parameter has a large effect size, but studies of more subtle effects on the human microbiota require carefully selected control populations and standardized protocols.Targeting our indigenous human microbial communities (microbiota) to prevent or treat disease is difficult due to their complexity, as well as their intra- and interpersonal variations (Lozupone et al. 2012b). Major efforts are underway to understand the predominant factors that shape the human gut microbiota and the inter-relationships between the organismal composition of the microbiota, its pool of microbial genes (microbiome), their expressed functions, and host physiologic and disease phenotypes. In the case of the gut, which contains the largest collection of microbes, these factors and interrelationships include diet (Muegge et al. 2011; Wu et al. 2011; Yatsunenko et al. 2012), host genetic and familial relationships (Turnbaugh et al. 2009; Hansen et al. 2011; Yatsunenko et al. 2012), varying cultural traditions and geography (De Filippo et al. 2010; Hehemann et al. 2010; Yatsunenko et al. 2012; Zupancic et al. 2012), age (Palmer et al. 2007; Biagi et al. 2010; Koenig et al. 2011; O''Sullivan et al. 2011; Yatsunenko et al. 2012), pregnancy (Koren et al. 2012), route of delivery (Huurre et al. 2008; Dominguez-Bello et al. 2010), obesity, metabolic syndrome, and type II diabetes (Ley et al. 2005; Turnbaugh et al. 2009; Qin et al. 2010; Graessler et al. 2012; Vrieze et al. 2012), cardiovascular disease (Wang et al. 2011), disturbances produced by antibiotics (Jakobsson et al. 2010; Dethlefsen and Relman 2011) including Clostridium difficile colitis (Chang et al. 2008; Khoruts et al. 2010; Gough et al. 2011), and other forms of inflammatory bowel diseases (Willing et al. 2010).Bacteria dominate our various microbial communities. The composition of these communities is typically evaluated by targeting the bacterial 16S rRNA gene as a phylogenetic marker. Trends in community-level diversity differences can be interrogated by computing the amount of diversity that is shared between samples (β-diversity), followed by clustering using an unsupervised multivariate statistical technique such as Principal Coordinates Analysis (PCoA). These techniques sometimes reveal clear associations between subject characteristics and overall diversity. Strong drivers of gut microbial community relatedness include age (Koenig et al. 2011; Yatsunenko et al. 2012), culture/geography (Yatsunenko et al. 2012), inflammatory bowel disease (IBD) (Willing et al. 2010), and kinship (Dicksved et al. 2008; Turnbaugh et al. 2009; Yatsunenko et al. 2012).Because different studies of the human microbiota often use the same 16S rRNA gene target, studies performed by different research groups can in principle be compared, and parallels among different disease, physiological, or developmental states discovered. Such comparative analyses have yielded key insights when applied to 16S rRNA gene libraries generated by different laboratories focusing on a number of environmental habitats; for example, these comparisons have revealed that salinity is an important factor structuring the bacterial diversity in free-living communities (Lozupone and Knight 2007) and that bacterial communities in the vertebrate gut are highly divergent from free-living communities (Ley et al. 2008). Analysis of sequences from different studies also allows comparisons to relevant control populations. For example, pregnant women differ in the composition of their gut microbiota between the first and third trimester (Koren et al. 2013), and third but not first trimester composition was shown to be distinctive from nonpregnant adults by comparison to the healthy reference data set sequenced by the NIH-sponsored Human Microbiome Project (HMP) (The Human Microbiome Project Consortium 2012).Particularly in comparisons restricted to a specific type of sample (e.g., only from human fecal samples), technical differences in experimental protocols between laboratories, including the manner in which samples are obtained and stored, DNA extraction methods, the selection of PCR primers for generating amplicons from bacterial 16S rRNA genes, the region of the 16S rRNA gene targeted for PCR, and the instruments used to determine the nucleotide sequences of these amplicons, might all produce variability that could outweigh biological differences (Mao et al. 2012). Here, we conducted meta-analyses to identify overall patterns that drive differences in the human microbiota and to ascertain the degree to which technical variability between studies impacts observed diversity.     

Sympatric chimpanzees and gorillas harbor convergent gut microbial communities

The gut microbial communities within great apes have been shown to reflect the phylogenetic history of their hosts, indicating codiversification between great apes and their gut microbiota over evolutionary timescales. But because the great apes examined to date represent geographically isolated populations whose diets derive from different sources, it is unclear whether this pattern of codiversification has resulted from a long history of coadaptation between microbes and hosts (heritable factors) or from the ecological and geographic separation among host species (environmental factors). To evaluate the relative influences of heritable and environmental factors on the evolution of the great ape gut microbiota, we assayed the gut communities of sympatric and allopatric populations of chimpanzees, bonobos, and gorillas residing throughout equatorial Africa. Comparisons of these populations revealed that the gut communities of different host species can always be distinguished from one another but that the gut communities of sympatric chimpanzees and gorillas have converged in terms of community composition, sharing on average 53% more bacterial phylotypes than the gut communities of allopatric hosts. Host environment, independent of host genetics and evolutionary history, shaped the distribution of bacterial phylotypes across the Bacteroidetes, Firmicutes, Proteobacteria, and Actinobacteria, the four most common phyla of gut bacteria. Moreover, the specific patterns of phylotype sharing among hosts suggest that chimpanzees living in sympatry with gorillas have acquired bacteria from gorillas. These results indicate that geographic isolation between host species has promoted the evolutionary differentiation of great ape gut bacterial communities.The compositions of the gut microbial communities harbored by great apes reflect the phylogeny of their hosts in a manner suggesting that host species and their gut microbiota have codiversified over evolutionary timescales (Ochman et al. 2010; Degnan et al. 2012). This pattern of codiversification between hosts and their gut microbiota could stem from both heritable factors, such as host genetics and the vertical, generation-to-generation transmission of gut microbes (Vaishampayan et al. 2010), and environmental factors, such as host diet and geography (Ley et al. 2008a,b; Turnbaugh et al. 2008, 2009; De Filippo et al. 2010; La Serre et al. 2010; Muegge et al. 2011; Claesson et al. 2012; Yatsunenko et al. 2012). However, because the great ape species sampled to date represent populations that are at once phylogenetically, ecologically, and geographically distinct, it has not been possible to separate the relative influences of heritable and environmental factors on the evolution of the great ape gut microbiota.One approach to parsing the effects of environmental factors on the gut microbiota from those of heritable factors is to compare sympatric (i.e., co-occurring) and allopatric (i.e., geographically separated) host populations. Gorillas diverged from the lineage leading to humans and chimpanzees/bonobos at least 6 million years ago (Glazko and Nei 2003; Langergraber et al. 2012), but since that time, the two groups have come into secondary contact throughout equatorial Africa. When living in sympatry, great ape species experience dietary convergence in addition to shared geography (Williamson et al. 1990; Tutin and Fernandez 1993; Shannon et al. 2006; Yamagiwa and Basabose 2006), but they do not mingle or interbreed, and their phylogenetic distinctiveness is maintained. Therefore, the effects on the gut microbiota of the environmental factors shared exclusively by sympatric chimpanzees and gorillas can be measured quantitatively as the degree of compositional convergence between the gut microbiota of sympatric populations relative to those of allopatric populations.To quantify the effects of shared environmental factors on the gut microbiota of sympatric great apes, we have investigated the gut microbiota of sympatric and allopatric populations of chimpanzees, bonobos, and gorillas from Tanzania, Cameroon, the Central African Republic (CAR), and the Democratic Republic of the Congo (DRC). We show that, while hosts of different species always maintain distinct gut microbiota (even when living in sympatry), the gut microbiota of sympatric Pan and Gorilla share significantly more bacterial phylotypes than do those of allopatric Pan and Gorilla. Moreover, the specific patterns of phylotype sharing indicate a history of transfer of gut bacteria between the two host species, with chimpanzees acquiring bacteria from sympatric gorillas. Recent analyses of human populations have shown how geographic factors can shape intraspecific variation in gut microbiota composition (De Filippo et al. 2010; Yatsunenko et al. 2012); our results broaden this principle to include a role for geographic isolation in maintaining differences in gut microbiota composition among closely related heterospecific hosts.     

Functional metagenomics to mine the human gut microbiome for dietary fiber catabolic enzymes

The human gut microbiome is a complex ecosystem composed mainly of uncultured bacteria. It plays an essential role in the catabolism of dietary fibers, the part of plant material in our diet that is not metabolized in the upper digestive tract, because the human genome does not encode adequate carbohydrate active enzymes (CAZymes). We describe a multi-step functionally based approach to guide the in-depth pyrosequencing of specific regions of the human gut metagenome encoding the CAZymes involved in dietary fiber breakdown. High-throughput functional screens were first applied to a library covering 5.4 × 10⁹ bp of metagenomic DNA, allowing the isolation of 310 clones showing beta-glucanase, hemicellulase, galactanase, amylase, or pectinase activities. Based on the results of refined secondary screens, sequencing efforts were reduced to 0.84 Mb of nonredundant metagenomic DNA, corresponding to 26 clones that were particularly efficient for the degradation of raw plant polysaccharides. Seventy-three CAZymes from 35 different families were discovered. This corresponds to a fivefold target-gene enrichment compared to random sequencing of the human gut metagenome. Thirty-three of these CAZy encoding genes are highly homologous to prevalent genes found in the gut microbiome of at least 20 individuals for whose metagenomic data are available. Moreover, 18 multigenic clusters encoding complementary enzyme activities for plant cell wall degradation were also identified. Gene taxonomic assignment is consistent with horizontal gene transfer events in dominant gut species and provides new insights into the human gut functional trophic chain.The human intestinal microbiome is the dense and complex ecosystem that resides in the distal part of our digestive tract. Its role in metabolizing dietary constituents (Sonnenburg et al. 2005; Flint et al. 2008; Ley et al. 2008) and in protecting the host against pathogens (Rakoff-Nahoum et al. 2004) is crucial to human health (Macdonald and Monteleone 2005; McGarr et al. 2005; Manichanh et al. 2006; Turnbaugh and Gordon 2009). It is mainly composed of commensal bacteria from the Bacteroidetes, Firmicutes, Proteobacteria, and Actinobacteria phyla (five), and of several archaeal and eukaryotic species. With up to 10¹² cells per gram of feces, the bacterial abundance is estimated to reach 1000 operational taxonomic units (OTUs) per individual, 70% to 80% of the most dominant ones being subject-specific (Zoetendal et al. 1998; Tap et al. 2009). However, only 20% of the bacterial species have been successfully cultured so far (Eckburg et al. 2005). Large-scale analyses of genomic and metagenomic sequences have provided gene catalogs and statistical evidence on protein families involved in the predominant functions of the human gut microbiome (Gill et al. 2006; Kurokawa et al. 2007; Flint et al. 2008; Turnbaugh et al. 2009; Qin et al. 2010), among which the catabolism of dietary fibers is of particular interest in human nutrition and health. Dietary fibers are the components of vegetables, cereals, leguminous seeds, and fruits that are not digested in the stomach or in the small intestine, but are fermented in the colon by the gut microbiome and/or excreted in feces (Grabitske and Slavin 2008). Chemically, dietary fibers are mainly composed of complex plant cell wall polysaccharides and their associated lignin (Selvendran 1984), along with storage polysaccharides such as fructans and resistant starch (Institute of Medicine 2005). Dietary fibers have been identified as a strong positive dietary factor in the prevention of obesity, diabetes, and cardiovascular diseases (World Health Organization 2003). Because of the wide structural diversity of dietary fibers, the human gut bacteria produce a huge panel of carbohydrate active enzymes (CAZymes), with widely different substrate specificities, to degrade these compounds into metabolizable monosaccharides and disaccharides. The functions and the evolutionary relationships of CAZyme-encoding genes of the human gut microbiome are being extensively studied through functional and structural genomics investigations (Flint et al. 2008; Lozupone et al. 2008; Mahowald et al. 2009; Martens et al. 2009), which are nevertheless restricted to cultivated bacterial species. CAZyme diversity has also been described in three metagenomics studies focused on this microbiome (Gill et al. 2006; Turnbaugh et al. 2009, 2010), and these revealed the presence of at least 81 families of glycoside-hydrolases, making the human gut metagenome one of the richest source of CAZymes (Li et al. 2009). However, the proof of function of annotated genes issued from metagenomes still constitutes a goal for enzyme discovery. This can be addressed by functional screening of metagenomic libraries, in order to retrieve genes of interest. Numerous studies have provided conclusive evidence on the potential of such an approach for the identification of novel glycoside-hydrolases from various ecosystems such as soil (Rondon et al. 2000; Richardson et al. 2002; Voget et al. 2003; Pang et al. 2009), lakes (Rees et al. 2003), hot springs (Tang et al. 2006, 2008), rumen (Ferrer et al. 2005; Guo et al. 2008; Liu et al. 2008; Duan et al. 2009), rabbit (Feng et al. 2007), and insect guts (Brennan et al. 2004; for review, see Ferrer et al. 2009; Li et al. 2009; Simon and Daniel 2009; Uchiyama and Miyazaki 2009). In all cases, the identification of the gene responsible for the screened activity was carried out by sequencing only a few kilobases of metagenomic DNA. Collectively these studies have established an experimental proof of function for 35 glycoside hydrolases (from eight families) issued from metagenomes (data from the CAZy database; http://www.cazy.org/), a number that is very small considering the known CAZy diversity. Here, we examined the potential of high-throughput functional screening of large insert libraries to guide in-depth pyrosequencing of specific regions of the human gut metagenome that encode the enzymatic machinery involved in dietary fiber catabolism.     

Country-specific antibiotic use practices impact the human gut resistome

Despite increasing concerns over inappropriate use of antibiotics in medicine and food production, population-level resistance transfer into the human gut microbiota has not been demonstrated beyond individual case studies. To determine the “antibiotic resistance potential” for entire microbial communities, we employ metagenomic data and quantify the totality of known resistance genes in each community (its resistome) for 68 classes and subclasses of antibiotics. In 252 fecal metagenomes from three countries, we show that the most abundant resistance determinants are those for antibiotics also used in animals and for antibiotics that have been available longer. Resistance genes are also more abundant in samples from Spain, Italy, and France than from Denmark, the United States, or Japan. Where comparable country-level data on antibiotic use in both humans and animals are available, differences in these statistics match the observed resistance potential differences. The results are robust over time as the antibiotic resistance determinants of individuals persist in the human gut flora for at least a year.When exposed to antibiotic compounds, bacteria evolve resistance mechanisms. These include polymorphisms in antibiotic targets that reduce vulnerability, as well as genes encoding efflux systems, drug modifiers, or proteins that fortify target sites (Wright 2007; Davies and Davies 2010). Resistance determinants can be transferred via mobile genetic elements, such as plasmids, prophages, or transposons, allowing horizontal transfer within and between bacterial species (Davies and Davies 2010), particularly in environments such as the gut microbiome (Salyers et al. 2004; Schjørring and Krogfeldt 2011; Smillie et al. 2011), and have collectively been dubbed the antibiotic resistome (Wright 2007; Marshall and Levy 2011). The transfer of resistance genes into the gut can come from diverse environments, for example, from soil bacteria (Forsberg et al. 2012). Previous studies have explored the pig gut resistome (Looft et al. 2012), as well as that of two human donors (Sommer et al. 2009), but population-scale studies are still lacking. Since antibiotics are widely used in medicine (Goossens et al. 2005) and food production (Barton 2000; Davies and Davies 2010; Marshall and Levy 2011; Aarestrup 2012), understanding the variation of the resistome within the population is crucial.     

The initial peopling of the Americas: A growing number of founding mitochondrial genomes from Beringia

Pan-American mitochondrial DNA (mtDNA) haplogroup C1 has been recently subdivided into three branches, two of which (C1b and C1c) are characterized by ages and geographical distributions that are indicative of an early arrival from Beringia with Paleo-Indians. In contrast, the estimated ages of C1d—the third subset of C1—looked too young to fit the above scenario. To define the origin of this enigmatic C1 branch, we completely sequenced 63 C1d mitochondrial genomes from a wide range of geographically diverse, mixed, and indigenous American populations. The revised phylogeny not only brings the age of C1d within the range of that of its two sister clades, but reveals that there were two C1d founder genomes for Paleo-Indians. Thus, the recognized maternal founding lineages of Native Americans are at least 15, indicating that the overall number of Beringian or Asian founder mitochondrial genomes will probably increase extensively when all Native American haplogroups reach the same level of phylogenetic and genomic resolution as obtained here for C1d.While debate is still ongoing among scientists from several disciplines regarding the number of migratory events, their timing, and entry routes into the Americas (Wallace and Torroni 1992; Torroni et al. 1993; Forster et al. 1996; Kaufman and Golla 2000; Goebel et al. 2003, 2008; Schurr and Sherry 2004; Wang et al. 2007; Waters and Stafford 2007; Dillehay et al. 2008; Gilbert et al. 2008a; O''Rourke and Raff 2010), the general consensus is that modern Native American populations ultimately trace their gene pool to Asian groups who colonized northeast Siberia, including parts of Beringia, prior to the last glacial period. These ancestral population(s) probably retreated into refugial areas during the Last Glacial Maximum (LGM), where their genetic variation was reshaped by drift. Thus, pre-LGM haplotypes of Asian ancestry were differently preserved and lost in Beringian enclaves, but at the same time, novel haplotypes and alleles arose in situ due to new mutations, often becoming predominant because of major founder events (Tamm et al. 2007; Achilli et al. 2008; Bourgeois et al. 2009; Perego et al. 2009; Schroeder et al. 2009). The scenario of a temporally important differentiation stage in Beringia explains the predominance in Native Americans of private alleles and haplogroups such as the autosomal 9-repeat at microsatellite locus D9S1120 (Phillips et al. 2008; Schroeder et al. 2009), the Y chromosome haplogroup Q1a3a-M3 (Bortolini et al. 2003; Karafet et al. 2008; Rasmussen et al. 2010), and the pan-American mtDNA haplogroups A2, B2, C1b, C1c, C1d, D1, and D4h3a (Tamm et al. 2007; Achilli et al. 2008; Fagundes et al. 2008; Perego et al. 2009).In the millennia after the initial Paleo-Indian migrations, other groups from Beringia or eastern Siberia expanded into North America. If the gene pool of the source population(s) had in the meantime partially changed, not only because of drift, but also due to the admixture with population groups newly arrived from regions located west of Beringia, this would have resulted in the entry of additional Asian lineages into North America. This scenario, sometimes invoked to explain the presence of certain mtDNA haplogroups such as A2a, A2b, D2a, D3, and X2a only in populations of northern North America (Torroni et al. 1992; Brown et al. 1998; Schurr and Sherry 2004; Helgason et al. 2006; Achilli et al. 2008; Gilbert et al. 2008b; Perego et al. 2009), has recently received support from nuclear and morphometric data showing that some native groups from northern North America harbor stronger genetic similarities with some eastern Siberian groups than with Native American groups located more in the South (González-José et al. 2008; Bourgeois et al. 2009; Wang et al. 2009; Rasmussen et al. 2010).As for the pan-American mtDNA haplogroups, when analyzed at the highest level of molecular resolution (Bandelt et al. 2003; Tamm et al. 2007; Fagundes et al. 2008; Perego et al. 2009), they all reveal, with the exception of C1d, entry times of 15–18 thousand years ago (kya), which are suggestive of a (quasi) concomitant post-LGM arrival from Beringia with early Paleo-Indians. A similar entry time is also shown for haplogroup X2a, whose restricted geographical distribution in northern North America appears to be due not to a later arrival, but to its entry route through the ice-free corridor (Perego et al. 2009). Despite its continent-wide distribution, C1d was hitherto characterized by an expansion time of only 7.6–9.7 ky (Perego et al. 2009). This major discrepancy has been attributed to a poor and possibly biased representation of complete C1d mtDNA sequences (only 10) in the available data sets (Achilli et al. 2008; Malhi et al. 2010). To clarify the issue of the age of haplogroup C1d and its role as a founding Paleo-Indian lineage, we sequenced and analyzed 63 C1d mtDNAs from populations distributed over the entire geographical range of the haplogroup.     

Widespread somatic L1 retrotransposition occurs early during gastrointestinal cancer evolution

Somatic L1 retrotransposition events have been shown to occur in epithelial cancers. Here, we attempted to determine how early somatic L1 insertions occurred during the development of gastrointestinal (GI) cancers. Using L1-targeted resequencing (L1-seq), we studied different stages of four colorectal cancers arising from colonic polyps, seven pancreatic carcinomas, as well as seven gastric cancers. Surprisingly, we found somatic L1 insertions not only in all cancer types and metastases but also in colonic adenomas, well-known cancer precursors. Some insertions were also present in low quantities in normal GI tissues, occasionally caught in the act of being clonally fixed in the adjacent tumors. Insertions in adenomas and cancers numbered in the hundreds, and many were present in multiple tumor sections, implying clonal distribution. Our results demonstrate that extensive somatic insertional mutagenesis occurs very early during the development of GI tumors, probably before dysplastic growth.Somatic mobilization of retroelements in the cancer genome has only recently been established as a widespread mutational phenomenon. In particular, Long INterspersed Element (LINE)-1 (L1)-mediated retrotransposition has been observed mostly in epithelial cancers. Somatic human-specific L1 (L1Hs) insertions are most abundant in these cancers, but L1-mediated Alu, SVA, and processed pseudogene insertions have also been detected (Iskow et al. 2010; Lee et al. 2012; Solyom et al. 2012; Ewing et al. 2013; Shukla et al. 2013; Cooke et al. 2014; Helman et al. 2014; Pitkanen et al. 2014; Tubio et al. 2014). L1s are autonomous mobile elements that comprise 17% of the human genome and retrotranspose by a “copy and paste” mechanism via an RNA intermediate. This process can lead to insertional mutagenesis and genetic instability (Goodier and Kazazian 2008). Potentially etiological L1 insertions have been reported in APC (Miki et al. 1992) and PTEN exons (Helman et al. 2014) in colorectal and endometrial cancer, respectively, and insertions of unknown significance have been found in numerous other cancer driver genes in a variety of malignancies (Iskow et al. 2010; Lee et al. 2012; Solyom et al. 2012; Ewing et al. 2013; Shukla et al. 2013; Cooke et al. 2014; Helman et al. 2014; Pitkanen et al. 2014; Tubio et al. 2014; Paterson et al. 2015). In a study of somatic retrotransposition during the evolution of prostate and lung cancer, Tubio et al. (2014) found evidence of insertions occurring during cancer development. Beyond this work, the timing of retrotransposition in cancer development has not been analyzed previously.     

Evolutionary assembly patterns of prokaryotic genomes

Evolutionary innovation must occur in the context of some genomic background, which limits available evolutionary paths. For example, protein evolution by sequence substitution is constrained by epistasis between residues. In prokaryotes, evolutionary innovation frequently happens by macrogenomic events such as horizontal gene transfer (HGT). Previous work has suggested that HGT can be influenced by ancestral genomic content, yet the extent of such gene-level constraints has not yet been systematically characterized. Here, we evaluated the evolutionary impact of such constraints in prokaryotes, using probabilistic ancestral reconstructions from 634 extant prokaryotic genomes and a novel framework for detecting evolutionary constraints on HGT events. We identified 8228 directional dependencies between genes and demonstrated that many such dependencies reflect known functional relationships, including for example, evolutionary dependencies of the photosynthetic enzyme RuBisCO. Modeling all dependencies as a network, we adapted an approach from graph theory to establish chronological precedence in the acquisition of different genomic functions. Specifically, we demonstrated that specific functions tend to be gained sequentially, suggesting that evolution in prokaryotes is governed by functional assembly patterns. Finally, we showed that these dependencies are universal rather than clade-specific and are often sufficient for predicting whether or not a given ancestral genome will acquire specific genes. Combined, our results indicate that evolutionary innovation via HGT is profoundly constrained by epistasis and historical contingency, similar to the evolution of proteins and phenotypic characters, and suggest that the emergence of specific metabolic and pathological phenotypes in prokaryotes can be predictable from current genomes.A fundamental question in evolutionary biology is how present circumstances affect future adaptation and phenotypic change (Gould and Lewontin 1979). Studies of specific proteins, for example, indicate that epistasis between sequence residues limits accessible evolutionary trajectories and thereby renders certain adaptive paths more likely than others (Weinreich et al. 2006; Gong et al. 2013; de Visser and Krug 2014; Harms and Thornton 2014). Similarly, both phenotypic characters (Ord and Summers 2015) and specific genetic adaptations (Conte et al. 2012; Christin et al. 2015) show strong evidence of parallel evolution rather than convergent evolution. That is, a given adaptation is more likely to repeat in closely related organisms than in distantly related ones. This inverse relationship between the repeatability of evolution and taxonomic distance implies a strong effect of lineage-specific contingency on evolution, also potentially mediated by epistasis (Orr 2005).Such observations suggest that genetic adaptation is often highly constrained, and the present state of an evolving system can impact future evolution. Yet, the aforementioned studies are limited to small data sets and specific genetic pathways, and a more principled understanding of the rules by which future evolutionary trajectories are governed by the present state of the system is still lacking. For example, it is not known whether such adaptive constraints are a feature of genome-scale evolution or whether they are limited to finer scales. Moreover, the mechanisms that underlie observed constraints are often completely unknown. Addressing these questions is clearly valuable for obtaining a more complete theory of evolutionary biology, but more pressingly, is essential for tackling a variety of practical concerns, including our ability to combat evolving infectious diseases or engineer complex biological systems.Here, we address this challenge by analyzing horizontal gene transfer (HGT) in prokaryotes. HGT is an ideal system to systematically study genome-wide evolutionary constraints because it involves gene-level innovation, occurs at very high rates relative to sequence substitution (Nowell et al. 2014; Puigbò et al. 2014), and is a principal source of evolutionary novelty in prokaryotes (Gogarten et al. 2002; Jain et al. 2003; Lerat et al. 2005; Puigbò et al. 2014). Clearly, many or most acquired genes are rapidly lost due to fitness costs (van Passel et al. 2008; Baltrus 2013; Soucy et al. 2015), indicating that genes retained in the long term are likely to provide a selective advantage. Moreover, not all genes are equally transferrable (Jain et al. 1999; Sorek et al. 2007; Cohen et al. 2011), and not all species are equally receptive to the same genes (Smillie et al. 2011; Soucy et al. 2015). However, differences in HGT among species have been attributed not only to ecology (Smillie et al. 2011) or to phylogenetic constraints (Popa et al. 2011; Nowell et al. 2014), but also to interactions with the host genome (Jain et al. 1999; Cohen et al. 2011; Popa et al. 2011). Indeed, studies involving single genes or single species support the influence of genome content on the acquisition and retention of transferred genes (Pal et al. 2005; Sorek et al. 2007; Iwasaki and Takagi 2009; Chen et al. 2011; Press et al. 2013; Johnson and Grossman 2014). For example, it has been demonstrated that the presence of specific genes facilitates integration of others into genetic networks (Chen et al. 2011), and genes are more commonly gained in genomes already containing metabolic genes in the same pathway (Pal et al. 2005; Iwasaki and Takagi 2009). However, to date, a systematic, large-scale analysis of such dependencies has not been presented. In this paper, we therefore set out to characterize a comprehensive collection of genome-wide HGT-based dependencies among prokaryotic genes, analyze the obtained set of epistatic interactions, and identify patterns in the evolution of prokaryotic genomes.