Spatial analysis of thickness changes in ten retinal layers of Alzheimer’s disease patients based on optical coherence tomography

The retina is an attractive source of biomarkers since it shares many features with the brain. Thickness differences in 10 retinal layers between 19 patients with mild Alzheimer’s disease (AD) and a control group of 24 volunteers were investigated. Retinal layers were automatically segmented and their thickness at each scanned point was measured, corrected for tilt and spatially normalized. When the mean thickness of entire layers was compared between patients and controls, only the outer segment layer of patients showed statistically significant thinning. However, when the layers were compared point-by point, patients showed statistically significant thinning in irregular regions of total retina and nerve fiber, ganglion cell, inner plexiform, inner nuclear and outer segment layers. Our method, based on random field theory, provides a precise delimitation of regions where total retina and each of its layers show a statistically significant thinning in AD patients. All layers, except inner nuclear and outer segments, showed thickened regions. New analytic methods have shown that thinned regions are interspersed with thickened ones in all layers, except inner nuclear and outer segments. Across different layers we found a statistically significant trend of the thinned regions to overlap and of the thickened ones to avoid overlapping.


Results
Demographic and clinical data. Data for patients with mild AD and age-matched controls are shown in Table 1. The two groups did not show statistically significant differences in age, gender composition, educational level or refractive error. In fact, the distributions of refractive error of both groups were similar and the difference was not statistically significant (Kolmogorov-Smirnov Z = 0.7632; p = 0.603). Educational level was categorized in terms of years of education as 1, <9 years; 2, 9-17 years; or 3, >17 years. All patients scored >17 on the MMSE (Mini-Mental State Examination), and the two groups showed a statistically significant difference in mean score.
Spatial distribution of thickness in each layer. Maps were generated showing mean thickness for a central square measuring 4.8 × 4.8 mm 2 in each retinal layer (Fig. 2). Thickness appears color-coded in Fig. 2a (patients) and 2b (controls). The search for regions with statistically significant thinning was restricted to a circle   Top row shows data of total retina, and subsequent rows show data for its 10 layers. Column (a) corresponds to the AD group and shows at each scanned point its color-coded thickness mean, as entered in the RFT-based statistical analysis, after preprocessing (i. e. after thickness correction, spatial normalization and Gaussian smoothing for signal to noise ratio enhancement; see Methods section for details). Column (b) shows the corresponding information for the control group. Column (c) displays the statistical parametric map (SPM) for thickness difference (AD-Control), thus giving for each scanned point the z-score corresponding to the difference in mean thickness between the two groups. Column (d) shows in red color the regions whose thinning in the AD group reached statistical significance; thus each red region is formed by the set of scanned points where the z-score of the difference (AD-Control) exceeded the critical value Z c ; a Z c value was obtained for each retinal layer from RFT analysis with the probability of family-wise type I error = 0.05, and is indicated at the bottom of the corresponding figure in the rightmost column.
www.nature.com/scientificreports www.nature.com/scientificreports/ of diameter 3.5 mm centered on the fovea and drawn in black in Fig. 2. A random field theory (RFT)-based search for regions thinner in patients than in controls led to the statistical parametric map 9 (SPM) in Fig. 2c, where green color means no difference; blue, thinner in patients; and yellow-red, thicker in patients. Figure 2d shows in red color on a reference anatomical image the regions of total retina and its layers where thinning in AD group reached statistical significance, based on a one-sided test with probability of family-wise type I error set at 0.05. The normality and homoscedasticity checks required by RFT techniques revealed statistically significant deviations from these assumptions in wide regions of the OPL, IS/OS, OSL and OPR layers, but these deviations did not alter our results, since similar results were obtained using a non-parametric random permutation test. The results showed a statistically significant thinning of the total retina of patients throughout its central region, including the fovea and parafovea, as well as thinning spread out in the upper nasal direction. Regions with statistically significant thinning were also found in the NFL, GCL, IPL, INL and OSL, at the locations shown in the corresponding rows of Fig. 2d. Thinning in the NFL was extrafoveal and in the upper side; as one moved toward outer layers, the thinning remained extrafoveal and appeared to rotate counterclockwise. Systematic but not statistically significant thinning was observed in regions of all other layers except RPE.
Next we analyzed the spatial distribution of AD-related changes in thickness in greater detail; differences in thickness between groups are shown in Fig. 3a. In patients, all layers were thinner to different extents and with different spatial distributions, with the exception of the RPE, which showed mainly thickening. For each retinal layer Fig. 2d offers a precise delimitation of the regions with a statistically significant thinning in the patient group: (1) NFL showed statistically significant thinning in the upper part of the parafovea; (2) GCL was thinner in a c-shaped area around the temporal side of the foveola, and this thinning reached statistical significance at the temporal side close to the foveola center; (3) IPL showed statistically significant thinning in the lower temporal quadrant of the parafovea and fovea, even within the foveola but not within the foveal pit; (4) INL was thinner all around the foveola, and this thinning was statistically significant in both lateral regions of the parafovea and in the lower region of both fovea and parafovea; (5) OPL showed only slight, not statistically significant thinning in the lower nasal quadrant of fovea and parafovea; (6) ONL showed thinning, which was not statistically significant, over the upper temporal half, including the entire foveal avascular zone; (7) IS/OS showed slight thinning at the fovea, without reaching statistical significance; (8) OSL showed a statistically significant decrease of perifoveal thickness in a comma form in the superior nasal area; (9) OPR showed only slight, not statistically significant thinning in the temporal fovea and parafovea; and (10) RPE also showed slight, not statistically significant thinning in the nasal area.
Some thickened regions were also clearly observed in the ONL and RPE, though the thickening did not achieve statistical significance. This thickening may be related to inflammatory processes and may have concealed thinning in previous studies where thickness was averaged among regions of different sizes and geometrical shapes. In NFL, a small blob of thickening covered the half-nasal part of the foveal avascular zone, outside the foveal pit. GCL showed some areas of thickening in the upper areas analyzed, as well as a small blob in the inferior macular area. ONL was thicker in the lower nasal para-and perifovea regions. RPE showed a wide thickened region dominating the temporal half of the scanned area. Thickening found in the NFL, GCL and RPE did not reach statistical significance.

