Gut microbial keystone taxa associate with genetic variation.

In total, we identified 31 distinct genetic variants associated (P < 5.0 × 10⁻⁸) with 39 microbial taxa related to identified keystone species as listed by Banerjee et al. (2018)^28,31, which included the Actinobacteria class, Helicobacter pylori, Bacteroides stercoris, Bacteroides thetaiotaomicron, Ruminococcus bromii, Klebsiella pneumoniae, Proteus mirabilis, Akkermansia muciniphila and the archaeon Methanobrevibacter smithii (Fig. 1c and Supplementary Table 1). Keystone species are defined as members of a microbial community exerting selective modulation and not broad effects on microbiome composition variation. Only one documented keystone species from Banerjee et al., Bacteroides fragilis, was not associated with genetic variation in our study²⁸. Although a lot of computationally identified keystone species remain to be experimentally verified^32,33, this observation suggests that they would generally associate with human genetic variation. This would indicate an intimate association with the human gut niche in line with their reported key ecological roles in microbiome modulation and functioning. Our work highlights human genotypes associating with keystone taxa (Supplementary Table 1), which could further improve our understanding of their ecology.

Combined effect of genetics and diet on LCT-associated taxa.

We compared the abundances of four bacterial taxa strongly associated with the LCT locus (Bifidobacterium genus, Negativibacillus genus, UBA3855 sp900316885 and CAG-81 sp000435795) in individuals with different rs4988235 genotypes and dairy diets (Fig. 2a). The abundance of Bifidobacterium in individuals producing lactase through adulthood (rs4988235:TT) was unaffected by dairy intake. However, lactose-intolerant individuals (rs4988235:CC) self-reporting a regular dairy diet had a significant increase in Bifidobacterium abundance (P = 1.75 × 10⁻¹³; Wilcoxon rank test). An intermediate genotype (rs4988235:CT) was linked to an intermediate increase in Bifidobacterium abundance (Fig. 2a). This trend did not seem to be affected by age³⁴ (Extended Data Fig. 4). Additionally, we observed a moderate negative correlation between Bifidobacterium abundances and age in rs4988235:CC individuals reporting a regular dairy diet (Spearman’s ⍴ = −0.17, P = 1.9 × 10⁻⁶) and in rs4988235:CC individuals reporting a low-lactose or lactose-free diet (⍴ = −0.19, P = 0.002). Furthermore, the Spearman correlation between the Bifidobacterium residual abundance and dairy diet was still significant (⍴ = −0.22, P = 2 × 10⁻¹²) in rs4988235:CC individuals. This indicated that the associations with age were consistent in individuals with and without regular dairy intake, and did not confound the association between Bifidobacterium and dairy diet.

Open in a separate window

Fig. 2 |

Interaction of human genotype, dairy diet and gut bacterial variation with the LCT locus.

a, The four panels present variation in microbial relative abundances (not CLR-transformed) for the four taxa associated at study-wide significance level with the LCT locus at P < 3.8 × 10⁻¹¹: Bifidobacterium, Negativibacillus, UBA3855 sp900316885 and CAG-81 sp000435795. Abundances are compared across stratified groups of individuals from the FR02 cohort according to LCT-MCM6:rs4988235 genotype and self-reported dietary lactose intake (red, regular dairy diet; blue, lactose-free diet). Sample sizes for groups of individuals self-reporting a regular dairy diet: rs4988235:TT (n = 1,786), CT (n = 2,413), CC (n = 736); self-reporting a nonregular dairy diet or lactose-free diet: TT (n = 150), CT (n = 198), CC (n = 245). All statistical comparisons denote the P values of Wilcoxon rank test on the distributions of untransformed relative abundances. Only significantly different comparisons (P < 0.05) are indicated. For all box plots, the central line, box and whiskers represent the median, interquartile range (IQR) and 1.5 times the IQR, respectively. Violin plots represent the distribution density of the data points.

b, Host genetics and gut microbes interact in the context of dairy intake and lactose intolerance.

