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

  1. Shauni Doms
  2. Hanna Fokt
  3. Malte Christoph Rühlemann
  4. Cecilia J Chung
  5. Axel Kuenstner
  6. Saleh M Ibrahim
  7. Andre Franke
  8. Leslie M Turner  Is a corresponding author
  9. John F Baines  Is a corresponding author
  1. Max Planck Institute for Evolutionary Biology, Germany
  2. Section of Evolutionary Medicine, Institute for Experimental Medicine, Kiel University, Germany
  3. Institute for Clinical Molecular Biology (IKMB), Kiel University, Germany
  4. Institute for Medical Microbiology and Hospital Epidemiology, Hannover Medical School, Germany
  5. Institute of Experimental Dermatology, University of Lübeck, Germany
  6. Sharjah Institute of Medical Research, United Arab Emirates
  7. Milner Centre for Evolution, Department of Biology & Biochemistry, University of Bath, United Kingdom

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.sa0

eLife 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).

Figure 1 with 6 supplements see all
Heritability estimates and their relationship to cospeciation rates.

(A–B) Lollipop chart showing the values of the chip heritability (dark red circles), narrow-sense heritability (green squares), PVE by the additive effect (blue triangles), and additive and dominance effect (yellow diamonds) of all significant single nucleotide polymorphisms (SNPs) within one taxon for DNA-based traits (A) and RNA-based traits (B) Only taxa with a non-zero narrow-sense heritability estimate are shown. Only the outline of non-significant heritability estimates (p<0.05, Restricted Likelihood Ratio test) are shown. Taxa without significant hits have no PVE reported. (C–D) Relationship between the narrow-sense heritability estimates for the relative abundance of bacterial taxa and cospeciation rate for the same genus calculated by Groussin et al., 2017. DNA level (C), and RNA level (D). The blue line represents a linear regression fit to the data and the grey area the corresponding confidence interval. p-Values and the correlation coefficient are calculated using the Spearman’s correlation test. The text labels on the y-axis (A–B) and the text labels in panels (C-D) are coloured according to taxonomic level: amplicon sequence variants (ASV) in black, genus in purple, family in light blue, order in red, class in green, and phylum in yellow.

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.

Figure 2 with 2 supplements see all
Heatmap of significant host loci from association mapping of bacterial abundances.

Karotype plot showing the number of significant loci found using a study-wide threshold, where (A) plots the significance intervals and (B) the significant single nucleotide polymorphisms (SNP) markers on the chromosomes. The position of the SNPs on panel (B) has been amplified by 0.5 Mb to visualise it. The position of the genes closest to SNPs with the lowest p-values (Tlr4, and Irak4 and Adamsts20) are indicated.

Table 1
Overview of mapping statistics.

Loci with a P-value below the study-wide threshold (P<1.29 × 10–8) are considered significant. ‘Significant loci total P’ are the loci having a significant P-value from the total model (additive and dominance effect), ‘Significant loci additive P’ have a significant additive effect, and ‘Significant loci dominance P’ have a significant dominance effect.

DNARNATotal
Mapped taxa101142153
Taxa with significant loci6593120
Median interval size (Mb)1.322.51.91
Total significant SNPs316596782
Total significant loci4437461184
Unique significant loci172305428
Significant loci total P83144204
Significant loci additive P144244351
Significant loci dominance P88171230
Median significant loci per trait558
Median unique significant loci per trait334
Median unique significant SNPs per locus222
Median number of genes per locus325443
Median protein coding genes per locus111714
  1. SNPs: single nucleotide polymorphisms.

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).

Genetic architecture of significant loci.

(A) Source of the allele with the highest phenotypic value. (B) Histogram of dominance values d/a of significant loci reveals a majority of loci acting recessive or partially recessive. (C) Histogram showing the percentage of variance explained (PVE) by the peak single nucleotide polymorphisms (SNP) for DNA (left) and RNA (right).

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).

Figure 4 with 4 supplements see all
High confidence protein-protein interaction (PPI) network of genes closest to single nucleotide polymorphisms (SNPs) significantly associated with bacterial abundances.

Network clusters are annotated using STRING’s functional enrichment (Doncheva et al., 2019). Nodes represent proteins and edges their respective interactions. Only edges with an interaction score higher than 0.9 are retained. The width of the edge line expresses the interaction score calculated by STRING. The colour of the nodes describes the expression of the protein in the intestine where yellow is not expressed and purple is highly expressed.

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.

Figure 5 with 2 supplements see all
Significant loci in this study previously found in quantitative trait loci (QTL) studies of the mouse gut microbiome.

The genes present in two repeatedly identified regions are depicted in boxes.

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.

