Abstract Purpose A long-term “memory” of hyperglycemic stress, even when glycemia is normalized, has been previously reported in endothelial cells. However, the molecular mechanism of “metabolic memory” (MM) remains unknown. In this report, we sought to screen at the whole transcriptome level the genes that participate in MM. Methods In the present research, RNA sequencing was used to determine the protein-coding mRNA expression profiles of human umbilical vein endothelial cells (HUVECs) under normal-glucose concentration (LG), high-glucose concentration (HG), and MM. A series of bioinformatic analyses was performed. HG-induced MM-involved up-regulated genes (up-HGMMGs) and HG-induced MM-involved down-regulated genes (down-HGMMGs) were identified. Afterward, based on up-HGMMGs and down-HGMMGs, the biological functions and signaling pathways were analyzed using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). In addition, several of the identified genes were validated by RT-qPCR. Results A total of 726 HGMMGs were identified, including 210 down- and 516 up-HGMMGs, which were enriched in the cell cycle (hsa04110), oocyte meiosis (hsa04114), p53 signaling pathway (hsa04115), and oxidative phosphorylation (hsa00190), among others. The protein–protein-interaction (PPI) network consisted of 462 nodes and 2656 connections, and four main modules were identified by MCODE. The cell cycle (hsa04110), oocyte meiosis (hsa04114), p53 signaling pathway (hsa04115), and oxidative phosphorylation (hsa00190), among others, could be potential therapeutic targets of HG-induced MM in endothelial cells. The real-time PCR results validated the RNA-seq data. Conclusion This study identified crucial mRNAs related to MM-persistent injury in endothelial cells even after switching the cells from high- glucose to normal glucose levels. Further research focusing on these mRNA may unravel new ways to modify MM in diabetes. Keywords: diabetes, high glucose, metabolic memory, messenger RNA, RNA-sequencing Introduction Type 2 diabetes mellitus (T2DM) is a chronic metabolic disease with a high blood glucose level leading to several complications associated with vascular disease, including vascular endothelial cell dysfunction, cardiomyopathy, and nephropathy.[38]^1^–[39]^3 Extensive research has shown that various factors play an important role in endothelial dysfunction in T2DM, including aging, obesity, hypertension, hyperlipidemia, hyperglycemia, low-grade inflammation, and insulin resistance.[40]^4^–[41]^6 The most recent investigations have been focused on effective approaches to prevent hyperglycemia-induced endothelial cell injury.[42]^7^–[43]^9 However, lowering glucose levels is not enough to turn off the oxidative damage in vascular endothelial cells induced by high glucose.[44]^10 This phenomenon can be partly explained by the “metabolic memory” (MM) theory, defined as the perpetuation of vascular damage despite the achievement of improved glycemic control.[45]^11 Some studies show that reactive oxygen species (ROS)-mediated vascular stress can maintain glucose normalization in DM animals and in cultured endothelial cells. Long-term maintenance of oxidative stress leads to the overactivation of pathways closely associated with DM and diabetic complications.[46]^10 However, previous investigations have failed to reveal a comprehensive profile of mRNA expression in MM induced by hyperglycemia in HUVECs, and we ignore the proteins that are the most relevant and could become targets of intervention. To comprehensively elucidate the MM molecular mechanisms, we applied an RNA-seq strategy. First, we performed RNA-seq in HUVEC cells under control, HG and MM conditions. After a series of bioinformatics analyses, up-HGMMGs and down-HGMMGs were identified. Then clusterprofiler, GO, and pathway annotation analyses of up-HGMMGs and down-HGMMGs were performed. Next, a PPI network was structured based on the STRING database, visualization was carried out by Cytoscape, and modules were identified by MCODE. Finally, some key HGMMGs were verified by qPCR. Materials and methods “Memory” experiments in cells HUVECs derived from healthy donors were purchased from the China Center for Type Culture Collection (CCTCC, Wuhan, China). Methods of cell culture and RNA preparation have been described in previous studies.[47]^12^–[48]^14 In general, the high glucose-induced HUVECs were cultured with 25 mM glucose for six days, the control HUVECs were cultured with 5 mM glucose and 20 mM mannitol for the same period, and the memory HUVECs were cultured with 25 mM glucose for three days, followed by 5 mM glucose and 20 mM mannitol for the next three days. After six days, HUVECs were collected, and each group was analyzed in triplicate. Total RNA preparation and qualification Total RNA was isolated and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions. The quantity and quality of the RNA samples in each group were determined by NanoDrop 2000 (Thermo, Wilmington, DE) and Agilent 2100 bioanalyzer (Agilent Technologies, USA). Library construction and RNA-sequencing TruSeq^® Stranded Total RNA Sample Preparation kit was used to produce libraries following the manufacturer’s instructions. The purified libraries were quantified by Qubit^® 2.0 Fluorometer (Life Technologies, USA) and Agilent 2100 bioanalyzer (Agilent Technologies, USA). The cluster was generated by cBot, the library being sequenced using the Illumina HiSeq 2500 (Illumina, USA). The sequencing was performed at Origin-Biotech Inc (Ao-Ji Biotech, Shanghai, China). Bioinformatics analysis FastQC was conducted for Quality control (QC) of RNA-Seq reads (v. 0.11.3) ([49]http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Trimming was performed by seqtk for known Illumina TruSeq adapter sequences, poor reads, and ribosome RNA reads ([50]https://github.com/lh3/seqtk). The trimmed reads were then mapped to the human reference genome (hg38) by the Hisat2 (version:2.0.4).[51]^15^,[52]^16 Stringtie (version:1.3.0) was performed for each gene counts from trimmed reads.[53]^16^,[54]^17 Gene counts were normalized by TMM,[55]^18 and FPKM was performed by Perl script.[56]^19 edgeR was performed for determining differentially expressed genes,[57]^20^,[58]^21 and threshold with p-value<0.05 and |fold change| >1.2.[59]^22^–[60]^24 Venny was applied to screen of up-HGMMGs (HG vs LG and MM vs LG comparisons were performed for up-regulated expression, and an HG vs MM comparison was performed for non-significant differential expression) and down-HGMMGs (HG vs LG and MM vs LG comparisons were performed for down-regulated expression, and an HG vs MM comparison was performed for non-significant differential expression). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched by clusterProfiler for a better understanding of the functions of the HGMMGs.[61]^25^,[62]^26 Protein-protein interaction (PPI) network construction and module analysis STRING is a database that provides comprehensive information about interactions between proteins, including prediction and experimental interactions data.[63]^27 In our study, the STRING tool was used to perform the PPIs among the DEGs, and interactions of a combined score of ≥0.4. Then, Cytoscape was performed to visualize the network.[64]^28 PPI network was used to filter modules based on the Molecular Complex Detection plug-in (MCOD) in Cytoscape[65]^29 with a standard set following degree cut-off=2, k-core=2, node score cut-off=0.2, and max depth=100. MCODE score ≥4 and node ≥10 were considered for functional enrichment analysis of the modules. Validation of differentially expressed mRNAs by qRT-PCR Total RNA was extracted and reversely transcribed into cDNA using SuperScriptTM III Reverse Transcriptase (Invitrogen). Six of the differentially expressed mRNAs related with high glucose-induced HUVECs were randomly selected to be quantified by FastStart Universal SYBR Green Master (Rox) with specific primers. PCRs were performed in triplicate according to the following temperature profile: denaturation at 95 °C for 10 min followed by an amplification composed of 40 cycles of 95 °C for 10 s and 60 °C for 1 min. The primers used for these amplifications are listed in [66]Table 1. Glyceraldehyde phosphate dehydrogenase (GAPDH) expression was used as an internal reference. The data were analyzed using the 2^−ΔΔCt method and the expression levels of each mRNAs were represented as fold change. Table 1. Primer sequences Gene Refseq Primer sequences GAPDH [67]NM_002046 F: CCTGGTATGACAACGAATTTG R: CAGTGAGGGTCTCTCTCTTCC CCNB1 [68]NM_031966 F: TACCTATGCTGGTGCCAGTG R: CAGATGTTTCCATTGGGCTT CDK1 [69]NM_001786 F: TAAGCCGGGATCTACCATACC R: TTTCATGGCTACCACTTGACC CCNA2 [70]NM_001237 F: TGAAGATGCCCTGGCTTTTA R: CACTCACTGGCTTTTCATCTTCT CDK6 [71]NM_001145306 F: TCCCAGGAGAAGAAGACTGG R: GGTCCTGGAAGTATGGGTGA CREB1 [72]NM_004379 F: ATCCGGGCCGTGAACGAAAGC R: CTGTGGCTGGGCTTGAACTGTCA RPL11 [73]NM_000975 F: AAAGGTGCGGGAGTATGAGTT R: TCCAGGCCGTAGATACCAATG [74]Open in a new tab Statistical analysis Although there were three sets of samples, pairs were used for the comparisons. Thus, the data were analyzed by Student’s t-test. Results are expressed as the mean ± SEM. Significance levels were set at the 0.05 threshold. Data analysis was carried out using Statistical Program for Social Sciences (SPSS) 22.0 software (SPSS, Chicago, IL, USA). Results Identification of HGMMGs To fully understand the multifaceted mechanism of HG-induced MM in HUVECs, we performed RNA-seq for the transcriptomes of LG, HG, and MM samples. By bioinformatics analysis, after filtration there were 516 up-HGMMGs and 210 down-HGMMGs. From the data in [75]Figure 1A, it can be seen that the 516 up-HGMMGs were screened from the intersection of 3568 up-regulated mRNAs (HG vs LG, fold-change >1.2, p-value<0.05), 1201 up-regulated mRNAs (MM vs LG, fold-change >1.2, p-value<0.05), and 10,742 non-significantly differentially expressed mRNAs (MM vs HG, p-value>0.05). Meanwhile, from the data in [76]Figure 1B, it can be seen that 210 down-HGMMGs were identified from the intersection of 2238 down-regulated mRNAs (HG vs LG, fold-change <0.833, p-value<0.05) and 708 down-regulated mRNAs (MM vs LG, fold-change <0.833, p-value<0.05), and 10,742 non-significantly differentially expressed mRNAs (MM vs HG, p-value>0.05). Hierarchical clustering of the HGMMGs can be seen in [77]Figure 2. Figure 1. [78]Figure 1 [79]Open in a new tab Transcriptome comparisons of the HG, LG, and MM groups. (A) up-HGMMGs (HG vs LG and MM vs LG comparisons were performed for up-regulated expression, and an HG vs MM comparison was performed for non-significant differential expression) (B) down-HGMMGs (HG vs LG and MM vs LG comparisons were performed for down-regulated expression, and an HG vs MM comparison was performed for non-significant differential expression). Abbreviations: HG, high glucose; LG, normal glucose; MM, metabolic memory; HGMMGs, HG-induced metabolic memory-involved regulated genes. Figure 2. [80]Figure 2 [81]Open in a new tab Heatmap of down-HGMMGs and up-HGMMGS with a total of 726 RNAs. LG1-LG3 refers to the control HUVECs; HG1-HG3 refers to the HG-induced HUVECs; MM1-MM3 refers to the HG-induced MM HUVECs. Abbreviations: HG, high glucose; LG, normal glucose; MM, metabolic memory. GO analysis of HGMMGs To further understand the genes associated with HG-induced Metabolic Memory in HUVECs, GO enrichment analysis was performed by clusterProfiler on the 516 up-HGMMGs and 210 down-HGMMGS. Enrichment analysis showed that a total of 421 GO-terms were significantly enriched in the 516 up-HGMMGs, 326 of which were biological processes (BP), 56 of which were cellular components (CC), and 39 of which were molecular functions (MF) with p<0.05. Mitosis, mitotic cell cycle, mitotic cell cycle process, cell division, cell cycle, chromosome segregation, cell cycle process, DNA packaging, chromosome organization, and organelle organization were the top10 enriched BP terms. The top10 enriched CC terms were chromosome, spindle, nucleosome, condensed chromosome, centromeric region, condensed chromosome kinetochore, chromosome, centromeric region, nuclear chromosome, protein-DNA complex, kinetochore, and condensed chromosome. The enriched MF terms were microtubule binding, tubulin binding, cytoskeletal protein binding, GTPase regulator activity, microtubule motor activity, protein heterodimerization activity, macromolecular complex binding, protein complex binding, GTPase activator activity, and kinase binding. These results show that cell proliferation is unbalanced, especially in chromosomal replication, segregation, and organization. The top30 GO-terms with highest enrich factor are shown in [82]Figure 3A. Figure 3. [83]Figure 3 [84]Open in a new tab Gene Ontology enrichment analysis for MM-maintained up-HGMMGs. Results of gene ontology enrichment analysis for MM-maintained up-HGMMGs (A) and down-HGMMGs (B) in HG-induced HUVECs. The 210 down-HGMMGs were significantly enriched in 83 GO-terms, of which 59, 12, and 12 were BP, CC, and MF (with p<0.05). The top10 enriched GO biological processes were detection of chemical stimulus involved in sensory perception of bitter taste, sensory perception of bitter taste, detection of chemical stimulus involved in sensory perception of taste, nuclear-transcribed mRNA catabolic process, mRNA catabolic process, regulation of cell size, mitochondrial ATP synthesis coupled electron transport, ATP synthesis coupled electron transport, RNA catabolic process, and T cell activation involved in immune response. Meanwhile, motile cilium, cell projection, axon, cell leading edge, mitochondrial respiratory chain, respiratory chain, cilium, cytosolic ribosome, sarcolemma, cell-cell junction, neuron projection, and cell junction were the enriched CC. The enriched MF were bitter taste receptor activity, taste receptor activity, oxidoreductase activity, acting on NAD(P)H, quinone or similar compound as acceptor, oxidoreductase activity, acting on NAD(P)H, calcium ion binding, growth factor receptor binding, structural constituent of ribosome, ubiquitin protein ligase activity, inorganic cation transmembrane transporter activity, protein binding, bridging, RNA polymerase II core promoter proximal region sequence-specific DNA binding, and transmembrane transporter activity. The top30 GO-terms with highest enrich factor are shown in [85]Figure 3B. KEGG analysis of HGMMGs Pathway enrichment analysis could provide further insights into the function of genes and their interaction. We performed the KEGG pathway enrichment analysis for up-HGMMGs, and found 221 pathway terms including 18 pathway terms with a p-value<0.05. The Top10 pathways with the greatest enrichment were Systemic lupus erythematosus (hsa05322), oocyte meiosis (hsa04114), cell cycle (hsa04110), alcoholism (hsa05034), progesterone-mediated oocyte maturation (hsa04914), p53 signaling pathway (hsa04115), viral carcinogenesis (hsa05203), glycosaminoglycan biosynthesis - chondroitin sulfate/dermatan sulfate (hsa00532), sulfur metabolism (hsa00920), and HTLV-I infection (hsa05166). The top30 enriched pathways are presented in [86]Figure 4A. Figure 4. [87]Figure 4 [88]Open in a new tab KEGG enrichment analysis for MM-maintained up-HGMMGs. Results of KEGG enrichment analysis for MM-maintained up-HGMMGs (A) and down-HGMMGs (B) in HG-induced HUVECs. KEGG enrichment analysis showed that a total of 143 pathway-terms were enriched with 210 down-HGMMGs, 12 of which were significant (p<0.05), including: oxidative phosphorylation (hsa00190), parkinson’s disease (hsa05012), glutathione metabolism (hsa00480), ribosome (hsa03010), mucin type O-Glycan biosynthesis (hsa00512), RIG-I-like receptor signaling pathway (hsa04622), chronic myeloid leukemia (hsa05220), insulin secretion (hsa04911), estrogen signaling pathway (hsa04915), glioma (hsa05214), epithelial cell signaling in Helicobacter pylori infection (hsa05120), and drug metabolism - cytochrome P450 (hsa00982). The top30 pathway-terms are shown in [89]Figure 4B. PPI network The significant HGMMGs were used to construct PPI network. The PPI network consists of 462 nodes and 2,656 interactions, as shown in [90]Figure 5. Dozens of gene nodes were high in connectivity degrees, including: CDK1 (cyclin-dependent kinase 1, degree=94, up-HGMMGs), CCNB1 (cyclin A2, degree=86, up-HGMMGs), CCNA2 (cyclin A2, degree=76, up-HGMMGs), AURKA (Aurora Kinase A, degree=76, up-HGMMGs), BUB1 (BUB1 mitotic checkpoint serine/threonine kinase, degree=74, up-HGMMGs), CDK6 (Cyclin-Dependent Kinase 6, degree=25, down-HGMMGs), RPL11 (Ribosomal Protein L11, degree=19, down-HGMMGs), and CREB1 (CAMP Responsive Element Binding Protein 1, degree=19, down-HGMMGs), among others. Details on the Top10 with the highest degree of down-HGMMGs and up-HGMMGs are included in [91]Table 2. Figure 5. [92]Figure 5 [93]Open in a new tab PPI network of the MM-maintained HGMMGs. The size of the nodes is positively correlated to the nodes’ degree; pink nodes denote the up-regulated genes, while green nodes denote the down-regulated genes. Abbreviations: PPI, protein–protein interaction; HGMMGs, HG-induced metabolic memory-involved genes. Table 2. Top10 of down-HGMMGs and up-HGMMGs details Gene_id Gene_name MM vs LG MM vs HG HG vs LG FC Pvalue Updown FC Pvalue Updown FC Pvalue Updown ENSG00000167083 GNGT2 0.255 0.004 DOWN 0.946 1.000 - 0.270 0.013 DOWN ENSG00000118260 CREB1 0.528 0.001 DOWN 0.892 0.266 - 0.592 0.014 DOWN ENSG00000063046 EIF4B 0.615 0.000 DOWN 0.983 0.740 - 0.626 0.001 DOWN ENSG00000162244 RPL29 0.718 0.027 DOWN 1.041 0.957 - 0.689 0.009 DOWN ENSG00000188846 RPL14 0.719 0.001 DOWN 0.956 0.445 - 0.751 0.012 DOWN ENSG00000105810 CDK6 0.739 0.000 DOWN 1.026 0.915 - 0.720 0.000 DOWN ENSG00000142676 RPL11 0.754 0.014 DOWN 0.989 0.680 - 0.762 0.027 DOWN ENSG00000115268 RPS15 0.775 0.030 DOWN 1.074 0.710 - 0.721 0.009 DOWN ENSG00000198727 MT-CYB 0.798 0.050 DOWN 1.119 0.458 - 0.713 0.005 DOWN ENSG00000167110 GOLGA2 0.803 0.008 DOWN 0.995 0.680 - 0.807 0.029 DOWN ENSG00000087586 AURKA 1.305 0.032 UP 0.881 0.175 - 1.481 0.002 UP ENSG00000145386 CCNA2 1.679 0.000 UP 1.002 0.726 - 1.675 0.000 UP ENSG00000169679 BUB1 1.889 0.000 UP 0.952 0.426 - 1.983 0.000 UP ENSG00000164109 MAD2L1 1.928 0.000 UP 1.236 0.095 - 1.560 0.000 UP ENSG00000080986 NDC80 1.950 0.000 UP 1.018 0.881 - 1.916 0.000 UP ENSG00000117399 CDC20 2.096 0.000 UP 0.929 0.220 - 2.256 0.000 UP ENSG00000157456 CCNB2 2.107 0.000 UP 1.119 0.444 - 1.883 0.000 UP ENSG00000134057 CCNB1 2.335 0.000 UP 1.137 0.289 - 2.053 0.000 UP ENSG00000170312 CDK1 2.464 0.000 UP 1.103 0.569 - 2.233 0.000 UP ENSG00000115163 CENPA 3.416 0.000 UP 1.290 0.069 - 2.647 0.000 UP [94]Open in a new tab A total of 17 modules were obtained using default criteria by MCODE. Modules were listed in descending order by MCODE score ([95]Table 3). Four modules with MCODE score ≥3 were named as module 1, module 2, module 3, and module 4. These four modules were selected for module network visualization ([96]Figure 6A–[97]D). Most of HGMMGs belonging to the four modules were up-regulated in the HG and MM groups. Table 3. List of modules exhibiting the HGMMG-related PPI Cluster Score Nodes Edges Node IDs 1 22.102 49 1083 DEPDC1, CDC20, CCNA2, CDK1, BUB1B, BUB1, NUF2, CCNB2, CDC25C, MAD2L1, RACGAP1, KIF23, NDC80, SPC25, PTTG1, CCNB1, CKS2, BIRC5, CENPE, AURKA, NEK2, CENPF, ECT2, CENPA, KIF2C, NCAPG, HMMR, SKA1, KIF4A, KIF15, CDCA5, SKA3, OIP5, KIF20B, PBK, GTSE1, NUSAP1, NCAPG2, SHCBP1, ANLN, RAD51AP1, KIAA0101, TRIP13, SPAG5, TROAP, FAM64A, POLQ, ARHGAP11A, FAM83D 2 4.5 12 54 HIST1H4J, HIST3H2A, HIST1H2BJ, HIST2H2BF, HIST3H2BB, BRCA2, HIST1H2BD, H3F3B, RBBP7, HIST1H3D, HIST2H2AA4, HIST2H3C 3 4.364 11 48 NMU, MAPK1, GNGT2, TAS2R13, TAS2R4, S1PR2, TAS2R14, TAS2R10, DRD4, OPRD1, ADCY7 4 3.5 8 28 UBE2S, ANAPC11, RNF6, TRAF7, FBXW5, TRAIP, RNF217, FBXL15 5 2.556 18 46 EIF4B, RPS15, BTF3, RPL29, UPF1, RPL14, PELP1, SRC, EPOR, SCN11A, RET, MRPL49, MRPL34, MRPS6, RPL11, CDKN2C, CCND3, EIF3J 6 2.333 6 14 NDUFA4, MT-ND3, MT-ND2, MT-ND4L, MT-ATP8, MT-CYB 7 1.5 4 6 CERK, PPAP2C, CERS1, UGCG 8 1.4 5 7 BORA, CKAP5, KIF22, NDE1, SGOL1 9 1.4 5 7 GALNTL1, GALNTL2, MUC20, B4GALT5, MUC16 10 1 3 3 RPS2, IMP3, PUSL1 11 1 3 3 NT5M, ITPA, APRT 12 1 3 3 TERT, MUS81, PIF1 13 1 3 3 CORO2A, KRT7, SCARF2 14 1 3 3 EFCAB6, KCNN3, KCNN4 15 1 3 3 MCRS1, GPS1, DDB2 16 1 3 3 HOMER3, SHANK2, SHANK3 17 1 3 3 NABP1, INTS5, NABP2 [98]Open in a new tab Figure 6. [99]Figure 6 [100]Open in a new tab PPI network of the MM-maintained HGMMGs involved in Module 1 (A), 2 (B), 3 (C) and 4 (D). Pink nodes denote the up-regulated genes, while green nodes denote the down-regulated genes; PPI, protein–protein interaction; HGMMGs, HG-induced metabolic memory-involved genes. The results of qRT-PCR verification of the HGMMGs As shown in [101]Figure 7, real-time PCR results revealed that CDK1, CCNB1, and CCNA2 were over-expressed in the HG group compared to the LG group. Meanwhile, CDK6, RPL11, and CREB1 showed lower-expression in the HG group compared to the LG group ([102]Figure 7, black bars). In the MM group, these genes were expressed at different levels ([103]Figure 7, gray bars). In addition, except for CCNB1, these genes were not significantly differentially expressed in the MM group compared with the HG group ([104]Figure 7, white bars). Overall, these results were consistent with the RNA-seq data. Figure 7. [105]Figure 7 [106]Open in a new tab Verification of the HGMMGs by qRT-PCR. The expression of six genes in HUVECs was detected by qRT-PCR and the data are shown as fold changes. GAPDH was used as control. Black bars denote the comparison of HG vs LG, and the # symbol denotes the p-value; Gray bars denote the comparison of MM vs LG, and the * symbol denotes the p-value; White bars denote the comparison of MM vs HG, and the $ symbol denotes the p-value; one, two, and three symbols mean 0.01