Does MHC heterozygosity influence microbiota form and function?

MHC molecules are essential for the adaptive immune response, and they are the most polymorphic genetic loci in vertebrates. Extreme genetic variation at these loci is paradoxical given their central importance to host health. Classic models of MHC gene evolution center on antagonistic host-pathogen interactions to promote gene diversification and allelic diversity in host populations. However, all multicellular organisms are persistently colonized by their microbiota that perform essential metabolic functions for their host and protect from infection. Here, we provide data to support the hypothesis that MHC heterozygote advantage (a main force of selection thought to drive MHC gene evolution), may operate by enhancing fitness advantages conferred by the host's microbiome. We utilized fecal 16S rRNA gene sequences and their predicted metagenome datasets collected from multiple MHC congenic homozygote and heterozygote mouse strains to describe the influence of MHC heterozygosity on microbiome form and function. We find that in contrast to homozygosity at MHC loci, MHC heterozygosity promotes functional diversification of the microbiome, enhances microbial network connectivity, and results in enrichment for a variety of microbial functions that are positively associated with host fitness. We demonstrate that taxonomic and functional diversity of the microbiome is positively correlated in MHC heterozygote but not homozygote animals, suggesting that heterozygote microbiomes are more functionally adaptive under similar environmental conditions than homozygote microbiomes. Our data complement previous observations on the role of MHC polymorphism in sculpting microbiota composition, but also provide functional insights into how MHC heterozygosity may enhance host health by modulating microbiome form and function. We also provide evidence to support that MHC heterozygosity limits functional redundancy among commensal microbes and may enhance the metabolic versatility of their microbiome. Results from our analyses yield multiple testable predictions regarding the role of MHC heterozygosity on the microbiome that will help guide future research in the area of MHC-microbiome interactions.

Introduction

MHC molecules display protein antigens (peptides) on the surface of cells and are central to the development of the vertebrate adaptive immune response. They are also important signaling molecules that coordinate immune cell development and effector function[1-6]. MHC genes are the most polymorphic loci known in vertebrates, and much of this polymorphism lies in the exons encoding the peptide-binding cleft of MHC molecules or the promoter regions of MHC genes[7-10]. MHC class I molecules are expressed on all nucleated cells and present intracellularly-derived peptide antigens to circulating cytotoxic T cells. These molecules are central to anti-viral and anti-tumor immunity. MHC class II molecules (MHCII) have more restricted expression and present extracellularly-derived peptides to circulating helper T cells. These molecules are central to humoral (i.e. antibody-mediated) immunity against extracellular microbes (primarily bacteria).

Several hypotheses have been put forth to explain the force of selection driving the evolution and maintenance of high allelic diversity at MHC loci and these can broadly be separated into two non-mutually exclusive mechanisms of diversifying selection; diversifying selection driven by antagonistic host-pathogen interactions and diversifying selection driven by inbreeding avoidance[8, 11, 12]. Heterozygote advantage (overdominance) is a pathogen-centric model of MHC evolution that is based on the argument that MHC heterozygotes have a more effective immune response than MHC homozygotes and consequently heightened resistance to infection. This phenomenon has been supported using epidemiological data from humans infected with viruses including Hepatitis B[13], HIV[14, 15], and HTLV[16]. Experimental infection studies in mice[17, 18] and studies across a diverse array of wild vertebrate species[19-23] also support the argument that MHC heterozygosity can enhance resistance to infectious disease. Because MHC genes are co-dominantly expressed enhanced resistance can arise because the immune system can mount responses against a more diverse array of pathogen epitopes, and because deleterious MHC alleles can be masked in the heterozygous condition[11, 18].

Hypotheses regarding the process of MHC gene evolution have focused on the fitness consequences associated with infection by pathogens. However, multicellular organisms are also stably colonized throughout their lives by a vast array of microbial species (collectively termed the ‘microbiota’) that provide essential metabolic services to their host[24] and enhance resistance to infection[25]. Moreover, previous empirical studies in animal models have demonstrated an association between MHC polymorphism and microbiota composition. Initial studies used liquid gas chromatography to compare short-chain fatty acid (SCFA) profiles (SCFAs are exclusively produced by bacterial symbionts in the gut) demonstrated differences in the composition of fecal flora among MHC congenic mouse strains[26, 27]. A subsequent study using bacterial 16S rRNA gene sequencing in stickleback fish, which are a classic model for studying MHC polymorphism, demonstrated associations between MHC alleles and microbial composition in the gut[28]. Importantly, this study was also the first to report an effect of MHC heterozygosity on the relative abundance of specific microbial species in the gut. A more recent study by Kubinak et al. (2015) used 16S rRNA sequencing to demonstrate that MHCII-/- mice had a distinct microbial community form in their gut compared to WT C57BL/6 control animals, that microbial communities were taxonomically distinct among MHC congenic mouse strains reared under a highly controlled environment, and that the effect of MHC on microbiota composition was most pronounced on mucosally-associated gut microbial communities[1].

MHC molecules regulate vertebrate adaptive immunity and are encoded by the most polymorphic genes known. Is MHC polymorphism favored because it helps establish and maintain a physiologically beneficial microbiome? Here, we provide a dataset that addresses the hypothesis that MHC heterozygosity influences taxonomic composition and functional gene content of the microbiome. Here, we aimed at exploring the link between the MHC-modulated microbiota and host fitness by comparing predicted metabolic pathways. We show that MHC heterozygotes develop microbiomes that are enriched in microbial functions positively associated with host health. Results described here suggest that MHC heterozygosity may confer a fitness advantage by promoting the development of gut microbiome that is less invasive/pro-inflammatory, more metabolically capable, more resistant to invasion by pathogens, and/or more functionally adaptive.

Materials and methods

Ethics statement

