Whole-blood transcriptome profiling reveals signatures of metformin and its therapeutic response

Metformin, a biguanide agent, is the first-line treatment for type 2 diabetes mellitus due to its glucose-lowering effect. Despite its wide application in the treatment of multiple health conditions, the glycemic response to metformin is highly variable, emphasizing the need for reliable biomarkers. We chose the RNA-Seq-based comparative transcriptomics approach to evaluate the systemic effect of metformin and highlight potential predictive biomarkers of metformin response in drug-naïve volunteers with type 2 diabetes in vivo. The longitudinal blood-derived transcriptome analysis revealed metformin-induced differential expression of novel and previously described genes involved in cholesterol homeostasis (SLC46A1 and LRP1), cancer development (CYP1B1, STAB1, CCR2, TMEM176B), and immune responses (CD14, CD163) after administration of metformin for three months. We demonstrate for the first time a transcriptome-based molecular discrimination between metformin responders (delta HbA1c ≥ 1% or 12.6 mmol/mol) and non-responders (delta HbA1c < 1% or 12.6 mmol/mol), that is determined by expression levels of 56 genes, explaining 13.9% of the variance in the therapeutic efficacy of the drug. Moreover, we found a significant upregulation of IRS2 gene (log2FC 0.89) in responders compared to non-responders before the use of metformin. Finally, we provide evidence for the mitochondrial respiratory complex I as one of the factors related to the high variability of the therapeutic response to metformin in patients with type 2 diabetes mellitus.


Introduction
Diabetes mellitus is a chronic disease affecting approximately 463 million people worldwide, which is nearly 9.3% of the global population [1]. Type 2 diabetes mellitus (T2DM) is the most common type of diabetes accounting for approximately 90% of all cases. The persistent hyperglycemia and insulin resistance, a characteristic of T2DM patients, is associated with an increased risk of serious microvascular and macrovascular complications, including nephropathy, retinopathy, neuropathy, myocardial infarction, and stroke, which may be reduced by early initiation of antidiabetic therapy [2][3][4]. Metformin is the first-line medication for treating hyperglycemia in T2DM with beneficial effects in the treatment of multiple nondiabetes related conditions, such as polycystic ovary syndrome, cancer, and neurodegenerative disorders [5][6][7]. Despite the pleiotropic effects of the drug, the variable efficacy and gastrointestinal side-effects observed cause a significant non-compliance and discontinuation of the therapy, justifying a need for studies exploring molecular mechanisms of metformin action, and biomarkers predicting both treatment response and tolerance of the drug [8].
The mechanism of metformin action is generally considered to involve modulation of the activity of mitochondrial complex I, activation of 5 0 AMP-activated protein kinase (AMPK)dependent mechanisms, and increased AMP concentrations, though some controversy remains since multiple studies are providing evidence for other indirect mechanisms, such as the significant contribution of the gut microbiome underlying the glucose-lowering effect of the drug [9][10][11].
RNA sequencing (RNA-Seq) is the state-of-the-art approach that may be used to profile drug response and efficacy biomarkers [12,13]. So far, transcriptome datasets obtained from cell cultures and tissue samples of animal models are extensively used in studies describing molecular mechanisms of metformin concerning various conditions, nevertheless, longitudinal data of in vivo studies in humans are still lacking. RNA-Seq has revealed various novel effects and therapeutic targets of metformin, such as enrichment of the transcriptional regulator forkhead box O3a (FOXO3a) in primary human fibroblasts [14], upregulation of activating transcription factor 3 (ATF3) in primary human hepatocytes [15], downregulation of cell division control protein 42 homolog (CDC42) in breast cancer cells [16], upregulation of krüppel-like factor 4 (KLF4) resulting in suppressed endothelial dysfunction [17], and even modulated alternative splicing in embryonic stem cells [18]. Moreover, multiple animal-based studies have reported metformin-specific signatures in gene expression profiles of rat arteries and mice epididymal fat, liver and muscle tissue, nevertheless, they still do not explain many beneficial effects of the drug [19][20][21].
Although the relatively high proportion (>30%) of patients failing to achieve glycemic control during metformin therapy has been partially explained by the contribution of genetic inheritance (allelic variants of organic cation transporters OCT1, OCT2, etc.) [22][23][24], recent studies report that heterogeneity of metformin response may be both patient and cell type-specific, suggesting the presence of yet unknown, non-genetic and selective manifestations of the drug [25,26]. The main objective of the study was to assess the systemic effect of metformin in T2DM patients and reveal potential biomarkers for accurate prediction of its therapeutic efficacy. We have previously reported direct evidence of the effects of metformin on the immediate and strong transcriptome changes in whole-blood samples of healthy subjects, though we considered the diabetic state as a significant confounding factor, therefore the study was continued in a well-characterized, prospective T2DM patient cohort with prolonged observational time, providing much wider applicability of the study results [27]. We believe that our strategy will promote the development of biomarker-based approaches for monitoring treatment outcomes and early identification of metformin responders, moving towards precision medicine.