Figure 6 with 1 supplement see all
Network of host candidate genes influencing bacterial traits using STRING (https://string-db.org).

The nodes represent proteins and are coloured according to a selection of enriched GO terms and pathways: G-protein coupled receptor (GPCR) signalling (red), regulation of the immune system process (blue), response to nutrient levels (light green), fatty acid metabolic process (pink), glucose homeostasis (purple), response to antibiotic (orange), regulation of feeding behaviour (yellow), positive regulation of insulin secretion (dark green), circadian entrainment (brown), and response to vitamin D (turquoise). The colour of the edges represents the interaction type: known interactions from curated databases (turquoise) or experimentally determined (pink); predicted interactions from gene neighbourhood (green), gene fusions (red), gene co-occurrence (blue); other interactions from text-mining (light green), coexpression (black), and protein homology (purple).

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

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
strain, strain background (Mus musculus musculus-Mus musculus-domesticus, males)(HZ)A-(HZ)HThis paperSecond generation wild-derived inter-crossed hybrid mouse line, Mus musculus musculus-Mus musculus domesticus originating from four breeding stocks captured in the wild hybrid zone around Freising, Germany.
biological sample (Mus musculus)Cecum tissue (mucosa)This paperCollected and stored in RNAlater overnight. After RNAlater removal tissue was stored in –20°C
biological sample (Mus musculus)Ear clipsThis paperFor genotyping
commercial assay or kitAllprep DNA/RNA 96-wellQiagen (Hilden, Germany)Cat. No.: 80,284DNA/RNA extraction
commercial assay or kitLysing matrix EMP Biomedical (Eschwege, Germany)SKU:116914050-CFLysing
commercial assay or kitHigh-capacity cDNA Reverse Transcription KitApplied Biosystems (Darmstadt, Germany)Cat. No.: 368,814cDNA transcription
sequence-based reagent27 Fhttps://doi.org/10.1016/j.ijmm.2016.03.004PCR primersForward primer V1-V2 hypervariable region
sequence-based reagent338 Rhttps://doi.org/10.1016/j.ijmm.2016.03.004PCR primersReverse primer V1-V2 hypervariable region
software, algorithmDADA2https://doi.org/10.1038/nmeth.386916S rRNA gene processing
software, algorithmphyloseqhttps://doi.org/10.1371/journal.pone.006121716S rRNA gene analysis
commercial assay or kitDNAeasy Blood and Tissue 96-wellQiagen (Hilden, Germany)Cat. No.: 69,504DNA extraction ear clips for genotyping
OtherGigaMUGANeogen, Lincoln, NEIllumina Infinium II array containing 141,090 SNP probes
software, algorithmplink 1.9https://doi.org/10.1186/s13742-015-0047-8Quality control genotypes
software, algorithmGEMMA (v 0.98.1)https://doi.org/10.1038/ng.2310Genetic relatedness matrix
software, algorithmlme4QTLhttps://doi.org/10.1186/s12859-018-2057-xSNP-based heritability, GWAS
software, algorithmexactLRT (RLRsim, v 3.1–6)Scheipl et al., 2008Significance heritability estimates
software, algorithminv.logit (Gtools, v 3.9.2)Grieneisen et al., 2021Inverse logistic transformation
software, algorithmmatSpDliteNyholt, 2019; https://doi.org/10.1038/sj.hdy.6800717Study-wide significance threshold
software, algorithmbiomaRt (mm10)Durinck et al., 2009Gene annotation
software, algorithmr.squaredGLMM (MuMIn, v 1.37.17)Kamil, 2020Percentage of variance explained
software, algorithmlocateVariants (VariansAnnotation, v 1.34.0)https://doi.org/10.1093/bioinformatics/btu168Nearest gene
software, algorithmSTRING (v 11)https://doi.org/10.1093/nar/gky1131Protein-protein interaction networks
software, algorithmCytoscape (v 3.8.2)https://doi.org/10.1101/gr.1239303Network analysis and visualisation
software, algorithmMCODEhttps://doi.org/10.1186/1471-2105-4-2Cytoscape plugin for identifying network clusters
software, algorithmstringApphttps://doi.org/10.1021/acs.jproteome.8b00702Cytoscape plugin for functional annotation
software, algorithmclusterprofiler (v 3.16.1)https://doi.org/10.1089/omi.2011.0118Enrichment analysis
software, algorithmpoverlapPedersen and Brown, 2013To determine significant overlap between studies
software, algorithmR (v 3.5.3)https://www.R-project.org/

Intercross design

Request a detailed protocol

We 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 in­stitutional 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 protocol

Mice 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 protocol

We 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 protocol

We 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 protocol

We 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 protocol

For 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 protocol

We 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:

Centered GRM=1p i=1p(xi1nx¯i)(xi 1nx¯i)T

The standardised relatedness matrix was calculated with the formula:

Standardized GRM=1p i=1p1vxi(xi1nx¯i)(xi 1nx¯i)T

where X denotes the n×p matrix of genotypes, xi as its ith column representing the genotypes of ith SNP, x¯i as the sample mean and 1n as a n×1 vector of 1’s, and vxi as the sample variance of ith SNP.

Heritability of microbial abundances

Request a detailed protocol

We 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:

h2=σg2σg2+σm2+σs2+σe2

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 protocol

Prior 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:

yi=μ+aiXija+diXijd+Wu+e

where yi is the phenotypic value of the jth individual; μ is the mean, Xija the additive and Xijd 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 (Xa ∈ 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 (Xd ∈ 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 protocol

We 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 protocol

We 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 protocol

The 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 protocol

To 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 protocol

From 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).

The following data sets were generated

References

    1. Abdi H
    (2007)
    The bonferonni and šidák corrections for multiple comparisons
    Encyclopedia of Measurement and Statistics, SAGE 1:e9.
    1. Castoldi A
    (2015) They must hold tight: Junction proteins
    Microbiota And Immunity In Intestinal Mucosa’, Current Protein & Peptide Science Curr Protein Pept Sci 16:655–671.
    https://doi.org/10.2174/1389203716666150630133141
  1. Book
    1. Falconer D. S
    (1996)
    Introduction to quantitative genetics
    Harlow, England: Prentice Hall).
  2. Software
    1. Nyholt DR
    (2019) matSpD, version matSpD
    Statistical and Genomic Epidemiology Laboratory.
    1. Rapp KG
    (1972)
    HAN-rotation, a new system for rigorous outbreeding
    Zeitschrift Fur Versuchstierkunde 14:133–142.
  3. Book
    1. Walsh B
    2. Lynch M
    (1998)
    Genetics and Analysis of Quantitative Traits
    Sunderland, MA: Sinauer.

Decision letter

  1. Jenny Tung
    Reviewing Editor; Duke University, United States
  2. Wendy S Garrett
    Senior Editor; Harvard T.H. Chan School of Public Health, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "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" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Wendy Garrett as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1. Demonstrate the robustness of the heritability analyses, including consideration of potential technical and biological confounding variables and ideally using complementary analyses (e.g., distinct methods for estimating heritability).

2. Demonstrate the robustness of the phylosymbiosis (heritability-cospeciation rate correlation) to controlling for taxon abundance.

3. Control for population structure are well as relatedness in the trait mapping analysis.

Reviewer #1 (Recommendations for the authors):

1. Heritability estimates in the study are quite high, including traits that are more highly heritable than even canonical high heritability traits in humans, like height. This could be accurate, particularly in a lab-controlled environment where genetic variation is ramped up by the crossing design and hybrid origins of the founder, but I think requires closer examination. In particular, I was surprised that there does not appear to be any adjustment for covariates in the heritability analysis: in my experience, batch and technical effects are omnipresent in genomic studies, including microbiome analysis, and can affect heritability estimates in unpredictable ways. Although sex and age are controlled for in the study design, are there differences in housing, diet, sample processing, library prep, or sequencing quality/read yield/lane/batch that account for measurable variance and need to be taken into account?

2. Related, part of the benefit of estimating heritability in a mixed effects models framework is the ability to flexibly control for other covariates. I don't see a benefit in using Mantel tests to correlate genetic relatedness and microbiome structure overall-as heritability is much more directly supported by direct estimation-and suggest removing this part of the analysis. I realize that the genetic relatedness-microbiome structure analysis might provide a "holistic" picture, but that could also be done by analyzing the top PCs of the microbiome data in a formal h2 analysis.

3. For the individual trait mapping analysis, I have some concerns that the models might not adequately control for background population structure. It's a good step to control for genetic structure with the K matrix, but I would like to see that the results are robust to including the top PCs of the genotype data in the model as additional fixed effects. Where possible, you could also consider using LD score regression to check that there is no unexpected inflation of test statistics.

4. In Figure 3 and lines 251 – 258, the results report PVE levels that are impossible (since more than 100% of variance can't be explained). While this is noted in the text, and some potential explanations are provided (lines 256 – 258), the current presentation isn't very useful because it's unclear how much trait variance is really being explained. A better alternative might be to construct the equivalent of polygenic scores for each trait and ask how much variance in each associated phenotype these composite scores explain (this approach would provide one value per mapped trait, upper-bounded by the overall trait heritability).

5. Finally, I lost the thread a bit when reading the second half of the results. All of these sections focus on one type of enrichment/annotation analysis or another. However, the hypotheses or models being tested in several of the analyses were not clear to me. E.g., why should we expect microbiome-associated genes to encode proteins that are networked to one another? Is the relationship to the GF versus conventional comparison based on the idea that if genes affect the microbiome, then changes in the microbiome should also affect the proteins those genes encode (not obvious to me)? Ultimately, what does it mean to be a "promising" candidate gene (promising for what kind of application), and how is that defined by these enrichment analyses?

Rather than having readers try to fill in the rationale for these analyses, it would be helpful to clarify the motivation in the paper itself-or even shorten these sections to focus on those where there is the clearest underlying model or hypothesis. This might also avoid diluting the more interesting, less speculative findings presented earlier.

Reviewer #2 (Recommendations for the authors):

Overall, this was an excellent study and I have few comments for improvement. A couple of points:

The association between heritability of taxa and their co-speciation rates is interesting, but I wonder whether this could be explained by an association between abundance and statistical power to detect heritability/co-speciation. Have the authors considered modeling heritability as a function of co-speciation rate and relative abundance? Such an analysis would be important to determine whether heritability and co-speciation rate are associated independently of abundance.

Relative abundances of microbes were used for all analyses (as opposed to absolute abundance quantification). Because findings based on relative abundances can be difficult to interpret, more discussion on this limitation would be helpful.

My understanding is that heritability estimates were generated only for ASV relative abundances. However, previous heritability studies have examined taxonomic levels from ASV to phylum. It would be interesting to extend the analyses presented to higher taxonomic levels beyond ASVs (eg species to phylum).

Reviewer #3 (Recommendations for the authors):

The main point of concern for me are the heritability estimates. The paper reports heritability results that may be a bit surprising considering the current knowledge and literature: the heritability values are very high, with several values around 90% or higher. This is unexpected, and I am not sure how this is reconciled with the expectation that most variation in the microbiome is environmental rather than genetic. Although some potential reasons are given in the Discussion (mice raised in a controlled environment, using cecal content, etc), it still makes me a bit uneasy to see such high heritability estimates. One potential way to approach this is to try a different statistical approach for calculating heritability. Another would be to compare heritability estimates from this study with estimates from other studies – by now there are quite a few studies that report microbiome-wide heritability estimates from humans, mice, and other host species. It could be useful to correlate the heritability estimates in the current study with those from these studies. Lastly, it would be good to compare microbiome heritability estimates from this study to heritability’s of other complex traits in the same system – are there other known phenotypes that have heritability estimates that are this high?

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "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" for further consideration by eLife. Your revised article has been evaluated by Wendy Garrett (Senior Editor) and a Reviewing Editor.

All reviewers appreciate the thorough revisions and re-analyses, which have resulted in a much improved manuscript. However, there are a few remaining issues that need to be addressed, as outlined below (the first item is the most essential):

1) Address the concern of Reviewer 3 below, which concern inflated heritability values in the main text. While the correlations in h2 estimates between methods are reassuringly high, Figure 8 in the response to reviewers also shows that the h2 reported in the paper (from lme4QTL) are typically several-fold (and sometimes an order of magnitude) higher than obtained from three alternative methods. This information would not be available to readers in the manuscript's current form. The response to reviewers indicates that lme4qtl was chosen because it had a high R2 with the PVE explained by the additive effects of significant SNPs (Figure 4 in the response to reviewers). However, both GEMMA and sommer produce similar R2 values (and in one case, a higher R2 value) but much lower heritability estimates. Do these lower estimates also correlate with co-speciation rates?

2) Consider the polishing revisions to the code repository suggested by Reviewer 3.

