Long-term exposure of human endothelial cells to metformin modulates miRNAs and isomiRs

Increasing evidence suggest that the glucose-lowering drug metformin exerts a valuable anti-senescence role. The ability of metformin to affect the biogenesis of selected microRNAs (miRNAs) was recently suggested. MicroRNA isoforms (isomiRs) are distinct variations of miRNA sequences, harboring addition or deletion of one or more nucleotides at the 5′ and/or 3′ ends of the canonical miRNA sequence. We performed a comprehensive analysis of miRNA and isomiR profile in human endothelial cells undergoing replicative senescence in presence of metformin. Metformin treatment was associated with the differential expression of 27 miRNAs (including miR-100-5p, -125b-5p, -654-3p, -217 and -216a-3p/5p). IsomiR analysis revealed that almost 40% of the total miRNA pool was composed by non-canonical sequences. Metformin significantly affects the relative abundance of 133 isomiRs, including the non-canonical forms of the aforementioned miRNAs. Pathway enrichment analysis suggested that pathways associated with proliferation and nutrient sensing are modulated by metformin-regulated miRNAs and that some of the regulated isomiRs (e.g. the 5′ miR-217 isomiR) are endowed with alternative seed sequences and share less than half of the predicted targets with the canonical form. Our results show that metformin reshapes the senescence-associated miRNA/isomiR patterns of endothelial cells, thus expanding our insight into the cell senescence molecular machinery.

Changes in the isomiR pattern associated with metformin treatment. Since studying miRNAs at the isomiR level could lend new insights into miRNA biology and function, we analysed isomiR modulation associated with metformin treatment of HUVECs during replicative senescence. IsomiRs result from a shift of the cutting site of Drosha/Dicer enzymatic activities during miRNA biogenesis 23,24 and can be classified into six categories according to the types of sequence modifications: (1) canonical miRNAs, (2) 3′ deletion isomiRs, (3) 3′ addition isomiRs, (4) 5′ deletion isomiRs, (5) 5′ addition isomiRs, and (6) mixed isomiRs, which represent a combination of the prior categories 25 . We also analysed the post-transcriptional addition of one or more uridines at the 3′ end of isomiRs and canonical miRNAs, namely uridylation. It has to be noted that the entire spectrum of isomiRs is covered by the standard miR-seq analysis. Figure 3a shows the contribution of different sequence isoforms to the total miRNA pool in SEN + M. On a total of 3,632,423 reads, the 43.1% was mapped to non-canonical isoforms. No statistically significant difference in the proportion of isomiR variations between SEN and SEN + M was observed (p = 0.103).
Of note, 39 of 133 metformin-modulated isomiRs showed a modification at the 5′ end (Fig. 3c), which leads to a shift of the seed sequence resulting in a change of the miRNA-target binding site 26 . Moreover, in some instances the seed sequence of the 5′ isomiR is identical to the seed sequence of another canonical microRNA. Indeed, miR-27b-3p|+ 3|0, miR-29a-3p|-1|-2, miR-34a-5p|+ 1|+ 1 and miR-423-5p|+ 2|0 share the same seed sequence of miR-5693, miR-5682, miR-6499-3p and miR-486-3p, respectively. Among the differentially regulated isomiRs, the most frequent modification was the 3′ deletion. Furthermore, 3′uridylation was extensively represented among all isomiR types, except those presenting a 5′ nucleotide addition (Fig. 3c). Interestingly, metformin affected the expression of 24 3′-uridylated miRNAs (Fig. 3b, www.nature.com/scientificreports/ Figure 4a shows the expression of the different isomiR variants of the 73 miRNAs including at least one isomiR differentially regulated by metformin. Notably, we observed a high variability in the proportion of isomiR variants among the evaluated miRNAs. Indeed, the canonical form, i.e. the one reported in the miRBase database, is not always prevalent (e.g. in the miR-30 family) and some miRNAs included non-canonical variants.
Moreover, metformin induced a significant redistribution of the isomiR variant proportions within 6 out of the 27 differentially regulated miRNAs (Fig. 4b). Only one isoform was detected for 3 miRNAs (miR-17-3p, miR-98-3p, and miR-92-1-5p), whereas no isoforms were detected for the remaining 12 miRNAs. The proportions of the different isomiR variation types between SEN and SEN + M are reported in Table 1.

MiRNA/isomiR expression trends in endothelial cells during replicative senescence. To gain
insight into the biological significance of miRNA/isomiR modulation induced by metformin treatment in senescent HUVECs (SEN, SEN + M) we used non-senescent HUVECs as control (Young, SA β-gal < 5%, Fig. 1b) 22 . This strategy allowed us to identify two different trends in miRNA modulation and to separate the 27 miRNAs according to their pattern of modulation. On one hand, 13 miRNAs were characterized by linear increasing or decreasing trend, when the Young/SEN/SEN + M sequence was examined (Fig. 5a, group 1). The most evident linear trends were observed for miR-100, -125b-5p, and -654-3p, showing the most abundant expression. Metformin treatment further increases the expression of these miRNAs in senescent cells. On the other hand, 14 miRNAs were characterized by a 'U-shaped'/ 'inverted U-shaped' trend of the Young/SEN/SEN + M sequence, suggesting that metformin induced a (partial) reversal of the miRNA expression induced by senescence (Fig. 5b, group 2). Among the 14 miRNAs belonging to this latter group, the 'U-shaped'/ 'inverted U-shaped' trend was confirmed by a significant likelihood ratio test (adjusted p value < 0.05) (in bold in Fig. 5b). The most relevant inverted U-shaped trends were observed for miR-217-5p, -216a-3p, and -216a-5p. Overall, metformin affected the expression of 18 SA miRNAs (Fig. 5c) and notably was able to rescue the miR-216a and miR-217-5p overexpression in senescent cells previously reported by our group 22 .

Metformin alters the miRNA and isomiR targetome of senescent endothelial cells.
To explore target genes and pathways affected by the 11 miRNAs showing significant U-shaped or inverted U-shaped trends, pathway enrichment analysis was performed using the miRPath v.3/Diana tool. Figure 6a lists the involved KEGG pathways (p < 0.01) ranked by the significance of the enrichment. The proportion of targeted genes over total genes for each pathway is also reported. A considerable number of pathways related to cell proliferation, i.e. TGFβ, ErbB, Wnt and MAPK pathways, is significantly enriched. Notably, the PI3K-Akt-mTOR pathway is the one containing the greatest amount of targeted genes (72.3%), in agreement with the inhibitory effects of metformin on mTOR signaling 27 .  www.nature.com/scientificreports/ Among these miRNAs, we focused on the only two miRNAs showing detectable levels of at least one 5′ isoform, i.e. miR-217-5p and miR-216a-3p. Interestingly, the canonical/3′ and the 5′ miR-217-5p isomiRs share only half of the predicted targets, while the other half is exclusive to either seed sequence. Regarding miR-216a-3p, the deletion of one or two nucleotides at the 5′ end leads to the generation of two alternative seed sequences. The 3 different seed sequences shared only a small pool (35) of predicted targets (Fig. 6b). The target genes of canonical and 5′isomiR seed sequences were evaluated also for those miRNAs including at least one 5′isomiR presenting  Figure S2). Notably, the canonical form and the 5′ isomiR of miR-100-5p shared no predicted target genes.