All mice were reared at the University of Utah School of Medicine, and all mouse work adhered strictly to University and Federal guidelines regarding the use of animals in research. Animals were euthanized via CO2 asphyxiation. All animal work was approved by the University of Utah Institutional Animal Care and Use Committee (Protocol#14–05009).

Animals

In this study, we utilized thirty 8–12 week old female BALB/c MHC congenic mice representative of three MHC homozygote (H2, H2, H2) and three MHC heterozygote (H2, H2, H2) genotypes (n = 5 female animals per genotype). We chose to limit our analysis to females due to the effect of sex on microbiota composition[29]. Multiple generations of each genotype used in this study were reared and maintained in the same facility under standardized environmentally-controlled conditions (IVC caging, automatic water, fed a standardized soy-free irradiated mouse chow (Envigo 2920X), 12:12 light:dark cycles). More specifically, in an effort to minimize the confounding effect of strain isolation, MHC homozygote animal purchased from Jackson Laboratories were used as founders to produce MHC heterozygote F1 progeny. F1 progeny were used to re-derive MHC homozygote F2 progeny. F2 MHC homozygote progeny were then used to produce the F3 progeny MHC homozygote and heterozygote animals used for this study (S1 Fig). To equalize cage effect, five female mice from each genotype were derived from a single breeding cage and came from the same litter. Thus, fifteen homozygote mice were derived from three separate cages and fifteen heterozygote mice were derived from three separate cages. Therefore, cage effect should not be a primary driver of observed differences between MHC homozygote and heterozygote mice. All mice were reared at the University of Utah School of Medicine, and all mouse work adhered strictly to University and Federal guidelines regarding the use of animals in research.

16S rRNA sequencing

Fecal samples were collected, immediately flash-frozen in liquid nitrogen, and stored at -80°C until downstream processing. DNA was extracted from fecal pellets using a bead-beating method following kit instructions (MoBio PowerFecal DNA Extraction Kit). Primers targeting the V3/V4 region of the bacterial 16S rRNA gene were used to generate PCR products for use in library prep. A detailed description of library preparation (including primer description) can be found in Kubinak et al.[1]. After quality filtering, all samples were rarified to a depth of 16,000 high quality sequence reads per sample for use in all of the analyses contained in this manuscript.

Sequence processing and analysis

QIIME 1.8.0 software was used for the upstream and downstream processing of the raw 16S rRNA gene sequences36. Followed by primer truncating and quality filtering, the sequences were assigned to operational taxonomic units (OTUs) at a minimum of 97% sequence similarity using closed reference OTU-picking method. The upstream processing steps include the removal of sequences with anonymous bases, chimera sequences using vsearch method82. In addition, phred quality score of at least 30 was used to maintain high quality sequences. The GreenGenes database41 was used to cluster the quality-filtered reads using Uclust algorithm40, and then RDP classifier42 was used to assign taxonomy. Before analysis, we performed following filtering strategies on the OTU-table: (1) singleton OTUs were removed to reduce the possibility of sequencing errors, (2) all samples were rarified to 16,000 sequences to avoid sequence effort bias. PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) analysis, which is a predictive metagenomics analysis that utilizes 16S rRNA sequencing data, was used to infer functional gene content37. A Bray-Curtis (BC) metric was used for all beta-diversity analyses comparing community diversity between MHC homozygotes and MHC heterozygotes. Phylogenetic (UniFrac) comparisons among genotypes are provided as supplemental data (S2 Fig). We use the terms "overall" and "core" to refer to all observed OTUs among mice or the subset of OTUs shared among all (100%) of samples, respectively. Alpha diversities were calculated to compare the amount of diversity contained within MHC heterozygotes and MHC homozygotes using Shannon index, a taxonomic metric and phylogenetic distance, a phylogenetic metric. A richness metric was used to quantify the total number of observed functional genes.

Network construction

Networks for each genotype were constructed to identify correlations among members of bacterial communities, where relative OTU abundance data was used to calculate Spearman correlation coefficients between all possible pairs of OTUs across mice. OTU-pairwise correlations between these OTUs were performed on within-group data. OTU-pairs with correlation coefficient values of either >0.8 or <-0.8, and corrected P ≤ 0.05 (Benjamini-Hochberg procedure) were used in the construction of networks. All OTU-pairs that did not satisfy these criteria were discarded from the analysis. The analysis was performed in R environment using “multtest” package38, and networks were plotted using Cytoscape v3.2.139. Interacting OTUs within networks represent significant correlations across samples, where the positive correlations indicate the similar abundance pattern of the interacting OTUs, while negative correlations signify opposite abundance pattern. Following construction, we characterized and compared networks by calculating few of their topological features: (1) number of nodes, which represents the total number of OTUs in the network; (2) number of edges, which represents the total number of significant pairwise correlations within a network; (3) connectivity, which is a ratio between number of edges and number of nodes in a network; (4) betweenness centrality, which represents the importance of a node by estimating the number of times it acts as a connector along the shortest path between two other nodes within a network. Therefore, higher betweenness centrality value of a node is indicative of its higher involvement in the overall connectivity with other nodes in the network, and vice versa.

Statistics

We employed DESeq2 algorithm80, an R package, to identify the fecal microbial taxa and PICRUSt-predicted functional KEGG Orthology (KO) genes that are differentially abundant between homozygotes and heterozygotes. Compared to the abundances of the microbial taxa and predicted KO genes in the fecal microbial community of homozygotes, this test identified whether they are differentially abundant in the heterozygote communities. The differences of their abundances were reported as log2FoldChange, where features (i.e., taxa and KO genes) enriched in homozygote are indicated by positive values, and heterozygote enriched features have negative values. To compare the alpha diversities, and pairwise taxonomic distances between samples, a nonparametric two-sided two sample t-test was carried out using 999 Monte Carlo permutations. In these statistical tests, Benjamini-Hochberg procedure was used to adjust the P-values, and statistical significance was considered at adjusted P < 0.05. Error bars represent standard errors of mean (S.E.M). Mantel test was performed to test the correlation of Bray-Curtis distances between taxonomic and functional gene contents within homozygote and heterozygote communities, where 999 permutations were used to calculate the P-values. Analysis of similarity (ANOSIM) test was used for multivariate comparisons shown in all PCoA plots. Area under curve (AUC) analysis was performed to compare the maintenance of taxonomic and functional gene content across MHC homozygote and MHC heterozygote animals. The areas under the curves were calculated in GraphPad Prism version 7, and a two-tail P-value was estimated as described elsewhere81. Unless otherwise stated, all statistical analyses were performed in QIIME 1.8.0 and 1.9.0 versions.

Results

Diversity patterns of fecal bacterial community

