Sequencing and microbiome data processing

For 2482 FGFP individuals, the V4 region of the 16S rRNA gene was amplified using the 515F/806R primer pair (GTGYCAGCMGCCGCGGTAA and GGACTACNVGGGTWTCTAAT, respectively), modified to contain a barcode sequence between each primer and the Illumina adaptor sequences to produce dual-barcoded libraries³². Size selection was performed using Agencourt AMPure to remove fragments below 200 bases. Sequencing was carried out on the Illumina HiSeq platform at the VIB Nucleomics core laboratory (Leuven, Belgium) with 500 cycles (sequencing kit HiSeq-Rapid SBS kit, version 2), producing 2x 250bp paired-end sequencing reads. After de-multiplexing with sdm as part of the LotuS pipeline³³ without allowing for mismatches, fastq sequences were further analyzed per sample using DADA2 pipeline (v. 1.6)¹⁵. In brief, after inspecting quality, sequences were trimmed to remove the primers and the first 10 bases after the primer, keeping only 200 bases and 130 for the R1 and R2 files, respectively. After merging paired sequences and removing chimeras, compositional matrices for each taxonomical level were carried out using the Ribosomal Database Project (RDP) training set ‘rdp_train_set_16’. Each sample was randomly down-sampled, also known as a rarefaction step, reducing the microbiome to a size of 10,000 reads. Classifications with low confidence at the genus level (<0.8) were organized in an arbitrary taxon of “unclassified_group”.

Microbiome trait preparation

The DADA2 pipeline yielded count data for 499 taxa across five levels of the microbiota phylogeny from phylum to genus (Extended Data Fig. 4a and 4b). Quality control on an individual level was performed by constructing an initial multi-dimensional scaling plot using Bray Curtis distances derived from the vegdist() function of the vegan package, and the argument method=”bray”, followed by a Kruskal’s nMDS using the function isoMDS() from the MASS package using rarefaction data from genera level counts. Two individual samples were identified as outliers in their genus-level microbiome profiles in this analysis (Supplementary Fig. 4), with the outlier cut-off set at greater than or less than five standard deviations (SDs) from the population mean of both nMDS axes. These two individuals were removed from all subsequent analysis including all association analyses.

α- and β-diversity statistics were estimated using the rarefaction data for all 288 genera level taxa counts. The α-diversity statistics used are (1) the number of genera observed, defined as the number of non-zero counts observed across all genus-level taxa, (2) Shannon diversity as calculated with the function diversity() in the vegan package, using the arguments index = “shannon”, MARGIN = 1, base = exp(1), and finally (3) Chao diversity estimated using the estimateR() function in the vegan package. β-diversity was estimated using a two-axis non-metric multidimensional scaling (nMDS) analysis using the function metaMDS() from the vegan package and the arguments, distance = “bray”, k=2, try=20, trymax=50, and trace=FALSE. The stress of the nMDS was 0.209 (Extended Data Fig. 4c). Enterotyping (or community typing) based on Dirichlet multinomial mixtures (DMM) as previously described³⁴ was performed using the DirichletMultinomial package in R. This analysis was performed on the FGFP-genus-level relative abundance matrix rarefied to 10,000 reads (Extended Data Fig. 4c).

For association analysis, we retained any taxa that met two criteria - (1) make up ≥5% of the reads for at least 1 individual and (2) have ≥15% of individuals with non-zero data (Extended Data Fig. 5a; Tables S15 and S16). In total, 139 taxa across all phylogenetic levels met these criteria. However, given that lower phylogenetic level count data can be exactly, or very close to, the same as count data at higher phylogenetic levels, we wanted to eliminate any statistical redundancies in the mGWAS. As such, we estimated Pearson correlation coefficients among all taxa, and any taxa pair with a correlation coefficient greater than 0.985 had the higher taxon level removed from association analyses. Seventy-three of the taxa exhibited such a correlation with at least one other taxa (Extended Data Figs. 5b and 5c). After removing higher-level taxa, 92 taxa remained for association analyses. The FGFP data set was used to identify these 92 taxa.

