Microbiome-wide PheWAS links gut microbial SNVs to health and exposures

Cohort description

Dutch Microbiome Project

The DMP consists of 8,719 individuals and is part of the Lifelines study, a multi-disciplinary prospective population-based cohort study that utilizes a unique 3-generation design to examine health and health-related behaviours in 167,729 people living in the northern Netherlands. Lifelines utilizes a broad range of investigative procedures to assess the biomedical, socio-demographic, behavioural, physical and psychological factors that contribute to health and disease, with a special focus on multi-morbidity and complex genetics⁴⁸. Microbiome data generation for the DMP is described in Gacesa et al.¹². In brief, fresh-frozen faecal samples were collected from participants of the Lifelines cohort. Microbial DNA was isolated using the QIAamp Fast DNA Stool Mini Kit (FSK, Qiagen) by the QIAcube automated sample preparation system (Qiagen). Metagenomic sequencing was performed at Novogene, China using the Illumina HiSeq 2000 sequencer.

Lifelines-DEEP cohort

LLD is a sub-cohort of Lifelines that was established to explore the risk factors of complex diseases. In LLD, 1,539 individuals were included, with multi-layers of omics data collected at baseline in 2013 (ref. ¹³). In the current study, high-quality metagenomic sequencing data were available for 1,135 individuals (474 male and 661 female participants). The average age of LLD participants was 45.04 years (18–81 years, s.d. = 13.60 years). The average BMI was 25.26 kg m⁻² (16.67–48.56 kg m⁻², s.d. = 4.18 kg m⁻²). Microbiome data generation for the LLD is described in Zhernakova et al.⁴⁹. In brief, fresh-frozen faecal samples were collected, and DNA was isolated with the AllPrep DNA/RNA Mini Kit (APK, Qiagen, 80204). Sequencing was performed using the Illumina HiSeq platform at the Broad Institute, Boston. 1,135 metagenomic samples passed quality control.

Follow-up APK samples of LLD were collected for 338 individuals (150 male and 188 female participants) at a second time point. The duration between the two time points ranged from 3.33 to 3.92 years (mean = 3.53 years, s.d. = 0.12 years). At baseline, the mean age of participants was 48.2 years (18–80 years, s.d. = 11.7 years) and their mean BMI was 25.4 kg m⁻² (17.6–43.3 kg m⁻², s.d. = 4.08 kg m⁻²). At follow-up, the mean age was 51.7 years (22–84 years, s.d. = 11.7 years) and the mean BMI was 25.6 kg m⁻² (16.1–37.6 kg m⁻², s.d. = 4.0 kg m⁻²). Two aliquots of each stool sample were processed and used for DNA isolation and library preparation using two kits: FSK (Qiagen) and APK (Qiagen, 80204). WGS was performed on the Illumina HiSeq 2000 platform.

500-Functional Genomics cohort

The 500FG is part of the Human Functional Genomics Project (HFGP) and consists of 534 individuals. The average age of 517 500FG-FSK participants was 28.47 years (18–75 years, s.d. = 13.69 years), and the average BMI was 22.75 kg m⁻² (15.09–35.43 kg m⁻², s.d. = 2.85 kg m⁻²). The metagenomic data generation was described previously^14,15. In brief, DNA was isolated from 517 and 471 faecal samples with the FSK (Qiagen) and APK (Qiagen, 80204), and libraries were sequenced on the Illumina HiSeq 2000 platform.

300-Obesity cohort

The 300OB cohort was established by Radboud University Medical Center, Nijmegen, the Netherlands⁵⁰. In total, 298 older individuals with metagenomic data (132 female and 166 male participants) with a high BMI > 27 kg m⁻² were enrolled in 300OB. The average age of 300OB participants was 67.07 (54–81 years, s.d. = 5.39 years). The average BMI was 30.73 kg m⁻² (26.28–45.47 kg m⁻², s.d. = 3.48 kg m⁻²). All participants were included between 2014 and 2016. The metagenomic sequencing of faecal samples was done using a similar protocol and analysis pipeline to LLD.

1000 Inflammatory Bowel Disease cohort

The 1000 Inflammatory Bowel Disease cohort (1000IBD) consists of patients with IBD recruited at the specialized IBD outpatient clinic of the University Medical Center Groningen (UMCG) in the Netherlands^16,51. IBD diagnosis was made based on accepted radiological, endoscopic and histopathological evaluation. The 1000IBD cohort is 60.70% female and 39.30% male. The mean age of participants was 42.61 (18–82 years, s.d. = 14.69 years) and their mean BMI was 25.44 kg m⁻² (15–49.4 kg m⁻², s.d. = 5.05 kg m⁻²).

