Identification of the Potential Key Genes and Pathways Involved in Lens Changes of High Myopia

Introduction

High myopia (HM) is estimated to affect near 10% of the global population in 2050, a leading cause of blindness worldwide.^1,2 Various mechanisms in myopia progression was reported in previous studies.^3–6 Similar to other complex diseases, myopia is considered result from the interaction of genetic and environmental factors.⁷ Several biological pathways, including ion transport, neurotransmission, retinoic acid metabolism, eye development and extracellular matrix remodelling, partially explained mechanism of a retina to sclera signalling cascade in myopia.⁸ In addition, these previous studies on HM mainly focused on the genome (eg genome-wide association studies and single nucleotide polymorphism).^9,10 However, the underlying molecular dysregulation of HM in transcriptome remain largely unknown.

Lens is the core refracting medium of human visual system and responsible for the full range of vision. A previous study¹¹ found that patients with HM may have higher rate of lens diseases. Yet despite the lens diameter of HM eyes was considered larger and associated with the TGF-β1-Smad signaling.¹²

It is not clear whether these lens changes are causal or merely markers for HM. It is urgently needed to propose more potential mechanisms. The development of high-throughput sequencing technology and bioinformatics has greatly accelerated our understanding of human diseases.^13–15 Studied of mining and re-analysis of these big-data also saves research expenditure. Interestingly, HM-related high-throughput data and bioinformatics research are few. Herein, we used a HM-related data set to perform bioinformatics analysis to reveal the potential key genes and pathway in lens epithelia of HM.

Methods and Materials

Data Processing

The HM-related data set GSE136701¹² was downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). GSE136701 contains gene expression profiles based on the GPL570 of lens epithelia from three human highly myopic eyes and three emmetropic control eyes and was based on the GPL570 platform. In our present study, the probes corresponding to multiple genes are deleted. If a gene corresponds to multiple probes, the average value of these probes is considered as the expression value of the gene. As these data are publicly available and open-access, ethical approval from the ethics committee of First Affiliated Hospital of Guangxi Traditional Chinese Medical University was not necessary for the present study.

Gene Set Enrichment Analysis (GSEA)

GSEA^16,17 was carried out using the GSEA java software to identify the HM-specific Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The canonical pathways gene sets derived from the KEGG pathway database were download from the Molecular Signatures Database (version 7.4).^17,18 Gene set with the nominal P value < 0.05 was considered significantly enriched in HM.

Gene Differential Expression Analysis and Cluster Analysis