16S rRNA microbial community profiling was performed to test whether MHC genotype impacts taxonomic composition of the microbiota. "Overall" (i.e. all observed OTUs) microbiota community diversity was significantly different between MHC homozygotes and heterozygotes (Fig 1A; ANOSIM: R = 0.23, P = 0.03). Overall microbiota composition was significantly more dissimilar among MHC heterozygous animals compared to homozygotes (Fig 1B; non-parametric two-sample t-test: Benjamini-Hochberg adjusted P < 0.01). There was no observed difference in alpha diversity measures (Shannon diversity (Fig 1C) and phylogenetic diversity (Fig 1D); adjusted P > 0.05 in both cases) between overall homozygote and heterozygote communities. Within the overall communities of 692 OTUs, we tested for differences in OTU abundance between samples in the homozygotes and heterozygotes, and estimated that 98 OTUs are differentially abundant in homozygous animals, whereas 91 OTUs are differentially abundant in heterozygous animals (S1 Table; DESeq2 test: Benjamini-Hochberg adjusted P < 0.05). Most of the differentially abundant taxa belong to the Firmicutes phylum.

Fig 1

MHC heterozygosity alters microbiota composition and drives community divergence among individuals.

A) PCoA of beta-diversity comparison of overall community taxonomic composition between MHC homozygote and MHC heterozygote animals. Bray-Curtis analysis, ANOSIM: R = 0.23; P = 0.03; no. of permutations = 99. B) Boxplots representing within-group taxonomic dissimilarity and taxonomic dissimilarity comparisons between MHC homozygote and MHC heterozygote overall community. Comparison of alpha diversities in overall communities between MHC homozygote and MHC heterozygote animals based on Shannon (C), and phylogenetic diversity (D) indices. E) PCoA of beta-diversity comparison of core community between MHC homozygote and MHC heterozygote animals. Bray-Curtis analysis, ANOSIM: R = 0.06; P = 0.09; no. of permutations = 99. F) Boxplots representing within-group taxonomic dissimilarity and taxonomic dissimilarity comparisons between MHC homozygote and MHC heterozygote core community. G) Comparison of alpha diversities in core communities between MHC homozygote and MHC heterozygote animals based on Shannon diversity index. H) Differential enrichment in abundance of specific taxa within the core microbiota of MHC homozygote and MHC heterozygote animals. The significance threshold in this plot was P <0.05, Benjamini-Hochberg corrected. X-axis shows the values of log2FoldChange, where positive values indicate the taxa enriched in homozygotes (black) and negative values indicate the taxa enriched in heterozygotes (grey). 1 = Firmicutes (Clostridiales); 2 = Actinobacteria; and 3 = Bacteroidetes. Symbols (*), (**), (***), and (****) indicate the Benjamini-Hochberg corrected significance values of P < 0.05, P < 0.01, P < 0.001, P < 0.0001 respectively, which were calculated using non-parametric two sample t-test.

The concept of a “core microbiota” is based on the idea that there is a subset of microbial species (or functions) within complex communities that play fundamental roles in ecosystem function. Their fundamental role in ecosystem function would predict that these species should disproportionately influence host fitness. Therefore, we wanted to focus a portion of our analyses on the influence of MHC heterozygosity on the functional capacity of the core microbiota. We strictly defined core microbiota membership as those species present in all of the mice used in this study. Of the 692 total OTUs identified in our analysis, approximately 12% (82 OTUs) were found in all mice. Predictably, there is no gross difference in taxonomic composition of the core microbiome between MHC heterozygote and MHC homozygote animals (Fig 1E; ANOSIM: R = 0.06, P = 0.09). However, we noticed, in comparison to the heterozygotes, that the core microbial community of homozygotes showed significantly lower pairwise taxonomic dissimilarity, indicating any two homozygote communities are more similar in taxonomic composition (Fig 1F; adjusted P < 0.0001). MHC homozygotes also had higher Shannon diversity, indicating that the relative abundance of core species was significantly less variable (uneven) among MHC homozygotes compared to heterozygotes (Fig 1G; adjusted P < 0.05). Within the core microbiota of MHC heterozygotes, specific species within the S24-7 family (phylum Bacteroidetes), Oscillospira genus, and Adlercreutzia genus are enriched in abundance (Fig 1H; DESeq2 test: Benjamini-Hochberg adjusted P < 0.05). Within the core microbiota of MHC homozygotes, specific species within the Lachnospriaceae family and Clostridiales Order were enriched in abundance (Fig 1H; DESeq2 test: Benjamini-Hochberg adjusted P < 0.05).

Diversity patterns of predicted functional gene profiles

We next performed PICRUSt analysis using 16S rRNA sequencing data as a method to infer whether functional gene diversity was influenced by MHC heterozygosity. In contrast to taxonomic diversity, the composition of functional genes contained within overall microbiota communities was not a discriminating factor between homozygotes and heterozygotes (Fig 2A; ANOSIM: R = 0.03, P = 0.23). Similarly, the total number of inferred functional genes was consistent between genotypes when considering the overall microbiome (Fig 2B; adjusted P > 0.05). The total number of observed functional genes were equivalent between overall and core microbiomes despite the core microbiota representing only 12% of the total observed species diversity observed (not shown). This result is due to bias in PICRUSt analysis that arises because quantification of inferred functions is based on an extremely limited number of known microbial genes. It is important to acknowledge this sampling bias here, because we do not want the reader to erroneously conclude that non-core species do not contribute novel functional genes. When comparing the degree of dissimilarity between MHC homozygote and MHC heterozygote microbial communities based on taxonomic and functional composition, we found that taxonomic dissimilarity between homozygotes and heterozygotes was significantly reduced when considering the core versus overall microbiota (Fig 2C; adjusted P < 0.0001). This makes sense based on how we have defined our core microbiota; the OTUs that are observed in all samples. Because the same species make up the core MHC homozygote and MHC heterozygote communities, this will decrease taxonomic dissimilarity between these communities. However, overall and core functional dissimilarity was the same between MHC homozygotes and heterozygotes despite exclusion of approximately 90% of the total observed species from the core microbiota. Again, this is explained by the sampling bias mentioned above and should not be interpreted as biologically meaningful. Finally, variability in gene content of the overall microbiome was significantly higher among MHC heterozygotes compared to homozygotes (Fig 2D; adjusted P < 0.0001).

Fig 2

MHC heterozygosity drives functional divergence among individuals.

A) PCoA of functional beta-diversity comparison between MHC homozygote and MHC heterozygote animals based on Bray-Curtis analysis of PICRUSt-annotated KEGG Orthology (KO) genes (ANOSIM: R = 0.03; P = 0.23; no. of permutations = 99). B) The number of functional genes is compared between MHC homozygote and MHC heterozygote animals. C) Boxplots representing within-group functional dissimilarity comparisons between MHC homozygotes and MHC heterozygotes based on Bray-Curtis analysis. D) Boxplots representing within-group predicted KO gene dissimilarity between MHC homozygote and MHC heterozygote animals based on Bray-Curtis analysis. Symbol (****) indicates Benjamini-Hochberg corrected significance values of P < 0.0001, which were calculated using non-parametric two sample t-test.

