Analysis of oral microbiome in glaucoma patients using machine learning prediction models

ABSTRACT

Purpose: The microbiome is considered an environmental factor that contributes to the progression of several neurodegenerative diseases. However, the association between microbiome and glaucoma remains unclear. This study investigated the features of the oral microbiome in patients with glaucoma and analyzed the microbiome biomarker candidates using a machine learning approach to predict the severity of glaucoma.

Methods: The taxonomic composition of the oral microbiome was obtained using 16S rRNA gene sequencing, operational taxonomic unit analysis, and diversity analysis. The differentially expressed gene (DEG) analysis was performed to determine the taxonomic differences between the microbiomes of patients with glaucoma and the control participants. Multinomial logistic regression and association rule mining analysis using machine learning were performed to identify the microbiome biomarker related to glaucoma severity.

Results: DEG analysis of the oral microbiome of patients with glaucoma revealed significant depletion of Lactococcus (P = 3.71e⁻³¹), whereas Faecalibacterium was enriched (P = 9.19e⁻¹⁴). The candidate rules generated from the oral microbiome, including Lactococcus, showed 96% accuracy for association with glaucoma.

Conclusions: Our findings indicate microbiome biomarkers for glaucoma severity with high accuracy. The relatively low oral Lactococcus in the glaucoma population suggests that microbial dysbiosis could play an important role in the pathophysiology of glaucoma.

KEYWORDS:

• Neurodegenerative
• oral microbiome
• glaucoma
• biomarker
• dysbiosis

Introduction

Glaucoma is a progressive optic neuropathy characterized by the degeneration of the retinal ganglion cells and their axons, with a corresponding visual field defect [1]. Intraocular pressure (IOP) is a major risk factor for glaucoma development. A high IOP exerts tremendous stress and strain on the lamina cribrosa, consequently blocking the axoplasmic flow [2]. However, glaucomatous axon loss also occurs in individuals with normal IOP, and the risk of glaucoma development might differ even in individuals with the same IOP; thus, glaucoma should be assessed as a complex neurodegenerative disease. Candidate factors associated with glaucoma pathogenesis include aging, inflammation [3], and immune responses, including heat shock proteins or local inflammatory responses [4–8]. As progressive retinal ganglion cell loss is a pathophysiological feature in patients, glaucoma shares signs and symptoms with other neurodegenerative diseases, including Alzheimer’s disease [9]. A population-based cohort study reported a significant association of glaucoma with Alzheimer’s disease, suggesting a potential link between neurodegenerative brain diseases and retinal neurodegenerative diseases [10,11].

Technological advances in next-generation sequencing (NGS) and bioinformatics tools have identified an association of human microbiota with neurodegenerative diseases and autoimmunity. The altered human microbiome, or microbial dysbiosis, can contribute to the pathogenesis of several diseases [12]. Although the ocular microbiome is related to eye diseases, non-ocular microbiome such as an oral microbiome could influence the pathogenesis of eye disease [3,6,13,14]. The microbiota-gut-brain axis involves a cross-talk between the gut and microbiota; enteric bacteria enters the circulation and crosses the blood-brain barrier, or to penetrates the brain via the vagus nerve [15,16]. Since the oral microbiome is a part of the gut microbiome family and 5% of oral species are found in the intestine, the microbiome-gut-retinal axis may have clinical importance. A previous study that investigated the relationship between microbiota and glaucoma reported high oral bacterial loads, such as Streptococcus species, in patients with glaucoma [7]. To date, only a few studies have investigated the influence of the gut microbiome on glaucoma. Therefore, this study aimed to investigate the differences in the oral microbiome in patients with glaucoma based on the microbiota-gut-retinal axis, and identify the candidate microbiome biomarkers related to glaucoma pathogenesis and severity using machine learning prediction models.

Material and methods

Participants

This multicentric study was conducted by researchers from Veterans Health Service (VHS) Medical Center, Daegu Veterans Hospital, Pusan National University Yangsan Hospital, and Pusan National University Hospital. The study protocol was approved by the Institutional Review Board of each participating center. The study was performed in accordance with the tenets of the Declaration of Helsinki and written informed consent was obtained from all participants before their enrollment.

During the study period, the research team enrolled patients with glaucoma who satisfied the inclusion criteria. Glaucoma was defined as having glaucomatous optic neuropathy, such as rim thinning, notching, retinal nerve fiber layer (RNFL) defect, or glaucomatous visual field defects. The control group was defined as participants with 1) a visual acuity of 20/40 or better; and/or those without 2) a glaucomatous optic nerve head or RNFL defect, visual field defect, or other retinal diseases; or 3) a history of IOP-lowering treatment. Participants with a history of ocular surgery (except uncomplicated cataract surgery), history of ocular trauma, or other diseases affecting the visual field (e.g. diabetic retinopathy) were excluded. Participants with a dental issue, those with probiotics use within the last 1 month, those with upper respiratory infection symptoms within the last two weeks, and those taking antibiotics or using topical antibiotic eye drops within the last one month were additionally excluded.