300-Tanzania Functional Genomics cohort

For analysis among non-Europeans, we included 300TZFG, a population cohort of 323 individuals from both rural and urban areas in the Kilimanjaro region of northern Tanzania. This study is part of the HFGP. Metagenomic data generation was described previously¹⁸. In brief, bacterial DNA was isolated using the APK (Qiagen), and libraries were sequenced on the Illumina NovaSeq 6000 platform. The average age of 300TZFG participants with metagenomic data was 33.08 years (17.76–66.43 years, s.d. = 11.25 years). The average BMI was 24.98 kg m⁻² (16.41–41.91 kg m⁻², s.d. = 4.93 kg m⁻²).

Quantification of metabolites

Bile acid quantification

The levels of 15 bile acids and C4 concentrations in fasting plasma were quantified by liquid chromatography–tandem mass spectrometry (LC–MS/MS) procedures, as previously described^52,53. In brief, for sample preparation, 25 µl plasma was mixed with 250 µl internal standard solution to precipitate proteins. Samples were centrifuged at 15,800g, and the supernatant poured into a clean glass tube. The fluid was evaporated under nitrogen at 40 °C. Before measuring, samples were reconstituted in 200 µl 50% methanol in water. For the quantitative determination of bile acids, we used a Nexera X2 Ultra High-Performance Liquid Chromatography system (Shimadzu), coupled to a SCIEX QTRAP 4500 MD triple quadrupole mass spectrometer (SCIEX) (UHPLC–MS/MS). The LC–MS/MS system is controlled by Analyst MD 1.6.2 software. Bile acids were separated with an ACQUITY UPLC BEH C18 Column (1.7 µm, 2.1 mm × 100 mm) equipped with an ACQUITY UPLC BEH C18 VanGuard Pre-Column (1.7 µm, 2.1 mm × 5 mm) (Waters). Separation was achieved in 28 min using 10 mM ammonium acetate in 20% acetonitrile (mobile phase A) and 10 mM ammonium acetate in 80% acetonitrile (mobile phase B), flow 0.4 ml min⁻¹.

The proportions of 15 bile acids (with suffix ‘_p’) were calculated by dividing by total bile acid concentration. Additionally, eight different total bile acid concentrations and ratios were calculated: (1) total bile acid (Total_BAs) = sum of all bile acid concentrations; (2) total primary bile acid (Total_primary_BAs) = sum of all primary bile acid concentrations; (3) total secondary bile acid (Total_secondary_BAs) = sum of all secondary bile acid concentrations; (4) ratio of secondary bile acids to primary bile acids (Secondary_primary_ratio) = Total_secondary_BAs/Total_primary_BAs; (5) ratio of CA to CDCA concentrations (CA_CDCA_ratio) = (CA + TCA + GCA)/(CDCA + TCDCA + GCDCA); (6) ratio of unconjugated bile acid to conjugated bile acid concentrations = (CA + CDCA + DCA + LCA)/(TCA + GCA + TCDCA + GCDCA + TDCA + GDCA + TLCA + GLCA); (7) ratio of dehydroxylated CA to deconjugated CA concentrations (CA_dehydro_deconju_ratio) = (DCA + TDCA + GDCA)/(CA + TCA + GCA); and (8) ratio of taurine-conjugated bile acid to glycine-conjugated bile acid concentrations (Taurine_glycine_ratio) = (TCA + TCDCA + TDCA + TLCA)/(GCA + GCDCA + GDCA + GLCA).

Short-chain fatty acid quantification

The stool level of short-chain fatty acids was measured using gas chromatography–mass spectrometry, as previously described^54,55.

Quantification of trimethylamine N-oxide and its precursors

Trimethylamine N-oxide and its precursors betaine, carnitine, deoxycarnitine, and choline were measured using LC–MS/MS, as previously described^56,57. In brief, an online cleanup of the sample was carried out using an analytic column containing a C₁₈ stationary phase, prior to its introduction into the mass spectrometer. This stationary phase did not retain the analytes, but it did retain significant matrix interferences, such as (phospho-)lipids⁵⁸.

Untargeted blood metabolome quantification