Given the ecological, observational count nature of 16S data, many individuals contain zero counts for some taxa. As such, a common feature of this data is zero-inflation, which can prove problematic for data transformation and linear modeling (Supplementary Fig. 1 and 2). To account for this possibility, we identified those taxa that should go through a two-step hurdle analysis that includes a presence/absence (P/A) association analysis and a zero-truncated abundance (AB) mGWAS. To do so, the proportion of individuals that were zero (absent) for each taxon was estimated, and those with greater than 5% zeros were pushed through the hurdle analysis. First, all non-zero counts were turned into 1’s for the binary P/A mGWAS and second, all zero counts were turned into NAs for the zero-truncated AB mGWAS. Sixty-two of the 92 retained taxa fit these criteria and were processed in this manner. The other 30 MTs were treated as simple abundance phenotypes and also denoted as AB. We note that the outcome of this procedure of course is a variety of different sample sizes in both the zero-truncated abundance phenotypes and between the number of absent individuals in P/A traits among taxa. This is an outcome that will introduce variability in power among the mGWAS performed here. Again, we note that only the FGFP data set was used to identify model type for each taxon.

In preparation for the association analysis, each abundance phenotype was rank normal transformed using the rntransform() function from the GenABEL³⁵ package and fit to a multivariate linear model, using the function lm() from the stats package, with the following covariates: the extraction type (drill or cut), the extraction year, the aliquot year (for 16S rRNA sequencing), the person performing the aliquot, the library preparation plate, genotype derived principle components 1-10, genotype predicted sex, and age. Residuals from this model were extracted using the function residuals() from the stats package and used in univariate linear modeling in the association analysis with genotypes, details below. Shapiro-Wilk W statistics for the raw and residualized data distributions can be found in Tables S17–S18. Analysis and preparation of the microbial trait data was carried out in R version 3.4.1 “Single Candle”³⁶.

Observational analysis

To identify biological phenotypes that may be influenced by gut microbiome variation in the FGFP data set (generalized) linear models, as described above, where fit with age, sex, and the top ten principle components as covariates along with each of the microbial traits (MT) analysed in the GWAS (results not included, but available on request). Human phenotypes include blood lipids, glycemic traits, anthropomorphic traits, diet and Bristol stool score. To identify laboratory batch variables that may have influenced 16S microbiome variation, we set all available variables as dependent variables in univariate analysis, with each MT set as the response variable to identify those that should be included as covariates in the GWAS. Batch variables that exhibited independent effects on at least one MT are the extraction type (drill or cut), the extraction year, the aliquot year (for 16S rRNA sequencing), the person performing the aliquot, and the library preparation plate. Further information and results from these analyses can be found in Supplementary Information (Supplementary Fig. 5).

Genotyping

A total of 2646 FGFP individuals were processed on two different arrays - the Human Core Exome v1.0 (N = 576 samples) and the Human Core Exome v1.1 (N = 2112 samples), which included repeat measurements. Allele calling was performed using GenomeStudio v2.0.4 following manufacturers default recommendations. While running GenomeStudio Log R Ratio (LRR) and B Allele Frequency (BAF) statistics were also extracted for copy number variant (CNV) calling with PennCNV³⁷. Unmapped and duplicate positions were removed, and the two batches were merged into a single data set resulting in 545,535 overlapping markers.