An inverse pattern was observed for the abundance distributions of Negativibacillus and uncultured CAG-81 sp000435795, for which abundances decreased in lactose-intolerant individuals reporting dairy intake, as compared with rs4988235:TT individuals consuming dairy products (Fig. 2a). Levels of UBA3855 sp900316885 were unaffected by a dairy diet in lactose-intolerant individuals but were surprisingly lower in rs4988235:TT individuals who reported dairy intake (P = 8.23 × 10⁻⁵). These opposite and contrasting effects of dairy on associated bacterial abundances in lactose-intolerant individuals could reflect competition for lactose in the gut. CAG-81 abundances were the most negatively correlated with those of the other LCT-associated taxa (Extended Data Fig. 5), which suggests that this competition could be strong and prevalent enough to drive coassociation at the LCT locus, possibly mediated by lactose intake (Fig. 2b).

Functional profiling of carbohydrate-active enzymes (CAZymes) in 11 Bifidobacterium species.

Of all 11 Bifidobacterium species prevalent enough in our study population to be included in the GWAS, only Bifidobacterium dentium was not associated with the LCT locus (P = 1.70 × 10⁻²), nor was it coabundant with any other Bifidobacterium species (Extended Data Fig. 6a). B. dentium has previously been suggested to have different metabolic abilities³⁵. A clustering of CAZyme profiles from reference genomes of all 11 Bifidobacterium species revealed that B. dentium clustered apart from the ten other species, which grouped consistently with their coabundance patterns (Extended Data Fig. 6b). B. dentium harbored more genes encoding CAZyme families with preferred fiber/plant-related substrates (GH94, GH26, GH53) than other Bifidobacterium species, which seemed to harbor more milk oligosaccharide-targeting CAZyme families (GH129, GH112) than B. dentium (Extended Data Fig. 6b), which could relate to the observed association differences. This suggests that bacterial metabolic abilities can be strong drivers of coabundance, and of association with human genetic variation.

Impacts of genotype and fiber intake on ABO-associated taxa.

A variety of bacteria metabolize blood antigens, with potential applications in synthetic universal donor blood production^36,37. Gut bacteria are particularly exposed to A- and B-antigens in the gut mucosa of secretor individuals³⁸. Our associations of F. lactaris (P = 1.10 × 10⁻¹²) and Collinsella (P = 2.59 × 10⁻⁸) with ABO suggest a possible metabolic link with blood antigens. A comparison of CAZyme profiles across a set of reference genomes revealed three CAZymes with blood-related activities in F. lactaris (GH110 (ref. ³⁹), GH136 (ref. ⁴⁰), CBM32 (ref. ⁴¹)), but none in any of nine Collinsella species (Fig. 3). More mucus-targeting and fewer fiber-degrading enzymes were found in F. lactaris than Collinsella, suggesting distinct functions in the gut.

Open in a separate window

Fig. 3 |

Functional profiling of reference genomes from two bacterial taxa associated with the ABO locus.

CAZyme distribution patterns in F. lactaris and Collinsella reference genomes (from the GTDB release 89 index used to classify metagenomes in this study). The heatmap indicates species abundance in corresponding CAZyme families, corresponding to the total count of detected families for each species divided by the number of reference genomes examined for the same species. Values < 1 (white to light blue) indicate that less than one copy per genome of the corresponding CAZyme family was detected for each species; values > 1 (light blue to dark blue) indicate that more than one copy per genome was detected. Preferred substrate groups are based on literature search and descriptions on CAZypedia.org.

As previously reported⁴, neither ABO blood types nor secretor status had an impact on alpha- and beta-diversity (Extended Data Fig. 7a). However, we observed that the effects of ABO genotypes on F. lactaris levels, underlying the association, were largely driven by secretor status, with increased abundances in secretor individuals from genotype groups rs545971:CT and rs545971:TT, A and AB blood type groups, but not in rs545971:CC genotype, or B and O blood type individuals (Fig. 4a). Levels in nonsecretors did not vary across ABO genotypes or blood types. Despite a slight increase in blood type A secretors, Collinsella only remained minimally affected by secretor status or blood group (Extended Data Fig. 7b). Taken together, this suggests that the secretion of soluble A- and B-antigens strongly affects F. lactaris in the gut, possibly through reduced opportunity to use them as substrate. Levels of both F. lactaris and Collinsella were significantly higher when individuals were predicted to secrete A-, B- and AB-antigens in their gut mucosa (Extended Data Fig. 7c).

Open in a separate window

Fig. 4 |

Effects of host genetics and dietary fiber intake on gut abundance variation of two bacterial taxa associated with the ABO locus.