Area of thinned and thickened regions in retinal layers. Differences in mean thickness between
patients and controls were studied at each scanned point within a 4.8 × 4.8-mm 2 central square extracted from the 6 × 6-mm 2 scanned area. Differences in thickness between the two groups are shown in Fig. 3a. The corresponding histogram in Fig. 3b shows the bias of the difference in the negative direction. This confirms the systematic dominance of thinned areas in all layers except RPE. Figure 3c shows regions that were classified as thinned (black) or thickened (white) in patients, based on thickness differences of any size between AD and control group. Figure 3d uses the same color scheme to show regions classified as thinned or thickened in patients when a thickness difference greater than 1 SD was required; regions with smaller differences are shown in gray. After identifying thinned and thickened regions for each layer of AD patients, we calculated two indices based on their corresponding areas. The first index, t 0 , was the percentage of thinned area when all differences were considered (regardless of their magnitude); the second t 1 , the percentage of area whose thinning exceeded 1 SD. Analogously, the percentages of thickened area under the same two conditions, T 0 and T 1 , were also calculated. Numeric values obtained for each index at each retinal layer are displayed under the corresponding images in Fig. 3c,d.
Graphical analysis of these indices (Fig. 4) shows that all layers were thinner in patients over an appreciable percent of their surface. In addition, all layers except INL and OSL contained areas that were thicker in patients. In neural layers, thinned area systematically predominated over thickened area: for all 9 neural layers, mean ± SD was 52.24 ± 30.98 for the difference t 0 -T 0 (Student's t = 5.06, df = 8, p < 0.001) and 38.91 ± 24.14 for the difference t 1 -T 1 (Student's t = 4.84, df = 8, p < 0.001). In RPE the thickened area was twice as large as the thinned area.

Volume of thinned and thickened regions in retinal layers.
To investigate AD-related retinal volume changes in all retinal layers, we defined, for each layer, v o as the ratio of the volume lost in thinned regions to the layer's mean volume in the AD group, and V 0 as the ratio of volume gained in thickened regions to the layer's mean volume in the AD group; these two ratios, v o and V 0 , took into account any thinning and thickening, regardless of its magnitude. We also defined the corresponding ratios v 1 and V 1 for the case in which only a thickening or thinning greater than 1 SD was considered. The results obtained for all layers (Fig. 5) paralleled those found for area changes. In neural layers, the ratio of volume lost in thinned regions was systematically greater than the ratio of volume gained in thickened regions: for all 9 neural layers, mean ± SD was 0.0528 ± 0.0546 for the difference v 0 -V 0 (Wilcoxon T = 0, p < 0.005) and 0.0476 ± 0.0531 for the difference v 1 -V 1 (Wilcoxon T = 0, p < 0.005). In RPE, the volume gained in the thickened regions exceeded that lost in the thinned ones. www.nature.com/scientificreports www.nature.com/scientificreports/

