Identification of the circRNA–miRNA–mRNA Regulatory Network in Pterygium-Associated Conjunctival Epithelium

To investigate the regulatory mechanism of pterygium formation, we detected differentially expressed messenger RNAs (DE-mRNAs) and differentially expressed circular RNAs (DE-circRNAs) in pterygium-associated conjunctival epithelium (PCE) and normal conjunctival epithelium (NCE). Genome-wide mRNA and circRNA expression profiles of PCE and NCE were determined using high-throughput sequencing. Bioinformatics analyses, including Gene Ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, gene set enrichment analysis (GSEA), and protein–protein interaction (PPI) analysis, were conducted. The microRNAs (miRNAs) interacting with the hub DE-mRNAs and DE-circRNAs were predicted and verified using real-time quantitative PCR (RT-qPCR). The data showed that there were 536 DE-mRNAs (280 upregulated and 256 downregulated mRNAs) and 78 DE-circRNAs (20 upregulated and 58 downregulated circRNAs) in PCE. KEGG enrichment analysis indicated that the DE-mRNAs were mainly involved in the following biological processes: IL-17 signalling pathway, viral protein interaction with cytokine and cytokine receptor, cytokine–cytokine receptor interaction, ECM-receptor interaction, and focal adhesion. The GSEA results revealed that the epithelial mesenchymal transition (EMT) process was significantly enriched in upregulated mRNAs. The pterygium-associated circRNA–miRNA–mRNA network was established based on the top 10 DE-circRNAs, 4 validated miRNAs (upregulated miR-376a-5p and miR-208a-5p,downregulated miR-203a-3p and miR-200b-3p), and 31 DE-mRNAs. We found that miR-200b-3p, as a regulator of FN1, SDC2, and MEX3D, could be regulated by 5 upregulated circRNAs. In addition, we screened out EMT-related DE-mRNAs, including 6 upregulated DE-mRNAs and 6 downregulated DE-mRNAs. The EMT-related circRNA–miRNA–mRNA network was established with the top 10 circRNAs, 8 validated miRNAs (upregulated miR-17-5p, miR-181a-5p, and miR-106a-5p, downregulated miR-124-3p, miR-9-5p, miR-130b-5p, miR-1-3p, and miR-26b-5P), and 12 EMT-related DE-mRNAs. We found that hsa_circ_0002406 might upregulate FN1 and ADAM12 by sponging miR-26b-5p and miR-1-3p, respectively, thus promoting EMT in pterygium. Briefly, the study provides a novel viewpoint on the molecular pathological mechanisms in pterygium formation. CircRNA–miRNA–mRNA regulatory networks participate in the pathogenesis of pterygium and might become promising targets for pterygium prevention and treatment.


Introduction
Pterygium is characterized by fibrovascular tissue hyperplasia from the bulbar conjunctiva towards the cornea, accompanied by foreign body sensation, dry eye, astigmatism, and visual impairment [1,2]. Once pterygium covers the pupillary area, it will cause a significant sight-threatening complication. The prevalence rates of pterygium were 8.8% in South Korea, 9.84% in China, and 38.7% in Northwest Ethiopia [3]. Currently, surgical excision is the main therapeutic method for pterygium. Nevertheless, recurrence is common. Excision combined with conjunctival autografts or amniotic membrane grafts is used to decrease postoperative recurrence [4]. However, the pterygium recurrence rate is still 3.3%-16.7% for conjunctival autografts and 6.4%-42.3% for amniotic membrane grafts [4]. Therefore, elucidating the pathogenesis and molecular mechanisms of pterygium is of great significance for preventing pterygium growth and recurrence.
Previous studies have demonstrated that inflammation, epithelial mesenchymal transition (EMT), DNA repair, cell proliferation, cell migration, and angiogenesis contribute to the pathogenesis of pterygium [5,6]. Recently, an increasing number of studies on noncoding RNAs have demonstrated that the mechanism of pterygium is complex [7,8].
Noncoding RNAs (ncRNAs), such as circular RNAs (cir-cRNAs) and microRNAs (miRNAs), are involved in the transcriptional regulation of gene expression. CircRNAs are covalently closed continuous loops. CircRNAs can inhibit miRNA functions by sponging the target miRNAs directly or indirectly. The miRNAs are single-stranded noncoding RNAs composed of 18-24 nucleotides. Previous studies have shown that miRNAs downregulate gene expression by suppressing target mRNA translation or promoting mRNA degradation [9]. As reported, miRNAs significantly contribute to pterygium development by interacting with pterygium-associated mRNAs [10][11][12][13][14]. Moreover, circRNAs acting as miRNA sponges are differentially expressed in pterygium in comparison to conjunctival tissues [15]. Therefore, we speculated that the endogenous RNA regulatory network may play a key role in pterygium formation and development. However, integrative analysis of the cir-cRNA-miRNA-mRNA regulatory mechanism in pterygium remains lacking. Ascertaining the circRNA-miRNA-mRNA regulatory network in pterygium is essential for preventing pterygium formation, growth, and recurrence.
In this study, as shown in Figure 1, we explored differentially expressed messenger RNAs (DE-mRNAs) and differentially expressed circRNAs (DE-circRNAs) in pterygiumassociated conjunctival epithelium (PCE) and normal conjunctival epithelium (NCE) using high-throughput sequencing. Afterwards, the circRNA-miRNA-mRNA networks were established to help us understand the pathogenesis and underlying molecular mechanisms of pterygium. epithelium in the pterygium. Three NCE samples were collected from the nasal conjunctival epithelium of donated eyes. The study was authorized by the Ethics Review Committee of the Affiliated Hospital of Nantong University and followed the Declaration of Helsinki. The participants signed the informed consent form. The samples were preserved at -80°C until further processing.

