Construction of lncRNA-related ceRNA regulatory network in diabetic subdermal endothelial cells

ABSTRACT Long non-coding RNAs (lncRNAs) were considered to be involved in vascular complications in diabetes mellitus, but still only limited knowledge in this regard has been obtained. Herein, we further explored the roles of lncRNAs and mRNAs in diabetic vasculopathy (DV) through conducting bioinformatics analysis using data set downloaded from GEO database. The differentially expressed lncRNAs and mRNAs were identified by edge package. GO enrichment analysis and KEGG pathway analysis were performed based on clusterprofiler package. The relationship between lncRNA and miRNA was predicted using starBase database, and the potential mRNAs targeted by miRNAs were predicted by TargetScan, miRTarbase and miRDB database. The string database was used to analyze the protein-protein interaction (PPI). As a result, a total of 12 lncRNAs and 711 mRNAs were found to be differentially expressed in the diabetic subdermal endothelial cells compared with normal controls. A ceRNA network was established, which was composed of seven lncRNA nodes, 49 miRNA nodes, 58 mRNA nodes and 183 edges, and MSC-AS1 and LINC01550 may serve as key nodes. GO function enrichment analysis showed enrichments of epithelial cell proliferation, intercellular junction, and cell adhesion molecule binding. KEGG pathway analysis revealed 33 enriched pathways. PPI protein interaction analysis identified 57 potential ceRNA-related proteins. Overall, this study suggests that multiple lncRNAs, specifically MSC-AS1 and LINC01550, may play an important role in DV development and they are like to be developed as the therapeutic targets for DV. However, further experiments in vitro and in vivo should be conducted to validate our results.


Introduction
Diabetes mellitus is a worldwide disease with high incidence and morbidity, and it is characterized with chronic hyperglycemia caused by multiple causes. The pre-diabetes changes and diabetes mellitus (DM) have affected 425 million people worldwide in 2017, and it was predicted to rise to 629 million by 2045 [1,2]. The long-term increase of blood glucose can lead to the damage of large blood vessels and micro vessels and endanger the heart, brain, kidney, peripheral nerve, eyes and feet. According to the statistics of the World Health Organization, there are more than 100 kinds of complications [3,4]. Diabetic vasculopathy (DV) is one of the most common and severe complications of both type 1 and type 2 diabetes mellitus. DV refers to a series of vascular damages and is related to many other complications, such as diabetic retinopathy, diabetic nephropathy and diabetic foot ulcer (DFU). DV and its relevant complications have been considered to influence the prognosis of patients with DM, so it is necessary to fully elucidate its mechanisms, so as to develop strategies to prevent or deal with it.
The damage of vascular endothelial cell is the most crucial step during DV development. Under physiological condition, insulin can interact with its receptor on endothelial cells to activate eNOS, upregulates NO release and reduces ET-1 secretion, consequently inducing the vasodilation of blood vessels [5]. However, the insulin and glucose uptake by endothelial cells decrease and even insulin resistance occurs ultimately when endothelial cells are damaged [6]. Additionally, oxidative stress has also been proposed to promote DV through overproduction of reactive oxygen species (ROS), which can induce endothelial cell apoptosis by negatively regulating Akt/NF-kappaB pathway [7]. Furthermore, evidence exists that mitochondrial energy metabolism dysfunction and advanced glycation end products (AGEs) participate in DV development as well, since they are able to disturb the physiological stability of endothelial cells, thereby inducing endothelial cell death [8,9]. To date, our knowledge of the damage of vascular endothelial cells in DM is still insufficient, and much effort should be paid to investigate its pathogenesis and further determine its roles in DV development.
Long non-coding RNA (LncRNA) is a class of noncoding RNA molecules with longer than 200 nucleotides that lack protein-coding functions. Recently, increasing studies suggested that multiple lncRNAs are dysregulated in DM and involved in the pathogenesis of DM. For instance, Ruan et al. revealed that lncRNA-p3134 contributes to glucose metabolism and insulin signal transduction in pancreatic B cells [10]. Feng et al. demonstrated that lncRNAs modulate inflammation in diabetic peripheral neuropathy as ceRNAs targeting miR-146a-5p [11]. Kong et al. reported that lncRNA LEGLTBC regulates glucolipotoxicityinduced INS-1 beta cell oxidative stress and apoptosis via the miR-34a/SIRT1 axis [12]. Nevertheless, there are few studies exploring the role of lncRNAs in DV.
To develop new strategies for preventing and treating DV, in this study, we tried to investigate the role of lncRNAs and their potential mechanisms in DV by the bioinformatics analysis using a high throughput sequencing dataset (GSE92724) and miRNA-predicted websites. As a result, we identified 12 differentially expressed lncRNAs in the diabetic subdermal endothelial cells. Then, the ceRNA regulatory network was established based on those lncRNAs and two key ceRNA subnetworks may be involved in the pathogenesis of DV. Overall, this study helps us to better understand the effects of lncRNAs in the dermal endothelial cells whose damage is involved in the pathogenesis of DV.

