Bioinformatics Analysis of Weighted Genes in Diabetic Retinopathy

Abstract

Purpose: Intricate signaling networks and transcriptional regulators translate pathogen recognition into defense responses. The aim of this study was to identify the weighted genes involved in diabetic retinopathy (DR) in different rodent models of diabetes.

Methods: We performed a gene coexpression analysis of publicly available microarray data, namely, the GSE19122 dataset from the Gene Expression Omnibus database. We conducted gene coexpression analysis on the microarray data to identify modules of functionally related coexpressed genes that are differentially expressed in different rodent models. We leveraged a richly curated expression dataset and used weighted gene coexpression network analysis to construct an undirected network. We screened 30 genes in the most closely related module. A protein-protein interaction network was constructed for the genes in the most related module using the Search Tool for the Retrieval of Interacting Genes. Gene Ontology enrichment analysis and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis were performed for the 30 genes.

Results: Five visual perception-related genes (Pde6g, Guca1a, Rho, Sag, and Prph2) were significantly upregulated. Based on the competing endogenous RNA hypothesis, a link between the long noncoding RNA metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) and visual perception-related mRNAs was constructed using bioinformatics tools. Six potential microRNAs (miR-155-5p, miR-1a-3p, miR-122-5p, miR-223-3p, miR-125b-5p, and miR-124-3p) were also screened.

Conclusions: MALAT1 might play important roles in DR by regulating Sag and Guca1a through miR-124-3p and regulating Pde6g through miR-125b-5p.

Diabetic retinopathy (DR) is the most common complication of diabetes mellitus and the major cause of vision loss globally. The treatment of DR remains challenging and most patients fail to obtain clinically significant visual function improvement.¹ To identify better therapies for DR, a thorough investigation of the molecular mechanisms of DR is critical. 

Only 1% to 2% of the genome encodes proteins, and the majority of the mammalian genome encodes a large amount of noncoding RNAs (ncRNAs), which regulate gene expression at the transcriptional and posttranscriptional levels.² MicroRNAs (miRNAs) are small ncRNA molecules consisting of 18 to 25 nucleotides that regulate a variety of cellular processes by binding to a specific target mRNA with a complementary sequence to induce its cleavage or degradation and thereby suppress the expression of target genes.³ Long ncRNAs (lncRNAs), a type of ncRNAs containing more than 200 nucleotides, play vital roles in modifying chromatin states and influencing gene expression.⁴ Circular RNAs (circRNAs), a novel class of ncRNAs, form a covalently closed loop, and their functions include acting as scaffolds in the assembly of protein complexes, modulating the expression of parental genes and RNA-protein interactions, and functioning as miRNA sponges.⁵ Emerging data suggest that ncRNAs are expressed and play critical roles in the pathogenesis of DR. Metastasis-associated lung adenocarcinoma transcript 1 (MALAT1), a highly conserved lncRNA, is abnormally expressed in DR.⁶ However, the involvement of this lncRNA in the complex molecular mechanisms of DR remains largely unknown. 

Recently, accumulating evidence has demonstrated that miRNAs play vital roles in DR by reducing the expression of their targets, which include mRNAs, lncRNAs, circRNAs, and pseudogenes.^7–9 Nevertheless, target overexpression can abolish the downregulatory effects of these miRNAs.^10–12 Moreover, multiple miRNA targets can function as competing endogenous RNAs (ceRNAs) and compete with each other to bind the miRNA; the overexpression of one ceRNA can also upregulate other ceRNAs.^13,14 This ceRNA crosstalk was first suggested by Poliseno et al.,¹⁵ who demonstrated that the tumor suppressor gene phosphatase and tensin homolog deleted on chromosome ten (PTEN) could be upregulated by phosphatase and tensin homolog pseudogene 1, which is a PTEN pseudogene. Furthermore, many subsequent studies have examined a plethora of other ceRNA phenomena correlated with DR, including the crosstalk between hypoxia-inducible factor 1-alpha and vascular endothelial growth factor (VEGF),¹⁶ as well as that between VEGF and angiogenin 2.¹⁷ 