a, ABO-associated F. lactaris relative abundances (not CLR-transformed) are compared across stratified groups of individuals from the FR02 cohort according to (left panel) ABO:rs4988235 genotype and predicted secretor status (blue, secretor status conferred by FUT2 rs601338:AG/AA genotype; red, nonsecretor status conferred by FUT2 rs601338:GG genotype), and (right panel) according to predicted A, AB, B and O blood types, and predicted secretor status. Sample sizes for compared groups: secretor status with rs545971:C/C (n = 1,538), C/T (n = 2,493), T/T (n = 1,050) and blood group A (n = 2,178), AB (n = 460), B (n = 900), O (n = 1,543); nonsecretor status with rs545971:C/C (n = 266), C/T (n = 437), T/T (n = 175) and blood group A (n = 383), AB (n = 80), B (n = 148), O (n = 267). b, ABO-associated F. lactaris and Collinsella sp. relative abundances, as well as compounded abundances from 13 mucin-degrading species from Tailford et al. (2015)⁴³, are compared across stratified groups of individuals from the FR02 cohort according to the predicted A/B/AB-antigen secretion status and dietary fiber intake. Secretion status was defined to segregate individuals. A/B/AB-antigen secretors were defined as secretor individuals from blood types A, AB and B. Non-A/B/AB-antigen secretors were defined as nonsecretor individuals and O-antigen secretors. Fiber intake was compared in individual groups from the top and bottom quartiles of total fiber score (Methods). Sample sizes for compared groups of individuals: A/B/AB-antigen secretors (n = 1,393) following a low-fiber diet (n = 723) or a fiber-rich diet (n = 670), or non-A/B/AB-antigen secretors (n = 952) following a low-fiber diet (n = 490) or a fiber-rich diet (n = 462). All statistical comparisons denote the P values of Wilcoxon rank test on the distributions of untransformed relative abundances. For all box plots (b and c), the central line, box and whiskers represent the median, IQR and 1.5 times the IQR, respectively. Violin plots represent the distribution density of the data points. c, Host genetics and gut microbes interact in the context of fiber intake, secretor status and blood types.

A high-fiber diet is thought to induce a metabolic transition from mucus-degrading to fiber-degrading activities in the colon, as carbohydrates from fiber are more easily metabolized⁴². The increase in F. lactaris abundances in A/B/AB-secretors (defined as secreting A-, B- and AB-antigens) compared with non-A/B/AB-secretors remained strongly significant irrespective of fiber intake (P = 1.15 × 10⁻⁹ in the low-fiber diet group, and P = 4.4 × 10⁻³ in the high-fiber diet group), suggesting that F. lactaris has a strong affinity for secreted A/B/AB-antigens, does not efficiently degrade dietary fiber or will not easily switch to it as an energy source (Fig. 4b). F. lactaris levels were increased in non-A/B/AB-secretors with a high-fiber diet compared with a low-fiber diet, implying a switch to fiber degradation or interaction with fiber-degrading bacteria (Fig. 4b). Collinsella variation in both A/B/AB-secretors and non-A/B/AB-secretors with high- and low-fiber diets was similar to the compounded abundances of 13 major mucin-degrading species in the human gut⁴³, suggesting a similar ecological response in stark contrast with F. lactaris (Fig. 4b,​,cc).

MED13L-associated E. faecalis as a putative link with colorectal cancer (CRC).

The allele frequency of the MED13L rs143507801 variant (A > G), associated with levels of E. faecalis (P = 7.26 × 10⁻¹¹), was low (minor allele frequency = 0.0111), consistent with reported allele frequencies in the gnomAD database⁴⁴. In our study population, 131 individuals carried rs143507801:G allele, 130 being heterozygous (GA) and only one being homozygous (GG). We observed that E. faecalis levels were increased in heterozygous rs143507801:GA individuals (Fig. 5). E. faecalis is a gut commensal, but also an opportunist pathogen believed to play a role in CRC development, possibly through direct damaging of colorectal cells^45–47. MED13L and MED13 encode for Mediator transcriptional coactivator complex modules associating with RNA polymerase II (ref.⁴⁸), and as such specifically interact with cyclin-dependent kinase 8 (CDK8) modules described for their oncogenic activation of transcription during colon tumorigenesis⁴⁹. Consequently, we observed slightly higher levels of E. faecalis (P = 0.014) in 14 individuals enrolled in FR02 with a history of CRC at the time of sampling (Fig. 5). Groups of individuals segregated by allelic variant and CRC status could not be compared robustly due to small sample size. Taken together, these results suggest a possible link between E. faecalis and CRC through the MED13 activation of CDK8 in colorectal tumors, which will need to be investigated further.

