A comprehensive competitive endogenous RNA network pinpoints key molecules in diabetic retinopathy.

Diabetic retinopathy (DR) is a severe microvascular complication of diabetes and the primary cause of vision loss in diabetic patients. Previous research has revealed that long non‑coding RNAs (lncRNAs) and microRNAs (miRNAs) play pivotal roles in the pathogenesis of DR. However, the roles of lncRNA‑miRNA‑mRNA interactions in DR are poorly understood. In the present study, we aimed to compute a global triple network of competitive endogenous RNAs (ceRNAs) in order to pinpoint essential molecules. We found that there were 802 nodes (121 lncRNA nodes, 17 miRNA nodes, and 664 mRNA nodes) and 949 edges in the ceRNA network. Further functional analysis suggested that some molecules were specifically related to DR. Surprisingly, these molecules were involved in visual perception, eye development, and lens development in camera‑type eye. In summary, our study highlighted specific lncRNAs and miRNAs related to the pathogenesis of DR, which might be used as potential diagnostic biomarkers and therapeutic targets for DR.

Introduction

Diabetic retinopathy (DR) is a severe microvascular complication of diabetes. It is the main cause of loss of vision in diabetic patients (1). Diabetic patients are known to suffer from hyperglycemia due to abnormalities in insulin signaling pathways; these abnormal levels of sugar can cause various pathogenetic processes, resulting in the following complications in diabetic patients: i) microvascular defects of the retina and ii) dysfunction and degeneration of the neuro-retina (2). The progression of DR depends upon two major types of diabetic mellitus. Research has shown that DR is the most common vision-threatening lesion in patients with type I diabetes. Diabetic macular edema seems to be more prevalent in patients with type II diabetes (3). Several studies have tried to determine the causes and progression of DR in diabetic patients (4–6); however, scientists have not yet been successful in completely elucidating the mechanisms associated with the pathogenesis of DR. Currently, diabetes mellitus cannot be prevented or treated effectively in clinical practice.

In the field of molecular biology, scientists have made major breakthroughs in understanding non-coding RNAs (ncRNAs). The prognosis of diabetic patients has improved tremendously because of various path-breaking studies conducted by researchers in the field of molecular biology (7). Recent studies have reported that long non-coding RNAs (lncRNAs) are involved in various biological processes (7), and it was found that lncRNAs are widely involved in signaling pathways, which regulate numerous aspects of life processes (7). In the context of oncology, pan-cancer molecular portraits of human tumors have been associated with essential lncRNAs (8). The activity of essential oncogenic pathways was altered with dysregulation in the quantities of lncRNAs (8). lncRNAs impact post-translation processes partly based on the competitive endogenous RNA (ceRNA) mechanism, where lncRNA interacts with microRNA (miRNA) through miRNA-binding sites (MREs) and hence regulate the expression of certain genes (9). As a result, ceRNA regulation is highly associated with several pathogenic processes.

lncRNAs and miRNAs were previously demonstrated to be associated with the pathogenesis of DR (10–12). Previous findings highlight that the development of DR is associated with various key molecules, such as miR-1273g-3p, miR-21 and miR-29a (13–15). Research has elucidated how miRNAs play a pivotal role in the diagnosis of DR (16). For example, Yan et al found that the expression of lncRNAs was aberrant in the early stages of DR. In other words, DR may develop as follows. The aberrantly expressed lncRNAs play pivotal roles in modulating target molecules, and they regulate multiple pathogenetic pathways (11). Presently, scientists have reported that disease-related lncRNAs are associated with function-related messenger RNAs (mRNAs) and miRNAs (17,18). In many other diseases except DR, ceRNA mechanism have been investigated to date; this has radically transformed the approach used to understand the pathogenesis of several diseases (19), including cardiac hypertrophy (20), rheumatoid arthritis (21) and papillary renal cell carcinoma (22). However, the role played by the lncRNA-miRNA-mRNA network in the pathogenesis of DR remains unclear. We have preliminary knowledge of the function and physiological significance of ceRNAs.