Discussion
In the present study, we investigated for the first time the miRNA landscape in endothelial cells (ECs) undergoing replicative senescence after a long-term treatment with metformin. Surprisingly, only 27 miRNAs on a total of 1706 detected by the small RNA-seq analysis were differentially regulated by metformin, despite the long duration of the exposure to a pharmacologically pertinent dose of the drug. To gain insight into the biological significance of these modulations, we used young proliferating HUVECs as reference group, in order to identify specific trends of modulation. We focused on the group of miRNAs characterized by a U-shaped/inverted U-shaped trend of expression in young vs. SEN vs. SEN + M, since this peculiar trend could reflect the ability of metformin to modulate the trajectories of senescence associated miRNAs (Fig. 7a). Increasing evidence suggests that a number of biomarkers of human aging followed non-linear trends when subjects representing the extreme phenotype of successful aging, i.e. the centenarians, are included in the analysis [28][29][30][31][32] . We therefore employed an in vitro cellular senescence model mimicking the gradual deterioration of endothelial function that accompanies human aging 33 , to unravel the ability of metformin to affect the senescence-associated miRNA/isomiR modulation. www.nature.com/scientificreports/ This approach allowed us to show that metformin can revert the SA trend of a number of miRNAs that were extensively studied in the context of cellular aging, including miR-216-3p, -216-5p, and -217-5p, which we previously identified among the most upregulated miRNAs in senescent HUVECs 22 . MiR-217-5p was proved to be involved in EC and human fibroblast senescence by targeting SIRT1 and DNMT1, respectively 34,35 . Furthermore, we recently demonstrated that the same pro-senescence effects of miR-217 can be spread through the exchange of small extracellular vesicles 22 . Similarly, miR-216a was shown to be involved in EC aging, in atherosclerosisrelated endothelial dysfunction by impairing the autophagy response to the accumulation of oxidized low-density lipoproteins 36 , and in macrophage pro-inflammatory M1 polarization by boosting the NF-κB pathway 37,38 . Among miRNAs showing a linear trend in Young vs. SEN vs. SEN + M, miR-100-5p was previously shown to be upregulated in senescent HUVECs 22,39 , while the metformin-mediated induction of miR-125-5p was consistent with previous reports on macrophages 40 and senescent ECs 12 .
Regarding the analysis of isomiRs, this is the first deep sequencing assessment of isomiRs in senescent HUVECs. One miRNA gene can potentially produce multiple distinct isomiRs, differing in length, sequence, or both 26 . Our results proved that the assessment of isomiRs can unravel complex modulations of the miRNA pool not detectable with standard miRNA analysis. Indeed, isomiR analysis allowed us to fully uncover the downregulating effects of metformin on the miR-17/92 cluster, which has been previously shown to be over-represented in a wide range of cancers and cardiovascular diseases and downregulated in physiological aging 41 . Therefore, further developments of isomiR analysis are warranted to increase our knowledge on miRNA modulation in a number of physiological and pathological processes. www.nature.com/scientificreports/ In agreement with previous reports 25, 42 , we observed a considerable presence of 3′ isomiRs, while more than half of the total reads was mapped to canonical miRNAs. It has to be noted, however, that the term 'canonical' refers to the sequence annotated in miRBase and do not necessarily indicate the most abundant miRNA isoform in a specific cell type or tissue or the primary product of pre-miRNA cleavage 43 .
On the other hand, only a small number of reads (about 3%) mapped to 5′ isomiRs, which are associated to a shifting of the seed sequence (Fig. 7b). For this reason, we evaluated the number of targets shared by the isoforms of miR-216a-3p and miR-217-5p, which were both modulated by metformin treatment and expressed 5′ isoforms. The inclusion of these additional seed sequences into the targetome analysis yielded a considerably greater number of target genes, most of which were not shared with the canonical miRNAs. As expected, the coexistence of more than one 5′ isomiR, as in the case of miR-216-3p, proportionally increased the amount of target genes. The ability of isomiRs of being loaded onto the RISC complex support their possible biological role 19,44 . Indeed, a previous report showed that the ratio between miR-411 and its 5′ isomiR in ECs is affected by acute ischemia and that only the 5′ isoform of miR-411 is capable of impairing angiogenesis by targeting a different subset of mRNAs 19 .
The 3′ end miRNA modifications are mostly related to post transcriptional deletion of nucleotides, i.e. trimming, or the addition of one or more nucleotides, i.e. tailing 45 . It has to be noted, however, that is quite challenging to distinguish templated nucleotides added during miRNA maturation from those added post-transcriptionally to the mature miRNA. In our study, we assessed isomiRs resulting from the untemplated nucleotide addition to the 3′ end of pre-miRNA or mature miRNA 46 . While these modifications are not associated with a shifting of the seed sequence, it has been demonstrated that 3′ uridylation enhances base-pairing between tailed miRNA and targets, a phenomenon named as tail-U-mediated repression (TUMR). Therefore, TUMR expands the miRNA target repertoires by producing novel miRNA-target binding sites in the presence of an incomplete seed-pairing 18 . Moreover, 3′ post-transcriptional modifications were shown to affect miRNA stability 47 , intracellular levels, and compartmentalization into extracellular vesicles 20 . Notably, miRNAs are not the sole substrates of the 3′ uridylation mediated by terminal uridyltransferases (TUTs). Indeed, 3ʹ-terminal uridylation of viral RNAs in mammalian cells has been recently identified as a conserved antiviral defense mechanism 48 . Interestingly, metformin affected the expression of 22 3′-uridylated miRNAs; therefore, it is straightforward to conceive a framework in which metformin could impact cellular senescence through the modulation of miRNA function, stability, and localization.
Our in vitro results support the role of isomiR assessment in biological samples as a useful tool to improve our knowledge on the aging process or discover new biomarkers of biological aging.
Nevertheless, several limitations need to be acknowledged. The study design does not allow to draw any mechanistical conclusion on the role of metformin on ECs or cellular senescence. In addition, some of the mechanisms of isomiR biogenesis are still unclear, implying the intrinsic difficulty to assess whether 3′ nucleotide addition occurs during or after miRNA transcription. Finally, qPCR validation of NGS assessment of isomiRs is still hampered by analytical challenges, such as the absence of dedicated protocols and reagents, e.g. probes www.nature.com/scientificreports/ and primers, and the relative inefficiency of the currently available techniques in differentiating highly similar sequences 49,50 . NGS studies on isomiRs paved the way to the exploration of novel non-canonical targets and allowed the identification of new regulatory mechanisms of miRNA expression and intracellular localization, adding an additional layer of complexity to the study of the epigenetic variations accompanying cell senescence, although further investigations are required to better understand the biological functions of the cellular isomiR pool.
Overall, we showed that long-term treatment with metformin is able to partly attenuate the complex miRNA/ isomiR remodeling observed during cellular senescence in ECs, supporting further exploration of the impact of metformin on the cellular epigenetic landscape as a possible mediator of the putative beneficial effect of this drug on the aging process.

