Screening of DEmRNAs
After high-throughput sequencing, 19467 mRNAs were detected between NCE and PCE samples (n=3). In addition, 536 mRNAs (log fold change>2.0, P<0.001) were identified to be differentially expressed. Hierarchical clustering analysis of gene expression showed the mRNA expression patterns (Fig. 1a). A volcano plot of the mRNA expression profile was also used to evaluate variations in mRNA expression between these two groups (Fig. 1b). A total of 536 DEmRNAs were detected; 280 of these mRNAs were upregulated, and 256 were downregulated in PCE samples compared with NCE samples. Moreover, the statistically significant DEmRNAs were identified (Fig. 1b).
Functional Annotation of the DEmRNAs
To explore the potential functions of DEmRNAs, we performed GO analysis and KEGG pathway enrichment analysis. The results showed that the significantly enriched BP terms among all DEmRNAs included the regulation of protein metabolism, vesicle-mediated transport, programmed cell death, tissue development, complement activation in the classical pathway, regulation of complement activation, and the response to wounding (Fig. 2a). The main CCs were the response to the immunoglobulin complex, circulating immunoglobulin complex, collagen-containing extracellular matrix, and endocytic vesicle lumen (Fig. 2b). The main MFs were the response to signalling receptor binding, antigen binding, structural molecule activity, immunoglobulin receptor binding, glycosaminoglycan binding, and organic acid binding (Fig. 2c).
The results of KEGG signalling pathway enrichment analysis of all DEmRNAs were shown in Fig. 2d. We found that the upregulated genes were mainly enriched in ECM-receptor interaction, the PI3K-Akt signalling pathway, the haematopoietic cell lineage, neuroactive ligand-receptor interaction, and metabolism of xenobiotics by cytochrome. In contrast, the downregulated genes were mainly related to viral protein interactions with cytokines, cytokine-cytokine receptor interactions, the TNF signalling pathway, and the IL-17 signalling pathway.
Functional Annotation of the DEGs Related to the Pathogenesis of Pterygium
To further explore the pathogenesis of pterygium in PCE tissues, we overlapped the DEmRNAs with genes correlated with the pathogenesis of pterygium. We mainly screened for genes related to inflammation, oxidative stress, autophagy, proliferation, DNA damage and repair, angiogenesis, lymph angiogenesis, apoptosis, autophagy and ferroptosis. We screened 145 DEGs related to the mechanism of pterygium using a Venn diagram and a two-dimensional histogram (Fig. 3a, 3b). In addition, we presented the six pterygium pathogenesis-related pathway genes using heat maps (Fig. 3c-g).
GSEA of mRNAs Related to Pterygium
GSEA was then performed to identify gene sets enriched in PCE samples. GSEA revealed that coagulation, EMT, Kras signalling and peroxisomes were significantly enriched among the upregulated DEGs (P<0.001, and NES ≥ 1) (Fig. 4). The pathways significantly enriched among the downregulated DEGs included TNFα signalling via NFκB, the inflammatory response, hypoxia, the interferon alpha response, allograft rejection, hedgehog signalling, the UV response and so on (P<0.05, and NES ≤ -1) (Fig. 4).
PPI Network Visualization
The DEmRNAs were analysed for gene-gene interaction networks via the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online tool (Fig. 5). The PPI network reflects the density distribution of gene betweenness in a biological process and contains 12 important pathways related to pterygium development. Subsequently, the DEmRNAs with a degree of connectivity ≥ 5 were selected as hub genes. In this study, the significant module was analysed by the MCODE plug-in (Fig. 5).
miRNA-mRNA Regulatory Network
To gain insight into the molecular mechanisms, we predicted the miRNAs of hub genes through miRbase. First, we constructed miRNA regulatory networks for EMT pathway-related DEmRNAs (Fig. 6). We also divided the target genes acting on miRNAs into two groups, in which one group was related to the pathogenesis of pterygium and the other group included DEGs (Fig. 7A). Subsequently, the associations between the 48 mRNAs targeted by 19 miRNAs were predicted. Finally, 19 miRNAs and 48 DEmRNAs were used to construct the miRNA-mRNA network (Fig. 7B).
Validation of Key miRNAs
In total, 11 predicted miRNAs (miR-5579-3p, miR-653-3p, miR-3617-5p, miR-4681, miR-4802-3p, miR-6844, miR-888-3p, miR-7106-3p, miR-200b-3p, miR-200a-3p and miR-543) were verified through RT-qPCR. Although these three miRNAs have been reported in diabetic retinopathy, they have not been reported in pterygium. The results showed that eight downregulated miRNAs (miR-653-3p, miR-3617-5p, miR-4681, miR-4802-3p, miR-888-3p, miR-7106-3p, miR-543, and miR-200b-3p) and two upregulated miRNAs (miR-5579-3p and miR-6844) were significantly differentially expressed (P < 0.001) (Fig. 8A). We screened three miRNAs studied in diabetic retinopathy from the top 20 miRNAs predicted to be related to the pathogenesis of pterygium. In our results, the expression of all three miRNAs was downregulated (Fig. 8b).