miRNAs and Genes Involved in the Interplay between Ocular Hypertension and Primary Open-Angle Glaucoma. Oxidative Stress, Inflammation, and Apoptosis Networks

Glaucoma has no cure and is a sight-threatening neurodegenerative disease affecting more than 100 million people worldwide, with primary open angle glaucoma (POAG) being the most globally prevalent glaucoma clinical type. Regulation of gene expression and gene networks, and its multifactorial pathways involved in glaucoma disease are landmarks for ophthalmic research. MicroRNAs (miRNAs/miRs) are small endogenous non-coding, single-stranded RNA molecules (18–22 nucleotides) that regulate gene expression. An analytical, observational, case-control study was performed in 42 patients of both sexes, aged 50 to 80 years, which were classified according to: (1) suffering from ocular hypertension (OHT) but no glaucomatous neurodegeneration (ND) such as the OHT group, or (2) have been diagnosed of POAG such as the POAG group. Participants were interviewed for obtaining sociodemographic and personal/familial records, clinically examined, and their tear samples were collected and frozen at 80 °C until processing for molecular-genetic assays. Tear RNA extraction, libraries construction, and next generation sequencing were performed. Here, we demonstrated, for the first time, the differential expression profiling of eight miRNAs when comparing tears from the OHT versus the POAG groups: the miR-26b-5p, miR-152-3p, miR-30e-5p, miR-125b-2-5p, miR-224-5p, miR-151a-3p, miR-1307-3p, and the miR-27a-3p. Gene information was set up from the DIANA-TarBase v7, DIANA-microT-CDS, and TargetScan v7.1 databases. To build a network of metabolic pathways, only genes appearing in at least four of the following databases: DisGeNet, GeneDistiller, MalaCards, OMIM PCAN, UniProt, and GO were considered. We propose miRNAs and their target genes/signaling pathways as candidates for a better understanding of the molecular-genetic bases of glaucoma and, in this way, to gain knowledge to achieve optimal diagnosis strategies for properly identifying HTO at higher risk of glaucoma ND. Further research is needed to validate these miRNAs to discern the potential role as biomarkers involved in oxidative stress, immune response, and apoptosis for the diagnosis and/or prognosis of OHT and the prevention of glaucoma ND.


