Profiles of microRNA in aqueous humor of normal tension glaucoma patients using RNA sequencing

We aimed to identify and compare microRNAs (miRNAs) from individual aqueous humor samples between normal-tension glaucoma (NTG) patients and normal controls. Aqueous humor (80 to 120 µl) was collected before cataract surgery. Six stable NTG patients and seven age-matched controls were included in the final analysis. RNA sequencing was conducted for RNA samples extracted from the 13 aqueous humor samples, and bioinformatics analysis was employed for the miRNA targets and related pathways. Two hundred and twenty-eight discrete miRNAs were detected in the aqueous humor and consistently expressed in all samples. Eight significantly upregulated miRNAs were found in the NTG patients compared to the controls (fold-change > 2, p < 0.05). They were hsa-let-7a-5p, hsa-let-7c-5p, hsa-let-7f-5p, hsa-miR-192-5p, hsa-miR-10a-5p, hsa-miR-10b-5p, hsa-miR-375, and hsa-miR-143-3p. These miRNAs were predicted to be associated with the biological processes of apoptosis, autophagy, neurogenesis, and aging in the gene ontology categories. The related Kyoto encyclopedia of genes and genomes pathways were extracellular matrix-receptor interaction, mucin-type O-glycan biosynthesis, biotin metabolism, and signaling pathways regulating the pluripotency of stem cells. The differentially expressed miRNA in the NTG samples compared to the controls suggest the possible roles of miRNA in the pathogenesis of NTG. The underlying miRNA-associated pathways further imply novel targets for the pathogenesis of NTG.

Differential miRNA expression in aqueous humor from normal-tension glaucoma using RNA sequencing. The miRNA targets were analyzed based on the data from miRWalk 2.0. 2588 miRNAs were tested by RNA sequencing. And a total of 228 mature miRNAs were identified in the aqueous humor (AH) of the NTG patients (Fig. 1A). Up-regulated miRNAs are shown to the upper region of the plot (red) and down-regulated miRNAs are shown to the lower region of the plot (green). Of these, 8 miRNAs were significantly up-regulated compared to the controls (fold-change > 2 or < -2, p < 0.05), in addition to the analysis of the fold-changes were selected based on the normalized array data (log2) > 4 (Fig. 1B). These eight significantly up-regulated miRNAs were hsa-let-7a-5p, hsa-let-7c-5p, hsa-let-7f-5p, hsa-miR-192-5p, hsa-miR-10a-5p, hsa-miR-10b-5p, hsa-miR-375, and hsa-miR-143-3p (Table 2). There was no significantly down-regulated miRNAs identified in NTG patients compared to control subjects according to the same significance criteria. microRNA validation and biological interpretation of differentially expressed miRNAs. RNA sequencing revealed eight significantly differentially expressed miRNAs between AH of NTG patients and control. The heatmap diagram shows those 8 miRNAs increased in AH of NTG patients (red, high relative expression, Z-score) compared to control (blue, low relative expression, Z-score) ( Fig. 2A). To verify the results of RNA sequencing from the AH of NTG patients, let-7c-5p was analyzed using quantitative PCR (qPCR) (n = 7), and similar results were acquired ( Table 2). The expression of let-7c-5p increased significantly in the AH of the NTG patients compared to control (44.32 ± 3.63-fold, p < 0.001) (Fig. 2B).   www.nature.com/scientificreports/ To explore the effects of the expressed miRNAs, gene ontology (GO) categories were used. Fifteen major GO categories were randomly selected among numerous GO pathways. The percentage of the total significant and number of genes with differences in expression among each GO-related gene are presented in Fig. 3. The percentage refers to the proportion of microRNA modified in AH of NTG patients compared to control in the total miRNAs identified by researches in each GO category. Cell death-related categories, including apoptosis (3.27%)  www.nature.com/scientificreports/ and autophagy (4.63%), occupied the major proportion. Categories related to neurogenesis (5.88%) and the inflammatory response (4.47%) also represented significant proportions. Categories related to cellular function, including proliferation, migration, and differentiation, which may occur in any pathological condition, occupied 7.83%. Those microRNAs involved in all of the gene ontology categories associated with biological processes of apoptosis, autophagy, and neurogenesis were hsa-let-7c-5p and hsa-miR-375. As significantly down-regulated miRNAs were not identified according to the same significance criteria with up-regulation (Fig. 1B), no green graph as down-regulated GO category is shown in Fig. 3. The leading Kyoto encyclopedia of genes and genomes (KEGG) pathways, which include the predicted gene targets of each miRNA, are presented in Table 3. The analysis of gene-annotation enrichment was performed with the database for annotation, visualization, and integrated discovery (DAVID). The pathways relating to extracellular matrix (ECM)-receptor interaction (17.03, enrichment score, − log10 (p-value)), mucin-type O-glycan biosynthesis (5.40), biotin metabolism (3.73) and signaling pathways regulating pluripotency of stem cells (2.39) were significantly associated with miRNAs increased in the AH of the NTG patients. Among them, ECM-receptor interaction showed the most significantly related KEGG pathway in NTG (Fig. 4). The pathways such as protein digestion and absorption (1.62), and PI3K-Akt signaling (1.50) were also associated with up-regulated miRNAs.
Since hsa-let-7c-5p and hsa-miR-375 were both involved in all of the gene ontology categories associated with biological processes of apoptosis, autophagy, and neurogenesis, we focused on these two microRNAs and inspected the KEGG pathways of each microRNA. The increase in hsa-let-7c-5p was the most associated with ECM-receptor interaction (18.1), like as shown in the results of KEGG pathway analysis in the AH of the NTG patients. Furthermore, the pathways such as amoebiasis (3.19), and mucin-type O-glycan biosynthesis (2.97) were also associated with hsa-let-7c-5p (Table 4). The increase in hsa-miR-375 was associated with biotin metabolism (8.34), the thyroid hormone signaling pathway (2.39), and the glutamatergic synapse (1.20) (Table 5).