Spatial correspondence between thinned and thickened regions of different retinal layers.
After observing the spatial location of thinned and thickened regions in each of the retinal layers (Fig. 3c,d), we examined whether these regions overlap across the layers. A Jaccard similarity index was computed for each pair of retinal layers when differences of any size were considered and when only differences >1 SD were considered ( Table 2). To assess the extent of overlap of thinned regions in any pair of layers L p and L q , the Jaccard index J pq was computed using the formula J pq = C/(A + B − C), where A is the number of thinned points present in L p , B www.nature.com/scientificreports www.nature.com/scientificreports/ is the number of thinned points present in L q , and C is the number of thinned points spatially coincident in both layers (i. e. the number of points where both layers appear thinned at the same spatial coordinates) 10 . Overlap of thickened regions was analyzed in an analogous manner.
The Jaccard index may be understood as the ratio of real to potential spatial coincidence of thinned (or thickened) regions in the two retinal layers being compared. A Jaccard value of 1 means maximum coincidence and 0, maximum incompatibility. Most values in Table 2 are statistically significant, leading to the conclusion that the location of thinned regions in any retinal layer correlated with the location of thinned regions in other layers, and that the locations of thickened regions also correlated across layers. But the sense of the correlation is different for thinned and thickened regions: most values below the main diagonal of the upper matrix in Table 2 are greater than randomly expected, while most values above the diagonal are lower than randomly expected. In other words, thinned regions co-located with one another to a statistically significant extent across layers, while thickened regions "avoided" co-localization to a statistically significant extent across layers (cf. OPL and ONL in Fig. 3c,d).

