Whole exome sequencing analyses reveal gene-microbiota interactions in the context of IBD

Abstract

Objective: Both the gut microbiome and host genetics are known to play significant roles in the pathogenesis of IBD. However, the interaction between these two factors and its implications in the aetiology of IBD remain underexplored. Here, we report on the influence of host genetics on the gut microbiome in IBD.

Design: To evaluate the impact of host genetics on the gut microbiota of patients with IBD, we combined whole exome sequencing of the host genome and whole genome shotgun sequencing of 1464 faecal samples from 525 patients with IBD and 939 population-based controls. We followed a four-step analysis: (1) exome-wide microbial quantitative trait loci (mbQTL) analyses, (2) a targeted approach focusing on IBD-associated genomic regions and protein truncating variants (PTVs, minor allele frequency (MAF) >5%), (3) gene-based burden tests on PTVs with MAF <5% and exome copy number variations (CNVs) with site frequency <1%, (4) joint analysis of both cohorts to identify the interactions between disease and host genetics.

Results: We identified 12 mbQTLs, including variants in the IBD-associated genes IL17REL, MYRF, SEC16A and WDR78. For example, the decrease of the pathway acetyl-coenzyme A biosynthesis, which is involved in short chain fatty acids production, was associated with variants in the gene MYRF (false discovery rate <0.05). Changes in functional pathways involved in the metabolic potential were also observed in participants carrying rare PTVs or CNVs in CYP2D6, GPR151 and CD160 genes. These genes are known for their function in the immune system. Moreover, interaction analyses confirmed previously known IBD disease-specific mbQTLs in TNFSF15.

Conclusion: This study highlights that both common and rare genetic variants affecting the immune system are key factors in shaping the gut microbiota in the context of IBD and pinpoints towards potential mechanisms for disease treatment.

Keywords: genetics; inflammatory bowel disease; intestinal microbiology.

Figure 1

Schematic overview of the study. (DATA part) We performed whole exome sequencing of the host genome and whole genome shotgun sequencing of faecal samples of 525 individuals (IBD) and 939 controls (LifeLines-DEEP). Nine covariates (age, sex, body mass index (BMI), smoking status, medication use (antibiotics, proton pump inhibitors (PPIs) or laxatives), disease location (in the IBD cohort) and sequencing read depth) were corrected for relative abundances of 242 taxa and 301 pathways. (ANALYSES WORKFLOW part) A four-step analysis was performed: step 1 includes a meta-analysis (p<6.83 $\times$ 10⁻⁷, corresponding to FDR<0.05) in which 73 164 exome-wide common variants with minor allele frequency (MAF) >5% were used for association analyses for microbial traits. Step 2 includes a meta-analysis (p<1.5 $\times$ 10⁻⁵, corresponding to FDR<0.05) using a targeted approach that only tested for 3010 variants located in IBD-associated genes known from IBD genome-wide association studies and PTVs with MAF >5%. Step 3 includes a meta-analysis (p<5 $\times$ 10⁻⁵, corresponding to FDR<0.05) using a gene-based burden test for 980 genes with rare PTVs (MAF <5%); a meta-analysis (p<1.87 $\times$ 10⁻⁴, corresponding to FDR<0.05) using a gene-based test for 267 genes with rare copy number variants (site frequency <1%). Step 4 includes joint analysis combining the two cohorts for disease and genetics interaction analyses. Step 4 focused only on single-cohort-significant microbial quantitative trait loci (mbQTLs) from steps 1 and 2 while adding a disease and a genetic interaction term into the model. All analyses were confined to non-zero values of taxa and pathways. All significance thresholds were set up by Bonferroni correction taking all variants/genes used into account.

Figure 2

Whole-exome-wide meta-analysis results from LifeLines-DEEP and IBD cohorts. Seventy three thousand one hundred and sixty-four common variants (minor allele frequency >5%), 242 taxa and 301 pathways (corrected for all covariates) were used in the association analyses. The discovery significance threshold was p<5 $\times$ 10⁻⁵ and the replication significance threshold was p<0.05. Manhattan plot displays −log10 p values for all association tests. Green and blue represent taxonomies and pathways, respectively. Red line indicates the whole-exome-wide association significance threshold: meta p<6.83 $\times$ 10⁻⁷, corresponding to exome-wide FDR<0.05 (n=73 164 variants, Bonferroni correction).

Figure 3

Microbial quantitative trait loci and eQTL analyses of MYRF. (A) Spearman correlation between genotype (TT, TC, CC) of rs2238001 in MYRF and the relative abundance of acetyl-coenzyme A (CoA) biosynthesis (IBD cohort, p=1.43 $\times$ 10⁻³, r=−0.19; LifeLines-DEEP (LLD) cohort, p=1.47 $\times$ 10⁻⁵, r=−0.20; meta p=7.50 $\times$ 10⁻⁸, FDR<0.05), the glyoxylate bypass and tricarboxylic acid cycle (TCA) MetaCyc pathways (IBD cohort, p=0.0149, r=−0.16; LLD cohort, p=1.04 $\times$ 10⁻⁵, r=−0.22; meta p=6.07 $\times$ 10⁻⁷, FDR<0.05). (B) The rs2238001 locus zoomed in on the IBD-associated region, including the IBD-associated genes MYRF, FADS2 and FADS3. P values are derived from meta-analyses between variants and the relative abundance of acetyl-CoA biosynthesis. (C) eQTL analysis between rs2238001 and MYRF gene expression in colon tissue from the GTEx database (n=246 tissues, p=2.46 $\times$ 10⁻⁷). r, Spearman correlation coefficient.

Figure 4

Associations between gene GPR151 and microbial pathways. (A) Meta p values based on burden test between 30 genes with rare protein truncating variants (PTVs) on chromosome 5 and relative abundance of MetaCyc pathway homolactic fermentation (top). Blue dot represents meta p value of gene GPR151. Lower panel shows the variants found along with the coding region in GPR151. Different colours indicate different variant categories. Red indicates two rare stop-gain mutations, rs114285050 and rs140458264. (B) Box plots for associations between the relative abundance of the homolactic fermentation (meta p=4.78 $\times$ 10⁻⁶, FDR<0.05), glucose xylose degradation (meta p=3.05 $\times$ 10⁻⁵, FDR<0.05) microbial pathways and GPR151, respectively. b, effect size. GPR151, without rare PTVs. GPR151*, with rare PTVs.