Metformin strongly affects transcriptome of peripheral blood cells in healthy individuals

Summary

Metformin is a commonly used antihyperglycaemic agent for the treatment of type 2 diabetes mellitus. Nevertheless, the exact mechanisms of action, underlying the various therapeutic effects of metformin, remain elusive. The goal of this study was to evaluate the alterations in longitudinal whole-blood transcriptome profiles of healthy individuals after a one-week metformin intervention in order to identify the novel molecular targets and further prompt the discovery of predictive biomarkers of metformin response. Next generation sequencing-based transcriptome analysis revealed metformin-induced differential expression of genes involved in intestinal immune network for IgA production and cytokine-cytokine receptor interaction pathways. Significantly elevated faecal sIgA levels during administration of metformin, and its correlation with the expression of genes associated with immune response (CXCR4, HLA-DQA1, MAP3K14, TNFRSF21, CCL4, ACVR1B, PF4, EPOR, CXCL8) supports a novel hypothesis of strong association between metformin and intestinal immune system, and for the first time provide evidence for altered RNA expression as a contributing mechanism of metformin’s action. In addition to universal effects, 4 clusters of functionally related genes with a subject-specific differential expression were distinguished, including genes relevant to insulin production (HNF1B, HNF1A, HNF4A, GCK, INS, NEUROD1, PAX4, PDX1, ABCC8, KCNJ11) and cholesterol homeostasis (APOB, LDLR, PCSK9). This inter-individual variation of the metformin effect on the transcriptional regulation goes in line with well-known variability of the therapeutic response to the drug.

Keywords:

Gene expression – Immune response – insulin – Gastrointestinal tract – Transcriptome analysis – Gene ontologies – RNA sequencing – CD coreceptors

Introduction

Metformin is the first-line antidiabetic agent used in pharmacotherapy of type 2 diabetes to improve glucose homeostasis[1]. Various additional therapeutic benefits beyond its antihyperglycaemic action have been highlighted lately, justifying the pleiotropic effect of the drug. In patients with type 2 diabetes metformin therapy is associated with reduced cardiovascular morbidity[2]. In addition, metformin exposure has a protective role against tumorigenesis in various types of cancers and it is proved to be beneficial in the preventive oncology regardless of the diabetic state[3, 4]. Furthermore, metformin therapy is often recommended for women with polycystic ovary syndrome to improve insulin sensitivity, facilitate menstrual regularity, induce ovulation, reduce circulating androgen levels and body weight[5, 6]. Likewise, metformin application in the treatment of such neurodegenerative disorders as Alzheimer’s and Parkinson’s diseases is currently under consideration[7–9].

Several possible molecular mechanisms of metformin action have been proposed, including inhibition of mitochondrial respiratory-chain complex 1, the reduction of cyclic adenosine monophosphate (AMP) levels, activation of AMP-activated protein kinase (AMPK) and recently described interaction with gut microbiota[10–13]. Although they partially explain major beneficial effects of the drug, exact mechanisms of metformin action remain unclear even after 60 years since its first clinical use.

Despite the widespread application of metformin, approximately 30% of type 2 diabetes patients using the drug are failing to achieve the adequate glycemic control[14]. Moreover, 20%-30% of type 2 diabetes patients suffer metformin-associated gastrointestinal adverse events and about 5% discontinue the therapy because of severe intolerance[15, 16]. Heritability of the glycaemic response to metformin has been suggested to depend on many allelic variants with small to moderate effects[17]. Contribution of inheritance to variation in metformin response has gained great interest in the past years, encouraging numerous targeted studies investigating genes coding for organic cation transporters OCT1, OCT2, OCT3 and multidrug and toxin extrusion proteins MATE1 and MATE2-K[18–22]. Moreover, Genome-Wide Association Studies have revealed multiple genetic variations within ATM, PRPF31, CPA6, and STAT genes associated with metformin response[23–25]. However, genetic alterations explain only a small proportion of the heterogeneous response to metformin therapy, therefore omics-based investigation of the pleiotropic mechanism of the drug is needed to promote the development of biomarkers for therapeutic efficacy[26].

