Expression profiling suggests the involvement of hormone-related, metabolic, and Wnt signaling pathways in pterygium progression

Background Pterygium is an ocular surface disease that can cause visual impairment if it progressively invades the cornea. Although many pieces of research showed ultraviolet radiation is a trigger of pterygium pathological progress, the underlying mechanism in pterygium remains indistinct. Methods In this study, we used microarray to evaluate the changes of transcripts between primary pterygium and adjacent normal conjunctiva samples in China. Then, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses. Moreover, we constructed protein-protein interaction (PPI) and miRNA-mRNA regulatory networks to predict possible regulatory relationships. We next performed gene set enrichment analysis (GSEA) to explore the similarities and differences of transcripts between Asian studies from the Gene Expression Omnibus database. Furthermore, we took the intersection of differentially expressed genes (DEGs) with other data and identified hub genes of the development of pterygium. Finally, we utilized real-time quantitative PCR to verify the expression levels of candidate genes. Results A total of 49 DEGs were identified. The enrichment analyses of DEGs showed that pathways such as the Wnt-signaling pathway and metabolism-related pathways were upregulated, while pathways such as hormone-related and transcription factor-associated pathways were downregulated. The PPI and miRNA-mRNA regulatory networks provide ideas for future research directions. The GSEA of selecting Asian data revealed that epithelial-mesenchymal transition and myogenesis existed in the pathology of pterygium in the Asian group. Furthermore, five gene sets (interferon-gamma response, Wnt beta-catenin signaling, oxidative phosphorylation, DNA repair, and MYC targets v2) were found only in our Chinese datasets. After taking an intersection between selecting datasets, we identified two upregulated (SPP1 and MYH11) and five downregulated (ATF3, FOS, EGR1, FOSB, and NR4A2) hub genes. We finally chose night genes to verify their expression levels, including the other two genes (SFRP2 and SFRP4) involved in Wnt signaling; Their expression levels were significantly different between pterygium and conjunctiva. Conclusions We consider hormone-related, metabolic, and Wnt signaling pathways may be important in the pathology of pterygium development. Nine candidate genes we identified deserve further study and can be potential therapeutic targets.


Introduction
Pterygium is a benign ocular surface disease that typically appears as a fibrovascular conjunctiva thickening from the nasal side and may progressively extend onto the cornea, causing visual impairment (1,2). The incidence of pterygium was 1.4% in the Asian population and 2.1% in the Chinese (3). Ultraviolet (UV) radiation, male gender, older age, ethnicity, outdoor occupation, and air pollution were associated with the occurrence and development of pterygium (3)(4)(5). Pterygium excision with conjunctival autograft and adjuvant therapy is the conventional treatment (6). Using fibrin glue in place of sutures for attaching conjunctival autograft during the surgery could further reduce the recurrence rate (0-4.5%) and surgical duration and help recover the cornea (7)(8)(9). Although many efforts have made progress, more surgical and non-surgical approaches are needed to be explored to decrease the recurrence further.
Currently, the underlying mechanism in pterygium remains indistinct. UV radiation is a well-recognized major cause of pterygium, and many research studies believe it triggers oxidative stress leading to DNA damage and subsequently stimulates multiple pathogenic factors (such as tumor suppressor P53, growth factors, pro-inflammatory cytokines, and extracellular matrix modulation) to induce a hyperproliferative state (10-15). On the other hand, hereditary factors and viruses (e.g., human papillomavirus) are also reported to contribute to the disease, so pterygium is considered an outcome of a combination of genetic and environmental factors (1,15). Moreover, pterygium is regarded as a limbal stem cell disorder with precancerous features (16). In short, despite many studies on pterygium, there are lots of unsolved puzzles that need to be unraveled.
Many studies have utilized transcriptome technologies to identify changes in transcripts of pterygium in recent years (17)(18)(19)(20)(21)(22)(23). They have greatly advanced the research on pterygium and identified several groups of abnormally expressed genes in pterygium, such as keratins (KRT3, KRT4, KRT6B, KRT13, KRT14, KRT16, and KRT24); S100 calcium-binding proteins (S100A8, S100A9, S100A11, and S100P); collagens (COL1A1, COL8A1, and COL10A1). Furthermore, other studies have reported the involvement of various matrix metalloproteinases in the progression of pterygium (MMP1, MMP2, MMP3, MMP9, and MMP13) (24)(25)(26). On the other hand, it has been reported that the Chinese have a higher incidence of the disease than the whole Asian population (3), so the transcriptomic information needs to be more abundant. Moreover, identifying the similarities and differences of transcripts between Asian studies is necessary to reveal the underlying mechanisms of pterygium.
In this study, we used microarray to evaluate the changes of transcripts between primary pterygium samples and adjacent normal conjunctiva samples in China and perform multiple bioinformatic analyses. Next, we took the intersection of differentially expressed genes with other data in Asian countries and identified hub genes. Finally, we verified the expression levels of candidate genes from microarray data. Our study hopes to provide more evidence of the progression of pterygium and find potential treatment targets.