Variant quality control (QC) steps included the removal of unmapped variants (n = 777), duplicated sites (n = 6899), variants with >5% missingness (n = 3445), those with Hardy-Weinberg equilibrium deviations p-values < 1×10^-05 - after accounting for relatedness (n = 404), those with ambiguous alleles (n = 12,095), and those that are tri-allelic or allele flip errors (n = 1026). A total of 509,886 variants remained after QC. Sample QC included a cross check between genetically predicted sex and available sex information (117 mismatches), removal of array failed samples (n = 5), samples with >5% variant missingness (n = 53), samples with heterozygosity ± three SDs from the population mean estimate (n = 33), the removal of cryptically related (relatedness > 0.025) samples (n = 262) using the function in rel-cutoff in plink 1.9³⁸, and those with genotypic discordance among replicates (n = 8). Data was then merged with that from phase three of the 1000 Genomes Project to identify those individuals exhibiting ancestry components from populations outside of Western Europe (n = 34), using principal component analysis (PCA, Supplementary Fig. 6). After QC 2293 individuals remained (Supplementary Fig. 7), 2257 of which were retained given the availability of microbiome data.

FGFP genotype data was phased using SHAPEIT3³⁹ and imputed with IMPUTE2⁴⁰ using UK10K and all 1000 Genome Project phase 3 samples as the reference panel⁴¹. Following imputation the 39,168,681 SNPs were filtered to retain only those sites with a minor allele frequency greater than or equal to 1% and with an imputation quality score (INFO) greater than or equal to 0.3, as estimated with qctool v2.0 -snp-stats (www.well.ox.ac.uk/~gav/qctool). In total 7,711,310 SNPs were retained for the FGFP mGWAS. A flowchart of this genotyping quality control steps is available in Extended Data Fig. 6a.

To acquire insertion deletion variants, the data was also phased and imputed to Genome of Netherlands reference panel using Impute2 v2.3.0⁴². All indels were isolated from this imputed data set and run in addition to the imputation data set from above.

Copy number variants (CNVs) were called with PennCNV v1.0.4³⁷ using the perl script detect_cnv.pl. Cleaning was performed with the perl script clean_cnv.pl, and filtered with the script filter_cnv.pl using the flags --numsnp 5 --length 250 -qclrrsd 0.35 -qcnumcnv 716. Unique CNVs were defined by unique base pair start and stop locations. In total 35,020 unique CNVs were identified across the FGFP sample; 949 CNVs were shared across 1% or more individuals. Global CNV burden was estimated for each individual as the number of CNVs that do not equal the copy number count of 2. Insertions (>2) and deletions (<2) were treated the same. Regional CNV burden was calculated in sliding windows of 200 kilo-bases and estimated following the same rules as for global burden estimation.

Heritability

Chip-based heritability was estimated for each microbiome presence/absence (62) abundance (92) and α-diversity metric (3) phenotype used in the association analysis, using the GCTA-GREML restricted maximum likelihood (REML) method, and a single genetic relationship matrix (GRM) as implemented in GCTA version 1.91.1beta⁴³. The GREML power calculator was used to estimate the power to detect genetic covariation in the FGFP data set (Extended Data Fig. 7)⁴⁴. For the abundance phenotypes, the residualized data used in the association analyses were used in the estimation of heritability. For binary phenotypes, the same covariates, mentioned above, were fit to the trait by GCTA. To produce the genetic relationship matrix (GRM) for running GCTA, we first identified all genotypes with an info (imputation quality) score greater than or equal to 0.9, a minor allele frequency greater than 0.05, and not deviating from Hardy-Weinberg equilibrium. Genotype probability score data was converted to hard call plink format data using qctools. SNP variation was linkage disequilibrium pruned using plink2 and the flag --indep-pairwise 50 5 0.45. Finally, the GRM was constructed using GCTA and the flags --grm-cutoff 0.025 --make-grm. In addition, for a more direct comparison with previously published studies⁷, we also performed box-cox transformations of the abundance phenotypes and regressed out our covariates.

Primary FGFP association analysis