3) Consider integrating the useful explanation in the response to reviewers about experimental design controls for technical and batch effects in the main manuscript, as future readers may have similar questions.

Reviewer #2 (Recommendations for the authors):

The authors addressed all of my comments from the previous round of review. The inclusion of multiple alternative calculations of heritability have improved the manuscript substantially.

All 16S rRNA gene sequence data have been deposited to NCBI.

Reviewer #3 (Recommendations for the authors):

I thank the authors for a comprehensive revision that addressed many of my concerns. It was good to see that the heritability values are correlated across different methods, and that several other methods produced heritability values smaller than those generated by lme4QTL. This, together with the fact that heritability values for some bacteria are higher than those for traits like length and weight, supports my notion that lme4QTL heritability values are overestimated. I am not sure what the reason is -- it's hard to know without spending time digging into the data, and I might not have the statistics background to advise on this -- but I am not confident that these values are robust. I would encourage the authors to investigate these analyses very thoroughly, making sure that all potential confounders are accounted for, the models are reasonable, and there are no artifacts in the data. I would suggest the text includes the results of heritability analysis with other approaches (in addition to lme4QTL) more prominently: Figure 1 should report and visualize heritability values from all methods used, and visualize the correlations between them. The text should describe these results, the methods used, the heritability values reported and their correlation. There should be a clear discussion about the possible reasons for the high heritability values (and why they are lower using other methods). I would also suggest including the analysis comparing heritability values across studies in the text, and include a visualization of these correlations, rather than just reporting the p-value.

