Introduction

Glaucoma is a chronic and progressive disease characterized by damage to the retinal nerve fiber layer (RNFL) and associated with morphological changes to the optic nerve1. The traditional view is that high levels of intraocular pressure (IOP) are not tolerated by the retina and optic nerve resulting in a deformed lamina cribrosa, causing axonal damage and consequent apoptotic death of retinal ganglion cells (RGCs). However, over the last decade, research has shown that glaucoma might not be solely an eye disease, but a complex neurodegenerative disease also involving processes within the brain tissue2,3,4. While the nature of the neurodegenerative component to glaucoma remains to be resolved, it is clear that damage in glaucoma is not restricted to the axons that form the optic nerve. In addition, numerous magnetic resonance imaging (MRI) studies have shown that the pre-geniculate (optic tract and chiasm), geniculate, and post-geniculate (the optic radiation (OR) and primary visual cortex) structures are affected as well in glaucoma, at least in later stages of the disease5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20. Moreover, changes in white matter (WM) tracts beyond the primary visual pathways such as the anterior thalamic radiation, corticospinal tract, superior longitudinal fasciculus, and forceps major17,18,21,22,23,24 and a network reorganization of the glaucomatous brain has been shown as well25. In line with this idea, it has also been proposed that glaucoma and Alzheimer’s Disease (AD) might share a common pathogenesis26. Indeed, both diseases are progressive, chronic, age-related diseases and cause RGC degeneration and irreversible neuronal cell loss. Furthermore, research suggests a genetic link between glaucoma and AD27,28,29. More recently, dysfunction of the glymphatic system (i.e. a cerebrospinal fluid transport system that facilitates clearance of neurotoxic molecules) was proposed as a novel and missing link between AD and glaucoma30,31,32,33,34,35.

At present, the most parsimonious explanation for changes to the brain tissue in glaucoma remains that these are a consequence of the decreased visual input due to retinal and optic nerve damage19,36. Yet, a complete picture of the alterations of the WM tissue in glaucoma is missing. Understanding the degree to which alterations in the brain tissue in glaucoma go beyond the early visual system, and the degree to which the effect of glaucoma outside of the visual system is similar to that of other diseases also reducing visual input, can help to clarify the pathophysiology of glaucoma. Answering this question is critical to advance our understanding of the disease and may, ultimately, inform future therapeutic approaches.

Our current study tested whether WM microstructure (WMM) alterations are unique to the glaucomatous brain or depend primarily on decreased visual input. To do so, the study used reproducible cloud computing methods to measure WMM in a much larger set of tracts37,38 than previously considered. We compared WMM between two independent sets of glaucoma patients (using two independently acquired and demographically different glaucoma datasets to replicate the pattern of results) and one monocular blind patients group. We note that blindness was defined as: “Light-perception negative for at least five years, due to perforation, enucleation, evisceration or ablatio retinae leading to blindness to the affected eye”. For the purpose of our study, monocular blindness acquired later in life was considered a reference model for the effect of reduced visual input to the brain. This is because due to unilateral sensory deafferentation half of the normal visual input is lost in these patients39.

Results

The goal of this study was twofold. First, we wanted to determine the extent of the WMM alterations in the glaucomatous brain outside of the visual WM40. Second, we wanted to establish whether the pattern of alterations across a large set of WM tracts was unique to glaucoma or was instead similar to those observed in a model of reduced visual input (i.e. monocular blindness). A total of 102 subjects participated in our study, including subjects with glaucoma (GL), monocular blindness (MBL) and healthy controls (HC). To validate our results in glaucoma, we used two independent datasets with major differences in the demographics and clinical characteristics of the subjects (one group from the Netherlands and another from Japan, see Methods: subjects). This resulted in a total of two glaucoma patients groups (GL1 and GL2) and three control groups (MBL, HC1 and HC2). The participant groups did not differ in mean age (GL1 vs HC1: p = 0.80, MBL vs HC1: p = 0.65, GL2 vs HC2: p = 0.90). The gender distribution did not differ in the NL subjects (GL1 vs HC1: p = 0.52, MBL vs HC1: p = 0.73), but differed in the JP subjects (GL2 and HC2 p = 0.04). The duration of the defect and the age at diagnosis/acquisition differed between the MBL and GL1 subjects (p = 0.02). All available demographics and clinical characteristics are summarized for the first cohort in Table 1 and the second cohort in Table 2.

Table 1 Characteristics of Cohort 1 (collection site at Netherlands; NL).
Table 2 Characteristics of Cohort 2 (collection site at Japan; JP).

Diffusion-weighted magnetic resonance imaging data was processed using an automated pipeline implemented as reproducible web-services on www.brainlife.io (see Table 3 for the full processing pipeline38,42,43). Using this pipeline, 37 major WM tracts were identified (Fig. 1a right-hand side; see also Table 4)37,44. After segmentation of all tracts, we computed WMM fractional anisotropy (FA) and mean diffusivity (MD) along the tract-length (Fig. 1a left-hand side; see also Methods: Tract profile generation)41,45.

Table 3 Description and web-links to the open source code and open cloud services used in the creation of this dataset.
Figure 1
figure 1

White matter anatomy segmentation. (a) Upper-left side. One representative subject white matter tract (the optic radiation) with overlayed microstructural estimate and representative optic radiation average tract profile for each participant group. Tract profiles were estimated using brainlife.io (brainlife.app.361;41). Bottom-right side. The 37 WM tracts identified using the brainlife.io white matter segmentation (brainlife.app.188). (b) White matter tract categories. Tracts were subdivided into 7 categories. These categories were defined by the degree of involvement in vision. Only WM tracts from the right hemisphere are visualized here for a representative subject.

Table 4 White matter tracts of interest and their classification.