Here, based on ceRNA theory, we attempted to explore the crosstalk between lncRNA-miRNA-mRNA complexities. Key molecules were identified by performing Gene Ontology (GO) analysis and Weighted Gene Co-expression Network Analysis (WGCNA) (Fig. 1). In this study, ceRNA theory was used for the first time in DR to compute a systematic profile of the lncRNA-miRNA-mRNA network. Our results may be useful in comprehending the role played by ceRNAs in the pathogenesis of DR.

Figure 1.

The reconstruction pipeline of the lncRNA-miRNA-mRNA network. lncRNA, long non-coding RNA; miRNA, microRNA; mRNA, messenger RNA; ceRNA, competitive endogenous RNA.

Materials and methods

Data availability

Using eight week-old C57BL/6 mice, a mouse model of streptozotocin (STZ)-induced diabetes was constructed. Two months after administering STZ injection to diabetic mice, total RNAs were isolated from the retinas of the diabetic mice. The mice were considered to be diabetic when their blood glucose levels exceeded 250 mg/dl. Using TRIzol reagent, wild-type mice were matched according to age and sex. Agilent Mouse Gene Expression Microarrays (Product Number G4852A; Agilent Technologies, Santa Clara, CA, USA) were used to generate microarray data from mouse mRNA and lncRNA. The raw data were provided by Yan et al, which has been described in a previous study (11), and the ethics statement was previously provided. To generate microarray data of miRNAs from STZ-induced diabetic rats, we referred to a previous study and its supplementary material (12).

Screening of differentially expressed lncRNAs, miRNAs and mRNAs

As for the differential expression analysis, the Benjamini-Hochberg method was applied to calculate the adjusted P-value to minimize the false discovery rate (FDR) as was previously described (12). |Fold change|>2 and P<0.05 were set as cut-offs.

Prediction of target lncRNAs and mRNAs of miRNAs

The target lncRNAs of miRNAs was predicted by referring to Miranda (23). The minimum free energy (MFE) was calculated under the following condition: Max energy ≤-20 and score >160 were considered as cut-offs. The target mRNAs of miRNAs was predicted on the basis of miRTarBase (24), miRecords (25) and starBase version 2.0 (26).

Construction of the ceRNA network

