Exploring the Molecular Mechanisms of Pterygium by Constructing lncRNA–miRNA–mRNA Regulatory Network

Abstract

Purpose: This research explores the aberrant expression of the long non-coding RNA (lncRNA), microRNA (miRNA), and messenger RNA (mRNA) in pterygium. A competitive endogenous RNA (ceRNA) network was constructed to elucidate the molecular mechanisms in pterygium.

Methods: We obtained the differentially expressed mRNAs based on three datasets (GSE2513, GSE51995, and GSE83627), and summarized the differentially expressed miRNAs (DEmiRs) and differentially expressed lncRNAs (DELs) data by published literature. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, protein-protein interaction (PPI), and gene set enrichment analysis (GSEA) analysis were performed. DEmiRs were verified in GSE21346, and the regulatory network of hub mRNAs, DELs, and DEmiRs were constructed.

Results: Overall, 40 upregulated and 40 downregulated differentially expressed genes (DEGs) were obtained. The KEGG enrichment showed the DEGs mainly involved in extracellular matrix (ECM)-receptor interaction, focal adhesion, and PI3K-Akt signaling pathway. The GSEA results showed that cornification, keratinization, and cornified envelope were significantly enriched. The validation outcome confirmed six upregulated DEmiRs (miR-766-3p, miR-184, miR-143-3p, miR-138-5p, miR-518b, and miR-1236-3p) and two downregulated DEmiRs (miR-200b-3p and miR-200a-3p). Then, a ceRNA regulatory network was constructed with 22 upregulated and 15 downregulated DEmiRs, 4 downregulated DELs, and 26 upregulated and 33 downregulated DEGs. The network showed that lncRNA SNHG1/miR-766-3p/FOS and some miRNA-mRNA axes were dysregulated in pterygium.

Conclusions: Our study provides a novel perspective on the regulatory mechanism of pterygium, and lncRNA SNHG1/miR-766-3p/FOS may contribute to pterygium development.

Pterygium is a prevalent ocular surface disease that occurs most frequently in tropical equatorial areas.^1,2 It comprises a wing-shaped progressively growing fibrovascular tissue usually located on the nasal side, and could lead to visual impairment, astigmatism, and dry eye.³ Over the years, various medical measures, such as anti-inflammatory eye drops and chemotherapeutic agents have been used in the treatment of pterygium. Surgical removal can be performed if the patient desires symptomatic or cosmetic improvement, but recurrence remains the main complication. The recurrence rate of different surgical techniques ranges from 0% to 88%.^4–6 Hence, elucidating the pathogenesis and molecular mechanisms of pterygium is crucial for improving surgical outcomes and decreasing the risk of recurrence. 

Previous clinical and laboratorial studies showed that immunologic mechanisms, extracellular matrix (ECM) modulation under ultraviolet radiation exposure, cell proliferation and hyperplasia, inflammation, angiogenesis, cholesterol metabolism modification, and hereditary factors attributed to the pathogenesis of pterygium.^7,8 Moreover, in recent years, the discovery of non-coding RNAs (regulatory RNAs) made the molecular etiology of pterygium more complex. 

MicroRNA (miRNA) is a group of single-stranded noncoding RNA (ncRNA) that downregulates gene expression at the post-transcriptional level by inhibiting translation or promoting degradation of target messenger RNA (mRNAs). Long non-coding RNA (lncRNA) also belongs as a member of the noncoding RNA family, which is longer than 200 nt in length and has little or no protein-coding ability. MiRNAs and lncRNAs are involved in many physiological and pathophysiological conditions. Recent studies have demonstrated that lncRNAs can work as competitive endogenous RNA (ceRNA) with miRNAs to compete with mRNAs for binding with miRNAs, thus affecting gene expression.⁹ However, the research of core RNAs and lncRNA-miRNA-mRNA regulatory networks in pterygium using bioinformatics analysis was still devoid. 

In the present study, we utilized mRNA microarray dataset from Gene Expression Omnibus (GEO) to obtain and analyze differentially expressed genes (DEGs) and pooled the expression profiling of differentially expressed miRNAs (DEmiRs) and lncRNAs (DELs) between pterygium and normal conjunctiva tissues. Afterward, gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, protein-protein interaction (PPI) network, and gene set enrichment analysis (GSEA) were performed. Furthermore, lncRNA-miRNA-mRNA regulatory networks were explored using StarBase and DIANA-LncBase, which helped us understand the pathogenesis and molecular mechanisms of pterygium. 