The 37 WM tracts were also assigned to 7 categories varying from tracts highly related to vision and tract more clearly related to cognition and motor control (Fig. 1b; see also Table 4). Categories were defined as follows: (1) Early vision: the optic radiation, OR40, (2) Ventral visual stream: inferior lateral fasciculus (ILF) and inferior fronto-occipital fasciculus (IFOF)46, (3) Occipital vertical tracts (associative tracts within the occipital lobe): The vertical occipital fasciculus (VOF) and forceps major (CC FMaj)47, (4) Dorsal stream tracts: superior longitudinal fasciculus (SLF) 1, 2 and 3 and the cingulum cingulate (Cing)46, (5) Posterior vertical association tracts: temporal parietal connection (TPC), posterior arcuate fasciculus (pArc), the middle longitudinal fasciculus projecting to the angular gyrus (MdLF Ang) and projecting to the superior parietal lobule (MdLF SPL)37, (6) Cross callosum tracts (projecting through the body or genu of the corpus callosum): CC forceps minor (CC FMin), anterior (CC Ant), middle (CC Mid), and parietal (CC Par)48 and (7) Other tracts: fronto-thalamic radiation (FTR), uncinate fasciculus (UF), corticospinal tract (CST), frontal aslant tract (FAT) and the arcuate fasciculus (Arc).