Discussion
Thickness changes in the 10 retinal layers were studied in patients with AD at a very early stage of disease development. We automatically segmented 10 retinal layers, evaluated their thickness at each retinal point scanned, and searched for spatial patterns of thickness differences between patients and controls. All 10 layers showed AD-related thinning over a relevant percent of their surface, and thinning reached statistical significance at various locations in NFL, GCL, IPL, INL, OSL and total retina. All layers except INL and OS also showed thickened regions, though their thickening did not reach statistical significance. In neural layers thinned regions showed a statistically significant increase in area over that of thickened ones, whereas the opposite was observed in the RPE. Volume lost in thinned regions of neural layers was greater than volume gained in thickened regions, the difference being statistically significant; the opposite was again observed in the RPE. Across different retinal layers, thinned regions appeared to co-localize while thickened regions appeared to "avoid" co-localization, in both cases to a statistically significant extent. Some layers affected by AD in our study differ from those found in some previous studies. Most studies reported a statistically significant thinning in inner layers [11][12][13][14] , but these results were obtained by analyzing all layers separately without any aggregation (e.g. GC-IPL) and by calculating only mean measures of an area. Our results are in agreement with those of Garcia-Martin et al. 13 . Their study and ours detected AD-related thinning of the NFL and ONL, but our method allowed us to go further by building a color-coded quantitative map of  www.nature.com/scientificreports www.nature.com/scientificreports/ thinned areas, which helped us delimit the location of macular thinning in mild AD with greater accuracy than in previous studies. Studies of AD patients with MMSE scores of 17.02 to 19.9 have suggested that inner retinal thinning parallels the decrease in MMSE score, reflecting disease progression [11][12][13][14][15][16] . Interestingly, we were able to detect macular thinning even at a quite early stage of AD in patients with an MMSE score of 23.42.
The thinning found in retinal layers may correspond to loss of ganglion, amacrine and bipolar cells, which may translate to visual function damage in early stages of AD. Loss of these cells means a loss in spatial sampling frequency of receptive fields tuned to all spatial frequencies, which may seriously disturb convolutional computation of all mono-and multi-channel processing. Indeed, this neuronal loss may help explain the observed nearly 40% loss of psychophysical contrast sensitivity in AD patients 6 . A possible explanation for the superior retinal thinning observed in our patients is that amyloid-β (Aβ) deposits and tau tangles predominate in the superior region of AD retinas [17][18][19][20][21][22][23][24] . Indeed, senile plaques and neurofibrillary tangles in AD are densest in the cuneal gyrus of the primary visual cortex, whose axons project from the superior retina 17,23,[25][26][27][28] .
An observation that caught our attention was that in our patients with mild AD, various retinal layers were thicker than the corresponding regions in controls. Our previous work in patients with early AD revealed thickening of the total retina, which we attributed to an early phase of neural inflammation prior to the degenerative process 4 . Another possible explanation for retinal thickening in mild AD is gliosis. Like García-Martín et al. 13 , we detected thickening in the INL, and also in GCL, IPL, OPL and ONL, although this never reached statistical significance in our study. The thickening in the GCL and IPL was distributed such that thicker voxels surrounded foveal voxels that were affected by statistically significant thinning, which highlights the ability of our method to gain detailed insights into the nature of retinal changes in AD.
Previous studies reported a statistically significant thinning of total retinal thickness in AD [11][12][13] , but those studies were limited because they could not identify which retinal layers were most affected. As a result, two of those three studies 11,12 included the ganglion cell layer. García-Martín et al. 13 have also reported on AD-related changes in mean thickness of each retinal layer in the macular area. Their patients had more advanced AD than ours, which may help explain why their observation of thinning in NFL, GCL, IPL and ONL differs from ours. Indeed, retinal thickness is known to correlate directly with MMSE score, with more advanced AD stage associated with more thinning 12,14,15,[29][30][31] . AD progression may also explain why previous studies of retinal layer thickness in AD patients reported thinning but not thickening of inner retinal layers [11][12][13][14][15][16] , whereas we observed thickening in our patients with early AD. Aside from disease differences, some technical factors may also contribute to this  Table 2. Spatial coincidence of thinned and thickened regions across retinal layers in AD patients. Jaccard coefficients for spatial coincidence of thinned regions of each pair of layers are shown below the main diagonals of the two correlation matrices; coefficients for thickened regions are shown above. The matrix on the top has been obtained when the classification of regions as thinned or thickened was based on thickness differences of any size between AD and control groups; and the matrix at the bottom when differences >1 SD were required. Results of statistical significance testing are indicated as follows: -, statistically non-significant; *p < 0.05; **p < 0.01; ***p < 0.001, where these p values represent the probability of family-wise type I error in two-tailed tests for the set of 45 Jaccard values obtained for thinned/thickened regions. discrepancy. One factor is the OCT software used: only the latest-generation OCT allows segmentation of the 10 retinal layers. Another factor is differences in spatial segmentation approaches used to analyze macular thickness. The use of wide regions (e.g. concentric rings), instead of the unconstrained search applied in the present work, may cause thickening and thinning within the same region to cancel each other out. This means that thickness changes may go undetected. In AD, appearance of Aβ in the brain and retina coincides with activation of surrounding microglia. The presence of Aβ in the retina can even precede deposits in the brain 27 . Microglia respond to Aβ deposition and accumulate within and around plaques 32,33 , where they may engage in phagocytosis. Microglia number and size correlate directly with plaque dimensions 18,[34][35][36] , and microglia in AD patients show statistically significant activation 37,38 . The microglial distribution in several retinal layers where we observed thickening leads us to postulate that the thickening could reflect microglial activation, at least in some cases. Activation triggers several changes in microglia that may thicken the retinal layer, including greater branching, process shortening and thickening, migration to the nearest nuclear layer, increase in number and increase in retinal area occupied [39][40][41][42] . In response to insults, microglia can take on macrophage-like morphology, becoming "ameboid microglia" lacking cell processes [43][44][45][46][47] . Ameboid microglia are commonly found in the vicinity of lesions in neuroinflammatory disorders 48 , and they have been reported in neurodegenerative diseases such as multiple sclerosis and AD 49 .
Although they may help to remove Aβ, microglia may actually contribute to degeneration by triggering inflammation. The interaction of Aβ and microglia leads to production of chemokines and neurotoxic cytokines destructive to the central nervous system 50 . Activated microglia can degrade the extracellular matrix, promote the retraction of dystrophic axons and destabilize synapses 51,52 . This so-called "synaptic stripping" 53,54 has been observed in neurodegenerative diseases, such as glaucoma, where retinal ganglion cell synapses are eliminated early in disease progression 41,55,56 . This can cause transneuronal degeneration, which may account for the observed thinning in the GCL, INL and OSL of our patients. Such transneuronal degeneration may help explain the observed thinning in the photoreceptor outer segments in our patients. Another possible explanation is photoreceptor degeneration induced by amyloid aggregation, as seen in animal models 22,57 .
During the last decades, in vivo analysis of the retina from optical cross sections in OCT has improved substantially due to hardware and software advances. Initially, only total thickness of the retina could be examined, from the vitreo-retinal interface to the RPE. Then, it became possible to discriminate between internal and external retina, and now the higher resolution provided by spectral domain OCT combined with software improvements allows separation of up to 11 retinal layers. A key difference between the methods commercially available and ours described here is that commercial methods provide the mean thickness of each retinal layer in geometrically predetermined areas. Our method, in contrast, explores a 6 × 6 mm 2 square at the pixel level without prior constraints, thus avoiding the risk that geometrically predetermined regions may lead to averaging that masks clinically relevant differences. Our spatially unconstrained analysis of retinal layers allows us to identify which layers are truly affected by the disease. Our analysis also applies spatial standardization and correction for distortions coming from OCT techniques and from individual differences, before data are combined for analysis. In fact, the present study appears to be the first to correct thickness values for local tilt of retinal layers. Moreover, RFT-based statistical analysis allows precise identification of the location, extension and shape of the retinal regions in which the two groups differ. Most of these approaches are well established in brain research, but the present work appears to be their first systematic application to OCT studies of retinal thickness. The methods described here may be also useful for other types of retinal studies and other disease contexts.
Our findings are consistent with histopathological studies of human eyes 36 suggesting that changes in retinal thickness may be a consequence of neuronal degeneration and inflammatory processes caused by AD-specific accumulation of retinal deposits that precede the appearance of deposits in the brain.
The detailed insights in this study were made possible by our novel method involving spatial normalization of OCT data and correction of thickness measurements for global, curvature and local tilt in each retinal layer. Moreover the present study is the first, to our knowledge, to apply local tilt correction. Our method is supported by custom-made software that performs statistical analyses based on Gaussian random field and random permutation theories. These statistical methods have four main advantages over classical ones: (1) they allow searches for thickness changes at pixel resolution, giving much finer spatial resolution than previous studies based on predefined grids 14,58 ; (2) they provide precise information about the spatial localization, shape and size of thinned (and thickened) regions in each layer; (3) they avoid masking and biasing results as a consequence of the shape and location of predefined search regions, such as ETDRS or rectangular grids; and (4) they are better suited for multiple comparisons than Bonferroni correction because they count on the spatial autocorrelation inherent in OCT data. Thus the proposed methods may contribute to solving existing controversies and providing more accurate and informative analyses of changes in normal and pathological retinas over time.