Following genotype and microbiome QC, 2257 individuals remained, and 2223 remained after accounting for data missingness among covariates. All microbiome α-diversity, abundance and presence/absence associations analyses were performed using snptest.2.5.0⁴⁵. All abundance traits were regressed on covariates (the aliquoting procedure, the extraction year, the aliquot year (for 16S rRNA sequencing), the person performing the aliquot, the library preparation plate, genotype derived principle components 1-10, genotype predicted sex, and age) and residuals were regressed on genotype probability data in a univariate fashion, assuming an additive genetic model and using the missing data likelihood score test in snptest (snptest flags: -frequentist 1 -method score and -use_raw_phenotypes). Presence/absence mGWAS were performed using the same covariates as those described above for abundance traits in a multivariate analysis again using the same snptest settings. These primary analyses were performed as a first pass signal detection step in order to determine signals to take forward for meta-analysis and to confirm the ability of score analyses to effectively rank the expectation–maximization (em) method (Supplementary Fig. 8). Association analyses for enterotype were run using a multinomial logistic regression for categorical traits as employed by snptest.2.5.4-beta3 and the flags -frequentist add -method newml and setting the -baseline_phenotype to “Bacteroides1”. Finally, associations for β-diversity (a two axis MDS) were run using a bespoke R script and the function manova() from the stats package, in a multivariate analysis using the same covariates stated above and genotype dosages as derived by qctool v2.0. A flow chart of the mGWAS is provided in Extended Data Fig. 6b.

Data preparation in German cohorts

The FoCus¹¹ and PopGen¹⁶ cohorts were genotyped using the Illumina Omni Express + Exome array and the Affymetrix Genome-Wide Human SNP Array 6.0, respectively. Genotyping QC and imputation in these two cohorts were performed following protocols defined here: https://github.com/alexa-kur/miQTL_cookbook#chapter-2-genotype-imputation. SNPs where filtered as a MAF of 0.01 and an INFO score of 0.3, as was done in the FGFP cohort. Microbiome census data for the Kiel based cohorts, targeting the 16S V1-V2 regions, was generated as described previously¹¹. Data processing was performed using DADA2¹⁵ modified for V1-V2 (https://github.com/mruehlemann/rep-cookbook/blob/master/scripts/Seq_dada2_V12_Kiel.R) following the same standardized workflow described in the above Microbiome trait preparation section. α- and β-diversity metrics, enterotypes, and abundance measures for association analysis were calculated as previously indicated. Three genera, Escherichia Shigella, Hespellia, and Methanobrevibacter were not present in the FoCus or PopGen cohorts (Supplementary Fig. 9). As such, in all three instances, their P/A and zero-truncated AB MTs were not available for association and inclusion in the meta-analysis. Association analyses were carried out as described below.

Meta-analysis

For the purposes of the meta-analysis, all outcomes were defined as described above in the “Primary FGFP association analysis” section, however, the expectation–maximization (em) method was used, rather than the score method, to account for genotype uncertainty and given the performance of “score” at low allele frequency and phenotype/trait group size (Supplementary Information).

The meta analyses were performed using the inverse-variance fixed effects method (method 1) as implemented in the software package META⁴⁶. The imputation quality threshold for each SNP was set at 0.3. To identify loci or unique genomic regions defined by shared linkage disequilibrium, we clumped all meta-supported markers using plink, the flag –clump and the p-values derived from the em meta-analysis. A locus or index tags was only identified if their meta-analysis p-value was < 0.0001.

Beta estimations (genotypic effects) for P/A traits are defined as an increase in the log odds ratio for each additional effect allele. For AB traits, beta is defined as a change in SD units for each effect allele carried. The study-wide p-value threshold was defined as a Bonferroni correction assuming 2 million independent genetic association tests across 159 mGWAS (0.05 / (2×10⁰⁶ × 159) = 1.57×10^-10).

Phylogenetic analysis

Representative 16S rRNA gene sequences of all the genera identified were retrieved from the RDP database. Multiple sequence alignments were performed for all taxa and for genera included in the GWAS analyses using MUSCLE v.3.8⁴⁷. The alignments were used to build maximum likelihood trees using FastTree v2.1.0⁴⁸ with default parameters. iTOL⁴⁹ was used for visualizing the trees with corresponding metadata, including the number of loci, prevalence of the MT, abundance of the MT and heritability of the AB and P/A trait(s) (Fig. 1 and Table S3).