The lncRNA-miRNA-mRNA network was reconstructed with ceRNA theory. The expression correlation of lncRNAs and mRNAs was determined with Pearson correlation coefficient (PCC). The mRNA-lncRNA pair was considered as the co-expressed one (P<0.05). The ceRNA network was constructed by assembling all the co-expressed competing triplets, and this network was visualized with Cytoscape software v3.5.1 (https://cytoscape.org/).

Function enrichment analysis

While performing GO analysis, these databases were used for annotation and visualization. Integration Discovery (DAVID) and Cytoscape plug-in BinGO were used to complete the GO analysis. To perform WGCNA analysis, we used R WGCNA package (https://cran.r-project.org/package=WGCNA). The cut height value was set to 0.992. We only considered modules in which the distance between two consecutive modules was <0.1. All the modules were then combined into a single module. The heat map was drawn by using pheatmap R package (https://cran.r-project.org/web/packages/pheatmap/).

Results

Screening of lncRNAs, mRNAs and miRNAs is associated specifically with DR

We compared expression profiles of lncRNAs, miRNAs and mRNAs, which were obtained from the retinas of different diabetic mice. To determine the key RNA molecules, differential expression analysis was performed (2). In this study, 305 differentially expressed lncRNAs and 17 differentially expressed miRNAs (diabetic vs. non-diabetic ones) were observed. Among them, 89 lncRNAs and 14 miRNAs showed upregulated expression, whereas 214 lncRNAs and 3 miRNAs showed downregulated expression in the diabetic retina samples compared with the controls (Fig. 2). The occurrence and development of DR was regulated by the differentially expressed genes.

Figure 2.

(A) Heat map of the differentially expressed lncRNAs for the DR group and control group. (B) Heat map of differentially expressed mRNAs for the DR group and control group. All the values were normalized to t values. In diabetic retina samples, 89 lncRNAs and 14 miRNAs were upregulated, whereas 214 lncRNAs and 3 miRNAs were downregulated. lncRNAs, long non-coding RNAs; DR, diabetic retinopathy; mRNAs, messenger RNAs; miRNAs, microRNAs.

Construction of lncRNA-miRNA and mRNA-miRNA co-expression network

lncRNAs can competitively interact with miRNAs, and they may function as ceRNAs. With the help of miRanda, the interactions between lncRNAs and miRNAs were predicted. The interactions were predicted by referring to the minimal free energy of miRNA-lncRNA duplexes (max energy ≤-20 and score >160). In total, we predicted 246 pairs of lncRNA-miRNA; 17 miRNAs and 121 lncRNAs were involved in these 246 pairs (Fig. 3A). Using miRTarBase, miRecords and starBase version 2.0, we identified the target mRNAs of miRNAs. Thus, we identified 664 target mRNAs of 14 miRNAs (Fig. 3B).

Figure 3.

(A) The lncRNA-miRNA co-expression network. In total, 246 lncRNA-miRNA pairs were predicted; 17 miRNAs and 121 lncRNAs were involved in 246 pairs. lncRNA, long non-coding RNA; miRNA, microRNA. (B) The mRNA-miRNA co-expression network. The computations of miRTarBase, miRecords, and starBase pinpointed the target mRNAs of miRNAs; moreover, we identified 664 target mRNAs of 14 miRNAs. miRNA, microRNA; mRNA, messenger RNA.

Construction of the competitive endogenous RNA network

Since lncRNAs were found to be specifically associated with DR, we further investigated them with Cytoscape software. Then, the results of the analyses were used to construct the ceRNA network. Based on previous results of the mRNA-miRNA and lncRNA-miRNA co-expression network, the interactions between the lncRNA-miRNA-mRNA network were visualized (Fig. 4). In total, 802 nodes (121 lncRNA nodes, 17 miRNA nodes, and 664 mRNA nodes) and 949 edges were found in this network. The key miRNAs were as follows: miR-223-3p, miR-34c-5p, and miR-200b-3p; these miRNAs were essential hubs in the entire network. Moreover, these miRNAs potentially play an important role in the pathogenesis of DR.

Figure 4.

The ceRNA network based on the mRNA-miRNA and lncRNA-miRNA co-expression pair network by using Cytoscape. ceRNA, competitive endogenous RNA; mRNA, messenger RNA; miRNA, microRNA; lncRNA, long non-coding RNA.

Functional evaluation of the diabetic retinopathy-specific lncRNAs

We observed that mRNAs were connected to lncRNAs in the network. To comprehend the functions of every lncRNA, the functions of connected mRNAs were assessed. Thus, we deciphered how gene products can be enriched in biological processes, cellular components, and molecular functions. Surprisingly, the results of GO analysis indicate that enriched GO terms were associated with the following events: the development of eyes and lens in camera-type eye, visual perception, structural constituents of eye lens, structural molecule activity and calcium-binding ions (Fig. 5A). We visualized gene function enrichment in the network by using BinGO (Fig. 5B).

Figure 5.

Gene Ontology (GO) function enrichment analysis of mRNAs: (A) GO terms displayed as an interaction network by using the plug-in BinGO of Cytoscape software. (B) The significant changes in the GO of lncRNA-related mRNAs. mRNAs, messenger RNAs; lncRNA, long non-coding RNA.

Key lncRNA-miRNA-mRNA sub-network based on the module analysis

To investigate the cross-talks between mRNAs and lncRNAs, key lncRNAs were extracted; moreover, the correlation between mRNAs and miRNAs was determined by WGCNA analysis. Thus, we constructed a novel sub-network of lncRNA-miRNA-mRNA. We discovered a module that was significantly related to DR (P<0.001, R=0.980). In this module, 33 lncRNAs and 10 mRNAs were involved (Fig. 6).

Figure 6.

(A) Heat map of key genes, including 33 lncRNAs and 10 mRNAs. (B) The interaction network of key genes, including 33 lncRNAs and 10 mRNAs. lncRNAs, long non-coding RNAs; mRNAs, messenger RNAs.

Discussion

In the past few decades, many studies have been conducted to determine the molecular mechanisms underlying diabetic retinopathy (DR). Protein-coding genes or miRNAs have been previously identified. Yet, very few studies have determined the function of lncRNAs which are associated with the pathogenesis of DR. Scientists have not been successful in elucidating the functions of the lncRNA-miRNA-mRNA network in DR to date. Here, we report certain lncRNAs which can be used as miRNA sponges, where these sponges are involved in the pathogenesis of DR. To the best our knowledge, this is the first study to determine whether ceRNAs play a pivotal role in the development of DR.

In the present study, we performed GO analysis to determine the possible functions of the differentially expressed lncRNAs. We found that the functions of lncRNAs were completely related to the functions of the connected mRNAs; this notion was used to determine gene product enrichment in biological processes, cellular components and molecular functions. By performing GO analysis of mRNAs, we found that differentially expressed genes were involved in the following ophthalmic functions: Eye development, lens development in camera-type eye, visual perception and structural constituents of eye lens. This indicates that disease-related lncRNAs play important roles in the pathogenesis of DR. Furthermore, we also performed KEGG analysis; however, the results indicated that disease-specific genes were not significantly associated with certain pathways.

It has been established that WGCNA is a widely-used bioinformatic method that describes correlation patterns among genes of microarray samples. We performed WGCNA analysis for the following purposes: i) to identify clusters (modules) of highly correlated genes, ii) to summarize clusters (module) by using eigen gene or an intramodular hub gene, iii) to relate modules with one another and external sample traits using eigen gene network methodology, and iv) to calculate module membership measures. Although the analysis was performed with utmost diligence, we obtained a paradoxical result: We found that only one module was specifically related to DR. This module included 33 lncRNAs and 10 mRNAs. These 33 lncRNAs might serve as potential diagnostic biomarkers and therapeutic targets for DR.