Sample collection and DNA extraction

The oral microbiome was collected by swabbing the right and left buccal surfaces and under the tongue four times using the Gene kit (Daeil pharm company, Korea) after washing the oral cavity with clean drinking water. Metagenomic DNA was individually extracted from swabs using the GeneAll Exgene^TM Blood/Clinic/Cell SV mini kit (GeneAll, Seoul, Korea) according to the manufacturer’s instructions. The DNA was stored at −80°C following DNA extraction, and the quantity and quality of the isolated DNA were measured using Trinean Dropsense 96 (Unchained Labs, Pleasanton, CA, USA) and PicoGreen (Thermo Fisher Scientific, Waltham, MA, USA).

Sequencing of 16S rRNA gene amplicon

A 16S rRNA sequencing library was generated according to the 16S metagenomic sequencing library preparation protocol (Illumina, San Diego, CA, USA), targeting the V3 and V4 hypervariable regions of the 16S rRNA gene, which were selected based on the results reported previously [17]. The KAPA HiFi HotStart ReadyMix (Kapa Biosystems, Wilmington, MA, USA) and Agencourt AMPure XP system (Beckman Coulter Genomics, Brea, CA, USA) were used for polymerase chain reaction (PCR) amplification and amplicon purification, respectively. The Illumina adapter overhang nucleotide sequences were added to the gene-specific sequences. The following full-length primer sequences were used to target these adapters:

16S Amplicon PCR Forward Primer:

5ʹ TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG

16S 3ʹ

Amplicon PCR Reverse Primer

5ʹ GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC 3ʹ

After magnetic bead-based purification of the PCR products, a second PCR with a limited cycle was performed using primers from the Nextera XT Index Kit (Illumina). Subsequently, the purified PCR products were visualized using gel electrophoresis and quantified with a Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific) on a Qubit 3.0 fluorometer. The samples were pooled and assayed using Agilent 2100 bioanalyzer (Agilent) for pre-sequencing quality analysis. The previously generated sequencing libraries were quantified by qPCR using the CFX96 Real-Time System (Bio-Rad, Hercules, CA, USA). After normalization, the libraries were sequenced on the MiSeq system (Illumina) with 300 bp paired-end reads.

Pre-processing

The adaptor sequences were removed from the original paired-end reads using CutAdapt v1.11. Next, merged reads were generated using the adaptor-trimmed paired-end reads using FLASH v1.2.11. The low-quality merged reads were filtered according to the following criteria: reads that had ≥2 ambiguous nucleotides, reads with average quality score <20, and reads with length <300 bp after trimming of low-quality bases. Finally, potential chimeric reads were removed using UCHIME v4.2.40 [18].

Calculation of Operational Taxonomic Units (OTUs) and diversity analysis

The pre-processed reads from each sample were used to calculate the number of OTUs, which were determined by clustering the sequences from each sample with a 97% sequence identity cut-off using QIIME software (v.1.9.0, www.quiime.org) [19,20]. The taxonomic abundance was counted with Ribosomal Database Project (RDP) Classifier v1.1, with a confidence threshold of 0.8 derived from the pre-processed reads of each sample [21]. The microbial composition in the samples was normalized using the value calculated by a division between the taxonomic abundance count and the number of pre-processed reads for each sample. Consensus sequences were clustered using cd-hit v4.6 using the following parameters – identify, >99%; and coverage, >80% – and then, aligned with the reference sequences from the National Center for Biotechnology Information (NCBI) database using the Megablast algorithm. Finally, taxonomic profiling was performed for the assembled genome using NCBI taxonomy information. The alpha-diversity for observed OTUs and Shannon index within the samples were calculated. Principal component analysis (PCoA) was performed using the result from the beta-diversity for the glaucoma and control groups, as well as for subgroup analysis.

Machine learning-based identification of microbiome biomarker for glaucoma severity

