microRNA Profiling in Glaucoma Eyes With Varying Degrees of Optic Neuropathy by Using Next-Generation Sequencing

Abstract

Purpose: To explore the microRNA (miRNA) profile and its putative role in glaucomatous optic neuropathy by using next-generation sequencing.

Methods: Aqueous humor (AH) samples were collected from 19 primary open-angle glaucoma (POAG) eyes and 17 cataract eyes before surgery. Next-generation sequencing was performed for RNA samples extracted from 18 AH samples, and the bioinformatics approach was applied for samples with adequate clean data output. The other 18 samples were used for quantitative PCR validation of sequencing results.

Results: In total, 12 (six POAG and six cataract controls) samples with sufficient clean data output after sequencing were used for further data analysis. Four hundred sixty-six and 480 mature miRNAs were detected in the POAG and cataract control groups, respectively. Among them, 164 miRNAs were detected in all POAG samples, and 96 miRNAs were detected in all cataract control samples. Furthermore, 88 miRNAs were identified as differently expressed between POAG and cataract control eyes. In addition, 16 miRNAs were differently expressed between POAG eyes with severe visual field damage and eyes with moderate visual field damage. This differential expression was predicted to regulate thiamine metabolism, purine metabolism, and transcriptional misregulation. Relative expression patterns of hsa-miR-184, hsa-miR-486-5p, and hsa-miR-93-5p were confirmed by quantitative PCR.

Conclusions: This study comprehensively demonstrated the miRNA expression profile in the AH of POAG eyes, especially the differential expression of miRNA in eyes with varying degrees of visual field damage, which, together with the underlying miRNA-related pathways, indicate new targets for the pathogenesis and progression of POAG.

MicroRNAs (miRNAs) are small, endogenous, noncoding RNA molecules consisting of 19 to 22 nucleotides (nt).¹ They are regulators of mRNA transcription and translation, which are stable and often specifically enriched in a particular tissue or during essential cellular processes.² More and more extracellular miRNAs from biofluids have been identified as biomarkers for cancer, cardiovascular disease, diabetes, ocular disease, and many other disorders.^3–11 

Tears, aqueous humor (AH), and vitreous humor are three sources of fluids that contain extracellular miRNAs in the eyes.^12–15 Among them, the AH is the most frequently collected ocular fluid sample for examinations of ocular diseases.¹⁶ The change of contents in the AH, including proteins, cells, and other substances such as miRNAs, which not only reflect the abnormality of anterior ocular segments but also indicate damage to the retina and choroid, serve as a good sample for finding biomarkers for eye diseases. 

Glaucoma is one of the leading causes of blindness worldwide.¹⁷ The common features of all forms of glaucoma are the loss of retinal ganglion cells, thinning of the retinal nerve fiber layer (RNFL), and cupping of the optic disc.¹⁸ With early detection and proper treatment, it is possible to slow the progression of this damage. However, despite controlling intraocular pressure (IOP), some patients' glaucoma continues to progress.¹⁹ There must be underlying mechanisms for this continued progression. Only a few studies about the miRNA in the AH of glaucoma eyes by using miRNA arrays have been reported.^12,20 However, miRNA arrays can detect only part of mature miRNAs in the miRBase database, whereas next-generation sequencing can detect all mature miRNAs and even novel miRNAs that have never been found. 

In this study, we used next-generation sequencing to analyze individual AH samples from primary open-angle glaucoma (POAG) patients with varying degrees of visual field damage and from cataract controls to identify specific miRNA expression characteristics in POAG eyes and certain biomarkers for the pathogenesis and progression of POAG. 

Subjects and Methods

Subjects

This research was conducted in accordance with the Declaration of Helsinki and with the approval of the Zhongshan Ophthalmic Center Institutional Review Board. The written informed consent was obtained from patients undergoing glaucoma and cataract surgery at Zhongshan Ophthalmic Center, Guangzhou, China. 