MHC heterozygosity is associated with less functional redundancy in the microbiome

We next quantified the relationship between taxonomic and functional diversity in the overall MHC homozygote and MHC heterozygote microbiomes. Using correlation analysis, we found a stronger linear relationship between taxonomic and functional diversity in MHC heterozygotes compared to MHC homozygotes (Fig 3A; Mantel test: r = 0.39, P < 0.01 for homozygotes; r = 0.71, P < 0.01 for heterozygotes). Moreover, when we considered the loss of functional diversity across individuals (i.e. functional decay) we also found that gene functions were more rapidly lost among MHC heterozygote than MHC homozygote animals (Fig 3B, right panel), while taxonomic diversity was lost at a similar rate between groups (Fig 3B, left panel). Collectively, these data suggest that the observed reduction in functional dissimilarity among MHC homozygote animals (Fig 2D) may be due to the fact that there is more functional redundancy among bacteria found in MHC homozygote compared to MHC heterozygote microbiomes.

Fig 3

MHC heterozygotes have a less functionally redundant microbiome.

A) Relationships of beta-diversity (Bray-Curtis) patterns between closed-referenced OTUs and predicted functional (KO) genes within MHC homozygotes and MHC heterozygotes. In this analysis, the beta diversity patterns are shown using the principal coordinates scores of the first axes of taxonomic and predicted KO gene datasets. Mantel tests conducted to estimate the significance of the correlation (Homozygotes: r = 0.39, P < 0.01; Heterozygotes: r = 0.71, P < 0.01). B) Decay plots illustrate the loss of taxonomic and functional diversity across MHC homozygote and MHC heterozygote animals. P-value in based on area under curve (AUC) analysis. Area under the curve of homozygote and heterozygote taxonomic datasets are estimated to be 6,759 and 6,974, respectively, where this difference is not statistically significant. On the other hand, AUCs of homozygote and heterozygote PICRUSt-annotated datasets are estimated to be 41,407 and 40,632, respectively, where this difference turns out to be significantly different (P < 0.01).

Genotype-level network analyses predict more complex interactions within heterozygote core bacterial communities

In an ecosystem, microorganisms do not remain isolated, they interact with each other and form a complex web of metabolic interactions. The functional capacity of a microbial community is influenced by the degree to which metabolically active species interact to form trophic cascades43,44. It is therefore assumed that microbial species that feed off of products produced by other species (syntrophic interactions), have positively correlated abundance patterns. Alternatively, when two microbes interact in an antagonistic manner (e.g. through the production of antimicrobials), they will have inversely correlated abundance patterns. Network analysis is a method used to characterize the general phenotype of species interactions in an ecosystem. We wanted to determine whether MHC heterozygosity could influence the structure of microbial networks. For this analysis, we again focused on the core microbiota for two reasons. First, these species are hypothesized to contribute most to the essential metabolic functions of the host. Second, peripheral species are transient (or rare) microbes that either may not form stable interactions with other species or may not provide essential functions to the host. Third, we wanted to reduce genotype-specific effects (i.e. the presence of bacterial species, and their associated functions, that are only found in certain homozygote or heterozygote genotypes) in order to focus on differences in more generalizable features of the microbiome across our six genotypes.

We constructed genotype-level networks (Fig 4A), and compared between them by their topological features (Fig 4A and 4B). The total number of nodes (i.e. number of species having one or more interactions), significant OTU-OTU pairs (r > 0.8 or r < -0.8, Benjamini-Hochberg adjusted P ≤ 0.05) (edges), and the average number of edges per node (connectivity) were all higher in MHC heterozygotes compared to homozygotes (Fig 4B). Later, we investigated node-wise feature in the networks by their betweenness centrality values, and identified few nodes (i.e., OTUs) in each network that, by definition, are most important in connecting other OTUs in their respective networks (Fig 4A). Each genotype has a set of these nodes, which, in general, have higher betweenness centrality values in heterozygotes compared to homozygotes (S2 Table), suggesting these nodes are more important in heterozygotes as they more frequently connect other nodes in heterozygote networks and contribute to forming more complex networks compared to homozygotes. For example, an unclassified member of Lachnospiraceae (OTU ID # 258202) is the most important node in both H2 and H2 networks, where the betweenness centrality values are 0.071, and 0.024, respectively. We also analyzed the relative proportion of positive and negative correlations in the networks. Positive interactions dominated in both MHC heterozygote as well as MHC homozygote core microbiota (~70% positive correlations versus ~30% negative), and there was no difference in this ratio between groups (Fig 4C). Taken together, higher connectivity within heterozygote networks might represent a more complex web of metabolic interactions within core bacterial community, where nodes with high centrality values are assumed to play major metabolic hubs. These data support higher diversity of functional (KO) genes in relation to taxonomic (OTU) diversity (Fig 3), suggesting MHC heterozygosity favors metabolic versatility of the microbiome.

Fig 4

MHC heterozygotes have a more connected microbial network in core communities.

A) Plots depict association between microbial species using Spearman correlation coefficient values of >0.8 for positive correlations, and <-0.8 for negative correlations with statistical significance P<0.01, Benjamini-Hochberg corrected. The microbial pairs that follow these criteria were considered to have strong association among them and included in network study. In the network analysis, nodes represent microbial taxa (OTUs), and lines (or edges) connecting the nodes represent significant pairwise correlations in species abundances. The size of the nodes is proportional to their betweenness centrality values. Nodes with highest betweenness centrality values (“hub” species; here we considered top seven OTUs across networks), along with their edges, are colored for comparative understanding of their connectedness across genotype-wise networks. These nodes are identified taxonomically in the accompanying legend. B) The comparative analysis of few network topologies between MHC homozygote and MHC heterozygote networks, number of nodes, number of edges, mean connection per node (i.e., ratio between number of edges and number of nodes in a network), and relative proportion of positive and negative correlations are illustrated.

Metabolic pathway analysis supports that MHC heterozygosity favors host-microbe symbiosis