Untargeted metabolic profiling of fasting plasma samples was conducted at General Metabolics using flow-injection–time-of-flight mass spectrometry on an Agilent 6550Q-TOF system, as previously described^59,60,61.

Quantification of biological age parameters

Telomere length measurement

Telomere length was assessed using an automated multicolour flow–fluorescence in situ hybridization (Flow-FISH) technique as previously described^62,63. In brief, white blood cells were obtained by lysing erythrocytes in whole blood with NH₄Cl. These white blood cells were then combined with bovine thymocytes, which have a known telomere length and serve as an internal reference. The cells underwent denaturation in formamide at 87 °C, followed by hybridization with a fluorescein-labelled (CCCTAA)₃ peptide nucleic acid probe specific to telomere sequences. They were then counterstained with the DNA dye LDS751. Flow cytometry was used to measure the fluorescence intensity in granulocytes, total lymphocytes, and lymphocyte subsets (identified by antibodies specific for CD20, CD45RA and CD57). This fluorescence intensity, compared to internal control cells and unstained controls, was used to determine the median telomere length from duplicate measurements.

sjTREC measurement

Signal joint excision circles produced during T cell development (sjTREC) expression was quantified using quantitative PCR (qPCR) as previously described⁶³. In brief, DNA was extracted from whole blood using proteinase K digestion, phenol extraction, and isopropanol precipitation. We then used TaqMan qPCR on a ViiA 7 system (Life Technologies) to quantify sjTREC expression. The qPCR was performed on 75–100 ng of DNA in a 12 μl reaction mix with 18 µM of sjTREC primers (5′-TCGTGAGAACGGTGAATGAAG-3′ and 5′-CCATGCTGACACCTCTGGTT-3′) and albumin primers (5′-TGAACAGGCGACCATGCTT-3′ and 5′-CTCTCCTTCTCAGAAAGTGTGCATAT-3′) and 5 µM probes for sjTREC (5′-(FAM) CACGGTGATGCATAGGCACCTGC-3′ (TAMRA)) and albumin (5′-(FAM) TGCTGAAACATTCACCTTCCATGCAGA -3′ (TAMRA)). The PCR protocol was 50 °C for 2 min, 95 °C for 15 min, followed by 45 cycles of 95 °C for 15 s and 60 °C for 30 s. Reactions were done in duplicate. Fluorescence from probe cleavage was measured at the threshold cycle (C_T). Valid results had a duplicate s.d. ≤ 1. Normalized sjTREC expression was calculated from the difference between average C_T values of albumin and sjTREC.

Metagenomic sequencing and quality control

We removed host genome-contaminated reads and low-quality reads from the raw metagenomic sequencing data using KneadData (v0.7.4), Bowtie2 (v2.3.4.3)⁶⁴ and Trimmomatic (v0.39)⁶⁵. In brief, the data-cleaning procedure includes two main steps: (1) filtering out human genome-contaminated reads by aligning raw reads to the human reference genome (GRCh37–hg19) and (2) removing adaptor sequences and low-quality reads using Trimmomatic with default settings (SLIDINGWINDOW:4:20 MINLEN:50).

Quantification of reference genome abundance

We estimated the abundances of reference genomes in metagenomic samples using CoverM (v0.7.0)⁶⁶ based on the identical reference genomes used in metagenomic SNV detection.

Metagenomic structural variant detection

Structural variants were defined as within-species highly variable genomic fragments (>1 kb) across microbial genomes. Structural variants of the gut microbiome were detected from cleaned metagenomic reads using SGVFinder (v2.0)⁴ based on the identical reference genomes used by metagenomic SNV detection.

Metagenomic SNV detection

We characterized the SNVs for metagenomic samples using GT-Pro (v1.0.1) based on the pre-constructed sck-mer database¹¹. In brief, a catalogue of 104 million biallelic SNVs in the core genomic region of 909 gut microbial species was initially constructed based on Unified Human Gastrointestinal Genome (UHGG) assemblies, the genotypes of the gut microbial SNVs were tagged by species-specific k-mers, and a k-mer-based gut microbiome metagenomic SNV panel was constructed. GT-Pro directly conducted genotyping for gut microbiota from metagenomic sequencing data by k-mer-searching and matching the genotype k-mers of the metagenomic SNV panel. GT-Pro output was further converted to SNV dosage format, with genotype indicating the frequency of the alternative allele.

Quality control of the SNV profile