Previous studies have demonstrated metformin-mediated changes at the transcriptome level in various animal tissues, nevertheless studies of metformin-related transcriptome profiles in humans are lacking. For instance, recent study by Guo et al. discovered metformin-induced alterations of the coding transcriptome profile and non-coding RNAs in the liver of high-fat diet induced mouse model of non-alcoholic fatty liver disease[27]. Likewise, microarray analysis of mice liver and muscle tissues revealed the ability of metformin to mimic the calorie restriction-like transcriptome[28]. Furthermore, a distinct gene expression profile related to cardiovascular benefits of metformin, was observed in a rat model of obesity and insulin resistance[29]. Meanwhile, cell culture studies of adrenal H295R cell and MCF7 breast cancer cell transcriptome have revealed the association of metformin with complex cellular processes related with energy metabolism, steroidogenesis and the immune system as well as glycolysis and cancer-related pathways[30, 31].

To identify the genes targeted by metformin, we performed the whole-transcriptome analysis with total RNA sequencing on whole-blood samples, obtained from the healthy individuals undergoing a seven-day course of metformin. To the best of our knowledge, this is the first study providing information about the immediate effects of metformin administration on global gene expression in healthy individuals.

Materials and methods

Study design

The study enrolled 25 healthy European descent volunteers with no history of chronic disease, meeting exclusion/inclusion criteria (S4 Table) set within the framework of the ongoing clinical trial ‘Pharmacodynamics of antidiabetic drug metformin’ (50 individuals to be included in total), protocol number MIKROMET16001, registration number of EU Clinical Trials Register: 2016-001092-74 (www.clinicaltrialsregister.eu) (Fig 1). Participants received twice-daily oral 850mg dose of metformin hydrochloride (Metforal 850mg, Berlin-Chemie AG) for 7 days. Medication adherence was reported by each participant at the end of the active period of the study. Fasting blood tests (e.g. measures of ALT, plasma glucose, creatinine levels) were performed in certified clinical laboratory 1–3 days before metformin administration in order to evaluate general hematological and biochemical parameters, and eligibility of volunteers (Table 1). RNA for transcriptome analysis was isolated from the whole blood samples collected in Latvian Biomedical Research and Study Centre at three time-points: (1) before administration of metformin (M0, morning, fasting state), 10 hours after the first metformin intake, but before the second dose (M10h, evening) and after 7 days long metformin course (M7d, morning, fasting state). The third blood sample was not collected from one out of 25 study subjects, due to the discontinuation of metformin treatment, nevertheless the rest of the blood samples collected from the particular subject were included in the data analysis. Longitudinal study design of the open-label trial was chosen as the most suitable method for global gene expression analysis with high inter-individual variability expected.

All individuals were concurrently involved in an ongoing analysis of gut microbiome and DNA methylation profiles according to the study protocol (S1 Text). The primary endpoint of the study was significantly shifted composition of the gut microbiome after administration of metformin. The secondary endpoint of the study was alterations in DNA methylation profiles and mRNA levels following metformin use.

Written informed consent was obtained from every participant and the study protocol was approved by the Central Medical Ethics committee of Latvia (1/16-05-12) and the State Agency of Medicines of the Republic of Latvia (17–1723). The research was conducted in accordance with the The Code of Ethics of the World Medical Association (Declaration of Helsinki) and International Conference on Harmonisation E6 (R2) Guidelines for Good Clinical Practice. All participants were included in the Genome Database of the Latvian Population[32].

RNA sample preparation and next generation sequencing