Therefore, ceRNA crosstalk is a crucial mechanism underlying the complex pathogenesis and multistep development of DR; it might represent a potential target for developing new therapies and should be further studied. Herein, we revealed the ceRNA network between ncRNAs and coding genes and identified several pivotal biological pathways closely associated with the onset and progression of DR, which could help build a promising new therapeutic system for this disease. 

Methods

Material and Data

The datasets used in the present study were downloaded from the National Center of Biotechnology Information Gene Expression Omnibus. The original gene expression profile was obtained from the GSE19122 dataset,¹⁸ which included 31 mouse retina samples from two mouse models of diabetes. Retinal transcriptomic responses of the both streptozotocin-induced and Ins2^Akita diabetic mice were compared after 3 months of hyperglycemia. The platform used for these data was the GPL6885 Illumina MouseRef-8 version 2.0 expression BeadChip (Illumina, San Diego, CA, USA). Background correction and quartile data normalization of the downloaded data were performed using the robust multiarray average algorithm.¹⁹ Probes without a corresponding gene symbol were filtered, and the average values of gene symbols with multiple probes were calculated. 

Weighted Gene Coexpression Network Analysis

Sample clustering was performed to demonstrate the relationship between expression profiles and clinical traits. After raw data preprocessing, weighted gene coexpression network analysis (WGCNA) was performed based on a previously described algorithm to identify significant gene modules.²⁰ Probe sets were first filtered based on the variance of expression values across all samples. Probe sets with duplicated gene symbols were removed based on expression variance. The R package WGCNA²¹ was used for this analysis. Briefly, Pearson's correlation coefficients were calculated for the selected genes in a pairwise manner, yielding a similarity matrix (Sij). The soft threshold (power) was set to 12. The matrix was transformed into an adjacency matrix (aij) by using a power function and the following formula: aij = Power (Sij, β) ≡ |Sij| β. Average linkage hierarchical clustering was then performed to identify modules of densely interconnected genes. Network interconnectedness was measured by calculating topologic overlap by using the Topological overlap matrix dissimilarity function with a signed TOM-Type. Average hierarchical clustering was performed using the hclust function to group the genes based on the topological overlap dissimilarity measure (1-TOM) of their connection strengths. Network modules were identified using a dynamic tree cut algorithm with minimum cluster size of 30 and a merging threshold function of 0.25. Genes that were not assigned to specific modules were assigned the color gray. 

Protein-Protein Interaction (PPI) Network Construction and Analysis

The PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes online database. PPI pairs with a combined score of >0.4 were used to construct the PPI network. Then, the regulatory relationship between genes was visualized using Cytoscape (version 3.4.0) and analyzed based on the topologic properties of the computing network, including the degree of the distribution of the network, using the CentiScaPe app.²² 

Gene Function Analysis