Methods
Subjects. In this cross-sectional study, participants were recruited from the database of the Memory Unit of the Geriatric Service of the Hospital Clínico San Carlos (Madrid, Spain). The study protocol followed the principles of the Declaration of Helsinki and was approved by the institutional Ethics Committee for Clinical Research of the Hospital Clínico San Carlos (code number 11/372-E). All participants gave written informed consent.
Review of the records of 2,635 patients allowed us to identify 87 with mild AD, defined as GDS 4 according to the NINCDS-ADRDA Alzheimer's Criteria. These patients underwent a full neurological examination and magnetic resonance imaging of the brain to rule out alternative diagnoses. Their ophthalmic medical records were also reviewed and those who had been previously diagnosed with an ophthalmological disease (glaucoma or suspected glaucoma, media opacity, and retinal diseases) were excluded. This left 29 patients with GDS 4 who were free of ocular disease and systemic disorders that might affect their vision. These 29 patients with mild AD and 37 age-matched healthy control subjects, who scored above 27  www.nature.com/scientificreports www.nature.com/scientificreports/ underwent a complete ophthalmological examination. Six patients and nine controls were subsequently excluded due to posterior pole pathology including macular degeneration, drusen, suspicion of glaucoma, glaucoma, epiretinal membrane, or cataracts that prevented ocular examination. Another four patients and four controls were excluded because signal intensity was inconsistent across the OCT scan. The remaining 19 patients and 25 controls were subjected to a complete ophthalmologic examination conducted by the same clinician. This examination involved assessment of visual acuity, refraction, slit-lamp analysis of the anterior and posterior segments of the eye, applanation tonometry (Perkins MKII tonometer, Haag Streit-Reliance Medical, Switzerland), dilated fundus examination and OCT. All participants showed an AREDS Clinical Lens Standards <2, a best-corrected visual acuity of 20/40, spherocylindrical refraction within ±5 diopters and intraocular pressure below 20 mmHg. After enrollment, one subject of the control group was excluded because automatic layer segmentation succeeded only for a quite small region. In the end, 19 patients and 24 controls were included in the study.
Only the right eye of each subject was studied; however data from the right eye of 4 patients were of insufficient quality, so their left eye was studied instead. In these cases, the OCT data were left-right flipped so that anatomical areas (temporal-nasal) became equivalent for all subjects. Optical coherence volumes were obtained after pupil dilatation using a spectral domain OCT (3D OCT-1000 Topcon, Japan). Three high-quality peripapillary and macular images were obtained in a raster pattern covering a 6 × 6-mm area with a scan density of 512 × 128 pixels in ∼2.5 sec (27000 A scans/sec) 59 . Voxel size was 11.71875 × 46.875 × 3.5 µm 3 (horizontal x vertical x depth), according to the calibration provided by the manufacturer. All OCT images were acquired by the same experienced technician, always keeping the light beam entry point centered on the pupil to avoid an oblique scanning artifact. Images were reviewed for quality, and the criteria for acceptable fundus images were as follows: (A) no large eye movements, defined as an abrupt shift completely disconnecting a large retinal vessel; (B) consistent signal intensity across the scan; and (C) no black bands (caused by blinking) throughout the examination. In addition, the criteria for acceptable scanning were (D) a signal-to-noise ratio >30 and acceptance of >95% of A-scans during fast nerve fiber layer (NFL) scanning. to allow thickness measurements of all patients or all controls to be combined before comparing the two groups with each other. (F) The mean thickness map of each retinal layer was compared between patients and controls using custom-designed software that performed parametric tests based on Gaussian RFT and non-parametric tests based on random permutations. These statistical techniques are well established in brain research and were adapted for the first time to the retina in the present study.

