Key features of the genetic architecture and evolution of host-microbe interactions revealed by high-resolution genetic mapping of the mucosa-associated gut microbiome in hybrid mice
Abstract
Determining the forces that shape diversity in host-associated bacterial communities is critical to understanding the evolution and maintenance of metaorganisms. To gain deeper understanding of the role of host genetics in shaping gut microbial traits, we employed a powerful genetic mapping approach using inbred lines derived from the hybrid zone of two incipient house mouse species. Furthermore, we uniquely performed our analysis on microbial traits measured at the gut mucosal interface, which is in more direct contact with host cells and the immune system. Several mucosa-associated bacterial taxa have high heritability estimates, and interestingly, 16S rRNA transcript-based heritability estimates are positively correlated with cospeciation rate estimates. Genome-wide association mapping identifies 428 loci influencing 120 taxa, with narrow genomic intervals pinpointing promising candidate genes and pathways. Importantly, we identified an enrichment of candidate genes associated with several human diseases, including inflammatory bowel disease, and functional categories including innate immunity and G-protein-coupled receptors. These results highlight key features of the genetic architecture of mammalian host-microbe interactions and how they diverge as new species form.
Editor's evaluation
This paper uses inbred hybrid mouse lines to estimate the heritability of the mucosa-associated microbiome and map variants in the mouse genome that are associated with the composition of the microbiome. The findings are of broad interest to microbiome researchers and improve on knowledge in the field, as the mapping design facilitates the identification of narrow association intervals and points to a novel correlation between heritability and cospeciation rates. The manuscript provides useful information about the approach to heritability estimation, allowing the results to be more readily placed in context. Congratulations on this important contribution to the literature.
https://doi.org/10.7554/eLife.75419.sa0eLife digest
The digestive system, particularly the large intestine, hosts many types of bacteria which together form the gut microbiome. The exact makeup of different bacterial species is specific to an individual, but microbiomes are often more similar between related individuals, and more generally, across related species. Whether this is because individuals share similar environments or similar genetic backgrounds remains unclear. These two factors can be disentangled by breeding different animal lineages – which have different genetic backgrounds while belonging to the same species – and then raising the progeny in the same environment.
To investigate this question, Doms et al. studied the genes and microbiomes of mice resulting from breeding strains from multiple locations in a natural hybrid zone between different subspecies. The experiments showed that 428 genetic regions affected the makeup of the microbiome, many of which were known to be associated with human diseases. Further analysis revealed 79 genes that were particularly interesting, as they were involved in recognition and communication with bacteria. These results show how the influence of the host genome on microbiome composition becomes more specialized as animals evolve.
Overall, the work by Doms et al. helps to pinpoint the genes that impact the microbiome; this knowledge could be helpful to examine how these interactions contribute to the emergence of conditions such as diabetes or inflammatory bowel disease, which are linked to perturbations in gut bacteria.
Introduction
The recent widespread recognition of the gut microbiome’s importance to host health and fitness represents a critical advancement of biomedicine. Host phenotypes affected by the gut microbiome are documented in humans (Ley et al., 2006; Turnbaugh et al., 2009; Lynch and Pedersen, 2016), laboratory animals (Bäckhed et al., 2004; Turnbaugh et al., 2008; Rolig et al., 2015; Rosshart et al., 2017; Gould et al., 2018), and wild populations (Suzuki, 2017; Roth et al., 2019; Suzuki et al., 2020a; Hua et al., 2020), and include critical traits such as aiding digestion and energy uptake (Rowland et al., 2018), and the development and regulation of the immune system (Davenport, 2020).
Despite the importance of the gut microbiome, community composition varies significantly among host species, populations, and individuals (Benson et al., 2010; Yatsunenko et al., 2012; Brooks et al., 2016; Rehman et al., 2016; Amato et al., 2019). While a portion of this variation is expected to be selectively neutral, alterations of the gut microbiome are on the one hand linked to numerous human diseases (Carding et al., 2015; Lynch and Pedersen, 2016), including diabetes (Qin, 2012), inflammatory bowel disease (IBD) (Ott et al., 2004; Gevers et al., 2014), and mental disorders (Clapp et al., 2017). On the other hand, there is evidence that the gut microbiome can play an important role in adaptation on both recent (Hehemann et al., 2010; Suzuki and Ley, 2020b) and ancient evolutionary timescales (Rausch et al., 2019). Collectively, these phenomena suggest that it would be evolutionarily advantageous for hosts to influence their microbiome.
An intriguing observation made in comparative microbiome research in the last decade is that the pattern of diversification among gut microbiomes appears to mirror host phylogeny (Ochman et al., 2010). This phenomenon, coined ‘phylosymbiosis’ (Brucker and Bordenstein, 2012a; Brucker and Bordenstein, 2012b; Lim and Bordenstein, 2020), is documented in a number of diverse host taxa (Brooks et al., 2016) and also extends to the level of the phageome (Gogarten et al., 2021). Several non-mutually exclusive hypotheses are proposed to explain phylosymbiosis (Moran and Sloan, 2015). However, it is likely that vertical inheritance is important for at least a subset of taxa, as signatures of cospeciation/codiversification are present among numerous mammalian associated gut microbes (Moeller et al., 2016; Groussin et al., 2017; Moeller et al., 2019), which could also set the stage for potential coevolutionary processes. Importantly, experiments involving interspecific fecal microbiota transplants indeed provide evidence of host adaptation to their conspecific microbial communities (Brooks et al., 2016; Moeller et al., 2019). Furthermore, cospeciating taxa were observed to be significantly enriched among the bacterial species depleted in early onset IBD, an immune-related disorder, suggesting a greater evolved dependency on such taxa (Papa et al., 2012; Groussin et al., 2017). However, the nature of genetic changes involving host-microbe interactions that take place as new host species diverge remains underexplored.
House mice are an excellent model system for evolutionary microbiome research, as studies of both natural populations and laboratory experiments are possible (Suzuki, 2017; Suzuki et al., 2019). In particular, the house mouse species complex is composed of subspecies that hybridize in nature, enabling the potential early stages of codiversification to be studied. We previously analysed the gut microbiome across the Central European hybrid zone of Mus musculus musculus and Mus musculus domesticus (Wang et al., 2015), which share a common ancestor ~0.5 million years ago (Geraldes et al., 2008). Importantly, transgressive phenotypes (i.e. exceeding or falling short of parental values) among gut microbial traits as well as increased intestinal histopathology scores were common in hybrids, suggesting that the genetic basis of host control over microbes has diverged (Wang et al., 2015). The same study performed an F2 cross between wild-derived inbred strains of M. m. domesticus and M. m. musculus and identified 14 quantitative trait loci (QTL) influencing 29 microbial traits. However, like classical laboratory mice, these strains had a history of rederivation and reconstitution of their gut microbiome, thus leading to deviations from the native microbial populations found in nature (Rosshart et al., 2017; Org and Lusis, 2018), and the genomic intervals were too large to identify individual genes.
In this study, we employed a powerful genetic mapping approach using inbred lines directly derived from the M. m. musculus–M. m. domesticus hybrid zone, and further focus on the mucosa-associated microbiota due to its more direct interaction with host cells (Fukata and Arditi, 2013; Chu and Mazmanian, 2013), distinct functions compared to the luminal microbiota (Wang et al., 2010; Vaga et al., 2020), and greater dependence on host genetics (Spor et al., 2011; Linnenbrink et al., 2013). Previous mapping studies using hybrids raised in a laboratory environment showed that high mapping resolution is possible due to the hundreds of generations of natural admixture between parental genomes in the hybrid zone (Turner and Harr, 2014; Pallares et al., 2014; Škrabar et al., 2018). Accordingly, we here identify 428 loci contributing to variation in 120 taxa, whose narrow genomic intervals (median <2 Mb) enable many individual candidate genes and pathways to be pinpointed. We identify a high proportion of bacterial taxa with significant heritability estimates and find that bacterial phenotyping based on 16S rRNA transcript compared to gene copy-based profiling yields an even higher proportion. Furthermore, these heritability estimates also significantly positively correlate with cospeciation rate estimates, suggesting a more extensive host genetic architecture for cospeciating taxa. Finally, we identify numerous enriched functional pathways, whose role in host-microbe interactions may be particularly important as new species form.
Results
Microbial community composition
To obtain microbial traits for genetic mapping in the G2 mapping population, we sequenced the 16S rRNA gene from caecal mucosa samples of 320 hybrid male mice based on DNA and RNA (cDNA), which reflect bacterial cell number and activity, respectively. After applying quality filtering and subsampling 10,000 reads per sample, we identified a total of 4684 amplicon sequence variants (ASVs). For further analyses, we established a ‘core microbiome’ (defined in Materials and methods), such that analyses were limited to those taxa common and abundant enough to reveal potential genetic signal. The core microbiome is composed of four phyla, five classes, five orders, 11 families, 27 genera, and 90 ASVs for RNA, and four phyla, five classes, six orders, 12 families, 28 genera and 46 ASVs for DNA. A combined total of 98 unique ASVs belong to the core, of which 38 were shared between DNA and RNA (Figure 1—figure supplement 1). The most abundant genus in our core microbiome is Helicobacter (Figure 1—figure supplement 2B), consistent with a previous study of the wild hybrid M. m. musculus/M. m. domesticus mucosa-associated microbiome (Wang et al., 2015).
Importantly, inspection of the taxonomic profiles of the mapping population confirms that key features of the native mouse microbiome were retained in our wild-derived lines, despite multiple generations of breeding in the laboratory. For example, the number, identity, and relative proportions of major bacterial orders are more similar to wild-caught mice than to rederived classical laboratory strains (Figure 1—figure supplement 2A – Order-level bar plot; Fig S4 from Rosshart et al., 2017). This is consistent with a previous study which showed that microbiomes of wild-derived strains maintained their distinctiveness over 10 generations of breeding in the laboratory (Moeller et al., 2018).
Heritability
Next, we estimated the narrow-sense heritability (h2) of bacterial traits using lme4QTL (Ziyatdinov et al., 2018). Of the 153 total core taxa, we identified 21 taxa for DNA and 30 taxa for RNA with significant heritability estimates (PRLRT<0.05, Restricted Likelihood Ratio test), with estimates ranging between 39 and 83% (Figure 1A–B and Supplementary file 1). The top values for bacterial abundances are similar to heritability estimates for body weight (87%) and body length (67%). ASV97 (genus Oscillibacter) followed by the genus Paraprevotella, and ASV7 (genus Paraprevotella) showed the highest heritability among DNA-based traits (80.8, 78.6, and 77.4%, respectively; Figure 1A), while ASV97 (genus Oscillibacter), followed by ASV36 (genus Oscillibacter), and ASV135 (unclassified order Bacteroidales) had the highest heritability among RNA-based traits (83.0, 80.2, and 79.3%, respectively; Figure 1B). The heritability estimates for DNA- and RNA-based measurements of the same taxa are significantly correlated (Spearman’s rho = 0.53, p=3.1 × 10–12, Figure 1—figure supplement 3).
Next, we estimated the ‘chip’ heritability (CH), i.e., the percentage of phenotypic variance in bacterial traits explained by genotyped SNPs (32,625 SNPs; Zhou et al., 2013). We find 23 DNA-based traits and 27 RNA-based traits with significant chip heritability estimates, ranging between 50.0 and 15.9% (Figure 1A, B and Supplementary file 1).
We compared these heritability estimates to estimates from previous studies in other mammals (Supplementary file 1), including mice (O’Connor et al., 2014; Org et al., 2015), humans (Davenport et al., 2015; Goodrich et al., 2016; Turpin et al., 2016; Xu et al., 2020; Ishida et al., 2020; Hughes et al., 2020; Kurilshikov et al., 2021), pigs (Chen et al., 2018), and primates (Grieneisen et al., 2021). DNA-based heritability estimates are positively correlated with DNA-based heritability estimates from male mice (Spearman’s rho = 0.60, p=0.049, n=11; Org et al., 2015; Figure 1—figure supplement 5A) and with DNA-based heritability estimates from one human study (Spearman’s rho = 0.38, p=0.049, n=28; Turpin et al., 2016; Figure 1—figure supplement 5B).
Heritability estimates are correlated with predicted cospeciation rates
In an important meta-analysis of the gut microbiome across diverse mammalian taxa, Groussin et al., 2017 estimated cospeciation rates of individual bacterial taxa by measuring the congruence of host and bacteria phylogenetic trees relative to the number of host-swap events. We reasoned that taxa with higher cospeciation rates might also demonstrate higher heritability, as these more intimate evolutionary relationships would provide a greater opportunity for genetic aspects to evolve. Intriguingly, we observe a significant positive correlation for RNA-based traits (h2, p=0.037, Spearman’s rho = 0.47; CH, p=0.012, Spearman’s rho = 0.55; Figure 1D; Figure 1—figure supplement 4D), but not for DNA-based traits (h2, Spearman’s rho = −0.062, p=0.80; CH, Spearman’s rho = −0.091, p=0.70; Figure 1C; Figure 1—figure supplement 4C) for both narrow-sense heritability and chip heritability estimates. To evaluate whether these results may be confounded by taxon abundance, we further used a multiple linear regression model incorporating both taxon abundance and cospeciation rate. The overall regression model was not significant (R2=0.27, F(2,17)=3.202, p=0.067). Furthermore, the cospeciation rate significantly predicts the heritability estimate (p=0.022), while the median abundance does not (p=0.92). Thus, heritability estimates and cospeciation rates are associated independent of taxon abundance. These results support the notion that cospeciating taxa evolved a greater dependency on host genes, and further suggest that bacterial activity may better reflect the underlying biological interactions.
Genetic mapping of host loci determining microbiome composition
Next, we performed genome-wide association mapping of the relative abundances of core taxa, in addition to two alpha-diversity measures (Shannon and Chao1 indices), based on 32,625 SNPs. We used a linear mixed model including both additive and dominance terms in the model to enable the identification of underdominance and overdominance in this hybrid mapping population. We included mating pair and the genotype-based genomic relatedness matrix (GRM) as random effects to control for maternal effects and relatedness, respectively (see Materials and methods). While we found no genome-wide significant associations for alpha diversity at either the DNA or RNA level (p>1.53 × 10–6), a total of 1030 genome-wide significant associations were identified for individual taxa (p<1.53 × 10–6, Supplementary file 2), of which 428 achieved study-wide significance (p<1.29 × 10–8). Apart from the X chromosome, all autosomal chromosomes contained study-wide significant associations (Figure 2). Out of the 153 mapped taxa, 120 had at least one significant association (Table 1). For the remainder of our analyses, we focus on the results using the more stringent study-wide threshold, and combined significant SNPs within 10 Mb into significant regions (Supplementary file 3). The median size of significant regions is 1.91 Mb, which harbour a median of 14 protein-coding genes. On average, we observe five significant mouse genomic regions per bacterial taxon.
Of the significant loci with estimated interval sizes, we find 69 intervals (16.1%) that are smaller than 1 Mb (Supplementary file 4). The smallest interval is only 231 bases and associated with the RNA-based abundance of an unclassified genus belonging to Deltaproteobacteria. It is situated in an intron of the C3 gene, a complement component playing a central role in the activation of the complement system, which modulates inflammation and contributes to antimicrobial activity (Ricklin et al., 2016).
The significant genomic regions and SNPs are displayed in Figure 2A and B, respectively. Individual SNPs were associated with up to 14 taxa, and significant intervals with up to 30 taxa. The SNPs with the lowest p-values were associated with the genus Dorea and two ASVs belonging to Dorea (ASV184 and ASV293; Figure 2—figure supplement 1). At the RNA level this involves two loci: mm10-chr4: 67.07 Mb, where the peak SNP is 13 kb downstream of the closest gene Tlr4 (UNC7414459, p=2.31 × 10–69, additive p=4.48 × 10–118, dominance p=1.37 × 10–111; Figure 2; Figure 2—figure supplement 1), and mm10-chr15: 94.4 Mb, where the peak SNP is found within the Adamts20 gene (UNC26145702, p=4.51 × 10–65, additive p=1.87 × 10–113, dominance p=1.56 × 10–105; Figure 2; Figure 2—figure supplement 1). Interestingly, the Irak4 gene, whose protein product is rapidly recruited after TLR4 activation, is also located 181 kb upstream of Adamts20. The five taxa displaying the most associations were ASV19 (Bacteroides), Dorea, ASV36 (Oscillibacter), ASV35 (Bacteroides), and ASV98 (unclassified Lachnospiraceae) (Figure 2—figure supplement 2).
Ancestry, dominance, and effect sizes
A total of 398 significant SNPs were ancestry informative between M. m. musculus and M. m. domesticus (i.e. represent fixed differences between subspecies). To gain further insight into the genetic architecture of microbial trait abundances, we estimated the degree of dominance at each significant locus using the d/a ratio (Falconer, 1996), where alleles with strictly recessive, additive, and dominant effects have d/a values of –1, 0, and 1, respectively. As half of the SNPs were not ancestry informative (Figure 3A), it was not possible to consistently have a associated with one parent/subspecies, hence we report d/|a| such that it can be interpreted with respect to bacterial abundance. For the vast majority of loci (83.79%), the allele associated with higher abundance is recessive or partially recessive (–1.25<d/|a|<–0.75; Figure 3B). On the basis of the arbitrary cutoffs, we used to classify dominance, only a small proportion of alleles are underdominant (0.23%; d/|a|<–1.25). However, for one-third of the significant SNPs, the heterozygotes display transgressive phenotypes, i.e., mean abundances that are either significantly lower (31% of SNPs) or higher (2% of SNPs) than those of both homozygous genotypes. Interestingly, the domesticus allele was associated with higher bacterial abundance in two-thirds of this subset (33.9 vs. 16.5% musculus allele; Figure 3A).
Next, we estimated phenotypic effect sizes by calculating the percent variance explained (PVE) by the peak SNP of each significant region. Peak SNPs explain between 3 and 64% of the variance in bacterial abundance, with a median effect size of 9.3% (Figure 3C). The combined PVE by the additive effects of all significant markers for each taxon ranged from 0.000018 to 41.6%, with an average of 12% (Figure 1A–B). As expected, the combined additive effects of significant loci are typically much lower than the h2, which is the upper bound as it represents the total additive genetic effect. Interestingly, there are several taxa for which the PVE by additive and dominance effects of all significant SNPs exceeds h2 (e.g. genus Odoribacter and ASV234 [unclassified Ruminococcaceae] for RNA-based traits), indicating there are strong dominance effects. For example, Odoribacter has two significant regions which show overdominance, and ASV234 has seven significant regions which show underdominance.
Functional annotation of candidate genes
Our mapping procedure involves testing for associations between a given microbial trait and a single SNP marker. However, the true genetic architecture is likely more complex. For example, multiple interacting genes in a common host regulatory pathway could influence a given taxonomic group and/or function, the latter of which could also be distributed among unrelated taxa (i.e. functional redundancy). Thus, in order to reveal potential higher-level biological phenomena among the identified loci, we performed pathway analysis to identify interactions and functional categories enriched among the genes in significant intervals. We used STRING (Szklarczyk et al., 2019) to calculate a protein-protein interaction (PPI) network of 925 protein-coding genes nearest to significant SNPs (upstream and/or downstream). A total of 768 genes were represented in the STRING database, and the maximal network is highly significant (STRING PPI enrichment p-value: 2.15×10–14) displaying 668 nodes connected by 1797 edges and an average node degree of 4.68. After retaining only the edges with the highest confidence (interaction score >0.9), this results in one large network with 233 nodes, 692 edges and ten smaller networks (Figure 4).
Next, we functionally annotated clusters using STRING’s functional enrichment plugin. The genes of the largest cluster are part of the G-protein-coupled receptor (GPCR) ligand binding pathway. GPCRs are the largest receptor superfamily and also the largest class of drug targets (Sriram and Insel, 2018). We then calculated the top ten hub proteins from the network based on Maximal Clique Centrality (MCC) algorithm with CytoHubba to predict important nodes that can function as ‘master switches’ (Figure 4—figure supplement 1). The top ten proteins contributing to the PPI network were GNG12, MCHR1, NMUR2, PROK2, OXTR, XCR1, TACR3, CHRM3, PTGFR, and C3, which are all involved in the GPCR signalling pathway.
Furthermore, we performed enrichment analysis on the 925 genes nearest to significant SNPs using the clusterprofiler R package. We found 14 KEGG pathways to be overrepresented: circadian entrainment, oxytocin signalling pathway, axon guidance, calcium signalling, cAMP signalling, cortisol synthesis and secretion, cushing syndrome, gastric acid secretion, glutamatergic synapse, mucin type O-glycan biosynthesis, inflammatory mediator regulation of TRP channels, PD-L1 expression and the PD-1 checkpoint pathway in cancer, tight junction, and the Wnt signalling pathway (Supplementary file 5, Figure 4—figure supplement 2 and Figure 4—figure supplement 3). Finally, genes involved in five human diseases are enriched, among them mental disorders (Figure 4—figure supplement 4).
Finally, due to the observation of a significant enrichment of cospeciating taxa among the bacterial species depleted in early onset IBD (Groussin et al., 2017) and the evidence that IBD is especially associated with a dysbiosis in mucosa-associated communities (Yang et al., 2020a; Daniel et al., 2021), we specifically examined possible overrepresentation of genes involved in IBD (Khan et al., 2021) among the 925 genes neighbouring significant SNPs. We found 14 out of the 289 IBD genes, which was significantly more than expected by chance (10,000 times permuted mean: 2.7, simulated p=00001, Fisher’s exact test; Supplementary file 6). Interestingly, SNPs in 5 out of the 14 genes are associated with ASVs belonging to the genus Oscillibacter, a cospeciating taxon known to decrease during the active state of IBD (Metwaly et al., 2020).
Comparison of significant loci to published gut microbiome mapping studies
Host loci that appear in multiple independent studies are more likely to represent true positive associations and/or less dependent on environmental perturbations. We therefore compiled a list of 648 unique confidence intervals of significant associations with gut bacterial taxa from seven previous mouse QTL studies (Benson et al., 2010; McKnite et al., 2012; Leamy et al., 2014; Wang et al., 2015; Org et al., 2015; Snijders et al., 2016; Kemis et al., 2019) and compared this list to our significance intervals for bacterial taxa at both the DNA and RNA level (341 intervals). Regions larger than 10 Mb were removed from all studies. We found 441 overlapping intervals, which is significantly more than expected by chance (10,000 times permuted mean: 372.8, simulated p=0.005, Fisher’s exact test, see Materials and methods). Several of our smaller significant loci overlapped with larger loci from previous studies and removing this redundancy left 190 significant loci with a median interval size of 0.83 Mb (Figure 5). The most frequently identified locus is located on chromosome 2 169–171 Mb where protein coding genes Gm11011, Znf217, Tshz2, Bcas1, Cyp24a1, Pfdn4, 4930470P17Rik, and Dok5 are situated.
Additionally, we collected genes within genome-wide significant regions reported in seven human microbiome GWAS (mGWAS) (Bonder et al., 2016; Turpin et al., 2016; Goodrich et al., 2016; Wang et al., 2016; Hughes et al., 2020; Rühlemann et al., 2021; Kurilshikov et al., 2021). However, no significant overrepresentation of genes was found within our significance intervals (p=0.16, Fisher’s exact test), nor within our list of genes closest to a significant SNP (p=0.62, Fisher’s exact test).
Proteins differentially expressed in germ-free vs. conventional mice
To further validate our results, we compared the list of genes contained within intervals of our study to a list of differentially expressed proteins between germ-free and conventionally raised mice (Mills et al., 2020). This comparison was made based on the general expectation that if a host gene influences microbial abundance, its own expression would be more likely to change according to differences in the microbiome. Thus, we examined the intersection between genes identified in our study and the proteins identified as highly associated (|π|>1) with the colonisation state of the colon and the small intestine (Mills et al., 2020). Out of the 373 overexpressed or underexpressed proteins according to colonisation status, we find 194 of their coding genes to be among our significant loci, of which 17 are the closest genes to a significant marker (Iyd, Nln, Slc26a3, Slc3a1, Myom2, Nebl, Tent5a, Fxr1, Cbr3, Chrodc1, Nucb2, Arhgef10l, Sucla2, Enpep, Prkcq, Aacs, and Cox7c). This is significantly more than expected by chance (simulated p=0.016, 10,000 permutations, Fisher’s exact test). Furthermore, analysing the PPI with STRING results in a significant network (STRING PPI enrichment <i>P-value = 1.73 × 10–14, and average node degree 2.4, Figure 5—figure supplement 1), with Cyp2c65, Cyp2c55, Cyp2b10, Gpx2, Cth, Eif3k, Eif1, Sucla2, and Rpl17 identified as hub genes (Figure 5—figure supplement 2).
Subsequently, we merged the information from Mills et al., 2020 and the seven previous QTL mapping studies discussed above to further narrow down the most promising candidate genes, and found 30 genes overlapping with our study. Of these 30 genes, six are the closest gene to a significant SNP. These genes are myomesine 2 (Myom2), solute carrier family 3 member 1 (Slc3a1), solute carrier family 26 member 3 (Slc26a3), nebulette (Nebl), carbonyl reductase 3 (Cbr3), and acetoacetyl-coA synthetase (Aacs).
Candidate genes influencing bacterial abundance
To compile a comprehensive set of promising candidate genes, we combined results from network analysis, overlap with previous mouse QTL studies, and differential expression in GF vs conventional mice (see Materials and methods 'Curation of candidate genes’ and ). Next, we used STRING to construct a PPI network with this curated gene set, which is highly significant (STRING PPI enrichment <i>P-value < 1.0 × 10–16, average node degree = 4.85). We identified genes with the highest connectivity and most supporting information (original network see Figure 6—figure supplement 1), resulting in a final set of 79 candidate genes (Figure 6 and Supplementary file 7). The G-protein, GNG12 and the complement component 3 C3, are the proteins with the most edges in the network (30 and 25, respectively), followed by MCHR1, CXCL12, and NMUR2 with each 18 edges. Of these 79 highly connected genes, 35 are associated with bacteria that are either cospeciating (cospeciation rate >0.5; Groussin et al., 2017) and/or have high heritability (>0.5) suggesting a functionally important role for these bacterial taxa.
Discussion
Understanding the forces that shape variation in host-associated bacterial communities within host species is key to understanding the evolution and maintenance of meta-organisms. Although numerous studies in mice and humans demonstrate that host genetics influences gut microbiota composition (McKnite et al., 2012; Leamy et al., 2014; Goodrich et al., 2014; Org et al., 2015; Davenport et al., 2015; Wang et al., 2016; Bonder et al., 2016; Goodrich et al., 2016; Kemis et al., 2019; Suzuki et al., 2019; Ishida et al., 2020; Hughes et al., 2020; Rühlemann et al., 2021), our study is unique in a number of important ways. First, the unique genetic resource of mice collected from a naturally occurring hybrid zone together with their native microbes yielded extremely high mapping resolution and the possibility to uncover ongoing evolutionary processes in nature. Second, our study is the first to perform genetic mapping of 16S rRNA transcripts in the gut environment, which was previously shown to be superior to DNA-based profiling in a genetic mapping study of the skin microbiota (Belheouane et al., 2017). Third, our study is one of the only to specifically examine the mucosa-associated community. It was previously reasoned that the mucosal environment may better reflect host genetic variation (Spor et al., 2011), and evidence for this hypothesis exists in nature (Linnenbrink et al., 2013). Finally, by cross-referencing our results with previous mapping studies and recently available proteomic data from germ-free versus conventional mice, we curated a more reliable list of candidate genes and pathways. Taken together, these results provide unique and unprecedented insight into the genetic basis for host-microbe interactions (Supplementary file 1).
Importantly, by using wild-derived hybrid inbred strains to generate our mapping population, we gained insight into the evolutionary association between hosts and their microbiota at the transition from within species variation to between species divergence. Genetic relatedness in our mapping population significantly correlates with microbiome similarity, supporting a basis for codiversification at the early stages of speciation. A substantial proportion of microbial taxa are heritable, and heritability is correlated with cospeciation rates. This suggests that (i) vertical transmission could enable greater host adaptation to bacteria and/or (ii) the greater number of host genes associated with cospeciating taxa could indicate a greater dependency on the host, such that survival outside a specific host is reduced, making horizontal transmission less likely.
By performing 16S rRNA gene profiling at both the DNA and RNA level, we found that 14% (DNA based) to 20% (RNA based) of bacterial taxa have significant heritability estimates, with values up to 83%. The proportion of heritable taxa and maximum value of heritability estimates are consistent with previous studies in humans and mice ( O’Connor et al., 2014; Org et al., 2015; Hughes et al., 2020). Several factors of our study design likely contribute to high heritability values for some taxa. First, mice were raised in a controlled common environment, and heritability estimates in other mammals were shown to be contingent on the environment (Grieneisen et al., 2021). Second, bacterial communities were sampled from cecal tissue (mucosa) instead of lumen/faeces (Linnenbrink et al., 2013). Third, genetic variation in our mapping population was higher than in typical mapping studies due to subspecies differences. Finally, not all studies discriminate between narrow-sense heritability estimates and chip heritability estimates. Chip heritability represents only the variance explained by genotyped SNPs, thus estimates are lower than narrow-sense heritability estimates, which include all additive genetic effects (Zhou et al., 2013). Chip heritability is particularly useful for phenotype prediction, for example, in the context of animal breeding and human disease, whereas narrow-sense heritability is useful for understanding genetic influence more broadly.
Notably, heritability estimates for RNA-based traits are significantly correlated with previously reported cospeciation rates in mammals (Groussin et al., 2017). This pattern, as well as the higher proportion of heritable taxa for RNA-based traits, suggest that host genetic effects are more strongly reflected by bacterial activity than cell number.
We found a total of 172 and 305 unique significant loci for DNA-based and RNA-based bacterial abundance, respectively, passing the conservative study-wide significance threshold. Taxa had a median of four significant loci, suggesting a complex and polygenic genetic architecture affecting bacterial abundances. We identify a higher number of loci in comparison to previous QTL and GWAS studies in mice (Benson et al., 2010; McKnite et al., 2012; Leamy et al., 2014; Wang et al., 2015; Org et al., 2015; Snijders et al., 2016; Kemis et al., 2019), which may be due to a number of factors. The parental strains of our study were never subjected to rederivation and subsequent reconstitution of their microbiota, and natural mouse gut microbiota are more variable than the artificial microbiota of laboratory strains (Kohl and Dearing, 2014; Weldon et al., 2015; Suzuki, 2017; Rosshart et al., 2017). Furthermore, as noted above, our mapping population harbours both within-subspecies and between-subspecies genetic variation. We crossed incipient species sharing a common ancestor ~0.5 million years ago, hence we may also capture the effects of mutations that fixed rapidly between subspecies due to strong selection, which are typically not variable within species (Walsh and Lynch, 1998; Barton and Keightley, 2002).
Importantly, our results also help to describe general features of the genetic architecture of bacterial taxon activity. For the majority of loci, the allele associated with lower relative abundance of the bacterial taxon was (partially) dominant. This suggests there is strong purifying selection against a high abundance of any particular taxon, which may help ensure high alpha diversity. For several bacterial taxa, the PVE by additive and dominance effects of significant SNPs exceeds the narrow-sense heritability, and the heterozygotes of one-third of significant SNPs displayed transgressive phenotypes. This is consistent with previous studies of hybrids (Turner et al., 2012; Turner and Harr, 2014; Wang et al., 2015), for example, wild-caught hybrids showed broadly transgressive gut microbiome phenotypes. This pattern can be explained by overdominance or underdominance, or by epistasis (Rieseberg et al., 1999). A curious observation is that domesticus alleles are associated with higher relative bacterial abundances twice as often as musculus alleles (Figure 3A). Although the biological explanation for this pattern is not discernible from our current results, future work incorporating quantitative profiling of bacterial load in hybrid and unadmixed musculus and domesticus individuals may help explain this phenomenon.
Notably, many loci significantly associated with bacterial abundance in this study were implicated in previous studies (Figure 5). For example, chromosome 2 169–171 Mb is associated with ASV23 (Eisenbergiella), Eisenbergiella and ASV32 (unclassified Lachnospiraceae) in this study, and overlaps with significant loci from three previous studies (Leamy et al., 2014; Snijders et al., 2016; Kemis et al., 2019). This region contains eight protein-coding genes: Gm11011, Znf217, Tshz2, Bcas1, Cyp24a1, Pfdn4, 4930470P17Rik, and Dok5. Another hotspot is on chromosome 5 101–103 Mb. This locus is significantly associated with four taxa in this study (Prevotellaceae, Paraprevotella, ASV7 genus Paraprevotella and Acetatifactor) and overlaps with associations for Clostridiales, Clostridiaceae, Lachnospiraceae, and Deferribacteriaceae (Snijders et al., 2016). Protein-coding genes in this region are: Nkx6-1, Cds1, Wdfy3, Arhgap24, and Mapk10. As previous studies were based on rederived mouse strains, identifying significant overlap in the identification of host loci suggests that some of the same genes and/or mechanisms influencing major members of gut microbial communities are conserved even in the face of community ‘reset’ in the context of rederivation. The identity of the taxa is however not always the same, which suggests that functional redundancy may contribute to these observations, if, for example, several bacterial taxa fulfill the same function within the gut microbiome (Moya and Ferrer, 2016; Tian et al., 2020).
A limitation of the current and previous genetic studies is that the phenotypes used for mapping are based on relative abundance rather than absolute, quantitative estimates. Specifically, an increase in a given taxon’s abundance will necessarily lead to a decrease in abundance among all other taxa, which can lead to a number of potential biases (Kathagen et al., 2017; Barlow et al., 2020). Thus, future studies incorporating absolute abundance estimates may improve the detection of host-microbe interactions.
Nevertheless, our results do display overlap with studies based on other independent methods, such as a list of proteins differentially expressed in the intestine of germ-free mice compared to conventionally raised mice (Mills et al., 2020). Furthermore, by analysing the functions of the genes closest to significant SNPs, we found that 12 of the 14 significantly enriched KEGG pathways were shown to be related to interactions with bacteria (Fonken et al., 2010 Thaiss et al., 2014; Neumann et al., 2014; Thaiss et al., 2015a; Thaiss et al., 2015b; Castoldi, 2015; Erdman and Poutahidis, 2016; Thaiss et al., 2016; Deaver et al., 2018; Wu et al., 2018; Peng et al., 2020; Nagpal et al., 2020; Hollander and Kaunitz, 2019; Supplementary file 5).
To improve the robustness of our results, we combined multiple lines of evidence to prioritise candidates, resulting in a network of 79 genes (Supplementary file 7). At the centre of this network is a set of 22 proteins involved in G-protein coupled receptor signalling (Figure 6, red nodes). MCHR1, NMUR2, and TACR3 (Figure 6, yellow) are known to regulate feeding behaviour (Saito et al., 1999; Cardoso et al., 2012; Smith et al., 2019), and CHRM3 to control digestion (Gautam et al., 2006; Tanahashi et al., 2009). Gut microbes can produce GPCR agonists to elicit host cellular responses (Cohen et al., 2017; Colosimo et al., 2019; Chen et al., 2019; Pandey et al., 2019). Thus, GPCRs may be key modulators of communication between the gut microbiota and host. Another interesting group of genes are those responding to nutrient levels (Bmp7, Cd40, Aacs, Gclc, Nmur2, Cyp24a1, Adcyap1, Serpinc1, and Wnt11) (Sethi and Vidal-Puig, 2008; Peier et al., 2009; Townsend et al., 2012; Yi and Bishop, 2015; Shi and Tu, 2015; Toderici et al., 2016; Yasuda et al., 2021; Gastelum et al., 2021), as gut microbiota affect host nutrient uptake (Chung et al., 2018). In addition, CYP24A1, BMP7, and CD40 respond to vitamin D. Previous studies identified vitamin D/the vitamin D receptor to play a role in modulating the gut microbiota (Wang et al., 2016; Malaguarnera, 2020; Yang et al., 2020b; Singh et al., 2020), and CD40 is known to induce a vitamin D dependent antimicrobial response through IFN-γ activation (Klug-Micu et al., 2013).
Another important category of candidate genes are those involved in immunity. Our most significant SNP was situated downstream of the Tlr4 gene and was associated with the genus Dorea and several Dorea species. Dorea is a known short chain fatty acid producer (Taras, 2002; Reichardt et al., 2018) and interacts with tight junction proteins Claudin-2 and Occludin (Alhasson et al., 2017). Tlr4 is a member of the Toll-like receptor family, and has been linked with obesity, inflammation, and changes in the gut microbiota (Velloso et al., 2015). These combined results reflect an important role for Dorea in fatty acid harvesting and intestinal barrier integrity, both of which could act systemically to activate TLR4 and to promote metabolic inflammation (Cani et al., 2008; Delzenne et al., 2011; Nicholson et al., 2012). Moreover, the SNP with the second lowest p-value was associated with the same taxa and situated 181 kb upstream of Irak4. IRAK4 is rapidly recruited after TLR4 activation to enable downstream activation of the NFκB immune pathway. Irak4 has previously been associated with a change in bacterial abundance using inbred mice (McKnite et al., 2012; Org et al., 2015).
Finally, we identified notable links between candidate genes and five human diseases (mental disorders, blood pressure finding, systemic arterial pressure, substance-related disorders, and atrial septal deficits; Figure 4—figure supplement 4). The connection to mental disorders is intriguing as involvement of the gut microbiota is suspected (Kelly et al., 2015; Foster et al., 2017; Cox and Weiner, 2018; Chen et al., 2019; Sarkar et al., 2020; Parker et al., 2020; Flux and Lowry, 2020). Taken together with our finding of an enriched set of GPCRs, this highlights the importance of host-microbial interplay along the gut-brain axis. Moreover, we also identify a significant overrepresentation of IBD genes (Khan et al., 2021) among the 925 genes nearest to significant SNPs (Supplementary file 6). Interestingly, SNPs in 5 out of 14 genes are associated with ASVs belonging to the genus Oscillibacter, a highly cospeciating taxon known to decrease during the active state of IBD (Metwaly et al., 2020).
In summary, our study provides a number of novel insights into the importance of host genetic variation in shaping the gut microbiome, in particular for cospeciating bacterial taxa. These findings provide an exciting foundation for future studies of the precise mechanisms underlying host-gut microbiota interactions in the mammalian gut and should encourage future genetic mapping studies that extend analyses to the functional metagenomic sequence level.
Materials and methods
Intercross design
Request a detailed protocolWe generated a mapping population using partially inbred strains derived from mice captured in the M. musculus–M. m. domesticus hybrid zone around Freising, Germany, in 2008 (Turner et al., 2012). Originally, four breeding stocks were derived from 8 to 9 ancestors captured from one (FS, HA, TU) or two sampling sites (HO), and maintained with four breeding pairs per generation using the HAN-rotation out-breeding scheme (Rapp, 1972). Eight inbred lines (two per breeding stock) were generated by brother/sister mating of the 8th generation lab-bred mice. We set up the cross when lines were at the 5th–9th generation of brother-sister meeting, with inbreeding coefficients of >82%.
We used power calculations to estimate the optimal cross design and sample size needed. We first set up eight G1 crosses, each with one predominantly domesticus line (FS, HO – hybrid index <50%; see below) and one predominantly musculus line (HA, TU – hybrid index >50%); each line was represented as a dam in one cross and sire in another (Figure 1—figure supplement 6). One line, FS5, had a higher hybrid index than expected, suggesting there was a misidentification during breeding (see genotyping below). Next, we set up G2 crosses in eight combinations (sub-crosses), such that each G2 individual has one grandparent from each of the initial four breeding stocks. We included 40 males from each sub-cross in the mapping population resulting in a total of 320 mice. Mice were kept together with male littermates until moved to single cages 1 week prior to sacrifice.
This study was performed according to approved animal protocols and institutional guidelines of the Max Planck Institute. Mice were maintained and handled in accordance with FELASA guidelines and German animal welfare law (Tierschutzgesetz § 11, permit from Veterinäramt Kreis Plön: 1401–144/PLÖ–004697).
Sample collection
Request a detailed protocolMice were kept in the same room and caged together with littermates after weaning before being separated into single cages 1 week prior to dissection (to minimize effects of social dominance on fertility traits for a related study). Mice were sacrificed at 91±5 days by CO2 asphyxiation. We recorded body weight, body length and tail length, and collected ear tissue for genotyping. The caecum was removed and gently separated from its contents through bisection and immersion in RNAlater (Thermo Fisher Scientific, Schwerte, Germany). After overnight storage in RNAlater at 4°C, the RNAlater was removed and tissue stored at –20°C.
DNA extraction and sequencing
Request a detailed protocolWe simultaneously extracted DNA and RNA from caecum tissue samples using Qiagen (Hilden, Germany) Allprep DNA/RNA 96-well kits. All samples were extracted together in the same extraction round and hence timepoint. We followed the manufacturer’s protocol, with the addition of an initial bead beating step using Lysing matrix E tubes (MP Biomedical, Eschwege) to increase cell lysis. We used caecum tissue because host genetics has a greater influence on the microbiota at this mucosal site than on the lumen contents (Linnenbrink et al., 2013). We performed reverse transcription of RNA with high-capacity cDNA transcription kits from Applied Biosystems (Darmstadt, Germany). We amplified the V1–V2 hypervariable region of the 16S rRNA gene using barcoded primers (27F-338R) with fused MiSeq adapters and heterogeneity spacers following the description in Rausch et al., 2016 and sequenced amplicons with 250bp paired-reads on the Illumina MiSeq platform. Individual sequencing libraries were prepared in parallel for all DNA and RNA samples, respectively, resulting in one MiSeq library each. Accordingly, all individual 16S rRNA gene profiles within a given DNA- or RNA-based mapping analysis were generated by a single sequencing run. Thus, only direct comparisons between DNA- and RNA-based traits could be possibly confounded by sequencing run, and all analyses were performed independently for these two categories.
16S rRNA gene analysis
Request a detailed protocolWe assigned sequences to samples by exact matches of MID (multiplex identifier, 10 nt) sequences processed 16S rRNA sequences using the DADA2 pipeline, implemented in the DADA2 R package, version 1.16.0 (Callahan et al., 2016). Processing of the raw reads with DADA2 was performed separately for the DNA- and RNA-based libraries. Briefly, raw sequences were trimmed and quality filtered with the maximum two ‘expected errors’ allowed in a read, paired sequences were merged, ASVs were inferred, and chimeras removed. The two libraries were merged and another round of chimera removal was performed. We classified taxonomy using the Ribosomal Database Project (RDP) training set 16 (Cole et al., 2014). Classifications with low confidence at the genus level (<0.8) were grouped in the arbitrary taxon ‘unclassified_group’. For all downstream analyses, we rarefied samples to 10,000 reads each. Due to the quality filtering, we have phenotyping data for 286 individuals on DNA level, and 320 G2 individuals on RNA level.
We used the phyloseq R package (version 1.32.0) to estimate alpha diversity using the Shannon index and Chao1 index, and beta diversity using Bray-Curtis distance (McMurdie and Holmes, 2013). We defined core microbiomes at the DNA- and RNA-level, including taxa present in >25% of the samples and with median abundance of non-zero values > 0.2% for ASV and genus; and >0.5% for family, order, class, and phylum.
Genotyping
Request a detailed protocolWe extracted genomic DNA from ear samples using DNAeasy Blood and Tissue 96 well kits (Qiagen, Hilden, Germany), according to the manufacturer’s protocol. We sent DNA samples from 26 G0 mice and 320 G2 mice to GeneSeek (Neogen, Lincoln, NE) for genotyping using the Giga Mouse Universal Genotyping Array (GigaMUGA; Morgan et al., 2015), an Illumina Infinium II array containing 141,090 SNP probes. We quality-filtered genotype data using plink 1.9 (Chang et al., 2015); we removed individuals with call rates < 90% and SNPs that were: not bi-allelic, missing in > 10% individuals, with minor allele frequency < 5%, or Hardy-Weinberg equilibrium exact test <i>P-values < 1e-10. A total of 64,103 SNPs and all but one G2 individual were retained. Prior to mapping, we LD-filtered SNPs with r2 > 0.9 using a window of 5 SNPs and a step size of 1 SNP. We retain 32,625 SNPs for mapping.
Hybrid index calculation
Request a detailed protocolFor each G0 and G2 mouse, we estimated a hybrid index – defined as the percentage of M. m. musculus ancestry. We identified ancestry-informative SNP markers by comparing GigaMUGA data from ten individuals each from two wild-derived outbred stocks of M. m. musculus (Kazakhstan and Czech Republic) and two of M. m. domesticus (Germany and France) maintained at the Max Planck Institute for Evolutionary Biology (L.M. Turner and B. Payseur, unpublished data). We classified SNPs as ancestry informative if they had a minimum of 10 calls per subspecies, the major allele differed between musculus and domesticus, and the allele frequency difference between subspecies was >0.3. A total of 48,361 quality-filtered SNPs from the G0/G2 genotype data were informative, including 8775 SNPs with fixed differences between subspecies samples.
Estimation of relatedness among individuals
Request a detailed protocolWe computed a centred and a standardised relatedness matrix using the 32,625 filtered SNPs with GEMMA (v 0.98.1; Zhou and Stephens, 2012). The centred relatedness matrix was calculated with the formula:
The standardised relatedness matrix was calculated with the formula:
where X denotes the n×p matrix of genotypes, as its ith column representing the genotypes of ith SNP, as the sample mean and as a n×1 vector of 1’s, and as the sample variance of ith SNP.
Heritability of microbial abundances
Request a detailed protocolWe calculated heritabilities for bacterial abundances using linear mixed models implemented in the lme4qtl R package (version 0.2.2; Ziyatdinov et al., 2018). We included mating pair nested within the subcross identifier (Figure 1—figure supplement 6) as random effects to control for maternal effects and population structure, respectively. The narrow-sense heritability (h2) is expressed as:
where σg2 is the genetic variance estimated by the centred GRM, σm2 variance of the mating pair component, σs2 the variance due to the subcross identifier, and σe2 the variance due to residual environmental factors.
We also estimated CH using the method from Zhou et al., 2013, which estimates the variance explained by genotyped SNPs. CH is estimated using the standardised GRM.
We determined significance of the heritability estimates using exact restricted likelihood ratio tests, following Supplementary Note 3 in Ziyatdinov et al., 2018, using the exactRLRT() function of the R package RLRsim (version 3.1–6; Scheipl et al., 2008). Correlation with cospeciation rates was calculated for taxa shared between studies using the Spearman’s correlation test.
Genome-wide association mapping
Request a detailed protocolPrior to mapping, we inverse logistic transformed bacterial abundances using the inv.logit function from the R package gtools (version 3.9.2; Grieneisen et al., 2021).
We performed association mapping in the R package lme4qtl (version 0.2.2; Ziyatdinov et al., 2018) with the following linear mixed model:
where is the phenotypic value of the jth individual; μ is the mean, the additive and the dominance genotypic index values coded as for individual j at locus i. a and d indicate fixed additive and dominance effects, W indicates random effects mating pair and centred kinship matrix, plus residual error e. For mapping, we used the centred GRM as the SNP effect size does not depend on its minor allele frequency. Moreover, the centred GRM typically provides better control for population structure in lower organisms and both GRMs perform equally in humans (Zhou and Stephens, 2012).
We estimated additive and dominance effects separately because we expected to observe underdominance and overdominance in our hybrid mapping population, as well as additive effects, and aimed to estimate their relative importance. To model the additive effect (i.e. 1/2 distance between homozygous means), genotypes at each locus, i, were assigned additive index values ( ∈ 1, 0,–1) for AA, AB, BB, respectively, with A indicating the major allele and B the minor allele. To model dominance effects (i.e. heterozygote mean – midpoint of homozygote means), genotypes were assigned dominance index values ( ∈ 0, 1) for homozygotes and heterozygotes, respectively.
We included mating pair as a random effect to account for maternal effects and cage effects, as male litter mates are kept together in a cage after weaning. We included kinship coefficient as a random effect in the model to account for population and family structure. To avoid proximal contamination, we used a leave-one-chromosome-out approach, that is, when testing each single-SNP association we used a relatedness matrix omitting markers from the same chromosome (Parker et al., 2014). Hence, for testing SNPs on each chromosome, we calculated a centred relatedness matrix using SNPs from all other chromosomes with GEMMA (v0.98.1; Zhou and Stephens, 2012). We calculated p-values for single-SNP associations by comparing the full model to a null model excluding fixed effects. Code for performing the mapping is available at https://github.com/sdoms/mapping_scripts (copy archived at swh:1:rev:d085e7782e9ac85e264fc6b70a5058a53fd7e9fe, Doms, 2022).
The chi-squared test statistics needed for the calculation of the genomic inflation factors were computed from the p-values assuming one degree of freedom. The genomic inflation factor was defined as the median of the observed chi-squared test statistic divided by the expected median of the corresponding chi-squared distribution, and was computed for each trait separately. For traits with a genomic inflation factor (λGC) above the proposed threshold of 1.05, we applied genomic control by dividing the chi-squared test statistic with λGC (Devlin and Roeder, 1999).
We evaluated significance of SNP-trait associations using two thresholds; first, we used a genome-wide threshold for each trait, where we corrected for multiple testing across markers using the Bonferroni method (Abdi, 2007). Second, as bacteria interact with each other within the gut as members of a community, bacterial abundances are non-independent, so we calculated a study-wide threshold dividing the genome-wide threshold by the number of effective taxa included. We used matSpDlite (Nyholt, 2019; Li and Ji, 2005; Qin et al., 2020) to estimate the number of effective bacterial taxa based on eigenvalue variance.
To estimate the genomic interval represented by each significant LD-filtered SNP, we report significant regions defined by the most distant flanking SNPs in the full pre-LD-filtered genotype dataset showing r2 >0.9 with each significant SNP. We combined significant regions less than 10 Mb apart into a single region. Genes situated in significant regions were retrieved using biomaRt (Durinck et al., 2009) and the mm10 mouse genome.
Percentage of variance explained
Request a detailed protocolWe estimated the percentage of variance explained by the lead SNP using a linear mixed model with the additive and dominance genotypes of the lead SNP included as a fixed effect and mating pair and kinship matrix as random effects in lme4QTL (Ziyatdinov et al., 2018). We used the r.squaredGLMM function from the MuMIn R package (v 1.37.17; Kamil, 2020) to calculate the marginal R_GLMM², which represents the variance explained by the fixed effects (i.e. the genotype effect). The total variance explained by all significant SNPs was calculated similarly with the exception of including all significant SNPs as fixed effects instead of only the lead SNP. This was performed using the additive genotypes only and both the additive and the dominance genotypes.
Dominance analyses
Request a detailed protocolWe classified dominance for SNPs with significant associations on the basis of the d/a ratio (Falconer, 1996) where d is the dominance effect, a the additive effect. As the expected value under purely additive effects is 0. As our mapping population is a multiparental-line cross, and not all SNPs were ancestry-informative with respect to musculus/domesticus, the sign of a effects is defined by the major allele within our mapping population, which lacks clear biological interpretation. To provide more meaningful values, we report d/|a|, such that a value of 1=complete dominance of the allele associated with higher bacterial abundance, and a value of –1=complete dominance of the allele associated with lower bacterial abundance. Values above 1 or below –1 indicate over/underdominance. We classified effects of significant regions the following arbitrary d/|a| ranges to classify dominance of significant regions (Burke et al., 2002; Miller et al., 2014): underdominant < −1.25, high abundance allele recessive between –1.25 and –0.75, partially recessive between –0.75 and –0.25, additive between –0.25 and 0.25, partially dominant between 0.25 and 0.75, dominant 0.75 and 1.25, and overdominant > 1.25.
Gene ontology and network analysis
Request a detailed protocolThe nearest genes upstream and downstream of the significant SNPs were identified using the locateVariants() function from the VariantAnnotation R package (version 1.34.0; Obenchain et al., 2014) using the default parameters. A maximum of two genes per locus were included (one upstream, and one downstream of a given SNP).
To investigate functions and interactions of candidate genes, we calculated a PPI network with STRING version 11 (Szklarczyk et al., 2019), on the basis of a list of the closest genes to all SNPs with significant trait associations. We included network edges with an interaction score >0.9, based on evidence from fusion, neighbourhood, co-occurrence, experimental, text-mining, database, and coexpression. We exported this network to Cytoscape v 3.8.2 (Shannon et al., 2003) for identification of highly interconnected regions using the MCODE Cytoscape plugin (Bader and Hogue, 2003), and functional annotation of clusters using the stringApp Cytoscape plugin (Doncheva et al., 2019).
We identified over represented KEGG pathways and human diseases using the clusterprofiler R package (version 3.16.1; Yu et al., 2012). p-values were corrected for multiple testing using the Benjamini-Hochberg method. Pathways and diseases with an adjusted p-value < 0.05 were considered overrepresented.
Calculating overlap with other studies and overrepresentation of IBD genes
Request a detailed protocolTo test for significant overlap with loci identified in previous mapping studies and for overrepresentation of IBD genes, we used the tool poverlap (Pedersen and Brown, 2013) to compare observed overlap to random expectations based on 10,000 permutations of significant regions. Regio`ping regions using the locateVariants() function from the VariantAnnotation R package (version 1.34.0; Obenchain et al., 2014).
Candidate gene curation
Request a detailed protocolFrom the total set of 11,618 genes situated within significant regions, we identified a set of high-confidence genes, which met one or more of the following criteria: (1) hub genes or nearest neighbours in the ‘closest gene network’ (Figure 4, Figure 4—figure supplement 1), (2) genes in regions overlapping with QTL from previous mouse studies (Figure 5), (3) genes differentially expressed between germ-free and conventional mice (Figure 5—figure supplement 1 and Figure 5—figure supplement 2; Mills et al., 2020). After filtering, the resulting set of 309 genes were given as input into STRING (Szklarczyk et al., 2019) to construct a PPI network (304/309 genes were represented in the database; Figure 6—figure supplement 1). Finally, we selected one top candidate gene per significant region on the basis of network properties (degree and number of nodes) and/or fitting the most of the above-mentioned criteria. In case of a tie, the gene with the highest intestinal expression score (source: STRING) is chosen. Nodes without any edges were removed.
Data availability
DNA- and RNA-based 16S rRNA gene sequences are available under project accession number PRJNA759194. Code is available at https://github.com/sdoms/mapping_scripts, (copy archived at swh:1:rev:d085e7782e9ac85e264fc6b70a5058a53fd7e9fe).
-
SRAID PRJNA759194. DNA- and RNA-based 16S rRNA gene sequences.
References
-
The bonferonni and šidák corrections for multiple comparisonsEncyclopedia of Measurement and Statistics, SAGE 1:e9.
-
Understanding quantitative genetic variationNature Reviews. Genetics 3:11–21.https://doi.org/10.1038/nrg700
-
The effect of host genetics on the gut microbiomeNature Genetics 48:1407–1412.https://doi.org/10.1038/ng.3663
-
Speciation by symbiosisTrends in Ecology & Evolution 27:443–451.https://doi.org/10.1016/j.tree.2012.03.011
-
The roles of host evolutionary relationships (genus: Nasonia) and development in structuring microbial communitiesEvolution; International Journal of Organic Evolution 66:349–362.https://doi.org/10.1111/j.1558-5646.2011.01454.x
-
Genetic analysis of sunflower domesticationGenetics 161:1257–1267.https://doi.org/10.1093/genetics/161.3.1257
-
Dysbiosis of the gut microbiota in diseaseMicrobial Ecology in Health and Disease 26:26191.https://doi.org/10.3402/mehd.v26.26191
-
Feeding and the rhodopsin family g-protein coupled receptors in nematodes and arthropodsFrontiers in Endocrinology 3:157.https://doi.org/10.3389/fendo.2012.00157
-
They must hold tight: Junction proteinsMicrobiota And Immunity In Intestinal Mucosa’, Current Protein & Peptide Science Curr Protein Pept Sci 16:655–671.https://doi.org/10.2174/1389203716666150630133141
-
Innate immune recognition of the microbiota promotes host-microbial symbiosisNature Immunology 14:668–675.https://doi.org/10.1038/ni.2635
-
Gut microbiota as a missing link between nutrients and traits of humanFrontiers in Microbiology 9:1510.https://doi.org/10.3389/fmicb.2018.01510
-
Gut microbiota’s effect on mental health: The gut-brain axisClinics and Practice 7:987.https://doi.org/10.4081/cp.2017.987
-
Ribosomal database project: Data and tools for high throughput rrna analysisNucleic Acids Research 42:D633–D642.https://doi.org/10.1093/nar/gkt1244
-
Microbiota signaling pathways that influence neurologic diseaseNeurotherapeutics 15:135–145.https://doi.org/10.1007/s13311-017-0598-8
-
Genetic variation shapes murine gut microbiota via immunityTrends in Immunology 41:1–3.https://doi.org/10.1016/j.it.2019.11.009
-
Circadian disruption changes gut microbiome taxa and functional gene compositionFrontiers in Microbiology 9:737.https://doi.org/10.3389/fmicb.2018.00737
-
Targeting gut microbiota in obesity: effects of prebiotics and probioticsNature Reviews. Endocrinology 7:639–646.https://doi.org/10.1038/nrendo.2011.126
-
Genomic control for association studiesBiometrics 55:997–1004.https://doi.org/10.1111/j.0006-341x.1999.00997.x
-
Softwaremapping_scripts, version swh:1:rev:d085e7782e9ac85e264fc6b70a5058a53fd7e9feSoftware Heritage.
-
Cytoscape stringapp: Network analysis and visualization of proteomics dataJournal of Proteome Research 18:623–632.https://doi.org/10.1021/acs.jproteome.8b00702
-
Microbes and oxytocin: Benefits for host physiology and behaviorInternational Review of Neurobiology 131:91–126.https://doi.org/10.1016/bs.irn.2016.07.004
-
Stress & the gut-brain axis: Regulation by the microbiomeNeurobiology of Stress 7:124–136.https://doi.org/10.1016/j.ynstr.2017.03.001
-
The role of pattern recognition receptors in intestinal inflammationMucosal Immunology 6:451–463.https://doi.org/10.1038/mi.2013.13
-
Adaptive changes in the central control of energy homeostasis occur in response to variations in energy statusInternational Journal of Molecular Sciences 22:2728.https://doi.org/10.3390/ijms22052728
-
The treatment-naive microbiome in new-onset Crohn’s diseaseCell Host & Microbe 15:382–392.https://doi.org/10.1016/j.chom.2014.02.005
-
Genetic determinants of the gut microbiome in uk twinsCell Host & Microbe 19:731–743.https://doi.org/10.1016/j.chom.2016.04.017
-
Microbiome interactions shape host fitnessPNAS 115:E11951–E11960.https://doi.org/10.1073/pnas.1809349115
-
Unraveling the processes shaping mammalian gut microbiomes over evolutionary timeNature Communications 8:14319.https://doi.org/10.1038/ncomms14319
-
The “leaky gut”: Tight junctions but loose associations?Digestive Diseases and Sciences 65:1277–1287.https://doi.org/10.1007/s10620-019-05777-2
-
Breaking down the barriers: the gut microbiome, intestinal permeability and stress-related psychiatric disordersFrontiers in Cellular Neuroscience 9:e00392.https://doi.org/10.3389/fncel.2015.00392
-
Wild-caught rodents retain a majority of their natural gut microbiota upon entrance into captivityEnvironmental Microbiology Reports 6:191–195.https://doi.org/10.1111/1758-2229.12118
-
An introduction to phylosymbiosisProceedings. Biological Sciences 287:20192900.https://doi.org/10.1098/rspb.2019.2900
-
The human intestinal microbiome in health and diseaseThe New England Journal of Medicine 375:2369–2379.https://doi.org/10.1056/NEJMra1600266
-
Vitamin D and microbiota: Two sides of the same coin in the immunomodulatory aspectsInternational Immunopharmacology 79:106112.https://doi.org/10.1016/j.intimp.2019.106112
-
The hologenome concept: Helpful or hollow?PLOS Biology 13:e1002311.https://doi.org/10.1371/journal.pbio.1002311
-
The mouse universal genotyping array: From substrains to subspeciesG3: Genes, Genomes, Genetics 6:263–279.https://doi.org/10.1534/g3.115.022087
-
Functional redundancy-induced stability of gut microbiota subjected to disturbanceTrends in Microbiology 24:402–413.https://doi.org/10.1016/j.tim.2016.02.002
-
Role of trp channels in shaping gut microbiomeMEDICINE & PHARMACOLOGY 9:653.https://doi.org/10.20944/journals/202007.0653.v1
-
Gut commensal bacteria and regional Wnt gene expression in the proximal versus distal colonThe American Journal of Pathology 184:592–599.https://doi.org/10.1016/j.ajpath.2013.11.029
-
Host-gut microbiota metabolic interactionsScience 336:1262–1267.https://doi.org/10.1126/science.1223813
-
Genetic and environmental control of host-gut microbiota interactionsGenome Research 25:1558–1569.https://doi.org/10.1101/gr.194118.115
-
Using the natural variation of mouse populations to understand host-gut microbiome interactionsDrug Discovery Today. Disease Models 28:61–71.https://doi.org/10.1016/j.ddmod.2019.08.003
-
The Gut Feeling: GPCRs Enlighten the WayCell Host & Microbe 26:160–162.https://doi.org/10.1016/j.chom.2019.07.018
-
HAN-rotation, a new system for rigorous outbreedingZeitschrift Fur Versuchstierkunde 14:133–142.
-
Analysis of factors contributing to variation in the C57BL/6J fecal microbiota across German animal facilitiesInternational Journal of Medical Microbiology 306:343–355.https://doi.org/10.1016/j.ijmm.2016.03.004
-
Complement component C3 - The “Swiss Army Knife” of innate immunity and host defenseImmunological Reviews 274:33–58.https://doi.org/10.1111/imr.12500
-
Transgressive segregation, adaptation and speciationHeredity 83 (Pt 4):363–372.https://doi.org/10.1038/sj.hdy.6886170
-
Gut microbiota functions: metabolism of nutrients and other food componentsEuropean Journal of Nutrition 57:1–24.https://doi.org/10.1007/s00394-017-1445-8
-
The role of the microbiome in the neurobiology of social behaviourBiological Reviews of the Cambridge Philosophical Society 95:1131–1166.https://doi.org/10.1111/brv.12603
-
Size and power of tests for a zero random effect variance or polynomial regression in additive and linear mixed modelsComputational Statistics & Data Analysis 52:3283–3299.https://doi.org/10.1016/j.csda.2007.10.022
-
Wnt signalling at the crossroads of nutritional regulationThe Biochemical Journal 416:e11–e13.https://doi.org/10.1042/BJ20082074
-
Acetyl-CoA and the regulation of metabolism: mechanisms and consequencesCurrent Opinion in Cell Biology 33:125–131.https://doi.org/10.1016/j.ceb.2015.02.003
-
Using the Mus musculus hybrid zone to assess covariation and genetic architecture of limb bone lengthsMolecular Ecology Resources 18:908–921.https://doi.org/10.1111/1755-0998.12776
-
Unravelling the effects of the environment and host genotype on the gut microbiomeNature Reviews. Microbiology 9:279–290.https://doi.org/10.1038/nrmicro2540
-
Links between natural variation in the microbiome and host fitness in wild mammalsIntegrative and Comparative Biology 57:756–769.https://doi.org/10.1093/icb/icx104
-
Host genetic determinants of the gut microbiota of wild miceMolecular Ecology 28:3197–3207.https://doi.org/10.1111/mec.15139
-
The gut microbiota and Bergmann’s rule in wild house miceMolecular Ecology 29:2300–2311.https://doi.org/10.1111/mec.15476
-
Multiple muscarinic pathways mediate the suppression of voltage-gated Ca2+ channels in mouse intestinal smooth muscle cellsBritish Journal of Pharmacology 158:1874–1883.https://doi.org/10.1111/j.1476-5381.2009.00475.x
-
Reclassification of eubacterium formicigenerans holdeman and moore 1974 as dorea formicigenerans genInternational Journal of Systematic and Evolutionary Microbiology 52:423–428.https://doi.org/10.1099/00207713-52-2-423
-
Deciphering functional redundancy in the human microbiomeNature Communications 11:6217.https://doi.org/10.1038/s41467-020-19940-1
-
Reduced male fertility is common but highly variable in form and severity in a natural house mouse hybrid zoneEvolution; International Journal of Organic Evolution 66:443–458.https://doi.org/10.1111/j.1558-5646.2011.01445.x
-
Light exposure influences the diurnal oscillation of gut microbiota in miceBiochemical and Biophysical Research Communications 501:16–23.https://doi.org/10.1016/j.bbrc.2018.04.095
-
Mucosal-associated microbiota other than luminal microbiota has a close relationship with diarrhea-predominant irritable bowel syndromeFrontiers in Cellular and Infection Microbiology 10:515614.https://doi.org/10.3389/fcimb.2020.515614
-
Genome-wide efficient mixed-model analysis for association studiesNature Genetics 44:821–824.https://doi.org/10.1038/ng.2310
Article and author information
Author details
Funding
Deutsche Forschungsgemeinschaft (261376515)
- John F Baines
Deutsche Forschungsgemeinschaft (TU 500/2-1)
- Leslie M Turner
Deutsche Forschungsgemeinschaft (390884018)
- John F Baines
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Diethard Tautz for generous support of mouse breeding and Camilo Medina and the MPI-Plön mouse team for performing mouse husbandry, and Katja Cloppenborg-Schmidt and Dr Sven Künzel for their excellent technical assistance. We thank Jason Wolf for constructive feedback on the statistical models and heritability estimation. We thank Mathieu Groussin for assistance with cospeciation rate data. Research funding for this project was provided by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Project-ID 261376515 – Collaborative Research Center 1182, ‘Origin and Function of Metaorganisms’, Project A2 (J.F.B. and A.F.), Cluster of Excellence 2167 ‘Precision Medicine in Chronic Inflammation (PMI)’ (grant no. EXC2167 to A.F. and J.F.B.) and TU 500/2–1 to L.M.T, and by the Max Planck Society (to D Tautz).
Ethics
This study was performed according to approved animal protocols and institutional guidelines of the Max Planck Institute. Mice were maintained and handled in accordance with FELASA guidelines and German animal welfare law (Tierschutzgesetz § 11, permit from Veterinä;ramt Kreis Plö;n: 1401-144/PLÖ;-004697).
Copyright
© 2022, Doms et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,857
- views
-
- 361
- downloads
-
- 13
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Evolutionary Biology
The majority of highly polymorphic genes are related to immune functions and with over 100 alleles within a population, genes of the major histocompatibility complex (MHC) are the most polymorphic loci in vertebrates. How such extraordinary polymorphism arose and is maintained is controversial. One possibility is heterozygote advantage (HA), which can in principle maintain any number of alleles, but biologically explicit models based on this mechanism have so far failed to reliably predict the coexistence of significantly more than ten alleles. We here present an eco-evolutionary model showing that evolution can result in the emergence and maintenance of more than 100 alleles under HA if the following two assumptions are fulfilled: first, pathogens are lethal in the absence of an appropriate immune defence; second, the effect of pathogens depends on host condition, with hosts in poorer condition being affected more strongly. Thus, our results show that HA can be a more potent force in explaining the extraordinary polymorphism found at MHC loci than currently recognized.
-
- Evolutionary Biology
Spatial patterns in genetic diversity are shaped by individuals dispersing from their parents and larger-scale population movements. It has long been appreciated that these patterns of movement shape the underlying genealogies along the genome leading to geographic patterns of isolation by distance in contemporary population genetic data. However, extracting the enormous amount of information contained in genealogies along recombining sequences has, until recently, not been computationally feasible. Here we capitalize on important recent advances in genome-wide gene-genealogy reconstruction and develop methods to use thousands of trees to estimate per-generation dispersal rates and to locate the genetic ancestors of a sample back through time. We take a likelihood approach in continuous space using a simple approximate model (branching Brownian motion) as our prior distribution of spatial genealogies. After testing our method with simulations we apply it to Arabidopsis thaliana. We estimate a dispersal rate of roughly 60km2 per generation, slightly higher across latitude than across longitude, potentially reflecting a northward post-glacial expansion. Locating ancestors allows us to visualize major geographic movements, alternative geographic histories, and admixture. Our method highlights the huge amount of information about past dispersal events and population movements contained in genome-wide genealogies.