Again, to focus on differences in microbiome functions generalizable across our six genotypes, we compared PICRUSt-predicted KEGG Orthology (KO) gene profiles of core microbial communities between homozygotes and heterozygotes and explored the alteration of specific functional genes and pathways caused by MHC heterozygosity. To check the accuracy of PICRUSt prediction, we calculated weighted Nearest Sequenced Taxon Index (weighted NSTI) scores for each sample that reflects the availability of sequenced genomes that are closely related to the OTUs in a given sample (S3 Table). A score of 0.2 indicates that functionality of microbial community in a given sample, on an average, can be predicted from reference genomes of similar species (80%). Since the close relatives of most of the core OTUs were not identified to the species level using the database used in our study (S2 Table), we obtained high NSTI scores, ranging 0.17 to 0.31.

Numerous differences were found in the relative abundance of inferred functional genes between MHC heterozygote and homozygote core microbiomes (Fig 5; S4 Table; DESeq2 test: Benjamini-Hochberg adjusted P < 0.05). In the core microbiome of MHC homozygotes, we observed increased abundance of genes associated with the metabolism of host-derived carbohydrates that comprise the mucus layer overlaying gut epithelium (alpha-L-fucosidase [EC:3.2.1.51], alpha-galactosidase [EC:3.2.1.22], alpha-mannosidase [EC:3.2.1.24], and beta-mannosidase [EC:3.2.1.25]). Additionally, MHC homozygotes have higher representation of predicted genes involved in the degradation of other glycans such as dietary fibers (N4-(beta-N-acetylglucosaminyl)-L-asparaginase [EC:3.5.1.26] and levanase [EC:3.2.1.65]). MHC homozygotes were also observed to be enriched for genes associated with D-lactate metabolism (D-lactate dehydrogenase [EC:1.1.1.28]), acetylaldehyde production (alcohol dehydrogenase (NADP+) [EC:1.1.1.2]), antimicrobial production such as streptomycin (Myo-inositol-1-phosphate synthase [EC:5.5.1.4]), bacterial toxin production (exfoliative toxin A/B), bacterial motility (Chemotaxis protein MotB), and oxalate synthesis (Formyl-CoA transferase [EC:2.8.3.16]).

Fig 5

MHC heterozygosity influences core microbiome function.

Major differentially abundant predicted KO genes encoding ECs involved in the KEGG functional categories between MHC homozygote and MHC heterozygote animals (see S4 Table for full list).X-axis shows the values of log2FoldChange, where positive values indicate the predicted KO genes enriched in MHC homozygotes, and negative values indicate the genes enriched in MHC heterozygotes. Bars are color coded by KEGG Level 2 functional groupings that are defined in the accompanying legend.

MHC heterozygote core microbiomes were enriched for genes associated with biosynthesis of essential amino acids like phenylalanine and tryptophan (3-phosphoshikimate 1-carboxyvinyltransferase [EC:2.5.1.19], 3-dehydroquinate synthase [EC:4.2.3.4], 3-deoxy-7-phosphoheptulonate synthase/chorismate mutase [EC:2.5.1.54 5.4.99.5], cinnamoyl-CoA:phenyllactate CoA-transferase [EC:2.8.3.17]), genes associated with the production and utilization of short-chain fatty acids (enoyl-CoA hydratase [EC:4.2.1.17], 2,4-dienoyl-CoA reductase (NADPH2) [EC:1.3.1.34], acetyl-CoA carboxylase [EC:6.4.1.2]), genes associated with the cellular respiration and ATP synthesis (pyruvate ferredoxin oxidoreductase, beta subunit [EC:1.2.7.1], succinate dehydrogenase flavoprotein subunit [EC:1.3.99.1], fumarate reductase flavoprotein subunit [EC:1.3.99.1], malate dehydrogenase (oxaloacetate-decarboxylating)(NADP+) [EC:1.1.1.40], pyrophosphatase PpaX [EC:3.6.1.1], flavoprotein-quinone oxidoreductase [EC:1.5.5.-], Protein NrfD; formate-dependent nitrate reductase complex), and genes associated with polyamine production (carboxynorspermidine decarboxylase [EC:4.1.1.-], putrescine carbamoyltransferase [EC:2.1.3.6]). Heterozygotes also have more genes associated with terpenoid backbone biosynthesis (1-deoxy-D-xylulose-5-phosphate synthase [EC:2.2.1.7], geranylgeranyl diphosphate synthase, type II [EC:2.5.1.1 2.5.1.10 2.5.1.29]).

Notably, we found one metabolic pathway that had opposing patterns of gene enrichment between MHC homozygotes and MHC heterozygotes, and it is associated with methane metabolism (Fig 6). MHC homozygotes were observed to be enriched for genes associated with methanogenesis (a byproduct of fermentation) (formylmethanofuran dehydrogenase subunit E [EC:1.2.99.5]), while MHC heterozygotes were enriched for genes associated with the utilization of methane as a carbon source (methanotrophy)(formate dehydrogenase, beta subunit [EC:1.2.1.2], formate dehydrogenase, gamma subunit [EC:1.2.1.2]).

Fig 6

MHC heterozygosity impacts methane metabolism in the gut.

Within the “core” microbial community of gut ecosystems, methanogenesis associated KO genes are enriched in MHC homozygotes, while MHC heterozygotes have higher abundance of KO genes associated with methanotrophic pathways. Green arrows indicate the genes enriched in MHC homozygote communities, while red arrows indicate the genes enriched in MHC heterozygote communities. Black arrows indicate either the observed genes are not differentially abundant between MHC homozygote or MHC heterozygote microbiomes (P > 0.05) or these genes are not included in PICRUSt database. The numbers in each step of the pathways indicates the EC number of associated enzymes.

Discussion

Several major observations from our sequencing study support the hypothesis that MHC heterozygosity significantly influences microbiome form and function. First, MHC heterozygotes have a unique taxonomic composition of their gut microbiota compared to that of MHC homozygotes. Second, while the microbiomes of MHC heterozygotes have similar numbers of total gene functions compared to MHC homozygotes, heterozygosity results in significant differential enrichment for certain functions relevant to host health; notably within the core microbiome that is hypothesized to be particularly important to host fitness. In particular, the core microbiomes of MHC heterozygotes appear to be enriched for beneficial functions and deficient in functions that have been associated with a variety of gastrointestinal diseases. MHC heterozygosity also results in greater individual variation in taxonomic and functional composition of the microbiome, which could be linked to selection favoring the addition of species that perform novel functions. All else being equal, novel functions are more likely to come from more phylogenetically distinct species. Interestingly, we also found that the core microbiome of MHC heterozygotes tends to reflect a more integrated network of species.

