Since the animal experiment was conducted for simultaneous gut microbiome data acquisition, the study was designed as a randomized block experiment, where each study arm was represented by three experimental units in specific pathogen-free conditions. An animal cage with three animals in it was considered as one experimental unit of the study due to an inevitable exchange of the microbiota between different animals within one cage. Experimental units/cages were organized in three blocks (24 animals per block) with time as the only blocking factor and a 2-week-long time shift between each block resulting in a 4-week-long extension of the experiment in total. In total, 72 animals (36 males and 36 females) were involved in the study compiling 24 experimental units. Half of the animals (n = 36) were allocated to the high-fat diet intervention, while the rest of the animals (n = 36) received the control diet after the first randomization procedure. The 20-week-long diet intervention was followed by a 10-week-long metformin treatment to which half of the animals (n = 16) of each diet group were randomly subjected. An equal animal sex distribution was ensured by performing the randomization procedures in males and females separately. Finally, there were eight study groups developed with three experimental units (nine animals) in each group based on three different factors (sex, diet, and therapeutic intervention): 1) high-fat diet-fed (HFD-fed) male mice receiving metformin, 2) HFD-fed male mice with no metformin intervention, 3) HFD-fed female mice receiving metformin, 4) HFD-fed female mice with no metformin intervention, 5) control diet-fed (CD-fed) male mice receiving metformin, 6) CD-fed male mice with no metformin intervention, 7) CD-fed female mice receiving metformin, and 8) CD-fed female mice with no metformin intervention. Nevertheless, for the RNA-Seq analysis, only one animal from each cage was randomly selected, resulting in the intestinal transcriptome data set derived from six animals (three males and three females) of each of the treatment groups: high-fat diet-fed animals (HFD), high-fat diet-fed animals receiving metformin (HFDmetf), control diet-fed animals (CD), and control diet-fed animals receiving metformin (CDmetf), representing all three experimental blocks. A schematic representation of the experimental design is provided in Figure 1 .

Four- to 5-week-old male and female C57BL/6N mice were obtained from the University of Tartu Laboratory Animal Centre. Animals of the same sex were housed in individually ventilated cages (Tecniplast, Milan, Italy) in groups of three. Mice were fed with either a high-fat diet (HFD: 60 kcal% fat, D12492, Research Diets, New Jersey, USA) or a control diet (CD: 10 kcal% fat, D12450J, Research Diets, New Jersey, USA) ad libitum depending on the experimental group assigned. The room temperature was controlled at 23°C ± 2°C with 55% humidity and a 12:12-h light–dark cycle.

For half of the animals, an obesity and insulin resistance phenotype was induced by HFD for 20 weeks. After week 20, metformin therapy was provided for randomly selected animals (half of both HFD-fed and CD-fed mice) with drinking water for 10 weeks in a concentration of 50 mg/kg body mass/day. The dose of metformin was chosen according to previous findings proving metformin dosage of 50 mg/kg/day as a clinically relevant dose retaining the property to improve hyperglycemia in high-fat diet-fed animals (17). The phenotype of experimental animals was monitored: body weight was measured once a week, and 6-h fasting plasma glucose was determined using an Accu-Chek Performa glucometer (Roche, Vienna, Austria) once in 2 weeks. In addition, at weeks 20 (before administration of metformin) and 30 (after 10 weeks of metformin therapy), fasting glucose and insulin levels were determined using the mouse glucose assay and mouse insulin ELISA kit (Crystal Chem, Zaandam, Netherlands) to calculate the homeostatic model assessment for insulin resistance (HOMA-IR) (18). Stool samples were collected from animals before cervical dislocation and analyzed using the Mouse Secretory Immunoglobulin A (sIgA) ELISA Kit (MyBioSource, USA) following the manufacturer’s instructions, and the absorbance was read at 450 nm. All the animals were euthanized by cervical dislocation without any other anesthesia to avoid potential bias of medications other than metformin. For tissue sampling from the distal part of the small intestine (ileum), one animal representing each experimental unit/cage was randomly selected (24 animals in total). Each tissue sample was stored at −80°C in 1.5 mL of RNAlater solution (Thermo Fisher Scientific, USA) until further processing.