Open in a separate window

Fig. 5 |

Effect of host genetics and prevalent CRC on gut levels of E. faecalis associated with MED13L variation across participants of the FR02 cohort.

Abundances are compared across individuals grouped according to (left panel) MED13L:rs143507801 genotype and (right panel) CRC prevalence according to the Finnish Cancer Registry. The comparison between E. faecalis variation and MED13L:rs143507801 reflects the GWAS results (Supplementary Table 1). The comparison of E. faecalis abundances in individuals with or without a history of CRC at the time of sampling was performed using a Wilcoxon rank test. Sample sizes for compared groups of individuals: rs143507801:A/A (n = 5,825), G/A (n = 130) (note: only 1 of 5,959 individuals in our cohort was G/G); with CRC (n = 14), without a history of CRC at baseline (n = 5,941). For all box plots, the central line, box and whiskers represent the median, IQR and 1.5 times the IQR, respectively. Violin plots represent the distribution density of the data points.

Cohort phenotype metadata and specific dietary information.

The phenotype data in this study comprised demographic characteristics, life habits, disease history, clinical measurements, laboratory test results and follow-up electronic health records. More specifically, baseline dietary factors were collected. Details of the method have been described previously⁹³. To broadly assess diet information within the cohort participants, a binary variable was used to indicate whether individuals were self-reporting to follow various possible dietary restrictions. Dietary consumption of specific food product categories was also reported. Habitual diet was assessed using a food propensity questionnaire which contained 42 food items or groups and had choices ranging 1–6 for consumption frequency, ranging from ‘Less than once a month’ to ‘Once a day or more often’. The consumption frequencies were converted to frequencies per month, ranging from 0.5 times per month to 30, 45 or 60 times per month. Food items that are rarely eaten more than once a day were given the value of 30 times per month. Food items that are often eaten multiple times a day such as fresh vegetables, breads, and so on were given a value of 60 times per month. Food items that fall in between these two groups were given 45 points.

Self-reporting of lactose-free diet and dietary fiber consumption.

Allelic distribution at the LCT-MCM6:rs4988235 variant responsible for lactase persistence in Europeans was as follows in our study population: 1,936 (35%) individuals had the T/T allele conferring a lactase persistence phenotype through adulthood, allowing them to digest lactose, while 981 (18%) individuals had the C/C allele conferring lactose intolerance. Most individuals (n = 2,611, 47%) had the intermediate allele, C/T, making them likely to be able to digest lactose. Most individuals reported a regular dairy intake in their diet (n = 5,002, 89%), while 706 (12.5%) individuals reported a regular lactose-free diet.

A total fiber consumption score was calculated from the questionnaires, reflecting the overall consumption of a combination of various fiber sources such as high-fiber bread, vegetables (vegetarian dishes, fresh vegetables, and boiled vegetables and legumes) and fruits, berries and natural juices. The resulting total fiber index values ranged from 9 (low dietary fiber intake) to 48 (high dietary fiber intake), with a median of 33. Comparisons of the effects of low- versus high-fiber diets were made between the 1st (n = 1,213) and 4th (n = 1,132) quartiles of the total fiber index.

Genotyping, imputation and quality control.

The genotyping was performed on Illumina genome-wide SNP arrays (the HumanCoreExome BeadChip, the Human610-Quad BeadChip and the HumanOmniExpress) and has been described previously⁹⁶. Stringent criteria were applied to remove samples and variants of low quality. Samples with call rate < 95%, sex discrepancies, excess heterozygosity and non-European ancestry were excluded. Variants with call rate < 98%, deviation from Hardy–Weinberg equilibrium (P < 1 × 10⁻⁶) and minor allele count < 3 were filtered. Data were prephased by using Eagle2 v.2.3 (ref. ⁹⁷). Imputation was performed using IMPUTE2 v.2.3.0 (ref. ⁹⁸) with two Finnish-population-specific reference panels: 2,690 high-coverage whole-genome sequencing and 5,092 whole-exome sequencing samples. To evaluate the imputation quality, we compared the sample allele frequencies with reference populations and examined imputation quality (INFO scores) distributions. Imputed SNPs with INFO > 0.7 were kept for analysis. Postimputation quality control was carried out by using plink v.2.0 (ref. ⁹⁹). Samples with >10% missing rate were removed. Individuals with extreme height or BMI values were further excluded (31 individuals with height < 1.47 m, 5 with BMI > 50 were removed). Both genotyped and imputed SNPs were kept for analysis if they met the following criteria: call rate > 90%, no significant deviation from Hardy–Weinberg equilibrium (P > 1.0 × 10⁻⁶) and minor allele frequency > 1%. SNP filtering was based on all individuals for which genotype information was available (n = 7,280), not on the 5,959 individuals selected subsequently for GWAS after quality control. The postquality control dataset comprised 7,967,866 SNPs.

