
==== Front
Sci Rep
Sci Rep
Scientific Reports
2045-2322
Nature Publishing Group UK London

72467
10.1038/s41598-024-72467-z
Article
Metformin increases gut multidrug resistance genes in type 2 diabetes, potentially linked to Escherichia coli
Kim Han-Bin 1
Cho Yong-Joon 23
Choi Sun Shim schoi@kangwon.ac.kr

1
1 https://ror.org/01mh5ph17 grid.412010.6 0000 0001 0707 9039 Division of Biomedical Convergence, College of Biomedical Science, Institute of Bioscience & Biotechnology, Kangwon National University, Chuncheon, 24341 Republic of Korea
2 https://ror.org/01mh5ph17 grid.412010.6 0000 0001 0707 9039 Department of Molecular Bioscience, Kangwon National University, Chuncheon, 24341 Republic of Korea
3 https://ror.org/01mh5ph17 grid.412010.6 0000 0001 0707 9039 Multidimensional Genomics Research Center, Kangwon National University, Chuncheon, 24341 Republic of Korea
14 9 2024
14 9 2024
2024
14 2148028 3 2024
9 9 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Metformin is the most commonly prescribed medication for treating type 2 diabetes (T2D). It is known that metformin can alter the gut microbiome, which influences the effectiveness of metformin treatment. We posited that if the gut microbiome, a reservoir of the resistome, is altered, then the resistome should change as well. To test this hypothesis, we reanalyzed microbiome data generated by Wu et al. (Nat Med 23(7):850–858, 2017), identifying antibiotic resistance genes (ARGs) and bacterial species. Through read-based analysis, we observed that the abundance of ARGs indeed changed in many samples treated with metformin. Moreover, the altered pattern was sufficiently heterogeneous across individual samples to allow subcategorization. We also found a strong correlation between the abundance of multidrug-resistant ARGs (MDR-ARGs) and the presence of E. coli. The contig-based analysis led to the same conclusion: an increase in MDR-ARGs due to metformin was associated with an increase in E. coli. In relation to this, we were able to confirm that the majority of MDR-ARGs are likely to originate from E. coli. These results suggest that metformin may have the potential side effect of increasing E. coli carrying ARGs, particularly MDR-ARGs, which could be a concern in T2D therapy that relies on metformin.

Keywords

Microbiome
Type 2 diabetes
Metformin
Antibiotic resistance genes
Resistome
Subject terms

Computational biology and bioinformatics
Microbiology
http://dx.doi.org/10.13039/501100003725 National Research Foundation of Korea RS-2024-00341909 RS-2023-00260267 Cho Yong-Joon http://dx.doi.org/10.13039/501100003653 Korea National Institute of Health 2024ER210600 Cho Yong-Joon issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Metformin (1,1-dimethylbiguanide hydrochloride) is the most commonly prescribed medication for treating type 2 diabetes (T2D), remaining one of the essential medicines for reducing blood glucose levels since its development in the early 1950s. Between 2000 and 2016, metformin was prescribed to 55–77% of T2D patients in the UK and the US1. Metformin promotes glucose uptake in peripheral tissues and reduces glucose release from the liver, thereby improving body tissue sensitivity to insulin. It is reported that metformin activates a transcriptional regulatory pathway, primarily involving two enzymes, liver kinase B1 (LKB1) and adenosine monophosphate-activated protein kinase (AMPK), which ultimately inhibits the synthesis of gluconeogenic enzymes, although its exact mechanism of action remains unknown2. Beyond T2D, its potential for cardiovascular protection, anti-inflammatory effects, and cancer treatment is being explored.

Recent studies have shown that metformin can selectively influence certain groups of bacteria within the gut microbiota in humans and mice3–5. Its effects on the compositional changes in the microbiome have been somewhat inconsistent across studies due to the high complexities, large individual differences in the human gut microbiome, and heterogeneous experimental designs2. However, there is consistency in the way that metformin can shift gut microbiomes. For example, Wu et al.6 showed that in a double-blind study design, groups treated with metformin had a significantly different gut microbiome composition compared to the placebo group not treated with metformin. In the fecal microbiome transfer experiment from human donors to germ-free mice, they found that fecal samples from four months of metformin treatment groups significantly improved glucose tolerance compared to those from donors before treatment. Additionally, some studies have shown that metformin is specifically associated with an increase in the abundance of several bacteria linked to positive metabolic effects, such as mucin-degrading A. muciniphila, and bacteria producing short-chain fatty acids (SCFA). This suggests that the beneficial effects of metformin on glucose metabolism may be associated with changes in gut microbiome2,7.

By contrast, a few other studies also attempted to explain the gastrointestinal side effects of metformin with changes in the gut microbiome. For instance, Maier et al.8 reported that approximately 24% of non-antibiotic drugs, including metformin, can inhibit the growth of at least one bacterial strain in vitro, suggesting that these drugs may have antibiotic-like side effects in humans. This study also indicated that antibiotic resistance mechanisms are shared among various drugs, by showing that antibiotic-sensitive strains, such as E. coli lacking an antibiotic resistance gene (ARG), are more susceptible to human-targeted drugs, while antibiotic-resistant strains, such as B. uniformis (HM-715), exhibit greater resistance to these drugs. Similarly, it has been reported that the abundance of ARGs, such as TolC and mdtC, increases with metformin intake in several disease cohorts6,9,10. However, few studies have investigated in detail which types of ARGs are increased by which types of bacteria. To explore these questions, we decided to reanalyze the microbiome data generated by Wu et al.6, which demonstrated changes in gut microbiomes by metformin. In the present study, using the same dataset, we tested whether metformin treatment in T2D patients can increase antibiotic resistance as a side effect. We identified ARGs belonging to different drug classes (DCs) and their proportions, analyzing the compositions at different time points after metformin intake, in comparison with placebo groups treated for the same periods.

Results

Obtaining metagenomic data to analyze the effect of metformin on the resistome in T2D