Data acquisition
The expression of dermal endothelial cells of type 2 diabetic patients processed by high throughput sequencing was downloaded from the Gene Expression Omnibus (GEO) database (https:// www.ncbi.nlm.nih.gov/geo/). The GSE92724 dataset included 4 diabetic patients and 6 control individuals. The mRNAs and lncRNAs were distinguished by using the annotation information obtained from Genecode database. Totally, 18,206 mRNAs and 2025 lncRNAs were obtained from the high throughput sequencing dataset.

Functional enrichment analysis of differentially expressed mRNAs and lnRNAs
The differentially expressed mRNAs and lnRNAs of dermal blood endothelial cell between diabetic patients and the healthy controls were screened using 'edgeR' [13] package in the R software (version 3.6.1), based on the cutoff values: |logFC| > 1 and the false discovery rate (FDR) < 0.05. The R package 'clusterProfiler' [14] was used to implement the functional annotation of Gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the differentially expressed ARGs. Adjusted P < 0.05 was considered statistically significant.

Protein-Protein Interaction (PPI) network analysis of mRNAs in ceRNA network
The PPI network of mRNAs involved in the ceRNA network was established through an online search tool, Search Tool for the Retrieval of Interacting Genes (STRING) database (https:// string-db.org) with a confidence score > 0.40 [21].

Conduction of Key ceRNA subnetwork
Network degree of each node in the ceRNA network was evaluated in Cytoscape. The lncRNA with node degree ≥6 was considered as a hub lncRNA. Then, the hub lncRNA-associated subnetworks were established by using the information between the relationship of hub lncRNA-miRNA and the relationship of miRNA-mRNA.

Identification of differentially expressed mRNAs and lnRNAs
In order to facilitate the understanding of the research, the entire workflow is presented in Figure 1. Raw data were downloaded from GSE92724 dataset in GEO database and 8206 mRNAs and 2025 lncRNAs were identified by using the annotation information downloaded from Genecode database. The edgeR package was performed to identify differentially expressed mRNAs and lncRNAs with a strict cutoff threshold of |logFC | > 1 and FDR <0.05. Compared with normal samples, 12 DElncRNAs and 711 DEmRNAs were differentially expressed, among which 9 lncRNAs and 614 mRNAs were downregulated as well as 3 lncRNAs and 97 mRNAs were upregulated in diabetic dermal blood endothelial cells (Table 1 and Supplementary  Table S1). The volcano plots, heat maps and PCA plot of DElncRNAs and mRNAs of dermal blood endothelial cell between diabetic patients and the healthy controls are shown in Figure 2.

Functional enrichment analysis of DEmRNAs
To determine the functions of the significantly upregulated and downregulated DEmRNAs, GO enrichment analysis and KEGG pathway analysis were conducted by using 'clusterProfiler' R package. The BP analysis showed that the DEmRNAs were mostly enriched in epithelial cell proliferation, regulation of epithelial cell proliferation, and epithelial cell development. The CC analysis showed that the DEmRNAs were significantly enriched in cell-cell junction, collagen−containing extracellular matrix, and basolateral plasma membrane. The MF analysis showed that the DEmRNAs were significantly enriched in cell adhesion molecule binding, integrin binding and glycosaminoglycan binding. The KEGG pathways enriched analysis results showed that major pathways, including the PI3K−Akt signaling pathway, Cytokine−cytokine receptor interaction, MAPK signaling pathway and MAPK signaling pathway, were involved in the DEmRNAs. The results of the GO and KEGG functional enrichment analyses are presented in Figure 3. The complete list of the GO and KEGG analyses is presented in Supplementary Table S2-S3.

Construction of ceRNA network
To establish the ceRNA network, the DElncRNAs available in the starBase database were further analyzed and the subcellular localizations of these lncRNAs were predicted in the lncLocator database. All of the 7 DElncRNAs predicted were located outside the nuclear ( Table 2). The potential DElncRNA targeted miRNAs were predicted using starBase database. Then, the DEmRNA-miRNA pairs were evaluated by the Targetscan, miRTarbase, and miRDB databases. When the DElncRNA-miRNA pairs and the DEmRNA-miRNA pairs contained a common miRNA, they were chosen to conduct the ceRNA network. In total, the ceRNA network was conducted with 183 edges and 114 nodes, including 7 lncRNAs, 49 miRNAs and 58 mRNAs (Figure 4).

Functional enrichment analysis of the ceRNA Network
To gain further insight into the functions of the targeted mRNAs in ceRNA network, GO enrichment analysis and KEGG pathway analysis were conducted ( Figure 5). The complete list of the GO and KEGG analyses is presented in Supplementary Tables S4-S5. The BP analysis showed that the targeted mRNAs were mainly enriched in extracellular structure organization, regulation of actin cytoskeleton organization, and extracellular matrix organization. The CC analysis showed that the targeted mRNAs were mainly enriched in collagen−containing extracellular matrix, basolateral plasma membrane, and basolateral plasma membrane. The MF analysis showed that the targeted mRNAs were mainly enriched in extracellular matrix structural   Figure 5.

Conduction of Key ceRNA subnetwork
To figure out the key lncRNA-associated ceRNA subnetwork, the node degree of each gene in this ceRNA network was evaluated by using a built-in 'cytoHubba' tool of Cytoscape. As listed in Table 3, a total of 14 genes with node degree no less than six were disclosed, including 2 lncRNAs (MSC-AS1, LINC01550). Then, the interactions of hub lncRNAs were extracted to conduct the related ceRNA subnetworks. As shown in Figure 6, the MSC-AS1-associated ceRNA sub-network was composed of 25 miRNAs, 26 mRNAs and 85 edges. The LINC01550-associated ceRNA subnetwork was composed of 13 miRNAs, 14 mRNAs and 49 edges.

Discussion
Diabetic vasculopathy is a defect in blood vessels caused by hyperglycemia, which causes severe organ-specific complications, such as nephropathy, retinopathy, myocardial infarction and cerebrovascular disease [22,23]. Evidence shows that the damage of vascular endothelial cells plays a crucial role in the occurrence and development of DV. Thus, to comprehensively understand the molecular mechanisms for the damage of vascular endothelial cells may help to find new strategies for the prevention and therapy of DV. In recent years, lncRNAs were found to regulate vascular endothelial cell functions, indicating their potential roles in DV development [24]. For example, lncRNA-MIAT can regulate microvascular dysfunction by functioning as a ceRNA [25]. In addition, Gong et al. suggested that LncRNA TDRG1 aggravates retinal microvascular endothelial cell dysfunction in diabetic retinopathy through upregulating VEGF [26]. With the development of  high-throughput sequencing, a large number of novel lncRNA were identified in DV. However, the association of these lncRNAs with DV development remains largely unknown.
In this study, we tried to further identify lncRNAs involved in DV and preliminarily explore the potential mechanisms by conducting the bioinformatics analysis of publica dataset (GSE92724). A total of 12 differently expressed lncRNAs in subdermal endothelial cells were identified between 4 diabetic patients and 6 control individuals. Then, a ceRNA network was finally established with 183 edges and 114 nodes, including 7 lncRNAs, 49 miRNAs and 58 mRNAs. Additionally, KEGG pathway analysis further revealed that 33 signaling pathways were enriched, such as PI3K/Akt pathway, Rap1 pathway and MAPK pathway.   [29]. Consistently, by GO analysis, we found that mRNAs in the ceRNA network were mainly enriched in extracellular structure organization, collagen-containing extracellular matrix and extracellular matrix structural constituent. Among seven differentially expressed lncRNAs, MSC-AS1 and LINC01550 that were markedly downregulated in diabetic subdermal endothelial cells were found to serve as key nodes. To date, no studies were performed to explore the roles of MSC-AS1 and LINC01550 in DM. However, numerous recent studies suggested that MSC-AS1 is dysregulated in multiple types of cancers and performs an oncogenic function in tumor progression [30,31], whereas LINC01550 was rarely studied. As mentioned above, our KEGG pathway analysis suggested that the differentially expressed lncRNAs may regulate DV development via the PI3K/Akt pathway. Coincidentally, dysregulation of MSC-AS1 expression was found to promote tumor progression in glioma and osteosarcoma through the PI3K-Akt pathway [32,33], which may further support the important role of MSC-AS1/PI3K-Akt pathway in DV to a certain extent.
There are some limitations in this study, and it should be seriously considered when interpreting our results. To begin with, the conclusion of this study was merely made based on the bioinformatic analysis. Therefore, in future, necessary experiments in vitro and in vivo should be performed to further validate our findings. Additionally, the publica data we used only included 4 diabetic patients and 6 control individuals, so the limitation of small sample size may lead to overestimated conclusion. Hence, large sample and multi-center clinical studies should be conducted to confirm our results.

Conclusion
In summary, 12 lncRNAs and 711 mRNAs were identified to differentially express in the diabetic subdermal endothelial cells versus normal controls. Based on these lncRNAs and mRNAs, a lncRNA-miRNA-mRNA interaction network was established, which was composed of seven lncRNA nodes, 49 miRNA nodes, 58 mRNA nodes and 183 edges, and MSC-AS1 and LINC01550 were found to serve as key nodes in the network. Overall, this study suggested that multiple lncRNAs, specifically MSC-AS1 and LINC01550, may play an important role in DV development and they are like to be developed as the therapeutic targets for DV. However, further experiments in vitro and in vivo should be conducted to validate our results.

Disclosure of potential conflicts of interest
No potential conflict of interest was reported by the author(s).