Profiling of Differentially Expressed MicroRNAs in Human Umbilical Vein Endothelial Cells Exposed to Hyperglycemia via RNA Sequencing

Hyperglycemia is the hallmark of diabetes mellitus that results in oxidative stress, apoptosis, and diabetic vascular endothelial dysfunction. An increasing number of microRNAs (miRNAs) have been found to be involved in the pathogenesis of diabetic vascular complications. However, there is a limited number of studies that characterize the miRNA profile of endothelial cells exposed to hyperglycemia. Therefore, this study aims to analyze the miRNA profile of human umbilical-vein endothelial cells (HUVECs) exposed to hyperglycemia. HUVECs were divided into two groups: the control (treated with 5.5 mM glucose) and hyperglycemia (treated with 33.3 mM glucose) groups. RNA sequencing identified 17 differentially expressed miRNAs between the groups (p < 0.05). Of these, 4 miRNAs were upregulated, and 13 miRNAs were downregulated. Two of the most differentially expressed miRNAs (novel miR-1133 and miR-1225) were successfully validated with stem-loop qPCR. Collectively, the findings show that there is a differential expression pattern of miRNAs in HUVEC following exposure to hyperglycemia. These 17 differentially expressed miRNAs are involved in regulating cellular functions and pathways related to oxidative stress and apoptosis that may contribute to diabetic vascular endothelial dysfunction. The findings provide new clues on the role of miRNAs in the development of diabetic vascular endothelial dysfunction, which could be useful in future targeted therapy.


Introduction
Cardiovascular diseases are one of the major complications of diabetes leading to mortality and morbidity in diabetic patients [1]. Types 1 and 2 diabetes mellitus, despite their different pathogenesis, are characterized by hyperglycemia [2]. Hyperglycemia is responsible for the development of diabetic vascular endothelial dysfunction that eventually leads to the macro-and microvascular complications of diabetes. Endothelial cells are vital for the regulation of vascular homeostasis and are sensitive to changes at the blood glucose level. A growing body of evidence demonstrates that hyperglycemia causes endothelial dysfunction, manifested as increased oxidative stress and apoptosis [2,3]. Endothelial cells exposed to hyperglycemia have increased hydrogen peroxide formation, which promotes oxidative damage to the vasculature [4]. A recent study showed Life 2023, 13, 1296 3 of 16 consent when they participated in the study. HUVECs were isolated from the umbilical veins by using an enzymatic technique as described previously [17]. Briefly, umbilical cords were cleaned and infused with 0.01% (w/v) collagenase (Worthington Biochemical Corporation, Lakewood, NJ, USA). Isolated HUVECs were cultured in an endothelial-cell medium supplemented with 5% fetal bovine serum, 1% endothelial cell growth, and 1% penicillin and streptomycin (ScienCell Research Laboratories, Carlsbad, CA, USA) at 37 • C in a humidified 5% carbon dioxide incubator. The medium was changed every other day until the cells had reached confluence. HUVECs were identified through the typical cobblestone morphology and positive endothelial cell markers (von Willebrand factor and CD31) in immunocytochemistry.
All the experiments were performed using passages 3-4 of HUVECs at 80% confluency. HUVECs were divided into two groups: the control (treated with 5.5 mmol/L D-glucose) and hyperglycemia (treated with 33.3 mmol/L D-glucose) groups. All treatments were given for 24 h. The dose of glucose and the treatment duration were based on the model of previous studies [18][19][20][21].

Total RNA Extraction
Following 24 h of treatment, total RNA was extracted using the miRNeasy Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. The extracted total RNA was then assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), for its integrity, purity, and concentration to guarantee high-quality total RNA for small RNA sequencing.

miRNA Expression Profiling
The RNA samples were further processed at Genewiz, New Jersey, where the complimentary DNA (cDNA) libraries were prepared and sequenced using Illumina Hiseq, 2 × 150 bp paired-end chemistry. Raw reads were preprocessed for adaptor trimming, quality, and size selection using Trimmomatic (V0. 30). Trimmed data were also subjected to quality assessment and assessed using FastQC (V0. 10.1). Sequences that passed the quality filtering were mapped to the reference genome, Homo sapiens (GRCh38.104), using Bowtie2 (V2.1.0). The clean reads were aligned to the miRbase database to identify the known miRNAs, and miRDeep2 software was used to predict the novel miRNAs. Raw reads were then normalized to transcripts per million (TPM). Differentially expressed miRNAs were identified using a negative binomial distribution model, available through DESeq2. The generated data in this study were submitted to Gene Expression Omnibus (GEO), accession number GSE229207.