Metagenomic sequencing from stool samples.

Stool samples were collected by participants and mailed overnight to the Finnish Institute for Health and Welfare for storing at −20 °C; the samples were sequenced at the University of California San Diego in 2017. No special arrangements were made regarding the temperature of the samples when they were shipped from the field clinics to the laboratory in THL but, as the survey was done in Finland during the winter months (January to March 2002), the average temperatures were well below 0 °C. Special care was anyway additionally taken to ensure that samples did not remain sitting in a post office over the weekend. The gut microbiome was characterized by shallow shotgun metagenomics sequencing with Illumina HiSeq 4000 Systems. We successfully performed stool shotgun sequencing in n = 7,231 individuals. The detailed procedures for DNA extraction, library preparation and sequence processing have been previously described⁹⁵. Adapter and host sequences were removed. To preserve the quality of data while retaining most of the disease cases, samples with a total number of sequenced reads lower than 400,000 were removed.

Taxonomic profiling, quality filtering and data transformation.

Taxonomic profiling of FR02 metagenomes was performed as follows: briefly, raw shotgun metagenomic sequencing reads were mapped using the k-mer-based metagenomic classification tool Centrifuge v.1.0.4 (ref. ¹⁰⁰) to an index database custom-built to encompass reference genomes that followed the taxonomic nomenclature introduced and updated in the GTDB release 89 (refs. ^87–89). This implies that unless specified otherwise, all taxonomic names in our study refer to their nomenclature in GTDB, which can be related to the original National Center for Biotechnology Information (NCBI) nomenclature using the GTDB database server: https://gtdb.ecogenomic.org/taxon_history/. The same profiling approach has also been used and described in recent studies from our consortium^94,95,101. Our study present results involving F. lactaris, which is named differently in NCBI and subsequent GTDB releases. A particular note on the evolution of this nomenclature can be found in the Supplementary Note.

Gut microbial composition was represented as the relative abundance of taxa. For each metagenome at phylum, class, order, family, genus and species levels, the relative abundance of a taxon was computed as the proportion of reads assigned to the clade rooted at this taxon among total classified reads. The relative abundance of a taxon with no reads assigned in a metagenome was considered as zero in the corresponding profile. For the purpose of this association study and because of reduced accuracy and power when considering rare taxa, we focused on common and relatively abundant microbial taxa, defined as prevalent in >25% of studied individuals, and defined with at least ten mapped reads per individual. For the purpose of association, and as previous studies have reported that only some microbial taxa are inheritable¹⁰², we also removed taxa with zero SNP heritability. This filtering resulted in a microbial dataset composed of a total of 2,801 taxa, including 59 phyla, 95 classes, 187 orders, 415 families, 922 genera and 1,123 species.

Taxonomic profiles derived from sequencing data are by nature compositional because of an arbitrary total imposed by the instrument¹⁰³. The compositional data of microbial taxa are not independent and can lead to inappropriate use of linear regression. To overcome this artificial bias, all relative abundance values were transformed by center-log-ratio (CLR)¹⁰⁴. More information about data transformation can be found in the Supplementary Note.

When visually comparing relative abundances in groups of individuals throughout the manuscript, we used untransformed relative abundances, for better interpretability. Alpha- (Shannon index) and beta- (Bray–Curtis distance) diversity were calculated at genus level used functions in the R package vegan v.2.5–6. We did not find a correlation between sequencing depth and Shannon diversity index (Spearman’s ⍴ = −0.001598, P = 0.90) in n = 5,959 samples (Extended Data Fig. 8). Additionally, to define CLR-transformed abundances of higher taxonomic levels than species, we summed the raw abundances of all taxa (for example, species) belonging to a specific higher taxonomic taxon (for example, genus), and then applied a CLR transformation. Additionally, we observed that Eastern and Western Finnish populations did not have different microbiome diversity, despite having overall slightly different lifestyles and mortality rates. To further investigate this, we visualized potential geographical effects using a Principal Coordinates Analysis (PCoA) plot on beta-diversity (Bray–Curtis dissimilarity) from metagenomic profiles of samples used in the GWAS from our study (n = 5,959; Extended Data Fig. 9).