Specimens preparation
This study approval was provided by the Ethics Committee of the second affiliated hospital of Anhui Medical University in accordance with the Helsinki Declaration of 1975. Sixteen nasal primary pterygium samples and fourteen adjacent normal conjunctiva samples were obtained from 16 participants (aged 47-70 years, five males and 11 females) through pterygium excision surgery. All the procedures were performed by an experienced ophthalmologist in the second affiliated hospital, Anhui Medical University. All of the participants in this study signed the informed consent before they were enrolled in the study. After pterygium excision, the whole part of the tissues was immediately placed in RNALater ™ RNA Stabilization Reagent for Animal Tissue (Beyotime, China, R0118) and infiltrated overnight at four degrees, then all of the tissues were stored at a −80°C refrigerator until RNA extraction.

RNA isolation
The total RNA of the tissues was isolated using TRIzol ™ Reagent (Invitrogen, USA, 15596026) according to the manufacturer's protocol. For microarray analysis, two pairs of pterygium and adjacent normal conjunctiva samples were used. The other tissues were used in real-time quantitative polymerase chain reaction (qPCR). The quality control and quantification of RNA were performed by electrophoresis on 1% agarose gels and NanoDrop 2000 (Thermo Fisher Scientific, USA).

Microarray analysis
After the quality test of RNA, gene expression profiling was conducted by Shanghai Baygene Biotechnology Co.Ltd (Shanghai, China) using GeneChip ® Human Transcriptome Array 2.0 (HTA 2.0). cDNA was synthesized, amplified, fragmented, and labeled for hybridization using the GeneChip ® WT PLUS Reagent Kit (Affymetrix, USA, 902280) from 2 ugs isolated RNA following the manufacturer's protocols. Then, GeneChip ® Expression Wash, Stain, and Scan were performed following protocol from GeneChip ® Hybridization Wash and Stain Kit (Affymetrix, USA, 900720). Probe cell intensity data (CEL) from microarrays are analyzed and normalized by Signal Space Transformation Robust Multichip Analysis (SST-RMA) algorithm in the Affymetrix ® Expression Console ™ software (v1.4). The normalized data were then subjected to R package limma (v3.42.2) to identify differentially expressed genes (DEGs) (27). Genes with |log2FC| >1 and p-value < 0.05 were identified as DEGs. Volcano plot visualization was exerted by ggplot2 (v3.3.0).
Next, the DEGs were submitted to execute the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses utilizing the R package ClusterProfiler (v3.14.3) (28). For Protein-protein interaction (PPI) network construction, we used the STRING database (https://string-db.org/cgi/input.pl) to recognize the possible connections between DEGs, the data exported from the STRING database was visualized by Cytoscape (v3.8.2), and the MCODE (v2.0.0) plugin of Cytoscape was applied to identify a key module of PPI network (29).

Data collection
For further bioinformatic analyses, we collected miRNA and microarray data from PubMed and the Gene Expression Omnibus database (GEO, www.ncbi.nlm.nih.gov/geo) based on our selection criteria. As for miRNA selection, we used "Pterygium" and "miRNA" as keywords to search on Pubmed and screen research from Asia. Next, we included 14 miRNAs that were validated by qPCR or could be found in at least two microarray data. The information on miRNA enrolled in our study is listed in Table 1 (30-37). The selection of microarray datasets for this study was based on screening the GEO database, and datasets from Asia were selected. The details of these microarray datasets are shown in Table 2.

The miRNA-mRNA regulatory network construction
All enrolled miRNAs were submitted to ENCORI (the encyclopedia of RNA Interactomes, https://starbase.sysu.edu. cn/) to predict and acquire their targeted mRNAs (38). The ENCORI could also connect to seven other target mRNApredicted programs (microT, miRanda, miRmap, PITA, RNA22, PicTar, and TargetScan); the mRNAs we selected were overlapped in at least four predicted programs. The miRNAs usually perform as negative regulators to affect the expression of the target mRNAs, upregulated miRNA led to downregulated target mRNA and vice versa. Considering such a regulatory pattern, we took the intersection of the acquired mRNAs with our microarray data based on the expression change direction. The visualization of the miRNA-mRNA r e g u l a t o r y n e t w o r k w a s c a r r i e d o u t b y T B t o o l s (v1.0986988) (39).

Gene set enrichment analysis
The gene set enrichment analysis (GSEA) based on Hallmark gene sets in the molecular signatures database (v7.5.1) was performed by GSEA software (v4.2.3). The results of GSEA were subjected to TBtools to visualize as a heat map, and data was clustered using Euclidean distances with complete hierarchical clustering.

Hub genes identification and enrichment analysis
To identify hub genes, we processed microarray datasets GSE83627, GSE51995, and GSE2513 by limma (v3.42.2). Besides, the microarray dataset GSE151872 only has one pterygium sample and is not applicable for limma's protocol, so we discarded this dataset to recognize hub genes. We took the intersection of DEGs from our microarray data and the above datasets, then submitted the results to the STRING database to perform further enrichment analysis. The visualizations were conducted by TBtools, Cytoscape, and ggplot2.

Real-time quantitative polymerase chain reaction
cDNA was synthesized from a total of 1 µg RNA with the ReverTra Ace ® qPCR RT Master Mix with gDNA Remover (Toyobo Co., Ltd., Japan, Code No.FSQ-301). Real-time quantitative polymerase chain reactions (qPCR) were carried out using SYBR ® Green Real-time PCR Master Mix (Toyobo Co., Ltd., Japan, Code No. QPK-201) with Roche LightCycler 480 (Roche Diagnostics Lid. Switzerland). Evaluating relative mRNA expression by using GAPDH as endogenous control via  the 2 −△CT methods. The experiments were repeated using three technical replicates. Table 3 lists all primers.

Statistical analysis
Analyses were performed and visualized by GraphPad Prism 8.3.0. After failing to pass the Shapiro-Wilk normality test, then comparisons of data were made using the Mann-Whitney nonparametric test. A p-value of <0.05 was considered a significant difference.

Results of microarray data analysis
3.1.1 The differential gene expression analysis The differential gene expression analysis identified a total of 49 DEGs (|log2FC| >1 and p-value < 0.05) on primary pterygium samples vs. adjacent normal conjunctiva samples. Among them, 21 were upregulated, and 28 were downregulated; all the DEGs and details are listed in Supplementary Table 1. The volcano plot (shown in Figure 1A) of the microarray data indicated the distribution of the DEGs. The upregulated DEGs were shown in yellow and downregulated were shown in green.

KEGG and GO pathway enrichment analysis
To get more bioinformatic signatures of the pterygium, we next performed enrichment analyses of the DEGs. Figure 1B presented the results of the KEGG pathway enrichment analysis. The upregulated pathways were related to infection and metabolism, and the downregulated were associated with inflammation, hormone, and cell differentiation. As for GO analysis, a total number of 195 GO terms of upregulated DEGs and 237 GO terms of downregulated DEGs were enriched (Supplementary Table 2 listed all the GO terms), the top 30 of them in different classifications (cellular component, molecular function, and biological process) were shown in Figures 1C, D for the up-and downregulated GO terms respectively. Negative regulation of the Wnt-signaling pathway, peroxide and oxygen related, and metabolism-related pathways could be found in upregulated GO terms ( Figure 1C). Additionally, hormone and skeletal muscle related and transcription factor associated pathways were enriched in downregulated GO terms ( Figure 1D).

PPI network construction
For protein-coding genes of the DEGs, we next put them into the STRING for a PPI network (Figure 2A). There were 32 nodes and 188 interacted edges in the PPI network. Then, we used MCODE to identify a key module for the whole network ( Figure 2B); this cluster scoring 8.444 had ten nodes and 76 edges and might play a role in the pterygium.

The miRNA-mRNA regulatory network construction
We searched on PubMed using the keywords "Pterygium" AND "miRNA" for all years, then selected the miRNAs in pterygium from Asian research based on the selection criteria. Table 1 summarizes the information on the selected miRNAs. Next, we put these miRNAs into ENCORI to predict the targeted mRNA and merged the results with our microarray data according to the screening criteria. Finally, a total of eight miRNAs and twelve mRNAs constructed 17 regulatory pairs. The miRNA-mRNA regulatory network was shown in Figure 3, and there existed a one-to-many regulatory relationship.

To compare the similarities and differences between Asian studies by gene set enrichment analysis
We found four datasets from Gene Expression Omnibus (GEO) in Asia (Singapore and Japan) to perform the following analysis (Table 2). We subjected the whole expression matrixes of our data and four datasets from GEO to GSEA software (v4.2.3), selecting Hallmark gene sets in the molecular signatures database (v7.5.1) to find similar and different biological processes of pterygium in Asia. As shown in Figure 4, only two gene sets were enriched in all datasets (epithelial-mesenchymal transition and myogenesis). Three gene sets were in the results of the other four datasets but not ours (inflammatory response, angiogenesis, and coagulation). Moreover, five gene sets were found only in our microarray data (interferongamma response, Wnt beta-catenin signaling, oxidative phosphorylation, DNA repair, and MYC targets v2). Furthermore, there were several gene sets in our data that were enriched in the other three data. Our microarray data was indicated as a separate clustering in cluster analysis based on Euclidean distance. The complete results of GSEA are demonstrated in Supplementary  Table 3. The bigger the NES, the stronger the effect. NES > 0 meant upregulated while NES < 0 represented downregulated.

Identifying the hub genes by taking intersection with other microarray data
To identify the hub genes in the development of pterygium, we took an intersection between our data and other data from GEO datasets (Figures 5A, B). There were only two upregulated and five  downregulated genes at the intersection of four microarray datasets ( Table 4). All of them were protein-coding genes. Strikingly, all five downregulated genes were in the key cluster that MCODE identified ( Figure 2B), so we next put these hub genes into STRING again and performed network enrichment analysis. Two upregulated genes showed no interaction, while five downregulated genes displayed a strong connection ( Figure 5C). Moreover, the network enrichment analysis (Figures 5D, E) of these five genes showed similarity with enrichment analysis of downregulated DEGs (e.g., Parathyroid hormone synthesis, secretion, and action; IL-17 The PPI network and module analysis.  The miRNA-mRNA regulatory network. Eight miRNAs and 12 mRNAs construct 17 regulatory pairs.

Validation of the mRNA expression levels of candidate genes by qPCR
Based on the previous analysis, we chose seven hub genes to verify their expression levels. We enrolled SFRP2 and SFRP4 in further analysis because they were DEGs involved in Wnt signaling. As illustrated in Figures 6A-I, the mRNA expression levels of candidate genes (SFRP2, SFRP4, SPP1, ATF3, FOS, EGR1, FOSB, NR4A2, and MYH11) were significantly different in primary pterygium samples compared to conjunctiva samples. Additionally, the correlation analysis of log2 foldchange between microarray and the qPCR quantitative data demonstrated a high correlation (R 2 = 0.8934, Figure 6J) which presented the accuracy of our microarray data.

Discussion
Pterygium is a common fibrovascular degeneration of conjunctiva, and ultraviolet radiation is the main cause. The best way of treating pterygium is to perform surgery with some adjuvant treatments (6,40). Currently, the concrete mechanisms of onset and progression of pterygium still remain obscure.   To explore the underlying mechanism of pterygium, we used microarray to compare changes in transcripts between primary pterygium samples and conjunctiva samples. 49 DEGs were identified in our study, and we next conducted multiple enrichment analyses to find the potential pathways of the disease. In the results, numbers of the metabolic and biosynthetic processes were upregulated (seven of the top ten GO BP terms), and these processes represent that the cells in pterygium were actively reacting to the compounds (e.g., hydrogen peroxide) involved in them. As part of reactive oxygen species, hydrogen peroxide is responsible for oxidative stress, while oxidative stress is considered a pathogenic mechanism of pterygium (12,41). Additionally, metabolites from healthy and diseased cells can be directly released to tear, making the difference in tear proteome between pterygium and other diseases (42). The Wnt-signaling pathway, which can control stem cell biology and growth, is also upregulated; the changes in gene expression associated with Wnt-signaling in pterygium have been reported; however, how this pathway affects the development of pterygium needs to be studied in depth (43,44). Downregulated pathways also provide us with pieces of evidence of disease progression. Hormone-related pathways are significantly downregulated; some research reported that hormone-related growth factors and receptors are implicated in pterygium (45,46). Even though some hormones can affect corneal morphology, physiology, and metabolism, there is few of evidence indicating that hormones directly participate in the development of pterygium (47,48). On the other hand, long-term UVA radiation to the eyes can affect the level of corticotropin-releasing hormone in the brain (49). So, we infer that the changes in hormone-related pathways are the response to UVA radiation. Furthermore, growth factors are also downregulated. Much evidence suggests that multiple growth factors are involved in pterygium, and their expression pattern may lead to a wide variation in the growth of pterygium (15). Other downregulated pathways (e.g., tumor necrosis factor (TNF) signaling, mitogen-activated protein kinase (MAPK) signaling, and transcription factor) have strong relationships with UV radiation or pterygium (50)(51)(52). The above results, to some extent, confirmed the reliability of our microarray data.
Nowadays, the PPI network and miRNA-mRNA regulatory network have been constructed in many studies (53)(54)(55). Since our new microarray data and new miRNA data from PubMed became available, we constructed new networks based on recent data. Our PPI network and miRNA-mRNA regulatory network hope to provide a new basis for further study of pterygium. Moreover, hsa-miR-199a-3p targets DUSP5 in our network have been confirmed in pterygium (33).
Besides performing analyses of our microarray data, we also compare our data with other Asian datasets from GEO. The hallmark gene sets include 50 gene sets, while 45 gene sets were enriched through GSEA. As displayed in Figure 4, epithelialmesenchymal transition (EMT) and myogenesis are found upregulated in all datasets. EMT, playing an essential role in wound healing and tissue remodeling, is a biological process that induces the transition of polarized and immotile epithelial cells into motile mesenchymal fibroblast-like cells (56). EMT and myogenesis in our results represent epithelial cells converted to myofibroblasts, serving as a key feature of pathological tissue repair and playing an essential role in pterygium progression (33,57). As for the differences between our data and other Asian data, we found five significant pathways only upregulated in our results: interferon-gamma response, Wnt beta-catenin signaling, oxidative phosphorylation, DNA repair, and MYC targets v2. As outlined previously, oxidative phosphorylation and DNA repair are related to oxidative stress and DNA damage which are downstream of UV radiation. Interferongamma response belongs to immune responses which have been associated with pterygium (58). MYC, as a proto-oncogene is also involved in the disease (59). The result of Wnt beta-catenin signaling is consistent with our previous GO analysis. On the other hand, three pathways (inflammatory response, angiogenesis, and coagulation) correlating with pterygium are not present in our data but exist in the other four datasets (15). Interestingly, three datasets show upregulation of the pathways named UV response downregulated. This result probably shows that the response to UV radiation is different between pterygium cells and healthy cells, while pterygium epithelium possibly has resistance to UV-induced apoptosis (14). In conclusion, the results of GSEA have offered a shred of potential evidence for the cause of differences in pterygium incidence in Asia. The pathways that are unique and absent in our datasets should be worth attention.
To further investigate the common mechanism of pathogenesis in pterygium, we took the intersection between our data and other GEO datasets we selected. A total of seven genes were identified as hub genes. The two upregulated genes show no interactions, and the network enrichment analysis shows no results. We later subjected them to ClusterProfiler, and the results were similar to previous GO terms that they involved (data not shown). SPP1 encodes osteopontin playing important roles in wound healing. Osteopontin is suggested to affect metalloproteinase (MMP) secretion and influence fibroblast proliferation dependent on the growth factors (60). This process is highly concordant with mechanisms that are reported in the progression of pterygium (25). MYH11 encoding myosin heavy chain 11 has been confirmed its involvement in several types of cancer and is probably related to metabolism-related genes (61). Five downregulated hub genes are included in the key module that MCODE identified and have interactions with each other, which means they are in the important section of disease development. These five hub genes are transcription factors. Network enrichment analysis shows they are also significantly associated with hormone-related pathways. Among them, FOS and FOSB, as members of the AP-1 transcription factor complex, were upregulated in corneal epithelial cell layers when under UV exposure for one to six hours (62). Furthermore, the c-fos mRNA level is intermediateearly induced by UV radiation in cultured pterygium epithelial cells but downregulated after 12h (52). Additionally, after applying oxidative stress to human retinal pigment epithelium, the changes of expression of ATF3 (another member of AP-1 transcription factors), FOS, and FOSB are in a dose-dependent manner (63). Moreover, overexpression of SPP1 in cells can downregulate the level of c-fos, and SPP1 plays a role in wound healing (64). So, since AP-1 transcription factors are downregulated in this study, we infer that their levels depend on the balance of effects of UV exposure and wound healing. EGR1 encodes early growth response, critically participating in neovascularization, tumor angiogenesis, and growth, together with ATF3, which can influence the limbal epithelial cell proliferation (65,66). NR4A2 is also associated with tumor proliferation, migration, and invasion (67). In summary, we speculate these hub genes are associated with pterygium proliferation and may response for resistance to UV exposure, these effects may be involved in inactivation of hormonerelated pathways.
Finally, we verified the mRNA levels of nine candidate genes through qPCR. Correlation analysis shows the qPCR data are consistent with microarray data, suggesting the credibility of our results. Besides seven hub genes, we also validated two upregulated gene expressions (SFRP2 and SFRP4). They are served as classical antagonists of Wnt signaling and can interact with Wnt protein directly through both autocrine and paracrine modes. In some cases, they also can promote Wnt signaling (68). Wnt signaling can control adult stem cell biology and growth (69); meanwhile, pterygium is thought to be a limbal stem cell disorder (16). The results of these two genes remind us that genes located upstream of the Wnt signaling pathway are involved in the progression of pterygium and probably can be used as therapeutic targets. Additionally, the Wnt beta-catenin signaling is the unique gene set in the results of GSEA in our data, so further in-depth studies are required, especially for the Chinese. Three hub genes (FOS, NR4A2, and SFRP2) are also in our miRNA-mRNA network; the predicted regulatory relationships are worthy of future research.
There were some deficiencies in our study. 1) Our sample size was insufficient and may cause a potential bias. 2) The validations of protein levels of our candidate genes were not conducted. The following study is currently ongoing to remedy the shortcomings of this study.
In conclusion, we consider hormone-related, metabolic, and Wnt signaling pathways may be important in developing pterygium besides its classical mechanisms. Nine candidate genes we identified deserve further study and can be potential therapeutic targets. In short, our findings shed important light on the further study of pterygium.

Data availability statement
The data presented in the study are deposited in the GEO repository, accession number GSE208384.

Ethics statement
The studies involving human participants were reviewed and approved by The Ethics Committee of The Second Affiliated Hospital of Anhui Medical University. The patients/participants provided their written informed consent to participate in this study.

Author contributions
JL and TT wrote this manuscript and performed bioinformatic analyses. YY, NX, and WD conducted data collection and experiment process. ZJ carried out a surgical procedure. MZ and LH provided the idea for this study, and LH revised this manuscript. JL and TT have contributed equally to this work. All authors contributed to the article and approved the submitted version.