Validation of the Selected miRNAs with Stem-Loop Reverse-Transcription-Quantitative Polymerase Chain Reaction
To validate the RNA sequencing data, two of the most differentially expressed miRNAs (novel miR-1133 and miR-1225) were validated via stem-loop quantitative polymerase chain reaction (qPCR). A total of 5 ng of small RNA-enriched total RNA was synthesized by using an MMLV Reverse Transcriptase 1st-Strand cDNA Synthesis Kit (Biosearch Technologies, Hoddesdon, UK) with an additional 0.1 µM of the stem-loop primer. The qPCR of the novel miR-1133 and miR-1125 was performed in 10 µL of total reaction consisting of 5 µL of 2× miRCURY SYBR Green Master Mix, 0.5 µM of each forward and reverse primer, and 3 µL of the synthesized cDNA (60× diluted). All the primer sequences are listed in Table 1. All reactions were prepared in a 96-well plate format, and qPCR was performed using a CFX96™ Real-Time PCR detection system (Bio-Rad, Hercules, CA, USA). The reactions were conducted using the following thermocycling conditions: initial denaturation at 95 • C for 10 min, followed by 40 cycles of amplification at 95 • C for 10 s, annealing at 60 • C for 30 s, and elongation at 72 • C for 10 s, with an additional elongation step at 40 • C for 1 s. The reaction kinetic of each primer set and protocol was verified with the melting profile. The amplification signals were acquired during the elongation step and recorded Life 2023, 13, 1296 4 of 16 by using Bio-Rad CFX Manager software version 3.1. The threshold cycle (C T ) value was used to calculate the expression of the miRNAs. The U6 gene was used to normalize the quantitative analysis. The relative miRNA expression was calculated with the 2 −∆∆CT method, with ∆∆C T = C T of miRNA − C T U6.

Bioinformatics Analysis
GOSeq (v1.34.1) was used to identify GO terms that annotated a list of enriched genes with a significance value of p < 0.05. The significant differential gene expression was enriched in KEGG pathways. Following adjustment with Benjamini and Hochberg's approach for controlling the false discovery rate, the p-value of miRNAs was set to <0.05, and the fold change of genes was set to a minimum of a twofold change to detect differentially expressed miRNAs.