Study design
The study was conducted within the framework of the ongoing observational, prospective and longitudinal study OPTIMED, which has been implemented since 2010 in collaboration with endocrinologists and general practitioners from the leading health care centers in Latvia, ensuring recruitment of newly diagnosed drug-naïve patients with ICD-10 diagnosis code E11 and follow-up data collection. Written informed consent was obtained from every participant after full explanation of the purpose and nature of all procedures used before their inclusion in the study, and the study protocol was approved by Central Medical Ethics Committee of Latvia (No. 01-29.1/22) and Committee of Ethics in Pauls Stradins Clinical University Hospital (No.3000610 -18L). The research was conducted in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki amended in Fortaleza, Brazil, October 2013) and The Convention for the protection of Human Rights and Dignity of the Human Being with regard to the Application of Biology and Medicine: Convention on Human Rights and Biomedicine. Management of patient recruitment, collection of samples and associated clinical data was ensured by Latvian Biomedical Research and Study Centre's core facility Genome Centre and the Genome Database of the Latvian Population following their standard procedures [28].
In total, 17 patients of European descent fulfilling the following inclusion criteria were enrolled: (1a) newly diagnosed type 2 diabetes mellitus (ICD-10 code E11) requiring oral antidiabetic therapy, or (1b) previously diagnosed type 2 diabetes mellitus but no oral antidiabetic therapy or insulin has been used for the last three years, or (1c) newly diagnosed type 2 diabetes mellitus and intensive insulin therapy initiated in a hospital for acute glycemic normalization; (2) the patient is not currently involved and is not planning to enroll in clinical trials during the OPTIMED study; (3) the patient has attained 18 years of age; (4) the patient is not pregnant at the time of application; (5) the patient meets the criteria for the diagnosis of type 2 diabetes mellitus: (a) fasting blood glucose level � 7 mmol/l, (b) a blood glucose level � 11.1 mmol/l for a two-hour glucose tolerance test with 75 g intake. The exclusion criteria of the study were as follows: (1) the patient is receiving oral antidiabetic therapy on a regular basis or has received the therapy during the last three years; (2) the patient is receiving insulin therapy at the time of application; (3) the patient is pregnant. According to the observational study design, the randomization procedure was not performed and metformin monotherapy (medication with metformin hydrochloride as the only active ingredient) was prescribed for each study participant by an endocrinologist for at least three months regardless of the research objectives. The drug manufacturer and dosage of metformin (varied from 850 mg to 2000 mg per day) were chosen by endocrinologists based on clinical experience, patient's health status, and manifestations of the disease. Enrollment in the study did not affect the choice of treatment strategy made by endocrinologists. Blood tests (e.g. measures of ALT, creatinine levels, HbA1c) were performed in a certified clinical laboratory before the administration of metformin and after metformin therapy for three months to evaluate general hematological and biochemical parameters and eligibility of the subjects (Table 1). Blood samples for RNA-Seq were collected at the same time points, hereinafter referred to as M0 (before administration of metformin) and M3m (after three months long metformin course). Considering the high inter- individual variability expected in RNA-Seq data, longitudinally repeated measurements were chosen as the most suitable approach for global gene expression analysis. HbA1c measurements were made for all participants at two consecutive time points analogous to the blood collection for RNA extraction (M0 and M3m). Study participants were further stratified in two subgroups based on metformin response, which was defined according to metformin-induced alterations in glycated hemoglobin (HbA1c) levels comparing measurements made in the time points M0 and M3m: responders (delta HbA1c � 1% or 12.6 mmol/ mol), non-responders (delta HbA1c < 1% or 12.6 mmol/mol) [29].