We calculated allele frequencies based explicitly on the average reference or alternative allele proportions across samples:

For a given SNV across n samples, let:

R_i = relative abundance (proportion) of the reference allele in sample i (for example, 0.7)

A_i = relative abundance (proportion) of the alternative allele in sample i (for example, 0.3)

These values typically sum to 1 per sample if only two alleles are considered.

Then:

Reference allele frequency (RAF):

Alternative allele frequency (AAF):

MAF:

We then performed quality control in each cohort separately, filtering out SNVs with a high missing rate (>90%) or low MAF (<0.05) to get the cohort-specific common SNVs. We next counted the number of SNVs retained for each species in each cohort and removed species with a low common SNV number (<1,000). The manipulation and quality control of genotype datasets were conducted using PLINK2 (v2.0.0-a.6.9LM)⁶⁷.

Annotation of microbial genes and SNVs

The genes of gut microbial reference genomes were predicted and annotated using Bakta⁶⁸ (v1.10.1, database v5.1). The functional effects of the common SNVs were annotated using SnpEff (v5.1f)²¹.

Between-sample genetic relatedness and population genetic structure

For each species, we calculated the species-specific between-sample genetic dissimilarity D_i,j measured by a Manhattan distance scaled by the number of shared SNVs within each sample pair using the parallelDist() function in R package parallelDist (v0.2.6). Sample pairs sharing <1,000 SNVs were skipped for the genetic dissimilarity calculation, and the corresponding D_i,j were marked as missing values. To maximize the size of the retained sample with no missing value, we removed the sample with the highest number of missing values from the matrix in each round until there are no missing values in the matrix. We further calculated the between-sample genetic relatedness, R_i,j = 1 – D_i,j, and the matrix M_R containing R_i,j was used as the representative of population genetic structure. We performed principal coordinate analysis on the genetic dissimilarity matrices M_D using the pcoa() function from R package ape (v5.7-1) with the parameter correction = “cailliez”.

Linkage disequilibrium calculation

To depict the linkage disequilibrium (LD) structure overview of each species, we randomly extracted 1,000 SNVs from the cleaned genotype dataset using plink2 (v2.0.0-a.6.9LM) with the parameters -thin-count 1000. We then calculated the LD as a squared correlation for all pairs of extracted SNVs using plink2 with the parameters -r2 inter-chr yes-really -ld-window-r2 0. We visualized the whole-genome LD structure using LDBlockShow (v1.40)⁶⁹ with the parameters -SeleVar 2 -NoShowLDist 100000000 -OutPdf -OutPng.

Quasi-phaseability assessment

The co-existence of multiple strains of a species within the host may bias the position of the sample in the population genetic structure and decrease the power of association analysis. For a given species, samples in which > 80% of the SNVs have the same allele in > 90% of the reads were considered quasi-phaseable, meaning the sample is very likely to be dominated by one strain of the corresponding species. The criterion and concept of quasi-phaseability were described before^38,70.

SNV-based metagenome-wide association study and meta-analysis

We first manually integrated and curated the phenotypes from different cohorts and datasets to ensure the consistency of units. Each phenotype was analysed with a LMM that contained: (1) fixed effects for corresponding species abundance, age, sex and sequencing depth, and (2) a random-effect term capturing within-species population structure, with the variance–covariance of this random effect specified by the sample-by-sample kinship matrix (M_R). When sex or age was the trait analysed, that variable was omitted from the covariate list. The LMM-based analysis was conducted using EMMAX (version beta-07Mar2010)⁷¹ with the parameters emmax -d 10 -t genotype.file -p phenotype.file -c covariate.file -k kinship.file -o output -Z. The association analyses were performed for each cohort and each species separately, with only the quasi-phaseable samples of the corresponding species included for the association analysis. For the SNV–phenotype pairs available in two or more independent cohorts (DMP, LLD1, 500FG, 300OB, IBD and 300TZFG), a follow-up meta-analysis was performed in the independent datasets (excluding the LLD2 and 500FG-APK datasets) using the Metal software (v2020-05-05)⁷² with default options (weighting cohort-based P values according to sample size). To ensure robustness and further reduce the chances of false discovery of identified associations, we conducted a permutational LMM with 999 times permutation following the Monte Carlo procedure within each cohort to assess the robustness of significant associations. Regional association plots showing LD structures on the gene regions covering the association signals were generated using LDBlockShow (v1.40)⁶⁹.

Determination of significance cut-off