Discussion
To our knowledge, the present study was the first to report the microRNAs significantly differentially expressed in individual aqueous humor samples of NTG patients compared to controls in a single ethnic group of Asians. The present study did not pool the aqueous humor samples despite the scanty volume of each sample. RNA sequencing was performed to detect the microRNAs in each aqueous humor sample. RNA sequencing identified a total of 228 microRNAs in all aqueous humor samples. We detected 8 significantly upregulated microRNAs in NTG patients compared to the controls. They were hsa-let-7a-5p, hsa-let-7c-5p, hsa-let-7f-5p, hsa-miR-192-5p, Figure 3. Percentage and number of microRNAs with significantly changed expression among gene ontology category-related microRNAs. Gene ontology (GO) category of microRNAs with relatively large expression changes are identified. Up-regulated miRNAs are shown as red graph and down-regulated miRNAs are shown as green graph. Cell death-related categories including apoptosis (3.27%) and autophagy (4.63%) occupied the major proportion. Categories related to neurogenesis (5.88%) and inflammatory response (4.47%) also mainly occupied the proportion. The percentage refers to the proportion of microRNA modified in AH of NTG patients compared to control in the total miRNAs identified by researches in each GO category. As significantly downregulated miRNAs were not identified according to the same significance criteria with up-regulation, no green graph as down-regulated GO category is shown in this figure. Graphs were presented by ExDEGA v1.2.1.0 software. www.nature.com/scientificreports/ hsa-miR-10a-5p, hsa-miR-10b-5p, hsa-miR-375, and hsa-miR-143-3p. Among them, the microRNAs involved in all the gene ontology categories associated with biological processes of apoptosis, autophagy, and neurogenesis were hsa-let-7c-5p and hsa-miR-375. Significantly downregulated microRNA was not identified in the present study according to the same significance criteria with up-regulation, although we also analyzed for downregulation in the process of RNA sequencing. miRNAs have a significant part in the post-transcriptional modulation of gene expression and also are concerned in cellular functions including differentiation, growth, and cell death 27 . The human genome has approximately 2500 miRNAs that modulate the level of expression of > 60% of all genes that code proteins 19,28 .
Hsa-let-7c-5p miRNA is placed on human chromosome 21q21.1 on the LINC00478 gene. It has decreased expression and functions as a tumor suppressor gene for many malignancies such as colorectal adenocarcinoma, prostate, and acute myeloid leukemia 29,30 . . KEGG pathway. Enrichment score was represented as − log10 (p-value). The higher the enrichment score value is, the more significant the pathway is. The KEGG pathways relating to ECM-receptor interaction (hsa04512; 17.03), mucin-type O-glycan biosynthesis (hsa00512; 5.40), biotin metabolism (hsa00780; 3.73) and signaling pathways regulating pluripotency of stem cells (hsa04550; 2.39) were significantly associated with miRNAs increased in the aqueous humor of NTG patients. Note that the pathway relating to ECM-receptor interaction showed the most significant association with up-regulated miRNAs in NTG patients compared to other pathways. Data was analyzed by DianaTools. KEGG Kyoto Encyclopedia Genes and Genomes, ECM extracellular matrix. www.nature.com/scientificreports/ It has been reported that let-7 was downregulated in several malignancies and regulated apoptosis through direct targeting the oncogenes Myc, Ras, HMGA2, CDK6, and CDC25A 31 . MicroRNA let-7c has been reported to have an apoptotic effect on endothelial cells by suppressing anti-apoptotic protein Bcl-xl 32 . As one of the let-7 family, let-7a-5p has been reported to crosstalk with BCL-xL and lead to cell death and toxic autophagy in A549 lung cancer cells via activating the PI3K-signaling pathway, independent from apoptosis 33 . Injury to RGCs leads to glaucoma and the subsequent characteristic clinical findings of cupping in the optic disc and loss of visual field. The results of our study suggest a possible role for hsa-let-7c-5p in the pathogenesis of NTG on apoptosis and autophagy, which has not been reported before.
Hsa-miR-375 was the microRNA most significantly differentially expressed in the NTG patients compared to controls among all 13 aqueous humor samples (fold-change (log2) = 8.19, p = 0.000). Hsa-miR-375 was originally detected to be profusely expressed in isles of pancreas, and have a practical role in modulating glucosestimulated insulin exocytosis 34,35 . In a study regarding mouse hippocampus, hsa-miR-375 was shown to affect dendrite formation 36 .
In a previous study, hsa-miR-375 was demonstrated to regulate neuronal cell death induced by ethanol 37 . Moreover, it was revealed that miR-375 upregulation was related with shrinkage of brain and defects in cognition in chronic alcoholism 37 . In another previous work, ketamine induced upregulation of hsa-miR-375 in hESCderived neurons and it was concentration dependent 38 . The combined results of these studies 37,38 suggest that miR-375 upregulation is associated with pathological cortical neuron environments. Several previous studies have shown the neuronal function of hsa-miR-375, including effects on dendrite formation 36 and the modulation of neuronal apoptosis 39 . Previous work demonstrated that hsa-miR-375 inhibition showed epigenetic protection in hESC-derived neurons from ketamine induced neurotoxicity, probably via the direct and inverse modulation of the brain-derived neurotrophic factor (BDNF) gene 38 .
In a previous study, 11% of aged marmosets demonstrated NTG-like degeneration of the optic disc, which was a similar rate with that of glaucoma patients. The marmosets did not show elevated IOP but demonstrated elevated oxidative stress levels, low BDNF and TrkB expression, and low cerebrospinal fluid (CSF) pressure in the retina, optic nerve head, and CSF 40 . The expression of BDNF and TrkB were decreased in the optic nerve and retina in glaucomatous marmosets, in accordance with the data assembled from observations of glaucoma patients 41 . Neurotrophins are one of therapeutic candidates for glaucoma and several studies have shown that BDNF eye drops rescued visual function in a mouse glaucoma model with high-IOP (DBA/2J mice) 42,43 . Considering the direct regulation of hsa-miR-375 effect on the BDNF gene, the results of our study may also support the BDNF gene as a potential target for glaucoma gene therapy, especially for NTG.
RGC dendrites are considered to respond and demonstrate compensatory changes after damage, for example, extending the dendritic perceptive fields and generating new branches [44][45][46] . A recent study by Park et al. reported that synaptic vesicle proteins in bipolar cells and RGCs were augmented following BDNF application in a chronic IOP elevation rat model 47 . In addition, BDNF improved the number of synapses connecting bipolar cells and the RGCs in the retina of glaucomatous change. The authors concluded that these findings concerning the ability of BDNF to encourage useful synaptic changes may assist in the development of neuroenhancement approaches for the treatment of disturbed synaptic function in glaucoma. Our results also suggest that dendrite formation and the related BDNF involvement could serve as novel treatment targets for glaucoma, especially for NTG, regarding the neuronal functions of hsa-miR-375. Moreover, RGC dendrite formation and the related BDNF action, which are gene ontology categories associated with neurogenesis, may be involved in the pathogenesis of NTG.
The inflammatory response was one of the gene ontology biologic process categories potentially affected by microRNAs in the aqueous humor of NTG patients (Fig. 3). The binding of hsa-miR-375 to toll-like receptor 4 (TLR4) was verified by a dual-luciferase reporter gene assay 48 . TLR4 activates nuclear factor-kappa B (NF-κB) through MyD88, stimulating the secretion of downstream inflammatory cytokines 49 . Studies of patients with glaucoma and glaucoma experimental models indicated that the inflammatory response is mediated partly by toll-like receptors (TLRs) 50 . Studies of anatomy have reported elevated levels of TLR2, TLR3, and TLR4 in microglia and astrocytes in the retina of glaucoma patients 51 . Another previous study demonstrated that the inhibition of astroglial NF-kB decreased the inflammatory conditions and enhanced the survival of RGCs following retinal ischemia 52 . The association between hsa-miR-375 and TLR and the inflammatory response may partially explain the results of our study in NTG patients.
Aging was also one of the gene ontology biological process categories potentially affected by microRNAs in the aqueous humor of NTG patients (Fig. 3). As glaucoma is an age-dependent disease and the prevalence of glaucoma increases with age as reported in many population-based studies worldwide and also in Asians 9 , the results of our study seem reasonable.
The other gene ontology categories of biologic processes affected by microRNA in the aqueous humor of NTG patients compared to control were cellular proliferation, migration, differentiation, and cycle, DNA repair, secretion, and angiogenesis. These categories may act in common in any pathological conditions and may not occur solely in glaucoma. In this regard, such categories are not discussed in detail in the present study.
Degeneration of optic nerve fiber originally takes place at the site of lamina cribrosa in glaucoma 53 . The lamina cribrosa is organized by the ECM and inactive astrocytes, and acts as a fibroelastic structure, offering biological and mechanical support for axons of optic nerve head 54,55 . Glaucoma brings about remodeling of ECM and the activation of astrocytes. Consequently, the reactive astrocytes express new ECM proteins, which are regarded to change its configuration or be neurotoxic to the RGCs 56 . A study that used bioinformatics analysis to detect differentially expressed glaucoma genes found that one of the most significantly enriched pathways was ECM-receptor interaction 57 . This finding is consistent with the results of the current study, which displayed that the most significant related KEGG pathway was ECM-receptor interaction in NTG patients (p-value (− log10) = 17.027). Compared to ECM-receptor interaction, other KEGG pathways showed p-value (− log10) of less than 5.402 and the least three were 1.378 (Table 3 and Fig. 4). www.nature.com/scientificreports/ When we analyzed the KEGG pathway of each differentially upregulated let-7c-5p and miR-375, the most significantly related KEGG pathway in let-7c-5p was ECM-receptor interaction (p-value (− log10) = 18.109) ( Table 4), which is in accordance with the results of all of differentially expressed microRNAs. We also confirmed the significant relative expression of let-7c-5p compared to the controls by qPCR (Fig. 2B). These validation results of confirmation of let-7c-5p which showed the same most significantly related KEGG pathway indicate the reliability and validity of our RNA sequencing results.
The pathogenesis of NTG is still under debate and it is a particular issue because NTG develops even within normal-range IOP. Studies using OCT imaging devices have reported differences in lamina cribrosa thickness and curvature in NTG patients [58][59][60] . The thickness of the lamina cribrosa was thinner in NTG than in POAG patients 58 . Moreover, the lamina cribrosa was thinner in VF-affected unilateral NTG eyes than in the fellow normal eyes of unilateral NTG patients or normal eyes 59 . In unilateral NTG patients, eyes with glaucomatous damage showed more steeply curved lamina cribrosa than the fellow healthy eyes 60 . These findings suggest that the lamina cribrosa undergoes significant remodeling in NTG eyes. Considering that the lamina cribrosa is the primary site of glaucomatous injury and NTG results in the remodeling of the lamina cribrosa, which is organized by the ECM, our microRNA analysis results indicate the important role of ECM alterations in the pathogenesis of NTG.
The other KEGG pathways possibly affected by microRNAs in the aqueous humor of NTG patients in our study included mucin-type O-glycan biosynthesis, biotin metabolism, and signaling pathways regulating the pluripotency of stem cells. In this current study, several gene ontology biological process categories and related KEGG pathways were associated with the differentially expressed microRNAs in the aqueous humor of the NTG patients. Confirmation of individual biological processes or KEGG pathways associated with NTG needs further clinical and experimental studies. However, our study is unique in that significantly differentially expressed microRNAs were identified in the individual sample of aqueous humor of NTG patients compared to the controls in a single ethnic group of Asians.
Since our study is the first to report the differentially expressed microRNAs in NTG patients and therefore, there is not yet a previous study to compare the results from. But we can compare our results from previous POAG patients with elevated IOP. Recent review article that investigated the mechanisms of microRNA in POAG described that microRNAs were involved in the regulation of IOP and related to optic nerve damage factors such as mechanical stress, hypoxia and inflammation 61 74,75 . These microRNAs investigated in the pathogenesis of POAG were not identified in the current study in NTG patients with normal range of IOP, which suggest that the differentially expressed microRNAs in the pathogenesis of NTG are different from POAG with high IOP. It may also partially indicate the reliability and validity of the present study to some extent.
The present exploratory study has limitations due to the relatively small sample size and the small volume of the aqueous humor samples, but it has a significant meaning in that it provides the potential for further research field of microRNA in NTG, which is prevalent in Asians. Since the volume of the aqueous humor samples was not enough for all miRNAs to undergo qPCR for validation, only hsa-let-7c-5p underwent qPCR, which showed a significant relative expression pattern consistent with the RNA sequencing results. Considering that hsa-let-7c-5p showed a much smaller fold-change (fold-change (log2) = 2.92) than that of hsa-miR-375 (fold-change (log2) = 8.19) in the NTG patients compared to the controls in RNA sequencing (Table 2), the consistent qPCR results (Fig. 2B) of hsa-let-7c-5p may indicate the reliability of our RNA sequencing results. The influence of hypotensive topical medications on microRNA expression within the aqueous humor of NTG patients has yet to be determined. The impact of using different hypotensive topical medications on our results is not known. Future studies with large numbers of samples would benefit from controlling the use of hypotensive topical medications.
In conclusion, we found significantly differentially expressed microRNAs in the individual aqueous humor of NTG patients compared to controls in a single ethnic group of Asians, which has not been investigated previously. The findings of the present study suggest a possible role of microRNA in the pathogenesis of NTG. Further studies with more cases should be performed to draw more definitive conclusions. In addition, microRNA in individual aqueous humor may have a further potential as biomarkers and novel targets for the pathogenesis of NTG.

