The Gut Microbiome Strongly Mediates the impact of Lifestyle combined variables on Cardiometabolic Phenotypes

Study cohort and sample acquisition

This study is based on the European MetaCardis cohort (http://www.metacardis.net/), initially designed to study the relationships between the different stages of the disease and the gut microbiome. The cross-sectional European MetaCardis study cohort covered a wide range of metabolic and cardiac phenotypes of 2214 participants that were recruited between 2013 and 2015 in Denmark, Germany and France. The study cohort, including protocol information, exclusion criteria, definition of the groups, biochemical analyses, collection of anthropometric and clinical data, were extensively described previously ^2,3,29,71.

Briefly, for the purpose of the present study we focused on a subset of 1,643 individuals, including healthy normal- weight individuals (BMI <25 kg/m²) as well as individuals with non-exclusive metabolic conditions such as type 2 diabetes, metabolic syndrome, obesity and cardiac disease, as summarized on the flow chart of Supplemental Figure 9 and described in the supplementary Table S1. The study protocol was approved by the Ethics Committee at the Medical Faculty at the University of Leipzig, Germany (application number: 047-13-28012013), the ethical committees of the Capital Region of Denmark (H-3-2013-145) and the ethics committee ‘Comité de Protection des Personnes’ (CPP) Ile-de-France III no.IDRCB2013-A00189-36 in France and was registered at https://clinicaltrials.gov/ (NCT02059538). The observational cohort design complied with all relevant ethical regulations, aligning with the Helsinki Declaration and in accordance with European privacy legislation. All participants provided written informed consent.

As a cohort of replication for people with overweight and obesity, we used data from the GutInside study as described in Table S2 ³⁰. Subject were recruited in France from September 2018 to January 2020 in different regions and the population was composed of 433 individuals with overweight or obesity (BMI ≥ 25 kg/m²). While there were involved in a dietary intervention program, we here only examined baseline data for the purpose of this work. Similar questionnaires were used to evaluate lifestyle aspects.

Clinical, anthropometric variables and definition of study groups

Groups were defined along international definitions of disease, with overweight and obesity assessed according to the WHO criteria (BMI ≥ 25 and BMI ≥ 30 kg/m², respectively). The metabolic syndrome (MS) was defined based on the International Diabetes Federation 2005 Consensus Worldwide Definition of the Metabolic Syndrome (http://www.idf.org/metabolic-syndrome), e.g. waist circumference> 94 cm in men and> 80 cm in women, and any two of the following four factors: i) triglycerides levels ≥ 1.7 mmol/L or treatment for lipid abnormality (statin and/or fibrate or ezetimibe); ii) high density lipoprotein (HDL) cholesterol< 1.03 mmol/L for European men and 1.29 mmol/L in European women and/or treatment for lipid abnormality; iii) abnormal blood pressure (BP): systolic BP ≥ 130 mmHg and/or diastolic BP ≥ 85 mmHg or patients taking anti-hypertensive treatment and, iv) abnormal fasting plasma glucose ≥ 5.6 mmol/L or prevalent T2D. T2D status was defined using the American Diabetes Association (ADA) definition: fasting glycemia> 6.9 mmol/l and/or 2h values during an oral glucose tolerance test> 11 mmol/l and/or hemoglobin A1c (HbA1c, glycated hemoglobin) ≥ 6.5% and/or use of any anti-diabetic treatment (American Diabetes Association 2018). Hypertension status was defined according to the American College of Cardiology and American Heart Association⁷². For obesity specifically, participants were recruited into two groups: group 2A, consisting of individuals with mostly grade II obesity (BMI≥ 35) none of whom had T2D or previous cardiovascular conditions, whereas group 2B consisted mostly of individuals with grade III obesity (BMI ≥40), who were eligible for bariatric surgery. T2D was not an exclusion criterion for this particular group (in contrast to group 2A) and patients from group 2B had overall more severe metabolic impairments than those in group 2A. Individuals with heart failure were defined according to the American College of Cardiology, American Heart Association and the Heart Failure Society of America ⁷³. We combined the group 2A and 2B in this study for descriptive purposes. Our group of participants in the present study with coronary heart disease included patients with first events of acute coronary heart diseases, chronic coronary artery disease with heart failure (defined as a left ventricular ejection fraction ≥45%) or without heart failure (<45%). We excluded participants suffering from heart failure related to an origin other than coronary artery disease.

Weight, height and waist circumference were assessed during the clinical inclusion visit according to standardized procedures using the same scales and units. Body fat mass and fat-free mass were measured via bioelectrical impedance analysis. Systolic and diastolic blood pressure were measured using a mercury sphygmomanometer (measures were taken three times on each arm, and the mean of the last 2 measurements in the right arm was used for analyses). Moreover, a detailed list of prescribed medications as well as the patient’s medical history were gathered. We have carefully cleaned and integrated each drug treatment according to its class of molecules, as described in this study ²⁹.

Dietary intake and lifestyle data assessment

Foods groups and nutrient content

The initial responses to the FFQ are frequency choices that indicate consumption units per month, week or day for each 160 food items. These frequency choices are then converted to daily frequencies and finally to food group intakes in g/d and to nutrient intakes in g/d, mg/d or mcg/d – based on information on portion size (i.e., grams contained in one portion) and nutrient content (i.e., nutrients contained in one portion). This information was taken from nationally relevant food composition databases.

Consumption of food groups was obtained by grouping the 160 food items into nutritional interest groups. The 22 initial groups of the MetaCardis study were broken down into 42 subgroups (shown in Supplementary Table S4). Dietary diversity scores were based on these 42 subgroups, according to ad hoc criteria. In sensitivity analysis, we constructed these same dietary diversity scores by removing “liquid” foods from the diet, i.e., hot drinks, fruit juices, sugary drinks, alcoholic drinks and milk because the quantities consumed can unbalance the diversity scores based on the uniformity of consumption (Simpson index).

Dietary “inflammatory potential” of the diet

The “inflammatory potential” of the diet was determined using the alternative Dietary Inflammatory Index (aDII) ^76,77, an alternative version of the original Dietary Inflammatory Index (DII) ⁷⁸, developed by Cavicchia et al. ⁷⁹ and following the methodology initially provided by van Woudenbergh et al. ⁸⁰. Briefly, the ADII included 34 food items, which were standardized after adjustment for energy intake using the residual method developed by Willett et al. ⁸¹ to reduce the variation in dietary intake resulting from differences in physical activity, body size, and metabolic efficiency ⁸². Each item was then multiplied by a proposed coefficient ⁷⁸, and items were summed to create the aDII score. Higher aDII values reflect a more proinflammatory diet.

Dietary quality scores

The Alternative Healthy Eating Index (AHEI) and Dietary Approaches to Stop Hypertension (DASH) scores were computed to evaluate overall dietary quality, two a priori scores, taking into account the existing correlations between the different components of the diet developed in order to overcome the limitations linked to approaches focusing only on food or nutrient groups ⁸³ and based on foods and nutrients predictive of chronic disease risk ^84,85. As a measure of healthy US-style eating, the AHEI-2010 assigns 0–10 points to each of the 11 dietary components based on the portion size ^83,84. The DASH score contains 8 components, each of which receives 1–5 points according to its consumption quintile ⁸⁶.

Dietary diversity and variety scores

Several dietary diversity and variety groups were constructed, depending on the groups retained in the method of construction (42 groups aggregated, with or without “liquid” foods, presented in Supplemental Table S4) and the method retained for the estimation (count-based or evenness-based estimation).

First, we used count measures of the diversity, whereby the number of consumed food items and food groups is recorded. Count measures focus on counting different food groups and subgroups; the distribution of consumed food quantities is not taken into account. Nevertheless, research studies have been limited by count-based measures and inconsistencies in the number and types of foods groups considered in these estimations ¹⁶. To overcome these limitations, we have systematically varied the number of foods groups in sensitivity analyses, in particular removing liquid foods from the diet. Secondly, we have used measures of the diversity in terms of number as well as distribution of different food items, using the Berry-Index, which is also known as the Simpson-Index⁸⁷, indexes also used to assess metagenomic or ecological diversity for each measure of the diversity considered in the count- based at the first step. This score increases if the food items are more equally consumed rather than being more concentrated.

We calculated the HFD score designed and published by Drescher et al. ⁴⁰. This index could be interpreted as a healthy food diversity indicator, i.e. the diversity of the healthy products quality (dietary diversity of “considered” health products) ⁴⁰. Finally, we also calculated the normalized composite dietary and variety score of the diet ^38,39, a score that combined diversity and variety (simple counts, normalized on a scale of 0 to 100).

All these scores account for different aspects of dietary diversity and variety. For example, some scores take into account the entire diet (identified by the word “diversity all diet” in Figure 1B), while others focus on specific groups of the diet, e.g., food items (identified by the word “variety” in Figure 1B). We used diversity and variety scores based on simple counts in consumption groups (identified by the word “count” in Figure 1C) and scores taking into account the evenness of consumption within groups (identified by the word “Simpson” in Figure 1C). Simple counts-based scores were calculated by assigning 0 points for each food group (diversity scores) or food considered in the group (variety scores) if the person did not consume the food or group in question and 1 point if she consumed it (consumption >0 grams). Thus, total diet diversity was a score ranging from 0 to 42 initially groups derived from the FFQ (see the initial 42 groups in Supplemental Table S4).

Simpson based scores were calculated using food groups (diversity) or foods in a given food group (variety) as different species, by analogy with (microbial) species used to calculate ecological indexes to measure the diversity of a given environment. Simpson indexes take into account the quantity of species present but also their abundance (equity).

R: richness: number of species considered, pi =relative abundance of the species (share of product i in the total amount of food consumed).

By subtracting this score from 1, we therefore obtain an index bounded between 0 and 1, 0 corresponding to the minimum diversity and 1 to the maximum (in equal shares of all products considered).

To overcome inconsistencies in the number and types of foods or food groups considered we have performed sensitivity analyses of these scores, varying the number of groups used in taking these calculations into account. In particular, when assessing total diet diversity, we created scores that did not take liquid food into account (See supplementary Table S4 for additional food groups excluding beverages and milk variables).

The QASD (Quality Activity Smoking Diversity) lifestyle score construction

A new lifestyle score was designed considering the quality of the diet, total physical activity, smoking status and dietary diversity. The score was calculated based on the tertile values of variables corresponding to: food quality (aHEI: 0 point for low tertile, 1 point for medium tertile and 2 points for high tertile), dietary diversity of the all diet, without beverages (Simpson Index) (0 point for low tertile, 1 point for medium tertile and 2 points for high tertile) and total physical activity in MET-h.week^-1 (0 point for low tertile, 1 point for medium tertile and 2 points for high tertile). For smoking status 0 point was given for smokers and 1 point for non-smokers. Taking together these variables gives a theoretical score between 0 and 7. We then grouped the modalities 0, 1 and 2 and 6 and 7 to correct for low numbers in some categories during the stratifications, giving from the lowest to the highest lifestyle score, 5 modalities: “0-1-2”, “3”, “4”, “5”, “6-7” for Figure 1,3 and Supplemental Figures 1-3 and 3 modalities “0-1-2-3”, “4”, “5-6-7” for Figures 2,4,5 and Supplemental Figures 4-8.

Extraction of fecal genomic DNA and gut microbiota sequencing

Gut microbiota sequencing in the MetaCardis cohort was performed as previously described^2,29,71. Briefly, participants collected fecal samples within 24 hours before enrollment. Samples were either stored immediately at −80°C or briefly conserved in home freezers, before transport to the laboratory where they were immediately frozen at −80°C following guidelines53. Total faecal DNA was extracted following the International Human Microbiome Standards (IHMS) guidelines (SOP 07 V2 H) and sequenced using ion-proton technology (ThermoFisher Scientific) resulting in 23.3 ± 4.0 million (mean ± SD) 150-bp single-end reads per sample on average. Reads were cleaned using Alien Trimmer (v0.2.4)⁸⁸ in order to remove resilient sequencing adapters and to trim low quality nucleotides at the 3’ side (quality and length cut-off of 20 and 45 bp, respectively). Cleaned reads were subsequently filtered from human and potential food contaminant DNA (using human genome RCh37-p10, Bos taurus and Arabidopsis thaliana with an identity score threshold of 97%). Filtered high-quality reads were mapped with an identity threshold of 95% to the 9.9-million-gene catalogue using Bowtie2 (v.2.3.4)⁸⁹ included in the METEOR v3.2 software⁹⁰ (https://forgemia.inra.fr/metagenopolis/meteor). A gene abundance table was generated by means of a two-step procedure using METEOR. First, the uniquely mapping reads (reads mapping to a single gene in the catalogue) were attributed to their corresponding genes. Second, shared reads (reads that mapped with the same alignment score to multiple genes) were attributed according to the ratio of their unique mapping counts. The gene abundance table was processed for rarefaction and normalization and further analysis using the R package MetaOMineR⁹¹. To decrease technical bias due to different sequencing depth and avoid any artefacts of sample size on low-abundance genes, read counts were rarefied. The gene abundance table was rarefied to 10 million reads per sample by random sampling of 10 million mapped reads without replacement. The resulting rarefied gene abundance table was normalized according to the FPKM (fragments per kilobase of transcript per million mapped reads) strategy (normalization by the gene size and the number of total mapped reads reported in frequency) to give the gene abundance profile table and binned by functional and phylogenetic categories as carried out within the MOCAT2 framework⁹². A total of 3463 metagenomic species (MGS; co-abundant gene groups with more than 50 genes corresponding to microbial species) were clustered from 1,267 human gut metagenomes used to construct the 9.9- million-gene catalogue (reference), as described previously⁹³. MGS abundances were estimated as the mean abundance of the 50 genes defining a robust centroid of the cluster (if more than 10% these genes gave positive signals). Microbial gene richness (gene count) was calculated by counting the number of genes that were detected at least once in a given sample, using the average number of genes counted in ten independent rarefaction experiments.

Quantitative metagenomic profiling of GutInside cohort

Gut microbiome sequencing of fecal samples from 433 individuals of the GutInside study (reference Biomedicines) was performed by GeneWiz (https://www.genewiz.com/) using Illumina HiSeq platform, resulting in 23.17±5.1 million read pairs (mean ± std. deviation) per sample. Raw reads were processed with NGLess⁹⁴ for quality trimming (minimum read quality=25; minimum read length=40), host contaminant removal vs. reference human genome (min. identity=90%, min. match size=45), alignment filtered reads over the 9.9-million gene catalogue (reference) (min. identity=95%, min. match size=45), and generation of gene abundance table with the dist1 metric of NGLess, equivalent to the counting strategy described above for the Metacardis cohort (first the genome abundances are computed from unique mapped reads and then are corrected by the multiple mapped reads weighted by the coverage of unique mapped reads). The gene abundance table was processed for rarefaction and RPKM normalization with MetaOMineR⁹¹ R package. The raw gene abundance table was rarefied to 10 million reads per sample and the rarefied gene abundance table was normalized according to the FPKM strategy, from which MGS abundances were computed as described above in Metacardis cohort. From NGLess filtered reads, genus-level abundance profiles in the mOTU space were generated with mOTU v2.6.1⁹⁵, from which enterotype classifications was performed after rarifying to 4000 counts per sample using the Dirichlet Multinomial Mixture (DMM) method⁵⁵ and implemented in the Dirichlet Multinomial R package (Dirichlet-Multinomial Mixture Model Machine Learning for Microbiome Data. R package version 1.28.0). The minimum Laplace metric was used as criteria to select the best stratification of the cohort⁵⁵, which was at 4 discrete microbiome compositions that corresponded to the Rum, Prev, Bact1 and Bact2 enterotypes described in the Metacardis cohort⁴¹.

Analyses of microbial diversity and community structure

The estimation of the explanatory power of nutritional features regarding relative microbiome profiles derived from MGS abundance data was performed using univariate and multi-variate stepwise distance-based redundancy analyses (dbRDA) as implemented in the R package vegan (Community Ecology Package. R package v.2.2–1)⁹⁶. Redundance filtering of collinear lifestyle variables were done previously to these analyses with the R package FMradio⁹⁷; threshold=0.9 on Spearman pairwise correlation matrix). Microbiome inter-individual variation was visualized by principal coordinates analysis using Bray–Curtis dissimilarity on the relative abundance matrix derived from MGS abundance data. Environmental fitting of nutritional covariates with significant impact on microbiome composition based on dbRDA analyses over PCoA ordination from Bray-Curtis inter-sample dissimilarity matrix was computed with envfit function of vegan R package. Associations of nutritional variables with enterotype status were tested with logistic regression analyses adjusted by age, gender, center, BMI, and the intake of metformin, statin and PPI.

Metabolomic profiling

We used untargeted UPLC–MS data generated by Metabolon, as described in previous study⁷¹. Briefly, serum samples were extracted and profiled by Metabolon using a UPLC–MS-based methodology98. Annotated metabolites and unknown features (denoted X-00000) were identified by comparing sample features with ion features in a reference database of pure chemical standards and previously detected unknowns, followed by detailed visual inspection and quality control, as reported99. For all metabolomic assays, we randomized the sample preparation order across the whole study so that each sample preparation batch included samples from all study groups. For MS untargeted assays, median batch correction was performed by adjusting batch-wise study sample variable medians according to a scalar derived from adjusting pooled reference sample medians, so that pooled reference sample medians are identical across all batches. The randomized sample preparation batches were also tested for association with study groups using univariate statistics (Fisher’s exact test or Kruskal–Wallis test), and P > 0.05 was observed across all methods (GC–MS, Fisher’s exact test, P = 0.23; UPLC–MS targeted, Fisher’s exact test, P = 0.12; and UPLC–MS untargeted (Metabolon), Fisher’s exact test, P = 0.65). In addition, NMR run order exhibited a Kruskal– Wallis P = 0.49.

Drug deconfounding analysis

The pipeline was used to assess to what extent observed differences among pairwise levels of the QASD score decomposed into 3 levels (high vs. medium, high vs. low, medium vs. low; Low=QASD 1,2,3; Medium = 4; High=5,6,7) between study participants in microbiome, metabolome and bioclinical feature abundance are confounded, in the sense of being consequences of other (treatment or risk factor) variables different among the groups more so than characteristic of the specific phenotype itself. We employed the post hoc filtering approach implemented in the R package metadeconfoundR (version 0.1.8; see https://github.com/TillBirkner/metadeconfoundR or https://doi.org/10.5281/zenodo.4721078) that was devised within the MetaCardis consortium^29,71. The pipeline has two steps. In the first step, all associations between -omics features and the set of independent variables (disease status, drug treatment status and risk markers, including age and smoking status) are determined under non-parametric statistics (Mann–Whitney U-test (MWU) or Spearman test, adjusted for multiple testing using the Benjamini–Hochberg method). In the second step, for each feature significantly (FDR < 0.1) associated with defined phenotype status, it is checked whether it has significant associations with any potential confounder. If not, it is considered trivially unconfounded (not confounded (NC)). If at least one covariate also has significant association with the feature, then, for each such covariate, a post hoc test for confounding is applied. This test takes the form of nested linear model comparisons (likelihood ratio test for P values), where the dependent variable is the feature (X), and the independent variables are the QASD group (A) and the tested covariate (B) versus a model containing only the covariate (B), thus testing whether QASD status (A) adds explanatory statistical power beyond the covariate (B). If this holds (likelihood ratio test (LRT) P < 0.05) and the confidence interval for predictive inference in the LRT test does not include zero for all covariates (B), then disease status is strictly deconfounded (SD) concerning its effect on feature X; it cannot be reduced to any confounding factor. For each covariate (B) where significance is lost, a complementary modeling test is performed of the complementary model pairs, predicting (X) as a function of (A) and (B) versus a model containing (A) alone, thus testing whether the covariate (B) in turn is equally reducible to (A). If for at least one such covariate (B), (B) has independent effect (LRT P < 0.05) on top of (A), then the feature (X) is considered confounded by (B). However, if in none of the pairwise tests the original significance holds, then (A) and (B) are considered so correlated that their relative influence cannot be disentangled. We consider these cases ambiguously deconfounded (AD), in the sense that, for these cases, clear confounding influence can neither be concluded nor ruled out. The R package was applied to the present dataset considering medication status (statin, metformin and PPI intake, all as binary variables) and Metacardis clinical group as covariates and with country of recruitmet as random effect, and was runned on metagenomic (MGS, mOTU, Gut Metabolic Modules (GMM)) and metabolomic features (serum, urine) presents in at least 20% of the samples. For metagenomic features, abundances were corrected for bacterial cell count by multiplying by an index factor calculated as the bacterial cell count of the sample divided by the mean value of this bacterial cell count over the dataset as a whole ⁷¹. Our deconfounding pipeline considers linear effects related to drug categories. Still, we were not able to control for every possibly lifestyle confounding factor, making a lack of full confounding adjustment a limitation of our study.

Statistical analyses and machine learning models

Data-management and statistical analysis were conducted using SAS version 9.3 and R (R statistical software, Vienna, Austria). The level of significance used was P<0.05, FDR<0.01, P for interaction<0.1.

Throughout the analyses, several adjustment strategies were used. First, for the analyzes concerning GMGR, a normally distributed variable, we systematically carried out and tested sensitivity analyses, taking into account the following confounding factors: age, sex, recruitment center, taking antibiotics, metformin and statins, BMI, which all are potential confounders identified in previous analyses from the MetaCardis consortium ^2,3,29,71. Then, in accordance with what is done for dietary intake, models were adjusted on energy in kcal and aHEI score when we wanted to isolate the independent effects of nutritional quality. Finally, for specific omics features, we used the analysis pipeline published and developed in the MetaCardis study, thus allowing a better comparison between studies in the literature (see previous Drug deconfounding analysis section).

Predomics R package version 1.01¹⁰⁰ was used to build interpretable predictive models of QASD status (high vs. medium, high vs. low, medium vs. low; High=QASD 1,2,3; Medium = 4; High= 5,6,7) based on MGS and metabolomic markers retrieved in the drug deconfounding analyses. Models were trained on three different algorithms including Random Forest as state-of-the-art (SOTA) method, and Binary and Ternary native predomics models describing simple ecological relationships in microbial ecosystems that were learned with Terbeam heuristics. Details of the different heuristics and BTR models can be found in Prifti et-al¹⁰⁰. Models were trained on 80% of the data and evaluated for accuracy and AUC on a 10 times 10-fold cross-validation schema, and the results of the best algorithm were further explored to extract a family of best models (FBM), described as models whose accuracy is within a given window of the best model’s accuracy. This window is defined by computing a significance threshold assuming that accuracy follows a binomial distribution (p<0.05). No hyperparameter optimization was performed. Features included in the FBM were further explored in terms of prevalence across models and feature importance, described as the mean decrease accuracy (MDA) of the model after feature removal. The best model in each prediction task was evaluated on the 20% of the data excluded during the training process (hold-out data).

Mediation analysis

Mediation analyses was used to infer causal mediation effects between MGS, serum metabolites and pairwise levels of the QASD score (High vs. Low, High vs. Medium, Medium vs. Low; Low=QASD 1,2,3; Medium = 4; High=5,6,7). To reduce the testing numbers, we focused these analyses on the MGS and metabolites with OK_sd status in the drug deconfounding pipeline and showing an absolute effect size greater than 0.1 on Cliff’s Delta between pairwise levels of the QASD score (features included in the Figures 3,4). In a first step, pairwise associations between metabolites and MGS were individually tested with linear regression analyses adjusted by country of origin, age, gender and intake of metformin, statin and protein pump inhibitors, and only significant pairs were retained for mediation analyses (FDR<0.05). MGS abundances were transformed with empirical normal quantile transformation method to ensure normality before the analyses. In the second step, similar regression framework with same adjustments was adopted to identify significant mediations under two possible hypothesis in terms of direction of the mediation using the mediate function of mediation V4.5.0 R package¹⁰¹. The first hypothesis delved into the role of MGS as mediators in shaping the relationship between the QASD score and the host’s metabolomic profile (Direction 1 in Figure 6A diagram). The second hypothesis aimed to explore whether the interplay between the QASD score and the host’s metabolome (mediator) could, in turn, influence the metagenomic landscape (Direction 2 in Figure 6A diagram). Significant mediations were considered those with FDR<0.05 for Average Causal Mediation Effect (ACME), Average Direct Effect (ADE) and Total effect. Similar analyses were carried out to infer causal mediation effects between levels of the QASD score, the microbial gene richness and the clinical status of the individuals in terms of corpulence (BMI) and metabolic status (glycated hemoglobin, HOMA-IR).