To investigate the impact of metformin on gut antibiotic resistome trends in patients with T2D, we obtained metagenomic data produced by whole genome shotgun sequencing (WGS) from the NCBI SRA, representing 40 patients with T2D. This cohort included a group that received metformin (M group, n = 22) and a placebo group (P group, n = 18), with samples collected before treatment (M0, P0), after two months (M2, P2), and after four months of treatment (M4, P4). A total of 102 ARGs were estimated from the WGS data derived from each sample after data cleaning (Supplementary Fig. S1). Samples missing any time point (n = 1 in M group, n = 1 in P group) were excluded from further analysis. We then estimated the values of ARG proportion (ARGP) for the remaining samples. Using these ARGP values (Supplementary Table S1), we performed a principal coordinates analysis (PCoA) combined with a permutational analysis of variance (PERMANOVA) test for the Canberra distance11 (Supplementary Fig. S2). As a result, we found that the ARG patterns in pre- and post-metformin groups were distinguishable (PERMANOVA P = 0.029 when performed on all groups, PERMANOVA P = 0.01 when performed on M2 and P2).

Notably, some samples at the initial time points (M0, P0) were far from the center in the PCoA plot (Supplementary Fig. S2). To determine how these samples differed from the others, we identified the DCs to which each ARG in the samples conferred resistance, and then calculated the total ARGPs in DC (ARGP-DC). Consequently, these five samples (42V1, 49V1, 23V1, 5V1, 35V1) contained an extremely high proportion (5–15%) of MDR-ARGs at the initial time points (M0, P0) (Supplementary Fig. S3). ARGs associated with peptide and phosphonic-acid DC were also found in greater quantities in these five samples.

In summary, certain samples exhibited distinct characteristics, with exceptionally high proportions of MDR-ARGs and other ARGs at the initial time point, unlike the majority of other samples. Hence, we decided to treat these outliers separately from the total of 33 filtered samples (n = 18 for M group; n = 15 for P group) for further analysis (Supplementary Fig. S1).

ARG proportions increased by metformin treatment

From the 33 filtered samples after preprocessing and cleaning the data, we identified 102 ARGs, which were then categorized into 11 DCs (Supplementary Table S2). We calculated the mean ARGP-DC (i.e., MARGP-DC). Our results indicated that MDR-ARGs significantly increased in MARGP-DC at M2 and M4 compared to M0 (Fig. 1a, Supplementary Table S3). A similar trend, albeit to a lesser extent, was observed in the placebo groups (P0, P2, P4) (Fig. 1a, Supplementary Table S4). This trend remained consistent when examined in terms of absolute abundance; the total read count of MDR-ARGs, normalized to FPKM (fragments per kilobase of transcript per million mapped reads), was notably higher at M2 and M4 than at M0, whereas no significant increase was observed in the P group (Fig. 1b). The ARGs associated with peptide DC (named peptide-ARGs) also showed a significant increase following metformin treatment.Fig. 1 Analysis of changes in abundance patterns in ARGs at different time points of metformin intake. (a) Stacked bar plot showing the percentage of MARGP-DC (i.e., mean of total ARG proportions in drug class). The X-axis represents each time point (M0: no metformin, M2: 2 months after metformin intake, M4: 4 months after metformin intake), and the Y-axis represents the MARGP-DC (%) value for each DC (i.e., drug class). Each color represents a specific DC as labeled at the bottom. (b) Paired bar plot showing the total normalized read count (i.e., FPKM; fragments per kilobase of transcript per million mapped reads) for ARGs associated with MDR and peptide DC at different time points after metformin intake. '**' means P < 0.01 according to the Nemenyi test. (c) Heatmap of different abundance measures of ARGs associated with different DC. ‘+’ in the cell indicates an increase in SPAD (i.e., sample proportions of ARG detection) values more than 30% compared to M0. '*' means P < 0.05, and ‘**’ means P < 0.01 (compared to M0) when performing Nemenyi post hoc test for MARGP (i.e., mean of ARG proportion) or MRCS (i.e., mean read count per sample). Refer to the abbreviations section for a detailed explanation of the acronyms.

We further assessed each ARG by estimating the Sample Proportion of ARG Detection (SPAD), finding that most MDR-ARGs increased by over 30% at M2 and M4 compared to M0. Similarly, the MARGP and Mean Read Count per Sample (MRCS) reinforced the conclusion that many MDR-ARGs significantly increased due to metformin intake (Fig. 1c). In line with Fig. 1b, some peptide-ARGs also showed an increase (Fig. 1c, Supplementary Table S5). In contrast, the P group exhibited only a minor increase in the average abundance of many genes, not as pronounced as in the M group. These findings suggest that metformin may lead to more significant changes in the ARG abundance pattern compared to placebo, particularly the increase in MDR-ARGs.

The effect of metformin on ARG abundance varies among individuals

The analysis of mean abundance in Fig. 1 does not fully represent the characteristics at each time point due to the heterogeneity of individual samples. In other words, the detailed pattern of ARG abundance change over time may differ between samples (Supplementary Fig. S4, Supplementary Table S1). To investigate this idea, we selected 13 ARGs with significant changes in the M group compared to the P group (Fig. 1c), and categorized individuals in the M group into sub-categories based on the changing pattern of these 13 ARG abundance in each sample. As a result, we identified a group of samples that increased at M2 and decreased at M4, a group that showed the most increase or remained stable at M4, and a group exhibiting no change, labeled as g1, g2, and g3, respectively (Fig. 2a, see “Methods”). Such categorization was not applicable to the P group, which served as a common control in this analysis. The characteristics of each subcategory, as determined by our criteria, became more prominent not only at the ARG level but also at the DC level. In g1, ARGs associated with MDR, peptide, and phosphonic-acid DCs significantly increased initially and then decreased, whereas in g2, after an initial significant increase, they further increased or were maintained. For g3, there were no significant changes (Fig. 2b). This result indicates that metformin can indeed increase the abundance of ARGs, especially MDR-ARGs, but not uniformly across all individuals treated with metformin. Notably, the outlier samples (i.e., two samples labeled P and three samples labeled M groups) with high baseline MDR abundance showed a decreased abundance of ARGs following metformin treatment (Supplementary Figs. S1, S5).Fig. 2 Analysis of ARG abundance patterns in different individuals. (a) Heatmap showing ARGP (i.e., ARG proportions) for those with significant changes in abundance in Fig. 1. The subcategories of samples, g1, g2, g3, were determined based on the individual differences in the abundance pattern of ARGs (see “Methods”). (b) Boxplots showing patterns of ARGP-DC (i.e., total ARG proportions in drug class) along with schematic drawings reflecting the ARG abundance at an individual level after the subcategorization of samples. Boxes represent the interquartile range (IQR) between the 25th and 75th percentiles, and the horizontal line inside the box denotes the median value. Whiskers represent the lowest and highest values within 1.5 times the IQR. The color inside the circle of the schematic drawings varies according to the magnitude of ARGP-DC value, labeled the ranges on the right side of the figure. For the P group, the average value by group for all time points is shown. "*" means P < 0.05, and "**" means P < 0.01 according to the Mann–Whitney U test. Refer to the abbreviations section for a detailed explanation of the acronyms.