Sample processing and RNA sequencing
For RNA isolation, 3 ml of whole blood were collected in Tempus™ Blood RNA Tubes (Thermo Fisher Scientific, USA) and further processed using Tempus™ Spin RNA Isolation Kit (Thermo Fisher Scientific, USA) according to manufacturer's instructions. The quantity and quality of extracted RNA and prepared libraries were determined by Qubit Fluorometer (Thermo Fisher Scientific, USA) and Agilent 2100 Bioanalyzer systems (Agilent, USA), respectively. The integrity of RNA was evaluated by RNA integrity number (RIN) within the Agilent 2100 Bioanalyzer system (Agilent, USA). For depletion of ribosomal RNA 500 ng of total RNA from each sample were processed using Low Input RiboMinus™ Eukaryote System v2 (Thermo Fisher Scientific, USA). Complementary DNA library preparation was performed with Ion Total RNA-Seq Kit v2 (Thermo Fisher Scientific, USA). Ion Proton™ System (Post-Light™ Ion Semiconductor Sequencing, Thermo Fisher Scientific, USA) and Ion PI™ Chip (Thermo Fisher Scientific, USA) was used for 200-base-read single-end sequencing, following the manufacturer's instructions. Since the shot-gun RNA-Seq is considered to be the most accurate and desirable method for the quantification of the individual transcript and gene expression, additional methods for technical validation were not applied in this study [30].

Data analysis
Trimmomatic 0.36 was used for read trimming applying window size 5 and quality threshold of 10. After trimming reads had to have a minimum length of 30 bp and an average quality of 10 to be included in subsequent analyses. Sequencing reads were mapped against human reference genome GRCh38 release 90 and per-gene read counts were calculated with STAR (v.2.5.3a.). STAR (v.2.5.3a.) outputs read quantification per gene while performing read mapping. The reads were quantified if they match only one gene. The obtained read counts were then normalized using trimmed mean normalization implemented in Bioconductor package edgeR in R (v.3.5.3). Differentially expressed genes (DEGs) were estimated using the Likelihood ratio test with added observation weights to reduce the influence of outliers, and sva (Surrogate Variable Analysis) package in R (v.3.5.3) was used for removing batch effects [31]. In order to evaluate metformin-induced alterations in the transcriptome profile (comparison of samples M3m vs M0), each sample was set as a factor considering longitudinal study design. When comparing responders against non-responders in each time point separately, sequencing run and baseline (M0) HbA1c levels were considered as covariates. FilterByExpr function was applied for gene filtering in edgeR, taking into account the sample library sizes [32]. Multiple testing correction was implemented using the Benjamini-Hochberg procedure and differential expression of the genes was determined using a false discovery rate (FDR) < 0.05 cutoff, regardless of the log 2 fold change of expression for each gene [33]. In order to identify key genes determining the metformin response, Partial least squares discriminant analysis (PLS-DA) was performed implemented in the mixOmics package of R (v.3.5.3). CPM values (obtained in edgeR and adjusted for the impact of sequencing run and baseline HbA1c levels) were used in PLS-DA. Key genes contributing to a separation of patients in both metformin response groups were identified by using a cutoff of variable importance of projection (VIP) score >1 obtained from PLS-DA [34]. Association between HbA1c levels and the log CPM expression values of each mitochondrial gene was performed with multiple linear regression using lm function in R (v.3.5.3). Sex and body mass index were included in the model to account for their potential confounding.
Genes showing the p-value <0.05 for differential expression were further used in the functional analysis. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were adopted as the functional terms. GO enrichment analysis was performed with R package Goseq (v.1.38.0), and KEGG pathway enrichment analysis was done using an online software Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8., the threshold value of enrichment was selected by a P-value <0.05 [35,36].
Heat maps were constructed with Matplotlib and SciPy. Hierarchical clustering with average linkage method implemented in SciPy was used for the clustering of genes according to their differences in CPM values [37,38]. Statistical analysis of anthropometric measures and biochemical data was performed in R (v3.5.3.) by applying the Wilcoxon rank-sum test and Pearson's chi-squared test with a p-value threshold < 0.05.