Venous blood samples (n = 74) were collected in Tempus Blood RNA Tubes, followed by total RNA isolation with PerfectPure RNA Blood Kit (5Prime GmbH, Germany), according to the manufacturer's instructions. The integrity of the extracted RNA was evaluated by RNA integrity number (RIN) within Agilent 2100 Bioanalyzer system (Agilent, USA). Ribosomal RNA depletion was done with Low Input RiboMinus Eukaryote System v2 (Thermo Fisher Scientific, USA) by processing 500ng of total RNA from each sample. For cDNA library preparation we used Ion Total RNA-Seq Kit v2 (Thermo Fisher Scientific, USA), sequencing was done on Ion Proton System and Ion PI Chip (Thermo Fisher Scientific, USA), following the manufacturer’s instructions. Since shot-gun RNA sequencing is considered to be the most accurate and desirable method for quantification of expression of individual transcripts and genes, additional methods for technical validation were not applied in this study[33].

Stool sample collection and detection of secretory IgA by ELISA

Within the framework of the clinical trial two aliquots of stool samples were collected from each study participant at three time points: before administration of metformin (M0), 24 hours after the first dose (M24h) and 7 days after the first intake of metformin (M7d), except for two participants for whom the third stool sample was not available (n = 73 in total). The samples were transferred at −80°C within 24 hours since the collection. The concentration of secretory immunoglobulin A (sIgA) in 100mg of each stool sample was determined by ImmuChrom ELISA Kit (ImmuChrom GmbH, Germany), according to the manufacturer’s instructions, the absorbance was read at 450nm and 620nm as the reference wavelength.

Bioinformatic analysis

The sequencing reads were trimmed with Trimmomatic 0.36 using window size 5 and quality threshold 10. After trimming reads had to have a minimum length of 30bp and average quality of 10 to be included in subsequent analyzes. Sequencing reads were mapped against human reference genome GRCh38 release 90 and read counts were calculated with STAR 2.5.3a. The obtained read counts were normalized using trimmed mean normalization implemented in Bioconductor package edgeR in R. Differentially expressed genes (DEGs) were estimated with two different methods. At first, likelihood ratio test with added observation weights was used to reduce the influence of outliers and to obtain a list of DEGs, where the counts per million (CPM) value had to be 1 or more in at least 24 samples for the gene to be included in the analysis (edgeR-robust)[34]. In order to account for subject-specific expression the quasi-likelihood F-test without any prior gene filtering was performed (edgeR-sensitive). Multiple testing correction was implemented using Benjamini-Hochberg procedure, significant DEGs were determined using false discovery rate (FDR) < 0.05 cutoff[35].

Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were adopted as the functional terms. GO and KEGG pathway enrichment analysis were performed with R package Goseq (1.30.0)[36]. KEGG pathways and GO terms with FDR < 0.05 were considered statistically significant.

Heat map was constructed with Matplotlib2 and SciPy3. Hierarchical clustering with average linkage method implemented in SciPy was used to cluster DEGs for contrasts M10h vs M0 and M7d vs M0 by their differences in read counts per million[37, 38].

Four-parameter algorithm in GraphPad Prism 8 was used for the calculation of sIgA concentrations in faecal samples and Wilcoxon signed rank test in R was applied for the evaluation of fluctuations in faecal sIgA levels among different time points. For the correlation analysis between sIgA levels and CPM values obtained in edgeR-robust analysis method, Spearman’s rank correlation test was performed in R, where significance level was set at 0.05.

Results

Differential global gene expression induced by administration of metformin

In order to reveal the target genes and pathways affected by metformin, we performed a transcriptome analysis in whole-blood samples of 25 healthy volunteers receiving metformin for one week (Table 1). Venous blood samples were obtained at three consecutive time-points, hereinafter referred to as M0 (before administration of metformin), M10h (10 hours after the first metformin intake/before the second dose) and M7d (after 7 days long metformin course) in order to observe both, acute and sustained effects of metformin. RNA-Sequencing produced an average of 24.6 ± 8.9 million reads per sample, 83.5% of the reads were mapped to the reference genome.