The average fractional anisotropy (FA) and mean diffusivity (MD) for each tract in each participant was estimated and used to compute Hedges’ g effect size (independently for FA and MD)49. To compute g the following participant groups were compared: GL1 and HC1, MBL and HC1 and GL2 and HC2. The average (and standard error) effect size across all 37 WM tracts is reported in Supplementary Table S1 (see also Methods: Statistical analysis). This procedure returned three sets of 37 estimated g. The variable group consisted of three different levels (GL1 vs HC1, MBL vs HC1 and GL2 vs HC2) and the variable category consisted of seven levels (early vision, ventral, occipital, dorsal, vertical, callosal and other tracts (Fig. 1 and Table 4). For all the subsequent analyses Two-Way Factorial ANOVAs were conducted independently for FA and MD. The ANOVA compared the main effects of the subjects' group comparisons (see above) and tract category and the interactions between group comparisons and tract category using the Hedges’ g.

Which white matter tracts show the strongest difference between group comparisons?

To gain insight into which WM tracts were most different in our individual group comparisons, and before discussing the specific statistical tests, we provide a visualization of the Hedges’ g, effect sizes for all tracts ranked by magnitude (Fig. 2). We further describe the strongest effect sizes (g <>|0.7|) based on FA and MD (see also Table S2). For the first glaucoma cohort (Fig. 2a), we found large effect sizes were observed for the right VOF (g = − 0.90), bilateral ILF (L g = − 0.72; R g = − 0.72) and bilateral OR (L g = − 0.72; R g = − 0.70). When comparing MBL vs HC1, the left IFOF (g = 0.77) revealed a large effect size. For the second glaucoma cohort, we found large effect sizes for the left OR (g = − 0.91), bilateral ILF (L g = − 0.76; R g = − 0.89) and bilateral IFOF (L g = − 0.82 R g = − 0.81). A similar analysis was performed for MD to examine effect sizes for each tract (Fig. 2b). For our first glaucoma cohort, MD did not reveal large effect sizes. On the other hand, when comparing MBL vs HC1, we found a large effect for the left IFOF (g = − 0.83) and left CST (g = − 0.71). Our second glaucoma cohort revealed large effect sizes for MD in the parietal CC (g = 0.75) and right OR (g = 0.73). As can be seen in Fig. 2a, when examining FA g estimates, tracts involved in visual processing (e.g. OR, ILF, VOF, IFOF) were negatively affected in all group comparisons suggesting decreased FA compared to controls. For monocular blindness, tracts less involved in vision, however, revealed a positive effect size indicating increased FA in comparison to controls. This pattern is uniquely different from glaucoma group comparisons: both glaucoma groups reveal a negative effect on the WMM for the majority of WM tracts, indicating FA is decreased in glaucoma compared to controls. For MD g estimates (Fig. 2b), the pattern is less apparent, but can still be observed. For monocular blindness, the majority of tracts revealed a negative effect size (i.e., decreased MD in comparison to controls). In opposition, the majority of tracts in both glaucoma cohorts revealed a positive effect size (i.e., increased MD). This suggests that in monocular blindness a decreased FA in vision tracts is accompanied by an increased FA in non-vision tracts (this could be interpreted as a plasticity and compensatory process). In glaucoma instead no plasticity was observed in non-vision tracts.

Figure 2
figure 2

Ordered effect sizes in glaucoma and monocular blindness for each white matter tract. (a) Fractional anisotropy (FA). (b) Mean diffusivity (MD). Effect sizes were calculated using Hedges’ g by comparing the average FA or MD values from each WM tract profile. For each white matter tract, FA (or MD) in the glaucoma patients of the first (GL1) or second (GL2) cohort was compared to the corresponding tract in healthy controls (HC1 or HC2 respectively.) Similarly, FA (or MD) in the monocular blind patients (MBL) was compared to the corresponding healthy controls (HC1). Tracts are ordered by effect size with positive effect size presented first (top). WM tract acronyms are fully described in the Supplementary Information.

Does glaucoma affect the WMM similarly to monocular blindness?

Next, we wanted to test whether the effect of glaucoma on the WMM is primarily due to reduced visual input. To do so, we set out to determine whether glaucoma and monocular blindness affect the WMM similarly. We remind the reader that monocular blindness was used in our study as a model of decreased visual input. When focussing on FA we found a main effect (F(2,90) = 59.29, p < 0.0001) indicating a difference between GL1 vs HC1 (g = − 0.25 SE = 0.04), GL2 vs HC2 (g = − 0.28, SE = 0.04) and MBL vs HC1 (g = 0.28, SE = 0.04). Also for MD the main effect was significant (F(2,90) = 62.52, p < 0.0001), indicating a difference between GL1 vs HC1 (g = 0.19 SE = 0.04), GL2 vs HC2 (g = − 0.08, SE = 0.04) and MBL vs HC1 (g = − 0.37, SE = 0.04). To determine which group comparisons differed from each other, post hoc comparisons were made using the Tukey HSD test. These tests indicated that the effect size of the glaucoma patients differed from that of MBL group (FA p < 0.0001). We repeated the same analysis for our second cohort of glaucoma patients and found that the effect size of the glaucoma patients differed from the MBL group as well (FA p < 0.0001). Furthermore, we compared effect size for FA and MD between the two glaucoma cohorts and found no difference overall for (FA: p = 0.8638 and a small but significant difference for MD: p = 0.0043). This suggests that the glaucomatous- and visual-deprived brain respond differently to the loss of visual input.

Are white matter tract categories similarly affected?

Next, we were interested in establishing whether the WM tract categories were similarly or different in the group comparisons used to estimate g. For FA, a significant main effect of category was found F(6,90) = 16.14, p < 0.0001, indicating a difference between tract mean g for each WM tract category – early vision (M = − 0.65 SE = 0.08), ventral stream (M = − 0.41, SE = 0.06), occipital (M = − 0.20, SE = 0.07), dorsal (M = − 0.01, SE = 0.05), vertical (M = − 0.0, SE = 0.04), callosal (M = − 0.02, SE = 0.06) and other tracts (M = − 0.04, SE = 0.04). For MD, a significant main effect of category was found F(6,90) = 6.69, p < 0.0001, indicating a difference between tract mean g for each WM tract category – early vision (M = 0.24, SE = 0.08), ventral stream (M = 0.10, SE = 0.06), occipital (M = − 0.12, SE = 0.06), dorsal (M = − 0.03, SE = 0.05), vertical (M = − 0.0, SE = 0.04), callosal (M = − 0.17, SE = 0.06) and other tracts (M = − 0.15, SE = 0.04). These results show that both for FA and MD there was a significant difference in effect size across WM tract categories. This difference in effect size across tract categories suggests that the impact of each disease to tracts beyond the early visual WM is not uniform.

Which white matter tract categories are different between group comparisons?

Finally, we were interested in determining which WM tract categories differed between group comparisons. We first plotted the Hedges’ g organized by tract category and group comparison for FA (Fig. 3a) and MD (Fig. 3b). We found a significant interaction effect between group comparisons and the tract categories (F(12,90) = 3.23, p < 0.0007). To determine which tract categories differed between group comparisons based on the g estimates of FA, post hoc analysis was performed using Tukey HSD. Our post hoc analysis revealed that both glaucoma group comparisons differed from monocular blindness in tracts of the ventral and dorsal stream and callosal, vertical and other tracts (all p < 0.0024; Fig. 3a and Table 5). The early vision and occipital WM structures did not reveal differences between any of the group comparisons. None of the categories differed between the two glaucoma group comparisons. The g estimates of MD did not reveal a significant interaction effect between group comparisons and tract categories (Fig. 3b).

Figure 3
figure 3

Effect size on the white matter microstructure across tract categories. (a) Fractional anisotropy (FA). (b) Mean diffusivity (MD). The effect of glaucoma and monocular blindness on the WMM compared to healthy controls is visualized for each tract category (abscissa) and participant groups (colors). Positive effects sizes indicate that the microstructural estimate (of either FA or MD) was increased in the clinical group as compared to the healthy controls. Asterisks indicate significance for the comparisons between conditions (colors) with Bonferonni corrected p values (p < 0.05/21 (category*group) = p < 0.0024). The color of the asterix indicates conditions (colors) were different.

Table 5 Significant differences between groups per category based on the effect size of fractional anisotropy.

When focussing on FA, Hedges’ g was strongest and negative for the OR (early vision), with a magnitude similar across group comparisons (Fig. 3a; green, orange and blue). This indicated that FA for the OR was reduced for all patient groups. Hedges' g for the rest of the tract categories was strikingly different between group comparisons. Notably, whereas the effect size was primarily positive for the comparison between the monocular blind patients and controls (Fig. 3a orange), g was primarily negative for the comparisons between the two glaucoma groups and matched controls (Fig. 3a green and blue). A positive Hedges' g for the monocular blind subjects can be interpreted as potential indication for WM reorganization processes beyond the visual system as result of prolonged visual input reduction. The opposing pattern of results in the two independent glaucoma group comparisons can be interpreted as an indication of a possible lack of reorganization in the two glaucoma patient groups. The pattern of results for MD is similar –but opposite in the sign of g. Namely, whereas the two glaucoma groups have similar and positive g values in the majority of the tract categories, the monocular blinds subjects have negative g for the majority of the tract categories. As in the case of FA also for MD we found no difference between the Hedges' g across the three group comparisons of the early visual WM (i.e., the OR).

We performed an additional analysis to explore whether the results are robust to the patient’s selection criteria (see Supplementary Analysis 1). A subset of patients in the GL2 cohort (i.e those with the largest visual field defect, VFMD > − 2 dB in the better eye) were compared versus HC2. A similar pattern of results to that reported in Fig. 3 was observed indicating that our analyses were not biased by our selection criteria (see Supplementary Fig. S1 and Supplementary Table S4).

Is there an association between white matter tracts microstructure and patients’ clinical outcomes?

We performed additional analyses to evidence any association between WM tract properties and a series of clinical outcomes available (Visual Field Mean Deviation, VFMD, in the worse eye, age of diagnosis/acquisition, duration of the disease, and retinal nerve fiber layer measured by NFI or OCT) in the three cohorts (GL1, GL2 and MBL). For all these analyses, we first selected the tracts with the strongest effect size in the previous analysis (Fig. 2; g <>|0.7|) because these tracts demonstrated the most meaningful group differences. This approach returned 5 tracts for the GL1 cohort (bilateral OR, bilateral ILF and right VOF), 5 tracts for GL2 cohort (left OR, bilateral ILF and bilateral IFOF) and 1 tract for MBL cohort (left IFOF). We then used a Spearman's rank correlation coefficient (ρ) to associate individual tracts microstructure and the clinical outcome across subjects. In general, the associations were small (with the highest ρ’s ranging between |0.5| and |0.57|) and only a few of them passed statistical significance at p < 0.05. Below, we report specific details of these additional analyses that returned significant results.

For the first patient cohort (GL1), we examined the correlation between VFMD in the worse eye, retinal nerve fiber layer measured by NFI in the worse eye, age at diagnosis and duration of the disease with the FA values of the bilateral OR, bilateral ILF and right VOF. We found a positive correlation between the VFMD of the worse eye and the right ILF (ρ = 0.57, p = 0.02; See Supplementary Analysis 2 and Supplementary Fig. S2), a negative correlation between the age at diagnosis and the right ILF (ρ = − 0.57, p = 0.02), left ILF (ρ = − 0.55, p = 0.02), and left OR (ρ = − 0.52, p = 0.05). For the second patient cohort (GL2), we examined the correlation between VFMD in the worse eye and retinal nerve fiber layer measured by OCT in the worse eye with the FA values of the left OR, bilateral ILF and bilateral IFOF. No statistical significant correlations were found. For the third patient cohort (MBL), we examined the correlation between age at acquisition and duration of the disease with the left IFOF. No statistical significant correlations were found.

Discussion

This study set out to investigate the effect of glaucoma on the human WM as it is an ongoing debate whether glaucoma is purely an eye disease or whether it affects the whole brain. To address this, we characterized whether the effect of glaucoma on the brain WM is comparable for WM tracts supporting communication within the visual system and for WM tracts serving brain structures beyond the visual system. We used two independent datasets to validate our results.

We used two independent glaucoma datasets and the overall pattern of results was consistent between the two datasets. The effect of glaucoma on WM was strongest in the early visual pathways and it decreased as the tracts of interest were categorized as less involved in serving the visual system. We also used an additional control group of monocular blind individuals. This group was used as a model of reduced visual input to the brain. We compared the pattern of results in the three participant cohorts. We estimated which tracts and categories were most affected in the glaucoma subjects and monocular blind patients to test whether any change in the WMM in the glaucoma brain could be explained primarily by a reduction of visual input. We found that whereas the pattern of effects within the two glaucoma patients group was very similar, the effect was strikingly different from that measured in the monocular blind subjects. The results indicated that in monocular blind patients processes of reorganization of the large set of tracts WMM: FA was increased and MD decreased for WM tracts outside of the visual system. This can be cautiously interpreted as plasticity and reorganization following long-term visual deprivation to possibly aid behavior (even though we have no direct evidence that the pattern of alterations in FA and MD relate to behavior). In the two glaucoma groups we found no evidence for a similar reorganization process in the WM. We found that in both glaucoma groups, FA was reduced and MD increases across all tract categories as compared to controls. Overall, these results show that glaucoma shows a clear difference in the majority of the WM tracts with our model of reduced visual input (i.e. MBL). Below, we will discuss the current findings, the relevance for understanding the aetiology of glaucoma and its clinical relevance. As it becomes apparent from our results from two independent cohorts − collected at different sites with different parameters − glaucoma consistently shows that WMM is affected strongest in early vision WM structures which then gradually decreases towards cognition structures. The involvement of the early vision structures is consistent with previous studies5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20. Our novel findings show a graded alteration of the WMM. The most severe alterations were found in early vision structures, followed by less severe alterations in tracts categorized as motor-, or cognitive-related (Fig. 3). This pattern of graded increasing alteration was specific to the glaucoma patients and differed in the monocular blind patients.

The pattern of results reported here indirectly support the network degeneration hypothesis (NDH). NDH suggests that neurodegenerative diseases primarily spread through distinct brain networks, characteristic of the disease (Tahmasian et al. 2016; Greicius & Kimmel 2012). In light of NDH hypothesis, the RGCs, the visual pathways and the vision-related WM tracts (such as the OR, VOF and ILF) are the primary network for glaucoma and the spread of degeneration among these structures could be part of its clinical manifestation. In regards to glaucoma as a more general neurodegenerative disease, our findings show similarities with involved WM tracts in posterior cortical atrophy (PCA)—a subtype of dementia also referred to as "the visual variant of AD" (Levine et al. 1993; Bokde et al. 2001; Lee & Martin 2004). PCA primarily manifests itself as degeneration of the occipito-temporal (ventral stream) and occipito-parietal (dorsal stream) thereby causing both visuospatial and visuoperceptual problems (Crutch et al. 2012). Furthermore, previous literature has suggested that glaucoma and AD are connected through the same underlying pathologic mechanism (Ghiso et al. 2013, Inoue, Kawaji & Tanihara 2013, Janssen et al. 2013, Sivak 2013). However, the latter notion is still being questioned since conflicting epidemiologic reports have been reported. Some of these epidemiologic studies find an increased prevalence of glaucoma in AD (Bayer, Ferrari & Erb 2002, Chandra, Bharucha & Schoenberg 1986, Helmer et al. 2013, Tamura et al. 2006), while other studies do not (Bach-Holm et al. 2012, Kessing et al. 2007, Ou et al. 2012). Other shared characteristics between AD and glaucoma are the observation of abnormal tau protein in the vitreous jelly of the eyes of glaucoma patients as amyloid beta protein in the retina which are both considered pathologic hallmarks of AD50.

An interesting result to emerge from the data is that glaucomatous WM alterations are significantly different from non-glaucomatous alterations of the WM. This suggests that dissimilarities between glaucoma and decreased visual input alterations likely reflect different underlying mechanisms. In both glaucoma group comparisons, on average, a negative effect size for FA and a positive effect size for MD was observed indicating increased and decreased values in glaucoma patients compared to healthy controls respectively, while the opposite effect is observed for monocular blindness. Glaucoma groups did not differ from each other, despite the fact that glaucoma groups had different pathology (NTG vs monocular POAG). This strengthens the results of comparable patterns of alterations. As little is known about mechanisms that play a role, the discrepancies between glaucoma and monocular blindness groups may favor the hypothesis that WM alterations are not caused by decreased visual input alone, likely reflecting different underlying mechanisms. Although the interpretation of diffusion measures is not straightforward51,52, FA is assumed to be highly sensitive for microstructural changes in general, other scalars are more specific for the type of microstructural changes. MD is thought to be sensitive to cellularity, necrosis and edema53,54,55,56. We speculate that on one hand, decreased FA and increased MD in glaucoma may reflect axonal degeneration and demyelination, respectively. On the other hand57, in monocular blindness the decreased MD and increased FA may reflect WM maturation and high myelination, respectively54,55,56. In support of this latter notion, decreases in MD and increases in FA have been observed after training, supporting the idea that this reflects neuroplasticity or reorganization57. Our findings corroborate a meta-analysis of 10 dMRI studies examining diffusion values (i.e., FA and MD) in patients with glaucoma, revealing FA decreases and MD increases in the early vision structures. Therefore, we hypothesize that glaucomatous changes are a consequence of a disease progress that involves neurodegeneration, whereas in monocular blindness WM alterations are a consequence of neuroplasticity and reorganization following loss of input.

Several studies have previously reported that glaucoma worsens with age10. An alternative proposal has been to consider glaucoma as a neurodegenerative disease that accelerates brain aging and also involves damage to the optic nerve18,58,59. Previous studies have suggested that visual WM structures such as the ILF or IFOF are particularly prone to age-related decline, showing the highest rate of decline as compared to other WM tracts60. Because our subjects were age-matched, we were not in the position to directly analyze the clinical outcomes, or WMM, as a function of the patients’ age. Instead, we looked at the age at which glaucoma was first diagnosed (Supplementary Fig. S2). In line with the proposal, we found a negative association between the age at which glaucoma was diagnosed and WMM in two visual pathways—the OR and ILF. More specifically, the older the age at which glaucoma was diagnosed, the lower the FA in cohort GL1 (Supplementary Fig. S2). The ILF FA was also positively associated with VFMD, indicating the greater the visual field defect the worse the WMM. These results are especially interesting when contrasted with those in the MBL cohort. Indeed, in contrast to the glaucoma patients, results with the MBL patients did not show an association between age of acquisition or duration of the visual defect and FA. Finally, several issues stand between a full interpretation of the results above and suggest caution. In most if not all glaucoma cases, visual field loss remains unnoticed until it becomes severe, so age of diagnosis is only a marginal proxy to look at a degenerative etiology of the disease. Given the difficulty to establish disease onset, we believe that future studies will be necessary to address how age of disease onset or progression rate affects WMM, VFMD and their association.

Our research advances the scientific understanding of glaucoma by using open science methods and data sharing. In addition to improving understanding WM alterations and their potential underlying mechanisms in glaucoma, we have used the open cloud computing platform brainlife.io for secure neuroscience data analysis. The brainlife.io platform allows publishing processed data and associated analyses steps (Apps) integrated into a single record referenced by a digital-object-identifier (DOI)38. This is a major advancement in supporting reproducibility efforts because it removes the complexities involved with modifying the computational environment for new software, increasing the ability of researchers to easily replicate or add to prior work. All of our data and methods are fully available to the wider scientific community on https://doi.org/10.25663/brainlife.pub.18.

Our research has several clinical implications. First, affected WM integrity shows that damage extends beyond the visual pathways in glaucoma. In addition, we show that monocular glaucoma and monocular blindness affects the WM differently. Determining the full extent of the alternation to the brain WM due to glaucoma could guide the future developments, for example suggesting the span of brain tissues to best target with neuroprotective treatments61,62,63,64. For glaucoma, such neuroprotective medication is currently under development (Doozandeh & Yazdani 2016). Neuroprotection may either prevent healthy RGCs from dying or “slow down” the process of already sick RGCs. Combining neuroprotection with IOP reduction—which is clinically proven to prevent RGC loss in glaucoma—may therefore optimize treatment.

As for all dMRI studies, the interpretation of differences between groups in diffusion values is not straightforward51,52. Differences may reflect WM abnormalities, but can also be affected by the anatomical features of a tract, such as its curvature, by partial volume effects with neighbouring fiber tracts due to crossing, and by merging or kissing fibers. Fibers originating from the occipital lobe can merge, kiss and cross with fibers from neighbouring tracts65. For example, WM alterations in the OR may consequently have caused (part of) the diffusion changes in neighbouring tracts such as the forcep major as well. However, changes in merging or crossing of fibers could also be part of clinical manifestation itself. Further, for our first cohort (NL), the T1-weighted image had a relatively low contrast ratio between grey- and WM compared to the second cohort (JP). The T1-weighted image is used for segmentation of ROIs and to create a GMWMI mask to constrain tractography. We resolved this by using different WM intensity level settings for each dataset. Although we visually checked all segmentations and masks, a lower contrast ratio between grey- and WM in the first cohort could result in misclassifications in WM segmentation. Additionally, susceptibility-induced distortions can only be corrected using a second diffusion scan acquired with opposing phase encoding directions. As common in clinical studies, the JP cohort did not include a second scan. As a result distortion correction for susceptibility artifacts was not applied on this dataset. Furthermore, our results may be influenced by differences in clinical outcomes for each patient. For example our subjects showed differences in disease duration, IOP treatment, glaucoma form (high- versus normal-tension glaucoma) and visual field defects (binocular or monocular glaucoma). Future studies would be necessary to either measure groups with more homogeneous clinical measures or with larger groups of patients so to allow measuring explicitly the relation between each clinical measure the the individual variability in the groups. We make all our data readily available online at DOI. We believe this is important for evaluating and reproducing our results but also for allowing future studies to build upon our current dataset. This will ultimately improve our understanding of the heterogeneous character of glaucoma.

In conclusion, glaucoma is associated with widespread WM reduction in fractional anisotropy and increase in mean diffusivity in both early vision as well as non-vision related WM tracts. Although the early vision and ventral tracts are affected in monocular blindness as well, the WM changes in glaucoma show a different and opposing pattern. Finding such widespread abnormalities in glaucoma suggests the presence of degeneration beyond what can be explained on the basis of decreased visual input. Together with degeneration of the RGCs, the degeneration of visual pathways and ventral WM tracts may be part of the manifestation of glaucoma.

Materials and methods

Data sources and App processing records

All data are de-identified and published online using the brainlife.io platform. Data and Apps are interlinked and each dataset is preserved with associated provenance information to allow other researchers to download the minimally processed or processed data and to reuse the Apps on the cloud platform on their own data. Data and Apps associated with the present project are published and archived at https://doi.org/10.25663/brainlife.pub.18.

Study subjects

Data from two sites were included in the current study. Data were collected among the ophthalmologic patient populations of the University Medical Center Groningen (Groningen, the Netherlands (NL)) and Jikei University hospital (Tokyo, Japan (JP)). This resulted in a total of two patient groups and three control groups. The first cohort included 17 NL subjects with primary open angle glaucoma (POAG) with a monocular visual field defect (GL1), 20 healthy age-matched NL subjects (HC1) and 14 NL age-matched subjects with non-glaucomatous monocular blindness (MBL). A second glaucoma cohort was included to validate our results. For the second cohort, we included 30 JP subjects with glaucoma with a high prevalence of normal tension glaucoma (NTG) (GL2) and 21 healthy age-matched JP subjects (HC2). The first cohort has been previously studied using morphometry of structural MRI data (T1-weighted images)66, whereas the second cohort has been previously studied using Tract-Based Spatial Statistics (TBSS)18.

The inclusion criteria can be found in the Supplementary Methods. Characteristics of the NL subjects can be found in Table 1 and of the JP subjects in Table 2. In addition to Tables 1 and 2 of the average clinical characteristics of our cohorts, all individual subject data and their disease characteristics are openly available online and accessible on https://doi.org/10.25663/brainlife.pub.18. This study conformed to the tenets of the Declaration of Helsinki. The study on the NL subjects was approved by the medical ethical committee (METC) of the University Medical Center Groningen while that on the JP subjects was approved by the METC of the Jikei University School of Medicine. All methods were performed in accordance with the relevant guidelines and regulations. All subjects gave their informed written consent prior to participation.

Clinical data acquisition

Visual acuity was measured for all NL subjects with a Snellen chart with optimal correction for the viewing distance.

Visual fields of the GL1, GL2 and HC2 subjects were assessed using a Humphrey Field Analyser (HFA; Carl Zeiss Meditec, Dublin, CA, USA) 30-2 Swedish Interactive Threshold Algorithm SITA) fast. An abnormal visual field was defined as a glaucoma hemifield test ‘outside normal limits’ for at least two consecutive fields; defects had to be compatible with glaucoma and without any other explanation; for a normal visual field the glaucoma hemifield test had to be within normal limits. Visual field loss is reported in terms of mean deviation (VFMD). The Frequency Doubling Technology perimeter (FDT; C20-1 screening mode) was used in the MBL and HC1 subjects. An abnormal test result was defined as at least 1 reproducibly abnormal test location at p < 0.01. Both eyes of the HC1 subjects and the contralateral eye of the MBL subjects had to be normal.