Identification of metformin-induced differential expression of genes in blood cells
We characterized the transcriptome profiles of whole-blood samples obtained from 17 drugnaïve T2D patients (characteristics of patients provided in Table 1) before any antidiabetic therapy (M0) and after three months of metformin monotherapy (M3m) by RNA-Seq technique to detect sustained transcriptional alterations in blood cells at the diabetic state due to the administration of the drug. The median of read counts per sample produced by RNA sequencing was 21.1 (IQR = 6.3) and 87.9% of the reads (median = 18.2; IQR = 5.1) were further mapped to the human reference genome.
Differential expression analysis performed by comparing transcriptome profiles of blood samples collected after metformin administration for three months against samples collected before the use of any antidiabetic therapy (M3m vs M0) revealed 28 DEGs (FDR < 0.05) from the pool of 9992 transcripts identified in total. Out of them, 20 genes were significantly downregulated and 8 genes showed significant up-regulation after administration of metformin (

Functional enrichment characteristic to metformin therapy
In order to describe the implementation of metformin-modulated transcriptome profiles in cell signaling pathways and core biologic functions, the KEGG pathway and GO enrichment analysis was performed. Over-representation of biological pathways (e.g. amino sugar and nucleotide sugar metabolism, antigen processing and presentation) and GO Terms (e.g. helper T cell chemotaxis and lipoprotein particle receptor activity) related to energy metabolism, immune responses and lipid metabolism were observed, see S1 and S2 Tables for the complete list of GO terms and enriched KEGG pathways.

Determining genes involved in differential metformin responsiveness
Considering the potential contribution of distinct molecular mechanisms mediating variable metformin response, all participants were stratified in two efficacy groups, according to metformin-induced changes in their HbA1c levels (responders: delta HbA1c �1% or 12.6 mmol/mol; non-responders: delta HbA1c <1% or 12.6 mmol/mol) ( Table 3).
Genes associated with variable metformin responsiveness were determined by differential gene expression analysis comparing the number of RNA-Seq reads from responders against non-responders in each time point of blood collection (M0 and M3m) separately. In total, 27 significant DEGs were identified contrasting responders against non-responders before the administration of metformin (M0). We observed a notable portion of downregulated genes coding for small nucleolar RNAs and upregulated insulin receptor substrate 2 (IRS2) gene in responders, revealing distinctive features between response groups already before the use of any antidiabetic therapy (Fig 2C), (Table 4).
PLS-DA was applied to understand whether patients could be clustered based on their gene expression profiles before the administration of metformin. PLS-DA revealed pronounced discrimination of samples, based on the obtained transcriptome profiles, explaining 13.9% of the variance in total (11% of variance explained by latent variable 1 and 2.9% of variance explained by latent variable 2) (Fig 3B). In total, 56 discriminatory genes showing the strongest separation of different metformin response groups were identified, based on variable importance for projection (VIP>1) value in two latent variables (components), generated by PLS-DA ( Fig  3A). See S5 Table for the full list and VIP values of discriminatory genes.
Comparison of gene expression profiles between metformin responders and non-responders after the use of metformin for three months (M3m) showed differential expression of 15 genes (12 up-regulated, 3 down-regulated), including 5 mitochondrial genes, which may be associated with mechanisms underlying variable efficacy of the drug (Table 5) (Fig 2D).  Log 2 FC, log 2 fold change; FDR, false discovery rate. † Genes showing significant differential expression due to the administration of metformin also in healthy individuals [27].

PLOS ONE
Moreover, a multiple linear regression model revealed a significant association between the expression of two out of five mitochondrial genes and HbA1c levels after metformin therapy for three months (MT-ND4, p-value = 0.047; MT-ND4L, p-value = 0.032) (S7 Table).

Biological functions contributing to metformin responsiveness
GO analysis was performed for each list of DEGs (responders against non-responders at time points M0 and M3m separately). RNA processing and negative regulation of fatty acid transport was revealed among the enriched biological processes using DEG list obtained by comparison of responders against non-responders before administration of metformin (Fig 2A).

PLOS ONE
Contrasting responders against non-responders after three months of metformin use identified enrichment of GO terms strongly related to the mitochondrial activity (e.g. mitochondrial respiratory chain complex I, ATP metabolic process) (Fig 2B). See S3 and S4 Tables for the full list of enriched GO terms identified. Although KEGG pathway analysis of DEG list, obtained from the baseline analysis (M0), did not reveal any significantly enriched cell signaling pathway, comparison of responders against non-responders after three months of metformin use showed differential expression of 5 mitochondrial genes (MT-ATP6, MT-ND2, MT-ND4, MT-ND4L, MT-ND6). All of these genes provided a significant enrichment of the following pathways: hsa00190: Oxidative phosphorylation (p-value = 4.47E-06), hsa05012: Parkinson's disease (p-value = 5.81E-06) and hsa01100: Metabolic pathways (p-value = 2.18E-02).

Discussion
The results of our study showed metformin-specific signatures in blood cell transcriptome profiles associated with some of its well-known properties: the ability to improve energy metabolism, influence immune responses, and inhibit cancer progression. This study has demonstrated a gene expression-based molecular discrimination between metformin responders Log2FC, log2 fold change; FDR, false discovery rate. † Genes showing significant differential expression comparing responders against non-responders also after metformin therapy for 3 months. https://doi.org/10.1371/journal.pone.0237400.t004 and non-responders, and suggested that mitochondrial respiratory complex I may be associated with metformin efficacy. To the best of our knowledge, this is the first longitudinal study focusing on metformin-induced transcriptional alterations of drug-naïve T2DM patients in vivo.
Our study revealed metformin-induced differential expression of genes involved in cholesterol homeostasis, such as solute carrier family 46 member 1 (SLC46A1), which is involved in the intestinal folate absorption affecting plasma high-density lipoprotein levels [39], and lipoprotein receptor-related protein 1 (LRP1), a crucial protein for cholesterol uptake. Metformininduced reduction of hepatic LRP1 expression level has been reported before [40], though here we report a similar effect in blood cells for the first time. Metformin-induced differential expression of both genes may serve as a contributing factor for the cholesterol-lowering effect of the drug. In our data, the same mechanism was supported by the enrichment of lipoprotein particle receptor activity among identified GO terms.
In addition, we found significant downregulation of multiple cancer-related genes coding for cytochrome P450 1B1 (CYP1B1), C-C chemokine receptor type 2 (CCR2), stabilin-1 (STAB1) and transmembrane protein 176B (TMEM176B). Some of the observed alterations have been previously identified in different tissue types, such as the downregulation of CCR2, which has been explained by the ability of metformin to block M2-like polarization of tumorassociated macrophages providing the anti-metastatic effect of the drug [41], or downregulation of CYP1B1 in breast cancer cells where due to its crucial role in estrogen metabolism, metformin has been suggested as a potential chemopreventive agent against carcinogenesis [42]. Multiple studies have reported that aberrant expression of STAB1 is related to the tumor progression in various types of cancer, serving as a potential molecular target for cancer therapy [43,44]. Similarly, elevated protein levels of TMEM176B are detected in multiple malignancies, moreover, inhibition of TMEM176B has already been proved to promote CD8 + T cell-mediated tumor growth control, enhancing the therapeutic efficacy of cancers [45,46]. Therefore, the observed downregulation of both, STAB1 and TMEM176B may highlight novel players in the anti-cancer activity of metformin.
Finally, the enrichment of antigen processing and presentation pathway and significantly reduced expression of genes coding for the cluster of differentiation 14 (CD14), which was already proved to be differentially expressed in metformin-treated monocyte cells, and the cluster of differentiation 163 (CD163) [47], a scavenger receptor which has been previously associated with insulin resistance in patients with T2DM, altogether may explain well-known participation of metformin in the inflammatory and immune responses [48].
According to our previous study, there are only three genes (HBB, PDZK1IP1, and RN7SL679P) showing metformin-induced differential expression in blood cells obtained from both, T2DM patients and healthy volunteers. Although both studies were longitudinal, the duration of the therapy (7 days for healthy volunteers and 3 months for T2DM patients) and the dose of metformin (850 mg twice a day for healthy volunteers and variable dose in the T2DM patient group), differed between both studies, which may explain the observed variability of metformin effects on the blood cell transcriptome profiles among both study groups. Moreover, the general health status including the presence of T2DM may serve as the major confounding factor to modulate the cell metabolism and effects of metformin in the patient group compared to healthy individuals [49].
One of the main findings of our study was the clear, transcriptome-based discrimination of study subjects into metformin responders and non-responders before the administration of any anti-diabetic therapy. Differential expression analysis before the use of metformin revealed significant upregulation of IRS2 gene coding for insulin receptor-2 (log 2 FC 0.89) in responders compared to non-responders. So far, multiple studies have demonstrated the ability of metformin to activate the hepatic insulin receptor and the IRS2/PI3K/Akt pathway resulting in reduced insulin resistance and increased glucose uptake [50,51]. The results of our study suggest that the activity of IRS2 in blood cells prior to the administration of metformin may be related to the efficacy of the therapy. Nevertheless, this hypothesis must be tested in a larger longitudinal cohort. PLS-DA analysis also provided strong evidence for transcriptome-based patient stratification which was mainly explained by a single latent variable (component) and defined by expression of 56 key genes (VIP>1). We found carcinoembryonic antigen-related cell adhesion molecule 1 (CEACAM1), Insulin receptor substrate 1 (IRS1), ABCC2 gene coding for multidrug resistance protein 2 and IGHA1 gene coding for the constant segment of immunoglobulin A heavy chain among the marker genes determining the distribution of patients in metformin response groups. So far, there are no studies reporting the role of multidrug resistance protein 2 in metformin efficacy, nevertheless, its homolog mitochondrial multidrug resistance protein 1 has been already associated with resistance to metformin in human malignant mesothelioma cells [52]. Interestingly, CEACAM1 is a mediator of insulin clearance in the liver [53] and IRS1 is playing a key role in glucose homeostasis [54], though no reports were linking their activity with metformin efficacy so far. In our previous studies, we have already reported the contribution of IgA, the most abundant intestinal antibody shaping the gut microbiome composition, in metformin action [27], and here we suggest a potential implication of IgA in microbiome-mediated metformin response [55]. Together these genes showed notable patient segregation in metformin response groups also when visualized in CPM-based heat map, therefore they may be considered as candidates for further studies of biomarkers for metformin response.
The comparison of blood cell transcriptome profiles from responders with transcripts from non-responders after the metformin therapy for three months revealed significant downregulation of mitochondrial genes (MT-ATP6, MT-ND2, MT-ND4, MT-ND4L, MT-ND6). Furthermore, MT-ND4 and MT-ND4L also showed a positive association with HbA1c levels according to the multiple regression model. Metformin directly acts on mitochondria and limits respiration by inhibiting mitochondrial respiratory-chain complex 1 (NADH: ubiquinone oxidoreductase), which is among the top targets of the drug, and it also catalyzes oxidative phosphorylation in mammalian mitochondria [56,57]. Thus, the observed downregulation of genes coding for NADH: ubiquinone oxidoreductase core subunits in responders and enrichment of oxidative phosphorylation seems rational. These data suggest altered mitochondrial complex I activity as one of the mechanisms linked to the variability of metformin response.
There are a few limitations in this study that could be addressed by future research. First, a small sample size due to the essential inclusion criterion of drug-naïve T2DM patients. Second, prolonged observation time (three months) that provides an advantage to evaluate long-term metformin effects, while causing an issue of possible accumulation of uncontrollable factors affecting gene expression and masking the true effects of the drug. Third, the lack of control arm allowing the elimination of these confounding factors. Nevertheless, we believe that RNA-Seq, a highly sensitive and accurate method for gene expression, together with the repeated measures within the longitudinal study design and strictly defined inclusion criteria ensures the absence of potential influence of other anti-diabetic therapies and improves the validity of our results.
In conclusion, the current study applying RNA-Seq for the discovery of transcriptional effects of metformin in drug-naïve T2DM patients provided detailed insight into potential molecular mechanisms underlying well-known beneficial effects of metformin. However, since there are no data confirming the accumulation of metformin in other blood cells than erythrocytes [58], it seems reasonable to assume that the observed effects of metformin on peripheral blood cells most likely are indirect and reflects the systemic consequence of the therapy. Nevertheless, we suggest that blood-derived transcriptome profiles may be used for evaluation of therapeutic efficacy and specific genes may be further applied in the development of biomarkers for metformin response.
Supporting information S1 Table.