Importantly, previous studies have reported that some key hub miRNAs, depicted in this study, play an important role in the pathogenesis of different diseases. For example, Bozec et al investigated the functions of miR-223-3p, reporting that this molecule is associated with the angiogenesis of head and neck squamous cell carcinoma (27). This indicates that miR-223-3p plays a pivotal role in the development of DR. In another study, miR-223-3p was found to be involved in the pathogenesis and progression of diabetic kidney disease (28). Another study reported that miR-34c-5p is an inflammation-related miRNA, which is associated with the vascular repair factor HGF and miR-574-3p. It is upregulated in type 2 diabetes monocytes, implying that anti-inflammatory cells play a pivotal role in adhesion, vascular repair and invasion (29). Moreover, miR-34c-5p is also a tumor suppressor (30), and it is associated with apoptosis and differentiation of cancer cells. Furthermore, miR-200b-3p is also a hub molecule in ceRNA network. It was found that miR-200b-3p is associated with several types of cancers, dysregulating monocyte/macrophage (31) and differentiating epithelial-to-mesenchymal transition of glioma cancer cells (32).

Our study has several limitations. A major concern is the insufficient data and a relatively small sample size. Owing to the limited sample size, the miRNA-lncRNA-mRNA network might be restricted. Moreover, the type of diabetes was not reported by study participants. Another major concern of this research study is that computational results might be noises and false-positive results; however, we set cut-offs to the most used values. Moreover, we used the widely accepted methods of computation. For example, the target mRNAs of miRNAs were predicted on the basis of three widely used databases (miRTarBase, miRecords and starBase). These results need further confirmation in terms of cell lines and animal models.

In summary, we constructed a competitive endogenous RNA network and identified several potential key molecules. We described the important role played by non-coding RNAs in the pathogenesis of DR. Our study highlighted specific lncRNAs and miRNAs related to the pathogenesis of DR, which might be used as novel diagnostic biomarkers and therapeutic targets for DR.