Functional annotation and enrichment

Annotation of variants of interest was carried out with the biomaRt R package⁵⁰ with additional linkage disequilibrium based annotation and enrichment analysis with DEPICT⁵¹. When using biomaRt, we referenced the (feb2014.archive.ensembl.org) Ensembl 75 archive for GRCh37/hg19 coordinates and identified all genes and the closest gene within 250 kilo-bases up- and down-stream of each polymorphism. Given that DEPICT utilizes pre-computed LD structure from genotypes derived from 1000 Genomes Project Phase 1 CEU, GBR and TSI and HapMap Project release 2 and 3 CEU data, a substantial proportion of our SNPs of interest were not represented. As such, we identified tag SNPs for our SNPs that are also present in the DEPICT data set, when possible. To do so, we extracted dosage data for all variants +/- 200kb of our variant, using qctools, and then computing r² using a bespoke R script. We kept all variants with an r² > 0.2 with our SNP of interest and then queried if any of those tag SNPs existed in the DEPICT data. If more than one was present, we kept the one with the highest r² value. Subsequently, the list of reference SNP identifiers (rs-ids) composed of our SNPs of interest, when they existed in the DEPICT data, or alternatively tag SNPs, when they were present, were run through the DEPICT framework. Each MTs was run individually using FGFP variants associated at a threshold of 1×10^-5. Results from these analyses, gene enrichment, tissue enrichment, and gene priority are available in Table S19, S20, S21, respectively.

The Gene-Tissue Expression portal (gtexportal.org) was used to determine if associated variants, specifically discussed in the text, were also expression quantitative trait loci (gtexportal.org).

