Introduction discovered to be prevalent in >25% of the cohort (Methods). Using a genome-wide 112 significance threshold (p<5.0×10 -8 ), a total of 478 distinct GTDB taxa, which represented 17% 113 of all tested taxa and included 11 phyla, 19 classes, 24 orders, 63 families, 148 genera and 213 114 species, were found to be associated with at least one genetic variant (Figure 1, Table S1 ). 115 Conditional analysis found 583 independent SNP-taxon associations at genome-wide 116 significance (Table S1) . Heritability across the 2,801 taxa ranged between h 2 =0.001 to 0.214, 117
with the highest values observed for taxa belonging to the Firmicutes and Firmicutes_A GTDB 118 phyla, both of which encompassed half (241/476, 50.4%) of all associated taxa with genetic 119 variation ( Figure S2 ). There were no differences in SNP heritability between groups of 120 associated or non-associated taxa at genome-wide significance (p=0.23). 121 122
Three loci were strongly associated with microbial variation at study-wide significance, as 123
shown on a Manhattan plot showing the lowest resulting p-value for each SNP tested against 124 each of the 2,801 taxa (Figure 1, Table 1 ). There was no evidence of excess false positive rate 125
in the GWAS (median λ GC =1.0051) ( Figure 1B) . After conditional analysis, the strongest 126 association by far (p=5.0×10 -35 ) involved members of class Actinobacteria and rs3940549, a 127 variant in the LCT-MCM6-ZRANB3 locus region which is in high LD (r 2 =0.87) with the well-128
described LCT variant rs4988235 causing lactase persistence in adults of European ancestry 129 ( Figure S3 ). In total, 29 taxa were associated with the LCT-MCM6 region, including 18 below 130 study-wide significance (Figure 1 , Table S1 ). These involved Bifidobacterium-related 131
Actinobacteriota and three taxa from the GTDB Firmicutes_A phylum which included 2 132 uncultured species defined from metagenome-assembled reference genomes (UBA3855 133 sp900316885 and CAG-81 sp000435795) ( Table 1) . The association of these three 134
Firmicutes_A with LCT was still genome-wide significant after adjusting for Bifidobacterium 135 abundances (Table S2) . A variant in ABO (rs545971), expressing the histo-blood 136
group ABO system transferase, was strongly associated (p=1.1×10 -12 ) with levels of 137
Faecalicatena lactaris. There was evidence for a second independent signal at ABO associated 138 with the Collinsella genus (chr9:133271182; p=2.5×10 -8 ) (Table S1, Figure 1 ). Rs187309577 139 and rs143507801 in MED13L, expressing the Mediator complex subunit 13L, were found to be 140 associated with genus Enterococcus (p=1.8×10 -12 ) and the Enterococcus faecalis species 141 (p=7.26×10 -11 ), respectively (Table S1, Figure 1 ).
Human gut microbiome keystone taxa are associated with genetic variation 144 145
In total, we identified 31 distinct genetic variants associated (p<5.0×10 -8 ) with 39 microbial 146 taxa related to identified keystone species as listed by Banerjee et al. (2018) 29 Akkermansia muciniphila 31 , and the archaeon Methanobrevibacter smithii 37,38 ( Figure 1C , 150 Table S1 ). Only one documented keystone species from Banerjee et al. 29 , Bacteroides 151 fragilis 39 , was not associated with genetic variation in our study. This observation suggests that 152 keystone species, although defined as exerting selective modulation and not broad effects on 153 microbiome composition variation, generally associates with human genetic variation, 154
suggesting an intimate association with the human gut niche, in line with their reported key 155 ecological roles in microbiome modulation and functioning. Our work highlights novel human 156 genotypes possibly associated with keystone taxa (Table S1) , which could further improve our 157
understanding of their ecology. 158 159
Combined effect of host genetics and dietary dairy intake on gut levels of LCT-associated 160 bacteria 161 162 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. .
We compared the abundances of 4 bacterial taxa strongly associated with the LCT locus 163 (Bifidobacterium genus, Negativibacillus genus, UBA3855 sp900316885 and CAG-81 164 sp000435795) in individuals with different rs4988235 genotypes and dairy diets (Figure 2A) . 165 The abundance of Bifidobacterium in individuals producing lactase through adulthood 166 (rs4988235:TT) was unaffected by dairy intake. However, lactose-intolerant individuals 167 (rs4988235:CC) self-reporting a regular dairy diet had a significant increase in Bifidobacterium 168 abundance (p=1.75×10 -13 ; Wilcoxon-rank test). An intermediate genotype (rs4988235:CT) was 169 linked to an intermediate increase (Figure 2A) . This trend did not seem to be affected by age 40 170 ( Figure S4 ). 171 172
An inverse pattern was observed for the abundance distributions of Negativibacillus and 173
uncultured CAG-81 sp000435795, for which abundances decreased in lactose intolerant 174
individuals reporting dairy intake, as compared to rs4988235:TT individuals consuming dairy 175 products (p=0.049 and p=0.041, respectively) (Figure 2A) . Levels of UBA3855 sp900316885 176
were unaffected by a dairy diet in lactose-intolerant individuals but were surprisingly lower in 177 rs4988235:TT individuals who reported dairy intake (p=8.23×10 -5 ) (Figure 2A) . These 178 opposite and contrasting effects of dairy intake on associated bacterial abundances in lactose-179 intolerant individuals could reflect competition for lactose in the gut. Genus CAG-81 180 abundances were the most negatively correlated with those of the other LCT-associated taxa 181 ( Figure S5 ), which suggests that this competition could be strong and prevalent enough to 182 drive co-association at the LCT locus, possibly mediated by lactose intake ( Figure 2B ).
Functional profiling of CAZymes in 11 Bifidobacterium species 185 186
Of all 11 Bifidobacterium species prevalent enough in our study population to be included in 187 the GWAS, only B. dentium was not associated with the LCT locus (p=1.70×10 -2 ), nor was it 188
co-abundant with any other Bifidobacterium species ( Figure S6A) . B. dentium has previously 189 been suggested to have different metabolic abilities 41 . A clustering of carbohydrate-active 190 enzymes (CAZyme) profiles from reference genomes of all 11 Bifidobacterium species 191 revealed that B. dentium clustered apart from the 10 other species, which grouped consistently 192 with their co-abundance patterns ( Figure S6B) . B. dentium harboured more genes encoding 193
CAZyme families with preferred fiber/plant-related substrates (GH94, GH26, GH53) than 194
other Bifidobacterium species, which seemed to harbour more milk oligosaccharide-targeting 195
CAZyme families (GH129, GH112) than B. dentium ( Figure S6B ), which could relate to the 196 observed association differences. This suggests that bacterial metabolic abilities can be strong 197 drivers of co-abundance, and of association with human genetic variation. 198 199
Functionally distinct ABO-associated bacteria are impacted differently by genotype and 200 dietary fiber intake 201 202
A variety of bacteria metabolize blood antigens, with potential applications in synthetic 203
universal donor blood production 42,43 . Gut bacteria are particularly exposed to A-and B-204 antigens in the gut mucosa of secretor individuals 44 . Our associations of Faecalicatena lactaris 205 (p=1.10×10 -12 ) and Collinsella (p=2.59×10 -8 ) with ABO suggest a possible metabolic link with 206 blood antigens. A comparison of CAZyme profiles across a set of reference genomes revealed 207
3 CAZymes with blood-related activities in F. lactaris (GH110 45 , GH136 46 , CBM32 47 ), but 208 none in any of 9 Collinsella species ( Figure 3A) . More mucus-targeting and less fiber-209 degrading enzymes were found in F. lactaris than Collinsella (Figure 3A) , suggesting distinct 210 functions in the gut. 211 212 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. .
As previously reported 5 , neither ABO blood types, nor secretor status had an impact on alpha 213
and beta diversity ( Figure S7 ). However, we observed that the effect of ABO genotypes on F. 214 lactaris levels, underlying the association, were largely driven by secretor status, with 215 increased abundances in secretor individuals from genotype groups rs545971:CT (p=3.6×10 -4 ) 216 and rs545971:TT (p=9×10 -4 ), A (p=1.24×10 -5 ) and AB blood type groups (p=1.24×10 -5 ), but 217 not in rs545971:CC genotype (p=0.4339), or B and O blood types individuals ( Figure 3B) . 218 Levels in non-secretors did not vary across ABO genotypes or blood types ( Figure 3B ).
Despite a slight increase in blood type A secretors, Collinsella only remained minimally 220 affected by secretor status or blood group ( Figure S8A) . Taken together, this suggests that the 221 secretion of soluble A and B-antigens strongly affects F. lactaris in the gut, possibly through 222 reduced opportunity to use them as substrate. Both levels of F. lactaris and Collinsella were 223 significantly higher when individuals were predicted to secrete A-, B-and AB-antigens in their 224 gut mucosa (p<2.2×10 -16 and p=1.3×10 -8 , respectively) ( Figure S8B ).
A high fiber diet is thought to induce a metabolic transition from mucus-degrading to fiber-227 degrading activities in the colon, as carbohydrates from fiber are more easily metabolized 48 . 228
The increase in F. lactaris abundances in A/B/AB-secretors (defined as secreting A-, B-and 229 AB-antigens) compared to non-A/B/AB-secretors remained strongly significant irrespective of 230 fiber intake (p=1.15×10 -9 in the low-fiber diet group, and p=4.4×10 -3 in the high-fiber diet 231 group), suggesting that either F. lactaris has a strong affinity for secreted A/B/AB-antigens, 232
does not efficiently degrade dietary fiber, or will not easily switch to it as an energy source 233 ( Figure 3C ). F. lactaris levels were increased in non-A/B/AB-secretors with a high fiber diet, 234
implying a switch to fiber degradation or interaction with fiber-degrading bacteria ( Figure  235 3C). Collinsella variation in both A/B/AB-secretors and non-A/B/AB-secretors with high-and 236 low-fiber diets was similar to the compounded abundances of 13 major mucin-degrading 237 species in the human gut 49 , suggesting a similar ecological response in stark contrast with F. 238 lactaris ( Figure 3C , Figure 3D ).
MED13L association with Enterococcus faecalis as a putative link with CRC development 241 242
The allele frequency of the MED13L rs143507801 variant (A>G), associated with levels of 243
Enterococcus faecalis (p=7.26×10 -11 ), was low (MAF=0.0111), consistent with reported allele 244 frequencies in the gnomAD database 50 . In our study population, 131 individuals carried 245 rs143507801:G allele, 130 being heterozygous (GA) and only one being homozygous (GG). 246 We observed that E. faecalis levels were increased in heterozygous rs143507801:GA 247 individuals (Figure 4 ). E. faecalis is a gut commensal, but also an opportunist pathogen 248 believed to play a role in colorectal cancer (CRC) development, possibly through direct 249 damaging of colorectal cells 51-56 . MED13L and MED13 encode for Mediator transcriptional 250 coactivator complex modules associating with RNA polymerase II 57 , and as such specifically 251 interact with cyclin-dependent kinase 8 (CDK8) modules described for their oncogenic 252 activation of transcription during colon tumorigenesis 58 . Consequently, we observed slightly 253 higher levels of E. faecalis (p=0.014) in 14 individuals enrolled in FR02 who had prevalent 254 CRC at the time of sampling (Figure 4) . Groups of individuals segregated by allelic variant 255
and CRC status could not be compared robustly due to small sample size. Taken together, these 256 results suggest a possible link between E. faecalis and CRC through the MED13 activation of 257 CDK8 in colorectal tumours, which will need to be investigated further. 258 259
Causal inference predictions between microbes and diseases highlight causal effect of 260
Morganella on MDD 261 262 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. .
Interpreting results of causal inference prediction using bacterial information entails to 263 particular caution, due to the possibility of multiple and unaccounted confounding factors 11 , 264 but can be useful to highlight potential focus for future research. Here, we predicted 96 causal 265 effects in both microbe to disease and disease to microbe directions using bidirectional 266
Mendelian Randomization (MR). Of these, 34 were from microbial levels as exposure to 267 disease as outcome, with a large proportion of causal effects in psychiatric and neurological 268 diseases (Table S5) sex and BMI ( Figure 5) . In our GWAS, Morganella variation in the study population 293 associated with a variant (rs192436108; p=6.16×10 -8 ) in the PDE1A locus, which has 294 previously been linked to depression 62,63 and psychiatric disorders 64 . Taken together, these 295 predicted links between Morganella and MDD suggest more efforts should be deployed into 296 exploring the possible roles of this bacterium as part of the brain-gut axis metabolic 297 modulation of health. 298 299
Discussion 300 301
Here, through GWAS and the subsequent investigation of functional and ecological factors 302
contributing to the most robust human-microbe associations, we present a diverse and global 303 picture of human-microbe interactions in a single cohort of ~6,000 European individuals. We 304 find 3 genetic loci to be strongly associated with gut microbial variation. Two of these loci, 305
LCT and ABO, are well-known and very segregated in human populations, possibly explaining 306 why our homogenous European cohort identified them as being associated so strongly. A third 307 more mysterious association with the MED13L locus highlights possible links with cancer 308
while predictive causal inference highlights several diseases as being causally linked to gut 309 microbes. 310 311
Lactase persistence as a recently evolved strong modulator of gut bacterial abundances 312 313 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. . Lactase persistence, or the continued ability to digest lactose into adulthood, is the most 314 strongly selected single-gene trait over the last 10,000 years in multiple human populations 65 , 315 believed to have spread amongst humans with the advent of animal domestication and the 316 culturally transmitted practice of dairying 66 . In our study, as in previous work 4,6,7,11,12 , the 317 association of LCT variants with Actinobacteria, more specifically Bifidobacterium, is by far 318 the most statistically significant, suggesting a profound interaction between Actinobacteria and 319 the human gut, in line with their reported keystone activities 30 . We reported a strong increase 320
of Bifidobacterium levels in genetically lactose intolerant people reporting a regular 321 consumption of dairy products 9 . This increase was not confounded by age in adults, despite 322
Bifidobacterium levels generally decreasing with age in our cohort. While self-reported dietary 323
information is not entirely reliable due to various social reasons 67,68 , our study population was 324 large, and the differences were significant enough to consider this a robust observation. These 325 observations can be explained by the evolutionary adaptation of Bifidobacterium species to 326 specifically use human and bovine milk oligosaccharides as an energy source 69 . In adults 327 unable to produce lactase in their small intestines, consumed lactose is likely to become 328 available for colonic bacteria as an energy source to compete for ( Figure 3A ). Hints of a 329 possible competitive relationship between Bifidobacterium and Negativibacillus, another LCT-330 associated taxon were revealed, which could be mediated by lactose intake and will need to be 331 investigated further in functional studies. 332 333
Two interesting questions stem from our findings. First, the genetic determinants of lactose 334
intolerance are known to vary across ethnicity 70 and cross-population heterogeneity in the 335
LCT-Bifidobacterium association was recently reported 12 . As more non-European-centric 336 genetic studies are conducted worldwide 12,71,72 , examining this combined interaction between 337 dairy diet and Bifidobacterium in different genetic backgrounds could bring new insights. 338
Secondly, despite recent progresses, lactose intolerance is still largely underdiagnosed, and 339 genetic prediction rates from large population studies exceed lactose intolerance prevalence 340 rates obtained using physical tests 70 . In our work, we lacked information on lactose 341 malabsorption symptoms in lactose intolerant individuals reporting a regular dairy diet. These 342 people could experience discomfort symptoms without knowingly implicating their own 343 lactose intake, but another possibility could be that the ability of Bifidobacterium to degrade 344 lactose may alleviate the perceived symptoms of discomfort associated with lactose 345 intolerance, therefore encouraging individuals to unknowingly continue consuming lactose that 346 they would otherwise not be able to digest 73 . This possible probiotic effect would be interesting 347
to investigate in controlled studies. 348 349
Blood antigen secretion can influence levels of specific gut microbial commensals 350 351
The ABO gene expresses a glycosyltransferase in many cell types, which determines the ABO 352 blood group of an individual by modifying the oligosaccharides on cell surface glycoproteins. 353
A comparison of humans and non-human primates has identified ABO (along with the MHC) 354
as harbouring ancient multiallelic polymorphisms that are maintained across species 74,75 . 355
Evolutionary selective pressures at this locus have been proposed to be linked to pathogen 356 infection. Indeed, many infectious diseases such as norovirus infection, bacterial meningitis, 357 malaria, cholera 76 , or even more recently SARS-CoV-2 77,78 are associated with host blood type 358 and secretor status 76 , suggesting that infection could be a driver of a strong balancing selection 359 that has maintained ABO polymorphisms. Furthermore, blood type variation has been 360 intriguingly linked to various chronic diseases 76 , such as heart and vascular diseases, gastric 361 cancers, diabetes, asthma or even dementia 76 . Many of these chronic diseases are also 362 associated with dysbiosis of the gut microbiota, which prompts interesting but largely 363 unexplored parallel between gut commensals, blood types and disease 44 . Our study confirms 364 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. . previous findings 5 that secretor status or blood types do not seem to globally affect gut 365 microbial alpha-or beta-diversity. It also confirms reports from two very recent studies: the 366 first of these studies, a meta-analysis across five German cohorts, using 16S rRNA sequencing 367
to characterize the gut microbiota, is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. . years after initial sampling. This observation supports previous experimental results showing 416 an increase of IgM and IgA-related immune response against Morganella secreted 417 lipopolysaccharide in major depression 88 . This finding potentially highlights the intimate 418
influence of the gut-brain axis on humans. 419 420
Our MR analysis suggested that known genetic risks of autoimmune and inflammatory 421 diseases could also influence gut microbes. One explanation could be that disease susceptibility 422 would affect host immunity and gut barrier integrity, which may favor an increase in some key 423 microbes. However, several studies have shown that manipulating gut microbial composition 424 could be a potential therapy for autoimmune and inflammatory diseases 89 , which would 425
suggest that composition variation in specific gut microbe maybe a requirement for the 426 penetration of a disease phenotype 90 . Further mechanistic studies are needed to untangle host-427 microbe interactions in disease, and further interpret these predictions. 428 429
The case for larger datasets and including uncultured novel species in metagenomic 430 studies 431 432
Our study highlights the benefits of increasing sample size to increase the statistical power for 433 discovery. Although the LCT locus has been reported multiple times to be associated with 434 bacterial taxa, our work is the first to report study-wide significant associations in a single 435 cohort, at the strongest significance ever reported. The association with Bifidobacterium in our 436 study was even stronger than the recent findings that used integrative data from 18,473 437 individuals in 28 different cohorts 12 , emphasizing the importance of standardized methodology 438 and homogeneity in participant ethnicity (especially when studying highly geographically 439 distributed traits such as lactose intolerance traits 91 ). ABO allelic variation is also notoriously 440 affected by geography 92 , which could explain why some meta-analyses in non-homogenous 441 populations could miss it or not. Importantly, metagenomic sequencing with standardized, 442 robust taxonomic definitions 93,94 can provide species-level characterization of microbial 443 profiles in the gut of individuals, which is challenging when using 16S rRNA-based studies. 444
An example from our work is the observation that Bifidobacterium dentium was prevalent but 445 not associated with the LCT locus like all other Bifidobacterium species in the population. 446
Observed difference in carbohydrate-active enzymes that are commonly found in other 447
Bifidobacterium species may explain this difference 41 . Furthermore, GTDB taxonomic 448 standardization results in greater taxon granularity, i.e. smaller, more discrete clades of similar 449 phylogenetic depth than commonly known lineages or species 93,94 . In theory, this would 450
increase overall accuracy 95 , as a weak association with a poorly-defined lineage may be caused 451 by a strong association with a well-defined subset of that lineage, defined as a coherent group 452 using GTDB 94 . Finally, a myriad of microbial taxa that are to date solely defined and 453
represented by uncultured metagenome-assembled genomes (MAGs) in the GTDB database 454
were found to be independently associated with various loci. Along with recent reports that the 455 more gut microbiome diversity is explored, the more novel, unknown species are 456 discovered 96,97 , this suggests that many discoveries are yet to be made in the field of human 457 microbiome studies. 458 459 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. .
Study population 462 463
The FINRISK study population has been extensively described elsewhere 98 A total fiber consumption score was calculated from the questionnaires, reflecting the overall 507 consumption of a combination of various fiber-rich foods such as high-fiber bread, vegetables 508 (vegetable foods, fresh and boiled) and berries (fruits, berries and natural juices). The resulting 509 total fiber index values ranged from 9 (low dietary fiber intake) to 48 (high dietary fiber 510 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this this version posted September 13, 2020. . intake), with a median of 33. Comparisons of the effects of low-vs. high-fiber diets were made 511 between the 1 st (n=1,213) and 4 th (n=1,132) quartiles of the total fiber index. 512 513
Genotyping, imputation and quality control 514 515
The genotyping was performed on Illumina genome-wide SNP arrays (the HumanCoreExome 516
BeadChip, the Human610-Quad BeadChip and the HumanOmniExpress) and has been 517 described previously 102 . Stringent criteria were applied to remove samples and variants of low 518 quality. Samples with call rate <95%, sex discrepancies, excess heterozygosity and non-519
European ancestry were excluded. Variants with call rate <98%, deviation from Hardy-520
Weinberg Equilibrium (p<1×10 -6 ), and minor allele count < 3 were filtered. Data was pre-521 phased by using Eagle2 v2.3 103 . Imputation was performed using IMPUTE2 v2.3.0 104 with two 522
Finnish-population-specific reference panels: 2,690 high-coverage whole-genome sequencing 523
and 5,092 whole-exome sequencing samples. library preparation and sequence processing have been previously described 101 . Adapter and 541
host sequences were removed. To preserve the quality of data while retaining most of the 542 disease cases, samples with a total number of sequenced reads lower than 400,000 were 543 removed. 544 545
Taxonomic profiling, quality filtering and data transformation 546 547
Taxonomic profiling of FR02 metagenomes has been described elsewhere 100,106 . Briefly, raw 548
shotgun metagenomic sequencing reads were mapped using the k-mer-based metagenomic 549 classification tool Centrifuge 107 to an index database custom-built to encompass reference 550 genomes that followed the taxonomic nomenclature introduced and updated in the GTDB 551 release 89 93-95 . This implies that unless specified otherwise, all taxonomic names in our study 552 refer to their nomenclature in GTDB, which can be related to the original NCBI nomenclature 553 using the GTDB database server: https://gtdb.ecogenomic.org/taxon_history/. 554 555
Gut microbial composition was represented as the relative abundance of taxa. For each 556 metagenome at phylum, class, order, family, genus and species levels, the relative abundance 557 of a taxon was computed as the proportion of reads assigned to the clade rooted at this taxon 558
among total classified reads. The relative abundance of a taxon with no reads assigned in a 559 metagenome was considered as zero in the corresponding profile. For the purpose of this 560 association study and because of reduced accuracy and power when considering rare taxa, we 561 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this this version posted September 13, 2020. . focused on common and relatively abundant microbial taxa, defined as prevalent in >25% 562 studied individuals, and defined with at least 10 mapped reads per individual. For the purpose 563 of association, and as previous studies have reported that only some microbial taxa are 564 inheritable 108 , we also removed taxa with zero SNP-heritability. This filtering resulted in a 565 microbial dataset composed of a total of 2,801 taxa, including 59 phyla, 95 classes, 187 orders, 566 415 families, 922 genera and 1,123 species. 567 568
Taxonomic profiles derived from sequencing data are by nature compositional because of an 569 arbitrary total imposed by the instrument 109 . The compositional data of microbial taxa is not 570 independent and can lead to inappropriate use of linear regression. To overcome this artificial 571 bias, all relative abundance values were transformed by centre-log-ratio (CLR) 110 . CLR 572 transformed data can vary in real space and better fit the normality assumption of linear 573 regression. To minimize the impact of zeros, the reads count profiles were shifted by +1 before 574 the transformation. This process was performed using the R package compositions. When 575 visually comparing relative abundances in groups of individuals throughout the manuscript, we 576 used untransformed relative abundances, for better interpretability. Alpha (Shannon index) and 577 beta (Bray-Curtis distance) diversity were calculated at genus level used functions in the R 578 package vegan. 579 580
Genome-wide association analysis 581 582
The protocol followed in this study was described elsewhere 111 included as covariates. Age, gender, and genotyping batch were adjusted. As no genetic variant 592 was reported to have large effect size on gut microbiota, statistic estimates were based on 593 infinitesimal model which assumes small non-zero effect for large number of genetic variants. 594
To identify independent associations, GCTA-COJO 114 was used to conduct approximate 595 conditional and joint analysis using individual genetic data. Window size was set to 10 Mb, 596
assuming SNPs on different chromosomes or more than 10 Mb distance are uncorrelated. The 597
resulting effect size (beta coefficient) indicated the number of standard deviation changes of a 598 taxon's CLR transformed abundance corresponding to one effective allele increase of SNP. 599 600
As microbes interact non-independently with each other in the gut, as part of larger ecological 601 and functional communities, matSpDlite 115,116 was used to estimate the number of independent 602 tests based on eigenvalue variance, the larger the eigenvalue variance the smaller the number 603 of effective tests. The number of independent tests was 1,328 for 2,801 tested taxa. We used 604 this information to calculate a Bonferroni-adjusted study-wide significant level for significant 605 associations, which was set to 5×10 -8 /1328=3.8×10 -11 . A genome-wide significant threshold 606 was set as 5×10 -8 . 607 608
Prediction of ABO blood groups and secretor status 609 610
SNP-based typing of ABO histo-blood group was performed. A combination of four SNPs 117 611 was used for the prediction, and a 98% concordance with phenotypically typed ABO histo-612 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. . blood group has been reported for this method 5 . For blood group allele A, the two different 613 types A1 and A2 were predicted by rs507666 and rs8176704 respectively. Blood group allele B 614 was inferred from rs8176746 and blood group allele O was predicted by rs687289. As the 615 combination of these SNPs are exclusive, no haplotype information was needed. To validate 616 the accuracy of prediction, we compared it with the prediction using a different combination of 617 SNPs 77 . The two predictions were highly consistent, with over 99.9% concordance. In addition, 618 the distribution of ABO groups was consistent with the population distribution found in public 619 database. Secretor status was predicted by the genotype of FUT2 variant rs601338, where AA 620 or AG genotypes are secretors and GG genotypes are non-secretors. An 100% concordance 621 between the variation in rs601338 and secretor status was reported in a study on Finnish 622 individuals 118 . 623 624
Bidirectional two-sample Mendelian randomization (MR) analysis 625 626
Causal relationships between diseases and gut microbiota were investigated at genus and 627 species levels only to maximise interpretability. In total, 213 species and 148 genera associated 628 with at least one variant at genome-wide significant level (p<1×10 -8 ) were included. GWAS 629 summary results were collected for 46 diseases from MR-Base 119 (Table S4) . These included 630 12 autoimmune or inflammatory diseases, 9 cardiometabolic diseases, 13 psychiatric or 631 neurological diseases, cardiovascular diseases, 4 bone diseases and 8 cancers. For disease with 632 more than one GWAS records, the record with the largest sample size was kept. 633 634
Bi-directional causal inference was performed as follows to infer causal effects of microbial 635 abundance variation (exposure) on disease risk (outcome), and of disease (exposure) on 636 microbial abundance levels (outcome). To select the SNP instruments for microbial exposures 637 in our study, we followed recommendations from a previous study showing that associated 638
SNPs below a significance threshold of p<1×10 -5 had the largest explained variance on 639 microbial features 120 . For each taxon, GCTA-COJO was used to perform a conditional analysis 640 to select independently associated SNPs at p<1×10 -5 . SNP instruments for disease exposures 641
were selected at genome-wide significant threshold (p<5×10 -8 ). Subsequently LD-clumping 642 with a strict threshold (r 2 <0.001 in 1000G EUR within 10 Mb windows) was conducted to 643 select independent instruments with the lowest p values for taxa and diseases, respectively. 644 645
Effective alleles of all genetic variants were oriented to the risk-increasing alleles of exposures. 646
For each inference, five different MR methods were used to estimate the causal effect: (1) 647 inverse variance weighted (IVW) 121 , (2) weighted median 122 , (3) simple mode 123 , (4) weighted 648 mode 123 and (5) MR-Egger 124 . IVW is the most sensitive method which requires all 649 instruments are valid. But in reality, it is hard to verify that no any genetic instrument violates 650 any instrumental assumptions. Weighted median only requires at least half of the instruments 651 are valid, making its inference robust to the cases where some instruments violating the 652 assumptions. Simple mode and weighted mode rely on the largest group of similar instruments, 653
reducing the effects of other instruments especially outliers. MR-Egger allows instruments 654 having non-zero pleiotropy and provides way to test and estimate the pleiotropy effect in 655 addition to causal estimate. As these methods are based on different assumptions, the 656 consistency among them indicates a credible estimate 125 , even if discrepancy in these methods 657
does not necessarily suggest the absence of causality. A predicted causal estimate was deemed 658
interesting in our study if: (1) it reached a nominal p<0.05 for at least three of the five tested 659 methods (Table S7) , (2) directionality testing supported the causal direction, and (3) no 660 significant casual effect in the reverse direction. In addition, MR-PRESSO 126 was used to 661 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. . formally detect and correct for the pleiotropic outliers. Analyses were conducted using the R 662
package TwoSampleMR 119 . 663 664
Cox proportional hazards regression 665 666
Cox proportional hazards regression was conducted to test the association between baseline 667 abundance of gut microbe and incident major depression (16 years follow-up, n=181 incident 668 events). Microbial abundances were CLR-transformed and standardized to zero-mean and unit-669
variance. The Cox models were stratified by sex and adjusted for age and log-transformed 670 BMI, with time-on-study as the time scale. Participants with prevalent major depression at 671 baseline were excluded. R function coxph() in the R package survival was used for this 672 analysis. 673 674
Profiling of carbohydrate-active enzymes (CAZymes) in bacterial genomes 675 676
The standalone run_dbCAN2 v2. We used GreenAlgorithms v1.0 131 to estimate that the main computational work in this study 692 had a carbon impact of at least 531.94 kg CO 2 e, corresponding to 560 tree-months. As a 693 commitment to the reduction of carbon emissions associated with computation in research, we 694 consequently funded planting of 30 trees through a local Australian charity, which across their 695 lifetime will sequester a combined estimated 8,040 kg CO 2 e, or 15 times the amount of CO 2 e 696 generated by this study. 697 698 699 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted September 13, 2020. . 1036 visualization of the taxonomic diversity of genome-wide associated microbial taxa. The central root of the tree represents the Bacteria domain, the first 1037 represents phylum, the second connected node class, the third order and the fourth family. Every node represents at least one associated taxa in the GWAS 1038 significance level. The three smaller trees on the right highlight all taxonomic groups containing at least one taxon identified as associated with the LCT-M 1039 MED13L loci (blue edges and nodes denote taxa associated at study-wide significance level and purple edges and nodes denote taxa associated at genome-w 1040 level). The main tree is annotated to indicate phyla harbouring >10 distinct genome-wide associated taxa, as well as previously described keystone taxa. Table S1 ). The mpared groups of (n=5,941).
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 