To identify the potential HM-specific genes, we screened the differentially expressed genes (DEGs) in lens epithelia of HM compared to emmetropic control using limma¹⁹ package. We considered genes with a |log2 (fold change)|>1 and a P value < 0.05 to be significant. Subsequently, we performed bidirectional hierarchical clustering to the optimal DEGs based on Euclidean distance and displayed the results as a heat map using the pheatmap [https://CRAN.R-project.org/package=pheatmap] package.

Protein–Protein Interaction (PPI) Network Analysis

The PPI networks of the DEGs were constructed in the STRING database (version 11.5)²⁰ and visualized using the Cytoscape software.²¹ The interactions with medium confidence (combine score > 0.4) were included in our present study. Then, we separately identified the hub genes using three methods (degree,²² Maximum Neighborhood Component (MNC),²³ and Maximal Clique Centrality (MCC)²⁴) using the cytoHubba²⁴ plug-in of Cytoscape. Additionally, we used the Cytoscape plug-in ClusterONE²⁵ to predict gene clusters based on a cohesion algorithm and nearest neighbor selection.

Functional Enrichment Analysis

In order to preliminary reveal the biological functions of each gene cluster, we performed functional enrichment analysis using the clusterProfiler²⁶ package. The P adjusted by Benjamini & Hochberg <0.05 was considered significant.

Results

High Metabolism Tendency Found in Lens Epithelia of HM

The GSEA results indicated that multiple metabolism-related pathways were significantly enriched in lens epithelia of HM, including galactose metabolism (Figure 1A), pyruvate metabolism (Figure 1B), insulin signaling pathway (Figure 1C), glycolysis gluconeogenesis (Figure 1D), neurotrophin signaling pathway (Figure 1E), and purine metabolism (Figure 1F).

Figure 1 Metabolism-related pathways enriched in lens epithelia of highly myopic. (A) Galactose metabolism pathway, (B) pyruvate metabolism pathway, (C) insulin signaling pathway, (D) glycolysis gluconeogenesis pathway, (E) neurotrophin signaling pathway, and (F) purine metabolism pathway.

Multiple Genes Differentially Expressed in HM Compared to Emmetropic

A total of 873 genes were considered as DEGs according to our cutoff criteria, including 426 up-regulated genes and 447 down-regulated genes (Figure 2A). The top 10 (ranked by FC) up-regulated genes in HM were USP9Y, CRYBB1, PALM2, HTATSF1P2, CTXN3, LGSN, TMEM130, SFTA3, CRYBA4, and CRYBA2, while the top 10 (ranked by FC) down-regulated genes in HM were LINC00314, RP1-28C20.1, CNGA4, RP11-327J17.2, CBLN1, LINC00682, LOC102724718, HIST1H2BK, BLACAT1, and LOC100506122. The results of clustering analysis indicated the expression patterns of these DEGs could accurately distinguish HM and emmetropic (Figure 2B).

Figure 2 Multiple genes differentially expressed in highly myopic compared to emmetropic. (A) The volcano of the differentially expressed genes. The top 10 up-regulated and down-regulated genes were highlighted. (B) Clustering analysis indicated the expression patterns of these differentially expressed genes accurately distinguish highly myopic and emmetropic. Abbreviations: NS, not significant; FC, fold change.

Identification of Potential Hub Genes Involving in HM from PPI Networks

A PPI networks of the DEGs, including 477 nodes and 1017 edges, were constructed using STRING database (Figure 3). This suggested that HM results from the synergy of multiple dysregulated genes. The ten hub genes identified by degree method were LCK, FGR, SHC1, EGF, CAMK2A, CD40, CD34, IL1A, CD40LG, and CXCL12 (Figure 4A). The ten genes identified by MCC method were IL33, CD34, IFNB1, CCL11, CD40, CXCR2, EGF, IL1A, CD40LG, and CXCL12 (Figure 4B). The ten genes identified by MNC method were VAV1, CD34, SHC1, CCL11, CD40, LCK, EGF, IL1A, CD40LG, and CXCL12 (Figure 4C). Thus, CD34, CD40, EGF, IL1A, CD40LG, and CXCL12, as the overlap hub genes in the three algorithms, maybe the potential key genes involved in HM.

Figure 3 The protein–protein interaction (PPI) networks of the differentially expressed genes.

Figure 4 Hub genes and the functional enrichment analysis. (A) Top 10 hub genes identified by degree method; (B) top 10 hub genes identified by maximal clique centrality method; (C) top 10 hub genes identified by maximum neighborhood component method; (D) Gene Ontology terms of the three gene clusters, (E) Kyoto Encyclopedia of Genes and Genomes pathway analysis of the three gene clusters.

Identification of Potential Key Pathways Involving in HM from PPI Networks

Three gene clusters were identified in the PPI networks of DEGs may possess according to the clusterONE analysis. The three gene clusters involved in distinct GO terms (Figure 4D) and pathways (Figure 4E). The gene cluster 2 were significantly involved in proliferation-related pathways, such as TGF−beta signaling pathway, and signaling pathways regulating pluripotency of stem cells. Genes of cluster 3 were significantly involved in immune-related pathways and NF−kappa B signaling pathway. We also found that gene clusters 1 and 2 are both involved in the MAPK signaling pathway and calcium signaling pathway, which indicated that these two pathways may play a critical role in HM. Subsequently, we separately highlighted the DEGs of clusters 1 and 2 in MAPK signaling pathway (Figure 5) and calcium signaling pathway (Figure 6).

Figure 5 Genes of clusters 1 and 2 in MAPK signaling pathway. Red represents up-regulated, and green represents down-regulated. GF includes EGF and FGF9 in cluster 2; CACN includes CACNA1B, CACNA2D1, and CACNA2D2 in cluster 1; PTP represents PTPN5 in cluster 1; RasGRP represent RASGRP1 in cluster 1, RTK represents FGFR4 in cluster 2.

Figure 6 Genes of clusters 1 and 2 in calcium signaling pathway. Red represents up-regulated, and green represents down-regulated. GF includes EGF and FGF9 in cluster 2; NCX represent SLC8A2 in cluster 1; CaV2 represents CACNA1B in cluster 1; RTK represents FGFR4 in cluster 2; RYR represents RYR3 in cluster 1; CAMK represents CAMK2A in cluster 1.

Discussion

HM is one of the blinding ocular diseases with an increasing prevalence worldwide. As the crystalline lens were the core refracting medium of the eye,²⁷ thus it may be a critical step to understand the changes of the crystalline lens to prevent HM. Previous studies suggested no difference in lens thickness^28–30 but larger in lens diameter size¹² in highly myopic eyes compared to emmetropic eyes. In our present study, various energy metabolic-related gene sets were enriched in the lens epithelia of HM eyes, which indicated more lens content is produced in HM eyes. This may partly explain the reason for the increase in the diameter of the lens. Yet further studies are required to validate this assumption.

In the current management of myopia, one of the reasons why few drugs are recommended³¹ is that the molecular mechanism of this disease is not well understood. Zhu et al reported that the aberrant TGF-β1 signaling activation by MAF underlies pathological lens growth in HM.¹² In our present study, compared to emmetropic eyes, we found multiple genes showed differentially expression patterns in lens epithelium of HM eyes. The potential key genes (CD34, CD40, EGF, IL1A, CD40LG, and CXCL12) were identified using different methods. EGF was considered associated with various ophthalmological diseases.^32–34 A previous study showed there is association between IL1A polymorphisms and primary open angle glaucoma.³⁵ CD40 and CD40LG were considered associated with pathogenesis of thyroid eye disease.³⁶ In addition, we also proposed that the imbalances of MAPK signaling pathway and calcium signaling pathway may be the two crucial steps of pathological lens growth in HM. In HM, up-regulated growth factors (EGF and FGF9) and down-regulated RASGRP1, PTPN5, CACNA1B, CACNA2D1, CACNA2D2, and FGFR4 maybe the mainly dysregulated in MAPK signaling pathway. However, this may be a compensatory result, as the energy metabolism-related pathway is enriched in HM compared to emmetropic control. A more active calcium signaling pathway was found in HM, behaving as the up-regulation of EGF, FGF9, SLC8A2, RYR3, and CAMK2A.

Though our present study may provide new insight into the pathological lens in HM, it has several limitations. Firstly, the sample size was small, however, this limitation cannot be addressed for a short time due to the lacking attention in this field. Secondly, this is bioinformatics-based study and further molecular experimental validation for these potential key genes and pathways is required. Thirdly, the key genes were identified based on transcriptomic level and PPI networks, whether these genes mutate or exhibit abnormal epigenetic changes is unknown.

In conclusion, the present study identified potential key genes in pathological lens growth of HM eyes and proposed that the imbalances of MAPK signaling pathway and calcium signaling pathway may be the two crucial steps of pathological lens growth in HM.