The diagnosis criteria of POAG was the same as in our previous study.²¹ All POAG patients underwent serial ophthalmic examinations by an experienced glaucoma specialist. The visual field was detected by the Humphrey computerized automated perimeter (Carl Zeiss Meditec, Jena, Germany). On the basis of the visual field data, the POAG patients enrolled in this study were divided into two groups according to Mills et al.²² Staging System: (1) POAG patients with moderate visual field defect (M-POAG) and (2) POAG patients with severe visual field defect (S-POAG). In addition, the RNFL thickness of POAG eyes was detected by the Heidelberg Spectralis (Heidelberg Engineering, Heidelberg, Germany). 

All the enrolled POAG patients underwent trabeculectomy. The drugs used before the operation were recorded. The control group was selected according to the following criteria: (1) underwent routine cataract surgeries (phacoemulsification plus intraocular lens implantation), and (2) their IOP was within the normal range, and they were free of other ocular or systemic diseases. After cataract surgery, the fundus of every enrolled patient was checked to exclude the glaucomatous neuropathy. 

AH Samples

Fifty to 100 μL AH were collected at the outset of surgery from 19 POAG eyes and 17 cataract control eyes in accordance with the procedure described in our previous study.²³ To avoid the influence of surgical traumatic damage to the blood-aqueous barrier, all samples were obtained prior to any conjunctival or intraocular manipulation. All samples were immediately stored at −80°C for further analysis. 

RNA Extraction

AH samples were centrifuged successively (3,000g for 10 minutes at 4°C and 12,000g for 10 minutes at 4°C) to remove cells and cell debris. Synthetic Caenorhabditis elegans miR-39 spike-in control was added to the supernatant, and RNA was extracted using the TRIzol reagent (TRIzol LS Reagent; Life Technologies Corp., Carlsbad, CA, USA) according to the manufacturer's instructions. 

Library Preparation and Sequencing

All extracted RNA from each individual eye was used to generate sequencing libraries with the NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England BioLabs, Inc., Ipswich, MA, USA). RNA segments of different sizes were separated by polyacrylamide gel electrophoresis. An 18 to 30 nt stripe was recycled and ligated with specific adapters to the 3′ and 5′ end. Adapter-ligated libraries were reverse transcribed and amplified by PCR, with primers containing index sequences specific for each sample. The purified PCR construct was recovered by polyacrylamide gel electrophoresis. Recycled products were dissolved in an ethidium bromide solution and labeled. An Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA) and an ABI StepOnePlus Real-Time PCR System (Life Technologies Corp., Carlsbad, CA, USA) were used to test quality and production of the prepared library. Individual libraries were pooled and sequenced on an Illumina sequencing platform (HiSeq 4000 platform; Illumina, Inc., San Diego, CA, USA). 

Sequencing Data Analysis

After obtaining raw data, we conducted bioinformatics analysis according to previous studies.^24–34 The 49 nt sequence tags from next-generation sequencing went through the data cleaning analysis, which included getting rid of the low-quality tags, or the 5′ adaptor contaminants from the 50 nt tags to obtain clean, credible tags. After filtering, the remaining clean tags were stored in the FASTQ format. Samples with adequate data output (clean tag count, >20,000,000) were used for further data analysis. Then, the length distribution of the clean tags and common and specific sequences between samples were summarized. Next, anchor alignment-based small RNA annotation was used to map clean reads to the reference genome and mirBase. The characteristic hairpin structure of the miRNA precursor can be used to predict novel miRNAs. We used miRDeep2 to predict novel miRNAs by exploring the secondary structure, the Dicer cleavage site, and the minimum free energy of the unannotated sRNA tags, which could be mapped to the genome. 

The miRNA expression level is calculated by using the Transcripts Per Kilobase Million method. Only miRNAs annotated as mature in mirBase data (mirbase version 21.0) were included. The NOISeq method was used to screen differently expressed miRNAs between two groups. First, NOISeq uses the sample's small RNA expression in each group to calculate log2(fold change) M and the absolute different value D of all pair conditions to build the noise distribution model. Second, for small RNA-A, NOISeq computes its average expression “Control_avg” in the control group and average expression “Treat_avg” in the treatment group. Then, the fold change M (M_A = log2[(Treat_avg)/(Control_avg)]) and absolute different value D (D_A = |Control_avg − Treat_avg|) is obtained. If M_A and D_A markedly diverge from the noise distribution model, gene A will be defined as a differently expressed gene. There is a probability value to assess how M_A and D_A both diverge from the noise distribution model. Finally, we screen differently expressed small RNA according to the following default criteria: fold change ≥4 and diverge probability ≥0.8. After getting the miRNA result, the target prediction for miRNAs and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway for target genes were analyzed. We used the tools RNAhybrid, miRanda, and TargetScan to find the possible miRNA targets and took intersection targets for further analysis. In order to predict the functions of these miRNAs in the regulation of disease pathogenesis, KEGG was used to perform pathway enrichment analysis. A P value of less than 0.01 and a false discovery rate of less than 0.01 were used as the threshold to select significant KEGG pathways. 