Genome-wide association analysis.

The protocol followed in this study was described elsewhere¹⁰⁵. Briefly, a linear mixed model (LMM) implemented in BOLT-LMM v.2.3.2 (ref. ¹⁰⁶) was used to search for genome-wide associations accounting for the individual similarity. Since BOLT-LMM only accepts <1 million SNPs in modeling the genetic relationship matrix, SNPs were pruned at the threshold of r² < 0.1 (plink2 (ref. ⁹⁹), command–indep-pairwise 1000 80 0.1), resulting in 106,201 independent SNPs. This list of independent SNPs was used to estimate heritability using BOLT-LMM. Additionally, BOLT-LMM automatically performs leave-one-chromosome-out analysis to avoid proximal contamination. Although the LMM accounts for the cryptic relatedness in individuals, there are still large population structures that cannot be addressed. Thus, the top ten genetic principal components (calculated by FlashPCA v.2.0 (ref. ¹⁰⁷) based on the pruned SNPs mentioned above) were included as covariates, in addition to age, sex and genotyping batch. We did not adjust for microbiome sequencing batch, as we observed that it had no effect on microbiome composition variation (Extended Data Fig. 9). As no genetic variant was reported to have a large effect size on gut microbiota, statistical estimates were based on infinitesimal model which assumes a small nonzero effect for a large number of genetic variants. To identify independent associations, GCTA-COJO v.1.91.3 (ref. ¹⁰⁸) was used to conduct approximate conditional and joint analysis using individual genetic data. Window size was set to 10 megabases (Mb), assuming SNPs on different chromosomes or more than 10-Mb distance apart are uncorrelated. The resulting effect size (beta coefficient) indicated the number of standard deviation changes of a taxon’s CLR-transformed abundance corresponding to one effective allele increase of SNPs. Additionally, for all but two reported SNPs (rs146740485 and rs2797225), the effect allele was the reference allele in the GWAS cohort.

As microbes interact nonindependently with each other in the gut, as part of larger ecological and functional communities, matSpDlite v.1.0 (refs. ^109,110) was used to estimate the number of independent tests based on eigenvalue variance—the larger the eigenvalue variance, the smaller the number of effective tests. The number of independent tests was 1,328 for 2,801 tested taxa. We used this information to calculate a Bonferroni-adjusted study-wide significance level for significant associations, which was set to 5 × 10⁻⁸/1,328 = 3.8 × 10⁻¹¹. A genome-wide significance threshold was set as 5 × 10⁻⁸. The identified SNPs were annotated using ANNOVAR v.2018Apr16 (ref. ¹¹¹) and grouped into genetic loci using 200-kilobase windows flanking the top SNPs.

We also examined whether antibiotic prescription before baseline sampling could be an important confounder of results. We obtained individual information on the prescription of any antibiotic up to 1 month before baseline fecal sampling, corresponding to 250 individuals of 5,959 (4.2%). We examined whether individual microbial profiles (via beta-diversity estimates using Bray–Curtis dissimilarity) were broadly affected by recent antibiotic prescription and observed a slight effect along PCoAs with significant variance explained (Extended Data Fig. 9c). After repeating the GWAS for all microbial taxa for which we initially had found at least one significantly associated locus, this time adjusting for previous antibiotic prescription status (‘yes’ versus ‘no’) (Supplementary Table 9), we found that recent antibiotic prescriptions had very minor effects on the GWAS association results. Adjusting for antibiotic prescription did not change any study-wide significant associations and only 32 of 567 genome-wide associations moved slightly above P = 5 × 10⁻⁸ (the largest P value was 3.2 × 10⁻⁷), which is likely by chance given inclusion of any additional covariate (Supplementary Table 9). In addition, the beta estimates with and without the adjustment of antibiotics usage were highly consistent (Pearson r = 0.9999487).