Methods

Subjects and Gene Information

The GEO is a national center for genetic information database, including a large number of gene chips, methylation, and sequencing data.¹⁰ The including criteria: 1. Samples were obtained from Homo sapiens. 2. The chip data included both pterygium and normal conjunctiva tissues. 3. Chip data belonged to different independent studies, and all chip data did not contain each other. In this study, three gene expression profiles (GSE2513, GSE51995, and GSE83627) were searched and selected from the GEO database. In brief, GSE2513 consisted of 4 conjunctiva samples and 8 pterygium samples, which were harvested from 7 Chinese and 7 non-Chinese with age distribution of 42 to 57 years. In GSE51995, four primary nasal pterygium and four uninvolved conjunctiva tissues were collected from the superior temporal quadrant of the same eye. GSE83627 contained four donor-matched pterygium and conjunctiva tissues without mentioning the other clinical information. 

Data Analysis and DEGs Screening

We used the GEO mirror of R packages to get the expression matrix of GSE2513, GSE51995, and GSE83627 from the GEO Dataset. Then the expression matrixes were normalized and differential genes were screened by limma package (http://www.bioconductor.org/packages/release/bioc/html/limma.html)¹¹ between pterygium and conjunctiva samples. The |log₂ Fold change| > 1.5 and adjust P value < 0.05 were used as the selection criteria. 

Enrichment Analysis of Gene Functions

GO¹² enrichment analysis and KEGG¹³ enrichment analysis are the two most widely used analysis strategies for gene functions. The basic unit of GO is “term,” which can be used for identifying cellular component (CC), molecular function (MF), and biological process (BP). The analysis of KEGG enrichment can show the main enrichment pathways of DEGs. In order to identify the GO annotations and pathways in which relevant DEGs were enriched, GO term and KEGG pathway enrichment analyses were performed with the Org.Hs.eg.db packages (http://www.bioconductor.org/packages/release/data/annotation/html/org.Hs.eg.db.html). The adjusted P value < 0.05 was served to distinguish significant enriched genes. 

PPI Network Visualization

PPI network of DEGs was simulated to screen out hub proteins, which played a key role in the progress of pterygium. The Search Tool for the Retrieval of Interacting Genes (STRING) database¹⁴ and Cytoscape¹⁵ were utilized for the visualization of PPI network. The STRING database covers a large number of information about proteins, including the results of protein interaction, three-dimensional structure, and related functional enrichment. Through this database, the network structure of multiple DEGs can be constructed. 

Analysis of GSEA

GSEA¹⁶ analysis was conducted in order to avoid missing the genes that actually play a crucial part during the process of screening out DEGs among three sets of data. The normal conjunctiva tissue was class A, and pterygium tissue was class B. Gene set permutations were performed 1000 times for each analysis. Absolute value of normalized enrichment score (NES) > 1 and nominal P value < 0.05 were considered as the threshold for statistical significance. 

Construction of lncRNA-miRNA-mRNA Regulatory Network

We identified miRNA through miRBase (http://www.mirbase.org),¹⁷ and constructed miRNA-mRNA regulatory network by screening Targetscan (http://www.targetscan.org/vert_72/)¹⁸ and miRDB (http://mirdb.org/).¹⁹ Both Targetscan and miRDB provided a wide-range of information on the interaction between miRNA and mRNA. Then, predicted lncRNAs interacted with miRNAs were constructed in downloaded databases StarBase version 2.0²⁰ and DIANA-LncBase version 2.0,²¹ both of which provided the experimentally validated lncRNA - miRNA interaction effect. DELs and DEmiRs were acquired by searching PubMed database for recent studies on pterygium. 

Validation of Key DEmiRs

We made a systemic search on PubMed and summarized the DEmiRs that met the criteria of P value < 0.05 and |log₂FC| > 1. The validation procedure proceeded in GSE21346, which consisted of three pterygium samples and three matched conjunctiva samples from patients diagnosed with primary pterygium. The dataset was based on GPL7723, from which we could match and compare the DEmiRs with those previously reported. The DEmiRs were considered as significant difference with P < 0.05. 

Results

Screening of DEGs

A total of 422, 1374, and 420 DEGs were identified from the GSE2513, GSE51995, and GSE83627 datasets, respectively. The specific filtering results of each dataset are shown in Table 1. A total of 40 upregulated genes were found in all 3 datasets and 40 downregulated genes were found in 2 datasets (Figs. 1A, 1B). The heatmap of 80 DEGs in GSE51995 was shown in Figure 2. 

Table 1.

 

Specific Filtering Results of Each Data Set

 

Figure 1.

 

Venn diagram DEGs were selected with a fold change > 1.5 and adjust P-value < 0.05 among the mRNA expression profiling datasets. (A) 40 up-regulated genes shared among GSE83627, GSE51995 and GSE2513. (B) Forty downregulated genes shared among GSE51995 and GSE2513.

Figure 2.

 

The heatmap of 80 DEGs in GSE2513, GSE51995, and GSE83627.

GO and KEGG Pathway Enrichment Analysis

GO analysis of individual DEGs and KEGG pathway enrichment analysis were performed by R software to obtain more insightful details into the diverse functions of particular DEGs. The main MF with significant enrichment involved with all DEGs were ECM structural constituent, peptidase regulator activity, serine-type endopeptidase inhibitor activity, sulfur compound binding, heparin binding, and endopeptidase regulator activity (Supplementary Table S1), and the main BP were response to steroid hormone, skin development, epidermis development, keratinocyte differentiation, and epidermal cell differentiation, etc. (Supplementary Table S2). The results of GO enrichment analysis of upregulated and downregulated DEGs were shown in Figures 3A, 3B. The upregulated genes were mainly related to ECM-receptor interaction, focal adhesion, PI3K-Akt signaling pathway, regulation of actin cytoskeleton, primary bile acid biosynthesis, and one carbon pool by folate, and the downregulated genes were mainly related to osteoclast differentiation, glycerolipid metabolism, and IL-17 signaling pathway (Table 2, Supplementary Tables S3A, S3B). 

Figure 3.

 

Gene ontology analysis of the DEGs. (A) GO enrichment of downregulated DEGs. (B) GO enrichment analysis of upregulated DEGs. BP, biological process; CC, cellular component; MF, molecular function. The X axis represents the number of DEGs involved in GO terms.

Table 2.

 

Mainly Pathways Involved With all DEGs

 

PPI Network Visualization

The 80 DEGs were introduced into the STRING online tool, and the proteins that did not interact with any other protein were removed to obtain the final protein interaction diagram (Fig. 4). A total of 58 nodes and 128 edges were selected to plot the PPI network, which consisted of 31 upregulated genes and 27 downregulated genes. Subsequently, 20 hub DEGs were screened out with degree ≥ 5, which might play a meaningful role in the development of pterygium (Table 3). With the help of MCODE plug-in, the top three sub-networks with most importance were analyzed (Figs. 5A–5C). 

Figure 4.

 

Protein-protein interaction network of DEGs. Note: A circle represents a protein, and an edge represents a degree. The higher the degree is, the more important the protein is in the network structure.

Table 3.

 

Top 20 Hub Genes With Degree ≥ 5

 

Figure 5.

 

Top three sub-networks in PPI.

GSEA Enrichment of all mRNAs in GSE2513

GSEA was applied to analyze the main gene sets enriched on pterygium. The results showed that cornification, keratinization, cornified envelope, peptide cross linking, extracellular structure organization, collagen fibril organization, ECM structural constituent conferring tensile strength, complex of collagen trimers, ECM component, and vascular smooth muscle contraction were top 10 gene sets with the largest NES, and the details were reported in Figure 6 and Supplementary Table S4. 

Figure 6.

 

Enrichment plots by GSEA.

DEmiRs and DELs in Pterygium

Through searching the results of 9 studies^22–30, there were 49 DEmiRs that met the criteria of P value < 0.05 and |log₂FC| > 1, all of these DEmiRs’ information was shown in Table 4. A total of 26 DEmiRs were upregulated and 23 DEmiRs were downregulated. From the research of Liu, 20 DELs were acquired (Table 5), among which 10 DELs were upregulated and 10 DELs were downregulated. 

Table 4.

 

DEmiRs of 9 Studies

 

Table 5.

 

DELs of 1 Research

 

lncRNA-miRNA-mRNA Regulatory Network

In the light of the DEmiR-DEL and DEmiR-DEG interactive pairs, the pterygium related lncRNA-miRNA-mRNA network was established in Figure 7, including 22 upregulated and 15 downregulated DEmiRs, 4 downregulated DELs, and 26 upregulated and 33 downregulated DEGs. 

Figure 7.

 

LncRNA-miRNA-mRNA regulatory network. Note: Circles represent DEGs, triangles represent DEmiRs, and chevrons represent DELs. Red means upregulated and blue means downregulated.

Validation of Key DEmiRs

In total, 6 upregulated DEmiRs (miR-766-3p, miR-184, miR-143-3p, miR-138-5p, miR-518b, and miR-1236-3p) and 2 downregulated DEmiRs (miR-200b-3p and miR-200a-3p) were verified in GSE21346, with a significant difference (P < 0.05; Fig. 8). The expression trends of these eight DEmiRs are consistent with our ceRNA network results. 

Figure 8.

 

Validation of key DEmiRs in GSE21346. Note: Green represents the control group, and purple represents the pterygium group. * P < 0.05.

Discussion

Pterygium is one of the most common ocular surface diseases. Epidemiological observations have suggested that environmental alteration, like ultraviolet radiation, is the most important factor contributing to this disease.³² However, its exact pathogenesis remains unknown. The development of high throughput microarray technology with high efficiency and high accuracy combined with bioinformatic algorithm allows us to identify key genes, which could provide a deeper understanding of molecular mechanisms on pterygium. 

In this study, we identified a total of 80 DEGs (40 upregulated and 40 downregulated mRNAs), and hub-genes with 3 subnetworks with most importance in PPI, including FN1, PI3, ERG1, SPRR1B, FOS, and FOSB, suggesting they may play a very important role in the pathogenesis of pterygium. SPRR1B belongs to keratinocyte protein, and their mRNA transcripts in conjunctival tissues increased in response to desiccating stress, which is associated with pterygium recurrence.³³ FOS and FOSB encode leucine zipper proteins that can dimerize with proteins of the JUN family, thereby forming the transcription factor complex AP-1. The encoded proteins have been implicated as regulators of cell proliferation, and experimental reports describe a tumor-suppressive function in various tumors.^34,35 EGR1 was found to suppress cell survival, proliferation, and activates expression of p53 and TGF-beta,³⁶ which was proved to play a pivotal role in the occurrence of pterygium.^8,37,38 FN1 and COL1A1 were key molecules in EMT, and were demonstrated to be involved in the pathogenic mechanism of pterygium.²⁸ Further studies were required to elucidate the complex interaction with these genes and clinical features. 

In this study, the results of functional enrichment analysis indicated that ECM-receptor interaction, focal adhesion, apoptosis, and PI3K-Akt signaling pathway were involved with significant DEGs. These results were in agreement with previous proteomics study by Hou,³⁹ in which they compared the protein expression from the conditioned medium of paired pterygium and normal conjunctival fibroblast cells from the same patients by iTRAQ-based proteomics strategy, and they found the differenced protein might serve as extracellular ligands to activate intracellular pathways. Aberrant ECM expressions were found to be a major characteristic feature of pterygium. ECM and its receptors, including fibronectin, versican, collagen III, and SPARC, were shown to be upregulated and remodeled in pterygium,^40–42 and they have turned out to be regulated by matrix metalloproteinase (MMP), an essential enzyme in local proteolysis of the extracellular matrix.⁴³ UV irradiation may lead to the imbalance of MMPs and its inhibitors TIMP, which enable the pterygial cells to dissolve corneal epithelium and Bowman's layer and invade the corneal stroma.^44–46 Besides, impression cytology specimens of pterygium show fewer apoptotic markers, and PI3K-Akt inhibitor has been proved to reduce the TGF-β-induced synthesis of FN in human pterygium fibroblasts.⁴⁷ Moreover, our GSEA-based GO analysis confirmed keratinization and cornification pathway participated in the pathogenesis of pterygium. Previous reports have confirmed that the expression of keratin markers is upregulated in the epithelium of pterygium, and takes effect in the epithelial abnormal differentiation.^48–50 Liu constructed pterygium-related mRNA libraries by using microarray and found upregulated mRNAs were closely related to proliferation and differentiation.⁵¹ If our results were combined with publications mentioned above, we could surmise the underlying mechanisms of pterygium that dysregulation of ECM proteins caused by increased expression of MMP might activate intracellular PI3K-Akt signaling pathway, leading to anti-apoptosis of fibroblast, abnormal matrix protein deposition, and hyperproliferation and keratinization of pterygial cells. 

It has been discovered that a group of miRNAs, including miR-200,²⁸ miR-145,^24,52 miR-122,³⁰ and miR-21,²⁷ were shown to contribute to the development of pterygium. In the current study, by summarizing 10 relevant studies in recent years, we found a total of 26 upregulated and 23 downregulated DEmiRs. Moreover, by establishing the miRNA-mRNA network, we found miR-215-3p significantly downregulated and directly regulated seven DEGs, including the hub gene FN. A previous research has shown that miR-215-3p takes part in regulating cell cycling and inhibiting proliferation of primary fibroblast cells from the ocular surface, and was downregulated in pterygium.²⁵ Because FN has been reported to be a key gene of epithelial mesenchymal transition (EMT) and was also identified as a target of miR-200b,²⁸ we speculated that downregulated miR-215-3p, along with miR-200b, might play important roles in pterygium by regulating EMT. 

Previous studies have identified the roles of various lncRNAs in regulating pterygium proliferation and apoptosis,³¹ but no study was reported about lncRNA-associated ceRNA network based on bioinformatics analysis in pterygium. Combining the results of DEmiR-DEL, DEmiR-DEG interactive network, and hub genes produced by the STRING, we found that the expression of lncRNA SNHG1, FOS, and FOSB were downregulated, whereas the expression of has-miR-766-3p was upregulated, which was consistent with the mechanism ceRNA hypothesis. lncRNA SNHG1 is a marker of tumor progression,⁵³ and it has been proved to negatively regulate p53, which is a tumor suppressor protein that could induce apoptosis and show strong expression in pterygium.³⁸ Downregulation of lncRNA SNHG1 could explain p53 reactivation in pterygium. MiR-766-3p has also been reported to act as a key gene in different tumors and immune diseases, and its expression was elevated in hepatocellular carcinoma, inflamed pulp, and acute promyelocytic leukemia.^53–55 This miRNA contributes to anti-inflammatory responses, cell proliferation, and apoptosis by targeting different downstream genes and signaling pathways, but the precise post-transcriptional mechanisms of miR-766-3p remain to be explored in pterygium. Further work will be needed with multiple clinical samples to clarify the ceRNA hypothesis that lncRNA SNHG1 regulates FOS and FOSB to act on the progression of pterygium through sponging miR-766-3p. 

There were also limitations in our study. The sample size of gene expression profile was not large, and no microarray data of dysregulated lncRNAs in pterygium was found in GEO and Array Express databases. In addition, the platform of GSE21346 was published 10 years ago with limited probe annotation for miRNA profile. All these may prevent us from finding a comprehensive noncoding RNA expression profile and related regulatory network. Although we proved that the expression of matched miRNAs on GSE21346 were consistent with previous studies, the mechanism and validation of these ncRNAs in pterygium still need further research in clinical and molecular biology experiments. 

To sum up, we identified core genes, and related crucial pathways, especially PI3K-Akt, keratinization, and cornification pathway involved in the pathogenesis of pterygium. Furthermore, we summarized all the studies on the miRNA and lncRNA related to pterygium from PubMed database and established lncRNA-miRNA-mRNA regulatory network. This study might deepen the understanding of potential molecular mechanism underlying pterygium and provide some new insights for use in further identification and development of new therapeutic targets for pterygium. 

Acknowledgments

Supported by the Fujian Provincial Natural Science Foundation (Grant numbers 2018J01238; 2018J01310; and 2019J01152), the Joint Funds for the Innovation of Science and Technology of Fujian province (Grant number 2018Y9035), the Startup Fund for Scientific Research of Fujian Medical University (Grant number 2017XQ1022), and the Youth Foundation of Fujian Provincial Health Commission (Grant number 2019-1-31). 

Disclosure: N. Xu, None; Y. Cui, None; J. Dong, None; L. Huang, None