RNA extraction from intestinal tissue was performed using AllPrep DNA/RNA Mini Kit (QIAGEN, Germany) according to the manufacturer’s instructions. For rRNA depletion, 200 ng of total RNA was processed using MGIEasy rRNA depletion kit (MGI Tech Co. Ltd., China). Complementary DNA library preparation was performed with MGIEsy RNA Directional Library Prep Set (MGI Tech Co. Ltd., China). The quantity and quality of extracted RNA and prepared libraries were determined by Qubit Fluorometer (Thermo Fisher Scientific, USA) and Agilent 2100 Bioanalyzer system (Agilent, USA), respectively. The integrity of RNA was evaluated by RNA integrity number (RIN) within the Agilent 2100 Bioanalyzer system (Agilent, USA). Next-generation sequencing was conducted on the DNBSEQ-G400RS sequencing platform (MGI Tech Co. Ltd., China), following the manufacturer’s instructions. Since shotgun RNA sequencing is considered to be the most accurate method for quantification of the expression of individual transcripts and genes, additional methods for technical validation were not applied in this study (19).

Isolation of microbial DNA from stool samples collected after metformin intervention from the same animals and the following shotgun metagenomic sequencing was done as described before (9). Briefly, for the DNA extraction, FastDNA® Spin Kit for Soil (MP Biomedicals, USA) was used according to the manufacturer’s instructions. The DNA libraries were prepared with MGIEasy Universal DNA Library Prep Kit (MGI Tech Co. Ltd., China) adding 1% PhiX control and sequenced on the DNBSEQ-G400RS next-generation sequencing platform (MGI Tech Co. Ltd., China) using the DNBSEQ-G400RS High-throughput Sequencing Set (FCL PE100) (MGI Tech Co. Ltd., China). The desired sequencing depth was at least 20 M 100 bp paired-end sequencing reads per sample.

For the trimming of sequencing reads, a sliding window size of 5 and a quality threshold of 25 were applied in Trimmomatic 0.39. Sequencing reads had to have a minimum length of 30 bp and an average quality of 25 to be included in the subsequent analyses. The reads were mapped to the mouse reference genome GRCm38, and per-gene read counts were calculated with STAR 2.5.3a. The obtained read counts were then normalized using the trimmed mean of M-values implemented in Bioconductor package edgeR (v.3.24.3) in R (v.3.5.3). For the gene filtering, filterByExpr() function was applied, taking into account the library sizes and the experimental design (20). Differentially expressed genes (DEGs) were identified with limma package (v.3.38.3) and voom method using a linear model with weighted least squares for each gene (21). Both diet and therapy were defined as factors to evaluate their transcriptional effects, while the animal sex and sequencing block were considered covariates. DEGs were identified by comparing animals representing different experimental groups in four contrasts: HFD-fed animals against CD-fed animals with no metformin therapy applied in both groups (HFD vs. CD), HFD-fed animals against CD-fed animals with metformin therapy applied for both groups (HFDmetf vs. CDmetf), HFD-fed animals with metformin therapy against HFD-fed animals with no metformin therapy (HFDmetf vs. HFD), and finally, CD-fed animals with metformin therapy compared with CD-fed animals with no metformin therapy applied (CDmetf vs. CD). Multiple testing correction was implemented using the Benjamini–Hochberg procedure, and significant DEGs were determined using a false discovery rate (FDR) <0.05 cutoff (22).

Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using the online software Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 applying a P-value threshold <0.05 (23, 24). Heatmaps were constructed with Matplotlib⁠ and SciPy⁠. Hierarchical clustering was performed with the average linkage method implemented in SciPy for the clustering of genes according to their differences in CPM values (25, 26). Cnetplot for GO terms was constructed using the clusterProfiler package in R (27).

For the correlation analysis, features (genes and taxonomies) with a median value of 100 reads were retained for further analyses. Read counts were normalized using a centered log-ratio (clr) transform using scikit-bio (0.5.8). Before clr transformation, imputation of zero values was performed with geometric Bayesian multiplicative replacement as implemented in R package zCompositions (1.4.0-1). Integration of metagenomic and transcriptomic data was performed with sparse partial least squares (sPLS) by selecting at most 30 variables from each contrast with mixomics (6.20.0). Pearson correlations of features were visualized as clustered image maps (heatmaps) using the complete linkage method as an agglomerative clustering algorithm and using Euclidean distances as a distance measure for both types of features. Feature stability was evaluated with 10-fold cross-validation and 100 repeats. Predictors with variable importance in projection value (VIP) >1 were considered relevant, while the frequency of occurrence >0.8 was accepted as stable.