We focused on DEGs at three contrasts comparing the samples collected at previously defined time-points: M10h vs M0, M7d vs M0 and M7d vs M10h. In total 561 unique genes showed significantly different expression levels among the analyzed contrasts (Fig 2A, 2B, 2C and 2D). The majority, 479 of DEGs appeared in the contrast M10h vs M0 (364 downregulated and 115 upregulated). Comparison of M7d vs M0 resulted in 82 DEGs (61 downregulated, 21 upregulated) and 120 DEGs were identified in the contrast M7d vs M10h with higher proportion of upregulated genes (17 downregulated, 103 upregulated) (Fig 3). The overlap of the two main contrasts (M10h vs M0 and M7d vs M0) consisted of 44 DEGs, including UBE2O, PHOSPHO1, MKRN1, possessing consistent metformin-evoked alterations in expression levels for 7 days. The complete lists of obtained DEGs are provided in S1 Table.

Functional annotation of the identified DEGs

In order to gain insights into the molecular mechanisms underlying the short-term effects of metformin the KEGG pathway enrichment analysis was performed. Three lists of DEGs corresponding to each contrast were submitted in Goseq package of R. DEGs obtained from the contrasts M10h vs M0 and M7d vs M10h were assigned to pathways related with immune responses, while list of DEGs obtained by comparing M7d samples with M0 samples showed no pathways enriched (Table 2). In order to gain a greater understanding of the biological implications of the obtained DEGs Gene Ontology (GO) enrichment analysis was performed, see S2 Table and S3 Table for the complete list of GO terms and KEGG pathways corresponding to DEGs obtained in all of the contrasts analyzed.

Elevated faecal sIgA levels during the administration of metformin

Considering the enrichment of pathways related to intestinal immune responses, sIgA concentration was determined in stool samples collected from the study participants at three consecutive time points, analogous to the course of blood sample collection: before administration of metformin (M0; median sIgA concentration = 7969.93μg/ml; IQR = 15587.55), within 24 hours after the first dose (M24h; median sIgA concentration = 6935.29μg/ml; IQR = 28953.46) and 7 days after the first intake of metformin (M7d; median sIgA concentration = 21207.64μg/ml; IQR = 36642.19), revealing significantly increased faecal sIgA levels comparing M7d samples vs M0 samples and M7d samples vs M24h samples (Fig 4).

The possible implication of specific immunity-related genes in intestinal IgA production in response to metformin administration was evaluated by Spearman’s rank correlation test between faecal sIgA concentration and expression level of DEGs involved in Intestinal immune network for IgA production and Cytokine-cytokine receptor interaction pathways, revealing 9 significant correlations (Table 3).

Subject-specific effects of metformin

In order to clarify subject-specific metformin effects on the transcriptome profile, we performed an additional data analysis (quasi-likelihood F-test without any prior gene filtering hereinafter referred to as edgeR-sensitive), which allowed us to identify 437 unique DEGs among all of the contrasts (S1 Table; S1A, S1B, S1C and S1D Fig). After a careful inspection of the generated heat map for the contrasts M10h vs M0 and M7d vs M0, we observed striking inter-individual variation in expression levels of most significant DEGs (S2 Fig). It was evident that the overall differences in expression levels for a number of gene clusters were influenced by extreme changes in expression observed only in few individuals from the study. Four such functional clusters of genes with similar expression profiles were distinguished after consideration of subject-specific effects. In each of the recognized clusters, genes sharing a common function were predominant, including upregulated genes, coding for small nuclear RNAs (snRNAs): SNORA20, SCARNA5, SNORA80E, SNORA3B, SCARNA22, SCARNA6, SNORD8, SNORA7B, SNORD9, SNORD83A, SNORA23, SNORD67, SNORD46, SNORA71B, SNORD71, SNORD11B, SNORD17, SNORD116-15, SNORD67; ribosomal genes and pseudogenes together with long non-coding RNAs (lncRNAs): Y_RNA, RNY4P10, RNA5-8SN1, RNA5-8SN5, FP671120.4, FP236383.2, FP671120.3, FP236383.3, RPL3P4, AC008038.1, AL122020.1, RPL13AP5, RPS2P46, RPL6P27, RPL13AP25, AL590867, RPS7P1, AP001324.1, AC034236.1, RPL37P2, AC004453.1, RPL21P16, AC007969.1, RPS23P8, RPL15P3, RPL10AP6; genes relevant to insulin production: HNF1B, HNF1A, HNF4A, GCK, INS, NEUROD1, PAX4, PDX1, ABCC8, KCNJ11 and downregulated genes involved in cholesterol homeostasis: APOB, LDLR, PCSK9.

