Pterygium is a common ocular surface condition frequently associated with irritative symptoms. The precise identity of its critical triggers as well as the hierarchical relationship between all the elements involved in the pathogenesis of this disease are not yet elucidated. Meta-analysis of gene expression studies represents a novel strategy capable of identifying key pathogenic mediators and therapeutic targets in complex diseases. Samples from nine patients were collected during surgery after photo documentation and clinical characterization of pterygia. Gene expression experiments were performed using Human Clariom D Assay gene chip. Differential gene expression analysis between active and atrophic pterygia was performed using limma package after adjusting variables by age. In addition, a meta-analysis was performed including recent gene expression studies available at the Gene Expression Omnibus public repository. Two databases including samples from adults with pterygium and controls fulfilled our inclusion criteria. Meta-analysis was performed using the Rank Production algorithm of the RankProd package. Gene set analysis was performed using ClueGO and the transcription factor regulatory network prediction was performed using appropriate bioinformatics tools. Finally, miRNA-mRNA regulatory network was reconstructed using up-regulated genes identified in the gene set analysis from the meta-analysis and their interacting miRNAs from the Brazilian cohort expression data. The meta-analysis identified 154 up-regulated and 58 down-regulated genes. A gene set analysis with the top up-regulated genes evidenced an overrepresentation of pathways associated with remodeling of extracellular matrix. Other pathways represented in the network included formation of cornified envelopes and unsaturated fatty acid metabolic processes. The miRNA-mRNA target prediction network, also reconstructed based on the set of up-regulated genes presented in the gene ontology and biological pathways network, showed that 17 target genes were negatively correlated with their interacting miRNAs from the Brazilian cohort expression data. Once again, the main identified cluster involved extracellular matrix remodeling mechanisms, while the second cluster involved formation of cornified envelope, establishment of skin barrier and unsaturated fatty acid metabolic process. Differential expression comparing active pterygium with atrophic pterygium using data generated from the Brazilian cohort identified differentially expressed genes between the two forms of presentation of this condition. Our results reveal differentially expressed genes not only in pterygium, but also in active pterygium when compared to the atrophic ones. New insights in relation to pterygium's pathophysiology are suggested.
These authors contributed equally: Juliana Albano de Guimarães and Bidossessi Wilfried Hounpke.
Pterygium is a common disease of the ocular surface, consisting of a fibrovascular tissue that arises from the conjunctiva and extends to the cornea. It is related to irritative symptoms, poor aesthetics, astigmatism and, sometimes, worsening visual acuity, being a condition frequently associated with decreased quality of life. The prevalence of pterygium varies according to the region studied, being higher in tropical regions, where the prevalence is approximately 22%, and lower in countries outside these areas, where it is close to 2%[
In view of the controversies on the multiple mechanisms of pterygium formation, studies were carried out using high-throughput genomic technologies, such as microarrays, to attempt to elucidate mRNA expression patterns that could underlie specific molecular events of interest in its pathogenesis[
Over the past decade, the availability of large public genomic datasets and robust bioinformatics tools enabled the development of new data mining approaches in biomedical research[
The aims of this study were to evaluate gene expression in pterygium samples, perform a meta-analysis with two databases available at the Gene Expression Omnibus public repository and then, correlate the findings through a miRNA-mRNA regulatory network.
The present study was carried out with the approval of the Research Ethics Committee of the University of Campinas (UNICAMP) and was conducted in accordance with the tenets of the Declaration of Helsinki and current legislation on clinical research. Written informed consent was obtained from all subjects after explanation of the procedures and study requirements.
Participants were selected by convenience, at the Ocular Surface outpatient clinic of the Hospital of University of Campinas. Participants that presented nasal pterygia, not previously operated, and wanted surgery were included. Individuals with conjunctival surface diseases such as previous chemical burns, pemphigoid and cicatricial diseases, or suspected ocular surface neoplasia were excluded.
The pterygium was classified according to the fibrovascular tissue extension towards the cornea (grade 1 when the lesion reaches the limbus, grade 2 when it covers the cornea at about 2 mm, grade 3 when it reaches the pupil margin and grade 4 when it exceeds the pupil) and according to its biomicroscopic aspect (atrophic when the visualization of structures immediately below is possible, or active when fibrovascular tissue prevents proper visualization of underneath structures). The latter classification is based on the classification proposed by Tan et al., in which the pterygia are divided according to tissue translucency. Grade T1 (atrophic) denotes a pterygium in which episcleral vessels underlying the body of the pterygium are unobscured and clearly distinguished. Grade T3 (active or fleshy) denotes a thick pterygium in which episcleral vessels underlying the body of the pterygium are totally obscured. All other pterygia that do not fall into these 2 categories are classified as Grade 2 (intermediate). An increased fleshiness or thickness of the fibrovascular component of the pterygium is believed to be associated with recurrence after excision[
Participants underwent pterygium excision followed by an upper bulbar conjunctival free autograft placed over the site of the original lesion. All pterygia specimens were nasally located. The entire pterygium tissue was collected, and the samples were immediately stored in RNAlater stabilization solution (Invitrogen) until RNA extraction.
Tissue samples were transferred to liquid nitrogen and ground into powder. Total RNA extraction was performed using TRIzol Reagent (Invitrogen) and purified with RNeasy RNA extraction kit (Qiagen) according to the manufacturer's protocol. RNA concentration and purity were assessed through NanoDrop 2000 spectrophotometer (Thermo Scientific) whereas RNA integrity was measured in the Bioanalyzer 2100 System (Agilent Technologies). Samples with an RNA integrity number ≥ 7 were subjected to microarray analysis.
Fifty nanograms of total RNA were used to generate ss-cDNA, which was hybridized on Human Clariom D Assay gene chip (Affymetrix). This array interrogates more than 540,000 human transcripts. Sample labeling, hybridization, washing and scanning were performed according to manufacturer's protocols at Molecular Core Facility (São Paulo, SP, Brazil).
The normalization of the transcriptomic data of 9 samples from our cohort was performed by SST-RMA summarization to generate gene level expression signals using Transcriptome Analysis Console (TAC 4.0, Applied Biosystems) software after the evaluation of quality control metrics. Affymetrix probes were annotated to human genome hg38 version using TAC software and the Clariom D annotation file Clariom_D_Human.r1.na36.hg38.a1.transcript.csv. Probes mapped to more than one gene were removed and logarithm transformed expression data was exported for further analysis[
This normalized expression data of the Brazilian cohort was later used for extracting the expression intensity of miRNAs and their respective targets to reconstruct the miRNA-mRNA regulatory network, that was based in data from the meta-analysis (described later).
To identify the divergent gene expression pattern between active and atrophic pterygium, a differential gene expression analysis was performed using a linear model fitted in limma package[
Microarray raw data were pre-processed using the Robust Multichip Average (RMA) method[
Gene expression datasets from microarray studies including human patients that underwent pterygium excision were searched in the Gene Expression Omnibus (GEO) public repository, by December 2020. Search was conducted using the term "pterygium". Datasets were only included if the studies analyze both affected patients and healthy tissues. All included datasets were from studies published in peer-reviewed journals. In addition, studies were excluded if samples were submitted to cell culture.
The meta-analysis was performed with RankProd package[
To facilitate the interpretation of the biological significance of the up-regulated genes identified by the meta-analysis, a functional gene set analysis (GSA) was performed using ClueGO[
To gain more insights into the biological processes associated with the down-regulated signature, these genes were used for an additional functional analysis based on the Functional Analysis of Individual Microarray/RNAseq Expression (FAIME) algorithm implemented in seq2pathway package[
The list of targeted genes included in this analysis was derived from the functional gene set analysis performed with the up-regulated genes identified by the meta-analysis. Briefly, we extracted the list of genes associated with the functional gene ontology terms and biological pathways predicted by ClueGO. The list of conserved miRNAs predicted to target these genes was identified in TargetScanHuman[
Nine participants were included, being 6 males (aged 29–70 years old) and 3 females (aged 27–48 years old). The pterygia were classified as grade 1 in 1 case, grade 2 in 4 cases, grade 3 in 3 cases and grade 4 in 1 case. According to their biomicroscopic aspect, there were 2 atrophic and 7 active pterygia. Demographic information of the participants is presented in Table 1.
Table 1 Demographic information of the participants.
Participants Age Gender Grade Biomicroscopic aspect 1 45 Female 2 Active 2 53 Male 2 Atrophic 3 34 Male 1 Active 4 70 Male 3 Atrophic 5 27 Female 2 Active 6 48 Female 3 Active 7 35 Male 2 Active 8 51 Male 4 Active 9 29 Male 3 Active
We evaluated the differential expression comparing active pterygium with atrophic pterygium using data generated from the Brazilian cohort. This analysis identified 219 up-regulated and 92 down-regulated genes. The pattern of gene expression generated with top 30 DE genes is presented as heatmap in Fig. 1. The top 10 up- and down-regulated genes and their main biological processes are presented in Table 2.
Graph: Figure 1 Heatmap plot of top 30 differentially expressed genes from the comparison of active pterygium with atrophic pterygium filtered by the fold change. The pattern of gene expression shows two different unsupervised clusters indicating a divergent expression pattern between these two clinical aspects. High expression is indicated in red and low expression is indicated in blue.
Table 2 List of top 10 up- and down-regulated genes in active vs atrophic pterygium and their associated biological processes.
DE Fold change AUC Main biological process 126.36 0.0003 0.86 Immune and inflammatory response 8.76 0.0065 0.86 Aldosterone synthesis and secretion 7.29 0.0002 1.00 Angiogenesis 7.04 0.0006 1.00 Degradation of the extracellular matrix 6.95 0.0006 1.00 (S)-reticuline biosynthesis and Tyrosine metabolism 6.78 0.0075 0.93 Cell growth. differentiation and apoptosis 6.28 0.0016 1.00 Modulation of intracellular signaling pathways 5.64 0.0007 1.00 NF-kappaB Signaling 5.44 0.0043 0.71 Targets extracellular lipids; anti-inflammatory and immunosuppressive functions 5.35 0.0002 1.00 Keratinization 0.03 0.0032 1.00 Circadian rhythm 0.16 0.0001 1.00 Transporter activity and glycerol channel activity 0.17 0.0091 1.00 * 0.18 0.0005 1.00 # 0.18 0.0027 1.00 Chemotactic and angiogenic properties 0.19 0.0083 1.00 Circadian rhythm 0.21 0.0024 1.00 Cell surface interactions at the vascular wall 0.22 0.0025 1.00 Integrin Pathway and Collagen chain trimerization 0.23 0.0009 1.00 # 0.23 0.0004 1.00 Ubiquitination and proteasome-dependent degradation
*Long non coding RNA; # miRNA.
Three studies (GSE2513, GSE51995 and GSE83627) fulfilled the inclusion and exclusion criteria described in methods section of which two datasets (GSE51995 and GSE83627) were identified as duplicated datasets by PCA analysis (Fig. 2; 8 samples in each dataset). An exploratory analysis also confirms the presence of 8 duplicates based on the comparison of the normalized expression of that datasets. Therefore, GSE83627 was removed from the meta-analysis. The studies of the final analysis include data from 8 conjunctiva samples and 12 pterygium samples.
Graph: Figure 2 Gene expression pattern of conjunctiva and pterygium from the meta-analysis. The heatmap was constructed using the top 30 differentially expressed genes (15 up- and 15 down-regulated). The dendrogram indicates two clusters stratified using hierarchical clustering. Expression pattern was rescaled using a list of human housekeeping genes from HRT Atlas database[
The meta-analysis identified 154 up-regulated and 58 down-regulated genes (Table 3). The top 10 up- and down-regulated genes are presented in Table 3. A heat map visualization of the expression pattern of the top 30 DE genes identified from the meta-analysis showed a stratification of conjunctiva and pterygium samples in two different clusters (Fig. 2).
Table 3 Top differentially expressed genes identified in the meta-analysis.
Genes Fold-change in individual studies (log2 FC) Meta-analysis results GSE2513 GSE51995 Ave log2(FC) FC FDR SPRR3 2.16 3.78 2.97 7.85 < 0.0001 SPRR1B 1.91 3.23 2.57 5.94 < 0.0001 OLFM4 0.23 4.61 2.42 5.35 < 0.0001 KERA 0.53 4.02 2.28 4.85 < 0.0001 TFF1 2.49 1.66 2.07 4.21 < 0.0001 ASPN 1.39 2.37 1.88 3.67 < 0.0001 GJA1 1.04 2.58 1.81 3.51 < 0.0001 KRT3 0.09 3.51 1.80 3.49 < 0.0001 S100P 1.15 2.44 1.79 3.47 < 0.0001 KRT16 0.74 2.76 1.75 3.36 < 0.0001 FOS − 3.84 − 3.21 − 3.53 0.09 < 0.0001 FOSB − 2.39 − 2.73 − 2.56 0.17 < 0.0001 WIF1 − 0.49 − 4.06 − 2.28 0.21 < 0.0001 NR4A2 − 2.21 − 1.68 − 1.94 0.26 < 0.0001 DUSP1 − 1.82 − 1.85 − 1.84 0.28 < 0.0001 TFPI2 − 0.84 − 2.47 − 1.65 0.32 < 0.0001 TNNI2 − 1.46 − 1.74 − 1.60 0.33 < 0.0001 CYP26A1 − 1.21 − 1.58 − 1.39 0.38 < 0.0001 IGFBP3 − 1.21 − 1.43 − 1.32 0.40 < 0.0001 IER2 − 1.12 − 1.47 − 1.29 0.41 < 0.0001 BTG2 − 1.17 − 1.38 − 1.28 0.41 < 0.0001
Genes were ranked according to the fold change. FC: Fold-change; Ave log2(FC): Base 2 logarithmic scale of average FC; FDR: False Discovery Rate. Ave log2(FC) is expressed as arithmetric mean log2(FC) across studies.
To gain further insights into the biological mechanisms associated with the expression pattern observed in the meta-analysis, a gene set analysis was performed with ClueGO using the full list of the up-regulated genes. This tool clusters biological pathways and gene ontology terms that participate in the same biological function, thereby showing the top significant non-redundant biological pathways and gene ontology terms (Fig. 3). Pathways associated with the remodeling of extracellular matrix were overrepresented in this gene set analysis. The following groups were also represented in the network: serine-type endopeptidase activity, formation of cornified envelope, estrogen signaling pathway, unsaturated fatty acid metabolic process, regulation of sensory perception of pain, positive regulation of calcium ion transport into cytosol and regulation of biomineral tissue development. ClueGO also enables the visualization of gene interactions between different gene ontology and biological pathways. Based on network connectivity, MMP2, FN1, COL11A1, LAMB3, THBS2 and YAP1 were predicted as hub genes and might play prominent role in the organization of this network. Interestingly, these genes show a tight relationship between extracellular matrix remodeling terms and other pathways and ontology terms (Fig. 3).
Graph: Figure 3 Functional analysis of up and down regulated genes in pterygium. The main pathways that are altered in patients with pterygium are indicated by colored nodes. GSA terms are interconnected with their shared genes. Close related terms are indicated by more than one color. GSA: Gene Set Analysis.
To further explore the mechanism of regulation underlying the biological processes identified in the meta-analysis, we reconstructed the miRNA-mRNA target prediction network based on the set of up-regulated genes presented in the gene ontology and biological pathways network (Fig. 3). This strategy allows the identification of 17 target genes negatively correlated with their interacting miRNAs from the Brazilian cohort expression data (Fig. 4). Based on network connectivity analysis we identified two main clusters (Fig. 4A,B). Once again, the first cluster is associated with 4 enriched GSA terms involving extracellular matrix remodeling mechanisms (Fig. 4A). The second cluster involves three biological processes: (
Graph: Figure 4 miRNA-mRNA targeting networks generated from the Brazilian cohort. miRNA-mRNA interaction was predicted using TargetScanHuman database. Pairwise correlation was computed between miRNA expression and the expression of the set of up-regulated genes presented in the gene ontology network generated in Fig. 2A. Only miRNA and genes which have negative correlation (less than − 0.50) were included. Regulatory networks associated with extracellular matrix remodeling are clustered (A). The second cluster (B) involves 3 different terms including the formation of the cornified envelope, the establishment of skin barrier and unsaturated fatty acid metabolic process. Genes are represented in grey dots and dot size is proportional to the level of expression; miRNA are represented in blue and edges are labeled with the correlation level between miRNA and genes.
In order to study the mechanisms underlying the down-regulated pattern observed in the meta-analysis, we next performed a functional analysis based on FAIME algorithm that attributes at a single sample level a score to each predicted term or biological process. This analysis shows mainly GSA terms associated with metabolic and regulation pathways (Fig. 5) ranked by their discriminatory capacity to classify pterygium and control in different clusters (AUC equal or more than 0.7). Gene interactions between different gene ontology and biological pathways are shown in Fig. 5.
Graph: Figure 5 Functional analysis of the down-regulated genes identified in the meta-analysis based on FAIME algorithm. Main down-regulated pathways are indicated by colored nodes. Pathways were selected based on the FAIME scores computed using the expression of down-regulated genes (grey dots). AUC was used to evaluate the discriminatory capacity of each enriched term for a binary classification (pterygium vs. control). FAIME: Functional Analysis of Individual Microarray/RNAseq Expression; AUC: Area Under the curve.
Despite being a common disease of the ocular surface, the pterygium exact etiopathogenesis is not yet fully understood. For this reason, a meta-analysis of transcriptomic datasets deposited in public repositories was performed in order to generate new biological insights about its pathogenic mechanisms. Furthermore, gene expression data was obtained from 9 Brazilian patients and later correlated with data from the meta-analysis.
In the Brazilian cohort, 66.7% of the patients were males, 44.4% had pterygium grade 2 and 77.8% had active pterygium. Demographic and epidemiological characteristics of this cohort are in accordance with previously described features of patients with pterygium in our region, as described by Artioli-Schelini et al., that found that most of the patients had active and grade 2 pterygia[
In our meta-analysis, we found 212 DEGs, being 154 up-regulated and 58 down-regulated genes. Amongst the top up-regulated genes, SPRR3, SPRR1B, OLFM4, KERA, TFF1, ASPN, GJA1, KRT3, S100P and KRT16 were the most significantly differentially expressed genes. On the other side, amongst the top down-regulated genes, FOS, FOSB, WIF1, NR4A2, DUSP1, TFPI2, TNNI2, CYP26A1, IGFBP3, IER2 and BTG2 were the most significantly differentially expressed ones. As in previous studies, we found that genes associated with wound healing response and extracellular matrix were the most predominantly overexpressed[
A meta-analysis performed by Xu and colleagues, based on three datasets (GSE2513, GSE51995 and GSE83627) available at the Gene Expression Omnibus public repository, disclosed upregulation of other genes in addition to the ones found in our meta-analysis, including FN1, a key molecule in the extracellular matrix, and PI3, involved in apoptosis, focal adhesion and extracellular matrix-receptor interaction[
Amongst the genes correlated with miRNA expression from the Brazilian cohort, GJA1 encodes a protein called connexin 43, a component of gap junctions, COL6A3 encodes one of the chains of type VI collagen, a component of connective tissues, and THBS4 and THBS2 encode for members of the thrombospondin family, which are adhesive glycoproteins that mediate cell-to-cell and cell-to-matrix interactions. YAP1 is a gene known to play a role in the development and progression of multiple cancers as a transcriptional regulator of a signaling pathway. LOXL1 encodes a protein essential to the biogenesis of connective tissue, while LHFPL2 is a member of the superfamily of tetraspan transmembrane protein encoding genes.
Also, regarding the genes correlated with miRNA expression obtained from the Brazilian cohort, pathways involved with unsaturated fatty acid metabolic processes were found to be differentially expressed in pterygium. This supports other studies findings in which pterygium fibroblasts were found to have an increased activity of intracellular cholesterol metabolism and high expression of LDL receptors[
Also reinforcing previous studies findings, increased cellular proliferation and reduced apoptosis appear to be involved in pterygium development[
OLFM4 consists of an antiapoptotic factor that promotes tumor growth and also an extracellular matrix glycoprotein that facilitates cell adhesion. In pterygium's pathogenesis, it could be associated an uncontrolled tissue growth. GJA1, a member of the connexin gene family, encodes the major protein of gap junctions and is related to cell signaling, apoptotic processes and chronic inflammatory processes. The FOS gene family consists of 4 members, which are FOS, FOSB, FOSL1 and FOSL2. The proteins encoded by these genes are leucine zipper proteins involved in regulation of cell proliferation, differentiation and transformation. In addition, these proteins play a role in cellular response to reactive oxygen species and in inflammatory response[
In Tong and colleagues' study, genes encoding for apoptosis (TGM2, IGFBP3 and DUSP1) and stress-inducible transcription regulator genes (ATF3, BTG2, EGR1, ERG2, FOS, JUN, NR4A1 and NR4A2) were down-regulated in pterygium, reinforcing the over-proliferative tendency in pterygium[
Finally, we could determine differences in gene expression between atrophic and active pterygium, leading to the conclusion that distinct biological pathways may be differently combined in terms of up and down-regulation to result in certain phenotypic characteristics of pterygium. In active pterygium, genes involved with immune and inflammatory response (CXCL9 and PLA2G2D), angiogenesis (SERPINB5 and BM4), keratinization (DSG1), response to UV light (TYR) and negative regulation of apoptotic process (PIM1) are up regulated in relation to atrophic pterygium. On the other side, genes involved with negative regulation of transcription (CIART), response to osmotic stress (AQP9), chemotaxis (CXCL6), intracellular signal transduction (TREM1), extracellular matrix organization (COL6A3) and protein ubiquitination (NEURL4) are down regulated in relation to atrophic pterygium.
Limitations of the study included the Brazilian sample size, which was not large, in addition to an absence of healthy conjunctival samples to compare with pterygium samples. Consequently, these samples were not included in the analysis of up and down regulated genes of the main meta-analysis, which was based on already available datasets. Brazilian data was used for the reconstruction of the miRNA-mRNA regulatory network and for a transcriptomic analysis in which active and atrophic pterygia were compared. Even though the sample size of this analysis was limited, the top ranked genes made it possible to stratify active and atrophic pterygium in different clusters using unsupervised clustering method as shown in Fig. 1.
In conclusion, reinforcing previous studies, we identified crucial pathways in the pathophysiology of pterygium, with the main ones being related to extracellular matrix remodeling and dysregulated wound healing response. Furthermore, different pathways involved in the pathogenesis of pterygium were evidenced to be related to different phenotypic characteristics of this condition, which may also provide some clarity in our understanding of this condition's pathophysiology. Thus, we believe this study enriches our understanding of molecular mechanisms involved with pterygium, providing insights that may contribute for the future development of new therapeutic targets for its management.
J.A.G. and B.W.H. wrote the main manuscript text. J.A.G., B.D., A.L.M.B. and M.G.M.V. collected data from participants and performed the experiments. J.A.G., B.W.H., L.C.B. and M.B.M. analyzed the results of the experiments and performed statistical analyses. B.W.H. prepared figures. M.A. and M.B.M. supervised the experiments and reviewed the manuscript. All authors reviewed the manuscript.
FAPESP Grant 2014/19138-5.
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
By Juliana Albano de Guimarães; Bidossessi Wilfried Hounpke; Bruna Duarte; Ana Luiza Mylla Boso; Marina Gonçalves Monteiro Viturino; Letícia de Carvalho Baptista; Mônica Barbosa de Melo and Monica Alves
Reported by Author; Author; Author; Author; Author; Author; Author; Author