Quantitative PCR (qPCR)

AH samples of 18 independent patients (including 10 AH samples from POAG eyes and 8 AH samples from cataract controls), obtained as above, were cleared of potential debris by centrifugation (3000g for 10 minutes at 4°C and 12,000g for 10 minutes at 4°C). A C. elegans miR-39 spike-in control was added to the supernatant, and small RNAs were extracted using an extraction kit (miRNeasy Plasma Kit; Qiagen, Hilden, Germany) in accordance with the manufacturer's instructions. Retrotranscription primers (stem loop) and qPCR primers (forward and reverse) for each miRNA were designed by RiboBio (Guangzhou, China).³⁵ Subsequently, a commercial kit (miDETECT A Track miRNA qRT-PCR Starter Kit, RiboBio, Guangzhou, China) was used for reverse transcription and qPCR detection of selected miRNAs in accordance with the manufacturer's instructions. 

General Statistical Analysis

Basic clinical data in this study were processed and analyzed by the IBM SPSS-20 statistical package (IBM SPSS Statistics for Windows, version 20.0; IBM Corp., Armonk, NY, USA). Demographic data as well as clinical measurements were tabulated. The significance of differences among groups was determined using the χ² test for categorical variables and ANOVA for normally distributed variables. P < 0.05 was considered statistically significant. 

The necessary sample size was also calculated. With adequate data output (clean tag count, >20,000,000), we assumed the same depth of coverage μ of the miRNA as 20,000, and within group coefficient of variation (σ) of 0.60 (for human samples, it will often be 0.4–1.0),³⁶ with the usual false positive rate of 0.05 (two-sided α) and a power (1−β) of 80%, in order to target a four-fold change (Δ). We determined that three samples per group were needed using the following formula³⁶:  