Regarding code availability, I want to thank the authors for enhancing the README file on the github repository, which now provides a nice description of the pipeline and analysis steps. However, I am not sure if this is sufficient for readers who want to reproduce the results: looking at the code itself, it seems like there are commands to load scripts that are not included in the repository (e.g. the snp_heritability_lme4qtl.R script loads the script function_for_gemma.r that I couldn't find anywhere), and these scripts might not be able to be run on other machines. I recommend amending the github repository and scripts so that anyone who wishes to do so is able to run the analysis and reproduce the results.

https://doi.org/10.7554/eLife.75419.sa1

Author response

Essential revisions:

1. Demonstrate the robustness of the heritability analyses, including consideration of potential technical and biological confounding variables and ideally using complementary analyses (e.g., distinct methods for estimating heritability).

2. Demonstrate the robustness of the phylosymbiosis (heritability-cospeciation rate correlation) to controlling for taxon abundance.

3. Control for population structure are well as relatedness in the trait mapping analysis.

Reviewer #1 (Recommendations for the authors):

1. Heritability estimates in the study are quite high, including traits that are more highly heritable than even canonical high heritability traits in humans, like height. This could be accurate, particularly in a lab-controlled environment where genetic variation is ramped up by the crossing design and hybrid origins of the founder, but I think requires closer examination. In particular, I was surprised that there does not appear to be any adjustment for covariates in the heritability analysis: in my experience, batch and technical effects are omnipresent in genomic studies, including microbiome analysis, and can affect heritability estimates in unpredictable ways. Although sex and age are controlled for in the study design, are there differences in housing, diet, sample processing, library prep, or sequencing quality/read yield/lane/batch that account for measurable variance and need to be taken into account?

We thank the reviewer for raising this important concern and agree that the heritability estimates appear high. We accordingly performed several analyses, which do confirm the robustness of the results. Details are given below in response to comment 3 and the major issue of Reviewer 3. In short, this resulted in the inclusion of a “sub-cross identifier” as an additional random effect (see Figure 1—figure supplement 4 for cross design). The current model includes the sub-cross identifier with the mating pair identifier nested within it (Author response image 1). We now clarify this in the text (see lines 834-836). The mating pair id accounts for both maternal effects and cage effects, as the mice were kept together with littermates after weaning and separated into single cages one week prior to dissection. All mice were kept in the same room.

Author response image 1
Visualization of the random effects: mating pair identifier (A) and mating pair identifier nested within the sub-cross ID (B).

We also agree that batch and technical effects are important to account for. Due to the setup of our extraction and sequencing library prep procedures, these factors were inherently accounted for at numerous levels. The samples were extracted using the Qiagen DNA/RNA Allprep kit in 96-well plate format, such that all samples were extracted in the same extraction round and hence timepoint. Individual sequencing libraries were subsequently 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.

2. Related, part of the benefit of estimating heritability in a mixed effects models framework is the ability to flexibly control for other covariates. I don't see a benefit in using Mantel tests to correlate genetic relatedness and microbiome structure overall-as heritability is much more directly supported by direct estimation-and suggest removing this part of the analysis. I realize that the genetic relatedness-microbiome structure analysis might provide a "holistic" picture, but that could also be done by analyzing the top PCs of the microbiome data in a formal h2 analysis.

Thank you for this suggestion, we agree that this result is not among the most important of our study and followed this advice as an opportunity to make the Results section more concise, also in reference to comment #5 of the same reviewer below.

3. For the individual trait mapping analysis, I have some concerns that the models might not adequately control for background population structure. It's a good step to control for genetic structure with the K matrix, but I would like to see that the results are robust to including the top PCs of the genotype data in the model as additional fixed effects. Where possible, you could also consider using LD score regression to check that there is no unexpected inflation of test statistics.

Thank you, we agree these suggestions would help show the robustness of the mapping analysis. Accordingly, we compared three different models and ran them on a representative subset of microbial traits (n=27 RNA-based genus level traits): (1) Current model with the mating pair and the kinship matrix as random effects, (2) Current model + PC1-PC3, (3) Current model + the mating pair id nested within the sub-cross identifier (see Figure 1—figure supplement 4 for cross design). We first compared the BIC scores of these models without a SNP effect (Author response table 1). We see that the BIC is, except for one trait (G_Odoribacter), the lowest (bold) in our current model.

Author response table 1
BIC scores models for RNA-based genera.
TaxaCurrent+PC1+PC2+PC3+PC4+PC5Current +subcross
G_Acetatifactor-2651.9-2637.7-2628.1-2613.8-2599.6-2584.5-2646.2
G_Alistipes-2714-2711.2-2701.5-2687.3-2673.9-2658.9-2712.8
G_Anaerostipes-3416.6-3400.9-3385.2-3369.2-3352.7-3336.1-3410.8
G_Bacteroides-2076.6-2065.9-2053.9-2045.2-2033.5-2021.5-2070.9
G_Butyricicoccus-3676.1-3661.1-3645.3-3628.4-3612.7-3596.6-3671.4
G_Clostridium_XlVa-2356.8-2354.2-2340.6-2327.2-2314.8-2301-2351.1
G_Dorea-2979.7-2965.4-2952-2936.6-2926.1-2916.2-2974.4
G_Eisenbergiella-1500.5-1492-1489.3-1479.8-1469.3-1458.2-1495.3
G_Fusicatenibacter-2932.8-2918.3-2902.8-2888.3-2873.5-2857.8-2927.1
G_Helicobacter-1145.4-1138.4-1131-1122.1-1114.8-1104.8-1140.3
G_Hungatella-2415.2-2402.3-2390.2-2376.7-2364.5-2351.1-2409.5
G_Intestinimonas-3286.5-3272.5-3256.3-3240.1-3228.2-3213.9-3282.4
G_Lactobacillus-2407-2395.6-2384.4-2380.5-2366.8-2352.6-2403.5
G_Marvinbryantia-2961.6-2950.7-2934.7-2918.8-2903.2-2887.3-2956.4
G_Mucispirillum-1938.8-1930.8-1918.9-1907.3-1896.7-1884.8-1933.1
G_Odoribacter-3183.2-3170.8-3156-3143.5-3134.7-3120-3186.1
G_Oscillibacter-2117.2-2105.8-2095.2-2084.3-2072.1-2060-2111.9
G_Paraprevotella-2230.2-2222.8-2219.2-2207.7-2195.4-2182.8-2228.4
G_Roseburia-2719-2704-2689.4-2674.8-2661.3-2647.6-2713.5
G_Staphylococcus-7666.6-7636.9-7607.6-7576.8-7546.5-7515.8-7660.9
unclassified_C_Deltaproteobacteria-3208.6-3197.8-3183.4-3167.3-3158.7-3142.1-3205.5
unclassified_F_Lachnospiraceae-1435.4-1433.5-1423-1413.9-1407.2-1396.8-1433.1
unclassified_F_Porphyromonadaceae-2997.7-2984.6-2970.3-2956.3-2943.4-2928.9-2991.9
unclassified_F_Prevotellaceae-4216.4-4198.5-4182-4175.1-4156.4-4136.9-4212.7
unclassified_F_Ruminococcaceae-2590.6-2578.8-2565.6-2554.4-2540.7-2525.9-2586.1
unclassified_O_Bacteroidales-2348-2339.5-2328-2317-2305.3-2291.7-2344.4
unclassified_O_Clostridiales-2721.2-2708-2694.8-2681-2667.4-2654.7-2715.4
unclassified_P_Bacteroidetes-3260.7-3252.8-3237.9-3222.6-3207-3194.2-3256.9

Next, we calculated the genomic inflation factors (Author response image 2) as well as the LD score regression intercept (Author response image 3) for this subset of traits. Running an LD score regression on data sets with less than 200k SNPs is however not advised (Bulik-Sullivan et al., 2015). Both the genomic inflation factor and the LD score regression intercept show the lowest inflation in the current model. Some traits did however show some inflation of the test statistics above the proposed threshold of 1.05. For this reason, we used genomic control for the traits with a genomic inflation factor above 1.05 (lines 887-895). Following this correction, the total number of study-wide significant loci changed from 443 to 428; all downstream analyses were accordingly repeated with the final updated list of loci.

Author response image 2
Genomic inflation factor for different P values (additive effect, left; dominance effect, middle; total model, right) calculated by the different models: current model including the mating pair id nested in the sub-cross ID (top), current model (middle), and the current model including PC1-PC3 as fixed effects (bottom).

Dashed line and values indicate the average genomic inflation factor.

Author response image 3
LD score regression intercept for P values calculated by the different models: current model including the mating pair id nested in the sub-cross ID (top), current model (middle), and the current model including PC1-PC3 as fixed effects (bottom).

Dashed line and values indicate the average genomic inflation factor.

4. In Figure 3 and lines 251 – 258, the results report PVE levels that are impossible (since more than 100% of variance can’t be explained). While this is noted in the text, and some potential explanations are provided (lines 256 – 258), the current presentation isn’t very useful because it’s unclear how much trait variance is really being explained. A better alternative might be to construct the equivalent of polygenic scores for each trait and ask how much variance in each associated phenotype these composite scores explain (this approach would provide one value per mapped trait, upper-bounded by the overall trait heritability).

We thank the reviewer for this suggestion. We used a linear mixed model that includes all the significant SNPs per trait to calculate a polygenic score equivalent. We adapted the manuscript Figure 3 accordingly and added a section to the Methods (line 908-917):

"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."

We used these "polygenic scores" to determine the choice of linear models and software/methods for calculating the heritability estimates. As seen in Author response image 4, using lme4QTL including the mating pair ID nested in the sub-cross ID (lme4qtl.fam.heritability) to estimate the heritability has a slope of 0.921, which is the closest to 1 of all the methods tested, and yields the highest R2 value (0.14).

Author response image 4
Correlation of heritability estimates calculated using different methods/programs and models to the percentage of phenotypic variance explained by the additive genotypes of all significant SNPs within the trait.

For Gaston, Sommer, Gemma 1PC, and lme4QTL.1PC, a linear mixed model was used with the first PC of the genotypes as a fixed effect and a genomic relatedness matrix (GRM) as a random effect. For Gemma.3PC, the first three PCs of the genotype data were used as fixed effects and a GRM as a random effect. For lme4QTL.MP, the mating pair ID was used as a random effect together with the GRM, while in lme4QTL.fam the mating pair ID was nested in the sub-cross ID as random effects together with a GRM to estimate the additive heritability.

5. Finally, I lost the thread a bit when reading the second half of the results. All of these sections focus on one type of enrichment/annotation analysis or another. However, the hypotheses or models being tested in several of the analyses were not clear to me. E.g., why should we expect microbiome-associated genes to encode proteins that are networked to one another? Is the relationship to the GF versus conventional comparison based on the idea that if genes affect the microbiome, then changes in the microbiome should also affect the proteins those genes encode (not obvious to me)? Ultimately, what does it mean to be a "promising" candidate gene (promising for what kind of application), and how is that defined by these enrichment analyses?

Rather than having readers try to fill in the rationale for these analyses, it would be helpful to clarify the motivation in the paper itself-or even shorten these sections to focus on those where there is the clearest underlying model or hypothesis. This might also avoid diluting the more interesting, less speculative findings presented earlier.

We thank the reviewer for this constructive input. We now added and/or re-worded text at the beginning of each respective subsection to more explicitly state the motivation/hypothesis behind these additional analyses, which we believe are ultimately important to draw higher-level conclusions/generate novel hypotheses to serve as a basis of future work (lines 379-383, lines 444-445, lines 476-480, lines 510-514).

Reviewer #2 (Recommendations for the authors):

Overall, this was an excellent study and I have few comments for improvement. A couple of minor points:

The association between heritability of taxa and their co-speciation rates is interesting, but I wonder whether this could be explained by an association between abundance and statistical power to detect heritability/co-speciation. Have the authors considered modeling heritability as a function of co-speciation rate and relative abundance? Such an analysis would be important to determine whether heritability and co-speciation rate are associated independently of abundance.

Thank you for this suggestion. We used a multiple linear regression model to test if the cospeciation rate and the median relative abundance significantly predict the heritability estimate. The overall regression model was not significant (R2=0.27, F(2,17)=3.202, P=.067). Further, the cospeciation rate significantly predicts the heritability estimate (P=.022), while the median abundance does not (P=.92). We thus conclude that the heritability estimates and the cospeciation rates are associated independent of the abundance of a taxon in question (lines 209-215).

Relative abundances of microbes were used for all analyses (as opposed to absolute abundance quantification). Because findings based on relative abundances can be difficult to interpret, more discussion on this limitation would be helpful.

We thank the reviewer for pointing out this important aspect. We added a specific passage in the Discussion section (lines 670-675):

“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 (Vandeputte et al., 2017; Barlow et al., 2020). Thus, future studies incorporating absolute abundance estimates may improve the detection of host-microbe interactions.”

My understanding is that heritability estimates were generated only for ASV relative abundances. However, previous heritability studies have examined taxonomic levels from ASV to phylum. It would be interesting to extend the analyses presented to higher taxonomic levels beyond ASVs (eg species to phylum).

We calculated heritability estimates with the updated model for all taxonomic levels (see reply to the comments of Reviewers 1 and 3). The labels on Figure 1A-B of the article are colored by taxonomic level.

Reviewer #3 (Recommendations for the authors):

The main point of concern for me are the heritability estimates. The paper reports heritability results that may be a bit surprising considering the current knowledge and literature: the heritability values are very high, with several values around 90% or higher. This is unexpected, and I am not sure how this is reconciled with the expectation that most variation in the microbiome is environmental rather than genetic. Although some potential reasons are given in the Discussion (mice raised in a controlled environment, using cecal content, etc), it still makes me a bit uneasy to see such high heritability estimates. One potential way to approach this is to try a different statistical approach for calculating heritability. Another would be to compare heritability estimates from this study with estimates from other studies – by now there are quite a few studies that report microbiome-wide heritability estimates from humans, mice, and other host species. It could be useful to correlate the heritability estimates in the current study with those from these studies. Lastly, it would be good to compare microbiome heritability estimates from this study to heritability’s of other complex traits in the same system – are there other known phenotypes that have heritability estimates that are this high?

We agree with the Reviewer that it is important to critically evaluate the robustness of heritabilities. Accordingly, we calculated the heritability estimates using gaston (Hervé and Claire, 2020), sommer (Covarrubias-Pazaran, 2018), and GEMMA (Zhou and Stephens, 2012), as well as lme4QTL (Ziyatdinov et al., 2018). To compare the heritability estimates between the different methods, we used the same model where the first genotype PC is included as a fixed effect and the genetic relatedness matrix (GRM) is included as a random effect, as GEMMA does not allow categorical covariates such as mating pair ID (Figure 9). We used identical GRMs in lme4QTL and GEMMA. However, sommer and gaston incorporate estimation of GRMs when calculating heritabilities.

Heritability estimates from each of the alternative methods are significantly correlated with estimates from lme4QTL (gaston R=0.52; GEMMA R=0.93; sommer R=0.79; P < 2.2 e-16, Spearman's correlation; Author response image 5). However, the heritability estimates calculated by lme4QTL are consistently higher than estimates from other methods. In the revised manuscript, we use lme4QTL to calculate the heritability’s with the mating pair nested in the sub-cross ID as random effects to control for maternal and population structure (see reply to comment 4 from Reviewer 1; Author response image 4).

Author response image 5
Heritability estimates calculated with lme4QTL compared to other methods.

The dashed line represents the identity line with a slope of 1. The colored lines are the corresponding linear regression lines.

This results in fewer taxa with significant heritability’s, and estimated values are lower. Importantly, the correlation of the cospeciation rate with the heritability estimates remains significant (P = .037, R=0.47)

Next, we compared our heritability estimates with other studies. Other microbiome QTL/GWAS studies in humans and mice also reported high heritability estimates for some taxa: 82 % in (O’Connor et al., 2014), 92% in (Org et al., 2015), and 99% in (Hughes et al., 2020). We calculated the correlation of our heritability estimates with seven human studies (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), two mouse studies (O’Connor et al., 2014; Org et al., 2015), one pig study (Chen et al., 2018), and one primate study (Grieneisen et al., 2021). We found our DNA-based heritability estimates to be positively correlated with DNA-based heritability estimates of male mice in (Org et al., 2015) (R=0.6047, P=.04872, n=11) and with the heritability estimates from one human study by (Turpin et al., 2016) (R=0.3760, P=.04867, n=28).

Finally, we estimated the heritability of other traits in our mapping population; estimates for body weight (87%) and body length (67%) are also high and comparable to previous studies in mice (Keenan et al., 2021). This has been added to results (lines 158-159, 193-200).

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

All reviewers appreciate the thorough revisions and re-analyses, which have resulted in a much improved manuscript. However, there are a few remaining issues that need to be addressed, as outlined below (the first item is the most essential):

1) Address the concern of Reviewer 3 below, which concern inflated heritability values in the main text. While the correlations in h2 estimates between methods are reassuringly high, Figure 8 in the response to reviewers also shows that the h2 reported in the paper (from lme4QTL) are typically several-fold (and sometimes an order of magnitude) higher than obtained from three alternative methods. This information would not be available to readers in the manuscript's current form. The response to reviewers indicates that lme4qtl was chosen because it had a high R2 with the PVE explained by the additive effects of significant SNPs (Figure 4 in the response to reviewers). However, both GEMMA and sommer produce similar R2 values (and in one case, a higher R2 value) but much lower heritability estimates. Do these lower estimates also correlate with co-speciation rates?

Thank you for the further input to improve our heritability results. After careful inspection, we determined that the analyses described in the previous reply (performed in response to original comments from Reviewer 1) were not appropriate. Including principle components of the SNP-based relatedness matrix as fixed effects in the model only makes sense if the kinship matrix included as a random effect is pedigree-based. This is not necessary because we include the full SNP-based relatedness matrix as a random effect in the model following the approach in GEMMA, because “using a lower-rank relatedness matrix compromises the ability of the linear mixed model to control for sample structure.” (Zhou and Stephens, 2012, Supplementary text).

We recalculated heritability estimates in GEMMA without including PCs, and address differences in heritability estimates in detail below, in response to Reviewer 3.

The previous RNA-based heritability estimates from GEMMA and sommer do significantly correlate with co-speciation rates (Author response images 6-8). In the revised manuscript, we report that both chip heritability (newly calculated as described below) and narrow-sense heritability estimates for RNA are correlated with co-speciation rates (lines 142 – 161 and 177-182).

Author response image 6
Correlation of heritability estimates, calculated with GEMMA using the first principal component as a fixed effect, with co-speciation rates reported in Groussin et al., 2017.
Author response image 7
Correlation of heritability estimates, calculated with GEMMA using the first three principal component as a fixed effect, with co-speciation rates reported in Groussin et al., 2017.
Author response image 8
Correlation of heritability estimates, calculated with sommer using the first principal component as a fixed effect, with co-speciation rates reported in Groussin et al., 2017.

2) Consider the polishing revisions to the code repository suggested by Reviewer 3.