Retinal nerve fiber layer (RNFL) was assessed for the GL1 patients using laser polarimetry (GDx; Carl Zeiss Meditec, Jena, Germany); outcome measure was the Nerve Fiber Indicator (NFI) (a summary value that indicates the likelihood of glaucomatous RNFL loss). RNFL thickness of the GL2 patients was measured by means of Optical Coherence Tomography (OCT; Stratus OCT 3000; Carl Zeiss Meditec, Dublin, CA, USA).

Neuroimaging data acquisition

Cohort 1 (collection site: Netherlands) NL subjects were imaged using a 3.0 T Philips Intera scanner (Philips, Eindhoven, The Netherlands) . An 8 channel head coil was used. Diffusion-weighted magnetic resonance imaging (dMRI) data were collected with two phase-encoding schemes, i.e. anterior–posterior (AP) and posterior-anterior (PA). The following parameters were used for the dMRI pulse sequence using DwiSE technique: 60 diffusion-weighting directions including an additional b0 image, resolution = 128 × 128, APP fat shift = 11.722 pixels, degree of angulation = 0.57 1.405 − 15.9, repetition time = 5485 ms, field of view = 240 × 102 × 240, echo time = 79 ms, diffusion b value = 800 s/mm2, EPI factor = 45 and voxel dimension = 1.875 × 1.875 × 2.000 mm. The second DW scan has similar parameters as the first, except for: APA fat shift = 2.035 pixels and repetition time = 5516 ms. One T1-weighted (T1w) anatomical image was acquired for each participant using T1W/3D/TFE-2 sequence and following parameters: flip angle = 8 degrees, TR = 8.70 ms, TE = 2.98 ms, matrix size = 256 × 256, field of view = 230 × 160 × 180, number of slices = 160 slices and voxel size = 1 × 1x1mm3.