Materials and methods
Cell culture and treatment. An in vitro model of endothelial replicative cell senescence was established using long-term cultured HUVECs. Cryopreserved HUVECs obtained from pool of donors were purchased from Clonetics (Lonza, Switzerland) and cultured in EGM-2 (CC-3162, Lonza) at 37 °C in a humidified atmosphere containing 5% CO 2 . Cells were seeded at a density of 5000/cm 2 and sub-cultured when they reached 70-80% confluence. All cells tested negative for mycoplasma infection. Before replating, harvested cells were counted using a hemocytometer. Population doublings (PDs) were calculated by the formula: (log 10 F -log 10 I)/ log 10 2, where F is the number of cells at the end of the passage and I is the number of seeded cells. Cumulative population doubling (cPD) was calculated as the sum of PD changes. Cells were cultured until the arrest of replication and classified based on SA β-galactosidase (β-gal) activity into young (SA β-gal < 5%) and senescent (SEN, SA β-gal > 80%) cells using Senescence Detection Kit (cat. no. K320, BioVision Inc., USA) as described previously 21 . Cells were treated with 20 μM metformin (cat. D150959, Sigma Aldrich, Italy) added at each medium replacement. mRNA expression level. CDKN2A mRNA expression was assessed as previously described 22  Small RNA sequencing analysis. Small RNA sequencing was performed in triplicate on Young and SEN cells, and SEN cells treated with metformin. TruSeq Small RNA Library PrepKit v2 (Illumina; RS-200-0012/24/36/48) was used for library preparation according to the manufacturer's indications. Briefly, 35 ng purified RNA was linked to RNA 3′ and 5′ adapters, converted to cDNA, and amplified using Illumina primers containing unique indexes for each sample. Each library was quantified using Agilent Bioanalyzer and High Sensitivity DNA Kit (cat. no. 5067-4626, Agilent Technologies, USA) and equal amounts of libraries were pooled together. Size selection allowed keeping 130-160 bp fragments. After ethanol precipitation, the library pool was quantified with Agilent High Sensitivity DNA Kit, diluted to 1.8 pM, and sequenced using NextSeq 500/550 High Output Kit v2 (75 cycles) (Illumina; FC-404-2005) on the Illumina NextSeq500 platform.