Please see reply below.

3) Consider integrating the useful explanation in the response to reviewers about experimental design controls for technical and batch effects in the main manuscript, as future readers may have similar questions.

Thank you for this advice. We added the explanation in the respective Methods sections (Sample collection, lines 614-616; DNA extraction and sequencing, lines 624-625 and 634-639; 16S rRNA gene analysis lines 643-655).

Reviewer #3 (Recommendations for the authors):

I thank the authors for a comprehensive revision that addressed many of my concerns. It was good to see that the heritability values are correlated across different methods, and that several other methods produced heritability values smaller than those generated by lme4QTL. This, together with the fact that heritability values for some bacteria are higher than those for traits like length and weight, supports my notion that lme4QTL heritability values are overestimated. I am not sure what the reason is -- it's hard to know without spending time digging into the data, and I might not have the statistics background to advise on this -- but I am not confident that these values are robust. I would encourage the authors to investigate these analyses very thoroughly, making sure that all potential confounders are accounted for, the models are reasonable, and there are no artifacts in the data. I would suggest the text includes the results of heritability analysis with other approaches (in addition to lme4QTL) more prominently: Figure 1 should report and visualize heritability values from all methods used, and visualize the correlations between them. The text should describe these results, the methods used, the heritability values reported and their correlation.