For the calculation of sIgA concentrations, the interpolation was done in GraphPad Prism by using the sigmoidal 4PL model. The p-values for the comparison of HOMA-IR, body mass, glucose, insulin, and sIgA levels were calculated using the Kruskal–Wallis test, followed by Dunn’s test for pairwise comparisons with the Bonferroni adjustment, and visualized with the ggplot2 package implemented in R.

In total, 24 animals were included in the study, with six animals representing each of the four treatment groups and equal gender distribution in each group (three males and three females). Four different contrasts were used in the statistical analysis of phenotypic as well as RNA-Seq data considering both the diet and metformin intervention (HFD vs. CD, HFDmetf vs. CDmetf, HFDmetf vs. HFD, CDmetf vs. CD); nevertheless, the greatest emphasis was placed on the contrasts reflecting the effect of metformin at the state of obesity and insulin resistance as well as the control phenotype (HFDmetf vs. HFD and CDmetf vs. CD).

To evaluate the intervention-induced phenotypic alterations, body mass and fecal sIgA levels were detected in addition to the HOMA-IR index, which was calculated for each animal based on fasting plasma glucose and insulin levels at the end of the experiment. Significantly higher body mass (HFD: median = 54.06 g, IQR = 52.38–54.26 g; HFDmetf: median = 50.96 g, IQR = 49.69–54.39 g; CD: median = 32.96 g, IQR = 29.56–39.60 g; CDmetf: median = 31.46 g, IQR = 30.86–35.69 g), insulin levels (HFD: median = 1,358.33 pmol/L, IQR = 652.05–1,918.17 pmol/L; HFDmetf: median = 1,072.35 pmol/L, IQR = 513.20–1,682.04 pmol/L; CD: median = 134.73 pmol/L, IQR = 78.65–174.68 pmol/L; CDmetf: median = 116.72 pmol/L, IQR = 78.83–134.70 pmol/L), and also HOMA-IR indexes (HFD: median = 9.39, IQR = 7.06–14.06; HFDmetf: median = 8.83, IQR = 5.06–15.07; CD: median = 1.21, IQR = 0.61–1.89; CDmetf: median = 0.82, IQR = 0.51–0.93) were detected in HFD-fed animals compared with CD-fed animals at the end of the experiment in all of the comparisons tested. The difference in blood glucose levels (HFD: median = 13.41 mmol/L, IQR = 9.87–15.37 mmol/L; HFDmetf: median = 11.47 mmol/L, IQR = 10.97–13.28 mmol/L; CD: median = 11.34 mmol/L, IQR = 11.06–11.95 mmol/L; CDmetf: median = 9.33 mmol/L, IQR = 9.27–9.67 mmol/L) did not reach statistical significance after adjusting the p-value. Similarly, although not reaching statistical significance, the body mass, HOMA-IR, glucose, and insulin levels tend to show lower levels in mice receiving metformin treatment against the untreated animals irrespective of the diet applied ( Figures 2A–C , Supplementary Table 5 ).

According to the previously reported property of metformin to elevate sIgA levels in human stool, additional ELISA-based measurements of fecal sIgA were performed after metformin intervention, showing notable though not significant variability in different treatment arms (HFD: median = 112.25 μg/mL, IQR = 63.01–176.24 μg/mL; HFDmetf: median = 48.92 μg/mL, IQR = 42.74–67.02 μg/mL; CD: median = 98.30 μg/mL, IQR = 59.38–138.35 μg/mL; CDmetf: median = 64.03 μg/mL, IQR = 55.33–107.92 μg/mL) ( Figure 2D , Supplementary Table 5 ). The data describing phenotypic variation among animals of different sexes are provided in Supplementary Table 5 and Supplementary Figure 7 ).

The RNA-Seq approach was applied to comprehensively identify transcripts responding to the diet and/or metformin treatment in samples obtained from the distal part of the small intestine of the 24 animals. Next-generation sequencing provided a median of 18.1 million sequencing reads per sample (IQR = 68.4), where 73.9% (IQR = 20.4) of the reads were mapped to the reference genome, ensuring the identification of 55,536 unique mouse genes in total.