Material and methods
Ethics statement. This study was conducted according to the tenets of the Declaration of Helsinki for research involving human subjects. The current study was approved by the Institutional Review Board of Gyeongsang National University Changwon Hospital, Gyeongsang National University, School of Medicine (GNUCH-2019-06-001-002). Written informed consent was obtained from all subjects enrolled in this study. All methods were carried out in accordance with relevant guidelines and regulations.
Patient selection and acquisition of aqueous humor samples. Samples of aqueous humor were obtained from patients who underwent uneventful phacoemulsification for routine cataract surgery after obtaining written informed consent. Six patients with NTG stably managed with only topical medication and seven age-matched control subjects agreed to take part in the present study. Approximately 80 to 120 µl of aqueous humor was obtained by corneal paracentesis with a 30-gauge needle before the initial main cataract incision at the beginning of surgery. Anterior chamber paracentesis was performed in the operating room under aseptic sterile conditions. The collection of aqueous humor was performed without trauma in all subjects, thus, removing any possibility of contamination with cellular remains or blood. All of collected samples were completely anonymized and immediately snap-frozen with liquid nitrogen, then transferred to the research laboratories. The clinical data was obtained from the electronic medical record and collected in a completely anonymized way. The clinical data collected were age, sex, eye laterality of right or left, baseline IOP, topical eye drops used, and ocular comorbidities.  76 . In brief, for library construction, 180 pg of total RNA from each sample was used to ligate 1ug of adaptors, and then cDNA was synthesized using reverse-transcriptase with adaptor-specific primers. PCR was performed for library amplification and the libraries were cleaned-up using a QIAquick PCR Purification Kit (Qiagen, Inc, Germany) and AMPure XP beads (Beckman Coulter, Inc., Pasadena, CA, USA). The yield and size distribution of the small RNA libraries were evaluated by the Agilent 2100 Bioanalyzer instrument for the High-sensitivity DNA Assay (Agilent Technologies, Inc., USA). High-throughput sequences were produced by the NextSeq500 system by single-end 75 sequencing (Illumina, San Diego, CA, USA).
microRNA validation by quantitative real-time PCR. cDNA synthesis and real-time PCR were conducted employing the miScript PCR system (Qiagen, Venlo, The Netherlands). cDNA was synthesized from 357 pg of RNA employing the miScript II RT Kit with HiSpec buffer following the manufacturer's instructions. The microRNA cDNA was amplified employing the following primer pair: hsa-let-7c-5p (Hs_let-7c_1, MS00003129) and internal control hsa-U6 (Hs_RNU6-2_11, MS00033740 Statistical analysis. Enrichment p-values were corrected for false discovery rate (FDR) 82 . Two sample t-test was employed for quantitative PCR validation. Data of microRNA validation are expressed as the mean ± standard error of the mean (S.E.M.). Statistical significance for the validation was determined using Unpaired Student's t-test (Prism 5; GraphPad Software, La Jolla, CA, USA). p < 0.05 was considered to indicate a statistically significant difference.