Discussion

This longitudinal study of metformin administration in healthy subjects demonstrated that metformin strongly affects the gene expression profiles in blood cells as estimated by RNA sequencing. To our knowledge this is the first study accessing the metformin-mediated changes in RNA expression in vivo in humans. Moreover, the significant transcriptomic changes were observed even after 10 hours as a result of single metformin dose, indicating the pronounced and immediate influence on the cell functions. The most striking finding was the strong evidence for metformin-induced enrichment of immunity-related pathways resulting in elevated faecal IgA levels.

Multiple DEGs identified in this study represent the main functional groups associated with previously described therapeutic targets of metformin. For instance, downregulation of gene coding for ubiquitin conjugating enzyme E2 O (UBE2O) impairs the tumorigenesis, moreover combined treatment of UBE2O inhibitors and AMPK agonists, such as metformin, has been suggested as promising treatment strategy for cancer already before[39]. Similarly, suppression of MKRN1 (makorin ring finger protein 1) activates AMPK, resulting in increased glucose consumption and reduced lipid accumulation, therefore MKRN1-mediated regulation of AMPK activity has been already considered as an attractive therapeutic approach for the treatment of metabolic disorders[40]. Finally, DNA methylation at the PHOSPHO1 (phosphoethanolamine/phosphocholine phosphatase 1) locus in blood cells has been previously linked to decreased type 2 diabetes mellitus risk, which goes in line with the well-known antidiabetic activity of metformin and metformin-evoked downregulation of PHOSPHO1 observed in this study[41].

A group of discovered DEGs correspond to the pathways involved in immune response or regulation of inflammation. The pathway enrichment analysis showed comprehensively decreased expression levels of genes related to immune responses, which confirms the anti-inflammatory effect as a universal property of metformin. Here we convincingly support the previously described ability of metformin to suppress inflammatory cytokines and their receptors irrespective of the diabetes status, specifying the occurrence of this process at the level of mRNA[42]. The downregulation of CXCL8 and CXCR4, coding for interleukin-8 and chemokine receptor type 4 respectively, has been previously attributed to the anticancer action of metformin in neoplastic cells, therefore the results of the present study describe the particular therapeutic effects as universal and characteristic also to the normal blood cells[43, 44].

Furthermore, RNA-Seq data revealed significant enrichment in the pathway of intestinal immune network for IgA production, which was further confirmed by elevated faecal sIgA in response to the metformin intervention. The observed correlation between faecal sIgA concentration and immunity-related gene expression levels points at transcriptional shift as a constitutive component of the intestinal immunity-related effects of metformin. IgA is responsible for bacteria-host interaction and is massively produced by mucosa in case of bacterial colonization; moreover, selective IgA deficiency is associated with mild intestinal dysbiosis and shifts in the microbial composition[45, 46]. Since the influence of metformin on gut microbiome is currently extensively studied[47], the potential implication of the intestinal immune network for IgA pathway in the metformin-microbiome interactions regarding the intestinal side effects should be considered.

In addition, we noticed subject-specific differences at gene expression levels, as number of individuals involved in study displayed overexpression or repression of specific functionally related gene sets. Thus for example insulin coding gene (INS) showed an apparent gain of the mRNA expression levels (from 0.52 CPM before the treatment to 282.97 CPM after one week of metformin administration) in a single individual. Based on the hierarchical clustering and functionality of these genes we categorized 4 gene clusters: ribosomal genes and their pseudogenes, small nuclear RNAs, genes relevant to insulin production and cholesterol homeostasis. Interestingly, the changes in RNA expression of those genes are highly subject-specific with strongly altered expression in only one or two participants of the study.