E. coli correlation with ARG trends depending on metformin intake

We next explored which changes in bacterial composition are related to alterations in ARG abundance due to metformin intake. After assigning taxonomic information to all reads using Kraken212, leading to the identification of a total of 444 bacterial species (Supplementary Table S6), we estimated bacterial compositions across different time points of metformin intake. We then assessed whether the abundance of ARGs, belonging to different DCs was associated with the species proportions per sample (SPS) for each of the significantly altered six bacterial species. We found that E. coli and B. wexlerae had significant correlations with many MDR-ARGs (Pearson correlation coefficient > 0.5 and P < 0.05) (Fig. 3a, Supplementary Table S7), indicating that the compositional changes in E. coli and B. wexlerae may be related to the composition of MDR-ARGs. Specifically, the composition of E. coli showed the strongest Pearson correlation with the abundance of each ARG (average correlation coefficient of 0.74) (Fig. 3a, Supplementary Table S7), followed by B. wexlerae with the next strongest correlation (average correlation coefficient of 0.5).Fig. 3 Analysis of bacterial abundance patterns correlated with ARGs. (a) Heatmap showing Pearson correlation between SPS (i.e., species proportions per sample) and ARGP (i.e., ARG proportions) for each ARG. The colors in the cell of heatmap indicate the Pearson correlation value between each bacterial species composition and ARG abundance. The "*" in each cell indicates Pearson correlation coefficient > 0.5 and p-value < 0.05. (b) Heatmap showing SPS (%) for bacteria with significant changes by metformin. The Y-axis is marked with species changed by metformin. The colors in the heatmap vary according to the magnitude of the SPS values. Note that B. wexlerae is marked with a different color range due to visual issues, indicated by an asterisk. Refer to the abbreviations section for a detailed explanation of the acronyms.

We further analyzed each species abundance patterns in the subcategorized samples as the same way that we did for Fig. 2, to see whether the sample heterogeneity pattern was also consistent with what was observed for the ARGs in Fig. 2. Most species did not match the ARG pattern within each subcategory, although showed significant differences in abundance between the time points of M groups. Specifically, the pathogenic bacterium R. ilealis13, significantly decreased in abundance in many samples following metformin intake (Fig. 3b, Supplementary Table S8), and a beneficial bacterium B. wexlerae showed an increasing trend due to metformin (Fig. 3b, Supplementary Fig. S6)14. However, contrary to previous reports6,7, A. muciniphila showed no change in the M2 and M4 groups (Supplementary Fig. S6).

Except for R. ilealis and B. wexlerae, the patterns of change in bacterial compositions did not differ between the M and P groups, suggesting that bacterial composition changes may be independent of metformin intake. In contrast, the composition of E. coli changed in both M and P groups, with the trend being substantially more pronounced in the M groups than in the P groups (Fig. 3b). Moreover, the composition of E. coli was strongly consistent with the changes in the MDR-ARGs in the subcategorized samples; in the g1 subgroup, E. coli significantly increased in M2 and then decreased in M4, while in g2, it increased in M2 (P < 0.05 in all three cases described above, Mann–Whitney U test) and then either increased further or remained the same in M4 (Fig. 3b). B. wexlerae also showed a pattern similar to E. coli, but the pattern was weaker (Fig. 3b). These findings suggest that metformin may specifically be related to the alterations in E. coli composition associated with MDR-ARGs (Supplementary Fig. S7).

The contig-based analysis confirms the effect of metformin on the increase of MDR-ARGs

We next examined if the findings from the read-based analyses were replicated in the contig-based analysis. Contig assembly and genome binning were conducted on all reads following conventional methods described by others15–17 (see “Methods”). As a result, a total of 134 ARGs belonging to 15 DCs were identified. Subsequently, after calculating Normalized Read Depth coverage (NRD_cov), we computed MARGP-DC as was done in Fig. 1 (see “Methods”, Supplementary Fig. S8). In general, the results were consistent with those derived from the read-based analysis presented in Fig. 1. Specifically, the MARGP-DC of the MDR-ARGs increased in M2 and M4, with more pronounced than P2 and P4 (Fig. 4a, Supplementary Table S4). Although the detection of ARGs in fewer samples compared to the read-based analysis made it challenging to ascertain statistical significance (Fig. 4b, Supplementary Table S10), the subcategorization analysis also confirmed the sample heterogeneity in the presence of ARGs as observed in the read-based analysis (Fig. 4b).Fig. 4 Contig-based Analysis of ARGs altered by metformin. (a) Stacked bar plot showing the MARGP-DC (i.e., mean of total ARG proportions in drug class) per sample for each DC. In this analysis, the ARGP (i.e., ARG proportions) is calculated based on the NRD_cov (i.e., normalized read depth coverage) of each ARG containing contig. (b) Heatmap of NRD_cov for each ARG with significant changes in abundance as shown in Fig. 1. The colors in the heatmap varies depending on the magnitude of NRD_cov value. (c) Heatmap of NRD_cov for each bacterial species which abundance was changed significantly in the read-based analysis. Only samples in which MDR-ARGs were detected are shown in this heatmap, and E. coli was detected only in these samples. (d) Bar plot showing Mean NRD_cov per sample for species-specific bins harboring MDR-ARGs. The pie chart on the right side showing the ratio of DC (based on Mean NRD_cov) observed in E. coli. The color of each DC is labeled in the middle panel of this figure. Refer to the abbreviations section for a detailed explanation of the acronyms.