corrections of raw thickness values obtained in the A-scan direction. When the A-scan direction
was not orthogonal to the layer due to retinal tilt and refraction effects, then raw measurement T' of retinal layer thickness (distance between the layer's delimiting surfaces in the A-scan direction) was greater than the true thickness value T.
As a result of retinal tilt with respect to the A-scan direction, raw thickness was not measured perpendicularly to the retinal layer, causing the measured thickness T' to exceed the true value T. There are at least three types of tilt, depending on its origin. One type is global tilt (Fig. 6a), which is mainly due to pupil eccentricity at the entry of the OCT beam 63 . This tilt is constant over the scanned area and reached 10.0° in our subject 24, a value quite close to 11° reported by Antony et al. 64 . A second type of tilt is spherical tilt (Fig. 6b), which is due to natural curvature of the eye. This curvature causes retinal tilt to increase radially, from zero near the center of the scanned area up to 7.0° at its corners, as in subject 10 in our study. The third type of tilt is local tilt (Fig. 6c), which is due to variation of the local orientation of each retinal layer. It varies at each point (x, y) of the scanned area and reached 14.4° in the NFL very close to the fovea center in subject 23 in our study.
Global tilt was modelled for each subject using the vector normal to a plane fitted to points defining the outer surface of the RPE. Spherical tilt was modelled using the vector normal to the sphere surface (radius) resulting from the revolution of the circle fitted to the 512 points of the OB_RPE corresponding to the B-scan crossing the center of the scanned area. Local tilt at point (x, y) of a retinal layer was modelled using the normal vector of the plane fitted to a 0.5 × 0.5 mm 2 region of the medial surface of the layer in the flattened OCT, centered at (x, y). The flattened OCT was obtained through z-axis displacement of all A-scans so that z-coordinates of the outer delimiting surface of RPE became zero (or some other constant). The medial surface of a retinal layer was defined as the surface dividing the layer into two equally thick sub-layers. The z-coordinate at any scanned point (x, y) was obtained by averaging the z-coordinates of the same point on the two delimiting surfaces. The angle α between the A-scan direction (z-axis in OCT) and the vector sum of the three normal vectors was used to compute the true thickness T at each point in a retinal layer as T(x, y) = T'(x, y) cos α(x, y). Applying an additional correction for refraction effects was considered unnecessary because of the automatic segmentation algorithm. Spatial normalization. Thickness maps in the x-y plane for each subject had to be first spatially normalized in order to allow thickness values from different subjects to be combined. This ensured that the values being averaged together referred to the same anatomical region. Since OCT lacks standardized references analogous to brain maps in functional MRI, all OCT data were spatially normalized through translation, rotation and scaling. The thickness matrix of each subject was displaced in the x-y plane so that the fovea center c m came to the central point www.nature.com/scientificreports www.nature.com/scientificreports/ of the image matrix (Fig. 7). The matrix was rotated so that the maculopapilar axis s of all subjects overlapped at an angle θ = 6.766°, and scaled so that all papillary centers c p overlapped at 4.377 mm from the common fovea center c m . These values of θ and s were the mean values in the entire study sample.

