Fecal microbiome predicts treatment response after the initiation of semaglutide or empagliflozin uptake

Study design and sample collection

Twenty subjects (10 women and 10 men) were enrolled in the study between November 2019 and January 2023 from the University of Tartu Hospital, Tartu, Estonia. The study followed approval from the Research Ethics Committee of the University of Tartu (290/T-20), and informed consent was obtained from each subject. All participants had been diagnosed with type 2 diabetes, had a BMI of 32 or more, glycohemoglobin (HbA_1c) <86 mmol/mol (<10%), and were taking 1.5 grams of metformin per day (Figure 1A). The main exclusion criteria included the following:

1. having positive glutamic acid decarboxylase (GAD 65) antibodies,

2. an earlier prescription of oral or intravenous antibacterial medications within the last two months,

3. a change in antidiabetic treatment regimen within the last 90 days,

4. previous use of GLP-1 receptor agonists or SGLT-2 inhibitors,

5. use of spironolactone within the last 60 days,

6. current usage of immunomodulating medications (glucocorticoids, cytostatic drugs, or biopharmaceuticals),

7. being pregnant or planning pregnancy,

8. consuming oral contraceptives or oral hormone replacement therapy medications,

9. having previously been diagnosed with a malignant tumor,

10. having a diagnosis of New York Heart Association (NYHA) class III and IV heart failure or a diagnosis of severe liver disease.

• Download figure
• Open in new tab

Figure 1.

Study design and clinical endpoints. A - overview of the study design. B - Changes in HbA_1c and BMI after 3 and 12 months following semaglutide (GLP-1-RA) and empagliflozin (SGLT-2i) initiation. After inclusion and exclusion criteria, 20 T2D patients were recruited to receive either semaglutide or empagliflozin. Subjects provided 4 samples: before treatment initiation and 1, 3, and 12 months after treatment initiation. T2D - type 2 diabetes; HRT - hormone replacement therapy; GAD65 - glutamic acid decarboxylase 65 antibodies, qd - daily (quaque die).

As a result, the study group taking semaglutide included 6 women and 4 men with ages ranging from 48 to 77 (mean 65.7 ± 8.4) and BMI ranging from 35.69 to 47.87 kg/m² (mean 42.2 ± 5.27). The group prescribed empagliflozin included 4 women and 6 men with age ranging from 52 to 75 (mean 61.7 ± 7.7) and with BMI ranging from 31.67 to 51.04 kg/m² (mean 37.02 ± 5.89) (Supplementary Table 1).