Next, we attempted to estimate bacterial composition by mapping data from single copy marker genes (SCGs), to investigate which bacterial species carry those ARGs. Bacterial species were identified by aligning contigs to a set of SCG sequences, which led to the creation of species-specific bins mapped to the same marker genes. A total of 7,201 bins were created, with 3795 bins assigned to 236 species (see “Methods”). The NRD_cov value was then calculated for each bin associated with bacteria in individual samples. E. coli was found in the same samples where MDR-ARGs had increased, and its pattern closely matched the pattern of MDR-ARGs observed in the subcategorization analysis (Fig. 4b,c). B. wexlerae also showed patterns similar to those MDR-ARGs demonstrated in the subcategorization analysis, although the trend was much weaker than that of E. coli (Supplementary Fig. S9). To identify which of the two species was more highly correlated with the trends in these MDR-ARGs, all contigs carrying ARGs were mapped to the species-specific bins. As a result, it was determined that E. coli had the highest correlation with the MDR-ARGs (Fig. 4d). Notably, a few MDR-ARGs detected in the P group were mapped to species other than E. coli (Supplementary Fig. S10).

Discussion

In the present work, we investigated the relationship between bacterial compositions and ARGs altered by metformin, using the metagenome data obtained from the study by Wu et al.6, based on two methods: read-based and contig-based analyses. We found that both analyses yielded consistent results, despite the limitation of identifying fewer ARGs mapped to E. coli bins than expected in the contig-based analysis. We concluded that metformin can indeed increase ARG abundance, but not uniformly across all individuals treated with metformin. Our results suggest that the ARG trends vary among individuals and that E. coli is likely the species most closely corresponding to the trends of ARGs, particularly those belonging to MDR genes, in patients with T2D.

Wu et al.6 demonstrated that metformin alters the microbiome composition in patients with T2D, suggesting that these changes could have beneficial effects on improving T2D. Consistent with this, our current study revealed that the beneficial bacterium B. wexlerae, known to transform the gut environment including SCFA composition and reduce the incidence of obesity and diabetes in Japanese adults14, showed increasing trends with metformin intake, while the pathogenic bacterium R. ilealis, prevalent in over 80% of obese patients and reported to worsen glucose metabolism13, decreased following metformin intake (Fig. 3b). However, A. muciniphila, a bacterium well-known for its beneficial effects, did not show significant changes with metformin intake, contrary to previous reports6,7. This inconsistency may be due to the differences in the number of identified bacterial species in each sample, as the number of bacterial species identified in the present study was much smaller than that reported by Wu et al.6.