The threshold of study-wide significance was determined by the total number of tested phenotype–SNV pairs (922,497,753) using the Bonferroni method: 0.05/no. of phenotype–SNV pairs = 5.42 × 10⁻¹¹. Similarly, suggestive significance was determined using Bonferroni correction based on the total number of tested SNVs in each species and for each phenotype. Thus the suggestive significance did not control for the number of species and phenotypes involved in the study.

Linkage disequilibrium-based grouping

The SNVs could be correlated due to LD in the genomes. To count the number of independent association signals, we further conducted a grouping analysis that grouped interrelated SNVs by LD and picked the most significant SNVs within each LD block as representatives or tags of the independent signals (between-SNV r² > 0.01, window size = 10,000 kbp) using PLINK2 (v2.0.0-a.6.9LM)⁶⁷. The primary and secondary significance cut-offs used in PLINK2 were study-wide significance (5.42 × 10⁻¹¹) and species-wide significance (species-specific) cut-offs. We visualized the local LD structure using LDBlockShow (v1.40)⁶⁹ with the parameters -SeleVar 2 -NoShowLDist 100000000 -OutPdf -OutPng.

Power analysis

We performed a post-hoc approximate power analysis for the tag SNVs of significant associations using the genpwr package (v1.0.4) in R. The proportion of phenotypic variance explained (R²) was estimated from the observed effect size, standard error, observed sample size and MAF using the pheno_var_explained() function. Power was then calculated with power.calc.linear() under an additive genetic model and a study-wide significance threshold of α = 5.42 × 10⁻¹¹. To evaluate the effect of sample size, power was also computed across fixed sample sizes (n = 100, 500, 1,000, 5,000 and 10,000) for the tag SNVs of significant associations.

Protein 3D structure prediction and protein structure-based analysis