Machine learning-based identification of glaucoma biomarkers consisted of three steps: significant taxon level identification by Random Forest (RF), identification of a significant taxon by Tag Count Comparison (TCC), and rule mining of the microbiome. RF model-based machine learning was performed with the R software (R Foundation for Statistical Computing, Vienna, Austria) using ‘ranger’ package [22]. Within each run, the number of trees (num.tree) was set to build the RF model [23–25]. The number of variables/taxa that could be selected in each splitting node (mtry) was set to a range of one-tenth of the maximum taxon number at each taxonomic level. The Gini index evaluated the importance of each taxon, and the prediction accuracy was determined by $\frac{TruePostive\left(TP\right)+TrueNegative\left(TN\right)}{TP+TN+FalsePositive+FalseNegative}$. Hundred runs were performed for each mtry parameter within each taxonomic hierarchy. We calculated the significant taxa using TCC (http://bioconductor.org/), which provides functions for differential expression analysis of genes (DEG) with normalization and multi-group comparison [25]. With TCC analysis, significant oral microbiome differences were analyzed in the glaucoma group, compared to the control group, and in the primary open-angle glaucoma (POAG) subgroups, compared to the control group.

Before rule mining, a multinomial logistic regression (MLR) was performed to analyze the significance of the identified microbiome biomarker in patients with glaucoma based on the clinical factors, such as age, baseline IOP level, RNFL thickness (assessed using optical coherence tomography [OCT]), and mean deviation value in the Humphrey visual field. Patients aged <60 years comprised the younger age group, and those ≥60 years comprised the older age group. The baseline IOP cut-off was 20 mmHg; normal baseline IOP, <20 mmHg; and high baseline IOP, ≥20 mmHg. Although arbitrary, this value is acceptable for statistical processing since the IOP criterion for normal-tension glaucoma is based on a cut-off of 21 mmHg, and the mean IOP of Asians is 1–2 mmHg lower than the other ethnicities. For glaucoma severity, the glaucoma group was divided into ≥80 μm and <80 μm subgroups based on the RNFL thickness. For functional glaucoma severity, the glaucoma group was stratified according to Anderson’s glaucoma stage criteria: ≥−6.0 dB, −6.0 to −12.0 dB, and <−12.0 dB.

The classification based on the predictive association rules (CPAR) algorithm was used to generate a complete ruleset [24]. This algorithm, which uses information metrics to generate rules, is more suited for bioinformatics applications than conventional support-confidence based measures for market basket data analysis [26]. The Laplace accuracy evaluated the accuracy of the CPAR-derived rules. Using rule r, the Laplace accuracy is defined as below:

$Laplaceaccuracy\left(r\right)=\frac{\left(Ng+1\right)}{\left(Ntotal+m\right)}$

where m, number of target groups; N_total, total number of examples that satisfy the body of the rules; and Ng, number of examples that belong to the predicted target group. Only the POAG subgroup was used for subgroup analysis using MLR and CPAR because of insufficient sample size in secondary glaucoma (SG) group.

Data analyses

All the data analyses were performed using the R Statistical Package, Version 3.6.2 (R Foundation for Statistical Computing, Vienna, Austria). Statistical significance was set at P < 0.05. For multiple comparison tests in DEG analysis, the false discovery rate was controlled using the Benjamini–Hochberg step-up procedure [27,28].

Results

Patient characteristics and metagenomic data

Ninety-six patients with glaucoma and 25 controls were recruited from March 2019 to December 2019. The glaucoma group had 62 POAG and 34 SG patients (SG cases were mostly uveitic glaucoma). There was no difference in age and sex between patients with glaucoma and the normal participants (Table 1). In the glaucoma groups, the baseline IOP was 18.49 mmHg, which was higher than that (13.76 mmHg) of the normal participants. RNFL thickness, a structural indicator of glaucoma, was low in glaucoma patients, and the visual field index, a functional indicator, was moderately advanced glaucoma. The number of the average lead counts/sample of the control group was higher than that of the glaucoma group (P = 0.038). Results of machine learning using the RF method showed that genus level identification was significant in 193 mtry (accuracy, 0.5714–0.7429) than those identified based on family, order, and phylum, and Lactococcus had the highest Gini indices (Figure 1, Table 2).

Table 1. Demographic and baseline characteristics of glaucoma subjects

Variables	Glaucoma n = 96	Control subject n = 25	P-value*
Age, years	61.48 ± 13.08	62.68 ± 1.73	0.385
Male/Female	66/30	18/7	0.753
BCVA, logMAR	0.21 ± 0.32	0.00 ± 0.20	<0.0001
Diagnosis, n (%)
POAG/Secondary glaucoma	62 (64.58)/34 (35.42)	-
Glaucoma medication period (yr)	2.74 ± 3.09	-
Baseline IOP, mmHg	18.49 ± 6.60	13.76 ± 1.01	<0.0001
RNFL thickness, µm	73.44 ± 14.52	105.36 ± 5.62	<0.0001
Visual field (MD), dB	−8.93 ± 8.44	0.03 ± 0.17	<0.0001
Visual field (VFI), %	76.97 ± 26.52	99.84 ± 0.47	<0.0001
Number of lead counts/sample	77,140.05 ± 23,087.21	83,310.64 ± 8,812.47	0.038

P-value* was calculated by independent T-test for continuous variables, Chi-square test was performed in case of numeric variable

IOP: intraocular pressure, BCVA: best-corrected visual acuity, logMAR: logarithm of the minimum angle of resolution, MD: mean deviation. VFI: visual field index, RNFL: retinal nerve fiber layer, POAG: primary open angle glaucoma

Table 2. Gini index for accuracy condition in genus levels using random forest machine learning

Taxon Level	Ntree	Ratio	mtry	Times	Accuracy	95% CI	P-Value^†
*genus	200	0.1	193	100	0.5714 ~ 0.7429	0.3935 ~ 0.8751	2.81e^−08 ~ 0.0003783
genus	200	0.2	385	100	0.5714 ~ 0.7143	0.3935 ~ 0.8536	1.869e^−07 ~ 0.0003783
genus	200	0.3	578	100	0.6 ~ 0.7143	0.4211 ~ 0.8536	1.869e^−07 ~ 0.0001039
genus	200	0.4	771	100	0.6 ~ 0.7143	0.4211 ~ 0.8536	1.869e^−07 ~ 0.0001039
genus	200	0.5	964	100	0.5714 ~ 0.7143	0.3935 ~ 0.8536	1.869e^−07 ~ 0.0003783
genus	200	0.6	1156	100	0.6 ~ 0.7143	0.4211 ~ 0.8536	1.869e^−07 ~ 0.0001039
genus	200	0.7	1349	100	0.6 ~ 0.6857	0.4211 ~ 0.8315	1.089e^−06 ~ 0.0001039
genus	200	0.8	1542	100	0.5714 ~ 0.7143	0.3935 ~ 0.8536	1.869e^−07 ~ 0.0003783
genus	200	0.9	1734	100	0.6 ~ 0.7143	0.4211 ~ 0.8536	1.869e^−07 ~ 0.0001039
genus	200	1	1927	100	0.5714 ~ 0.7429	0.3935 ~ 0.8751	2.81e^−08 ~ 0.0003783

*Best accuracy condition in genus level, P-Value^† was calculated by machine learning using random forest model

95% CI: confidential Interval

Figure 1. Metagenome analysis using machine learning and statistical models

a: Schematic of machine learning for the study of algorithms and statistical models [23–25] TCC, tag count comparison; CPAR, Classification based on predictive association rules

b: Using Random Forest method, the significant taxa and genus levels in the oral microbiome were identified.

c: The top 10 species of the oral microbiome as identified from the Ribosomal Database Project (RDP) Data Base; Lactococcus was the most significant species in glaucoma patients and control subjects.

Diversity analysis in glaucoma patients and control subjects

On comparisons of the alpha-diversity of oral microbiome between glaucomatous and control subjects, the observed OTUs (P = 0.003 and P = 0.011, respectively) and the Shannon index (P = 6.1e⁻⁶ and P = 3.6e⁻⁵, respectively; Figure 2) were significantly higher in the control participants than that in the glaucoma group and glaucoma subgroups. PCoA of beta-diversity for the oral microbiome showed that the microbial clusters in the glaucoma and control groups differed in their taxonomy (Figure 3). On subgroup analysis, the microbial clusters in POAG and SG subgroups were similar, but that in the control group was different. The beta-distance calculated from the RDP classifier showed significant differences, but the beta-distance from the NCBI method was not different between the glaucoma and control groups (P < 2.2e⁻¹⁶, P = 0.52, respectively; Figure 4), whereas the beta-distance of subgroup analysis showed a significant difference in the oral microbiome (P < 2.2e⁻¹⁶).

Figure 2. Comparison of alpha-diversity in glaucoma patients as compared with the control subjects

(a) Observed OTUs was higher in the control group than in the glaucoma (P = 0.003). It was also higher in the control group than in the two subgroups (P = 0.011)

(b) Shannon index was higher in the control group than in the glaucoma group (P = 6.1e⁻⁰⁶). It was also higher in the control group than in the two subgroups (P = 3.6e⁻⁰⁵). POAG, primary open-angle glaucoma; SG, secondary glaucoma

Figure 3. Differences in the oral microbiome of glaucoma patients compared with control subjects using principal component analysis for beta-diversity

Beta-diversity was calculated in two ways – RDP and NCBI library. PCoA shows the different clusters in glaucoma vs. control groups, and POAG and SG subgroups vs. control group.

RDP, Ribosomal Database Project Classifier; NCBI, National Center for Biotechnology Information; POAG, Primary open-angle glaucoma; SG, secondary glaucoma

Figure 4. Beta-Diversity differences in the microbiome of glaucoma patients and control subjects

Beta-diversity (estimated by RDP method) was significantly higher in glaucoma group than in control groups (P < 2.2e⁻¹⁶); however, the beta-diversity using NCBI was not different between control and glaucoma groups (P = 0.52).

RDP, Ribosomal Database Project Classifier; NCBI, National Center for Biotechnology Information; POAG, Primary open-angle glaucoma; SG, secondary glaucoma

DEG analysis in glaucoma group and subgroup

For TCC-based taxon significance analysis, the MA plot revealed the different quantities of some species and the DEG analysis showed various taxa such as Lactococcus and Atopobium were significantly depleted (P = 3.71e⁻³¹, P = 6.69e⁻²¹, respectively), whereas Faecalibacterium was enriched in the glaucoma group (P = 9.19e⁻¹⁴) compared with control participants (Figure 5, Supplementary Table S1). On subgroup analysis, the volcano plot of the oral microbiome showed that Lactococcus, Candidatus Pelagibacter, and Isobaculum were significantly depleted P = 1.67e⁻²⁰, P = 9.37e⁻²⁰, and P = 1.43e⁻¹⁸, respectively, whereas Faecalibacterium and Lachnospiracea incertae sedis were enriched (P = 7.21e⁻¹⁶ and P = 6.02e⁻¹², respectively) in the POAG subgroup compared with the control group (Figure 6 and Supplementary Table S2).

Figure 5. MA and volcano plots generated using analysis with tag count comparison for oral microbiome of glaucoma patients and control subjects

a: MA plot of glaucoma vs control subjects (DE number of the upper panels 166 vs. that of the lower panels 171; P < 0.05)

b: Volcano plot of DEG in glaucoma vs control subjects Lactococcus and Candidatus Pelagibacter were downregulated in the glaucoma group DEG, differential expression of genes

Figure 6. MA and volcano plots generated using subgroup analysis with tag count comparison for oral microbiome in the POAG subgroups compared with the control group

a: MA plot of POAG vs control subjects (DE number of the upper panels 175 vs. that of the lower panels 142; P < 0.05)

b: Volcano plot of DEG in POAG vs control subjects DEG, differential expression of genes; POAG, primary open-angle glaucoma

MLR and association rules mining analysis for biomarkers of glaucoma

The Lactococcus, Candidatus Pelagibacter, and Atopobium were the significant taxa for factors age and baseline IOP. In addition, they were significantly associated with glaucoma severity based on RNFL thickness and visual field index (Table 3). Association rule generated by CPAR showed two ‘glaucoma’ rules and four ‘normal’ rules generated with a high accuracy in the 0.9–0.96 range (Table 4). ‘Glaucoma’ rule 1 is the most accurate, which states that if the microbial composition is Atopobium ≤0.347, Candidatus Pelagibacter ≤0.001, Endobacter ≤0.015, and Lactococcus ≤0.013, the sample patients is glaucomatous (accuracy, 96.0%). Similarly, ‘glaucoma’ rule 2 states that if the composition is Lactococcus, ≤0.013, the patient is ‘glaucoma’ (accuracy, 90.0%). Contrarily, ‘normal’ rule 1 states that if the composition is Candidatus Pelagibacter, 0.001–0.014 and Lactococcus, 0.013–0.889, the patient is ‘normal’ (accuracy, 92.0%).

Table 3. Multinomial logistic regression analysis of factors related taxon in glaucoma group compared with control subjects

	coefficient	SD	Odd ratio	P-value	coefficient	SD	Odd ratio	P-value*
Taxa	Age < 60				Age ≥ 60
Lactococcus	−210.481	55.386	3.88e^−92	0.0001	−172.977	45.470	7.53e^−76	0.0001
Endobacter	0.223	1.750	1.25e^+00	0.8984	−27.032	15.330	1.82e^−12	0.0778
Rikenella	89.032	237.423	4.63e^+38	0.7077	−80.714	338.227	8.83e^−36	0.8114
Candidatus Pelagibacter	−560.729	205.488	3.01e^−244	0.0064	−921.479	259.819	0.00e^+00	0.0004
Atopobium	−3.561	1.322	0.02842	0.0071	−4.122	1.244	0.01620	0.0009
	Baseline IOP < 20 mmHg				Baseline IOP ≥ 20 mmHg
Lactococcus	−177.477	45.896	8.37e^−78	0.0001	−191.253	52.588	8.70e^−84	0.0003
Endobacter	−13.841	12.697	9.75e^−07	0.2757	0.094	1.730	1.10e^+00	0.9567
Rikenella	14.974	47.738	3.18e^+06	0.7538	33.475	46.674	3.45e^+14	0.4732
Candidatus Pelagibacter	−509.789	171.543	3.99e^−222	0.0030	−319.291	145.031	2.16e^−139	0.0277
Atopobium	−4.134	1.252	1.60e^−02	0.0010	−3.561	1.310	2.84e^−02	0.0065
	RNFL thickness ≥ 80 μm				RNFL thickness < 80 μm
Lactococcus	−207.229	53.820	1.00e^−90	0.0001	−173.639	45.928	3.89e^−76	0.0002
Endobacter	−13.957	15.484	8.68e^−07	0.3674	−0.536	1.711	5.85e^−01	0.7542
Rikenella	−38.900	295.676	1.28e^−17	0.8953	61.128	142.806	3.53e^+26	0.6686
Candidatus Pelagibacter	−679.649	250.852	6.80e^−296	0.0067	−782.464	226.057	0.00e^+00	0.0005
Atopobium	−4.731	1.622	8.81e^−03	0.0035	−3.545	1.118	2.89e^−02	0.0015
									coefficient	SD	Odd ratio	P-value
	Early glaucoma				Moderate glaucoma				Advanced glaucoma
Lactococcus	−155.716	46.085	2.36e^−68	0.0007	−231.783	65.006	2.18e^−101	0.0004	−218.896	60.662	8.60e^−96	0.0003
Endobacter	−25.875	17.045	5.79e^−12	0.1290	−6.534	10.847	1.45e^−03	0.5469	0.488	1.756	1.63e^+00	0.7813
Rikenella	−35.096	155.035	5.73e^−16	0.8209	32.121	63.202	8.91e^+13	0.6113	43.848	62.576	1.10e^+19	0.4835
Candidatus Pelagibacter	−1695.298	445.654	0.00e^+00	0.0001	−1370.249	457.311	0	0.0027	−1342.329	430.780	0	0.0018
Atopobium	−5.753	1.710	3.17e^−03	0.0008	−15.354	5.126	2.15e^−07	0.0027	−2.135	0.941	1.18e^−01	0.0232

P-value* was calculated using multinomial logistic regression analysis

IOP: intraocular pressure, RNFL: retinal nerve fiber layer, SD: standard deviation

Table 4. Association rules generated by association rule mining for glaucoma using oral microbiome data

Glaucoma rules
Rule 1	Atopobium = [1] Candidatus Pelagibacter = [1] Endobacter = [1] Lactococcus = [1] → Glaucoma 0.96
Rule 2	Lactococcus = [1] → Glaucoma 0.9
Normal rules
Rule 1	Candidatus Pelagibacter = [2] Lactococcus = [2] → Normal 0.92
Rule 2	Lactococcus = [2] → Normal 0.85
Rule 3	Candidatus Pelagibacter = [2] → Normal 0.7
Rule 4	Atopobium = [2] → Normal 0.65
Rule explanations (Rule 1)
Original rule	Atopobium = [1] Candidatus Pelagibacter = [1] Endobacter = [1] Lactococcus = [1] → Glaucoma 0.96
Rule explanations	If the compositions of Atopobium ≤ 0.347 Candidatus Pelagibacter ≤ 0.001 Endobacter ≤ 0.015 Lactococcus ≤ 0.013; the sample is Glaucoma with 96.0% accuracy
Applied category for rule generation
Category	Atopobium	Candidatus Pelagibacter	Endobacter	Lactococcus	Rikenella
[1]	≤ 0.347	≤ 0.001	≤ 0.015	≤ 0.013	≤ 0.001
[2]	0.347–2.908	0.001–0.014	0.015–1.106	0.013–0.889	0.001–0.159

CPAR (classification based on predictive Association Rules) algorithm used [24].

MLR and association rules mining analysis for biomarkers of POAG

In patients with POAG aged <60 years, Lactococcus was the significant taxon, whereas Lactococcus, Candidatus Pelagibacter, and Atopobium were the significant taxa in patients with POAG aged ≥60 years (Table 5). In the POAG group with baseline IOP <20 mmHg, Lactococcus, Candidatus Pelagibacter, and Atopobium were the significant taxa, whereas Lactococcus and Atopobium were the significant taxa in the POAG group with baseline IOP≥20 mmHg (Table 5). According to glaucoma severity based on RNFL thickness and visual field index in the POAG group, Lactococcus, Candidatus Pelagibacter, and Atopobium were the significant taxa. CPAR-derived association rule shows one ‘POAG’ rule and four ‘normal’ rules derived with high accuracy (Table 6). The ‘POAG’ rule 1 states that if the composition is Atopobium, ≤0.433; Candidatus Pelagibacter, ≤0.001; Endobacter, ≤0.015; and Lactococcus, ≤0.015, the patient has POAG (accuracy, 93.0%). The ‘normal’ rule 1 states that if the composition is Candidatus Pelagibacter, 0.001–0.014, the patients is ‘normal’ (accuracy, 93.0%).

Table 5. Multinomial logistic regression analysis of factors related taxon in primary open-angle glaucoma

	coefficient	SD	Odd ratio	P-value	coefficient	SD	Odd ratio	P-value*
	Age < 60				Age ≥ 60
Lactococcus	−175.539	65.322	5.81e^−77	0.0072	−157.432	45.664	4.25e^−69	0.0006
Endobacter	−10.362	14.712	3.16e^−05	0.4813	−57.292	24.940	1.31e^−25	0.0216
Rikenella	95.277	195.573	2.39e^+41	0.6261	−97.172	373.910	6.29e^−43	0.7950
Candidatus Pelagibacter	−328.189	215.986	2.95e^−143	0.1286	−637.095	220.119	2.06e^−277	0.0038
Atopobium	−3.874	2.006	0.0208	0.0534	−3.826	1.297	0.02179	0.0032
	Baseline IOP < 20 mmHg				Baseline IOP ≥ 20 mmHg
Lactococcus	−159.538	46.999	5.17e^−70	0.0007	−163.518	56.609	9.66e^−72	0.0039
Endobacter	−27.835	15.960	8.15e^−13	0.0812	−52.336	33.014	1.87e^−23	0.1129
Rikenella	18.174	46.865	7.81e^+07	0.6982	39.064	45.859	9.23e^+16	0.3943
Candidatus Pelagibacter	−602.133	230.776	3.14e^−262	0.0091	−435.431	232.158	7.85e^−190	0.0607
Atopobium	−5.446	1.773	4.31e^−03	0.0021	−2.531	1.203	7.96e^−02	0.0353
	RNFL thickness ≥ 80 μm				RNFL thickness < 80 μm
Lactococcus	−204.775	58.817	1.17e^−89	0.0005	−140.536	46.567	9.24e^−62	0.0025
Endobacter	−44.285	25.725	5.85e^−20	0.0852	−28.076	16.803	6.41e^−13	0.0947
Rikenella	−41.052	327.684	1.48e^−18	0.9003	64.161	125.357	7.32e^+27	0.6088
Candidatus Pelagibacter	−323.115	164.296	4.71e^−141	0.0492	−401.957	172.983	2.71e^−175	0.0201
Atopobium	−4.243	1.712	1.44e^−02	0.0132	−3.599	1.341	2.74e^−02	0.0073
									coefficient	SD	Odd ratio	P-value
	Early glaucoma				Moderate glaucoma				Advanced glaucoma
Lactococcus	−125.291	45.677	3.86e^−55	0.0061	−229.962	79.424	1.35e^−100	0.0038	−216.259	72.337	1.20e^−94	0.0028
Endobacter	−31.014	18.134	3.39e^−14	0.0872	−29.294	25.434	1.89e^−13	0.2494	−45.306	32.995	2.11e^−20	0.1697
Rikenella	−28.362	105.738	4.82e^−13	0.7885	29.192	42.020	4.76e^+12	0.4872	39.907	41.176	2.14e^+17	0.3325
Candidatus Pelagibacter	−1122.831	394.907	0.00e^+00	0.0045	−854.987	435.199	0	0.0495	−777.340	372.864	0	0.0371
Atopobium	−5.434	1.875	4.36e^−03	0.0038	−19.236	7.884	4.42e^−09	0.0147	−1.926	1.071	1.46e^−01	0.0722

P-value* was calculated using multinomial logistic regression analysis

IOP: intraocular pressure, RNFL: retinal nerve fiber layer, SD: standard deviation

Table 6. Association rules generated by association rule mining for primary open-angle glaucoma

POAG rules
Rule 1	Atopobium = [1] Candidatus Pelagibacter = [1] Endobacter = [1] Lactococcus = [1] → POAG 0.93
Normal rules
Rule 1	Candidatus Pelagibacter = [2] → Normal 0.93
Rule 2	Lactococcus = [2] → Normal 0.8
Rule 3	Atopobium = [2] → Normal 0.73
Rule 4	Endobacter = [2] → Normal 0.53
Rule explanation (Rule 1)
Original rule	Atopobium = [1] Candidatus Pelagibacter = [1] Endobacter = [1] Lactococcus = [1] → POAG 0.93
Rule explanations	If the compositions of Atopobium ≤ 0.433 Candidatus Pelagibacter ≤ 0.001 Endobacter ≤ 0.015 Lactococcus ≤ 0.015; the sample is POAG with 93.0% accuracy
Applied category for rule generation
Category	Atopobium	Candidatus Pelagibacter	Endobacter	Lactococcus	Rikenella
[1]	≤ 0.433	≤ 0.001	≤ 0.015	≤ 0.015	≤ 0.001
[2]	0.433–2.908	0.001–0.014	0.015–0.475	0.015 − 0.889	0.001–0.159

CPAR (classification based on predictive Association Rules) algorithm used [24].

Discussion

In this study, we demonstrated that the oral microbiome in patients with glaucoma was different from that of the elderly without glaucoma, and identified candidate microbiome biomarkers for glaucoma severity with high accuracy using machine learning. Our study findings are consistent with those of previous studies in that we suggested a possible relationship between microbiome and glaucoma [3–5,7,29]. This relationship can be interpreted in several ways. As an immunity problem, the microbiome acts as a possible antigen, which may activate T-helper cell development and increased the production of antibodies. The trabecular meshwork cells can recognize the bacterial antigen, a pathogen-associated molecular pattern, such as the toll-like receptor, to generate imbalanced signaling of inflammation [30]. Although this mechanism is not fully understood, it is considered that this cascade results in the anterior chamber-associated immune deviation, which induces immune cell activation as well as homing to the ocular region [29].

A previous study on the association of glaucoma and the oral microbiome has shown that altered commensal microbiome induces changes in cytokine signaling and complement activation [7]. In addition, one possible mechanism by which the microbiome could aid the management of glaucoma is by modulating brain-derived neurotrophic factor levels, as shown in an animal model [31]. From the perspective of neuroprotection, microbiome research suggests that the pathogenesis of glaucoma may be identified. According to the study by Astafurov et al., oral bacterial load, such as of Streptococcus spp., in glaucoma cases was increased compared with the control [7]. However, in our study, the bacterial diversity was not significantly high in glaucoma patients, but there was a significant difference in the composition of the bacteria, such as Faecalibacterium, which was overexpressed in patients with glaucoma. Faecalibacterium modulated the butyrate metabolite in the gut, which was known to be related to anti-inflammatory effect [32]. Mental health disease such as major depression negatively correlated with Faecalibacterium, which was different in Alzheimer's disease and Parkinson’s disease [32]. However, it should be judged in consideration of the overall microbiome composition. Additionally, another study on the intestinal microbiota profile and POAG showed that Megamonas and Bacteroides load was decreased in POAG [33]. The study found negative correlations between visual acuity and Blautia and between the mean deviation of visual field and Faecalibacterium, and a positive correlation between RNFL thickness and Streptococcus. This study also showed that the gut microbiome might differ between glaucoma and control subjects, which is consistent with the results of our study, although taxa expression level was different due to differences in the targeted gut microbiome.

The reduced Lactococcus load in the oral microbiome suggests that dysbiosis in the microbiota could play a role in the pathophysiology of glaucoma and can be used as a surrogate marker for glaucoma. It may be due to competition with other taxa in the oral microbiome or decreased antimicrobial activity of H₂O₂ produced by downregulation of the lactic acid bacteria. Gastrointestinal permeability is altered depending on diet and microbiota. In this regard, diet enriched with a high sugar content may cause a decrease in the number of beneficial microbes [34], which may directly or indirectly affect the development of glaucoma. The carbohydrate ingestion in POAG patients showed systemic autonomic dysregulation [35], whereas ketogenic diets and modified gut microbiome and permeability may play a role in neuroprotection in glaucoma patients [36,37].

Several similarities may exist in glaucoma and neurodegenerative disease, as the optic nerve is a part of the central nerve system [38,39]. Hence, the microbiome study of Parkinson’s disease and Alzheimer’s disease may provide new insights into glaucoma research, suggesting that microbial dysbiosis can induce or exacerbate the disease via toxic effects from direct invasion or epigenetic changes [40]. One example of commensal pathogenic bacteria is H. Pylori, which has been associated with the risk of several neurodegenerative diseases, including Parkinson’s disease, Alzheimer’s disease, and POAG [41–43].

A major strength of our study is the inclusion of a relatively large number of glaucoma samples from multiple centers using metagenomic analysis. In addition, our study used recently developed bioinformatics technology and machine learning for analysis. However, there are several limitations to this study. First, it was difficult to determine causality in this study; we could not determine whether the microbiome biomarker alone is a bystander or culprit. This study is worthwhile as an early study to suggest biomarkers that could be related to the identification of the severity of glaucoma. Additional verification is needed using cell-based and animal model studies for mechanisms. The validation study by PCR would provide for more reliable results, although the NGS technology has improved. Another limitation was the small sample size of the control group compared to that of glaucoma patients. A large number of control subjects (over 50 control subjects) and other replication cohorts will be required for further validation. The sample size for metagenomics research is usually 30; as such, our sample size was not very small. However, larger-scale prospective studies in the future are necessary to support our study results. This study was conducted as a multi-center study, and it could potentially be biased by the inclusion of people living at various locations of Korea. The microbiome may vary according to different time points and physical conditions. The repeatability of microbiome results can be established by obtaining different samples from the same individuals. Finally, the taxa preferentially selected for analysis were limited by the computing resources for machine learning modeling.

In conclusion, the oral microbiome in patients with glaucoma was different from that in control participants that indicated disturbing environmental influences on the pathogenesis of glaucoma via the microbiota-gut-retina axis. Our findings highlight the identification of microbiome biomarkers for glaucoma severity with high accuracy. The reduced Lactococcus count in the oral cavity suggests that dysbiosis in microbiota could play a role in the pathogenesis of glaucoma.

Ethical approval

The study protocol was approved by the Institutional Review Board (IRB) of each participating center—Veterans Health Service Medical Center, Korea (IRB No. 2018-12-015 and IRB No. 2020-06-025) Pusan National University Hospital (IRB No. H-1904-001-077); Daegu Veterans Health Service Medical Center (IRB No. 2019-001); and Pusan National University Yangsan Hospital (IRB No.05-2019-022).

Informed consent

Written informed consent was obtained from all subjects before their enrollment in the study.

Financial disclosure

This study was supported by a VHS Medical Center Research Grant, Republic of Korea (grant number: VHSMC 19043)

Conflicts of interest

The authors have no conflicts of interest to declare.

Acknowledgments

We thank Theragen Bio for providing technical support for genome sequencing and metagenome analysis.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed here.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This study was supported by VHS Medical Center Research Grant (grant number: VHSMC 19043), Republic of Korea.