We agree the differences between estimates from different methods needed further explanation. We thoroughly scrutinized the publications describing each method, together with input from a colleague who is an expert in quantitative genetics (Jason Wolf, University of Bath).

The differences in estimates are due to the way in which the genetic relatedness matrix (GRM) is calculated. Several methods use standardized GRMs, which standardize the variance across SNPs by accounting for allele frequencies (see manuscript methods for GEMMA approach, lines 685-713; sommer calculates the GRM following (VanRaden, 2008)). This approach assumes causative alleles are rare in the mapping population, as might be expected when mapping disease traits in humans or mapping in a sample from a natural population. As our mapping population was produced by intercrossing eight inbred strains, causative alleles are not expected to be rare, even for deleterious phenotypes. In this case, GEMMA recommends using a centered GRM for mapping (GEMMA manual 4.4.2), which we followed. However, when calculating heritability estimates, GEMMA always accounts for allele frequency differences – by multiplying the genetic variance component by a standardization factor to (Zhou et al., 2013, supplementary text). This is because the ‘chip heritability’ calculated in GEMMA only estimates the heritability that can be explained by genotyped SNPs (Zhou and Stephens, 2012) – in contrast to narrow-sense heritability (h2), which is the proportion of phenotypic variance explained by all additive genetic effects. Chip heritability is useful for genotype prediction, such as in animal breeding and predicting disease risk in humans. Typically chip heritability values are much lower than narrow-sense heritability estimates from other methods (e.g. for height in humans h2 ~0.8, chip heritability is 0.41; Zhou et al., 2013). This explains why the heritability estimates from GEMMA (and other methods using GRMs adjusting for allele frequencies) were lower than estimates from lme4QTL (narrow-sense heritability).