Prediction and Analysis of the Potential Target Genes of Differentially Expressed miRNAs
The potential target genes of differentially expressed miRNAs were identified using in silico prediction analysis through independent prediction miRDB algorithms (http: //mirdb.org, accessed on 1 March 2023). To date, the only miRNA-mRNA prediction tool available for novel miRNAs is the miRDB algorithm. Since novel miRNAs were part of the differentially expressed miRNAs in this study, the miRDB algorithm was used to predict the target genes and subsequently validate the downstream target of miRNAs.
The predicted target genes in miRDB v6.0 were downloaded using the default parameters with prediction scores ranging from 50 to 100. The predicted target genes of the differentially expressed miRNAs were mapped into the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (v6.8) (https://david.ncifcrf.gov/, accessed on 26 March 2023). The functional annotation clustering of differentially expressed miRNAs was performed against a GO term with medium stringency settings (kappa similarity term overlap of three or more, similarity threshold of 0.50, minimal three group members ab initio with EASE scores of ≥1) to characterize the overall role and function of differentially expressed miRNAs target genes. Fisher's exact and Chi-square tests were used to determine the significance of the GO terms. Only the GO terms with p-value of <0.05 were selected. Then, the protein-protein interaction (PPI) network of target genes was constructed using the STRING database (http://string-db.org, accessed on 28 March 2023). Subsequently, the hub genes in the network were identified by analyzing the degree of connectivity using Cytoscape software (version 3.6.1).

Statistical Analysis
GraphPad Prism version 9 software was utilized to analyze the qPCR data. The Shapiro-Wilk test was used to determine the normality of the data. Results are expressed as mean ± SEM. Differences between the groups were evaluated using unpaired t-test. p-values < 0.05 were considered statistically significant.

RNA Sequencing Identified Differentially Expressed miRNAs in Hyperglycemia-Induced HUVECs
The deep RNA sequencing of HUVECs following exposure to normal and high glucose for 24 h identified a total of 814 known miRNAs, of which 580 miRNAs were detected in both the control and hyperglycemia groups. We also found a total of 1229 novel miRNAs in these samples ( Figure 1). Known miRNAs refer to miRNAs that have been identified and documented in miRNA databases. On the other hand, unknown or novel miRNAs are miRNAs that have not been identified or documented in previous studies or established miRNA databases. However, MiRDeep2 prediction is able to determine possible novel miRNAs on the basis of the secondary structure of miRNA precursors [22].
NAs are miRNAs that have not been identified or documented in previous studies or established miRNA databases. However, MiRDeep2 prediction is able to determine possible novel miRNAs on the basis of the secondary structure of miRNA precursors [22].

Validation of Novel miR-1133 and miR-1225 via Stem-Loop qPCR
To validate our RNA sequencing data, we selected one of the most significantly upregulated miRNAs (novel miR-1133) and one of the most significantly downregulated miRNAs (novel miR-1125). These miRNAs were further quantified with stem-loop qPCR. The results confirmed that novel miR-1133 was significantly upregulated by 8.7-fold in HUVECs exposed to hyperglycemia compared to the control HUVECs (p < 0.001). Meanwhile, novel miR-1225 was significantly downregulated by 1.62-fold in HUVECs following exposure to hyperglycemia (p < 0.01) ( Figure 3).

GO Analysis
To better understand the potential implications of the 17 differentially expressed miRNAs, a total of 5655 target genes of the 17 known and unknown, significantly expressed miRNAs were annotated using GO analysis, and their functions were categorized. The GO functional classification of all differentially expressed miRNAs is presented in the histogram and shown as the -log p-value with the specification of the relevant biological process, cellular component and molecular function (Figure 4). Among all, 6, 18, and 6 GO terms were significantly involved in molecular function, biological process, and cellular component, respectively, on the basis of the p-value.
The   regulated miRNAs (novel miR-1133) and one of the most significantly downregulated miRNAs (novel miR-1125). These miRNAs were further quantified with stem-loop qPCR. The results confirmed that novel miR-1133 was significantly upregulated by 8.7-fold in HUVECs exposed to hyperglycemia compared to the control HUVECs (p < 0.001). Meanwhile, novel miR-1225 was significantly downregulated by 1.62-fold in HUVECs following exposure to hyperglycemia (p < 0.01) (Figure 3).

GO Analysis
To better understand the potential implications of the 17 differentially expressed miRNAs, a total of 5655 target genes of the 17 known and unknown, significantly expressed miRNAs were annotated using GO analysis, and their functions were categorized. The GO functional classification of all differentially expressed miRNAs is presented in the histogram and shown as the -log p-value with the specification of the relevant biological process, cellular component and molecular function (Figure 4). Among all, 6, 18, and 6 GO terms were significantly involved in molecular function, biological process, and cellular component, respectively, on the basis of the p-value.

KEGG Pathway Analysis
The KEGG pathway database was used to identify diabetes-, oxidative stress-, and apoptosis-related genes involved in the pathways using differentially expressed genes. The the most significant pathways were Type 2 diabetes mellitus, retrograde endocannabinoid signaling, the GnRH signaling pathway, the Wnt signaling pathway, the MAPK signaling pathway, dopaminergic synapse, hippo signaling pathway-fly, insulin secretion, RIG-I-like receptor signaling pathway, prolactin signaling pathway, pancreatic cancer, adipocytokine signaling pathway, hepatitis B, neurotrophin signaling pathway, epithelial cell signaling in Helicobacter pylori infection, Ras signaling pathway, ErbB signaling pathway, flavone and flavonol biosynthesis, glycosaminoglycan biosynthesis-keratan sulfate, and calcium signaling pathway with p-values between 6.55 × 10 −143 and 6.10 × 10 −21 ( Figure 5).
Life 2023, 13, x FOR PEER REVIEW 10 of 17 Figure 5. Scatter plot of differentially expressed gene's KEGG enrichment. The differentially expressed genes were mainly enriched in flavone and flavonol biosynthesis, and Type 2 diabetes mellitus KEGG pathways. The size of the dot is positively correlated with the number of differential microRNA target genes in the pathway. Color code is to indicate different q-value ranges. The larger the rich factor is, the greater the degree of enrichment. The smaller the q-value is, the more significant the enrichment.

Target Prediction and Analysis of Candidate for Differentially Expressed miRNAs
As shown in Table 3, there were a total of 1384 and 6572 predicted targets of the upregulated and downregulated miRNAs, respectively. For the upregulated miRNAs, hsa-miR-10526-3p had the most target genes with the number of 591 target genes. For the downregulated miRNAs, hsa-miR-4429 possessed the most targets of 1047 genes (see Supplementary Materials Table S1 and Figure S1 for details). . Scatter plot of differentially expressed gene's KEGG enrichment. The differentially expressed genes were mainly enriched in flavone and flavonol biosynthesis, and Type 2 diabetes mellitus KEGG pathways. The size of the dot is positively correlated with the number of differential microRNA target genes in the pathway. Color code is to indicate different q-value ranges. The larger the rich factor is, the greater the degree of enrichment. The smaller the q-value is, the more significant the enrichment.

Discussion
Hyperglycemia is the hallmark of diabetes mellitus that results in oxidative stress, with subsequent damage to cellular components and apoptosis [24]. The role of miRNAs has been reported in the pathogenesis of various diseases including diabetes mellitus. In this study, we conducted a comprehensive profiling of 17 differentially expressed miRNAs in HUVECs exposed to hyperglycemic state, which is a hallmark of diabetes. Our subsequent integrated bioinformatic analyses identified potential pathways associated with hyperglycemia-induced endothelial oxidative stress.
RNA sequencing revealed that a number of miRNAs were present in both hyperglycemiainduced and control HUVECs. A total of 17 differentially expressed miRNAs were found comprising 13 downregulated and 4 upregulated miRNAs. We subsequently verified the expression of the most upregulated and downregulated miRNAs in the samples: novel miR-1133 and miR-1125 by stem-loop qPCR. The stem-loop qPCR results of novel miR-1133 and miR-1125 expression were in agreement with the results of high-throughput sequencing, which indicates that the high-throughput RNA-sequencing for screening the differentially expressed miRNAs in our study was accurate and reliable.
We then assessed the GO and pathway enrichment analysis that showed the functions and pathways targeted by the 17 differentially expressed miRNAs. In total, we discovered 3992 notable GO terms and 308 pathways. The wide target range suggests that these miRNAs may play pivotal roles in the biological processes of HUVECs in hyperglycemic state. The functions of the 17 differentially expressed miRNAs involve the regulation of apoptotic signaling pathway, response to water deprivation, cellular hyperosmotic response, histone lysine methylation, cellular response to organic substance, neuropeptide catabolic process, immune response-regulating signaling pathway, protein phosphorylation, and regulation of insulin secretion. Most of these functions are closely related to oxidative stress and apoptosis, and could significantly contribute to the pathogenesis of diabetic vascular complications. In addition, the miRNAs targeted multiple signaling pathways. One of the most enriched KEGG pathways of the 17 differentially expressed miRNAs was the Type 2 diabetes mellitus pathway, which is related to the study.
miRNA-mRNA network analysis connects the 17 differentially expressed miRNAs to their predicted target genes and the related biological functions, allowing for a deeper understanding of the critical role of dysregulated miRNAs. In the networks, certain miR-NAs function as network hubs, such as hsa-miR-10526-3p, mir-90 and miR-70 (upregulated miRNAs); and hsa-miR-4429, miR-28 and hsa-miR-4709-3p (downregulated miRNAs). Of these miRNAs, hsa-miR-10526-3p and hsa-miR-4429 were the most active miRNAs in the networks, targeting 591 and 1047 genes, respectively. Only one study has investigated hsa-miR-10526-3p. In silico analysis has shown that hsa-miR-10526-3p is 1 of the 14 human miRNAs with the highest number of nucleotide matches with the nucleocapsid region of SARS-CoV-2. Thus, inhibiting these miRNAs is a potential way to treat SARS-CoV-2 infection [25]. However, hsa-miR-10526-3p had not been reported to be involved in oxidative stress and apoptosis.
Nevertheless, hsa-miR-4429 has attracted much attention in recent years. A study showed that cadmium exposure downregulates hsa-miR-4429 in adenocarcinomic human alveolar basal epithelial cells [26]. Cadmium is a potent cell poison that induces oxidative stress through the generation of ROS and by depleting the intrinsic antioxidant reserves [26]. Another study revealed that advanced glycation end products downregulate hsa-miR-4429 to promote apoptosis and fibrosis in pelvic organ prolapse [27]. Moreover, circulating hsa-miR-4429 levels were reduced following ischemic stroke [27], and hsa-miR-4429 was downregulated in the retinal pigment epithelial cells of diabetic mice [28]. Although the specific mechanisms of function of hsa-miR-4429 are not well-understood, findings from previous studies demonstrated that the downregulation of has-miR-4429 exerts a modulatory effect on oxidative stress and apoptosis. Therefore, further studies are needed to elucidate the precise role hashsa-miR-10526-3hasnd hsa-miR-4429 in hyperglycemiainduced endothelial oxidative stress.
In the PPI network of predicted target genes, we found that EGFR had the highest overlap degree of interaction in both upregulated and downregulated miRNAs. EGFR is critically involved in regulating a variety of cellular functions, such as cell growth, proliferation, motility, and survival [29]. The activation of EGFR signaling has been implicated in numerous pathophysiological processes, including diabetes [30]. Emerging studies reported that EGFR can play a dual role. In physiological conditions, EGFR acts as a beneficial factor for the normal development and function of the cardiovascular system. However, when its regulation is disrupted, EGFR acts as a detrimental factor that contributes to cardiovascular pathologies. EGFR inhibition reduced systemic oxidative stress in diabetic mice [30]. Hyperglycemia also increases the ROS that disturbs the EGFR pathway, resulting in delayed corneal epithelial wound healing [31]. While it may not be surprising to see ubiquitous hub genes such as EGFR among the highly ranked genes in both upregulated and downregulated gene lists, it is important to consider their potential functional implications. The fact that EGFR is highly ranked in both lists suggests its involvement in multiple regulatory mechanisms or pathways affected by hyperglycemia. Hence, further study is required to validate the function of significant miRNAs on their target genes such as EGFR.
PI3K acts as an effector of EGFR in diabetes-induced vascular dysfunction [30]. We identified the PI3KCA gene among the highest degrees of PPI network interaction in both upregulated and downregulated miRNAs. PI3K/AKT signaling is a common target downstream of epidermal growth factor receptor/epidermal growth factor receptor family (EGFR/ErbB) signaling. This pathway can activate endothelial nitric oxide synthase to produce nitric oxide (NO), which plays an important role in regulating endothelial function. Impaired PI3K/AKT signaling pathway leads to endothelial dysfunction in diabetic mice as evidenced by the reduction in vasorelaxation and NO production [32]. The PI3K-Akt signaling pathway is also involved in enhancing insulin sensitivity [33]. Moreover, PI3K expression was decreased in oxidative stress-induced HUVECs [34]. Thus, further research is required to validate these target genes that are involved in hyperglycemia-induced oxidative stress in endothelial cells.
Although we performed a comprehensive analysis on the miRNA profile of HUVECs exposed to hyperglycemia, there are a few limitations in this study. First, the sample of the RNA sequencing was relatively small with three samples, hindering identifying less predominantly differentiated miRNAs due to the low statistical power. Second, despite the high similarity of miRNA expression in endothelial cells isolated from various locations, some diversity of the miRNA profile could be present among the different types of endothelial cells [35]. Nevertheless, the scope of this study was limited to the miRNA profile of hyperglycemia-induced oxidative stress in HUVECs. In order to identify the groups of miRNAs that cause vascular complications in a particular target organ, it is necessary to conduct further research on the endothelial cells derived from the blood vessels of the target organ of interest. Lastly, this study was limited by the absence of biological validation to confirm the functional impact of the significant miRNAs, including their effects on the target genes such as EGFR, hence needing further functional study in the future.

Conclusions
Collectively, this study showed that 17 miRNAs are differentially expressed in HU-VECs following exposure to hyperglycemia. Bioinformatic analyses predict the functions, signaling pathways, target genes, miRNA-gene, and protein networks related to oxidative stress and apoptosis, which may contribute to diabetic vascular endothelial dysfunction. The findings provide new clues on the role of miRNA in the development of diabetic vascular endothelial dysfunction, which could be useful in future targeted therapy. Further functional study is needed to validate the role of these miRNAs in regulating oxidative stress and apoptosis in endothelial cells exposed to hyperglycemia.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/life13061296/s1. Supplementary Table S1: list of target genes for differentially expressed miRNAs. Supplementary Figure S1: predicted target genes of (a) upregulated and (b) downregulated miRNAs, highlighting the presence of a higher number of connected genes. Supplementary  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
Publicly available datasets were analyzed in this study. These data can be found in GEO, accession number GSE229207.