Cohort 2 (collection site: Japan) JP subjects were imaged using a 3.0 T Siemens MAGNETOM Trio A Tim System scanner (Siemens, Erlangen, Germany) with a 32 matrix head coil. Diffusion-weighted magnetic resonance imaging (dMRI) data were collected. The following parameters were used for the dMRI pulse sequence using single shot echo planar sequence. The scanning parameters were: 64 diffusion-weighting directions including an additional b0 image, TR = 8800 ms, TE = 87 ms, field of view = 210 ~ 230, matrix = 140, slice distance factor = 5%, voxel size = 1.5 × 1.5 × 1.575 mm, b value = 1000 s/mm2 and acceleration factor = 2 (GRAPPA). One T1w anatomical image was acquired for each participant using the following T1W 3D MPRAGE sequence using the following parameters: flip angle = 9 degrees, TR = 2300 ms, TE = 2.98 ms, matrix size = 256 × 256, field of view = 256, number of slices = 176 slices and voxel size = 1 × 1x1mm3. The brain images were corrected for intensity inhomogeneity using the built-in console.

Data processing using open cloud computing services on brainlife.io

All data was processed using the reproducible, open cloud computing services available at brainlife.io38,67,68. The workflow used consisted of a series of 20 brainlife.io Apps (Table 3) and is described in detail below. All the Apps used for this study can be accessed at https://doi.org/10.25663/brainlife.pub.18. Table 3 reports the DOIs to reproducible web services on brainlife.io and the URLs to the source code on GitHub.com.