Gene Ontology (GO) enrichment analysis of mRNAs was implemented via DAVID (https://david.ncifcrf.gov/, in the public domain). Briefly, gene identifiers were first converted into their H. sapiens Entrez gene identifiers (IDs) by using the most recent database. If multiple identifiers corresponded to the same Entrez gene ID, they were considered a single Entrez gene ID in downstream analyses. For each given gene list, pathway and process enrichment analyses were performed using the following ontology sources: Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, GO biological process, GO molecular function, and GO cellular component. All genes in the genome were used as the enrichment background. More specifically, P values were calculated based on accumulative hypergeometric distribution, and q-values were calculated using the Benjamini-Hochberg procedure to account for multiple testing. Kappa scores were used as the similarity metric when performing hierarchical clustering on the enriched terms, and subtrees with a similarity of >0.3 were considered a cluster. The most statistically significant term within a cluster was chosen as the term representing the cluster. 

Prediction of miRNA/mRNA/lncRNA Interactions

The interaction between MALAT1 and miRNAs was predicted using LncBase version 2.²³ The interaction between miRNAs and mRNAs was predicted using TargetScan,²⁴ and only conserved miRNAs were included in this study. Cytoscape (version 3.40) was used to visualize the miRNA regulatory network. 

Results

Construction of Coexpression Network and PPI Network

After filtering out the probe sets with no significant variance in expression across all samples (Fig. 1A), we used the R package WGCNA to generate 9 modules from 6424 probe sets (Figs. 1B–C). All uncorrelated genes were assigned a gray module. The trait in this study is the different mouse models. Among the nine modules, the magenta module was significantly associated with the Ins2^Akita diabetic mouse model (r = 0.66, P = 6E−5) (Fig. 1D). Therefore, we chose the 30 genes in the magenta module to construct the PPI network and coexpression network. Ultimately, a PPI network containing 15 nodes and 27 edges (Fig. 2A; Supplementary Table S1) and a coexpression network containing 27 nodes and 139 edges (Fig. 2B; Supplementary Table S2) were obtained. 

Figure 1

 

(A) Sample clustering to detect outliers. All samples are represented in the clusters, and all samples passed quality control. (B) Analysis of network topology for various soft-thresholding powers. The left panel presents the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). The right panel displays the mean connectivity (degree, y-axis) as a function of the soft-thresholding power (x-axis). (C) Hierarchical clustering dendrogram of the 6424 probe sets used in the analysis. Each line represents an individual probe. Genes were clustered based on a dissimilarity measure (1-TOM). The branches correspond to modules of highly interconnected groups of genes. The tips of the branches represent genes that are the least dissimilar and thus share the most similar network connections. Genes within gray-shaded boxes were not assigned a module. (D) Correlation matrix of module eigengene values obtained for different mouse models of diabetes. WGCNA groups mRNAs into modules based on the patterns of their coexpression. Each of the modules was labeled with a unique color as an identifier. Nine modules were identified; the eigengene value of each module was tested for correlation with the trait. Within each cell, the upper values represent correlation coefficients between the module eigengene and the trait, and the lower values are the corresponding P values. The scale bar represents the correlation pattern (1 = strong positive correlation, 0 = no correlation, −1 = strong negative correlation).

Figure 2

 

(A) PPI network for the genes grouped in the magenta module. Nodes represent mRNAs, and edges represent interactions between proteins. The nodes with red borders represent mRNAs enriched for visual perception terms, and the nodes filled with orange represent mRNAs enriched for phototransduction. (B) Coding-non-coding gene co-expression (CNC) network for the genes grouped in the magenta module. Nodes represent genes, and edges represent coexpression between genes. The nodes with a blue border represent mRNAs enriched for visual perception terms.

GO and KEGG Enrichment Analysis of the Selected mRNAs

GO and KEGG enrichment analyses were performed for the 30 selected mRNAs to investigate the biological function of these genes, and the results are presented in the Table. In the GO and KEGG enrichment analyses, all the results were ranked using an enrichment score (−log [P value]). In the KEGG enrichment analysis, phototransduction was enriched. In the biological process analysis, visual perception was the most enriched term. Moreover, visual perception-related genes (guanylate cyclase activator 1A [Guca1a], phosphodiesterase 6G [Pde6g], peripherin 2 [Prph2], rhodopsin [Rho], and S-antigen visual arrestin [Sag]) were upregulated in the Ins2^Akita diabetic mouse model (Fig. 3). 

Table

 

GO and KEGG Enrichment Analysis of Select mRNAs in the Magenta Module

 

Figure 3

 

Expression of visual perception-related genes in different mouse models of diabetes. (A) Guca1. (B) Pde6g. (C) Prph2. (D) Rho. (E) Sag. (F) Heatmap of five genes related to visual perception.

Construction of a MALAT1-miRNA-mRNA Regulatory Network

According to a previous report, MALAT1 is an epigenetic regulator of inflammation and is upregulated in DR.⁶ To investigate the relationship between MALAT1 and visual perception, we constructed a ceRNA network between MALAT1 and visual perception-related genes using bioinformatics tools. Ultimately, we obtained six potential miRNAs (miR-155-5p, miR-1a-3p, miR-122-5p, miR-223-3p, miR-125b-5p, and miR-124-3p) that provide a link between MALAT1 and visual perception-related genes (Fig. 4; Supplementary Table S3). 