RNA extraction.
Raw base-call data generated by the Illumina NextSeq 500 system were demultiplexed using Illumina Bas-eSpace Sequence Hub (https ://bases pace.illum ina.com/home/index ) and converted to FASTQ format. After a quality check with FastQC (https ://www.bioin forma tics.babra ham.ac.uk/proje cts/fastq c/), sequence reads were quality trimmed using the cutadapt tool 51 . Sequence reads were aligned to the miRBase version 21.0 database 52 using the STAR algorithm 53 . Standard miRNA quantification (including the canonical form and all isoforms) was obtained as previously detailed 22 . Quantification of miRNA isoforms. Sequence reads were quality trimmed using the cutadapt tool, and mapped unambiguously using SHRIMP2 54 to the human genome assembly GRCh38. During the mapping, no insertions or deletions, and at most one mismatch was permitted. IsomiRs were identified as done previously 16,17,[55][56][57] . The isomiR nomenclature used is based upon the one used previously in Loher et al. 17 . For example, the isomiR whose 5′ terminus begins one position to the right (+ 1) of the archetype's 5′ terminus and whose 3′ terminus ends two positions to the left (− 2) of the archetype's 3′ terminus is labeled " + 1|− 2". The archetype isomiR, the sequence found in public databases, is labeled as "0|0".
IsomiR abundances were quantified in reads per million (RPM). Only reads that passed quality trimming and filtering and could be aligned exactly to miRNA arms were used in the denominator of this calculation. The abundance of a miRNA arm is calculated as the sum of normalized abundances of all isomiRs from the arm.