We are interested in estimates of narrow-sense heritability, which is the upper bound of what could be explained if all additive effects are known. However, we decided to report both chip heritability and h2 in the manuscript (Results lines 143 – 161, Figure 1A-B, Supplementary File 1), for ease of comparison with some previous studies, which report estimates of chip heritability.

It is not possible to include categorical covariates in GEMMA, therefore we used lme4QTL to calculate both chip heritability and h2 (Methods, lines 685-708), after first verifying that using the standardized GRM without covariates in lme4QTL produces the same results as GEMMA.

In addition, we added a comparison of h2 estimates to PVE by significant SNPs in Results lines 279 – 287, Figure 1A-B.

There should be a clear discussion about the possible reasons for the high heritability values (and why they are lower using other methods).

We added an explanation of the difference between chip heritability and narrow-sense heritability estimates to the Discussion (lines 461-468).

I would also suggest including the analysis comparing heritability values across studies in the text, and include a visualization of these correlations, rather than just reporting the p-value.

We agree. We added a plot (Figure 1—figure supplement 5) showing the correlation of the narrow-sense heritability estimates of this study (lme4QTL) with Org et al., 2015 and Turpin et al., 2016. Supplementary File 1 includes heritability estimates from previous studies, along with our estimates.

Regarding code availability, I want to thank the authors for enhancing the README file on the github repository, which now provides a nice description of the pipeline and analysis steps. However, I am not sure if this is sufficient for readers who want to reproduce the results: looking at the code itself, it seems like there are commands to load scripts that are not included in the repository (e.g. the snp_heritability_lme4qtl.R script loads the script function_for_gemma.r that I couldn't find anywhere), and these scripts might not be able to be run on other machines. I recommend amending the github repository and scripts so that anyone who wishes to do so is able to run the analysis and reproduce the results.

Thank you for spotting that this script was missing. It has now been added. We however respectfully maintain that further development of our documented pipeline would amount to the creation of a new analysis package, which we feel would go beyond the scope of the current study. The github repository was rather intended as a basis for researchers to adapt these approaches to their individual datasets.

References:

VanRaden, PM (2008), ‘Efficient methods to compute genomic predictions.’, J Dairy Sci, 91 (11), 4414-23.

Zhou, Xiang and Matthew Stephens (2012), ‘Genome-wide efficient mixed-model analysis for association studies’, Nat. Genet., 44 (7), 821-24.

Zhou, Xiang, Peter Carbonetto, and Matthew Stephens (2013), ‘Polygenic modeling with Bayesian sparse linear mixed models’, PLoS Genet., 9 (2), e1003264.

https://doi.org/10.7554/eLife.75419.sa2

Article and author information

Author details

  1. Shauni Doms

    1. Max Planck Institute for Evolutionary Biology, Plön, Germany
    2. Section of Evolutionary Medicine, Institute for Experimental Medicine, Kiel University, Kiel, Germany
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing – review and editing, Writing – original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2129-4138
  2. Hanna Fokt

    1. Max Planck Institute for Evolutionary Biology, Plön, Germany
    2. Section of Evolutionary Medicine, Institute for Experimental Medicine, Kiel University, Kiel, Germany
    Contribution
    Resources
    Competing interests
    No competing interests declared
  3. Malte Christoph Rühlemann

    1. Institute for Clinical Molecular Biology (IKMB), Kiel University, Kiel, Germany
    2. Institute for Medical Microbiology and Hospital Epidemiology, Hannover Medical School, Hannover, Germany
    Contribution
    Formal analysis, Software, Validation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0685-0052
  4. Cecilia J Chung

    1. Max Planck Institute for Evolutionary Biology, Plön, Germany
    2. Section of Evolutionary Medicine, Institute for Experimental Medicine, Kiel University, Kiel, Germany
    Contribution
    Resources
    Competing interests
    No competing interests declared
  5. Axel Kuenstner

    Institute of Experimental Dermatology, University of Lübeck, Lübeck, Germany
    Contribution
    Formal analysis, Methodology, Validation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0692-2105
  6. Saleh M Ibrahim

    1. Institute of Experimental Dermatology, University of Lübeck, Lübeck, Germany
    2. Sharjah Institute of Medical Research, Sharjah, United Arab Emirates
    Contribution
    Conceptualization, Methodology
    Competing interests
    No competing interests declared
  7. Andre Franke

    Institute for Clinical Molecular Biology (IKMB), Kiel University, Kiel, Germany
    Contribution
    Conceptualization, Funding acquisition, Supervision
    Competing interests
    No competing interests declared
  8. Leslie M Turner

    Milner Centre for Evolution, Department of Biology & Biochemistry, University of Bath, Bath, United Kingdom
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Supervision, Validation, Visualization, Writing – review and editing, Resources
    Contributed equally with
    John F Baines
    For correspondence
    l.m.turner@bath.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5105-3546
  9. John F Baines

    1. Max Planck Institute for Evolutionary Biology, Plön, Germany
    2. Section of Evolutionary Medicine, Institute for Experimental Medicine, Kiel University, Kiel, Germany
    Contribution
    Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – review and editing, Investigation, Resources
    Contributed equally with
    Leslie M Turner
    For correspondence
    Baines@evolbio.mpg.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8132-4909

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).