Statistical analysis.
To investigate which areas of each retinal layer were altered in patients relative to controls, matrices of raw thickness values with dimensions of 128 × 512 were resized to 512 × 512 using bilinear interpolation, then smoothed with an isotropic Gaussian filter whose variance was chosen based on theoretical and practical considerations. The theoretical consideration was that the filter's full width at half maximum (FWHM) should be >4 times the voxel size in its lower-resolution dimension 65 , which in our case meant >187.5 µm. The practical consideration was that we wished to improve on the spatial resolution of previous studies, but without gaining a resolution so fine that might be spurious or inaccurate; retinal signals with a high spatial frequency may not be detected with sufficient signal-to-noise ratio with our OCT equipment. These considerations led us to choose FWHM = (1, 1) mm, giving a resel size of 1 × 1 mm 2 and 36 resels for the entire scanned area.
The area searched for differences in retinal thickness between patients and controls was restricted to a central circle of the scanned region, because data were not available outside that region for all subjects. This happens mainly because OCT rotation and scaling during spatial normalization generated empty peripheral regions and because the automatic layer segmentation failed in some peripheral regions of a small number of subjects, thus rendering layer thickness unavailable at those regions. The area searched for differences was therefore restricted to a fovea-centered square of 4.8 × 4.8 mm 2 . The search region for RFT-based analysis was restricted to a fovea-centered circle of diameter 3.5 mm. This region contained data from all subjects and covered the fovea and parafovea. Restricting the search region to this circle implied the resel count to be set at 17 resels 66 .
The family-wise null hypothesis can be tested in a variety of ways. We chose "height thresholding", because voxels inside the search region with above-threshold signal allow us to conclude the existence of statistically significant thinning in those voxels 67 . Other hypothesis-testing methods, in contrast, do not provide spatial localization of effects. Family-wise error (FWE) probability under the null hypothesis was set to α = 0.05, and one-sided statistical tests were adopted. The total number of subjects was 43. The statistical parametric map for the macular region was built on the z statistic using a standard Gaussian distribution, and its critical value was computed for each layer according to RFT. For those layers whose data failed to satisfy the requirements of RFT, Bonferroni correction was discarded because nearby voxels could not be treated as mutually independent; we used instead the non-parametric permutation test 68 . This statistical testing was coded into in-house software to identify thinned regions in each retinal layer. Two versions of the software were independently developed and then cross-checked to ensure that results were reliable.
Spatial overlap of thinned/thickened regions of different layers was evaluated using Jaccard coefficients and hypothesis-tested using exact probability tests 10,69 . A Jaccard value was considered statistically significant when the probability of obtaining equal or more extreme values was below 0.05/(2 × 45) in a two-sided test with Bonferroni correction for the 45 layer pairs being compared. Given the extremely large integer numbers generated by combinatorial coefficients in probability computations, we used the Symbolic Math, Variable Precision Integer Arithmetic and High Precision Floating Point Arithmetic toolboxes in Matlab R2018a 70 . Determining the statistical significance of each Jaccard coefficient required in most cases approximately 19 h of computing time on a dedicated server with an Intel Xeon processor E5-2690 v3, 12 core, 2.6 GHz, 35 MB, and 32 GB RAM. Figure 7. Parameters for spatial normalization: center of the fovea c m , center of the papilla c p , length s and angle θ of maculopapillar axis. All OCTs were moved, rotated and scaled so that their macular and papillar centers overlapped. Data on macula and papilla were acquired independently in two OCT sessions, so their coordinate systems could not be directly related to each other. Thus, the two datasets had to be co-registered based on shared retinal regions before maculopapilar axis length and angle could be computed.