One important association in our study involved F. lactaris abundance and variants in the ABO locus. We observed that the distribution of F. lactaris abundance in our GWAS cohort (n = 5,959) was slightly bimodal (Extended Data Fig. 10). To investigate whether a logistic model gives the same result for this taxon, we arbitrarily coded F. lactaris abundance as ‘1’ if the relative abundance was higher than 5 × 10⁻⁴ (n = 2,866), and ‘0’ if smaller (n = 3,093). Akaike information criterion (AIC) value was smaller for logistic than for linear models (AIC = 8,196 versus AIC = 12,463, respectively), and the strongest association was also observed in the same top SNP (rs545971, P = 5.5 × 10⁻¹⁸) as when using linear regression (rs545971, P = 1.1 × 10⁻¹²).

Replication of previously reported associations.

To evaluate the reproducibility of our results with previously reported associations, we collected GWAS summary results from eight studies published in peer-reviewed journals at the time of this work^{3,6–10,102,112}. These studies reported associations between 548 SNPs and microbial features. ANNOVAR was used to annotate the reported SNPs to the hg38 human reference genome¹¹¹ and we used plink2 (ref. ⁹⁹) to identify a further 15,427 SNPs in high LD (r² > 0.8, within 5 Mb) with any of these 548 SNPs. To assess replication, we first examined whether previously reported associations could be matched in our results to identical or linked SNPs, with an association below the Bonferroni-corrected suggestive significance threshold, which was set to 0.05/548 = 9.124 × 10⁻⁵. More details about the replication methods and the use of the GTDB taxonomic system can be found in the Supplementary Note.

Prediction of ABO blood groups and secretor status.

SNP-based typing of ABO histo-blood group was performed. A combination of four SNPs¹¹³ was used for the prediction, and a 98% concordance with phenotypically typed ABO histo-blood group has been reported for this method⁴. For blood group allele A, the two different types, A1 and A2, were predicted by rs507666 and rs8176704, respectively. Blood group allele B was inferred from rs8176746 and blood group allele O was predicted by rs687289. As the combinations of these SNPs are exclusive, no haplotype information was needed. To validate the accuracy of prediction, we compared it with the prediction using a different combination of SNPs⁶⁸. The two predictions were highly consistent, with over 99.9% concordance. In addition, the distribution of ABO groups was consistent with the population distribution found in public databases. Secretor status was predicted by the genotype of FUT2 variant rs601338, where AA or AG genotypes are secretors and GG genotypes are nonsecretors. A 100% concordance between the variation in rs601338 and secretor status was reported in a study on Finnish individuals¹¹⁴.

Bidirectional two-sample MR analysis.

Causal relationships between diseases and gut microbiota were investigated at genus and species levels only to maximize interpretability. In total, 213 species and 148 genera associated with at least one variant at genome-wide significant level (P < 1 × 10⁻⁸) were included. GWAS summary results were collected for 46 diseases from MR-Base¹¹⁵ (Supplementary Table 4). These included 12 autoimmune or inflammatory diseases, 9 cardiometabolic diseases, 13 psychiatric or neurological diseases, 4 bone diseases and 8 cancers. For diseases with more than one GWAS record, the record with the largest sample size was kept.

Bidirectional causal inference was performed to infer causal effects of microbial abundance variation (exposure) on disease risk (outcome), and of disease (exposure) on microbial abundance levels (outcome). To select the SNP instruments for microbial exposures in our study (Supplementary Table 7), we followed recommendations from a previous study showing that associated SNPs below a significance threshold of P < 1 × 10⁻⁵ had the largest explained variance on microbial features¹¹⁶. For each taxon, GCTA-COJO was used to perform a conditional analysis to select independently associated SNPs at P < 1 × 10⁻⁵. F statistics were calculated to estimate the strength of instruments for each bacterial exposure, and were found to be >10 for all exposures (Supplementary Table 5). SNP instruments for disease exposures were selected at genome-wide significance threshold (P < 5 × 10⁻⁸). Subsequently LD-clumping with a strict threshold (r² < 0.001 in the 1000 Genomes European data within 10 Mb windows) was conducted to select independent instruments with the lowest P values for taxa and diseases, respectively.

Details about the precise methods used for MR inference can be found in the Supplementary Note.

Cox proportional hazards regression.

Cox proportional hazards regression was conducted to test the association between baseline abundance of gut microbes and incident major depression (16 yr of follow-up, n = 181 incident events). Microbial abundances were CLR-transformed and standardized to zero-mean and unit-variance. The Cox models were stratified by sex and adjusted for age and log-transformed BMI, with time-on-study as the time scale. Participants with prevalent major depression at baseline were excluded. R function coxph() in the R package survival v.3.1–8 was used for this analysis.