Our results complement other studies that have linked MHC polymorphism with taxonomic differences in microbiota composition[26-28]. Most current efforts to understand the physiologic role of the microbiota have relied on taxonomic compositional descriptions like this, but the functions that these species perform, are what impacts host fitness directly. Therefore, we wanted to infer the functional relevance of observed taxonomic differences between MHC heterozygotes and homozygotes. While taxonomic diversity drove significant differences in overall microbiota composition (Fig 1A), functional diversity did not (Fig 2A). MHC heterozygotes and homozygotes contained similar numbers of functional genes (Fig 2B). This finding is consistent with another study demonstrating largely overlapping functionality of the microbiomes of diverse mammal species, despite significant taxonomic divergence[30]. However, despite possessing largely overlapping functional profiles, we observed differential enrichment of specific functional genes between MHC homozygotes and heterozygotes (discussed below). Overall, several observations support the hypothesis that MHC heterozygotes may be more functionally versatile and metabolically capable.

The relationship between taxonomic and functional diversity within microbiomes is influenced by the degree of functional redundancy that exists among microbial species[31]. This is an important relationship because it determines the degree to which species diversity impacts host fitness and can be acted upon by natural selection. Functional redundancy can be adaptive by increasing the stability of a microbiome (i.e. chance of losing a critical function is reduced), or maladaptive by decreasing functional versatility of the host (i.e. the capacity to acquire novel functions). One of the most interesting results from our analysis was the significant difference in the correlation between taxonomic and functional diversities between MHC homozygotes and heterozygotes (Fig 3A and 3B). Microbiomes among MHC heterozygotes (beta diversity) were more variable in taxonomic composition (Fig 1B) and functional gene content (Fig 2D) compared to MHC homozygotes. This is in contrast to the fact that MHC heterozygotes did not have on average more taxonomic (Fig 1C and 1D) or functional diversity (Fig 2B) within their communities (alpha diversity). We can infer two points from this observation. First, this means that increased inter-individual variation observed in MHC heterozygotes is not simply a consequence of heterozygotes being colonized by a wider variety of species. Second, it suggests that whatever mechanism drives the effect of MHC heterozygosity on microbiota composition, it is a diversifying force of selection that favors the addition of novel functions as microbial species are added to a community. In contrast, that same force of selection appears to drive communities towards a taxonomically and functionally more homogenous composition in MHC homozygotes. Throughout life, microbes with novel functions are ingested through the act of feeding. If these bacteria confer enhanced ability to extract nutrients from a particular food source, then we would hypothesize that the ability to integrate novel functions into an existing microbiome is a fitness enhancing trait. Future studies measuring fitness differentials between MHC homozygote versus MHC heterozygote animals fed simple to increasingly complex diets, or fed the same diet versus a diverse array of diets will be useful for addressing this.

Microorganisms in the gut form a complex network of metabolic interactions. Understanding the processes that shape these interactions is important because they are deterministic factors governing the adaptability and resilience of microbiomes in response to environmental variability. The shape of these interactions (i.e. network topology) may also be a diagnostic feature of certain disease states. Positive interactions between microbial species occur through nutritional cross feeding; the synthesis of molecules that enhance the fitness of other microbes. In contrast, negative interactions occur through the synthesis of antimicrobial molecules and environmental modifications (nutritional or immunological modifications) that inhibit the growth of other microbes. Here, we have provided for the first time an analysis of the potential impact of MHC polymorphism on gut microbiome network topology. Our data suggests that MHC heterozygotes have more potentially interacting species, more total interactions between species, and more interactions per species within their core microbiomes compared to MHC homozygotes (Fig 4). Network modularity (i.e. species clustering) is also more complex in MHC heterozygotes versus homozygotes, with those species possessing the highest betweenness centrality scores (i.e. metabolic “hub” species) represented by OTUs, closest relatives of which are known beneficial microbial taxa such as Lachnospiraceae, Dorea, Oscillospira, and Ruminococcus (Fig 4A). Modularity is thought to enhance the adaptability of the microbiome through enhanced horizontal gene transfer between species[32] or by promoting the acquisition of novel functions by making the overall metabolic network more robust to the addition of novel functions[33, 34]. Modularity has been positively correlated with environmental variability. For commensal organisms in the gut, host factors as well as other bacteria or their metabolites can provide the necessary environmental variation to drive module formation. Given the highly controlled nature of the environment and diet these mice experienced, we believe variable host factors, and in particular antibodies) may facilitate module formation. Higher connectance and/or modularity is also thought to increase resistance of networks to biological invasion, which can be extrapolated to include pathogen colonization of the gut. These emergent properties of the microbial network could make MHC heterozygotes more resistant to enteric infection, which has been demonstrated previously in MHC heterozygous mice infected with Salmonella[1, 17, 18]. Finally, less connected and less modular microbiomes have been associated with disease states in humans including obesity, inflammatory bowel disease, and alcoholism[35-37].

Multiple gene functions were differentially abundant between MHC homozygotes and MHC heterozygotes. In an effort to describe how the effect of MHC heterozygosity on microbiome form and function might impact host fitness, we have focused our discussion on those genes whose functions are known to be important for host health. Carbohydrate metabolism is an essential function provided by members of the microbiota that is key for nutrient extraction from complex polysaccharides derived from plant tissues[38, 39]. However, the mucus lining of the gut epithelium is also covered in host-derived carbohydrates in the form of glycan, and this glycan layer is an important element of barrier defense in the gut that limits inflammation[40-43]. Multiple genes encoding enzymes that degrade host-derived glycans were enriched in MHC homozygotes compared to MHC heterozygotes (Fig 5). Microbial digestion of host-derived carbohydrates that comprise the mucus layer in the gut have recently been shown to enhance intestinal permeability and susceptibility to enteric infection in mice[44], and mucus-deficient mice develop spontaneous colitis[40]. Thus, barrier integrity could be reduced in MHC homozygotes compared to MHC heterozygotes because of reduced thickness of the mucus lining of the gut, and this may reduce fitness by pre-disposing individuals to acute (e.g. septicemia) and chronic (e.g. liver dysfunction) diseases. Additionally, MHC homozygotes have higher representation of predicted genes involved in the degradation of other glycans such as dietary fibers including N4-(beta-N-acetylglucosaminyl)-L-asparaginase [EC:3.5.1.26], and levanase [EC:3.2.1.65]. This observation probably explains why their gut microbiomes have higher abundance of genes involved in metabolism of other carbohydrates such as cellulose (cellulose synthase (UDP-forming) [EC:2.4.1.12]), sucrose (trehalose-6-phosphate hydrolase [EC:3.2.1.93]), and galactose (galactonate dehydratase [EC:4.2.1.6]). Higher potential access to carbohydrates from diet and host may favor a specific group of microorganisms in the colon with the ability to perform D-lactate fermentation, explaining why MHC homozygote microbiomes have increased D-lactate dehydrogenase gene. Increased accumulation of D-lactate in the blood causes D-lactate acidosis that can promote neurological disease and is linked to a gastrointestinal disorder called short-bowel syndrome[45, 46]. These data also imply that MHC homozygotes may be more susceptible to diseases associated with diets high in sugar like diabetes and colitis, which has been reported before in MHC congenic mouse models as well as humans[47-51].