Introduction
Glaucoma is a neurodegenerative disease and a leading cause of irreversible blindness, affecting over 60 million people worldwide [1]. The number of people with glaucoma will increase to 111.8 million by 2040 [2]. These estimates are important in guiding the designs of glaucoma screening, diagnosis and treatment, research milestones, and related public health strategies.
Primary open-angle glaucoma (POAG) is the most prevalent type of glaucoma, typically characterized by adult onset, chronic intraocular pressure (IOP) elevation, IOPdependent progressive apoptotic retinal ganglion cell (RGC) death, and visual field loss [1,2]. Main clinical features of glaucomatous optic nerve degeneration include the following: optic disc deepening, papillary hemorrhages, and specific defects of the retinal nerve fiber layer (RNFL) [3].
Diagnosis of POAG underlies a variety of clinical hallmarks such as the IOP elevation (the major risk factor) and the optic nerve head changes, as reflected by the structural/functional imaging techniques to appropriately establish the glaucoma stage [4][5][6]. Some decades ago, the term ocular hypertension (OHT) arose to catalogue the eyes with elevated IOP but not displaying optic disc damage or altered visual field. In 1977, Shaffer warned professionals about the false sense of security that may involve the managing of these patients that only have OHT, but no signs of neurodegeneration (ND) [7], suggesting that the term glaucoma suspects is preferred to ocular hypertensives. Glaucoma suspects define the group of patients with increased IOP and borderline optic discs (mildto-moderate alterations), RNFL anomalies and/or visual field changes as well as glaucoma family history and the occurrence of other POAG risk factors. These patients have a higher risk of undergoing optic nerve degeneration (OND) than the normal population [8]. Therefore, classic and emerging technologies are essential for the early detection of POAG damage, but the identification of glaucomatous pre-perimetric changes continues to be a challenging issue for ophthalmologists and researchers. As the elevated IOP is the main risk factor, current knowledge on the etiopathogenic mechanisms for HTO and POAG remains incomplete. Among the cellular and molecular processes underlying these diseases, the following have largely been considered: oxidative/nitrosative stress, mitochondrial failure, inflammation and immune response, autophagy/mitophagy, apoptosis, neurotoxicity, ND, etc. [3][4][5][6][7][8][9]. It is imperative that the above processes are individually and integrally addressed.
The regulation of gene expression and gene networks as well as the multifactorial pathways involved in glaucoma disease are high-priority landmarks for ophthalmic research. Recent data have pinpointed potentially interesting routes to mediate glaucomatous RGC dysfunction [9]. Meanwhile, hypotensive therapy (medical, laser, surgical) is the only way to fight against the elevated IOP [10,11]. Despite experimental advances in neuroprotection [12][13][14], there is no definitive cure for glaucoma ND. Despite great advances in glaucoma, there is still no reliable biomarker that can pre-clinically identify subjects at risk of POAG initiation and progression. Molecular-genetic diagnostic challenges for POAG are needed to complete the current knowledge in disease pathogenesis as well as to design new diagnostic and therapeutic strategies gathered under the recently proposed term glaucoma theranostics [15] for better eye care.
MicroRNAs (miRNAs) are a family of small endogenous non-coding, single-stranded RNA molecules (18-22 nucleotides) that regulate gene expression by either inhibiting mRNA translation or by degrading mRNA [16,17]. miRNAs are involved in the pathological processes of numerous diseases including eye pathologies [18,19]. Specific miRNAs have also been proposed to regulate IOP [20,21] as well as for use as noninvasive biomarkers for the diagnosis of glaucoma [22][23][24]. Our group has been widely utilizing tear samples for ophthalmic research with satisfactory results [25,26]. Because of this, we focused on analyzing tear samples that are relatively easy to collect, store, and process, to identify specific miRNAs (and its genetic targets) that differentially express themselves in the clinically silent interface between OHT and glaucoma, and to support their presumable interest as diagnostic biomarkers for individuals at risk of glaucoma ND.