\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\unicode[Times]{x1D6C2}}\)\(\def\bupbeta{\unicode[Times]{x1D6C3}}\)\(\def\bupgamma{\unicode[Times]{x1D6C4}}\)\(\def\bupdelta{\unicode[Times]{x1D6C5}}\)\(\def\bupepsilon{\unicode[Times]{x1D6C6}}\)\(\def\bupvarepsilon{\unicode[Times]{x1D6DC}}\)\(\def\bupzeta{\unicode[Times]{x1D6C7}}\)\(\def\bupeta{\unicode[Times]{x1D6C8}}\)\(\def\buptheta{\unicode[Times]{x1D6C9}}\)\(\def\bupiota{\unicode[Times]{x1D6CA}}\)\(\def\bupkappa{\unicode[Times]{x1D6CB}}\)\(\def\buplambda{\unicode[Times]{x1D6CC}}\)\(\def\bupmu{\unicode[Times]{x1D6CD}}\)\(\def\bupnu{\unicode[Times]{x1D6CE}}\)\(\def\bupxi{\unicode[Times]{x1D6CF}}\)\(\def\bupomicron{\unicode[Times]{x1D6D0}}\)\(\def\buppi{\unicode[Times]{x1D6D1}}\)\(\def\buprho{\unicode[Times]{x1D6D2}}\)\(\def\bupsigma{\unicode[Times]{x1D6D4}}\)\(\def\buptau{\unicode[Times]{x1D6D5}}\)\(\def\bupupsilon{\unicode[Times]{x1D6D6}}\)\(\def\bupphi{\unicode[Times]{x1D6D7}}\)\(\def\bupchi{\unicode[Times]{x1D6D8}}\)\(\def\buppsy{\unicode[Times]{x1D6D9}}\)\(\def\bupomega{\unicode[Times]{x1D6DA}}\)\(\def\bupvartheta{\unicode[Times]{x1D6DD}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bUpsilon{\bf{\Upsilon}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\def\iGamma{\unicode[Times]{x1D6E4}}\)\(\def\iDelta{\unicode[Times]{x1D6E5}}\)\(\def\iTheta{\unicode[Times]{x1D6E9}}\)\(\def\iLambda{\unicode[Times]{x1D6EC}}\)\(\def\iXi{\unicode[Times]{x1D6EF}}\)\(\def\iPi{\unicode[Times]{x1D6F1}}\)\(\def\iSigma{\unicode[Times]{x1D6F4}}\)\(\def\iUpsilon{\unicode[Times]{x1D6F6}}\)\(\def\iPhi{\unicode[Times]{x1D6F7}}\)\(\def\iPsi{\unicode[Times]{x1D6F9}}\)\(\def\iOmega{\unicode[Times]{x1D6FA}}\)\(\def\biGamma{\unicode[Times]{x1D71E}}\)\(\def\biDelta{\unicode[Times]{x1D71F}}\)\(\def\biTheta{\unicode[Times]{x1D723}}\)\(\def\biLambda{\unicode[Times]{x1D726}}\)\(\def\biXi{\unicode[Times]{x1D729}}\)\(\def\biPi{\unicode[Times]{x1D72B}}\)\(\def\biSigma{\unicode[Times]{x1D72E}}\)\(\def\biUpsilon{\unicode[Times]{x1D730}}\)\(\def\biPhi{\unicode[Times]{x1D731}}\)\(\def\biPsi{\unicode[Times]{x1D733}}\)\(\def\biOmega{\unicode[Times]{x1D734}}\)\begin{equation}{\rm{n}} = 2{\left( {{{\rm Z}_{1 - {\alpha \over 2}}} + {{\rm Z}_\beta }} \right)^2}{{({1 \over \mu } + {\sigma ^2})} \over {({{\log }_e}{\Delta ^2})}}\end{equation}

The tests were conducted using the statistical package R, version 3.1.2.  

The specific R package is also available via Bioconductor (available in the public domain, http://bioconductor.org/packages/release/bioc/html/RNASeqPower.html). 

Results

Basic Clinical Characteristics of the Subjects

In total, RNA samples from nine POAG eyes and nine cataract controls were sequenced. However, three POAG eyes and three cataract controls were excluded for data analysis due to low data output (clean tag count, <20,000,000). The baseline characteristics of the remaining six POAG eyes and six cataract controls are summarized in Table 1. The age of the subjects in the cataract controls (n = 6), POAG (n = 6), S-POAG (n = 3), and M-POAG (n = 3) groups were 66.7 ± 7.6, 66.5 ± 7.3, 66.0 ± 8.2, and 67.3 ± 8.7 years, respectively, and there were no significant differences between the groups (P > 0.05). The mean IOPs of the four groups were 13.3 ± 3.9, 20.0 ± 6.0, 19.3 ± 7.5, and 20.7 ± 5.5 mm Hg, respectively. The IOP of the POAG patients was higher than that of the cataract eyes (P = 0.045). POAG patients were using a mean of 2.5 topical pressure-lowering agents. The patients of the M-POAG and S-POAG groups had visual field defects of −8.9 ± 2.9 dB mean deviation (MD) and −30.6 ± 1.3 dB (MD), respectively (P < 0.001). The average RNFL thickness of the M-POAG groups was 55.0 ± 13.7 μm and that of the S-POAG groups was 38.7 ± 7.8 μm (P = 0.211). 

Table 1

 

Baseline Demographics and Characteristics of Participants

 

Characteristics of Next-Generation Sequencing and miRNA Expression

The mean RNA yield of the entire AH sample was 0.38 ± 0.15 ng. To determine the miRNA expression profile in POAG, next-generation sequencing was performed on independently pooled RNA samples with index sequences specific for each sample from the POAG and cataract groups. After sequencing, samples with adequate data output (clean tag count, >20,000,000), including six AH samples of POAG and six samples of cataract, were used for data analysis. On average, 21,867,979 (±1,963,465) clean reads were obtained per sample among the remaining six POAG and six cataract samples. When restricted to sequences of 20 nts and longer, 2,798,234 (±1,981,378) of the reads were mapped to miRBASE entries. A total of 164 mature miRNAs were detected in all six AH samples from POAG patients and 96 miRNAs in all six cataract samples. Of these, 84 miRNAs were detected in both groups. The top 20 miRNAs with high expression in all samples are listed in Figure 1, and miR-184 was the most abundantly expressed miRNA in our study. Moreover, in total, 426 novel miRNAs were predicted in the AH of all samples (Supplementary Sheet S1). The expression mode of these predicted miRNAs was also analyzed (Supplementary Sheet S2, S3). 

Figure 1

 

The top 20 miRNAs with high expression in AH of POAG and cataract eyes.

To evaluate the significance of differently expressed miRNAs, we calculated the probability of being differently expressed using NOISeq scores, M-value, and D-value. Based on these calculations, 88 miRNAs showed significantly different expression with a fold change of ≥4 and a diverge probability of ≥0.8 when comparing the POAG group with the cataract group (Table 2). Of these, 73 miRNAs were upregulated (those with a positive M-value) and 15 were downregulated (those with a negative M-value), as shown in the volcano plot (Fig. 2). 

Table 2

 

Differently Expressed miRNA Between POAG and Cataract Groups

 

Figure 2

 

Volcano plot showing differentially expressed miRNAs. Eighty-eight miRNAs were identified as differentially expressed, with a fold change greater than 4 and a probability score greater than 0.8 between POAG and cataract eyes. Red dots represent significant upregulated miRNAs in POAG eyes, and green dots are significant downregulated miRNAs.

When comparing the S-POAG with the M-POAG group, 16 miRNAs were significantly upregulated, with a fold change of ≥4 and a diverge probability of ≥0.8, as summarized in the volcano plot (Table 3; Fig. 3). After intersection analysis, 12 of 16 differently expressed miRNAs between the S-POAG and M-POAG groups were also differently expressed between the POAG and cataract groups, as shown in the Venn diagram (Fig. 4). 

Table 3

 

Differently Expressed miRNA Between S-POAG and M-POAG Groups

 

Figure 3

 

Volcano plot showing differentially expressed miRNAs. Sixteen miRNAs were differentially expressed, with a fold change greater than 4 and a probability score greater than 0.8 between POAG eyes with severe visual field damage (PD < −30dB) and with modest visual field damage (−6Db < PD < −12dB). Red dots represent significant upregulated miRNAs in POAG eyes.

Figure 4

 

Twelve of 16 differentially expressed miRNAs between S-POAG and M-POAG groups were differentially expressed simultaneously between POAG and cataract groups. Red: differentially expressed miRNAs between POAG and cataract groups. Green: differentially expressed miRNAs between S-POAG and M-POAG groups. Intersection: miRNAs differentially expressed in both comparisons (POAG vs. cataract, S-POAG vs. M-POAG).

Target Genes and Functional Analysis

To explore the various targets and their associated pathways, a putative target list was obtained. We performed pathway enrichment analysis of the differently expressed sequences' target genes based on the KEGG database³³ and generated a report for the differently expressed sequences' target genes in each pairwise (POAG vs. cataract, S-POAG vs. M-POAG). The top three KEGG pathways of these two comparisons are shown in Tables 4 and 5, respectively. In addition, we generated a scatter plot for the top 20 of KEGG enrichment results (S-POAG vs. M-POAG), as shown in Figure 5. 

Table 4

 

Significant Pathways Predicted to be Highly Regulated by Differently Expressed miRNAs (POAG vs. Cataract)

 

Table 5

 

Significant Pathways Predicted to be Highly Regulated by Differently Expressed miRNAs (S-POAG vs. M-POAG)

 

Figure 5

 

The top 20 of KEGG enrichment results (S-POAG vs. M-POAG).

miRNA Target Gene Regulation Network

To better understand the role of miRNAs in varying degrees of glaucomatous optic neuropathy, we constructed a miRNA target gene regulation network. The top three KEGG pathways (thiamine metabolism, purine metabolism, and transcriptional misregulation in cancer) of comparison between S-POAG and M-POAG are shown in Figure 6. The blue nodes are putative target genes. MiR-16-5p, miR-205-5p, miR-206, miR-200a-3p, and miR-200b-3p are differently expressed miRNAs that can regulate these three pathways. 

Figure 6

 

The miRNA target gene regulation network for the top three KEGG pathways of comparison between S-POAG and M-POAG.

Verification of Next-Generation Sequencing Data by Real-Time PCR

To reassess possible differences in miRNA abundance suggested by next-generation sequencing, three miRNAs were analyzed by qPCR in 18 independent AH samples. According to our sequencing results, we chose to study miR-184 as the most abundant in sequencing, miR-486-5p as of intermediate abundance, and miR-93-5p as of low abundance. The expression of hsa-miR-184, hsa-miR-486-5p, and hsa-miR-93-5p detected using qPCR was accordant with the sequencing results (Fig. 7). 

Figure 7

 

The relative expression levels of miRNAs were assayed by qPCR and normalized to levels of a C. elegans miR-39 spike-in control. Then, the expression fold changes of hsa-miR-184 and hsa-miR-486-5p were normalized to that of hsa-miR-93-5p. The expression mode of these miRNAs detected using qPCR was accordant with the sequencing results.

Discussion

In the present study, we identified miRNAs profiles in AH from POAG and cataract patients by using next-generation sequencing. Eighty-eight miRNAs showed significantly different expressions when comparing the POAG group with the cataract group; this differential expression may regulate glycosaminoglycan biosynthesis and the longevity regulating pathway. When comparing the S-POAG with the M-POAG group, 16 miRNAs were significantly differently expressed, and predicted pathways of thiamine metabolism, purine metabolism, and transcriptional misregulation may play roles in the pathogenesis and progression of POAG. Moreover, 426 novel miRNAs were predicted to exist in POAG and cataract AH, but these novel miRNAs must be further studied. 

miRNAs in AH can not only be influenced by damaging of the anterior ocular segment but also by damaging of the posterior segment. For instance, the concentration of inflammation factors was found to be especially higher in the AH, such as in diabetes retinopathy and age-related macular degeneration eyes, suggesting that retinal damage could be reflected in the AH.^37–39 Cell fragments from retina-like photoreceptor outer segments can also be found in AH, which suggests that miRNAs in the AH may also come from the posterior segment of the eye.⁴⁰ 

In our study, miR-184 was the most abundantly expressed miRNA among all 12 AH samples of POAG and cataract patients (Fig. 1). In accordance with our study, several prior studies regarding aqueous miRNAs of glaucoma or cataract patients have detected miR-184 with high expression.^{12,14,20,41,42} However, the differently expressed miRNAs found in the AH between glaucoma and cataract patients in three studies, including ours and two other previous studies using miRNA arrays, had no overlap with each other.^12,20 This may have resulted from the inherent inconformity of miRNA expression detection among different platforms, including miRNA microarray systems, next-generation sequencing, and qPCR, which were systematically analyzed by Anna Git and colleagues.⁴³ The different subtypes of glaucoma patients (including POAG and primary angle closure glaucoma) enrolled in prior studies may also lead to this discordance. Therefore, it is essential to conduct follow-up studies of these differently expressed miRNAs. 

POAG is a currently irreversible optical neurodegenerative disease. Using IOP lowering, the progression of the disease can be halted in some patients. However, some patients suffer continuous visual field damage even with radical IOP lowering.^44,45 There must be underlying mechanisms for this disparity. The expression profile of miRNAs in AH from POAG eyes with varying degrees of visual field defects may offer certain clues. Several miRNAs were found to be differently expressed in the POAG, cataract, S-POAG, and M-POAG groups (Tables 2, 3). Eighty-eight miRNAs showed significantly different expressions when comparing the POAG group with the cataract group (Table 2). Sixteen miRNAs were significantly upregulated in S-POAG patients compared with M-POAG patients. After intersection analysis, 12 of the 16 differently expressed miRNAs between the S-POAG and M-POAG groups were also differently expressed between the POAG and cataract groups, and these miRNAs are most likely implicated in the optic neuropathy of POAG. 

Some of the differently expressed miRNAs we found have already been reported in other studies as indicative of retinal damage. For instance, the differently expressed miRNAs in our study, miR-200a-3p and miR-200b-3p (S-POAG vs. M-POAG, POAG vs. cataract), were found to be critical regulators of POAG induced by the mutation of the OPTN (E50K) gene, which participates in neurodegeneration by induction of apoptosis of retinal ganglion cells in transgenic mice models and progressive retinal degeneration exclusively in the peripheral region of the retinas.⁴⁶ Drewry et al.⁴⁷ also found that miR-200a-3p, miR-200b-3p, miR-200c-3p, miR-141-3p, and miR-429 were differently expressed in normal human ciliary body, cornea, and trabecular meshwork through next-generation sequencing. These five miRNAs were all validated to target the ZFPM2 gene, also known as FOG2, which was glaucoma related.⁴⁸ ZFPM2 was observed to be upregulated after injury to the optic nerve⁴⁹ and may have a role in ocular development.⁵⁰ The differently expressed miRNA in our study, miR-16 (S-POAG vs. M-POAG, POAG vs. cataract), was significantly dysregulated in the glaucomatous retina of rat eyes with advanced nerve damage induced by elevated IOP.⁵¹ Moreover, hsa-miR-183-5p and hsa-miR-96-5p were significantly differently expressed between POAG and cataract eyes in our study. Stephen et al.⁵² found that the knockout of miR-183/96/182 cluster in mice results in early-onset and progressive synaptic defects of the photoreceptors, leading to abnormalities of scotopic and photopic electroretinograms, with decreased b-wave amplitude as the primary defect and progressive retinal degeneration. In addition, inactivation of the miR-183/96/182 cluster resulted in global changes in retinal gene expression, and these changes play an important role in postnatal functional differentiation and synaptic connectivity of photoreceptors. All these studies concerning specific miRNA functions in vivo support our results of differential miRNA profiles in POAG eyes. However, some of the differently expressed miRNAs between the POAG and cataract groups may be related to cataract development rather than POAG. For example, hsa-miR-125b-5p, which was differently expressed when comparing the POAG group with the cataract group in our study, was reported to inhibit lens epithelial cell apoptosis by targeting gene p53 in age-related cataracts.⁵³ 

Genes usually interact with each other to play roles in certain biologic functions. We performed pathway enrichment analysis for target genes of altered miRNA expression based on the KEGG database.²⁵ The observed increased expression of miR-16-5p, miR-205-5p, miR-206, miR-200a-3p, and miR-200b-3p in S-POAG compared with M-POAG is predicted to regulate thiamine metabolism, purine metabolism, and transcriptional misregulation in cancer. Thiamine (or vitamin B1) deficiency is relevant to various neuropathies. Neurologic features of thiamine deficiency include central manifestations, such as Wernicke–Korsakoff encephalopathy and optic neuropathy,⁵⁴ and peripheral manifestations such as axonal peripheral neuropathy.^55,56 Therefore, the dysregulation of these miRNAs related to thiamine metabolism may play an important role in glaucomatous optic neuropathy. Purine metabolism has also been implicated in various neuropathies and neuronal differentiation.⁵⁷ For instance, inborn errors of purine metabolism can cause a series of pediatric neurologic syndromes.⁵⁸ The levels of purine, including adenosine, inosine, uric acid, and xanthine, may be useful for monitoring the progression of diabetic retinopathy and evaluating the treatment.⁵⁹ Our previous study found that the balance of purines might determine the life or death of retinal ganglion cells, as A3 adenosine receptors prevent loss following P2X7 receptor stimulation.⁶⁰ Therefore, regulating these miRNAs related to purine metabolism may find new therapy strategies in glaucomatous optic neuropathy. 

However, the predicted KEGG pathway of dysregulated miRNAs when comparing the POAG group with the cataract group may be more relevant to cataract pathogeneses. The top three predicted pathways were glycosaminoglycan biosynthesis of chondroitin sulfate and dermatan sulfate, glycosaminoglycan biosynthesis of heparan sulfate and heparin, and longevity regulating pathway. The pathway of glycosaminoglycan biosynthesis has long been reported to be implicated in cataract pathogeneses.^61,62 When comparing cataract with POAG patients, the interest of cataract on aqueous miRNA expression cannot be ignored, so part of the dysregulated miRNAs found in previous studies may be related to cataract rather than POAG. 

To get more comprehensive information, we compared our results with other studies with a similar background. Previous studies of Jayaram et al.⁵¹ and our group⁶³ on miRNA expression in the rat retina during acute ocular hypertension, using miRNA microarrays or panels, both found that the dysregulated miRNAs in retina are involved in pathways of inflammation and regulation of apoptosis. These findings indicated that acute glaucoma might induce neuropathy mainly by inflammatory responses, a different mechanism compared with POAG. In addition, Romano et al.⁶⁴ predicted miRNAs that target the reported genes from previous gene-association studies in glaucoma, age-related macular degeneration, and Alzheimer's disease, and found nine miRNAs common to the three neurodegenerative diseases. They predicted that those nine miRNAs might regulate pathways involved in inflammation. However, we can see in their study that the functions of the reported genes themselves are not definite in these three diseases, then the miRNAs are predicted based on these related genes, rather than being found directly in the diseases, and related pathways were finally predicted according to the predicted miRNAs. Therefore, we think that the results in their study may provide a hint, but should be treated with caution and further studies are urgently needed. 

Ghanbari et al.⁶⁵ used data from the recent genome-wide association studies to investigate the association of miRNAs with POAG endophenotypes. They found variants in the miR-612 precursor and in the miR-4707 seed region that are significantly associated with vertical cup-to-disc ratio and cup area. The first variant was demonstrated to increase the miR-612 expression, and the second variant reduced the binding of miR-4707-3p to CARD10, a glaucoma-related gene. These two miRNAs were not found differently expressed in our study, and we speculate the differential expression of miRNAs found in our study might not result from miRNA-related genetic variants. 

A previous study²⁰ comparing miRNA profiles using PCR arrays in AH between POAG and cataract patients revealed that pathways in cancer, chemokine signaling, and Wnt signaling might be regulated by the differently expressed miRNAs. Pathways in cancer are also predicted to be regulated by the differently expressed miRNAs (S-POAG vs. M-POAG) in our study. Another study on the AH of glaucoma patients by using miRNA microarrays discovered several dysregulated miRNAs relating broad pathways, such as RNA posttranscriptional modification, developmental disorders, and hereditary disorders.¹² However, the subjects of the glaucoma group enrolled in that study contained not only POAG, but also pseudoexfoliation glaucoma and primary angle closure glaucoma. In addition, the control group had patients of both cataract and epiretinal membrane patients. There were no overlaps of dysregulated miRNAs and predicted pathways between their results and ours, as expected, which may be due to the large differences in inclusion. 

These informatics analyses in our study and others provide us with specific research directions, but further in vitro and in vivo experiments should be done to prove these hypotheses. This study is limited by the relatively small sample size. In addition, the relatively low RNA yield of each AH sample led to insufficient data output of one-third of the samples (6 of 18 samples), which had to be ruled out for further analysis. The impact of topical medications to lower IOP on miRNA expression within AH has yet to be determined, and future studies can enroll patients with different use of topical glaucoma medications. Moreover, follow-up POAG patients with or without progression of visual field defects can be included to further study the underlying mechanism of different resistibility to IOP in the future. 

To our knowledge, this is the first report on the comprehensive human aqueous miRNA expression profile of POAG patients with varying degrees of visual field damage by using next-generation sequencing. Several dysregulated miRNAs and predicted pathways of thiamine metabolism and purine metabolism may play roles in the development of glaucomatous optic neuropathy. Further insights in understanding their roles in glaucoma may help design new therapeutic strategies. 

Acknowledgments

Supported by the National Natural Science Foundation of China (81670847, 81600728), and the Cultivation Projects for Young Teaching Staff at Sun Yat-sen University (17ykpy75). 

Disclosure: Y. Liu, None; Y. Chen, None; Y. Wang, None; X. Zhang, None; K. Gao, None; S. Chen, None; X. Zhang, None