Senior Editor

  1. Wendy S Garrett, Harvard T.H. Chan School of Public Health, United States

Reviewing Editor

  1. Jenny Tung, Duke University, United States

Publication history

  1. Preprint posted: September 28, 2021 (view preprint)
  2. Received: November 9, 2021
  3. Accepted: June 14, 2022
  4. Version of Record published: July 19, 2022 (version 1)

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

  • 248
    Page views
  • 66
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. Shauni Doms
  2. Hanna Fokt
  3. Malte Christoph Rühlemann
  4. Cecilia J Chung
  5. Axel Kuenstner
  6. Saleh M Ibrahim
  7. Andre Franke
  8. Leslie M Turner
  9. John F Baines
(2022)
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
eLife 11:e75419.
https://doi.org/10.7554/eLife.75419

Further reading

    1. Evolutionary Biology
    2. Microbiology and Infectious Disease
    Pramod K Jangir et al.
    Research Article

    Bacterial pathogens show high levels of chromosomal genetic diversity, but the influence of this diversity on the evolution of antibiotic resistance by plasmid acquisition remains unclear. Here, we address this problem in the context of colistin, a ‘last line of defence’ antibiotic. Using experimental evolution, we show that a plasmid carrying the MCR-1 colistin resistance gene dramatically increases the ability of Escherichia coli to evolve high-level colistin resistance by acquiring mutations in lpxC, an essential chromosomal gene involved in lipopolysaccharide biosynthesis. Crucially, lpxC mutations increase colistin resistance in the presence of the MCR-1 gene, but decrease the resistance of wild-type cells, revealing positive sign epistasis for antibiotic resistance between the chromosomal mutations and a mobile resistance gene. Analysis of public genomic datasets shows that lpxC polymorphisms are common in pathogenic E. coli, including those carrying MCR-1, highlighting the clinical relevance of this interaction. Importantly, lpxC diversity is high in pathogenic E. coli from regions with no history of MCR-1 acquisition, suggesting that pre-existing lpxC polymorphisms potentiated the evolution of high-level colistin resistance by MCR-1 acquisition. More broadly, these findings highlight the importance of standing genetic variation and plasmid/chromosomal interactions in the evolutionary dynamics of antibiotic resistance.

    1. Evolutionary Biology
    Milo S Johnson, Michael M Desai
    Research Article Updated

    As an adapting population traverses the fitness landscape, its local neighborhood (i.e., the collection of fitness effects of single-step mutations) can change shape because of interactions with mutations acquired during evolution. These changes to the distribution of fitness effects can affect both the rate of adaptation and the accumulation of deleterious mutations. However, while numerous models of fitness landscapes have been proposed in the literature, empirical data on how this distribution changes during evolution remains limited. In this study, we directly measure how the fitness landscape neighborhood changes during laboratory adaptation. Using a barcode-based mutagenesis system, we measure the fitness effects of 91 specific gene disruption mutations in genetic backgrounds spanning 8000–10,000 generations of evolution in two constant environments. We find that the mean of the distribution of fitness effects decreases in one environment, indicating a reduction in mutational robustness, but does not change in the other. We show that these distribution-level patterns result from differences in the relative frequency of certain patterns of epistasis at the level of individual mutations, including fitness-correlated and idiosyncratic epistasis.