Numerous bacterial metabolites are utilized by gut epithelial cells to enhance barrier integrity. Commensal-derived essential amino acids (EAAs) (e.g. tryptophan) have recently been shown to promote barrier integrity by upregulating IL-22 expression in gut epithelial cells[52]. In addition to the genes important for the biosynthesis of these EAAs, heterozygote microbiomes are enriched for genes associated with polyamine production (carboxynorspermidine decarboxylase [EC:4.1.1.-], putrescine carbamoyltransferase [EC:2.1.3.6]). Polyamines such as spermidine, putrescine are some of the most important molecules produced by the microbiota that are essential in regulating not only barrier integrity[53, 54], but also anti-inflammation, anti-mutagenicity, and autophagy[55]. In contrast, homozygotes are enriched in genes important for the synthesis of compounds known to decrease barrier integrity. For example, alcohol dehydrogenase (NADP+) [EC 1.1.1.71] catalyzes the conversion of ethanol to acetaldehydes, which disrupts the mucosal barrier in the intestine[56]. In addition, homozygote microbiomes have higher abundance of genes associated with the production of bacterial toxins (exfoliative toxin A/B) that might damage epithelial cells, and bacterial motility (chemotaxis protein MotB). Together, these observations suggest that MHC homozygotes may have a gut barrier deficiency that may be more permissive to leakage of pro-inflammatory microbial metabolites or whole bacteria from the gut into deeper host tissues.

MHC heterozygotes were enriched for genes whose functions are known to be important for dietary energy harvest and biosynthesis of essential nutrients. Indeed, twice as many genes associated with energy metabolism and amino acid synthesis were enriched in MHC heterozygotes compared to MHC homozygotes. MHC heterozygotes have higher abundance of microbial genes associated with terpenoid (also known as isoprenoid) biosynthesis via the methylerythritol phosphate (MEP) pathway (1-deoxy-D-xylulose-5-phosphate synthase [EC:2.2.1.7], geranylgeranyl diphosphate synthase, type II [EC:2.5.1.1 2.5.1.10 2.5.1.29]), which is important for the biosynthesis of vitamins essential for host health, including B1 (thiamine) and B6 (pyridoxine)[57]. The TCA cycle is an essential component of cellular respiration process through the oxidation of acetyl-CoA, with central importance in energy harvest, and biosynthesis of essential nutrients including amino acid, SCFAs. SCFAs are entirely produced by the microbiota[58, 59], which are particularly important for enterocyte health[60, 61], and have anti-carcinogenic and anti-inflammatory properties[62]. In addition, the microbiomes of MHC heterozygotes were enriched for genes that provide substrate for the TCA cycle (pyruvate ferredoxin oxidoreductase, beta subunit [EC:1.2.7.1]) as well as enzymes directly associated with the TCA pathway (succinate dehydrogenase flavoprotein subunit [EC:1.3.99.1], fumarate reductase flavoprotein subunit [EC:1.3.99.1] and malate dehydrogenase (oxaloacetate-decarboxylating)(NADP+) [EC:1.1.1.40]). Another striking difference between MHC homozygote and MHC heterozygote core microbiomes was the differential abundance of genes associated with methane metabolism (Fig 6). While MHC homozygote microbiomes had higher gene abundances associated with the synthesis of methane, heterozygote microbiomes had higher gene abundances for the utilization of it. KEGG pathways analysis showed that methanogens in the homozygotes employ reductive acetyl-CoA pathway in methanogenesis, whereas the methanotrophic pathway in MHC heterozygote contributes to the synthesis of acetyl-CoA and pyruvate. Production of methane is a major pathway for removing fermentation products including H2 from gut ecosystem[63], which may enhance fermentation efficiency such as D-lactate synthesis. In addition, previous studies showed the links between methanogenesis and a pathological increase of bacteria in the small intestine known as small intestinal bacterial overgrowth (SIBO)[64], which may contribute to irritable bowel syndrome (IBS)[65]. Two other important genes enriched in MHC heterozygotes are enoyl-CoA hydratase [EC:4.2.1.17], which catalyzes the synthesis of acetyl-CoA (and energy) from fatty acids, and acetyl-CoA carboxylase [EC:6.4.1.2], which catalyzes the conversion of acetyl-CoA into malonyl-CoA for fatty acid synthesis. While acetyl-CoA can be funneled into TCA cycles, malonyl-CoA is a precursor of fatty acid biosynthesis, suggesting that MHC heterozygote microbiomes may maintain a better tradeoff in fatty acid metabolism critical for host health.