To control for the effects of the diet, at first, the transcriptome profiles of animals receiving diets with distinct fat content (HFD vs. CD) and no metformin were compared, revealing an HFD-induced differential expression of 585 genes (FDR < 0.05). The majority of DEGs (395 genes) were uniquely differentially expressed in the HFD vs. CD contrast, including the ones functionally associated with lipid metabolism, such as Bco1 (logFC = −3.08; FDR = 1.69E−03) and Acaa1b (logFC = 1.56; FDR = 2.07E−02) ( Table 1 ). In total, 253 genes were downregulated, and 332 genes appeared to be upregulated due to increased dietary fat content solely ( Supplementary Table 1 , Supplementary Figures 1, 2 ). The comparison of the intestinal gene expression levels obtained from HFD-fed animals against CD-fed animals both receiving metformin indicated the modulatory effect of metformin on the diet-related intestinal transcriptome profiles exerted as the differential expression of 580 genes (333 downregulated and 247 upregulated). Out of those, 103 genes showed differential expression also in the HFD vs. CD contrast, while for the remaining majority (476 DEGs unique to the HFDmetf vs. CDmetf contrast), the diet-related alterations in expression levels were provoked directly by the presence of metformin ( Figure 3 , Supplementary Table 3 , Supplementary Figures 3, 4 ).

The KEGG pathway analysis of the list of identified DEGs in the HFD vs. CD comparison showed significant enrichment of 29 pathways, including butanoate metabolism among the uniquely enriched pathways (not represented in other contrasts) ( Table 2 ), while the GO analysis revealed the involvement of cellular response to the fatty acid and lipid metabolism exhibited by statistically significant enrichment of 197 GO terms in total. Additionally, 112 significantly enriched GO terms and 18 KEGG pathways were discovered by the functional annotation analysis of the DEGs identified in the comparison of both diet types together with metformin intervention (HFDmetf vs. CDmetf). Significant enrichment of 11 unique immunity-related signaling pathways was found representing the modulatory effects of metformin ( Table 2 , Supplementary Tables 1, 3 ).

We considered the HFDmetf vs. HFD as the most representable contrast reflecting the host transcriptomic responses to metformin treatment in the presence of obesity and insulin resistance. The metformin-related comparison revealed apparent differences in the expression of 46 genes (23 genes downregulated and 23 genes upregulated in the intestinal tissue of metformin-treated mice), with the most prominent upregulation of the ighv1-7 gene (logFC = 5.76; FDR = 3.64E−02) showing differential expression also in the contrast of CDmetf vs. CD (logFC = 5.09; FDR = 7.13E−03) ( Figures 4A, B ). The KEGG pathway analysis of 46 DEGs revealed the enrichment of three pathways: primary immunodeficiency, intestinal immune network for IgA production contrasts, and cytokine–cytokine receptor interaction ( Table 2 ); nevertheless, none of these reached statistical significance, and all three pathways were targeted also in other contrasts. Eight significantly enriched GO terms were identified in total, all falling within the ontology group of biological processes, and the majority of them related to immune functions, e.g., adaptive immune response ( Supplementary Table 2 , Figure 5 ).

Similarly, the comparison of the intestinal transcriptome of CD-fed animals receiving metformin with equally fed animals not treated with metformin was done resulting in 271 DEGs (56 genes downregulated and 215 upregulated in metformin-treated animals), with 91 of the genes appearing exclusively in the particular contrast (CDmetf vs. CD), including Ighv1-81 (logFC = 2.52; FDR = 7.13E−03) and Ighv14-3 (logFC = 2.76; FDR = 8.09E−03) ( Figure 2 , Supplementary Table 4 , Supplementary Figures 5, 6 ). The comparison of CDmetf vs. CD revealed 15 immunity-related KEGG pathways enriched, including leukocyte transendothelial migration, NF-kappa B signaling pathway, and the three pathways that were affected by metformin treatment in HFD-fed animals. In total, 98 GO terms were significantly enriched due to the administration of metformin in CD-fed animals with only two terms (adaptive immune response and positive regulation of interferon-gamma production) overlapping with the results of the HFD-fed animal analysis above ( Supplementary Table 4 ).

The compositional data of the animal gut microbiome were obtained at the end of the metformin intervention from the same animals and used for the integration with expression levels of DEGs identified in both contrasts reflecting the transcriptional effects of metformin (HFDmetf vs. HFD and CDmetf vs. CD). The correlation analysis of both datasets provided 31 metagenomics features of the highest relevance (VIP > 1) involving two genes from the comparison of HFD-fed animals and 16 genes derived from the contrast of CD-fed animals ( Figure 6 , Supplementary Tables 6, 7 ).