Metformin-induced overexpression of Maturity onset diabetes of the young (MODY)—related genes (INS, PDX1, PAX4, HNF4A, HNF1A, HNF1B, NEUROD1, GCK) coding for transcription factors and regulators of β-cell function was observed in both main comparisons (M10h vs M0, M7d vs M0), but not in the contrast M7d vs M10h, suggesting that these alterations are likely to be associated with metformin intervention rather than discrepancies between fasting and feasting states at the time points of blood collection. Metformin-induced differential expression of MODY genes has been previously reported in the liver of spontaneously hypertensive rats, overlapping several homologues of human genes[48]. Moreover, subject-specific activation of the MODY genes may be the reason why the metformin ability to induce insulin secretion has not been observed before in human trials. One may speculate that metformin exerts insulin secretagogue ability only in subgroup of metabolically compromised individuals, however, to prove this additional research in patients with metabolic syndrome and diabetes is needed.

Very similar to our detection of the MODY cluster our study revealed a considerable downregulation of the genes coding for apolipoprotein-B (APOB), low-density lipoprotein receptor (LDLR) and proprotein convertase subtilisin/kexin type 9 (PSCK9) in one person from the study group. All of these genes are previously associated with cholesterol homeostasis and phenotype of familial hypercholesterinaemia[49–51]. To date, several reports have described a dose-dependent cholesterol lowering effect of metformin[52, 53]. Taken together, it would be reasonable to argue that the downregulation of LDLR, PCSK9, APOB might serve as potential mechanism of action, underlying beneficial cardiovascular properties of metformin, yet in a case-specific manner.

The following are some study limitations. These limitations include small sample size and only one week of the intervention time due to the fact that study was performed on healthy individuals. Another limitation is that this is an exploratory study without a placebo control arm that would be needed to draw definitive conclusions on the causality of observed transcriptional changes. Also a larger study group of diabetic subjects and longer observation time would provide information on the factors that may explain subject-specific differences in expression levels and relation of these differences to the treatment response including glucose control.

These limitations are compensated by the longitudinal design of the study in which the first sample from an individual was the control for further samples. We believe that this design and short-term observation should have minimized false associations and conclusions arising from unaccounted factors playing important role in human studies, meanwhile making reasonable the interpretation of observed inter-individual variability of gene expression profiles.

In conclusion, we were able to provide, for the first time, direct evidence of the effects of metformin on the immediate and strong transcriptome changes in whole-blood samples. Our results have pinpointed some important targets that need further investigation. First, the ability of metformin to induce extensive immune responses may be executed at the level of transcription and serve as the basis of common therapeutic effect of metformin. Second, the induction of IgA pathway may explain the widely discussed interaction of metformin with the gut microbiome. Third, the subject-specific response may explain the large percent of unresponsiveness to the metformin therapy. Altogether these results may serve the ground for development of expression based biomarker sets to predict and/or monitor the treatment outcomes.

Supporting information

S1 Table [xlsx]
List of differentially expressed genes.

S2 Table [xlsx]
Gene ontology analysis.

S3 Table [xlsx]
KEGG pathway enrichment.

S4 Table [pdf]
Inclusion and exclusion criteria of the clinical trial ‘Pharmacodynamics of antidiabetic drug metformin’.

S5 Table [pdf]
CONSORT Checklist.

S1 Text [pdf]
Protocol of the clinical trial ‘Pharmacodynamics of antidiabetic drug metformin’.

S1 Fig [a]
Metformin-induced alterations in gene expression profiles, considering subject-specific effects.

S2 Fig [tiff]
Heat map and hierarchical clustering of 404 significant DEGs in the contrasts M10h vs M0 and M7d vs M0, accounting for subject-specific effects.