Preparation of Sequencing
Libraries and Sequencing of mRNAs and circRNAs. Briefly, total RNA was purified from the samples with the RNeasy Kit (Qiagen, Duesseldorf, Germany). After qualification of RNAs using the NanoDrop One system (Waltham, MA, USA), the sequencing libraries were established by the Whole RNA-seq Lib Prep Kit (ABclonal, Shanghai, China). The sequencing libraries were prepared and quantified via a Qubit 3.0 fluorometer (Invitrogen, Carlsbad, USA). Finally, the library sequencing was performed by an Illumina NovaSeq sequencer (Illumina, CA, USA).
The Arraystar Super RNA Labeling Kit was used in cir-cRNA amplification and fluorescent cRNA transcription. Next, the acquired cRNAs were transcribed into the Human circRNA Array (8 ×15K, Arraystar). Then, the circRNA expression level was detected with "mapped back-splicing junction reads per million mapped reads" after circRNA verification.

DE-mRNAs and DE-circRNAs Identification.
The lowquality reads for base quality lower than 20 were eliminated by Cutadapt software (v1.9.1) [16]. Gene annotation and reference genome files were obtained from online databases (UCSC and NCBI). Next, HISAT2 software (v2.0.1) was used to align the clean data to the reference genome [17]. The transcripts were annotated and indexed for expression analysis. Then, we used an annotated file as a reference gene set and estimated the related gene expression using RSEM software (v1.2.15) [18]. The differential expression was analysed using the Bioconductor package DESeq2 [19]. The mRNAs or circRNAs with P < 0:05 and jlog 2 ðfoldchangeÞj ≥ 1 were identified as DE-mRNAs or DE-circRNAs.      BioMed Research International  BioMed Research International molecular function (MF), and enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were performed to reveal the functions of DE-mRNAs and DE-circRNAs. Enrichment analysis of GO categories was performed by R "cluster Profiler" package, and enrichment analysis of KEGG pathways was tested upon hypergeometric distribution by R "phyper" function [20]. Those GO categories with adjusted P < 0:05 and KEGG pathways with P < 0:05 were considered as significantly enriched. Meanwhile, GO categories or KEGG pathways with less than 3 DEGs were discarded.