Anatomical data (T1w) processing

The T1w anatomical images were preprocessed using the FMRIB Software Library (FSL) function fsl_anat (brainlife.app.273). Using this app, the T1w anatomical images were reoriented and linearly-aligned to the MNI152 1 mm atlas using FLIRT, and non-linearly aligned to the same atlas using FNIRT. Subsequently, the anatomical images were segmented into multiple tissue types (e.g. cortical grey matter, subcortical grey matter, WM, cerebrospinal fluid and pathological tissue) using MRtrix3’s 5ttgen function69 (brainlife.app.239). A mask representing all of the voxels of grey matter-white matter interface (GMWMI) was created and used to seed tractography (see Methods: White matter microstructure modelling). Additionally, the aligned T1w anatomical image was used for segmentation and surface generation using Freesurfer’s recon-all function70 (bl.app.0). For the NL dataset, we used a WM intensity level of 105 using the flag -seg-wlo. Next, the bilateral lateral geniculate nuclei were segmented using Freesurfer’s developmental segment_thalamic_nuclei function71 (brainlife.app.222). To be able to reconstruct regions of interest (ROIs) for tractography of the OR, population receptive field (pRF) maps were fit to the cortical surfaces generated from Freesurfer using a retinotopy atlas (brainlife.app.187). This app performs a retinotopic mapping and segmentation of visual areas V1, V2, and V3, as well as higher-order visual areas, from from the surface topology of the T1w anatomical images using an open source library (github.com/noahbenson/neuropythy)72.