As AlphaFold3 supports structure prediction for protein complexes, the three-dimensional (3D) structures of the OPP complex were modelled using AlphaFold3 (ref. ⁷³) and visualized in ChimeraX (v1.10.1)⁷⁴ and PyMOL (v2.5.7). However, according to the terms of use of the AlphaFold3 server, its output cannot be used in docking or screening tools (https://alphafoldserver.com/about). We therefore used local AlphaFold2 for monomer protein structure prediction and conducted in silico docking and other downstream analyses. To explore interactions between OppA/SarA_2 and glutathione, the protein structure of OppA/SarA_2 was first predicted using AlphaFold2 (v2.3.1)⁷⁵ for both the wild-type and mutated protein sequences under monomer settings. The default reference databases and a template cut-off date of 14 May 2020 were used to guide the structure generation. For each experiment, five models were produced and further refined using graphics processing unit-based AMBER relaxation. The best model, based on pLDDT confidence scores (mean ≥ 90), was selected for downstream analysis, excluding the low-confidence N-terminal residues. The glutathione structure was downloaded from PubChem (https://pubchem.ncbi.nlm.nih.gov/) and converted into PDB format. Ligand-binding pockets in the OppA/SarA_2 structure were predicted using P2Rank (v2.5)⁷⁶, focusing on residues with high ligand-binding probabilities (>0.6) from top-scoring pockets (probability > 0.85). These residues were then used as active sites in HADDOCK (v2.4)^77,78, a docking platform that generated 200 structural complexes for each variant. The top four docked poses from the best-scoring cluster were selected and energy-minimized in water. Binding free energies (ΔG) were estimated using PRODIGY-LIG (v1.1.4)^79,80, and PLIP (v2025)⁸¹ was used to identify non-covalent interactions. As a docking benchmark, a literature-derived^82,83 tripeptide panel representing basic (KAK, KWK, KKK), neutral (AAA, GFG, GAG) and acidic (EEE, DED) chemistries was docked to the same OppA binding site using the same set of active residues defined for the glutathione docking in HPEPDOCK (v2.0)⁸⁴. All the resulting complexes were refined through water minimization in HADDOCK (v2.4)⁷⁷, and binding free energies were estimated with PRODIGY-LIG (v1.1.4)⁷⁹.

The protein structural property analysis was conducted using an internal pipeline (https://github.com/GArcho44/snvs-annotation-pipeline). In brief, to interpret effect of SNVs at the protein structure level, only missense variants in coding regions were retained. Functional gene annotations from Bakta (v1.10.1)⁶⁸ and UniProt (https://www.uniprot.org/)^85,86 enabled linkage of genes to protein sequences within UniRef100 and consequently to structure entries in the AlphaFold Database (AlphaFoldDB, v2022)²⁷. Recognizing that many genes lacked direct matches in UniRef100, we applied a homology-based approach to expand structural coverage. Gene coding sequences were extracted using GffRead (v0.12.7)⁸⁷ and translated to proteins with EMBOSS transeq (v6.6.0)⁸⁸. The resulting sequences were queried against AlphaFoldDB and ESM Metagenomic Atlas (version 1 November 2022, https://esmatlas.com/)²⁶ using DIAMOND BLASTP (v2.0.15)⁸⁹, retaining only full-length matches with ≥98% identity. When multiple matches occurred, the best match was selected by identity score, preferring AlphaFoldDB entries when there was a tie. Protein accessions (UniProtKB or MGnify) were then used to retrieve high-confidence structural models (mean pLDDT ≥80) via custom application programming interfaces. Multiple per-residue properties were measured on these curated protein models to functionally annotate the structural properties of SNVs. Relative solvent accessibility (RSA) was calculated using DSSP (v3.0.0)^90,91, normalizing surface exposure per residue against empirical maxima. Ligand-binding pockets were detected using fpocket (v4.2.1)⁹² and rescored with fpocket_PRANK to boost detection sensitivity. Binding residues were defined from pockets with druggability score > 0.5. Each residue’s distance to ligand (DTL) was computed based on the minimal Euclidean distance to the closest pocket centre. Protein–protein interaction interfaces were predicted using PeSTo (model version i_v4_1)⁹³. Residues with interface confidence scores > 0.5 were grouped into clusters based on spatial proximity (Cα–Cα < 10 Å), and only high-confidence clusters (average score > 0.8, pLDDT > 0.7 and ≥6 residues) were retained. Each residue’s distance to protein-protein interface (DTPPI) was calculated based on Euclidean distance to closest residue of predicted interface.

The 3D structure of the marR homologue dimer was modelled using AlphaFold3 (ref. ⁷³) and visualized with ChimeraX (v1.10.1)⁷⁴. Protein-protein interface analysis, specifically in silico alanine scanning, was performed using the mCSM-PPI2 web server⁹⁴.

Identification and annotation of glutathione transport system gene cluster

The original gene annotation was conducted using Prokka (v1.13.3)⁹⁵, and the age-association signal was identified in the gene annotated as sarA (UniProtKB: P31306), which is located next to oppB, oppC, oppD and oppF. The sarA gene encodes a periplasmic oligopeptide-binding protein, SarA, and contains an oppA gene family domain. According to the literature^96,97, OppABCDF is a typical oligopeptide permease (OPP) or ABC transport system complex. Therefore, sarA, oppB, oppC, oppD, and oppF could form a gene cluster whose products can assemble into an oligopeptide ABC-transport system complex.

To confirm whether the gene organization of these five genes is conserved as a gene cluster across different bacterial genomes within the F. prausnitzii taxonomic clade, we conducted a gene cluster search using cblaster (v1.3.20)⁹⁸ in 17,496 assemblies of F. prausnitzii from publicly available databases such as UHGG (v2.0.2)⁹⁹, SPIRE¹⁰⁰, and our in-house F. prausnitzii MAG and isolate genomes⁶. We found that all 5 genes co-occur in 2,185 F. prausnitzii genomes, situated close to each other (maximum distance between any 2 hits in a cluster <20 kb). The F. prausnitzii species with the complete OPP operon are Faecalibacterium butyricigenerans (223), Faecalibacterium duncaniae (3), Faecalibacterium longum (1,568), F. prausnitzii _I (372), Faecalibacterium sp900539945 (1), Faecalibacterium sp934514285 (15) and 3 unclassified genomes. Therefore, the SarA gene and the other four genes (oppB, oppC, oppD and oppF) do form a gene cluster or operon that is conserved in its gene organization within the F. prausnitzii genome. We also conducted 3D structure prediction using AlphaFold3 for the protein complex, compared it to the experimentally derived structures, and found that the predicted structure is similar to experimentally validated ABC transporter structures in M. tuberculosis (PDB: 8JA7)²⁴ and E. coli MetNI (PDB: 6CVL)²⁵. Therefore, the SarA with the OppA domain can function as OppA and can assemble with the other four subunits to form an OPP complex.

To further identify the potential ligand of this OPP complex, we conducted gene function annotation with Bakta (v1.11.0)⁶⁸. We found that one of the five genes in the OPP gene cluster was annotated as an ABC-type glutathione transport system ATPase component (gsiA). We then compared the protein sequences of SarA (OppA), OppB, OppC, OppD and OppF to 3,239 existing protein sequences of the ABC-type glutathione transport system using blastp (v2.16.0 + ). We found that all five proteins match existing protein sequences of the ABC-type glutathione transport system (E-value < 1.00 × 10⁻¹⁰). The OPP complex may therefore function as a glutathione transport system.

Identification of nicotine degradation genes

To identify the genes from reference genomes potentially related to nicotine degradation, we obtained sequences of 92 proteins involved in two nicotine-degradation pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) (KEGG entries: M00810 and M00811) and searched for them in the genes of 909 reference genomes using tblastn (v2.14.1-gompi-2023a)¹⁰¹. In total, we identified 2,434 genes similar to known nicotine-degradation genes (E-value < 1.00 × 10⁻¹⁰).

Gut microbial SNV-based age prediction model

We constructed a host age prediction model based on light gradient-boosting machine implemented in the R package lightgbm (v3.3.5). To ensure generalizability of the model, we used tidymodels (v1.1.0) to conduct nested cross-validation with a five-times-fivefold split for the outer loop and a fivefold split for the inner loop to assess the predictive power of gut microbial SNVs on host age. To avoid information leakage for the samples collected from individuals with close microbiome kinship, samples from the same lineage group were assigned to the same split dataset (either internal train dataset or internal test dataset) using the group_vfold_cv() function during the nested cross-validation procedure. We defined lineage groups using all SNVs used for age model training, pooled across species. Pairwise sample distances (scaled Manhattan distance) D were computed, and kinship was then defined as K = 1 – D. From D = 1 – K, we built a neighbour-joining tree and rooted it to orient the cut. Lineage groups were assigned by a single tree-height cut at h = (1 – c)/2 with kinship cut-off c = 0.9 (thus h = 0.05, equivalent to distance cut-off 0.1). Each descendant clade crossing this cut was then treated as one lineage group. In all cross-validation runs, we used a set of predetermined hyperparameters: trees = 2000, tree_depth = 15, min_n = 30, sample_size = 0.6, mtry = round(0.2 * total_feature_number), and learning_rate = 0.01. We estimated the importance of features using Shapley additive explanation, which assigns each feature an importance value for a particular prediction. After 5 × 5 rounds of fivefold cross-validation, we calculated the median importance value of each SNV feature, sorted the SNV features from high to low median importance values, cumulatively summed the median importance values until we reached 95% of the total sum of importance values, and took the counted SNV features as the predictive features for the final model. Because lineage groups were defined using SNVs pooled across species, this approach captures overall microbiome kinship but may not fully prevent information leakage for individual species whose species-specific kinship structure differs from the overall pattern.

Partial correlation analysis

The partial correlation test was conducted using the pcor() function of R package ppcor (v1.1). For the partial correlations between gut microbial SNV age and other biological age indices and between gut microbial SNV age and host intrinsic phenotypes, we used partial Spearman’s correlation, adjusting for chronological age. For the partial correlations between gut microbial SNV age and untargeted blood metabolome, we used a partial Pearson correlation, adjusting for chronological age, sex and blood creatinine level.

Ethical approval

The Lifelines study was approved by the Ethics Committee of the UMCG (METc 2007/152). All participants signed an informed consent form prior to enrolment. Additional written consents were signed by the DMP participants or legal representatives for children aged under 18 years. The LLD study was approved by the Institutional Ethics Review Board of the UMCG (ref. M12.113965), the Netherlands. The 300OB study was approved by the Institutional Review Board (IRB) CMO Regio Arnhem-Nijmegen (nr. 46846.091.13). The 500FG study was approved by the Ethical Committee of Radboud University Nijmegen (NL42561.091.12, 2012/550). The inclusion of volunteers and experiments was conducted according to the principles expressed in the Declaration of Helsinki. All volunteers provided written informed consent before any materials were collected. The 1000IBD study was approved by the IRB of the UMCG, Groningen (ref. M12.1139652008.338), the Netherlands. The 300TZFG study was approved by the Ethics Committee of the KCMC University (KURERC) (No. 936) and the National Institute for Medical Research (NIMR/HQ/R.8a/Vol. IX/2290) in Tanzania.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.