Finally, Commensal microbes are important for the detoxification of xenobiotic compounds[66, 67], which could impact the host health in a variety of ways. The fecal microbiomes of MHC heterozygotes are enriched for genes whose functions are important for the detoxification of respiration-associated reactive oxygen species such as catalase, chloride peroxidase [EC:1.11.1.10], glutathione production (6-phosphogluconate dehydrogenase [EC:1.1.1.44], glucose-6-phosphate 1-dehydrogenase [EC:1.1.1.49]. These enzymes are important for scavenging/neutralizing free radicals, peroxides, and lipid peroxides formed during cellular respiration, suggesting that MHC heterozygotes may be less susceptible to diseases associated with oxidative stress.

It is important to acknowledge major caveats to the observations we describe in this study. The first stem from technological limitations of predictive meta-genomic approaches based on 16S rRNA sequence data. Many of the operational taxonomic units (OTUs) identified from gut ecosystem remain undefined based on our current knowledge. When defined, the biology of many of them is not fully understood, which challenges us to explore and compare the complete functional roles between any two systems. As a consequence of this, we cannot rule out that an observed functional advantage contributed by a microbial species found in heterozygotes cannot be rescued by a currently undefined microbial species found in MHC homozygotes. Additionally, PICRUSt annotated metabolic pathways are based on the KEGG reference database. This database contains all currently known microbial functions, but is only a small subset of the total possible functions (based on coding sequence). Therefore, many physiologically important functions will go unappreciated based on this analysis method. Also, because PICRUSt is based on 16S rRNA genes, eukaryotic microbial and viral contributions are ignored in this analysis. Finally, it is not within the scope of this manuscript to compare/contrast all of the unique features of the microbiota/microbiomes observed among our six genotypes. To simplify the story we have focused on core (i.e. shared) species and functions among MHC homozygotes and MHC heterozygotes and it must be explicitly restated that non-core species excluded from our analyses will contribute many novel functions to the microbiome that we do not have the resolution to discuss here. Finally, as stated above, more MHC homozygote and MHC heterozygote comparisons are needed to assess the generality of these effects. Given these limitations, we believe this study provides a framework for the generation of novel hypotheses regarding the influence of MHC-microbiota interactions on host physiology.

Given the role MHCII antigen presentation plays in the generation of antibody responses, we believe that MHC polymorphisms differentially regulate anti-commensal antibody responses to drive differences in community structure/function. If true, then MHC-mediated deficiencies in this response could result in unstable microbiomes that can enhance susceptibility microbiota-dependent diseases. Indeed, genes within the MHC, and MHCII specifically, have recently been strongly linked to IgA deficiency in humans[68-71] and IgA deficiency has been associated with increased susceptibility to a variety of autoimmune, inflammatory, and infectious diseases in humans. Additionally, it has also been shown that MHCII expression by group 3 innate lymphoid cells is important for inducing tolerance against commensal antigens[3, 72], and that non-classical MHC class I-like genes have been identified that can present unconventional antigens derived from the microbiota to influence immune development and microbiota composition. For example, a recent study has demonstrated that a non-classical MHC class I molecule (CD1d), that presents lipid antigens and is central to development and activation of natural killer-like T (NKT) cells in the gut[73], may also be important for sculpting microbiota composition[74]. Another non-classical MHC class I-like molecule, MR1, has been shown to bind and present B vitamins, which are exclusively produced by the microbiota, to drive the activation of mucosa-associated invariant T (MAIT) cells[75]. Thus, accumulating evidence suggests that classical MHC genes can operate through unconventional means to influence the fitness-consequences of host-microbiota interactions, and that non-classical MHC-like molecules represent novel mechanisms facilitating molecular crosstalk between the microbiota and various aspects of the mucosal immune system.

MHC polymorphisms can be favored evolutionarily if they promote the establishment of a beneficial microbiome. The benefits of our microbiomes arise primarily from the metabolic services and protection from infection they provide us. Results from our experiments provide support for the hypothesis that MHC heterozygosity is favored because it positively influences microbiome functionality. Future studies directly testing the physiologic predictions arising from our results are underway. Multiple studies in humans have revealed links between the presence of specific HLA alleles and susceptibility to autoimmune, infectious, and inflammatory diseases including rheumatoid arthritis[76], multiple sclerosis[77], inflammatory bowel disease[78-81], lupus[82, 83], and diabetes[84], in which the microbiome has been implicated as a contributing factor[85-87]. To what degree (and how) an individual’s MHC genotype contributes to molding the community of symbiotic microbes and how this influences host health are not understood. Collectively, our data complement the studies demonstrating that MHC polymorphisms can influence microbiota composition. However, we also make an effort to describe the impact of MHC heterozygosity on microbiome function, and provide preliminary evidence to support the hypothesis that MHC heterozygosity may be evolutionarily favored by promoting the establishment of a more beneficial microbiome. Given the life-long association between vertebrate hosts and their microbiota and the pervasive impact these symbionts have on host physiology, these interactions may have been important drivers of MHC gene evolution. In conclusion, we make no allusion that the data generated from six genotypes of inbred mice (maintained in ventilated cages, eating the same irradiated chow, housed under specific-pathogen-free conditions, etc) reflects a natural condition. The true biological relevance of such interactions need to be studied by Evolutionary Ecologists in the field. Excitingly, this has already begun. It is our hope that this descriptive/speculative study provokes conversation within such groups on the possibility that MHC polymorphism may regulate fitness through modulation of the microbiome.

MHC congenic BALB/c mice were purchased from Jackson Laboratories representing three different MHC genotypes (H2, H2, and H2 MHC genotypes).

These animals were used as founders to re-derive MHC homozygote and MHC heterozygote animals for use in this study.

(TIFF)

Click here for additional data file.

A) PCoA of beta-diversity comparison of overall community between MHC homozygote and MHC heterozygote animals based on Unweighted Unifrac analysis (ANOSIM: R = 0.18; P = 0.01; no. of permutations = 99). B) PCoA of beta-diversity comparison of core community between MHC homozygote and MHC heterozygote animals based on Weighted Unifrac analysis (ANOSIM: R = 0.06; P > 0.05; no. of permutations = 99).

(TIFF)

Click here for additional data file.

Differential enrichment of specific taxa within the overall microbiota of MHC homozygote and MHC heterozygote animals, which are represented as log2FoldChange values with P < 0.05, Benjamini-Hochberg adjusted.

The positive log2FoldChange values indicate that the taxa enriched in MHC homozygotes, while negative values indicate the taxa enriched in MHC heterozygotes.

(DOCX)

Click here for additional data file.

Distribution of betweenness centrality values of nodes with their closed relatives across genotype-level networks.

If an OTU is not calculated as a node in a network, its value is reported as N/A.

(XLSX)

Click here for additional data file.

Summary of the weighted Nearest Sequenced Taxon Index (weighted NSTI) scores of each of the samples used in this study.

The score represents the extent to which OTUs in a given sample are related to their sequenced genomes.

(DOCX)

Click here for additional data file.

Differential enrichment of predicted KO genes involved in the KEGG metabolic pathways within the core microbiota of MHC homozygote and MHC heterozygote animals, which are represented as log2FoldChange values with P < 0.05, Benjamini-Hochberg adjusted.

The positive log2FoldChange values indicate that the taxa are enriched in MHC homozygotes, while negative values indicate that the taxa are enriched in MHC heterozygotes.

(XLSX)

Click here for additional data file.