Our study suggests that metformin can alter bacterial compositions in ways that may impact human health both positively and negatively. Specifically, an increase in B. wexlerae could be associated with beneficial effects, potentially reducing the risk of obesity and T2D. However, an increase in ARGs associated with E. coli could be harmful. Increases in E. coli following metformin intake have been observed in T2D patients18,19, yet this finding is not without contradiction: an in vitro experiment indicated that metformin did not directly affect the growth of E. coli6. Most of the MDR genes identified in this study (Fig. 2) as increased by metformin treatment are chromosomally encoded, mainly encoding efflux pump components and their transcriptional regulators. Therefore, the rise in MDR gene abundance merely due to the increased prevalence of E. coli carrying these genes may not directly cause an enhancement in antibiotic resistance. Instead, it could involve the spontaneous or mutational upregulation of these MDR genes in E. coli. Another potential cause of resistance could be the increased prevalence of E. coli strains carrying MDR genes in their genomes. According to the CARD database, approximately 61.6% of E. coli NCBI WGS samples carry TolC (https://card.mcmaster.ca/ontology/36376), indicating that about 40% of E. coli may not carry MDR genes in their genomes. Further experimental studies are needed to clarify this issue.

According to Magruder et al.’s study based on meta-analysis, patients with T2D are twice as likely to experience urinary tract infections (UTIs), and most bacterial infections originating in the urinary tract are often presumed to originate from the gut. Indeed, one study demonstrated that E. coli strains from urine samples were most similar to E. coli in fecal samples from the same subjects20. This leads us to speculate that the increase in E. coli during metformin intake could elevate the risk of UTIs or potentially lead to the risk associated with the spread of MDR in the gut.

Finally, the observation that the effects of metformin vary individually as shown in Fig. 2 may provide novel insight into how patients with T2D can be distinguished and treated more effectively with future metformin-based therapies. Understanding the relationship between the abundance of B. wexlerae or E. coli and clinical blood test data will be essential. Nonetheless, caution is warranted in interpreting the clinical implications of these findings, given the limited sample sizes in this study. Extensive, large-scale research is necessary to confirm how metformin influences individuals differently and to elucidate the mechanisms behind metformin-induced increases in E. coli, which likely contribute to the rise in multi-drug resistant antibiotic resistance genes.

Conclusions

Our study suggests that metformin treatment in patients with T2D can affect the gut antibiotic resistome, particularly by increasing the abundance of MDR-ARGs, with E. coli playing a central role in this process. The findings underscore the complexity of metformin's effects on the gut microbiome and its potential implications for antibiotic resistance.

Methods

Data preparation

In this study, we utilized publicly available data primarily produced by the research conducted by Wu et al.6. Briefly, WGS data from fecal samples of T2D patients (n = 40), divided into a group that took metformin (M group, n = 22) and a group that took a placebo drug (P group, n = 18). Fecal samples from the patients were collected before medication (M0, P0), two months after taking the medication (M2, P2), and four months after taking the medication (M4, P4). According to Wu et al.’s6, the cohort consists of patients who did not have any systemic or metabolic diseases other than T2D, had not an infection within a month and had not taken medications such as antibiotics for three months. The metagenome data was downloaded from NCBI SRA database (ID: PRJNA361402). Data analysis was performed with FASTQ format by converting the files in SRA format using the fastq-dump-orig tool of the sratoolkit (v3.0.0) (https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software).

Data preprocessing

FastQC (v0.11.9) (https://www.bioinformatics.babraham.ac.uk/) and MultiQC (v1.11)21 were used for quality control. Next, sequencing reads with a quality score below 20 and a minimum read length of less than 35 bp were removed using fastp (v0.23.0)22, with all other thresholds set to default values. After using fastp, it was verified using MultiQC that the adapter content of all samples was less than 0.1%. All remaining reads were then aligned to the human reference genome FASTA file (GRCh38/hg38) using bowtie2 (v2.4.5)23. Finally, only the unaligned reads from the bowtie2 output were used for read-based analysis, also used when performing contig-assembly.

Contig-assembly and binning of contigs

ATLAS (v2.18.1)24 was used for contig-based analyses. Within ATLAS, the built-in tools called metaSPAdes25 and MetaBAT226 were used for contig-assembly and binning method, respectively. Only contigs with a minimum length of 1500 bp were only selected for further analysis, with default settings for cutoffs and commands implemented in ATLAS except –skip-qc option. About 2 million contigs were produced from the assembly process, which were then used as input for ARG detection and taxon identification. Approximately 1.5 million contigs were used for binning to identify taxon identification.

ARG identification and filtering

In read-based analysis, all reads processed through the steps mentioned above were aligned to the reference database, CARD (v3.2.5)27, using RGI bwt (v6.0.0)27 integrated with KMA (v1.4.9)28. Only reads with an identity of 90% or higher and mapping quality (MAPQ) of 30 or above were used for analysis. The breadth coverage value, indicating the extent to which the entire sequence of a specific ARG is covered by reads mapped to it within a sample, was calculated using mpileup from bcftools (v1.11)29. If the breadth coverage of a specific ARG in a sample was less than 60%, that ARG was considered not identified in that sample. To minimize noise, any ARG with SPAD below 20% across all samples was removed from the analysis. We also excluded ARGs that confer antibiotic resistance due to SNP occurrence (for example, E. coli gyrA conferring resistance to fluoroquinolones), because the RGI tool has yet to implement the function to detect such genes (https://github.com/arpcard/rgi/blob/dev/README.rst). In the end, 102 ARGs were identified. In contig-based analysis, all contigs processed through the steps mentioned above were aligned to the CARD using RGI main (v6.0.0) integrated with DIAMOND (v0.8.36)30. Only contigs with an identity of 90% or higher and breadth coverage of 60% or above were used for analysis. Finally, 134 ARGs were identified.

Abundance calculation

In read-based analysis, the raw read count of ARGs was normalized using the FPKM calculation formula, and ARGP was calculated based on the relative abundance obtained from normalized read counts.ARG read count=ARG raw read count per sampleARG length∗sample total read count∗109

In contig-based analysis, we calculated NRD_cov which explain read mapping level for taxon-contig in each bins or ARG-contig. To calculate average depth across all sequences in contig, we used ‘jgi summarize bam contig depths’ from MetaBAT226 built into ATLAS (v2.18.1)24, and used –percentIdentity 95 option for consider only reads mapped to target contig by identity 95%. Additionally, when calculate depth about ARG-contig, we considered only read aligned to the ‘specific-one ARG location’ estimated by rgi, not the entire location of contig to which can multiple ARG is assigned. Reads that align to the edge parts of a ‘specific-one ARG’ location are also be calculated for depth (More than 60% of each read's length must be aligned). Finally, we calculated NRD_cov following formula described by other study31.NRD_cov=average depth in sequences∗150length(ARG or taxon contig)∗sample total read count∗106

Assigning ARGs into DC

We manually curated which DC each ARG confers resistance to. The process began with searching for all ARGs in the CARD database27. After verifying the "Drug Class," "Definition," and "Publications" annotations for each ARG in the database, we referred to the curation methods of other studies32 and reassigned the DC related to each ARG using the following rules: MLS (if at least one of streptogramin, streptogramin A, streptogramin B, lincosamide, or macrolide is listed in "Drug Class"), beta-lactam (if at least one of cephalosporin, cephamycin, monobactam, penam or carbapenem is listed in "Drug Class"), other classes (following CARD's curation), MDR (if it belongs to at least three different curated classes or if the "Definition" or "Publications" list mentions that the gene is associated with multidrug resistance). As a result of applying these curation rules, all ARGs used in the analysis were assigned to a total of 11 DCs in read-based analysis, and 15 DCs in contig-based analysis.

Sample subcategorization criteria

Sample subcategorization was conducted according to the following criteria; g1: Samples carrying ARGs which abundance increased by ≥ 3.5-fold in M2 and decreased by ≥ 3.5-fold in M4 compared to M2; g2: Samples carrying ARGs which abundance increased by ≥ 3.5-fold in M2 and either increased by ≥ 3.5-fold or decreased to < 1-fold in M4 compared to M2. An additional criterion was whether the number of ARGs in g1 or g2 subcategories in a specific sample exceeded 60% of target ARGs number. The target ARGs refer to ARGs that show significant changes in M group compared to P group (Fig. 1c). If a sample did not meet any of these conditions, it was assigned to the g3 category.

Taxon identification

For the read-based analysis, each read from all the samples that had completed sample cleaning process was assigned to one of 5213 species-level taxa using Kraken2 (v2.1.2)12. After filtering out the species with the relative abundance of < 0.01% as described in Refs.12,33, the cleaned reads were ultimately assigned to a total of 444 species (Supplementary Table S9). For the contig-based analysis, after gene prediction within all assembled contigs of each sample using Prodigal (v2.6.3)34, using the gene-predicted contig sequences as queries and the bac120 SCG set35 from the gtdb (release 214.1)36 as the subject, we conducted BLASTN (v2.9.0)37. Subsequently, we assigned the best hit SCG to each gene-predicted contig sequence (alignment identity = 100%, breadth coverage ≥ 60%). Then, we counted the number of each SCG present in all contigs within each bin. We designated species name by aligning the bin sequence with the SCG in the gtdb database, when the proportion of a specific SCG within a specific bin was over 80%. Consequently, a total of 3795 bins across all samples were assigned to 236 species.

Data analysis and visualization

PCoA was performed using the vegan (v2.5.7), ape (v5.8) package in R. All datasets targeted for analysis were first subjected to the Shapiro–Wilk test to assess normality. Datasets not demonstrating normality were analyzed with the Friedman test and the Nemenyi post hoc test or Mann–Whitney U test to identify statistically significant differences between time points. Statistical analyses were conducted using scipy (v1.10.1), scikit-posthocs (v0.7.0), numpy (v1.24.2) and pandas (v2.0.0) library in Python (v3.11.2). Visualization was carried out using the ggplot2 (v3.3.2) and complexheatmap (v2.2.0) packages in R.

Supplementary Information

Supplementary Figures.

Supplementary Tables.

Abbreviations

ARG Antibiotic resistance gene

ARGP ARG proportions; this term represents the concept of ARG relative abundance.

ARGP-DC Total ARG proportions in drug class; e.g., “ARGP-DC of peptide class = 3%” means “the abundance of ARGs conferring resistance to a peptide class total 3% of all ARG abundance in that sample”.

DC Drug class; a term used to indicate the drug class to which a particular ARG confers resistance.

MARGP Mean of ARG proportions in a specific time point (e.g., M2)

MARGP-DC Mean ARGP-DC; this term refers to the mean of ARGP-DC values for samples at a specific time point (e.g., M2).

MRCS Mean read count per sample; this term refers to the average count of ARG reads at a specific time point (e.g., M2). The count of ARG reads was normalized using an FPKM-based calculation.

NRD_cov NoRmalized Read Depth coverage; this measure represents abundance in assembly-based analysis. For detailed calculation formulas, refer to the “Methods” section

SPAD Sample proportions of ARG detection; this term refers to the detection frequency of a specific ARG (e.g., “SPAD of TolC gene = 20% at time point M2” means “the TolC gene was detected in 20% of all samples from M2”)

SPS Species proportions per sample; this term represents the concept of species relative abundance

Supplementary Information

The online version contains supplementary material available at 10.1038/s41598-024-72467-z.

Author contributions

SSC supervised the research and wrote the manuscript. HBK conducted the analysis and also contributed to writing the manuscript. YJC participated in discussion related to study design and assisted in the writing process.

Funding

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (RS-2024-00341909), Ministry of Science and ICT, Korea (RS-2023-00260267), and the Korea National Institute of Health research project (#2024ER210600).

Data availability

We downloaded and used data from the SRA with accession number PRJNA361402.

Competing interests

The authors declare no competing interests.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Ahmad E Where does metformin stand in modern day management of type 2 diabetes? Pharmaceuticals 2020 13 12 427 10.3390/ph13120427 33261058
Ahmad, E. et al. Where does metformin stand in modern day management of type 2 diabetes?. Pharmaceuticals 13(12), 427 (2020).33261058 10.3390/ph13120427
2. Zhang Q Hu N Effects of metformin on the gut microbiota in obesity and type 2 diabetes mellitus Diabetes Metab. Syndr. Obes. 2020 13 5003 5014 10.2147/DMSO.S286430 33364804
Zhang, Q. & Hu, N. Effects of metformin on the gut microbiota in obesity and type 2 diabetes mellitus. Diabetes Metab. Syndr. Obes. 13, 5003–5014 (2020).33364804 10.2147/DMSO.S286430
3. Silamiķele L Metformin strongly affects gut microbiome composition in high-fat diet-induced type 2 diabetes mouse model of both sexes Front. Endocrinol. 2021 12 626359 10.3389/fendo.2021.626359
Silamiķele, L. et al. Metformin strongly affects gut microbiome composition in high-fat diet-induced type 2 diabetes mouse model of both sexes. Front. Endocrinol. 12, 626359 (2021).10.3389/fendo.2021.626359
4. He D Metformin reduces blood glucose in treatment-naive type 2 diabetes by altering the gut microbiome Can. J. Diabetes 2022 46 2 150 156 10.1016/j.jcjd.2021.08.001 35148952
He, D. et al. Metformin reduces blood glucose in treatment-naive type 2 diabetes by altering the gut microbiome. Can. J. Diabetes 46(2), 150–156 (2022).35148952 10.1016/j.jcjd.2021.08.001
5. Ezzamouri B Metabolic modelling of the human gut microbiome in type 2 diabetes patients in response to metformin treatment NPJ Syst. Biol. Appl. 2023 9 1 2 10.1038/s41540-022-00261-6 36681701
Ezzamouri, B. et al. Metabolic modelling of the human gut microbiome in type 2 diabetes patients in response to metformin treatment. NPJ Syst. Biol. Appl. 9(1), 2 (2023).36681701 10.1038/s41540-022-00261-6
6. Wu H Metformin alters the gut microbiome of individuals with treatment-naive type 2 diabetes, contributing to the therapeutic effects of the drug Nat. Med. 2017 23 7 850 858 10.1038/nm.4345 28530702
Wu, H. et al. Metformin alters the gut microbiome of individuals with treatment-naive type 2 diabetes, contributing to the therapeutic effects of the drug. Nat. Med. 23(7), 850–858 (2017).28530702 10.1038/nm.4345
7. De La Cuesta-Zuluaga J Metformin is associated with higher relative abundance of mucin-degrading Akkermansia muciniphila and several short-chain fatty acid–producing microbiota in the gut Diabetes Care 2017 40 1 54 62 10.2337/dc16-1324 27999002
De La Cuesta-Zuluaga, J. et al. Metformin is associated with higher relative abundance of mucin-degrading Akkermansia muciniphila and several short-chain fatty acid–producing microbiota in the gut. Diabetes Care 40(1), 54–62 (2017).27999002 10.2337/dc16-1324
8. Maier L Extensive impact of non-antibiotic drugs on human gut bacteria Nature 2018 555 7698 623 628 10.1038/nature25979 29555994
Maier, L. et al. Extensive impact of non-antibiotic drugs on human gut bacteria. Nature 555(7698), 623–628 (2018).29555994 10.1038/nature25979
9. Nagata N Population-level metagenomics uncovers distinct effects of multiple medications on the human gut microbiome Gastroenterology 2022 163 4 1038 1052 10.1053/j.gastro.2022.06.070 35788347
Nagata, N. et al. Population-level metagenomics uncovers distinct effects of multiple medications on the human gut microbiome. Gastroenterology 163(4), 1038–1052 (2022).35788347 10.1053/j.gastro.2022.06.070
10. Vich Vila A Impact of commonly used drugs on the composition and metabolic function of the gut microbiota Nat. Commun. 2020 11 1 362 10.1038/s41467-019-14177-z 31953381
Vich Vila, A. et al. Impact of commonly used drugs on the composition and metabolic function of the gut microbiota. Nat. Commun. 11(1), 362 (2020).31953381 10.1038/s41467-019-14177-z
11. Lance GN Williams WT Computer programs for hierarchical polythetic classification (“similarity analysis”) Computer J. 1966 9 1 60 64 10.1093/comjnl/9.1.60
Lance, G. N. & Williams, W. T. Computer programs for hierarchical polythetic classification (“similarity analysis”). Computer J. 9(1), 60–64 (1966).10.1093/comjnl/9.1.60
12. Wood DE Lu J Langmead B Improved metagenomic analysis with Kraken2 Genome Biol. 2019 20 1 13 10.1186/s13059-019-1891-0 30606230
Wood, D. E., Lu, J. & Langmead, B. Improved metagenomic analysis with Kraken2. Genome Biol. 20, 1–13 (2019).30606230 10.1186/s13059-019-1891-0
13. Rodrigues RR Transkingdom interactions between Lactobacilli and hepatic mitochondria attenuate western diet-induced diabetes Nat. Commun. 2021 12 1 101 10.1038/s41467-020-20313-x 33397942
Rodrigues, R. R. et al. Transkingdom interactions between Lactobacilli and hepatic mitochondria attenuate western diet-induced diabetes. Nat. Commun. 12(1), 101 (2021).33397942 10.1038/s41467-020-20313-x
14. Hosomi K Oral administration of Blautia wexlerae ameliorates obesity and type 2 diabetes via metabolic remodeling of the gut microbiota Nat. Commun. 2022 13 1 4477 10.1038/s41467-022-32015-7 35982037
Hosomi, K. et al. Oral administration of Blautia wexlerae ameliorates obesity and type 2 diabetes via metabolic remodeling of the gut microbiota. Nat. Commun. 13(1), 4477 (2022).35982037 10.1038/s41467-022-32015-7
15. Liao J Microdiversity of the vaginal microbiome is associated with preterm birth Nat. Commun. 2023 14 4997 10.1038/s41467-023-40719-7 37591872
Liao, J. et al. Microdiversity of the vaginal microbiome is associated with preterm birth. Nat. Commun. 14, 4997 (2023).37591872 10.1038/s41467-023-40719-7
16. Gálvez EJ Distinct polysaccharide utilization determines interspecies competition between intestinal Prevotella spp Cell Host Microbe 2020 28 6 838 852 10.1016/j.chom.2020.09.012 33113351
Gálvez, E. J. et al. Distinct polysaccharide utilization determines interspecies competition between intestinal Prevotella spp. Cell Host Microbe 28(6), 838–852 (2020).33113351 10.1016/j.chom.2020.09.012
17. Chevalier C Warmth prevents bone loss through the gut microbiota Cell Metab. 2020 32 4 575 590 10.1016/j.cmet.2020.08.012 32916104
Chevalier, C. et al. Warmth prevents bone loss through the gut microbiota. Cell Metab. 32(4), 575–590 (2020).32916104 10.1016/j.cmet.2020.08.012
18. Mueller NT Metformin affects gut microbiome composition and function and circulating short-chain fatty acids: A randomized trial Diabetes Care 2021 44 7 1462 1471 10.2337/dc20-2257 34006565
Mueller, N. T. et al. Metformin affects gut microbiome composition and function and circulating short-chain fatty acids: A randomized trial. Diabetes Care 44(7), 1462–1471 (2021).34006565 10.2337/dc20-2257
19. Zhernakova A Population-based metagenomics analysis reveals markers for gut microbiome composition and diversity Science 2016 352 6285 565 569 10.1126/science.aad3369 27126040
Zhernakova, A. et al. Population-based metagenomics analysis reveals markers for gut microbiome composition and diversity. Science 352(6285), 565–569 (2016).27126040 10.1126/science.aad3369
20. Magruder M Gut uropathogen abundance is a risk factor for development of bacteriuria and urinary tract infection Nat. Commun. 2019 10 1 5521 10.1038/s41467-019-13467-w 31797927
Magruder, M. et al. Gut uropathogen abundance is a risk factor for development of bacteriuria and urinary tract infection. Nat. Commun. 10(1), 5521 (2019).31797927 10.1038/s41467-019-13467-w
21. Ewels P Magnusson M Lundin S Käller M MultiQC: Summarize analysis results for multiple tools and samples in a single report Bioinformatics 2016 32 19 3047 3048 10.1093/bioinformatics/btw354 27312411
Ewels, P., Magnusson, M., Lundin, S. & Käller, M. MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32(19), 3047–3048 (2016).27312411 10.1093/bioinformatics/btw354
22. Chen S Zhou Y Chen Y Gu J fastp: An ultra-fast all-in-one FASTQ preprocessor Bioinformatics 2018 34 17 i884 i890 10.1093/bioinformatics/bty560 30423086
Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34(17), i884–i890 (2018).30423086 10.1093/bioinformatics/bty560
23. Langmead B Salzberg S Fast gapped-read alignment with Bowtie 2 Nat. Methods 2012 9 4 357 359 10.1038/nmeth.1923 22388286
Langmead, B. & Salzberg, S. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9(4), 357–359 (2012).22388286 10.1038/nmeth.1923
24. Kieser S Brown J Zdobnov EM Trajkovski M McCue LA ATLAS: A Snakemake workflow for assembly, annotation, and genomic binning of metagenome sequence data BMC Bioinform. 2020 21 1 8 10.1186/s12859-020-03585-4
Kieser, S., Brown, J., Zdobnov, E. M., Trajkovski, M. & McCue, L. A. ATLAS: A Snakemake workflow for assembly, annotation, and genomic binning of metagenome sequence data. BMC Bioinform. 21, 1–8 (2020).10.1186/s12859-020-03585-4
25. Nurk S Meleshko D Korobeynikov A Pevzner PA metaSPAdes: A new versatile metagenomic assembler Genome Res. 2017 27 5 824 834 10.1101/gr.213959.116 28298430
Nurk, S., Meleshko, D., Korobeynikov, A. & Pevzner, P. A. metaSPAdes: A new versatile metagenomic assembler. Genome Res. 27(5), 824–834 (2017).28298430 10.1101/gr.213959.116
26. Kang DD MetaBAT 2: An adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies PeerJ 2019 7 e7359 10.7717/peerj.7359 31388474
Kang, D. D. et al. MetaBAT 2: An adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ 7, e7359 (2019).31388474 10.7717/peerj.7359
27. Alcock BP CARD 2020: Antibiotic resistome surveillance with the comprehensive antibiotic resistance database Nucleic Acids Res. 2020 48 D1 D517 D525 31665441
Alcock, B. P. et al. CARD 2020: Antibiotic resistome surveillance with the comprehensive antibiotic resistance database. Nucleic Acids Res. 48(D1), D517–D525 (2020).31665441
28. Clausen PTLC Aarestrup FM Lund O Rapid and precise alignment of raw reads against redundant databases with KMA BMC Bioinform. 2018 19 307 307 10.1186/s12859-018-2336-6
Clausen, P. T. L. C., Aarestrup, F. M. & Lund, O. Rapid and precise alignment of raw reads against redundant databases with KMA. BMC Bioinform. 19(307), 307 (2018).10.1186/s12859-018-2336-6
29. Li H The sequence alignment/map format and SAMtools Bioinformatics 2009 25 16 2078 2079 10.1093/bioinformatics/btp352 19505943
Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25(16), 2078–2079 (2009).19505943 10.1093/bioinformatics/btp352
30. Buchfink B Reuter K Drost HG Sensitive protein alignments at tree-of-life scale using DIAMOND Nat. Methods 2021 18 4 366 368 10.1038/s41592-021-01101-x 33828273
Buchfink, B., Reuter, K. & Drost, H. G. Sensitive protein alignments at tree-of-life scale using DIAMOND. Nat. Methods 18(4), 366–368 (2021).33828273 10.1038/s41592-021-01101-x
31. Xiong W Antibiotic-mediated changes in the fecal microbiome of broiler chickens define the incidence of antibiotic resistance genes Microbiome 2018 6 1 1 11 10.1186/s40168-018-0419-2 29291746
Xiong, W. et al. Antibiotic-mediated changes in the fecal microbiome of broiler chickens define the incidence of antibiotic resistance genes. Microbiome 6(1), 1–11 (2018).29291746 10.1186/s40168-018-0419-2
32. Zhang Z Assessment of global health risk of antibiotic resistance genes Nat. Commun. 2022 13 1 1553 10.1038/s41467-022-29283-8 35322038
Zhang, Z. et al. Assessment of global health risk of antibiotic resistance genes. Nat. Commun. 13(1), 1553 (2022).35322038 10.1038/s41467-022-29283-8
33. Tan CC No evidence for a common blood microbiome based on a population study of 9,770 healthy humans Nat. Microbiol. 2023 8 5 973 985 10.1038/s41564-023-01350-w 36997797
Tan, C. C. et al. No evidence for a common blood microbiome based on a population study of 9,770 healthy humans. Nat. Microbiol. 8(5), 973–985 (2023).36997797 10.1038/s41564-023-01350-w
34. Hyatt D Prodigal: Prokaryotic gene recognition and translation initiation site identification BMC Bioinform. 2010 11 1 11 10.1186/1471-2105-11-119
Hyatt, D. et al. Prodigal: Prokaryotic gene recognition and translation initiation site identification. BMC Bioinform. 11, 1–11 (2010).10.1186/1471-2105-11-119
35. Parks DH Recovery of nearly 8000 metagenome-assembled genomes substantially expands the tree of life Nat. Microbiol. 2017 2 11 1533 1542 10.1038/s41564-017-0012-7 28894102
Parks, D. H. et al. Recovery of nearly 8000 metagenome-assembled genomes substantially expands the tree of life. Nat. Microbiol. 2(11), 1533–1542 (2017).28894102 10.1038/s41564-017-0012-7
36. Parks DH A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life Nat. Biotechnol. 2018 36 10 996 1004 10.1038/nbt.4229 30148503
Parks, D. H. et al. A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat. Biotechnol. 36(10), 996–1004 (2018).30148503 10.1038/nbt.4229
37. Camacho C BLAST+: Architecture and applications BMC Bioinform. 2009 10 1 9 10.1186/1471-2105-10-421
Camacho, C. et al. BLAST+: Architecture and applications. BMC Bioinform. 10, 1–9 (2009).10.1186/1471-2105-10-421
