ORIGINAL RESEARCH article
Comparative Transcriptomic Analysis to Identify the Important Coding and Non-coding RNAs Involved in the Pathogenesis of Pterygium
- 1Shenzhen Eye Hospital, Shenzhen Key Laboratory of Ophthalmology, Affiliated Shenzhen Eye Hospital of Jinan University, Shenzhen, China
- 2Department of Ophthalmology, The Affiliated Hospital of Guizhou Medical University, Guiyang, China
Pterygium is a common ocular surface disease characterized by abnormal fibrovascular proliferation and invasion, similar to tumorigenesis. The formation of tumors is related to a change in the expression of various RNAs; however, whether they are involved in the formation and development of pterygium remains unclear. In this study, transcriptome analysis of messenger RNAs (mRNAs), long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs) of paired pterygium and normal conjunctiva was performed to explore key genes regulating the development of pterygium. In total, 579 mRNAs, 275 lncRNAs, and 21 circRNAs were differentially expressed (DE) in pterygium compared with paired conjunctival tissues. Functional enrichment analysis indicated that DE RNAs were associated with extracellular matrix organization, blood vessel morphogenesis, and focal adhesion. Furthermore, through protein-protein interaction network and mRNA-lncRNA co-expression network analysis, key mRNAs including FN1, VCAM1, and MMP2, and key lncRNAs including MIR4435-2HG and LINC00968 were screened and might be involved in the pathogenesis of pterygium. In addition, several circRNAs including hsa_circ_0007482 and hsa_circ_001730 were considered to be involved in the pterygium development. This study provides a scientific basis for elucidating the pathogenesis of pterygium and will be beneficial for the development of preventive and therapeutic strategies.
Pterygium is a common ocular surface disease characterized by benign subconjunctival fibrous connective tissue hyperplasia and vascular invasion into the limbus, causing eye dryness, foreign body sensation, and even visual loss. Presently, there is no effective drug for the treatment of this disease, and surgery is the best treatment option; however, recurrence is common. Repeated recurrence may cause symblepharon adhesion and affect eye movement. Most researchers have reported that the pathogenesis of pterygium is related to an exposure to UV radiation, and occurrence of anti-apoptosis mechanisms, extracellular matrix remodeling, and other mechanisms, but the exact cause of pterygium remains unclear (Liu et al., 2017b). Pathologically, pterygium is characterized by proliferation, fibrosis, corneal infiltration, and angiogenesis, which are similar to those observed in tumor tissues (Chui et al., 2008). Epithelial-mesenchymal transition (EMT), a recent focus in the study of tumor pathogenesis, was reported to be involved in pterygium generation (Kato et al., 2007; Meshkani et al., 2019). The expression levels of EMT markers, such as FN1 and VIM, were significantly higher in pterygium than those in normal conjunctiva (Engelsvold et al., 2013; Hou et al., 2014). Additionally, angiogenesis, an important pathological mechanism of disease development, participates in the occurrence and development of pterygium (Livezeanu et al., 2011). The formation of pterygium is extremely complex, and its complete pathogenesis remains to be further explored.
In recent years, non-coding RNAs, including long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs), have been reported to play crucial roles in transcriptional regulation of gene expression and post-transcriptional levels (Jiang et al., 2015; Gao et al., 2017). lncRNAs are autonomously transcribed non-coding RNAs greater than 200 nt in length (Wilusz et al., 2009). circRNAs are widespread endogenous non-coding RNAs that forms covalently bonded closed-loop structures due to the absence of 5' and 3' ends (Hao et al., 2020). Non-coding RNAs are involved in several biological processes, including cell proliferation, differentiation, and apoptosis, and their abnormal expression is directly related to the occurrence of many diseases, such as cancer and retinitis pigmentosa (Beermann et al., 2016; Donato et al., 2020c). Multiple lncRNAs and circRNAs have been found to be differentially expressed (DE) in pterygium compared with conjunctival tissues via microarray analysis (Liu et al., 2016; Li et al., 2018a). Thus far, the data available on pterygium are insufficient, and the complete profile and roles of coding and non-coding RNAs in the pathogenesis of pterygium remain unknown. In this study, we compared the expression patterns of coding and non-coding RNAs in pterygium tissues and paired normal conjunctival tissues by RNA-seq to clarify the pathogenesis of pterygium, and to provide new clues for the development of strategies for its prevention and treatment, and to reduce postoperative recurrence.
Materials and Methods
Sample Preparation From Patients
The samples used in this study were from six patients with primary pterygium from the Shenzhen Eye Hospital. All patients were in good health, aged between 55 and 70 years, and had undergone pterygium excision combined with conjunctival autograft. Pterygium tissue was collected from the nasal limbus, and a small piece of conjunctival tissue was collected as a control. Pterygium tissues and paired normal conjunctival tissues collected from three patients were placed in normal saline to wash the residual blood on the surface, and then quickly placed in liquid nitrogen for rapid freezing, and finally stored at −80°C for total RNA extraction. In addition, samples from the other three patients were fixed in 4% paraformaldehyde (PFA) and stored at 4°C for histological examination. This study was approved by the Ethics Committee of the Shenzhen Eye Hospital. All patients signed an informed consent form and voluntarily participated in the study.
Three pairs of paraffin-embedded pterygium tissues and normal conjunctival tissues were prepared for immunohistochemistry (IHC) studies according to the manufacturer’s instructions of the SignalStain®DAB Substrate Kit (Cell Signaling Technology, Beverly, MA, United States). The primary antibody used was rabbit anti-human CD31 (Abcam, Cambridge, United Kingdom) and the secondary antibody used was Alexa fluor 488-conjugated goat anti-rabbit IgG (Thermo Fisher Scientific, Waltham, MA, United States). The percentage of CD31-positive cells was calculated using the Image Proplus 6.0 software.
RNA Extraction, Library Construction, and RNA-Seq
Total RNAs were extracted from three paired pterygiums and normal conjunctivas using TRIzol reagent, respectively (Invitrogen, Carlsbad, CA, United States). RNA degradation and contamination were monitored on 1% agarose gels, and RNA purity was checked using the NanoPhotometer® spectrophotometer (Implen, Westlake Village, CA, United States). Additionally, RNA concentration and integrity were determined using a Qubit® RNA Assay Kit with a Qubit® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, United States) and an RNA Nano 6,000 Assay Kit with the Bioanalyzer 2,100 system (Agilent Technologies, Santa Clara, CA, United States), respectively. Ribosomal RNAs from each sample were removed using Ribo-Zero Gold Kits (Epicenter, Madison, WI, United States). Six chain-specific libraries were constructed using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ispawich, MA, United States) following the manufacturer’s recommendations, and Qubit 2.0 and quantitative real-time PCR (qRT-PCR) were used to check the quality of the library. Finally, RNA-seq was performed using an Illumina NovaSeq 6,000 instrument (Illumina, San Diego, CA, United States). The raw data obtained had been submitted to the NCBI database (accession number: PRJNA669964).
Reference Genome Mapping and Transcriptome Assembly
To obtain clean reads, low quality reads and reads containing adapter contamination or ploy-N, were removed. Simultaneously, the Q20, Q30, and GC contents of the clean data were calculated to assess their quality. All downstream analyses were based on high-quality clean data. The software HiSAT2 was used to map the sequence data to the human reference genome (GRCh38; Kim et al., 2015).
Identification of lncRNAs and circRNAs
Based on the mapped results, transcriptome assembly and reconstruction were performed using StringTie software (Pertea et al., 2015). Candidate lncRNAs were identified by the following workflow: (1) Transcripts with exons ≥2 and lengths >200 bp were obtained; (2) the known transcripts including mRNAs and other non-coding RNAs were filtered by comparing with annotation files using Cuffcompare (Trapnell et al., 2010); and (3) coding potential analysis were performed by coding potential calculator (CPC), coding-non-coding index (CNCI), and Pfam (Sonnhammer et al., 1998; Kong et al., 2007; Sun et al., 2013). Finally, the remaining transcripts with no coding potential were defined as candidate novel lncRNAs. To improve the accuracy of circRNA identification, a conjoint analysis of the two software, CIRI and find_circ, was used with default parameters (Sekar et al., 2019).
Principal Component Analysis
Principal component analysis (PCA) is performed to obtain an overview of the expression profile of mRNAs, lncRNAs, and circRNAs using OmicShare tools.1
Analysis of Differentially Expressed mRNAs, lncRNAs, and circRNAs
The fragments per kilobase of transcript per million reads mapped (FPKM) value was used to determine the expression levels of mRNAs and lncRNAs, while the expression levels of circRNAs in each sample were calculated using transcripts per kilobase million mapped (TPM; Mandric et al., 2017). Significance analysis of the two groups was performed using the DEseq2 R package (Love et al., 2014), and an FDR-adjusted p < 0.05 and log2 fold change > 1.0 were used to screen for DE genes.
Functional Analysis of mRNAs, lncRNAs, and circRNAs
The function of the DE lncRNAs was predicted according to the co-location and co-expression correlation of lncRNAs and mRNAs. For determining the mechanism of cis regulation, we set the co-location threshold to 100 kb upstream and downstream of lncRNA. For determining the mechanism of trans-regulation, Pearson’s correlation coefficient ≥ 0.9 was defined as the target correlation. The function of DE circRNAs was revealed via functional analysis of their parental genes. Gene ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using the clusterProfiler R package (Yu et al., 2012). The GO terms and KEGG pathways with p < 0.05 were indicated to be significantly enriched, and the top 15 enriched terms are shown in the figures.
Protein-Protein Interaction Network Construction
To investigate the interactive relationships of DE mRNAs, a protein-protein interaction (PPI) network was constructed using the STRING database2 with a combined score >0.4 as the cutoff criterion. PPI network visualization was carried out using Cytoscape software (version 3.5.1).
Co-expression Network Construction
A co-expression network of DE lncRNAs with their target-associated DE mRNAs was constructed using Cytoscape software (version 3.5.1) to explore the function of key lncRNAs.
Gene Expression Validation by Quantitative Real-Time PCR
To verify the sequence of circRNAs, PCR reactions were performed using Takara Ex Taq polymerase (Takara, Dalian, China) following the manufacturer’s instructions. PCR products of the circRNAs were detected by conducting 2% agarose gel electrophoresis, and Sanger sequencing was performed at TSINGKE Biotechnology Co., Ltd. (Beijing, China). In the RNase R digestion assay, the RNAs extracted from pterygium and normal conjunctiva were divided into two parts: one part was incubated with RNase R (Geneseed, Guangzhou, China) for 30 min, and the expression levels of linear RNA and circular RNA were detected by qRT-PCR; for the other part, the RNAs were directly used for the detection of the expression of linear RNA and circular RNA by qRT-PCR (control group, Mock).
Gene expression levels were verified using qRT-PCR. First-strand cDNA was synthesized from three pairs of total RNA samples from pterygium and conjunctive groups using the PrimeScript™ Master mix (Takara), following the manufacturer’s instructions. qRT-PCR was performed using SYBR Premix Ex Taq™ (Takara) and the Step-OnePlus™ Real-time PCR System (Applied Biosystems, Foster City, CA, United States). The reaction conditions were as follows: 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. Each qRT-PCR assay was performed in triplicate, and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene was used as the normalization control. All primers used in the qRT-PCR are shown in Supplementary Table S1. The relative gene expression level was calculated using the 2−ΔΔCT method.
Experiments were performed independently at least three times, and each independent test was conducted using three biological replicates. Statistical analyses of the data from CD31 immunohistochemistry and qRT-PCR were performed using GraphPad Prism 6.0. The Student’s t-test was used to assess the differences between pterygium and normal conjunctival groups. All data are expressed as mean ± standard deviation. Values of p < 0.05 were considered statistically significant.
Differences in Morphology Between Pterygium and Normal Conjunctiva
As shown in Figures 1A,B, compared with normal conjunctiva, pterygium was richer in blood vessels and connective tissues. The positive expression of CD31 in pterygium connective tissue was significantly higher than that in normal conjunctiva, according to CD31 immunohistochemical staining results (Figures 1C–E).
Figure 1. Morphologic characteristics of pterygium and normal conjunctiva. Morphology of normal conjunctiva (A) and pterygium (B) were observed by hematoxylin-eosin staining (×200), and the lumen has been indicated using a red arrow. The positive expression of CD31 in normal conjunctiva (C) and pterygium (D) was detected by immunohistochemistry (×200). (E) Percentage of CD31-positive cells in normal conjunctiva and pterygium were calculated. Data are represented as means ± SDs, n = 3 per group. *p < 0.05.
Summary of Raw Sequence Data
In this study, 657,824,860 clean reads were retained for the conjunctival and pterygium tissues after excluding low-quality reads. The average Q30 and GC content of the clean data were 93.60 and 48.04%, respectively, and more than 95.72% of the clean reads were mapped to the human reference genome (GRCh38; Table 1).
Expression Pattern of mRNA, lncRNA, and circRNA
After gene mapping, 19,775 mRNAs, 106,231 lncRNAs (including 2,884 novel lncRNAs), and 7,878 circRNAs (including 3,063 novel circRNAs) were identified in the conjunctiva and pterygium (Supplementary Table S2). Based on the matrix analysis by PCA, the expression profiles of all obtained mRNAs, lncRNAs, and circRNAs in this study could clearly distinguished pterygium from conjunctival tissue (Supplementary Figure S1). Furthermore, the first two components can explain 54.6, 44.3, and 44.3% variance of the mRNA, lncRNA, and circRNA dataset, respectively.
DE Profiles of mRNAs, lncRNAs, and circRNAs
The novel lncRNAs were classified into four types, with the maximum proportion being sense_overlapping lncRNAs (Figure 2A). There were three types of circRNAs, of which exon circRNAs accounted for 94.87%, followed by intron circRNAs which counted for 3.82% (Figure 2B). Compared with conjunctival tissue, 579 mRNAs, 275 lncRNAs, and 21 circRNAs were differentially expressed in pterygium (Supplementary Table S3). Among them, 452 mRNAs, 149 lncRNAs, and 11 circRNAs were upregulated, and 127 mRNAs, 126 lncRNAs, and 10 circRNAs were downregulated in pterygium (Figure 2C). Cluster analysis showed that more than 75% of the DE mRNAs were upregulated, and nearly 50% of the DE lncRNAs and circRNAs were downregulated (Figures 2D–F).
Figure 2. Gene expression characterization. (A) The type and proportion of novel long non-coding RNAs (lncRNAs) and (B) circular RNAs (circRNAs). (C) The count of differentially expressed (DE) messenger RNAs (mRNAs), lncRNAs, and circRNAs between pterygium and normal conjunctiva. Hierarchical clustering analysis of DE mRNAs (D), lncRNAs (E), and circRNAs (F) in pterygium and normal conjunctiva.
Functional Annotation of DE mRNAs, lncRNAs, and circRNAs
To identify the mechanism of pterygium formation, GO functional enrichment analysis and KEGG pathway analysis of DE genes were conducted. As shown in Figure 3A; Supplementary Table S4, the significantly enriched GO terms for upregulated mRNAs were mainly associated with extracellular matrix organization, muscle system process, and blood vessel morphogenesis, while response to cAMP, response to calcium ion, and regulation of neuron death were the main enriched terms for downregulated mRNAs (Figure 3C; Supplementary Table S4). Similarly, the upregulated lncRNAs were mainly involved in extracellular matrix organization, blood vessel morphogenesis, and homophilic cell adhesion via plasma membrane adhesion (Figure 4A; Supplementary Table S4), while the downregulated lncRNAs were mainly related to positive regulation of extrinsic apoptotic signaling pathway via death domain receptors, hormone-mediated signaling pathway, and phosphate ion homeostasis (Figure 4C; Supplementary Table S4). Furthermore, the significantly enriched GO terms for upregulated circRNAs were mainly centered on positive regulation of protein kinase B signaling, epithelial tube formation, and positive regulation of apoptotic process (Figure 5A; Supplementary Table S4), while downregulated circRNAs were mainly involved in columnar/cuboidal epithelial cell differentiation, phosphatidylinositol-mediated signaling, and angiogenesis (Figure 5C; Supplementary Table S4). KEGG pathway analysis indicated that the upregulated mRNAs and lncRNAs were mainly involved in ECM-receptor interaction, focal adhesion, and vascular smooth muscle contraction (Figures 3B, 4B; Supplementary Table S5), while downregulated mRNAs were mainly enriched in Toll-like receptor signaling pathway, MAPK signaling pathway, and melanogenesis (Figure 3D; Supplementary Table S5); downregulated lncRNAs were related to adherens junction, phenylalanine metabolism, and tyrosine metabolism (Figure 4D; Supplementary Table S5). Furthermore, upregulated circRNAs were related to axon guidance (Figure 5B; Supplementary Table S5), and downregulated circRNAs were associated with mucin type O-Glycan biosynthesis (Figure 5D; Supplementary Table S5).
Figure 3. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DE mRNAs in pterygium vs. normal conjunctiva. (A,B) The top 15 enriched biological processes and KEGG pathways for upregulated mRNAs. (C,D) The top 15 enriched biological processes and KEGG pathways for downregulated mRNAs.
Figure 4. Gene ontology and KEGG enrichment analysis of target genes of DE lncRNAs in pterygium vs. normal conjunctiva. (A,B) The top 15 enriched biological processes and KEGG pathways for target genes of upregulated lncRNAs. (C,D) The top 15 enriched biological processes and KEGG pathways for target genes of downregulated lncRNAs.
Figure 5. Gene ontology and KEGG enrichment analysis of parental genes of DE circRNAs in pterygium vs. normal conjunctiva. (A,B) The top 15 enriched biological processes and KEGG pathways for parental genes of upregulated circRNAs. (C,D) The top 15 enriched biological processes and KEGG pathways for parental genes of downregulated circRNAs.
PPI Network Analysis of DE mRNAs
To further explore the pathogenesis of pterygium and to identify the key regulatory factors, the 579 DE genes were submitted to the STRING online database to analyze their protein interactions. The results showed that 489 DE mRNAs established interactions with each other, forming 2,528 edges, with a combined score >0.4 (Supplementary Table S5). Of these, the two most significant modules were detected using MCODE in Cytoscape. Module 1 was composed of 31 nodes and 252 edges (Figure 6A), while module 2 contained 30 nodes and 142 edges (Figure 6B). Furthermore, the top 10 hub genes, including FN1, MMP2, PECAM1, VWF, ENG, CXCL12, EGF, CD34, VCAM1, and CDH5, were identified using the cytoHubba plugin in the Cytoscape software by MCC (Figure 6C).
Figure 6. Identification of key modules and hub genes. (A,B) Two key modules identified by MCODE in Cytoscape. (C) Ten hub genes identified by CytoHubba in Cytoscape. Pink- and blue-colored terms represent up and downregulated mRNAs, respectively.
Co-expression Network of mRNA-lncRNA
To identify the key lncRNAs related to pterygium formation, 10 hub mRNAs and 72 target-associated DE lncRNAs were screened to construct an mRNA-lncRNA co-expression network using Cytoscape. As shown in the network, one lncRNA was associated with multiple mRNAs. For instance, MIR4435-2HG and CDKN2B-AS1 were co-expressed with six mRNAs and three mRNAs, respectively (Figure 7; Supplementary Table S6).
Figure 7. The co-expression network of lncRNAs and mRNAs. The co-expression network of 10 hub mRNAs and their trans-acting lncRNAs. Pink- and blue-colored terms represent up and downregulated genes, respectively. Circle and triangles represent mRNAs and lncRNAs, respectively.
Validation of DE Genes
To validate the RNA-seq data, the expression levels of eight mRNAs (PECAM1, VCAM1, FN1, VIM, MMP2, CD34, EGR1, and NR4A2) related to blood vessel morphogenesis, focal adhesion, and EMT progression, and five co-expressed lncRNAs (ZFPM2-AS1, LINC00968, LINC00884, ZNF295-AS1, and MUC20-OT1), were detected using qRT-PCR (Figures 8A,D). These genes showed considerably similar expression trends based on qRT-PCR results and RNA-seq data.
Figure 8. Validation of the DE genes. (A) Validation of the expression of DE mRNAs by qRT-PCR. Data from the quantitative real-time PCR (qRT-PCR) are shown as columns and are placed on the left Y-axis, while the data from RNA-seq are shown as lines and are placed on the right Y-axis. (B) Sanger sequencing-based validation for the back-splice site sequence of circRNA. Figures with small letters (a–e) represent five circRNAs in the following order: hsa_circ_0007482, hsa_circ_0023988, hsa_circ_0004846, novel_circ_0000151, and novel_circ_0012935, respectively. (C) RNase R digestion assays were performed to evaluate the stability of circRNAs. The X-axis indicates circRNAs, and the Y-axis indicates the relative expression levels of glyceraldehyde-3-phosphate dehydrogenase (GAPDH), β-actin, and circRNAs. (D) Validation of the expression of DE lncRNAs and circRNAs between pterygium vs. normal conjunctiva by qRT-PCR. Data are represented as mean ± SD, n = 3 per group. *p < 0.05; **p < 0.01.
Five circRNAs (hsa_circ_0007482, hsa_circ_0023988, hsa_circ_0004846, novel_circ_0000151, and novel_circ_0012935) were selected to explore their splice sites according to Sanger sequencing results, and the RNase R digestion assay was used to assess the stability of circRNAs. The results of Sanger sequencing showed that the splice sites were consistent with those predicted by the software (Figure 8B). Furthermore, RNase R digestion assay results indicated that the five circRNAs were resistant to RNase R digestion compared to linear GAPDH and β-actin RNAs (Figure 8C). qRT-PCR results showed that the trends of five circRNAs expression were consistent with that of RNA-seq results (Figure 8D).
Pterygium is a common conjunctival degenerative disease, and pterygium excision with conjunctival autograft is the recommended treatment to reduce recurrence (Chu et al., 2020). Compared with the normal conjunctiva, pterygium possesses more vascular and fibrous tissues, and the morphology of these vessels suggests the existence of active angiogenesis in the upper subcutaneous connective tissue (Livezeanu et al., 2011; Engelsvold et al., 2013). Previous studies have performed the mRNA, miRNA, and lncRNA expression comparisons of pterygium and conjunctiva tissues via microarray or high-throughput sequencing technology (Liu et al., 2016, 2017a; He et al., 2019; Yue and Gao, 2019; Chen et al., 2020). However, studies on lncRNA and circRNA in pterygium development based on multiple types of RNA interaction analysis are limited. In this study, we first determined the expression profiles of lncRNA and circRNA in pterygium and conjunctiva by RNA-seq, and identified the biological processes and pathways related to the formation of pterygium. Abundant mRNAs, lncRNAs, and circRNAs that might be important in the pathogenesis of pterygium were identified through bioinformatics analysis and biochemical analysis.
Based on strict screening conditions, 579 DE mRNAs, 275 DE lncRNAs, and 21 DE circRNAs were obtained in this study. Functional enrichment analysis revealed that upregulated mRNAs were mainly related to extracellular matrix organization, blood vessel morphogenesis, and focal adhesion, while the downregulated mRNAs were mainly involved in cell death, which were consistent with the results of previous studies (He et al., 2019; Chen et al., 2020). Furthermore, the two most significant modules and 10 hub genes, including FN1, MMP2, PECAM1, VWF, ENG, CXCL12, EGF, CD34, VCAM1, and CDH5, were identified by PPI network analysis. Coincidently, several of the above mRNAs including FN1, MMP2, VWF, and CXCL12 were identified as the hub genes associated with pterygium development in previous studies (He et al., 2019; Yue and Gao, 2019; Chen et al., 2020). Moreover, we validated several of these genes using qRT-PCR, and compared with RNA-seq data, similar expression trends were observed. These results suggest that angiogenesis and the changes in the extracellular matrix are involved in the formation of pterygium, which is consistent with the morphological characteristics of pterygium (Engelsvold et al., 2013).
Angiogenesis is a physiological process in human development and an important pathological mechanism in the development of diseases such as cancer, obesity, and ocular disease (Cristofanilli et al., 2002; Nijhawans et al., 2020; Scimone et al., 2020). Neovascularization and fibrous connective tissue are the main components of the pterygium, and an imbalance occurring between pro- and anti-angiogenic factors results in abnormal vascular proliferation (Livezeanu et al., 2011; Tonthat and Di Girolamo, 2014). Compared with the conjunctiva, several key factors related to angiogenesis, such as PECAM1, CD34, VWF, and ENG, were highly expressed in the pterygium. PECAM1, an important endothelial cell marker, is mainly expressed in vascular endothelial cells (Murray et al., 2019; Tan et al., 2019). Consistent with our results, immunostaining of PECAM1 revealed much richer vascularization in pterygium (Aspiotis et al., 2007; Livezeanu et al., 2011). Additionally, CD34, VWF, and ENG, endothelial cell surface markers, and their positive microvascular counts were markedly higher in pterygium than those in conjunctival tissue (Marcovich et al., 2002; Lee et al., 2007; Zhang et al., 2011). During angiogenesis, vascular endothelial cells must migrate through the extracellular matrix, and this process is strictly controlled by the action of cell adhesion molecules (Bischoff, 1997). VCAM-1 is a member of the immunoglobulin superfamily and an important cell adhesion molecule, the positive rate of which in pterygium is higher than that observed in conjunctiva (Tekelioglu et al., 2006). Moreover, CDH5 and CDH6, two important epithelial cell adhesion factors, were highly expressed in pterygium tissues investigated in this study. These factors promote neovascularization, provide nutrition for the proliferation of pterygium, and accelerate pterygium growth.
Degradation and reconstruction of the extracellular matrix are key steps for angiogenesis, and these steps are controlled by the action of the members of the matrix metalloproteinases (MMPs) family (Sang, 1998; Reid and Dushku, 2010). Vascular endothelial cells secrete a variety of MMPs, which can degrade the extracellular matrix individually or in concert, and promote angiogenesis (Sang, 1998). In this study, multiple types of MMPs (MMP2, MMP3, MMP8, MMP9, MMP11, MMP16, and MMP28) were highly expressed in pterygium. The expression of MMP1, an important factor in corneal collagen degradation, in pterygium head fibroblasts is higher than that in the body and conjunctiva, which creates conditions for degradation of the exposed elastic layer by MMP2 (Li et al., 2001). MMP2 and MMP9 encode two different forms of gelatinase, are produced after the pterygium cells invading over Bowman’s layer, which contribute to the dissolution of Bowman’s layer (Dushku et al., 2001). Previous studies have shown that MMP2 and MMP9 were expressed in advanced-stage pterygium, but they were not detected in early-stage pterygium tissues and fibroblasts (Yang et al., 2009).
Extracellular matrix remodeling enables pterygium-mediated invasion of the cornea, while EMT is another important mechanism responsible for the promotion of the migration and invasion of pterygium epithelial cells (Oh et al., 2014; Yue and Gao, 2019). Kato et al. (2007) found that newly synthesized collagen fibrils accumulated between the epithelium and Bowman’s layer at the pterygium head. FN1, an important factor involved in cell adhesion and migration processes, was enriched in the PPI network in this study (Cai et al., 2017). Transforming growth factor-β is a master switch for myofibroblasts, which can induce the synthesis of FN1 and α-SMA in primary human pterygium fibroblasts and primary human conjunctival fibroblasts, and the expression levels of TGF-β2 and TGF-β3 were both upregulated in this study (Di Girolamo et al., 2004; Chen et al., 2019). VIM is a major cytoskeletal component of mesenchymal cells. Expression of VIM and α-SMA was detected in epithelial keratins from pterygium tissue but not in normal corneal epithelium (Kato et al., 2007). Upregulation of these EMT markers is an important mechanism of pterygium fibrosis and the main cause of the high recurrence rate of pterygium.
Previous reports have shown that multiple lncRNAs were DE in pterygium compared with normal conjunctiva, but the function of lncRNAs in the pathogenesis of pterygium has not been studied extensively (Liu et al., 2016). The potential functions of important lncRNAs were predicted based on the mRNA-lncRNA co-expression network created in this study. The results showed that these DE lncRNAs were mainly involved in angiogenesis and focal adhesion, which might play a key role in the pathogenesis of pterygium. Through the analysis of 10 hub mRNAs and their co-expressed lncRNAs, a considerable number of key lncRNAs that might be involved in the pathogenesis of pterygium were screened. MIR4435-2HG was co-expressed with FN1 and MMP2 and was highly expressed in pterygium. It can promote the proliferation and metastasis of various cancers by upregulating TGFβ-1 expression or by promoting the occurrence of EMT (Yang et al., 2019; Zhang et al., 2020). Moreover, LINC00968, which is co-expressed with ENG and VCAM1, can promote the proliferation and migration of endothelial cells and accelerate the proliferation and fibrosis in diabetic nephropathy (Li et al., 2018b; Wang et al., 2018). However, another study showed that LINC00968 could inhibit breast cancer cell proliferation, migration, and angiogenesis (Sun et al., 2019). Additionally, ST3GAL6-AS1 was co-expressed with a variety of mRNAs, including PECAM1, VCAM1, CXCL12, CD34, CDH5, and VWF, and was associated with colorectal cancer progression (Hu et al., 2019). In-depth functional analysis of these lncRNAs would enable the identification of key genes regulating the formation of pterygium, which may provide new insights for the development of strategies for the treatment of pterygium.
In addition to lncRNAs, circRNAs play an important role in the pathogenesis of pterygium (Li et al., 2018a). By GO functional analysis of parental genes of circRNAs, we first revealed that the upregulated circRNAs were mainly related to protein kinase B signaling, epithelial tube formation, and positive regulation of apoptotic process, while downregulated circRNAs were associated with epithelial cell development, angiogenesis, and positive regulation of mesenchymal cell proliferation. Interestingly, few DE circRNAs have been reported to be involved in tumor progression and proliferative diseases. For instance, hsa_circ_0007482 expression was upregulated in pterygium in this study, and its expression level in keloid tissue was upregulated compared with that in normal skin tissue (Shi et al., 2019). Furthermore, hsa_circ_001730 expression was upregulated in pterygium and was reported to promote the proliferation and invasion of glioma cells (Lu et al., 2019). These DE circRNAs may be involved in the regulation of abnormal proliferation and postoperative recurrence of pterygium, and it is necessary to further explore their molecular function.
Furthermore, several studies have shown that UV radiation and ion channel dysfunction are also involved in the pathogenesis of pterygium (Cimpean et al., 2013; Garreis et al., 2016; Liu et al., 2017a). Recent reports have revealed that DNA damage and mutations in ion channel-related genes are also closely related to the occurrence and progression of certain ophthalmic diseases, including retinal degeneration and retinal dystrophy (Donato et al., 2020a,b). In this study, the genes related to response to oxidative stress, including FOS, JUN, and DUSP1, were significantly downregulated, and the genes involved in ion transmembrane transport, such as CALM1, CALM2, and FXYD4, were also significantly downregulated. Therefore, oxidative stress injury caused by UV exposure and ion channel dysfunction resulting in abnormal cell proliferation and apoptosis, which may also be important factors in the occurrence of pterygium. In general, pterygium is a multifactorial disease, and further studies are warranted to elucidate its pathogenesis in the future.
In this study, we revealed that several mRNAs, lncRNAs, and circRNAs might participate in pterygium formation by regulating cell adhesion, EMT, and angiogenesis. Through the construction of the PPI network and mRNA-lncRNA co-expression networks, we screened for important mRNAs and lncRNAs that might be related to the pathogenesis of pterygium. This study provides a comprehensive differential expression profile of different types of RNA in pterygium, which helps us to explore the molecular pathomechanism of pterygium.
Data Availability Statement
RNA-seq data were submitted to NCBI SRA (accession number: PRJNA669964).
The studies involving human participants were reviewed and approved by the Ethics Committee of the Shenzhen Eye Hospital. The patients/participants provided their written informed consent to participate in this study.
XL, JZ, XhL, and JW designed the experiment and wrote and modified the manuscript. XL, DN, and JT performed the experiment and bioinformatics analysis. KZ, HH, LS, and LP collected the samples. All authors contributed to the article and approved the submitted version.
This study was supported by the Shenzhen San Ming Project (SZSM201812091).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.646550/full#supplementary-material
Supplementary Figure S1 | Principal component analysis (PCA) of mRNA, lncRNA, and circRNA expression profiles. PCA of mRNA (A), lncRNA (B), and circRNA expression profiles (C) identified based on the analysis of six samples with three biological replicates.
Aspiotis, M., Tsanou, E., Gorezis, S., Ioachim, E., Skyrlas, A., Stefaniotou, M., et al. (2007). Angiogenesis in pterygium: study of microvessel density, vascular endothelial growth factor, and thrombospondin-1. Eye 21, 1095–1101. doi: 10.1038/sj.eye.6702495
Beermann, J., Piccoli, M. T., Viereck, J., and Thum, T. (2016). Non-coding RNAs in development and disease: background, mechanisms, and therapeutic approaches. Physiol. Rev. 96, 1297–1325. doi: 10.1152/physrev.00041.2015
Cai, X., Liu, C., Zhang, T. N., Zhu, Y. W., Dong, X., and Xue, P. (2017). Down-regulation of FN1 inhibits colorectal carcinogenesis by suppressing proliferation, migration and invasion. J. Cell. Biochem. 119, 4717–4728. doi: 10.1002/jcb.26651
Chen, K., Lai, K., Zhang, X., Qin, Z., Fu, Q., Luo, C., et al. (2019). Bromfenac inhibits TGF-β1-induced fibrotic effects in human pterygium and conjunctival fibroblasts. Invest. Ophthalmol. Vis. Sci. 60, 1156–1164. doi: 10.1167/iovs.18-24743
Chui, J., Di Girolamo, N., Wakefield, D., and Coroneo, M. T. (2008). The pathogenesis of pterygium: current concepts and their therapeutic implications. Ocul. Surf. 6, 24–43. doi: 10.1016/S1542-0124(12)70103-9
Di Girolamo, N., Chui, J., Coroneo, M. T., and Wakefield, D. (2004). Pathogenesis of pterygia: role of cytokines, growth factors, and matrix metalloproteinases. Prog. Retin. Eye Res. 23, 195–228. doi: 10.1016/j.preteyeres.2004.02.002
Donato, L., Scimone, C., Alibrandi, S., Abdalla, E. M., and Sidoti, A. (2020a). New omics-derived perspectives on retinal dystrophies: could ion channels-encoding or related genes act as modifier of pathological phenotype? Int. J. Mol. Sci. 22:70. doi: 10.3390/ijms22010070
Donato, L., Scimone, C., Alibrandi, S., Pitruzzella, A., Scalia, F., D'Angelo, R., et al. (2020b). Possible A2E mutagenic effects on RPE mitochondrial DNA from innovative RNA-seq bioinformatics pipeline. Antioxidants 9:1158. doi: 10.3390/antiox9111158
Donato, L., Scimone, C., Alibrandi, S., Rinaldi, C., Sidoti, A., and D'Angelo, R. (2020c). Transcriptome analyses of lncRNAs in A2E-stressed retinal epithelial cells unveil advanced links between metabolic impairments related to oxidative stress and retinitis pigmentosa. Antioxidants 9:318. doi: 10.3390/antiox9040318
Dushku, N., John, M. K., Schultz, G. S., and Reid, T. W. (2001). Pterygia pathogenesis: corneal invasion by matrix metalloproteinase expressing altered limbal epithelial basal cells. Arch. Ophthalmol. 119, 695–706. doi: 10.1001/archopht.119.5.695
Engelsvold, D. H., Utheim, T. P., Olstad, O. K., Gonzalez, P., and Raeder, S. (2013). MiRNA and mRNA expression profiling identifies members of the miR-200 family as potential regulators of epithelial-mesenchymal transition in pterygium. Exp. Eye Res. 115, 189–198. doi: 10.1016/j.exer.2013.07.003
Garreis, F., Schröder, A., Reinach, P. S., Zoll, S., Khajavi, N., Dhandapani, P., et al. (2016). Upregulation of transient receptor potential vanilloid Type-1 channel activity and Ca2+ influx dysfunction in human pterygial cells. Invest. Ophthalmol. Vis. Sci. 57, 2564–2577. doi: 10.1167/iovs.16-19170
Hao, Z., Zhou, H., Hickford, J. G. H., Gong, H., Wang, J., Hu, J., et al. (2020). Identification and characterization of circular RNA in lactating mammary glands from two breeds of sheep with different milk production profiles using RNA-Seq. Genomics 112, 2186–2193. doi: 10.1016/j.ygeno.2019.12.014
He, S., Sun, H., Huang, Y., Dong, S., Qiao, C., Zhang, S., et al. (2019). Identification and interaction analysis of significant genes and MicroRNAs in pterygium. Biomed. Res. Int. 2019:2767512. doi: 10.1155/2019/2767512
Hou, A., Lan, W., Law, K. P., Khoo, S. C., Tin, M. Q., Lim, Y. P., et al. (2014). Evaluation of global differential gene and protein expression in primary pterygium: S100A8 and S100A9 as possible drivers of a signaling network. PLoS One 9:e97402. doi: 10.1371/journal.pone.0097402
Hu, J., Shan, Y., Ma, J., Pan, Y., Zhou, H., Jiang, L., et al. (2019). LncRNA ST3Gal6-AS1/ST3Gal6 axis mediates colorectal cancer progression by regulating α-2, 3 sialylation via PI3K/Akt signaling. Int. J. Cancer 145, 450–460. doi: 10.1002/ijc.32103
Jiang, Q., Wang, J., Wu, X., Ma, R., Zhang, T., Jin, S., et al. (2015). LncRNA2Target: a database for differentially expressed genes after lncRNA knockdown or overexpression. Nucleic Acids Res. 43, D193–D196. doi: 10.1093/nar/gku1173
Kato, N., Shimmura, S., Kawakita, T., Miyashita, H., Ogawa, Y., Yoshida, S., et al. (2007). Beta-catenin activation and epithelial-mesenchymal transition in the pathogenesis of pterygium. Invest. Ophthalmol. Vis. Sci. 48, 1511–1517. doi: 10.1167/iovs.06-1060
Kong, L., Zhang, Y., Ye, Z. Q., Liu, X. Q., Zhao, S. Q., Wei, L., et al. (2007). CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–W349. doi: 10.1093/nar/gkm391
Li, X. -M., Ge, H. -M., Yao, J., Zhou, Y. -F., Yao, M. -D., Liu, C., et al. (2018a). Genome-wide identification of circular RNAs as a novel class of putative biomarkers for an ocular surface disease. Cell. Physiol. Biochem. 47, 1630–1642. doi: 10.1159/000490982
Li, D. Q., Lee, S. B., Gunja-Smith, Z., Liu, Y. Q., and Tseng, S. C. G. (2001). Overexpression of collagenase (MMP-1) and stromelysin (MMP-3) by pterygium head fibroblasts. Arch. Ophthalmol. 119, 71–80. doi: 10.1016/S0002-9394(00)00720-0
Li, Z., Yu, Z., Meng, X., and Yu, P. (2018b). LncRNA LINC00968 accelerates the proliferation and fibrosis of diabetic nephropathy by epigenetically repressing p21 via recruiting EZH2. Biochem. Biophys. Res. Commun. 504, 499–504. doi: 10.1016/j.bbrc.2018.08.048
Liu, J., Ding, X., Yuan, L., and Zhang, X. (2016). Identification of pterygium-related long non-coding RNAs and expression profiling by microarray analysis. Int. J. Mol. Med. 38, 529–536. doi: 10.3892/ijmm.2016.2641
Lu, Y., Deng, X., Xiao, G., Zheng, X., Ma, L., and Huang, W. (2019). circ_0001730 promotes proliferation and invasion via the miR-326/Wnt7B axis in glioma cells. Epigenomics 11, 1335–1352. doi: 10.2217/epi-2019-0121
Mandric, I., Temate-Tiagueu, Y., Shcheglova, T., Al Seesi, S., Zelikovsky, A., and Mandoiu, I. I. (2017). Fast bootstrapping-based estimation of confidence intervals of expression levels and differential expression from RNA-seq data. Bioinformatics 33, 3302–3304. doi: 10.1093/bioinformatics/btx365
Marcovich, A. L., Morad, Y., Sandbank, J., Huszar, M., Rosner, M., Pollack, A., et al. (2002). Angiogenesis in pterygium: morphometric and immunohistochemical study. Curr. Eye Res. 25, 17–22. doi: 10.1076/ceyr.220.127.116.1159
Meshkani, S. E., Kooshan, N., Moghadam, A. B., Falanji, F., Adli, A., Baghbani-Arani, F., et al. (2019). Signaling roadmap to epithelial-mesenchymal transition in pterygium, TWIST1 centralized. J. Cell. Physiol. 234, 18146–18155. doi: 10.1002/jcp.28447
Murray, G. P., Post, S. R., and Post, G. R. (2019). ABO blood group is a determinant of von Willebrand factor protein levels in human pulmonary endothelial cells. J. Clin. Pathol. 73, 347–349. doi: 10.1136/jclinpath-2019-206182
Oh, S. J., Shin, J. H., Kim, T. H., Lee, H. S., Yoo, J. Y., Ahn, J. Y., et al. (2014). β-Catenin activation contributes to the pathogenesis of adenomyosis through epithelial-mesenchymal transition. J. Pathol. 231, 210–222. doi: 10.1002/path.4224
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122
Scimone, C., Alibrandi, S., Scalinci, S. Z., Trovato Battagliola, E., D'Angelo, R., Sidoti, A., et al. (2020). Expression of pro-angiogenic markers is enhanced by blue light in human RPE cells. Antioxidants 9:1154. doi: 10.3390/antiox9111154
Shi, J., Yao, S., Chen, P., Yang, Y., Qian, M., Han, Y., et al. (2019). The integrative regulatory network of circRNA and microRNA in keloid scarring. Mol. Biol. Rep. 47, 201–209. doi: 10.1007/s11033-019-05120-y
Sonnhammer, E. L., Eddy, S. R., Birney, E., Bateman, A., and Durbin, R. (1998). Pfam: multiple sequence alignments and HMM-profiles of protein domains. Nucleic Acids Res. 26, 320–322. doi: 10.1093/nar/26.1.320
Sun, X., Huang, T., Zhang, C., Zhang, S., Wang, Y., Zhang, Q., et al. (2019). Long non-coding RNA LINC00968 reduces cell proliferation and migration and angiogenesis in breast cancer through up-regulation of PROX1 by reducing hsa-miR-423-5p. Cell Cycle 18, 1908–1924. doi: 10.1080/15384101.2019.1632641
Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41:e166. doi: 10.1093/nar/gkt646
Tan, A., Lam, Y. Y., Pacot, O., Hawley, A., and Boyd, B. J. (2019). Probing cell-nanoparticle (cubosome) interactions at the endothelial interface: do tissue dimension and flow matter? Biomater. Sci. 7, 3460–3470. doi: 10.1039/C9BM00243J
Tekelioglu, Y., Turk, A., Avunduk, A. M., and Yulug, E. (2006). Flow cytometrical analysis of adhesion molecules, T-lymphocyte subpopulations and inflammatory markers in pterygium. Ophthalmologica 220, 372–378. doi: 10.1159/000095863
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621
Wang, X., Zhao, Z., Zhang, W., and Wang, Y. (2018). Long noncoding RNA LINC00968 promotes endothelial cell proliferation and migration via regulating miR-9-3p expression. J. Cell. Biochem. doi: 10.1002/jcb.28103 [Epub ahead of print]
Yang, M., He, X., Huang, X., Wang, J., He, Y., and Wei, L. (2019). LncRNA MIR4435-2HG-mediated upregulation of TGF-β1 promotes migration and proliferation of nonsmall cell lung cancer cells. Environ. Toxicol. 35, 582–590. doi: 10.1002/tox.22893
Yang, S. F., Lin, C. Y., Yang, P. Y., Chao, S. C., Ye, Y. Z., and Hu, D. N. (2009). Increased expression of gelatinase (MMP-2 and MMP-9) in pterygia and pterygium fibroblasts with disease progression and activation of protein kinase C. Invest. Ophthalmol. Vis. Sci. 50, 4588–4596. doi: 10.1167/iovs.08-3147
Zhang, S., Li, C., Liu, J., Geng, F., Shi, X., Li, Q., et al. (2020). Fusobacterium nucleatum promotes epithelial-Mesenchymal transition through regulation of the lncRNA MIR4435-2HG/miR-296-5p/Akt2/SNAI1 signaling pathway. FEBS J. 287, 4032–4047. doi: 10.1111/febs.15233
Zhang, J., Zhang, M., Li, X., Zheng, T., Mu, G., Liu, W., et al. (2011). Correlation of vascular endothelial growth factor and CD105-microvascular density in primary pterygium. J. Huazhong Univ. Sci. Technol. Med. Sci. 31, 560–564. doi: 10.1007/s11596-011-0490-4
Keywords: messenger RNA, long non-coding RNA, circular RNA, RNA-seq, pterygium
Citation: Liu X, Zhang J, Nie D, Zeng K, Hu H, Tie J, Sun L, Peng L, Liu X and Wang J (2021) Comparative Transcriptomic Analysis to Identify the Important Coding and Non-coding RNAs Involved in the Pathogenesis of Pterygium. Front. Genet. 12:646550. doi: 10.3389/fgene.2021.646550
Edited by:Peter Igaz, Semmelweis University, Hungary
Reviewed by:Luigi Donato, University of Messina, Italy
Shardul Kulkarni, Pennsylvania State University (PSU), United States
Copyright © 2021 Liu, Zhang, Nie, Zeng, Hu, Tie, Sun, Peng, Liu and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
†These authors have contributed equally to this work