Diffusion data (dMRI) processing

The DWI images were preprocessed using the procedure dwiprepoc provided by MRtrix3 (bl.app.68). The app is largely based on the DESIGNER pipeline (Diffusion parameter EStImation with Gibbs and NoisE Removal), developed to identify and minimize common sources of methodological variability. The steps consisted of performing PCA denoising, correction on Gibbs ringing artifact, motion and eddy current correction with FSL, bias correction with ANTs N4 and removal of Rician background noise. If AP and PA diffusion scans are fed into this App, FSL’s TOPUP73,74 is additionally run as well. This was the case for the NL data and as such, AP and PA diffusion scans were combined into a single corrected one by estimating the susceptibility-induced off-resonance field using TOPUP. This step was omitted for the JP data automatically by the brainlife.io App because only one diffusion scan was provided. As a result susceptibility-induced distortions were not corrected in the JP data. DWI slices whose average intensity is at least four standard deviations lower than the expected intensity were considered as outliers and replaced with predictions made by the Gaussian Process. The DWI data and gradients were aligned to the ACPC aligned T1-weighted anatomical scan using FSL’s BBR. The final isotropic voxel dimension in mm was set to 1.875 and 1.5 for the NL and JP data, respectively.

Region-of-interest (ROI) generation

A multiple region of interest (ROIs) approach was used to be able to reconstruct the OR using the roiGenerator app (brainlife.app.223). This app uses AFNI's 3dROIMaker functionality75 to generate ROIs in dMRI space by taking the Freesurfer’s parcellation volume, thalamic nuclei segmentation and pRF visual area segmentation, to generate exclusion ROIs, the lateral geniculate nuclei (LGN) and V1, respectively. Exclusions ROIs used from the Freesurfer’s aseg segmentation were as follows: the bilateral cerebral WM, cerebral cortex, lateral ventricle, cerebellum WM and cerebellum cortex.

White matter microstructure modelling (tractography)

Probabilistic whole-brain fiber tracking was performed using the MRtrix3 app76 (bl.app.101) which uses MRtrix3’s Anatomically Constrained Tractography (ACT) algorithm69. Constrained spherical deconvolution (CSD) was used to estimate local fiber orientation distributions (FOD) (Tournier et al. 2007). CSD estimates were generated with a maximum spherical harmonic order of Lmax = 8. The GMWMI mask was used as a seed mask to constrain fiber tracking. The following parameters were used: a maximum curvature angle of 45° between successive steps, a step size of 0.2 mm, a maximum track length of 250 mm and a minimum length of 29 mm. A total of 1,500,000 streamlines per parameter combination was produced.

White matter microstructure modelling (tracking the human optic radiation (THORA))

The OR was reconstructed using the THORA App (brainlife.app.347). This app will generate ORs using MRtrix3’s iFOD2 algorithm with the capability of using the ACT framework and performs backtracking between the bilateral LGN and bilateral V1 as an end ROI. Previously generated exclusion ROIs were used to exclude fibers that are not part of the OR. The following parameters were used: a maximum spherical harmonic order of Lmax = 8, number of repetitions, minimum and maximum length of streamlines of 10 and 180 mm, respectively, a step size of 0.5 mm, a streamline count of 1000, a minimum FOD amplitude of 0.05 as a cutoff and a maximum curvature angle of 35°.

White matter microstructure modelling (segmentation)

To examine WM differences in early vision, ventral stream, dorsal stream, occipital, vertical, callosal and other tracts, we identified and segmented a total of 37 WM tracts (see Table 4) using the app WM anatomy segmentation37 (brainlife.app.188). This app Automatically segments a tractogram into the major WM tracts using the Freesurfer 2009 parcellation and the ACT tractography. For an overview of the segmented tracts see Fig. 1.

Remove tract outliers

To remove spurious fibers, tracts were cleaned by removing fibers using a brainlife.io wrapper created for the mbaComputeFibersOutliers algorithm. It takes the tract classification and prunes classified fibers that are unlike other fibers within the same tract (brainlife.app.195). Tracts that deviate more than 3 standard deviations (SDs) from the core of the tract (i.e. maximum distance) as well as fibers that are more than 4 SDs above the mean length of the tract (i.e. maximum length) were removed.

Tract profile generation

Tract profile generation is performed after the reconstruction of tracts via tractography and subsequent segmentation. For each of the segmented tract, plots of diffusion metrics (i.e. FA, MD, RD, AD) were created, known as Tract Profiles41; bl.app.43). The app utilizes the Compute_FA_AlongFGcommand provided by Vistasoft. For each tract the core representation is estimated by applying a gaussian weight to each streamline based on distance away from the “core”. Each tract is resampled into 100 equidistant nodes along their trajectories and the weighted average metric at each node for each participant can be extracted.

Quality assurance

All raw data was manually inspected by two expert observers. Processed images were evaluated for quality of processing and data quality using the streamlined method available on brainlife.io. To check the quality of our analyses, we used a series of Quality Assurance (QA) Apps provided by brainlife.io that allows evaluating the quality of the initial data and the processing steps. Using these Apps, we generated images of each subject of the T1-weighted scan, tissue type masks, DWI scans and DWI scan overlaid on T1-weighted scan. We further plotted the response function, computed SNR and calculated statistics associated with average streamline characteristics of each segmented tract (i.e. count, volume occupied, average length and length distribution). Table 3 provides an overview of the Apps used. The output of these QA Apps was visually inspected for each participant by two expert observers. The output is openly available on https://doi.org/10.25663/brainlife.pub.18 (see tag “QA_image”) and readers can appreciate it by visualizing the data from the platform to judge the quality of the processing achieved with the pipeline. Additional information about how to visualize data using brainlife.io can be found here: https://brainlife.io/docs/user/tutorial/#visualize-dataset.