When evaluating enrichment for all meta-supported variants as a unit, we extracted the closest gene (within +/- 250kb) for each variant and used GENE2FUNC (http://fuma.ctglab.nl/), and its integrated hypergeometric test, to identify tissue and pathways functional enriched, as reported in Table S12. The background set of genes, included in hypergeometric test, was 19,283 protein coding genes. The default parameter of GEN2FUNC in the FUMA platform²⁹.

All microbial traits analyzed were also analyzed under the GARFIELD⁵² framework for identification of regulatory and tissue enrichments. This analysis uses the complete GWAS data set of a trait, so pooling data from multiple traits, as done above, was not carried out. Further, the FGFP results only were used in these analyses. The results of these analyses are available in Table S22.

Finally, variants with associations that met our study-wide and genome-wide confidence thresholds as well as those that may be deemed as replicated from other studies were passed through PhenoScanner V2, an online platform to screen for genotype-to-phenotype associations, expression quantitative loci, and methylation quantitative loci from previously published genome-wide -omics association analysis²⁸. All of these results are provided in Table S11.

Applied analyses

We undertook two-sample, bi-directional Mendelian randomization³⁰ analyses to estimate potentially causal relationships between gut MTs and 11 metabolic health, inflammatory and neurological traits. These were selected a priori as they have all been repeatedly been associated with variation in the gut microbiome and have been the focus of credible and accessible GWAS studies. They include waist circumference, waist-hip ratio, BMI, and type 2 diabetes; Crohn’s disease, inflammatory bowel disease, ulcerative colitis and rheumatoid arthritis; and Alzheimer’s disease, Parkinson’s disease and major depressive disorder.

MR analyses interrogating the role of the gut microbiome on each of these outcomes were restricted to the gut MTs that had the greatest evidence of a host-genetic contribution - where independent meta-derived genetic variants reaching a genome-wide threshold of p<2.5×10^-08 were used as “instruments”. In order to assess causality in relationships from microbiome to outcome, summary statistics for these genetic variants were obtained from publicly available genome-wide summary-level data for the 11 metabolic health, inflammatory and neurological traits. For analyses assessing causality in relationships from each trait to microbiome, independent genetic variants reaching a genome-wide threshold of p<5×10^-08 in each respective GWAS were used as “instruments” for the relevant trait. Summary statistics were obtained from the current mGWAS meta-analysis.

Once all summary-level data was obtained, causal effect estimates were derived using the inverse variance weighted (IVW)⁵³ method (or the Wald ratio, if only 1 genetic variant was available) alongside sensitivity analyses including the weighted median⁵³, weighted mode⁵⁴ and MR-Egger⁵⁵ tests (if ≥3 genetic variants were available). All exploratory MR analyses were conducted using the TwoSampleMR package (https://github.com/MRCIEU/TwoSampleMR) in R version 3.4.1 with RStudio, created and provided by MR-Base (www.mrbase.org/)⁵⁶, a large-scale database of GWAS summary-level data and automated pipeline for two-sample MR analyses. Summary-level GWAS results for the 11 selected metabolic health, inflammatory and neurological traits were obtained from the following publications (indicated with pubmed ID): waist circumference and waist-to-hip ratio (25673412); BMI (25673413); type 2 diabetes (22885922); Crohn’s disease, inflammatory bowel disease and ulcerative colitis (26192919); rheumatoid arthritis (24390342); Alzheimer’s disease (24162737), Parkinson’s disease (19915575) and major depressive disorder (22472876).

In analyses interrogating the impact of each continuous (AB) MT on each outcome, effect estimates represent the SD change for continuous outcomes (waist circumference, waist-to-hip ratio and BMI) or risk ratio for binary outcomes (type 2 diabetes, Crohn’s disease, inflammatory bowel disease, ulcerative colitis, Alzheimer’s disease, major depressive disorder, Parkinson’s disease and rheumatoid arthritis) for a SD unit of AB MT phenotype. For analyses interrogating the impact of each binary (P/A) MT on each outcome, effect estimates represent the SD change for continuous outcomes or risk ratio for binary outcomes for a doubling of the genetic liability to presence (vs. absence) of each P/A MT phenotype. Results reaching a p-value threshold of (p<0.05) are presented (Table 2).

FGFP procedures were approved by the medical ethics committee of the University of Brussels–Brussels University Hospital (approval 143201215505, 5/12/2012). A declaration concerning the FGFP privacy policy was submitted to the Belgian Commission for the Protection of Privacy. Written informed consent was obtained from all participants. The FGFP was funded with support of the Flemish government (IWT130359), the Research Fund–Flanders (FWO) Odysseus program (G.0924.09), the King Baudouin Foundation (2012-J80000-004), FP7 METACARDIS HEALTH-F4-2012-305312, VIB, the Rega Institute for Medical Research, and KU Leuven. We acknowledge the contribution of Flemish GPs and pharmacists to data and sample collection. Finally, we thank all FGFP volunteers for participating in the project. RB is funded by the Research Fund–Flanders (FWO) through a Postdoctoral Fellowship (1221620N). SVS is supported by Marie Curie Actions FP7 People COFUND–Proposal 267139 (acronym OMICS@VIB). NJT is a Wellcome Trust Investigator (202802/Z/16/Z), is supported by the University of Bristol NIHR Biomedical Research Centre (BRC-1215-20011) and works within the CRUK Integrative Cancer Epidemiology Programme (C18281/A19169). SR and NJT work within a unit that receives funding from the MRC and University of Bristol (MC_UU_12013/3, MC_UU_00011/6). DH is supported by the Wellcome Investigator Award (202802/Z/16/Z). KHW is supported by the Elizabeth Blackwell Institute for Health Research, University of Bristol and the Wellcome Trust Institutional Strategic Support Fund (204813/Z/16/Z). MR was funded by the German Research Foundation (DFG) Collaborative Research Center 1182 (CRC1182; Origin and Function of Metaorganisms).