BioMed Research International
2.5. Protein-Protein Interaction (PPI) Analysis. The pathway database (http://www.pathwaycommons.org/) was used to analyse the PPI network in DE-mRNAs [21]. We selected direct interactions among DE-mRNAs. The R igraph package was applied to identify whether the genes were significant according to three attributes (event, betweenness, and degree) in the PPI network [22]. Genes involved in identical biological processes were grouped. The ggplot2 package was used to plot the distribution of gene attributes [23].
2.6. Gene Set Enrichment Analysis (GSEA). To avoid omission of key pathways and hub mRNAs in DE-mRNAs, GSEA was performed. The analysis involved 1,000 permutations of the gene set. The following criteria were applied to

Verification of Hub miRNAs and Hub circRNAs.
The expression of hub miRNAs that could interact with hub mRNAs and circRNAs was detected in the samples by realtime quantitative PCR (RT-qPCR). Total RNA was purified from those samples. In the case of miRNA PCR, the miRNA RT-qPCR Starter Kit (RiboBio, Guangzhou, China) was used to reverse-transcribe miRNAs following the manufacturer's instructions. Total miRNAs were reverse-transcribed to cDNAs. Then, PCR was performed using SYBR Green Master Mix (Vazyme, Nanjing, China) with miDETECT A Track miRNA primers (RiboBio, Guangzhou, China) in QuantStudio 5 Real-Time PCR Systems (Applied Biosystems, Foster City, CA, USA). Relative miRNA expression was normalized to U6. In the case of circRNA PCR, the reaction was carried out using the RiboBio circRNA qRT-PCR Starter Kit (Ribo-Bio, Guangzhou, China) and SYBR Green Master Mix (Vazyme, Nanjing, China). Relative circRNA expression was normalized to GAPDH expression. Significant miRNAs and circRNAs were identified depending on previous reports and predicted results.
2.9. Statistical Analysis. The limma package of R 3.6.3 was applied for data standardization [26]. The DE genes were screened following a standardized expression matrix. Moreover, when jlog 2 ðfoldchangeÞj ≥ 1 and P < 0:05, the expression levels of mRNAs and circRNAs were deemed to be significantly different. The pheatmap package and ggplot2 package in R 3.6.3 were used in volcano plots and heatmap construction [27]. P < 0:05 was considered to indicate statistical significance. Further analysis was performed using GraphPad Prism software 8.0.

DE-mRNA Screening.
The results showed that 19467 mRNAs were identified by high-throughput sequencing in PCE and NCE. A total of 536 mRNAs (jlog 2 foldchangej ≥ 1, P < 0:05) were deemed to be differentially expressed. The mRNA expression patterns were shown using hierarchical clustering analysis (Figure 2(a)). Differences in mRNA expression between PCE and NCE were evaluated by volcano plot analysis (Figure 2(b)). Among those DE-mRNAs, 280 mRNA were upregulated, and 256 mRNAs were downregulated in PCE compared with NCE ( Figure 2(b)).

Functional Analysis of DE-mRNAs. GO and KEGG analyses were used to reveal the underlying roles of DE-mRNAs.
Among all DE-mRNAs, the significantly enriched BPs included complement activation classical pathway, regulation of complement activation, response to external stimulus, positive regulation of response to stimulus, and vesicle-  (Figure 3(a)). The main enriched CCs included immunoglobulin complex, immunoglobulin complex circulating, extracellular matrix, and endocytic vesicle lumen (Figure 3(a)). The main enriched MFs included antigen binding, immunoglobulin receptor, signalling receptor binding, glycosaminoglycan binding, and chemokine receptor binding (Figure 3(a)). The main enriched KEGG pathways of all DE-mRNAs are shown in Figure 3(b). The DE-mRNAs were primarily involved in the following pathways: IL-17 signalling pathway, viral protein interaction with cytokine and cytokine receptor, cytokine-cytokine receptor interaction, ECM-receptor interaction, focal adhesion, tumour necrosis factor (TNF) signalling pathway, chemokine signalling pathway, and PI3K-Akt signalling pathway (Figure 3(b)).
To further study the pathogenesis of pterygium, we overlapped the DE-mRNAs with pterygium-associated genes. We primarily screened for crucial genes involved in pterygium-associated biological processes (oxidative stress, DNA injury, DNA repair, inflammation, cell proliferation, autophagy, apoptosis, ferroptosis, angiogenesis, and so on). Briefly, 145 pterygium-associated DE-mRNAs were shown using a Venn diagram (Figure 4(a)). The two-dimensional histogram shows the numbers of DE-mRNAs in the 14 BioMed Research International processes of cell death, DNA repair, EMT, angiogenesis, inflammation, and cell proliferation (Figure 4(b)). Moreover, we showed pterygium-associated DE-mRNAs which were involved in the above pathophysiological processes (Figures 4(c)-4(h)).

PPI Network
Analysis. The PPI network of DE-mRNAs was established ( Figure 6). The PPI network revealed that the DE-mRNAs were involved in 12 key pterygiumassociated pathways. The MCODE plugin was used to analyse the significant module ( Figure 6).

Screening and Functional Annotation of circRNAs.
Cir-cRNAs can act as miRNA sponges to modulate gene expression [28]. To explore the potential pterygium-related circRNAs, we detected the circRNA expression profiles of PCE and NCE using high-throughput sequencing. The DE-circRNAs were clustered into two groups (Figure 7(a)). A volcano plot was used to filter DE-circRNAs (jlog 2 ðfold changeÞj ≥ 1, P < 0:05), revealing 20 upregulated and 58 downregulated DE-circRNAs in PCE compared with NCE (Figure 7(b)). GO analysis of DE-circRNA host genes was performed in three dimensions, including BP, CC, and MF. In the BP analysis, the main enriched categories included organic substance catabolic process, nuclear export, response to temperature stimulus, regulation of neuron projection development, and proteolysis (Figure 7(c)). In the CC analysis, the main enriched categories included the endoplasmic reticulum membrane and cytoplasmic vesicle membrane (Figure 7(c)). In the MF analysis, the main enriched categories included GTPase regulator activity, Rab GTPase binding, endopeptidase activity, phosphoric ester hydrolase activity, and phospholipid binding (Figure 7(c)).

Discussion
Pterygium is a common ocular disease characterized by hyperplasia of conjunctival tissues on the cornea. The pterygium is mainly composed of vessels, fibroblasts, and epithelium. The pathogenesis and molecular mechanism of pterygium are still unclear. It was reported that EMT plays a key role in the pathogenesis of pterygium [29,30].
Previous studies have focused on pterygium tissues (including pterygium epithelium and stroma); however, pterygium-associated EMT mainly occurs in pterygium epithelium. Therefore, we carefully dissected and sequenced the pterygium epithelium instead of the entire pterygium tissue. Moreover, previous studies have identified the roles of the lncRNA-miRNA network in pterygium, but very few studies have explored the circRNA-miRNA-mRNA regulatory network [7,8]. In this study, we established the circRNA-miRNA-mRNA network in pterygium epithelium and provided deeper understanding of pterygium. The workflow of study design is shown in Figure 1.
Functional enrichment analysis demonstrated that DE-mRNAs were significantly enriched in the PI3K-Akt signalling pathway, ECM receptor interaction process, and haematopoietic cell lineage process, corresponding to previous studies [8,37]. Bioinformatics analysis also indicated that EMT-related processes were significantly enhanced in PCE. It has been demonstrated that pterygium formation is closely connected with EMT of the conjunctiva [35]. Conjunctival epithelial cells can transform into pterygium stromal cells via the EMT process [38]. Moreover, EMT is also an important mechanism contributing to the migration and invasion of pterygium cells [39]. In this study, we found that EMT-related gene expression was significantly increased in the pterygium epithelium. To further investigate the pathway, we assessed the enrichment of DE-mRNAs in the EMT-related pathways. Among the EMTrelated DE-mRNAs, MMP-3, FN1, TNC, and NNMT were significantly upregulated in PCE.
Studies have demonstrated that noncoding RNAs are involved in regulating pathogenic genes in pterygium [8,40]. In the present research, we found that 20 DE-circRNAs were upregulated and 58 DE-circRNAs were downregulated in PCE compared with NCE. After predicting miRNAs interacting with circRNAs, we established the circRNA-miRNA-mRNA network, including the top 10 DE-circRNAs and top 10 DE-mRNAs. In addition, we verified 2 upregulated DE-miRNAs (miR-376a-5p and miR-208a-5p) and 2 downregulated DE-miRNAs (miR-203a-3p and miR-200b-3p) in the regulatory network using qRT-PCR. It has been previously discovered that miR-200b-3p was downregulated in pterygium as a regulator of FN1 [7,8]. MEX3D and SDC2 were also important targets of miR-200b-3p. It was reported that MEX3D could promote proliferation of cervical carcinoma and transformation of prostatic epithelium in prostate cancer [41,42]. SDC2 is involved in the EMT process of prostate cancer [43]. Moreover, we found that miR-200b-3p was negatively associated with five upregulated circRNAs, including hsa_circ_ 0002406, hsa_circ_0002564, hsa_circ_0072688, 6:32610387-32713849, and 19:11941432-12058122. We speculated that these highly expressed circRNAs may increase the expression of FN1, SDC2, and MEX3D by sponging miR-200b-3p and further promote pterygium formation and growth.
Since EMT is a key mechanism of pterygium, we also constructed an EMT-related circRNA-miRNA-mRNA network involving 10 DE-circRNAs (5 upregulated and 5 downregulated DE-circRNAs), 8 validated DE-miRNAs, and 12 DE-mRNAs. We found that the downregulation of miR-1-3p and miR-26b-5p could simultaneously induce the upregulation of FN1 and ADAM12. Previous studies have demonstrated that ADAM12 can promote tumour invasion and EMT [44,45]. FN1 has been reported in pterygium, as mentioned above. In addition, hsa_circ_ 0002406 was significantly upregulated and found to target seven miRNAs, including miR-1-3p and miR-26b-5p. Therefore, we speculated that hsa_circ_0002406 might promote the EMT process by sponging miR-1-3p/miR-26b-5p and then upregulating FN1 in pterygium epithelium. The functions of hsa_circ_0002406 in pterygium have not been reported previously and need to be further studied.
In conclusion, we identified DE-mRNAs and DE-circRNAs in PCE using high-throughput sequencing and then established circRNA-miRNA-mRNA regulatory networks. We found that circRNAs may regulate EMT-related genes by targeting miRNAs and play an important role in pterygium. The limitation of this study is that the content of circRNA databases is not sophisticated enough, and we could not fully elucidate the functions of all DE-circRNAs. In the future, we will further verify the regulatory functions of filtered key DE-circRNAs and elucidate the underlying molecular mechanisms in pterygium occurrence and growth.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Ethical Approval
The study followed the tenets of the Declaration of Helsinki and was authorized by the Ethics Committee of the Affiliated Hospital of Nantong University.