Materials and Methods
We performed an analytical, observational, case-control study including 42 patients recruited from the ophthalmological department of the University Hospital Dr. Peset (Valencia, Spain), who agreed to participate in the study and signed the informed consent. Sample size calculation was performed using the ssize.fdr R package (R Core Team, Vienna, Austria) to detect a 1.5-fold change and achieve an 80% statistical power with a false discovery rate of 15% and an estimated proportion of non-differentially expressed miRNAs of 0.85. The study adhered the Declaration of Helsinki (Edinburgh, 2000) and the Ethics Committee standards of the study center (no. 81/16). All requirements for clinical research to maintain the privacy of the data obtained were met. Two ophthalmologists from the glaucoma section performed a systematized examination of the suitable study participants to ensure their appropriated status (Table 1), which were distributed into two groups: (1) patients diagnosed of POAG (n = 20), and (2) patients with OHT (n = 22).  [27].
Sampling was done by collecting reflex tears from the inferior meniscus of the eye without instilling anesthetics, as described in our previous work [25,26], using microhe-matocrite capillary tube, that were appropriately labeled, and immediately transferred into microcentrifuge tubes and stored in an ultra-freezer at −80 • C. On the day of processing, samples were defrosted and prepared for RNA extraction using the miRCURY RNA Isolation Kit-Biofluids (EXIQON Inc., Woburn, MA, USA). This kit is designed to isolate all RNAs sized less than 1000 nucleotides, from mRNA and tRNA to microRNA and small interfering RNA. We carried out the RNA extraction according to the manufacturer's instructions. Briefly, the purification is based on spin column chromatography using a proprietary resin as the separation matrix. Small RNAs are separated from other cellular components such as proteins without the use of phenol or chloroform.
The quality and quantity of total RNA obtained from tears was assessed using a Bioanalyzer 2100 (Agilent ® Technologies, Inc., Santa Clara, CA, USA), and the RNA 6000 Nano Kit (Agilent ® Technologies, Inc.). RNA libraries were prepared using NEBNext ® Multiplex Small RNA Library Prep Set for Illumina ® (#E7300 y #7580; New England BioLabs ® , Inc., Ipswich, MA, USA), according to the manufacturer's protocol https://international.neb.com/ protocols/2018/03/27/protocol-for-use-with-nebnext-small-rna-library-prep-set-for-illu mina-e7300-e7580-e7560-e7330 (accessed on 17 May 2021). According to the guidelines for low RNA concentration samples, the adapters and RT primers were diluted 1:2 with nuclease-free water and 15 cycles were used for the amplification by PCR. The indexed libraries were purified using the QIAquick ® PCR Purification Kit (#28104, QIAGEN ® , Hilden, Germany). Library quality control was assessed using a 4200 TapeStation (Agilent ® Technologies, Inc.) and High Sensitivity D1000 Kit (Agilent ® Technologies, Inc.). The miRNA fraction of each library (120-200 bp) was collected using the Pippin Prep System (Sage Science, Inc., Beverly, MA, USA) following the manufacturer's guidelines and using 3% agarose dye free gel cassettes with internal standards (Marker P) (Sage Science # CDP3010). miRNAs were quantified using a 4200 TapeStation (Agilent ® Technologies, Inc.) and High Sensitivity D1000 Kit (Agilent ® Technologies, Inc.) prior to normalization and pooling. Sequencing was performed on a NextSeq 500 System (Illumina, Inc., San Diego, CA, USA) with a Mid-Output flow cell for 150-cycle reads obtaining about 3.5 million reads per sample. FASTQ file quality was assessed using FASTQC tool (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 21 May 2021). Adapters and low liability reads were removed. Non-coding RNAs previously described in the ENSEMBL database were selected and characterized. Statistical analyses (normalization, differential expression, and significance) were performed using Limma and edgeR packages deposited in Bioconductor (www.bioconductor.org). A predictive analysis based on receiver operating characteristic (ROC) curves was performed to select those miRNAs showing an area under the curve (AUC) greater than 0.75. Subsequently, an analysis of the main components (PCA) was performed. Later, the target genes of the selected miRNAs were determined using the R miRFA pipeline [28], which is supported by data from the DIANA-TarBase v7, DIANA-microT-CDS and TargetScan v7.1 databases. Additionally, UniProt and Gene Ontology databases have been used to search for terms associated with the trabecular meshwork functions and the glaucoma pathology [29]. To select all genes needed to build a network of metabolic pathways, we considered only those genes that appeared in at least four of the following databases: DisGeNet, GeneDistiller, MalaCards, OMIM PCAN, UniProt, and GO [30]. These genes also need to have at least seven terms in any of the UniProt categories above-mentioned. Two approaches have been used to build the networks: (1) enrichment of metabolic pathways using the g: Profiler, GSEA, Cytoscape, and EnrichmentMap tools [31], and (2) gene function prediction using the GeneMANIA gene integration tool [32], which can function both as an independent server and as an application in Cytoscape.
For statistical analysis, comparison of two categorical variables was performed using the Pearson Chi square test. We used the Shapiro-Wilk test to check the distribution quantitative variables. The comparison of two means was analyzed by means of the Student t-test for independent samples (normal variables) or the Mann-Whitney U test (non-normal variables). The comparison of more than two means was carried out by the analysis of variance (ANOVA, normal variables) or the Kruskal-Wallis test (non-normal variables). Statistical proceedings were performed using the statistical package ((BM SPSS Statistics for Windows, Version 24.0. Armonk, NY, USA: IBM Corp).

Sociodemographic and Ophthalmologic Data
Mean age of participants was 64.5 ± 1.4 years for the POAG group and 61.1 ± 2.4 years for the OHT group. Comparison between groups did not show statistically significant differences (p > 0.05).
Regarding the distribution by gender, both groups showed a greater proportion of women. No statistically significant differences were observed between groups (p > 0.05).
Sociodemographic characteristics of the study participants are listed in Table 2. No significant differences were found in any of the study variables. The OHT subjects showed elevated IOP but no damage at OCT examination, normal visual fields, and normal ocular fundus. However, POAG patients showed IOP elevation, an increase of optic disc excavation, optic nerve damage, and/or altered visual fields. Ophthalmologic parameters of the study participants are shown in Table 3.  The next step was to compare the expression levels of those 95 miRNAs in both study groups (POAG patients vs. OHT subjects). From these 95 miRNAs, 87 showed no significant differences between groups (p > 0.05). However, we found six upregulated and two downregulated miRNAs in tears from the POAG patients (Table 4). At this point, we build again the PCA plot, taking into consideration only those eight miRNAs that differentially expressed between groups. The PCA showed a better separation of the miRNAs present in tears from the OHT vs. POAG patients ( Figure 2).  In addition, a ROC curve was performed to analyze whether these eight miRNAs could predict POAG in the HTO individuals. Table 5 shows the AUC for each miRNA. The AUC of four of these miRNAs was greater than 0.75, so they could be considered as good predictors of POAG.

Bioinformatic Processing
Using the databases mentioned in the methodology section and combining their predictions, we identified a total of 14,379 potential target genes of the eight miRNAs. Then, these 14,379 genes were searched in four databases (DisGeNet, GeneDistiller, MalaCards, and OMIM PCAN) that collected information on the association between genes and pathologies, finding the presence of 390, 183, 145, and 7 genes in each of them, respectively.
Eleven categories related to trabecular meshwork functions and the glaucoma pathology were selected and the number of genes that had terms in each of them are shown in Table 6. The biological functions marked in bold represent those with the highest number of genes (extracellular matrix and apoptosis processes). In addition, genes related to neurodegeneration, inflammation, and oxidative stress, among others, were also identified. We performed the Gene Ontology (GO) Biological Process analysis and found 201 items to which the eight identified miRNAs were significantly associated. Table 7 shows the top 40 biological processes of the above miRNAs. In addition to this GO analysis, we identified 36 reactomes (REACT) significantly associated with the eight miRNAs (Table 8), among them, genes related to apoptosis and extracellular matrix conditions. Additionally, those genes involved in aqueous humor homeostasis, oxidative stress, immune response regulation (TGF-β) and neurodegeneration were identified.
Following the screening criteria detailed in the methodology section, a total of 114 genes were selected to build a network of metabolic pathways (Figure 3). This network shows the clusters involved in apoptosis (blue box), extracellular matrix and collagen concerns (red box), endopeptidase activity (green box) as well as in development pathways (upper left side) among other clusters. Network data are included as Supplementary Materials, as both CYS (Cytoscape session file) and SIF (simple interaction file) formats.
By using the GeneMania program, different maps were created. In these maps, each circle is a node representing a gene. Black nodes belong to the original set of 114 genes studied, while grey circles represent those not present as genes in the original set, but very closely related to them according to the selected network criteria. These diagrams also give information about the number of interactions of each gene (the size of the circle is proportional to this number), and the edges indicate a relationship between genes, meaning that the thicker the line, the more intense the relationship.
With these proceedings, we finally built three maps: the first is a map of genetic interactions between these genes ( Figure 4); the second is a map showing the physical interactions ( Figure 5); and the third is a map of the metabolic pathways ( Figure 6). In these figures, grey circles represent genes not present in the original 114 set, but closely related to them. Circle size is proportional to the number of interactions of each gene. Edges indicate a relationship between genes(a thicker line the more intense relationship). Network data from Figures 4-6 are also included in the Supplementary Materials as both CYS (Cytoscape session file) and SIF (simple interaction file) formats.
People more predisposed to glaucoma must be early identified to undertake measures for avoiding optic nerve irreversible damage and vision-threatening consequences [1][2][3][4][5][6]. Over the years, this concern has continued to shape the thinking of ophthalmologists and researchers worldwide because hypotensive treatment has to be promptly initiated to appropriately control the elevated IOP. In this context, biomarkers might be key diagnostic tools to identify individuals at higher risk for glaucoma such as HTO individuals, just like we did in this report. Our differential expression analysis carried out by NGS led us to identify the tear fingerprint of eight miRNAs that significantly expressed between the POAG patients and HTO individuals. Then, we comprehensively explain the above miRNAs as well as the underlying miRNA-related pathways involved in HTO and POAG, as reflected in Figures 3-6 as well as in the Supplementary Materials.
Specific miRNAs involved in IOP homeostasis and glaucomatous OND were reported by different authors [21][22][23][24]. Others have also conducted similar studies to ours in different ophthalmic disorders [33] using different biological samples [34][35][36]. Liu et al. [37] analyzed the miRNA signature in the aqueous humor of POAG patients and its role in the etiopathogenic mechanisms of glaucomatous OND, reporting 88 miRNAs that differentially expressed between POAG and cataracts. Among them, the miR-151-a-3p was identified in the aqueous humor by these authors, in a similar manner to us in the course of the present work, in which the same miRNA was downregulated in tears of our study participants. Ertekin et al. [38] found miR-26b-5p in plasma from patients with wet macular degeneration. The authors associated this miRNA with oxidative stress occurring in the retinal pigment epithelial cells of the eyes with this disease. We also found the miR-26b-5p in tears of our GPAA patients and OHT individuals. Oxidative/nitrosative stress plays a pivotal role in POAG pathogenesis [12,15,39]; we may hypothesize that in the course of glaucoma, the pro-oxidants may induce the expression of miR-26b-5p, as reflected in this work [39][40][41]. Interestingly, hsa-miR-26b-5p has recently been involved in pseudoexfoliative glaucoma through interaction with the TGFR1 and TGR2 part of SMAD2 [42].
We also described herein that the miRN-27-a-3p showed a significant differential expression in tears from the POAG patients and the OHT participants. This miRNA has been implicated in cell proliferation and apoptosis in different types of cancer [43,44]. Overexpression of this miRNA promotes proliferation and inhibits apoptosis of cancer cells. Likewise, it has been shown that miR27a-3p induces ischemia-reperfusion injury by triggering oxidative stress [45]. In previous reports, it has been shown that the expression of this miRNA was increased in stretched versus control TM cells [46,47]. Consistent with these findings, we suggest that the miR-27a-3p overexpression in tears from POAG patients may be related to both the apoptotic and oxidative stress processes occurring in the POAG eyes [39][40][41].
Most target genes corresponding to the eight miRNAs showing significant differential expression profile in tears of the OHT individuals and POAG patients were related to extracellular matrix concerns and apoptosis processes (see the Tables 6-8).
Desjarlais et al. [48] reported the upregulation of miR-152-3p in the choroid during the vasodegenerative phase of an oxygen-induced retinopathy model compared to the control animals. However, these authors did not delve into the role of this miRNA, but just described the differential expression profile between groups. However, Wang et al. [49] pointed out that miR-152-3p may be an interesting molecular target for keloid treatment as a relevant regulator of cell proliferation, invasion, and extracellular matrix expression through targeting FOXF1 in keloid fibroblasts [49]. In this regard, we found the miR-152-3p upregulation in tears from the POAG group, supporting data previously reported by other researchers, emphasizing the importance of this miRNA in the trabecular meshwork changes and elevated IOP occurring in OHT and POAG [47].
The miR-30e-5p was differentially expressed in tears from OHT and POAG participants during the present work. However, only one article dealt with its role at the ocular level, highlighting its possible function as a predictor of the myasthenia gravis through the main ocular manifestation of the disease, the blepharoptosis [50]. No publications on the relationship between miR-30e-50 and glaucoma and/or other degenerative eye diseases have been found. Meng et al. [51] reported significant differential expression of the miR-30e-5p in patients with Alzheimer's disease compared to healthy controls. Additionally, Hardeland et al. [52] reported a negative connection between melatonin availability and neurodegenerative disorders, and in this regard, Scuderi et al. [53] reported that melatonin could protect the eye tissues from free radicals and pro-inflammatory mediators. Taking all these reports into consideration, we may hypothesize that miR-30e-5p may play a role in the ethiopathogenic mechanisms of glaucoma by its negative connection with melatonin.
Toro et al. [54] carried out a recent study in the vitreous humor of patients with different degrees of diabetic proliferative vitreoretinopathy (PVR) in which six miRNAs whose expression significantly increased with disease progression were reported, among them miR-224-5p. Dysregulation of its expression in PVR patients suggests a possible association with certain retinal degenerative processes. In the present study, we also found the upregulation of this miRNA in the POAG patients, pointing to a potential involvement of miR-224-5p in the onset and/or progression of the glaucoma ND.
The miR-151a-3p seems to be associated with a high risk of uveal melanoma [55]. Bianciotto et al. [56] reported that OHT is a risk factor for the development of uveal melanoma. In our OHT participants, a significantly higher expression of this miRNA was seen, suggesting that its downregulation in the POAG patients may be related to the onset of glaucoma ND. Nevertheless, the function of this miRNA in glaucoma needs further investigation.
Oxidative stress, inflammation, and apoptosis are relevant mechanisms involved in POAG [3][4][5][6][7][8]12,15,39]. In this regard, oxidative stress can trigger a wide spectrum of transcription factors such as p53, Nrf2, NF-κB, AP-1, PPAR-γ, and β-catenin/Wnt, among others [57][58][59][60]. When the above molecules undergo activation, they can induce the expression of hundreds of genes such as those related to pro-inflammatory cytokines and chemokines, cell adhesion molecules, metalloproteinases, cell cycle regulators, proapoptotic molecules, etc. [60][61][62][63][64]. In the case that the oxidative/inflammatory environment will be prolonged, this can chronically affect the ocular cells and tissues leading to apoptotic cell death and glaucoma ND. Given the target genes of the eight miRNAs with different expression between POAG and OHT, and following the methods explained above, we conducted a bioinformatic work to build several maps that showed the interactions between these genes as well as a connection network between the different metabolic pathways associated with these genes and miRNAs, as reflected in Figures 3-6. Target genes of the eight miRNAs that were identified in tears such as those showing a significant differential expression profile between the OHT and the POAG groups were mainly related to apoptosis and extracellular matrix functions as well as aqueous humor homeostasis, focal adhesion, oxidative stress, inflammation signaling and ND, as reflected in Figures 3-6.
Up to today, genetic testing can help identify people at risk for early-onset glaucoma types in almost 30% of cases [65]. A useful gene-based screening test for adult-onset POAG forms is not yet available [66], but a new test in blood or saliva by multi-trait analyses have found the risk of a person developing glaucoma [67]. Nevertheless, no way exists to predict the risk of a person with OHT to develop glaucoma.
We propose that the miRNAs identified in the present work and their target genes/ signaling pathways will allow for a better understanding of the molecular and genetic bases of glaucoma and, in this way, gain knowledge that leads us to develop a better diagnosis strategy for properly identifying HTO at a higher risk of glaucoma ND.
Individuals predisposed to glaucoma such as those suffering OHT should be closely followed to identify early changes to promptly start hypotensive therapy. Top advances in glaucoma genetics have provided an open window for stratifying the risk of glaucoma in HTO cases based on molecular-genetic predictions. Here, we show for the first time a new approach and future avenues dealing with miRNA expression in tears to better diagnose HTO individuals and identify those at highest risk of POAG as well as to better eye/vision care and blindness prevention.

Conclusions
In summary, eight miRNAs: the miR-26b-5p, miR-152-3p, miR-30e-5p, miR-125b-2-5p, miR-224-5p, miR-151a-3p, miR-1307-3p and the miR-27a-3p were found to be differentially expressed in the tears of patients in the interface between OHT and POAG. We demonstrated that some of these miRNAs, their target genes, and signaling pathways have already been related to apoptosis, extracellular matrix functions, inflammation, and oxidative stress as well as aqueous humor homeostasis and neurodegeneration. Further studies are needed to validate our results and deepen the knowledge of these miRNAs to discern its potential as biomarkers for the diagnosis and/or prognosis and future biotherapies for OHT and glaucomatous NDG.