Statistical analyses

Extraction of white matter properties of tracts For each WM tract, we computed diffusion along the tract by resampling each tract into 100 equidistant nodes. At each node, FA MD were calculated. Individual tracts were excluded from the analysis if they had less than 50 streamlines (see Supplementary Table S3 for overview of excluded tracts per group). Additionally, the first and last 10 nodes of each tract were excluded from the analyses since the ends of the tract may contain ambiguous voxels of the grey-white matter junction. The average FA and MD of each tract was calculated.

Effect size calculation To examine the effect of glaucoma and monocular blindness on the WMM compared to healthy controls, we calculated the effect size for both FA and MD separately between groups using the average tract values. We compared subjects of GL1 vs HC1 and MBL versus HC1. To validate results, we additionally compared GL2 vs HC2 in our second cohort. Due to heterogeneity of our groups (neuroimaging parameters and disease characteristics, effect sizes were calculated using within-cohort group comparisons only. This is an approach similar to that used by ENIGMA, the largest multi-site neuroimaging study of clinical populations. The approach has been benchmarked by the consortium and it has been shown to allow comparison of very heterogeneous measurements and clinical populations77.

The effect sizes were computed using Hedges’ g for each tract using The Measures of Effect Size (MES) Toolbox49 (github.com/hhentschke/measures-of-effect-size-toolbox). Hedges’ g was calculated using the following equation (1): \(Hedges' g = \frac{{\bar{y}}_{1}- {\bar{y}}_{2}}{{s}_{p}}\), with \({\bar{y}}_{1}\),\({\bar{y}}_{2}\) and \({s}_{p}\) denoting the mean of sample 1, the mean of sample 2, and the pooled standard deviation, respectively. The equation for the pooled standard deviation is \({s}_{p}=\surd \frac{({{n}_{1}}-1){s}_{1}^{2}+({n}_{2}-1){s}_{2}^{2}}{ {n}_{1}+{ n}_{2}-2}\), with \({s}_{1}\) and \({n}_{1}\) denoting the standard deviation and number of observations for sample 1, respectively, and \({s}_{2}\) and \({s}_{2}\) denoting the standard deviation and number of observations for sample 2, respectively. Additional subsets of GL2 were also compared versus HC2 and are included in the Supplementary Analysis 1.

Two-Way Factorial ANOVAs were conducted independently for FA and MD to compare the main effects of participant group and tract category and the interactions between group and category on the g estimates. The variable group consisted of three different levels (GL1 vs HC1, MBL vs HC1 and GL2 vs HC2) and the variable category consisted of seven levels (early vision, ventral, occipital, dorsal, vertical, callosal and other tracts.

What are the similarities between glaucomatous WM alterations and those caused by decreased visual input? To examine whether the effect of glaucomatous alterations is similar to decreased input (as measured by monocular blindness), we investigated the main effect of group for Hedges’ g FA and MD, separately, using a Two-Way Factorial ANOVA. The main effect of group consisted of three levels (GL1 vs HC1, MBL vs HC1, GL2 vs HC).

Are alterations in the visual WM similar to those in the rest of the WM?

To examine whether alterations in visual WM are similar to the rest of the WM, all WM tracts were categorized based on their involvement in vision. Categorizations were made as follows: (1) early vision tracts are tracts part of the visual pathways40, (2) ventral stream tracts are associated with ventral stream processing46, (3) dorsal stream tracts are associated with dorsal stream processing46, (4) occipital tracts are defined as tracts originating from or projecting to the occipital lobe47, (5) vertical tracts are defined as vertically-oriented tracts37, (6) callosal tracts are tracts projecting through the body or genu of the corpus callosum48) and (7) other tracts are defined as WM tracts that are not categorized in previous categories. An overview of these categorizations of WM tracts can be found in Table 4 and Fig. 1. Next, the Hedges’ g value for FA and MD of each tract was averaged per tract category to create an average effect size measure of WM property changes in glaucoma groups and non-glaucomatous blindness compared to healthy controls. To examine whether the effect of glaucomatous alterations in visual WM is similar to other WM tracts), we investigated the main effect of WM tract categories for Hedges’ g FA and MD, separately, using a Two-Way Factorial ANOVA.

Interaction effect group and category To examine whether WM tract categories are differently affected between groups, we examined the interaction effect between group and categories using Two-Way Factorial ANOVA. If significant at p < 0.05, then post hocs using Tukey HSD were performed to examine which categories differed between the groups. If significant at p < 0.05, then post hocs were performed to examine which categories differed within the group. Post hoc results were Bonferroni corrected for multiple comparisons and significance was set at p < 0.05/21 (category*group) = p < 0.0024.

Investigating the relation between clinical outcomes and the white matter tracts microstructure. We investigated the association between WM tract microstructure properties and a series of clinical outcomes available for the difference datasets. More specifically, VFMD in the worse eye, age of diagnosis/acquisition, disease duration, and RNFL measured by NFI or OCT were used. For each of these analyses, we first selected the tracts with the strongest effect size in the previous analysis (Fig. 2; g <>|0.7|) because these tracts demonstrated the most meaningful group differences. Five tracts passed the threshold of g <>|0.7|, more specifically for the GL1 cohort: bilateral OR, bilateral ILF and right VOF. Five tracts passed the threshold of g <>|0.7| for the GL2 cohort: left OR, bilateral ILF and bilateral IFOF. Only one tract passed the threshold of g <>|0.7| for the MBL cohort: the left IFOF. Spearman's rank correlation coefficient was used to associate individual tracts microstructure and the clinical outcome across subjects. More specifically for the first cohort we correlated VFMD in the worse eye, age of diagnosis/acquisition, duration of the disease, and RNFL measured by NFI with bilateral OR, bilateral ILF and right VOF. For the second cohort we correlated VFMD in the worse eye and RNFL measured by OCT in the worse eye with left OR, bilateral ILF and bilateral IFOF and for the last cohort we correlated age of acquisition and disease duration and with the left IFOF.

Demographics Group differences in age, duration of defect, age at diagnosis/acquisition were tested using two‐sample t-tests. Gender differences were tested using Fisher's exact tests.