Figure 4

 

MALAT1-miRNA-mRNA regulatory network. Circles represent the mRNAs, arrowheads represent miRNAs, and the octagon represents the lncRNA MALAT1.

Discussion

Many lncRNAs have important roles in diverse biological processes. Previous studies have implicated the upregulation of the lncRNA MALAT1 in the pathogenesis of diabetes-related microvascular disease and DR.²⁵ MALAT1 upregulation represents a critical pathogenic mechanism for diabetes-induced microvascular dysfunction. MALAT1 knockdown could clearly ameliorate DR in vivo, as demonstrated by pericyte loss, capillary degeneration, microvascular leakage, and retinal inflammation.⁶ 

MALAT1 inhibition may serve as a potential target for antiangiogenic therapy for diabetes-related microvascular complications.²⁵ However, the role of MALAT1 in visual perception remains elusive. Here, we demonstrate that five visual perception-related mRNAs (Pde6g, Guca1a, Rho, Sag, and Prph2) are significantly upregulated in diabetic mice, and four of these mRNAs are regulated by MALAT1 through six miRNAs. 

Given that the MALAT1-miRNA-mRNA regulatory network was constructed based on the ceRNA hypothesis and that both MALAT and mRNAs are upregulated in DR, the expression of the six miRNAs should be downregulated. Huang et al.²⁶ demonstrated that miR-155 expression is more than five-fold greater in kidney samples from diabetic nephropathy patients than in those from controls and that miR-155 expression gradually increases during disease induction and progression in rat models of type 1 and type 2 diabetic nephropathy. Loscher et al.²⁷ and Anasagasti et al.²⁸ reported that miR-1a was upregulated in mouse models. In addition, miR-122 is a highly abundant, hepatocyte-specific miRNA, and miR-122 affects lipid metabolism in the liver.^29,30 

Another study assessing cardiomyocyte glucose metabolism found that miR-223 was upregulated in the left ventricle of type 2 diabetes mellitus patients. Moreover, miR-223 overexpression increased glucose transporter 4 (GLUT4) protein levels in cardiomyocytes, and miR-223 inhibition in vivo significantly reduced GLUT4 expression.³¹ Indeed, based on previous studies, these four miRNAs might not represent a link between MALAT1 and the mRNAs. 

Gong et al.³² demonstrated that miR-124/miR-125b expression decreased with DR progression in vitro and in vivo; miR-124-3p and miR-125b-5p expression gradually decreased as the exposure time to high glucose increased, and these effects occurred concomitantly with the early development of DR. The tumor suppressor miR-124 is downregulated in cancers and plays a role in regulating cell proliferation and apoptosis by targeting different genes in distinct types of cancer.^33,34 Gong et al. also reported that the downregulation of miR-125b correlates with the upregulation of SP1, a key transcription factor involved in various types of diseases that plays vital roles in angiogenesis, inflammation, cell proliferation, migration/invasion, and survival.^35,36 

In summary, based on these previous studies and our results, we hypothesized that MALAT1 might play important roles in DR by regulating Sag and Guca1a through miR-124-3p and regulating Pde6g through miR-125b-5p. 

Acknowledgments

Supported by National Natural Science Foundation of China (Grant numbers 81760176 and 81860175), Jiangxi Provincial Training Program for Distinguished Young Scholars (Grant number 20171BCB23092), Jiangxi Provincial Key R&D Program (Grant number 20171BBG70099), Jiangxi Provincial Natural Science Foundation for Youth Scientific Research (Grant number 20171BAB215032), Scientific Research Foundation of Jiangxi Education Department (Grant number GJJ170095), and Youth Scientific Research Foundation of the Second Affiliated Hospital of Nanchang University (Grant number 2014YNQN12011). 

Disclosure: Z.-P. You, None; Y.-L. Zhang, None; B.-Y. Li, None; X.-G. Zhu, None; K. Shi, None