The study workflow is shown in Figure 1. Each participant donated a venous blood sample, an overnight urine sample, and a fecal sample for microbiome analysis at baseline before treatment initiation, after 1 month, 3 months, and 12 months of treatment. This means that each participant served as their own control, allowing us to minimize inter-individual variability. In addition, body anthropometric measures were recorded. The dose of empagliflozin was 10 mg once a day. A once-weekly subcutaneous (s/c) semaglutide was started at 0.25 mg dose and then up-titrated to 0.5 mg once a week at week 5 and 1 mg once a week at week 9. Treatment efficacy was assessed by measuring glycohemoglobin (HbA_1c), fasting blood glucose, body mass index (BMI), waist circumference, and inflammation markers (neutrophils, lymphocytes, white blood cells, and neutrophil-to-lymphocyte ratio at Month 3 and Month 12. Several other parameters, like glomerular filtration rate (eGFR), blood lipids, and transaminases, were also analyzed.

DNA extraction and sequencing

For DNA extraction, study participants collected fecal samples collected with a special stool sampling kit with instructions. Samples were collected using an RNALater sample stabilizing solution and taken to the University of Tartu Hospital, Tartu, Estonia. Immediately after arrival, samples were stored at -20°C and then taken to the Institute of Genomics, University of Tartu, Estonia, where samples were stored at -80°C before follow-up analyses.

DNA extraction of fecal samples was carried out with Qiagen DNeasy PowerSoil^® Pro extraction kit (Qiagen, Venlo, Netherlands). Samples were first brought out to room temperature at 25°C for complete dethawing, and 100 µl of starting material was taken from samples. In the following step, 100 µl of Qiagen PowerBead Pro Tube silica gel beads were added to the sample material and heated at 37°C to ensure complete dissolution of salt crystals in RNALater. Next, 800 µl of solution CD1 was transferred to samples following incubation at 65°C for 10 minutes to activate lysis processes of microbial and human intestine cells. The sample disruption step was completed using a tissue homogenizer Precellys^® 24 (Bertin Instruments, Montigny-le-Bretonneux, France) at 2500 rpm for 2 x 30 seconds with a 30-second pause between two cycles. The following extraction process was conducted according to the manufacturer’s protocols. All extracted DNA samples were stored at -20°C.

The sequencing process of genomic DNA was conducted twice at the Core Facility of Genomics, Institute of Genomics, University of Tartu, Tartu, Estonia, first in February 2022 (N = 52 samples) and then in January 2023 (N = 26 samples). The concentrations and quality of extracted DNA samples were first assessed using a Qubit^® 2.0 Fluorometer (Invitrogen, Grand Island, USA). Extracted genomic DNA was then amplified in the PCR amplification reaction with forward and reverse primers designed for the approximately 460 bp long V3-V4 hypervariable region of the prokaryotic 16S ribosomal RNA gene – 16S_F (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG-3′) and 16S_R (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAAT CC -3)²¹.

Microbial community analysis

Raw sequences were demultiplexed with Illumina bcl2fastq2 Conversion Software v2.20. Bioinformatic analyses of genomic DNA reads were conducted with the open-source software platform QIIME2 (version 2021.11)²². For importing data into QIIME2, the script q2-tools was used in the format of PairedEndFastqManifestPhred33. The first analysis involved denoising sequences with DADA2 software as part of the QIIME2 platform²³. In addition, forward and reverse sequences were joined at this step. As a result of the first quality control, each base pair position was given a quality score according to whether reads were trimmed and truncated in the following analysis. Due to the better overall quality of reads, the remaining analysis was conducted using forward reads trimmed and truncated at positions 15 and 288. Reads marked as chimeric were removed with the consensus filter. Following the quality control, reads were aligned to amplicon sequence variants (ASV) with software MAFFT²⁴. The taxonomic classification was carried out using the QIIME2 script q2-feature-classifier with pre-trained naïve Bayes classifier against the SILVA 16S V3-V4 v132_99 with a similarity threshold set to 99%²⁵.

Statistical analysis

Statistical analyses of microbial ASVs were conducted in R (version 4.3.1). Packages used in the analyses included mia (v.1.8.0), dplyr (v.1.1.3), tidyr (v.1.3.0), and purrr (v.1.0.2). Figures were obtained with the package miaViz (v.1.8.0), pheatmap (v.1.0.12), and ggplot2 (v.3.4.3). Additional customizations were done using vector graphics software Inkscape (v.1.3).

Initially, 78 samples from 21 patients were included in the microbiome analysis, but four samples from one individual were removed due to the missing baseline sample. Additionally, one sample was excluded due to the low number of reads. Thus, the final number of samples was 73 samples from 20 patients.

Alpha diversity was assessed using the mia package from microbial genus-level data. The alpha diversity indices included the number of observed genera, the Shannon diversity index, and Pielou’s evenness index. Beta diversity was assessed using five principal components (PC1-PC5) of the principal component analysis on the Aitchison distance. In addition to the diversity analysis, the abundance of genera present in at least 10 samples per study group was analyzed. As a result, we investigated changes in CLR-values of 128 genera in the semaglutide study group and 147 genera in the empagliflozin study group.

To study the effect of two diabetes medications on the gut microbiome, we analyzed the changes in the microbial diversity and the abundance of the prevalent genera after treatment initiation. Changes in the microbiome were assessed with repeated measures ANOVA. Benjamini-Hochberg false discovery rate (FDR) method was used to account for multiple testing²⁶. In addition, two-sided dependent Welch t-tests were performed as post-hoc tests to study changes in the abundance of genera between two timepoints. Post-hoc tests were carried out for genera that showed a difference in the abundance of genera in the repeated measures ANOVA analysis.

To understand whether baseline gut microbiome could predict the efficacy of semaglutide and empagliflozin in antidiabetic treatment, baseline microbiome values were correlated with the change of clinical markers by the third treatment month. This is the timepoint where we detected the largest treatment effect of semaglutide and empagliflozin therapy. The baseline gut microbiome values included alpha and beta diversity indices and CLR-transformed values of the prevalent genera. Pearson correlation coefficient was calculated, and the Benjamini-Hochberg procedure for adjusting the p-values was used to account for multiple testing (FDR < 0.05).

Cohort overview and study design

In this study, we recruited 20 subjects with T2D who had been previously using metformin as a part of the treatment regime and had a BMI of at least 32. As the study was conducted in a routine clinical setting, the patients who had more severe obesity and/or diabetes compensation were more likely to be prescribed semaglutide, and subjects with heart and kidney failure or T2D complications were more likely prescribed empagliflozin as a supplementary treatment. Participants were asked to donate microbiome fecal samples four times – at baseline, 1 month (M1), 3 months (M3), and 12 months (M12). Fasting blood samples, overnight urine samples as well as clinical anthropometric measurements were collected at baseline, 3 months (M3) and 12 months (M12) (Figure 1A). The microbiome was profiled from fecal samples by sequencing the V3-V4 regions of the 16S rRNA gene. In total, 73 fecal samples were analyzed, with 10 participants taking semaglutide (35 fecal samples) and 10 participants taking empagliflozin (38 fecal samples).

In the semaglutide study group, mean glycohemoglobin (HbA_1c) values decreased by -11.37 mmol/mol (1.11%) between treatment initiation and the third month (Figure 1B, Supplementary Table 1). The decrease for empagliflozin users was -13.88 mmol/mol (0.88%). At M12, a decrease in HbA_1c constituted -8.2 mmol/mol (1.4%) and -18.25 mmol/mol (0.48%) for semaglutide and empagliflozin users, respectively.

Compared to baseline, BMI (kg/m²) for the semaglutide users decreased by 1.69 units by the third treatment month. The BMI in the empagliflozin study group decreased by 0.88 units by M3 (Figure 1C, Supplementary Table 1). At the 12-month timepoint, the BMI on average decreased by 5.61 units in the semaglutide group and increased by 0.7 units in the empagliflozin group.

Drug initiation shows no direct effects on the changes in the fecal microbiome

To elucidate the effect of antidiabetic drug usage on the gut microbiome, we compared the baseline gut microbiome with the microbiome composition of M1, M3, and M12 (Figure 1A) (Methods). First, we analyzed changes in the alpha diversity (observed richness, Shannon index, Pielou index) using repeated measures analysis of variance. No significant changes were detected before and after correction for multiple testing (Figure 2A, B; Supplementary Table 2). Second, we studied shifts in the overall composition (beta diversity) by analyzing changes in the first 5 principal components of the genus-level CLR-abundance profile. We detected nominally significant changes in PC1 in the semaglutide group and PC5 in the empagliflozin group (repeated measures ANOVA, p <= 0.05). Post-hoc tests indicated a nominally significant alteration in PC2 between baseline and M3 for empagliflozin users. However, similarly to alpha diversity, no changes remained statistically significant after applying FDR correction (Figure 2B, Supplementary Table 2).

• Download figure
• Open in new tab

Figure 2.

Effect of semaglutide and empagliflozin on the gut microbial diversity in T2D patients. Microbial alpha diversity (observed genera, Shannon diversity, Pielou’s evenness) and three principal components (beta diversity) at four timepoints. Pairwise t-test p-values comparing the microbiome diversity between baseline and later timepoints are shown. A - semaglutide (GLP-1-RA) study group. B - empagliflozin (SGLT-2i) study group. BL - baseline, M1 - 1^st month of treatment, M3 - 3^rd month, M12 - 12^th month.

Next, to characterize the changes in the microbiome in greater detail, we analyzed alterations in the abundance of bacterial genera prevalent in at least 10 samples per study group (128 genera for semaglutide, 147 genera for empagliflozin). First, we aimed to identify genera with altered abundance in between any of the four timepoints. In the semaglutide group, we identified the abundance to be changed between six genera (repeated measures ANOVA, p <= 0.05). These included genera Actinomyces, Candidatus Soleaferrea, Fusobacterium, Haemophilus, Pseudobutyrivibrio, and Veillonella (Figure 3A, Supplementary Table 2). However, none of these associations remained significant after applying the FDR correction. Similarly, in the empagliflozin group, the abundance of nine bacterial genera was identified to be different between the timepoints (repeated measures ANOVA, p <= 0.05). These genera included Agathobacter, Akkermansia, Coprococcus 2, DTU089, Parasutterella, Romboutsia, Ruminococcaceae UCG-014, Turicibacter, and uncultured_bacterium 8 from the order Mollicutes RF39 (Figure 3B, Supplementary Table 2). However, no significant associations were detected after p-value correction.

• Download figure
• Open in new tab

Figure 3.

Changes in the CLR-transformed abundances of bacterial genera in the semaglutide and empagliflozin study group. A - semaglutide study group (GLP-1-RA), B - empagliflozin study group (SGLT-2i). Genera detected as nominally significant in the repeated measures ANOVA analysis are shown. The p-values on the plot refer to post-hoc analysis using pairwise t-tests to compare the changes between the baseline timepoint and later timepoints. BL - baseline, M1 - 1st month of treatment, M3 - 3rd month, M12 - 12th month.

Next, for the nominal associations identified, we carried out further analyses to identify the timepoints between which the abundance of the genera was changed (post-hoc pairwise t-tests). Interestingly, we observed that in both study groups, the changes in the abundances of the bacteria were not evident when comparing the baseline with the latter timepoint. Instead, the change was more evident between the microbiome collected at Month 3 and Month 12 (Figure 3, Supplementary Table 2). Thus, the results suggest that semaglutide and empagliflozin do not directly affect the gut microbiome. The changes in the microbiome observed between the latter timepoints are likely due to changes in metabolic health or lifestyle.

Gut microbiome correlates with treatment efficacy

Next, we evaluated the ability of the baseline gut microbiome to predict changes in the major clinical endpoints, indicating the efficacy of the antidiabetic therapy. We calculated the Pearson correlation between the CLR-transformed abundance of the prevalent genera, alpha and beta-diversity metrics, and changes in clinical markers between BL and M3. As for clinical parameters, we focused on BMI, glycohemoglobin (HbA_1c), measures of lymphocytes and neutrophils, neutrophil-to-lymphocyte ratio (NLR), and amount of white blood cells (WBC). In the empagliflozin study group, glomerular filtration rate (GFR) was additionally included to assess how the baseline microbiome influences the drug’s effect on kidney function. In the semaglutide study group (GLP-1-RA), 55 out of 816 associations were nominally significant (p < 0.05), with the most association with the change in white blood cells and neutrophil-to-lymphocyte ratio (NLR) (Supplementary Table 2). After the FDR correction, the correlation between the abundance of Alistipes and neutrophil and NLR change remained significant (Supplementary Table 2). A higher abundance of Alistipes was associated with a higher decrease in neutrophils and NLR. Most importantly, six genera were found to be associated with the change in HbA_1c (Figure 4A, Supplementary Table 2). In the empagliflozin group, 48 associations among 1085 were detected as nominally significant (p < 0.05) (Figure 4B, Supplementary Table 2). However, none of these correlations were significant after correcting for multiple testing. Similarly to the semaglutide study group, the associations between six bacterial genera and the change of Hb_A1c were found to be nominally significant, including genera Prevotellaceae NK3B31 group, Candidatus Soleaferrea, and Escherichia-Shigella (Figure 4B, Supplementary Table 2). For example, a higher abundance of Escherichia-Shigella indicates a higher reduction in HbA_1c. Therefore, our results suggest that the fecal microbiome could be used to predict the treatment outcome of semaglutide and empagliflozin.

• Download figure
• Open in new tab

Figure 4.

Three nominally most significant associations between the baseline microbiome and change in glycohemoglobin (HbA_1c) by the third treatment month. A - semaglutide study group (GLP-1-RA), B - empagliflozin study group (SGLT-2i). Associations were estimated using Pearson correlation coefficients.