Profiling of CAZymes in bacterial genomes.

The standalone run_dbCAN2 v.2.0.11 tool¹¹⁷ (https://github.com/linnabrown/run_dbcan) was used to scan for the presence of CAZyme genes in public assembled bacterial genomes taken from the GTDB release 89 reference. We used a CAZyme reference database taken from the CAZy database¹¹⁸ (31 July 2019 update). In total, we scanned 327 Bifidobacterium sp., 2 F. lactaris and 15 Collinsella sp. reference genomes included in GTDB release 89. Three methods were compared as part of the run_dbCAN2 procedure (HMMER, DIAMOND and Hotpep). We considered it a positive detection result when all three methods agreed on a CAZyme family identification. Identification of preferred reported substrates for the various CAZyme families was done manually from key publications^42,119, from literature searches and from the CAZypedia website¹²⁰. Certain CAZyme families have a broad range of substrates, many of which are still unknown, which results in our reported preferred substrates to be as accurate as possible, but nonexhaustive.

Carbon impact and offsetting.

We used GreenAlgorithms v.1.0 (ref. ¹²¹) to estimate that the main computational work in this study had a carbon impact of at least 2,660 kg of CO₂ emissions (CO₂e), corresponding to 233 tree-years. As a commitment to the reduction of carbon emissions associated with computation in research, we consequently funded planting of 30 trees through a local Australian charity, which across their lifetime will sequester a combined estimated 8,040 kg of CO₂e, or three times the amount of CO₂e generated by this study.

The study protocol of FINRISK 2002 was approved by the Coordinating Ethical Committee of the Helsinki and Uusimaa Hospital District (Ref. 558/E3/2001). All participants signed an informed consent. The study was conducted according to the World Medical Association Declaration of Helsinki on ethical principles. All necessary patient/participant consent has been obtained and the appropriate institutional forms have been archived. We thank all participants of the FINRISK 2002 survey for their contributions to this work. The FINRISK surveys are mainly funded by budgetary funds from the Finnish Institute for Health and Welfare with additional funding from several domestic foundations. Y.Q. was partially supported by The Albert Shimmins Fund (Faculty of Science Postgraduate Writing-Up Award 2020). M.I. was supported by the Munz Chair of Cardiovascular Prediction and Prevention and the NIHR Cambridge Biomedical Research Centre (BRC-1215-20014). V.S. was supported by the Finnish Foundation for Cardiovascular Research. L.L. was supported by the Academy of Finland (decision 295741) and EU/H2020 (FindingPheno; grant 952914). T.N. was supported by the Emil Aaltonen Foundation, the Finnish Medical Foundation, the Paavo Nurmi Foundation and the Academy of Finland (grant no. 321351). A.S.H. was supported by the Academy of Finland, grant no. 321356. R.L. receives funding support from NIEHS (grant no. 5P42ES010337), NCATS (grant no. 5UL1TR001442), NIDDK (grant nos. U01DK061734, R01DK106419, P30DK120515, R01DK121378, R01DK124318) and DOD PRCRP (grant no. W81XWH-18-2-0026). S.C.R. is funded by a BHF Programme Grant (RG/18/13/33946). This study was supported by the Victorian Government’s Operational Infrastructure Support (OIS) program, and by core funding from the British Heart Foundation (grant no. RG/13/13/30194; RG/18/13/33946) and the NIHR Cambridge Biomedical Research Centre (BRC-1215-20014). National Institute for Health Research (Cambridge Biomedical Research Centre at the Cambridge University Hospitals NHS Foundation Trust) (The views expressed are those of the authors and not necessarily those of the NIHR or the Department of Health and Social Care). This work was supported by Health Data Research UK, which is funded by the UK Medical Research Council, the Engineering and Physical Sciences Research Council, the Economic and Social Research Council, the Department of Health and Social Care (England), the Chief Scientist Office of the Scottish Government Health and Social Care Directorates, the Health and Social Care Research and Development Division (Welsh Government), the Public Health Agency (Northern Ireland), the British Heart Foundation and Wellcome. We thank Dr Annalisa Buniello (EMBL-EBI, Cambridge, UK) for valuable help with GWAS Catalog submissions.