Longitudinal Analysis of Retinal Ganglion Cell Damage at Individual Axon Bundle Level in Mice Using Visible-Light Optical Coherence Tomography Fibergraphy

Abstract

Purpose: We developed a new analytic tool based on visible-light optical coherence tomography fibergraphy (vis-OCTF) to longitudinally track individual axon bundle transformation as a new in vivo biomarker for retinal ganglion cell (RGC) damage.

Methods: After acute optic nerve crush injury (ONC) in mice, we analyzed four parameters: lateral bundle width, axial bundle height, cross-sectional area, and the shape of individual bundles. We next correlated the morphological changes in RGC axon bundles with RGC soma loss.

Results: We showed that axon bundles became wider and taller at three days post ONC (pONC), which correlated with about 15% RGC soma loss. At six days pONC, axon bundles showed a significant reduction in lateral width and cross-sectional area, followed by a reduction in bundle height at nine days pONC. Bundle shrinking at nine days pONC correlated with about 68% RGC soma loss. Both experimental and simulated results suggested that the cross-sectional area of individual RGC axon bundles is more sensitive than bundle width and height to indicate RGC soma loss.

Conclusions: This study is the first to track and quantify individual RGC axon bundles in vivo after ONC injury.

Translational Relevance: Recognizing RGC loss at its earliest stage is crucial for disease diagnosis and treatment. However, current clinical methods to detect the functional and structural changes in the inner retina are not sensitive enough to directly assess RGC health. In this study, we developed vis-OCTF-based parameters to track RGC damage, making possible to establishing a quantifiable biomarker for glaucoma.

Introduction

Retinal ganglion cell (RGC) loss is a hallmark of optic neuropathies, such as glaucoma,^1–3 and neurodegenerative diseases that affect vision, such as Parkinson's and Alzheimer's disease.^4,5 Thus recognizing RGC loss at its earliest stage is crucial to prevent further irreversible vision loss.^6–8 However, current clinical diagnostic methods to detect the functional and structural changes in the inner retina are often not sensitive or specific enough to directly assess RGC health.^9,10 An essential parameter for glaucoma diagnosis is visual field (VF) loss; however, studies suggest that more than 30% of RGCs may be lost before detection on VF testing.^11–15 Characteristic changes in the optic nerve head (ONH) and optic disc are also used for glaucoma diagnosis, but identification of damage can be subjective, and grading varies between observers.¹⁶ Among the clinical noninvasive imaging modalities,^17–20 optical coherence tomography (OCT) is most widely used for diagnosing and monitoring optic neuropathies. OCT's cross-sectional imaging capabilities enable measurement of the retinal nerve fiber layer (RNFL) and ganglion cell–inner plexiform layer (GCIPL) thickness in vivo as indirect indicators of RGC health.^14,21,22 However, inconsistent axial resolution (5 µm to 10 µm) and segmentation algorithm discrepancy among different clinical devices make estimating RGC health less reliable.^23,24 In conventional RNFL analysis, algorithm failure due to edge smoothing may also cause underestimation of local areas of thinning when they do occur over time.²⁵ 

Importantly, the clinical parameters measured by OCT, such as the RNFL or GCIPL thinning, are not specific indicators for RGC damage. First, a significant variability in RGC density and RNFL thickness exists among individual patients. For example, the RNFL thickness of healthy human subjects varies from about 50 µm to 120 µm as measured by OCT.²⁶ Overlaps of RNFL ranges were also observed between healthy subjects and glaucoma patients.²⁶ Shin and colleagues²⁷ followed up 292 eyes from 192 patients with primary open-angle glaucoma. They found that only 72 eyes (24.7%) showed progressive GCIPL thinning, and among the 72 eyes, 41 eyes showed progressing visual field (VF) loss. The different patterns of disease progression thus emphasize the need to monitor individual patients longitudinally.^22,27 Second, the ganglion cell layer (GCL) contains non-RGC cells, called displaced amacrine cells, which vary from 30% to 80% in the GCL and are largely unaffected by glaucoma.^2,28 We observed 38% RGC axon loss at three to five days (d) post optic nerve crush (pONC) injury, but at the same time, we found no significant reduction in overall RNFL thickness in mice.^29,30 In addition, the inner plexiform layer (IPL) contains synaptic connections among different retinal cell types in rodent and primate eyes.^28,31,32 In other words, the RNFL or GCIPL thinning may not be sensitive enough to detect RGC damage specifically. Therefore there is a need for more accurate and sensitive biomarkers for RGC damage following the disease insult. 

We recently applied visible-light OCT (vis-OCT), which operates from 510 nm to 610 nm and reaches an axial resolution of 1.3 µm, in the mouse retina.^33,34 Because optical scattering contrast exponentially decays with increasing wavelength, the visible-light spectral range offers a much higher tissue contrast than the near-infrared (NIR) range.³⁴ As a result, vis-OCT provides unique anatomical and functional imaging capabilities that facilitate RGC damage evaluation.³³ Recently, we developed vis-OCT fibergraphy (vis-OCTF) to visualize and quantify changes in individual RGC axon bundles in mice.^34–36 In this study, we applied vis-OCTF to track the changes of axon bundles after ONC injury by quantifying four parameters: (1) lateral bundle width, (2) axial bundle height, (3) cross-sectional area, and (4) bundle shape. We correlated the in vivo changes in RGC axon bundle morphology with RGC soma loss after ONC injury. Based on the experimental data, we created a numerical simulation for each parameter to determine the sensitivity and the floor value—a threshold at which no further change is observable.^37,38 

Methods

Healthy adult (3–6 months) male and female wildtype C57BL/6 mice were used for this study. All animal protocols were approved by the University of Virginia institutional animal care and use committee and complied with the guidelines of National Institutes of Health and the Association for Research in Vision and Ophthalmology. 

Optic Nerve Crush Surgery

The ONC procedure was performed as described previously.^35,39 Briefly, mice were anesthetized with an intraperitoneal injection of a mixture of ketamine (100 mg/kg, Kataset, Zoetis; NADA no. 043-304) and xylazine (8 mg/kg, AnaSed, Akorn; NADA no. 139-236). A small incision was made in the superior and lateral conjunctiva, and the optic nerve was exposed by gentle dissection. The optic nerve was then gently clamped with a pair of forceps approximately 1 mm behind the globe for seven to 10 seconds. After surgery, moxifloxacin (0.5%, NDC 60505-0582-4; Apotex Corp., Toronto, Canada) was applied to the crushed eyes to prevent infection. Mice were kept on a heating pad until fully recovered. 

Vis-OCT Animal Preparation

Before imaging, mice were anesthetized by intraperitoneal injection of a mixture of ketamine (114 mg/kg, Kataset; Zoetis, Parsippany-Troy Hills, NJ, USA) and xylazine (17 mg/kg, AnaSed; Akorn Pharmaceuticals, Gurnee, IL, USA). The pupils were dilated using tropicamide drops (1%; Henry Schein Animal Health, Covetrus, Portland, ME, USA). During imaging, mice were kept warm with an infrared heat lamp and given polyvinyl alcohol artificial tear drops (1.4%; Rugby Laboratories, Inc., Hempstead, NY, USA) to prevent corneal dehydration. After imaging, mice were recovered on a heating pad and monitored until alert and active. 

Visible-Light Optical Coherence Tomography

Mice were imaged using a small-animal vis-OCT system (Halo 100; Opticent Health, Evanston, IL, USA) as previously reported.^34,35 For each eye, four to six volumes were acquired with the ONH aligned in each corner of the field of view (FOV) to cover the entire retina.³⁵ Each volume took ∼10.5 seconds to acquire, requiring ∼5–7 minutes to capture the whole eye including time for repositioning and artificial tear application between acquisitions. 

Vis-OCT Fibergram Processing

Vis-OCT fibergrams were extracted from each image volume.³⁴ We used an intensity-based threshold method to detect the surface of the retina and cropped the RNFL by selecting the first ∼15 µm in depth. We then calculated the mean intensity projection along the axial (z) direction to generate the fibergram composed of RGC axon bundles and surrounding vasculature. The fibergrams were then montaged, creating a final FOV of approximately 1.2 mm × 1.2 mm.³⁵ 

Each vis-OCT volume was digitally resampled to generate a circumpapillary B-scan with a radius of 425 µm. To do so, we manually marked the ONH in the en face image and plotted a ∼15 µm-thick arc around the ONH. The pixels were then sorted as a function of the angle measured between each sampled A-line and the nasal direction with the ONH as the vertex. Adjacent A-lines within a 0.1° sector were averaged to reduce speckle noise while preserving spatial density. 

Individual Axon Bundle Size Quantification

We used the blood vessel pattern and the ONH as reference points to identify and track individual axon bundles in each retina (Supplementary Fig. S1). We measured the axial height (H_b) and lateral width (W_b) of the individual RGC axon bundles using MATLAB. Bundle width measurements were recorded as previously reported. Briefly, the center axis of each bundle was manually marked, and the mean intensity profile along the center axis was plotted as shown in Figure 1C. The intensity profile was normalized between 0 and 1, and the bundle width was recorded as the profile width at 1/e². Bundle height was similarly recorded by extracting the axial intensity profile of the bundle and recording the height value at 1/e², as shown in Figure 1D. Measurement values were reported for individual bundles. Blood vessels were excluded from analysis by identifying the dark shadows in the B-scan images and uniquely distinguishable branching structures compared with surrounding axon bundles in fibergram images. 

Figure 1.

 

In vivo identification and quantification of RGC axon bundle morphology. (A) A fibergram from a single OCT volume of a wildtype mouse. Two RGC axon bundles (1) and (2) were labeled at the radius of 425 µm from the ONH as indicated by the red arc. (B) Circumpapillary B-scan image reconstructed along the red arc in A shows the cross-sectional image of the retina. The blue arrow in A indicates the leftmost A-line in B. The green dashed lines indicate the bundles (1) and (2). (C, D) Intensity profile of the bundle width (C) and the bundle height (D) of the bundles (1) and (2). The intensity profile width is measured at 1/e² decay, as indicated by the black (1) and red (2) double headed arrows. (E) Table of the lateral width, bundle height, cross-sectional area, and the shape indicators of bundles (1) and (2). INL, inner nuclear layer; OPL, outer plexiform layer; ONL, outer nuclear layer; RPE, retinal pigment epithelium.

Bundle cross-sectional area was approximated by the area of an ellipse as follows:  

\begin{eqnarray*}{A_b} = \frac{\pi }{4}{W_b}{H_b},\end{eqnarray*}

where A_b is the cross-sectional bundle area. The elliptical area approximation was used throughout this study instead of pixel-based measurements because of segmentation errors that occur when neighboring bundles are packed too closely together. To justify this approximation, we performed pixel-based measurements on isolated bundles and compared the results with the approximated area measurements. First, isolated bundles were manually selected from the resampled B-scans. Cropped images were contrast enhanced and median filtered to reduce noise. The filtered images were then binarized and morphologically opened. A watershed transform was applied and basins outside the selected bundle were removed. The pixel area of the segmented cross section was then calculated. Finally, we compared the cross-sectional area measurements approximated by the area of an ellipse with pixel-base measurements of the same bundles and found no significant difference between the two methods (See Results). 

We also defined a shape parameter to describe the width-to-height ratio of individual bundles using a dimensionless value normalized between −1 and 1. The shape parameter is defined as:  

\begin{eqnarray*}{S_b} = \left( {1 - \left( {\frac{{\min \left( {{W_b},{H_b}} \right)}}{{\max \left( {{W_b},{H_b}} \right)}}} \right)} \right){\rm{*}}\left( {\frac{{{W_b} - {H_b}{\rm{\;}}}}{{\left| {{W_b} - {H_b}} \right|}}} \right),\end{eqnarray*}

where S_b is the shape parameter of a single RGC axon bundle. A positive S_b suggests a wider axon bundle elongated laterally, and a negative S_b suggests a taller bundle elongated axially. 

To calculate the axon bundle density, we tracked the same area of each retina before and after ONC injury. We manually counted the number of bundles located at the radius of ∼425 µm in the fibergram images and calculated the bundle density per mm. We also measured the GCIPL thickness, which was a measure of the top edge of the vitreous/RNFL and the bottom edge of the IPL.³⁵ 

Immunohistochemistry and Confocal Imaging

After acquiring vis-OCT images, mice were euthanized with euthasol (600 mg/kg, Euthasol, Virbac ANADA, no. 200-071; Virbac, Carros, France) and perfused with paraformaldehyde (4%, ChemCruz, sc-281692; Santa Cruz Biotechnology, Dallas, TX, USA). Eye cups were dissected, post-fixed in paraformaldehyde for 30 minutes, washed with phosphate-buffered saline solution containing Triton-X detergent (PBST, 0.5% Triton X-100), and then blocked for one hour in blocking buffer (1% BSA and 10% normal donkey serum, 0.5% Triton X-100; Sigma-Aldrich, St. Louis, MO, USA). Primary antibodies, diluted using blocking buffer, included rabbit anti-rbpms (ab194213, 1:500; Abcam, Cambridge, MA, USA), mouse anti-neurofilament H antibody (MCA1321GA, 1:250; Bio-Rad Life Science, Hercules, CA, USA), mouse anti-Tuj1 (gift from Tony Spano, University of Virginia, 1:250), and rat anti-Icam-II (553325, 1:500; Pharmingen, BD Bioscience, Franklin Lakes, NJ, USA). Secondary antibodies, including donkey anti-mouse immunoglobulin G conjugated to Alexa Fluor 594 dye (A-21203, RRID: AB_141633; Invitrogen, Carlsbad, CA, USA) and donkey anti-rabbit immunoglobulin G conjugated to Alexa Fluor 488 dye (A-21206, RRID: AB_2535792; Invitrogen) were diluted at 1:1000 in blocking buffer and incubated overnight at 4°C. After immunostaining, retinas were flat-mounted and cut into four quadrants: temporal (T), nasal (N), inferior (I), and superior (S). The blood vessel pattern was used as a landmark to align vis-OCTF and confocal images of flat-mounted retinas (Supplementary Fig. S1). 

Confocal images were taken using a Zeiss LSM800 confocal microscope (Zeiss, Thornwood, NY, USA). Z-stack images covering the depth of the outer nuclear layer to the GCL (approximately 50–80 µm) were acquired. Lower-magnification (5×) pictures were captured for the whole retina using the tiling/stitch function in Zen (Zen 3.2; Zeiss, Oberkochen, Germany). For cell counting, individual images were captured at magnification of 10×, covering an area of 0.408 mm². To show morphological changes in degenerating axons, individual z-stack images covering an area of 0.0037 mm² were taken using a water immersion objective 63×. 

RGC Soma Quantification

We followed our published protocols to quantify RGC soma density.^40–43 In brief, the nasal side of the eye was marked during eye dissection for orientation. To quantify RGC density, mouse retinas were immunostained with anti-rbpms antibody. For each retina, 16 enface z-stack images covering the depth of the GCL were captured (20–30 µm of total depth). Four images were acquired for each of the four quadrants to ensure broad coverage of the entire retina. Three rectangles covering no less than 0.03 mm² area were randomly drawn on the images, avoiding overlaps with large blood vessels. All rbpms-positive cells within the rectangles were manually counted in Zen (Zen 3.2; Zeiss).⁴³ Cell density was calculated using total cell counts divided by area for the four quadrants in each retina. All data analyses were cross-examined by two independent observers. 

Statistical Analyses

All statistical analyses were performed using MATLAB and Prism. We used a linear mixed-effects model for all width, height, area, and shape comparisons to remove influence from individual subjects. Two sets of ONC experiments were performed. One set was to track changes in the axon bundle parameters following the ONC of the same retina. The other set was to correlate changes in axon bundle morphology with RGC soma loss using different retinas. For both sets of experiments, we performed one-way analysis of variance (ANOVA) followed by Dunnett test for multiple comparisons. A significance level of 0.05 was used, and the P values of the Dunnett tests were reported unless otherwise stated. All results were reported as mean ± standard deviation. 

Results

Establishing a New Analytic Tool for Identifying and Tracking Individual Axon Bundles In Vivo

Taking advantage of the improved axial resolution offered by vis-OCTF and associated data processing methods,³⁴ we established an analytic tool for tracking change in individual RGC axon bundle size rather than bulk RNFL thickness (Fig. 1). Figure 1A shows one of the four processed vis-OCT fibergrams (vis-OCTF) from an adult wildtype mouse. We measured the RGC axon bundle lateral width using the vis-OCTF image, which has been validated by confocal imaging.^34–36 As shown in Figure 1C, we plotted the signal intensity profile of the bundle width, which gives 24.9 µm for bundle 1 and 11.7 µm for bundle 2. Next, we measured the RGC axon bundle height using the resampled circumpapillary B-scan image.^35,36 The red arc shows the path reconstructed at 425 µm radius from the ONH (Fig. 1B). The axial bundle height was also measured using the signal intensity profile, which gave 19.7 µm for bundle 1 and 12.8 µm for bundle 2, respectively (Fig. 1D). 

We quantified the cross-sectional area using the pixel-based cross-sectional area and cross-sectional area approximated by the area of an ellipse (See Methods). The mean of 57 individual axon bundles measured by pixel area was 158.1 ± 69.1 µm², consistent with the mean approximated by the area of an ellipse (160.0 ± 82.3 µm², P = 0.90, Student's t-test). Thus we used the area of an ellipse to approximate bundle cross-sectional area throughout this study. As summarized in Figure 1E, the area of the axon bundles 1 and 2 are 386 µm² and 117 µm², respectively. In addition, we developed a dimensionless indicator for the shape of bundle (S_b), which normalizes the bundle width to height ratio between −1 and 1 (see Method section for more details). A wider axon bundle, such as bundle 1, has a positive shape value (S_b1 = 0.21); whereas a taller axon bundle, such as bundle 2, has a negative shape value (S_b2 = −0.08). 

In Vivo Tracking of Morphological Changes in RGC Axon Bundles After ONC

We applied vis-OCTF to examine the structural changes at the single axon bundle level following ONC injury. Figure 2A shows example fibergrams and circumpapillary B-scans of a mouse retina acquired at baseline (before ONC) and three, six, nine, and 15 days post ONC (pONC). The red arc indicates the path of the resampled circumpapillary B-scan at the radius of 425 µm from the ONH. The blue boxes in the middle panel of Figure 2A show the magnified views of the fibergrams with five axon bundles tracked over time. The right panels of Figure 2A show the circumpapillary scans resampled at the same location at different time points with bundle profiles traced in green. 

Figure 2.

 

In vivo tracking of individual axon bundles after ONC injury. (A) In vivo images of the same retina at baseline and 3-days, 6-days, 9-days and 15-days pONC. Left panel shows the vis-OCTF images at different time points. Middle panels show the magnified views of highlighted regions in left panels (blue boxes). Right panel is the resampled circumpapillary B-scans at the radius of 425 µm (red arc). The green points (1-5) indicate the same five axon bundles tracked over time. V₁ and V₂ indicate blood vessels. (B) Quantification plots of the width, height, area, and shape (S_b) of the five tracked axon bundles over time.

At baseline, the five axon bundles exhibit diverse width, height, area, and shape (Fig. 2B). After ONC injury, the width of bundle 1 (blue) increased from 23.5 µm (before ONC) to 25.0 µm at three days, and then progressively decreased at six days pONC (19.9 µm). In contrast, the width of bundle 4 (green) decreased from 23.6 µm (before ONC) to 21.3 µm at three days and continued to decrease at six days pONC (20.4 µm). The height also changed following different patterns. For example, the height of bundle 1 increased from 11.6 µm (before ONC) to 20.3 µm at three days and 20.6 µm at six days pONC and then decreased at nine days pONC. The height of bundle 4 increased from 10.3 µm (before ONC) to 16.8 µm at three days pONC and then decreased to 15.7 µm at six days pONC. The cross-sectional bundle area presented an overall clear pattern of increase-to-decrease from bundles 1 to 4, although the peak point for each bundle is different. Bundle 5, on the other hand, remained stable until nine days pONC (Fig. 2B). The shape, however, did not show a clear trend for the selected bundles. 

We identified and tracked 141 axon bundles from three mice after ONC injury. We plotted the histograms of the four parameters (Supplementary Fig. S2), and their fitted distributions are shown in Figure 3. The RGC axon bundle width significantly increased from baseline (black, mean: 13.7 ± 4.6 µm) to three days pONC (dark red; 14.5 ± 4.5 µm; P = 4.6e-2, one-way ANOVA), followed by a decrease at six days pONC (red; 12.8 ± 4.3 µm, P = 5.9e-2) and nine days pONC (light red; 11.5 ± 3.8 µm, P < 1e-4, Fig. 3A). Similarly, the height distribution plot shown in Figure 3B indicates a significant increase in bundle height at three days (black; baseline mean: 12.5 ± 3.9 µm; dark blue; three days pONC: 13.5 ± 4.5 µm, P = 1.5e-2). However, the height returned to baseline at six days after ONC (blue; 13.1 ± 4.3 µm, P = 2.5e-1), followed by a significant decrease at nine days pONC (light blue; 11.1 ± 4.0 µm, P = 5.0e-4). 

Figure 3.

 

Distribution of the changes in RGC Axon Bundle morphology after the ONC injury. (A–D) Smoothed distributions of the lateral width (A), axial height (B), cross-sectional area (C), and shape (D) for 141 axon bundles (n = 3 mice) tracked over time pONC. Shaded arrows point to the mean value of the distribution curve. (E) Percent change of RGC axon bundle width (red), height (blue), and area (purple) with respect to the baseline values. (F) Percentage of axon bundles exhibiting increased (gray) or decreased (orange) cross-sectional area at different times compared to baseline.

Figure 3C shows the distribution of bundle area. At three days pONC (dark purple), axon bundles presented larger size compared to baseline (black; baseline: 136.6 ± 71.7 µm²; three days pONC: 158.6 ± 82.9 µm²; P = 9.0e-4). At six days pONC, the area returned to baseline (purple; 6-d pONC: 133.8 ± 71.3 µm²; P = 9.1e-1), whereas at nine days pONC the area of the axon bundles was significantly smaller compared to the baseline area (light purple; 9-d pONC: 103.8 ± 58.9 µm², P < 0.0001). We also observed that large bundles tend to elongate axially immediately after the ONC injury (blue arrows in Supplementary Fig. S2). The distribution curves of the axon bundle shape show a progressive shift toward the negative values from baseline (black; 0.065 ± 0.3) to three days (dark green; 0.067 ± 0.3, P = 9.9e-1), six days (green; 0.011 ± 0.3, p = 2.8e-2), and nine days pONC (light green; −0.037 ± 0.3, P = 4.8e-2, Fig. 3D), suggesting that axon bundles naturally have a broader shape that begin to shrink to a more circular shape after injury. 

We calculated the percent change of width, height, area, and shape of the 141 tracked axon bundles (Figs. 3E, 3F). The width of the axon bundles (red) was 14% higher at three days pONC and only 0.5% higher at six days pONC compared to baseline. The height of the axon bundles (blue) was 13% higher at three days and 9% higher at six days pONC compared to baseline. By nine days pONC, the width and height of RGC axon bundles decreased by 10% and 7%, respectively. In other words, the width dropped below baseline close to six days pONC, whereas height dropped below baseline between six and nine days pONC. The cross-sectional bundle area, which enhances the subtle changes in the axon bundle morphology, reveals a clear trend (Fig. 3E). The area increased by 30% at three days, and 8% at six days, and then decreased by 16% at nine days pONC, compared to baseline. We also observed a reduction in bundle shape from baseline by 7% at three days, then to 32% at six days followed by 13% at nine days (Fig. 3E). Figure 3F shows changes in size in the overall axon bundle population. About 60% of bundle population presented a size increase at three days pONC, and 70% of the bundle population showed a size decrease at nine days pONC. At six days pONC, about half of the bundle population increased (51%) and half of the bundle population decreased (48%). 

We next compared the GCIPL measurements with our new axon bundle parameters (Fig. 4). In mice, the GCIPL includes RNFL, GCL, and IPL as indicated by the blue arrows in Figure 4A. We quantified the GCIPL thickness from the circumpapillary B-scan reconstructed at the radius of 425 um (Fig. 4A). The smoothed distribution of the GCIPL measurements (n = 120) shows that the GCIPL thickness was not affected at three days pONC compared to the baseline (baseline, black, 60.3 ± 5.2 µm; three days pONC, dark blue, 61.7 ± 6.1 µm; P = 0.26; Fig. 4C). Whereas at six days pONC (blue) and nine days pONC (light blue), the GCIPL thickness was significantly decreased to 56.0 ± 4.8 µm (P < 0.0001) and 51.8 ± 4.5 µm (P < 0.0001), respectively, compared to the baseline. 

Figure 4.

 

No significant change in the GCIPL height was detected at three days pONC by vis-OCT imaging. (A) Resampled circumpapillary B-scans from the same location from at different time points. Blue arrows indicate GCIPL axial height. (B) In vivo fibergram images of the same retina at baseline and three, six, and nine days pONC. V₁ and V₂: blood vessels. Green dots indicate axon bundles between V₁ and V₂ at the radius of 425 µm (red arc). (C) Smoothed distribution of the GCIPL axial height measurements. (D) Average of the number of RGC axon bundle per mm at baseline and three, six, and nine days pONC.

In addition, we tracked the change in bundle density (Fig. 4D). Before ONC, the density of RGC axon bundles was 60 ± 9 bundles/mm. After ONC, there was not a significant change in density (three days pONC, 56 ± 9 bundle/mm; six days pONC, 63 ± 10 bundles/mm; and nine days pONC, 58 ± 10 bundles/mm; P > 0.5, one-way ANOVA). 

Confocal Microscopy Imaging Confirmed Axon Bundle Damage After ONC

We performed confocal microscopy to validate the axon bundle damage following ONC injury. The retina was dissected and double-immunostained with mouse anti-tubulin beta 3 (Tuj1) and rat anti-neurofilament H (Supplementary Figs. S3, S4).⁴⁴ We confirmed that somas of degenerating RGCs became labeled by anti-neurofilament H at three days pONC (top panel, Supplementary Fig. S3). Confocal images in Supplementary Fig. S4 were taken at the ONH region close to the lesion site of the ONC surgery. In the control retina, axon bundles were well organized and directly converged toward the ONH (Supplementary Fig. S4A). At three days pONC, some axon bundles became disorganized with loose or splitting fibers (yellow arrows, Supplementary Fig. S4B). We also observed retraction bulbs, the non-growing counterparts of growth cones, at the tip of the lesioned axons (white arrows).⁴⁵ At nine days pONC, the axon bundles near the ONH became entangled, with some retraction bulbs and lesioned axon tips pointing away from the ONH (Supplementary Fig. S4C). At 15 days pONC, the overall bundle structure degenerated at the ONH (Supplementary Fig. S4D). 

Establishing the Correlation of Morphological Changes in Axon Bundles With RGC Soma Loss

We next determined the correlation between the in vivo vis-OCTF parameters and the RGC soma loss by ex vivo confocal imaging. For this set of experiments, we acquired vis-OCT images from healthy wild type mice and mice at three, nine, and 15 days pONC, respectively. Immediately after vis-OCT imaging, mice were sacrificed and perfused, and the retinas were dissected and immunostained with rbpms, an RGC marker.⁴³ Flat-mounted retinas were then imaged by confocal microscopy.^35,43 As demonstrated by the flattened retina schematic in Figure 5A, we divided the retina into superior (S), inferior (I), nasal (N), and temporal (T) quadrants. Examples of in vivo vis-OCTFs (left panel) and their corresponding ex vivo confocal microscopy images of flattened retinas (middle panel) stained with rbpms are shown side by side in Figures 5B–E. Magnified views of retinal confocal images show a steady decrease in RGC soma density compared to control at three, nine, and fifteen days pONC (right panel in Figs. 5B–E). 

Figure 5.

 

Confocal imaging of flat mounted retinas for quantification of RGC loss after vis-OCT imaging. (A) Schematic representation of the flat mounted retina with vis-OCT (blue box) and confocal microscopy (orange boxes) FOVs overlaid. (B–E) Vis-OCTF (left panel) followed by confocal (middle panel) images of RGCs labeled by rbpms for estimating RGC density in control (B), three days (C), nine days (D), and 15 days (E) pONC eyes. Right panels show the magnified views of highlighted regions in middle panels (orange boxes). Green dashed lines in vis-OCTF divided the field of view into four regions: superior (S), nasal (N), temporal (T), and inferior (I).

We plotted the rbpms⁺ RGC density of each retina as a function of each of the four axon bundle parameters quantified from the vis-OCT images (Figs. 6A–D). Each data point represents the average reading per retina. The mean density of rbpms⁺ RGCs in controls was 4095 ± 209 cells/mm² (n = 6 retinas). At three days pONC, the mean density was decreased to 3475 ± 343 cells/mm² (n = 9 retinas), a 15% reduction on average. The mean density of rbpms continued to decrease with time: 1305 ± 104 cells/mm² (68.1% reduction, n = 8 retinas) at nine days pONC, and 611 ± 37 cells/mm² (85.1% reduction, n = 3 retinas) at fifteen days pONC. 

Figure 6.

 

Correlation of RGC soma density and axon bundle size measurements. (A–D) RGC axon bundle cross-sectional area (A), width (B), height (C), and shape (D) plotted as a function of rbpms + RGCs density for each retina. Black, dark gray, light gray, and white dots represented data from control (N = 6), three days (N = 9), nine days (N = 8), and 15 days (N = 3) pONC, respectively. Second-order polynomial regression (dotted line) were fitted to the data for A–C, and linear regression were fitted to D. The equation and R² values are labeled on the figure.

The overall relationship between axon bundle height, width, and area with RGC density was non-linear (Figs. 6A–C), because of a swelling phase that immediately followed ONC. For example, at three days pONC (Fig. 6A), a 15% of RGCs lost (P < 0.001) correspond to 36% increase of the axon bundles area compared to the control (control: 92.4 ± 13.6 µm², n = 6; three days pONC: 126 ± 17.8 µm², n = 9; P < 0.01). At nine days pONC, 68% of RGCs lost (P < 0.001, compared to the controls), correspond to a 12% reduction in the axon bundle area (81.9 ± 17.8 µm², n = 8, P = 0.5). At 15 days pONC, the axon bundle area decreased by 37% (58.4 ± 10.4 µm², n = 3; P = 0.03), whereas RGC density suffered an 85% reduction (P < 0.001, Fig. 6A). 

To estimate where the correlation switches from negative to positive for each parameter, we used the second-order polynomial regression model to fit the bundle width plot (R² = 0.55, P < 0.001), height plot (R² = 0.45, P < 0.001), and area plot (R² = 0.55, P < 0.001) (dashed lines in Figs. 6A–C; also see Table 1). The second-order polynomial regression model is a descriptive model that it provides a snapshot of what happens to the bundles as soma density varies. We found that the axon bundle area peaked at 152 µm² with an RGC density of 3780 cells/mm² (8% cell loss), whereas the width and height both peaked at a lower RGC density of 3000 cells/mm² (27% cell loss) and 2875 cells/mm² (30% cell loss), respectively. Our results reinforced the idea that the area might represent a more valuable parameter to track subtle changes in axon bundle morphology after injury. We did not observe a clear trend of shape changes in axon bundle changed after ONC injury (Fig. 6D), to which we fit a linear regression model (P = 0.47, One-way ANOVA; Fig. 6D). The R² values, significances, and estimates of peak values are listed in Table 1. 

Table 1.

 

Regression Analysis of RGC Axon Bundle Morphology and Soma Loss (See Figure 6)

 

Location-Dependent Changes in Axon Bundle Morphology After ONC

Because of the unpredictability of location and the extent of RGC damage caused by the crush surgery, we compared the relationship between bundle morphology and RGC density in different retinal quadrants after ONC. In Figure 7, the four axon bundle parameters for the superior (S), inferior (I), nasal (N), and temporal (T) regions of the retina were plotted as a function of RGC density for each orientation. Second-order polynomial regression models were fit to the data for each region (dashed lines, Figs. 7A–C). The R² values, significances, and estimates of peak values for all four regions are listed in Table 2 (all regressions were statistically significant, except for temporal width regression). For example, at three days pONC, the average density of RGCs had significantly decreased by 20% in the superior leaflet (3314 ± 356 cells/mm² for control, and 2681 ± 337 cells/mm² at three days pONC, P < 0.01) and significantly decreased by 13% in the nasal leaflet (4371 ± 334 cells/ mm² for control, and 3826 ± 385 cells/mm² at three days pONC, P < 0.01). The 20% RGC reduction in the superior region was correlated with a 48% significant increase of axon bundle area in the same region (control: 94 ± 14.9 µm², three days pONC: 140 ± 26.5 µm², P = 0.01). The 13% reduction in the nasal RGC density correlated with a 36% increase in the RGC axon bundle area (control: 91 ± 11.7 µm², three days pONC: 124 ± 31.6 µm²), although the change was not significant (P = 0.07, Fig. 7A). At 15 days pONC, axon bundles from both nasal (15 days pONC: 56 ± 15.7 µm², P = 0.2) and superior regions (15 days pONC: 54 ± 14.8 µm², P = 0.1) showed a reduction in area, which correlated with an 84% (nasal: 15 days pONC: 712 ± 52cells/ mm², P < 0.001, Student's t test) and 87% reduction in RGC density, respectively (superior: 15 days pONC: 414 ± 20 cells/ mm², P < 0.001, Student's t test, Fig. 7A). 

Figure 7.

 

Regional differences were detected between RGC density and axon bundle size measurements. (A–D) RGC axon bundle cross-sectional area (A), width (B), height (C), and shape (D) plotted as a function of rbpms + RGCs density from each quadrant of each retina (S, superior; I, inferior; T, temporal; N, nasal). Second-order polynomial regressions (dashed lines) were fitted to each region separately for A–C, and linear regressions were fitted for D.

Table 2.

 

Region-Dependent Polynomial Regression Analysis of RGC Axon Bundle Morphology and Soma Loss (See Figure 7)

 

A linear regression model was used to fit the shape parameter plot for each region (Fig. 7D). Although no regression was significant for superior, temporal, and nasal regions, a significant linear relationship was found between axon bundle shape and RGC density in the inferior region (R² = 0.25, P < 0.01), suggesting uneven changes in axon height and width after the ONC. 

Numerical Simulation Suggests Bundle Area is Most Sensitive to RGC Damage

We performed a numerical simulation to determine which bundle parameter is most sensitive to detect RGC damage. To model the change in RGC soma density as a function of days after ONC, we fit our experimental soma density data with a logistic function, as shown in Figure 8A. Next, we simulated soma density values normally distributed along the modeled function within the recorded standard deviation. A total of 10 density values were simulated every 0.25 days, as shown in Figure 8B. Simulated density values were used as input values for the bundle parameter models described in Figure 6. Simulated parameter values were normally distributed along the model for each parameter using the experimentally recorded standard deviation values. The simulated parameter values were plotted as a function of time and an unpaired t test (P > 0.05) was used to determine at which time the parameter become significantly different from the baseline. Figures 8C–F show the change over time for simulated cross-sectional area (C), width (D), height (E), and shape (F). We repeated the simulation 50 times and the recorded average P value as a function of time for each parameter. Figure 8G shows the simulation results for determining the most sensitive parameter (black dashed line indicates significance threshold): RGC density reaches the significance threshold at 1.3 days (blue), followed by bundle area at 9.6 days (purple), width at 10.1 days (orange), and height at 10.4 days (yellow). Unlike the other bundle parameters, bundle shape (green) never crosses the significance threshold, indicating that it does not have high enough sensitivity to detect a significant difference within 25 days of the ONC procedure. 

Figure 8.

 

Simulation using experimental data to determine which RGC axon bundle size parameter is most sensitive to RGC damage and estimate parameter floor values. (A) RGC density (ρ_RGC, blue data points) plotted as a function of time to establish relationship modeled by a logistic function (dashed red line). (B) Example simulated RGC density values (black data points) as a function of time. (C–F) Example simulated RGC axon bundle size parameters: cross-sectional area (A_B) (C), width (W_B) (D), height (H_B) (E), and shape (S_B) (F). (G) P values as a function of time for the simulated size parameters to determine which parameters are most sensitive to RGC loss. (H) P values as a function of time for the simulated size parameters to determine parameter floor times (t_f).

The floor value is defined as the time point beyond which no further change in the respective parameter is detected.^37,38 Using our simulated data, we estimated the floor value for each parameter. An unpaired t-test with a significance level of 0.05 was used to compare each simulated time point with the final time point. The time at which the p-value of the parameter crossed above the significance threshold was defined as floor time (t_f). The mean of the parameter value beyond t_f was treated as the floor value. The simulation result for determining t_f for each parameter is shown in Figure 8H: bundle width is the first to cross the significance threshold at 12.9 days (orange), followed by height at 13.4 days (yellow), area at 13.9 days (purple), and RGC density at 15.4 days (blue). The bundle shape parameter (green) did not cross the significance threshold and, thus, does not have a parameter floor value within 25 days of the ONC procedure. RGC density had a floor value of 131 cells/mm², area had a floor value of 35.0 µm², width had a floor value of 6.4 µm, and height had a floor value of 7.1 µm. The parameter floor model suggests that bundle area remains sensitive to size changes for the longest amount of time after crush because it reaches its t_f later than bundle width or height. 

Discussion

RGC Axon Bundle Morphology as a New In Vivo Biomarker for RGC Damage

Glaucoma can progress without having easily identifiable symptoms until reaching an advanced stage of vision loss. Early diagnosis and intervention are crucial to slow glaucoma progression.⁴⁶ Clinical OCT systems for the diagnosis and monitoring of optic neuropathies operate using NIR illumination. By shifting the illumination wavelengths to the visible light spectrum (510 nm to 610 nm), vis-OCT has an improved axial resolution of 1.3 µm in the retina compared to 4 µm in the best clinical NIR OCT devices.³³ In addition, vis-OCT has greater contrast between retinal layers due to the higher backscattering properties of biological tissues in the visible light wavelength range.³³ Taking advantage of this improved resolution and increased contrast sensitivity, we developed vis-OCTF to analyze individual RGC axon bundles in vivo.^34–36,47 Thus far, we have demonstrated vis-OCTF in the mouse retina to visualize RGC axon bundles in vivo and validated these structures using confocal microscopy ex vivo.³⁴ We then applied vis-OCTF to visualize changes in RGC axon bundle structure in the case of increased RGC population using BAX^−/− mice and decreased RGC population using ONC mice.³⁵ In the present study, we developed new analytic tools for extracting RGC axon bundle size parameters from vis-OCTF images to determine which parameter is most sensitive to RGC damage. 

Because every RGC extends one axon in the RNFL, we seek to examine whether directly quantifying changes of individual RGC axon bundles can be a more sensitive and accurate indicator for RGC damage than bulk RNFL or GCIPL thickness. In this study, we used an acute ONC injury model to longitudinally track morphological changes in single RGC axon bundles. We measured four bundle size parameters: (1) lateral width, (2) axial height, (3) cross-sectional area, and (4) shape. First, we found bundle width from vis-OCTF is more sensitive to early damage compared to bundle height. The reduction of lateral width was detected between three and six days pONC (Figs. 3A-3E), which was earlier compared to the reduction in the bundle height (six to nine days pONC, Figs. 3B–E). 

Second, the bundle height measurements obtained after ONC matched the fact that RNFL thinning is not sensitive to the early neuropathic damage.^22,23 We first observed a significant thinning of the RGC axon bundle at 9 days pONC (Figs. 3B–E). Importantly, we detected morphological changes in the RGC axon bundles but no significant change in the overall height of the GCIPL at three days pONC (Fig. 4), supporting the notion that RGC axon bundle morphology could serve as a more sensitive indicator than the GCIPL thickness. 

Finally, we introduced two novel parameters to measure and track the size and shape of single RGC axon bundles. By combining changes in both dimensions (lateral and axial), the cross-sectional area has shown to be an accurate indicator of RGC damage. At three days pONC, we found that 60% of all axon bundles showed a cross-sectional area increase of 30%, corresponding to about a 14% increase in width and a 15% increase in height (Figs. 3E, 3F). Swelling of the RNFL⁴⁸ and retina³⁰ was also observed in early response to optic nerve injuries. However, in vivo changes in individual RGC axon bundles have not been demonstrated previously. Combining vis-OCTF with the new analytic tools we developed, we detected and quantified in vivo, for the first time, the swelling of individual axon bundles after acute ONC injury. 

RGC and Axon Degeneration

Dendritic and axonal dysfunction is an early event in animal glaucoma models and may precede RGC soma degeneration.^46,49 The ONH is hypothesized to be one of the most vulnerable structures to disease insult by glaucoma.^46,50,51 In our study, we observed obvious axon bundle swelling at three days pONC (Figs. 2, 3, 5, 6), an early time point when only 15% of RGC somas have degenerated. This suggests that morphological changes in axon bundles occur early in disease progression, which agrees with previous studies.⁵⁰ In both DBA-2J mice, a genetic model of glaucoma⁵¹ and a rat ocular hypertension model,⁵⁰ axonal cytoskeleton damage was observed at early disease stages. Although individual axon damage could lead to RNFL swelling, other mechanisms may also contribute, including inflammatory responses such as macro and microglial proliferation in the RNFL.^52–54 

Interestingly, although axons were affected early in the disease, the rate of axon shrinkage was not as fast as that of RGC soma loss. At nine days pONC, we observed 68% of RGC soma loss whereas the axon bundles showed about a 12% decrease in cross-sectional area. At 15 days pONC, most RGC somas had degenerated (85%), whereas a substantial amount of axon bundles remained (30% decrease in axon bundle area). This desynchronization of soma and axon loss suggested a compartmentalized degeneration in RGCs, which also agrees with previous studies.^55,56 Future studies are needed to investigate the molecular and cellular mechanisms underlying RGC loss, axon bundle degeneration, as well as the changes in RNFL. 

We also show regional differences in bundle changes and RGC soma degeneration (Fig. 6). In wildtype retinas, the distribution of general RGCs varies across different regions.^31,57,58 For example, the density of rbpms-positive RGCs in the superior region of the retina is 24% lower than the nasal retina in control mice. This difference in the overall RGC density could affect axon bundle organization and change how they react to disease insults. Furthermore, RGC types distribute unevenly across the retina,^31,57,58 and different types of RGC respond differently to injury and disease.^40,58–60 For example, one type of RGCs, the ipRGC, is involved in circadian functions,⁶¹ and at least one subtype, the M4 ipRGCs, has higher distribution in the superior and temporal retina.⁵⁷ It has been found that ipRGC's are more resistant to chronic elevation of IOP⁴³ and ONC injury⁶² than general RGCs. The retained circadian functions of the mice after chronic IOP elevation also suggest a normal function of ipRGCs. Therefore it is likely that ipRGC axons suffered less damage and morphological changes than other RGCs in the case of ONC injury, potentially contributing to the regional difference in axon bundle morphology changes observed in our study. More studies remain needed to investigate whether the axons of RGC are organized into axon bundles based on function or location. 

In summary, we established vis-OCTF parameters to track axon bundle morphology in vivo after the acute ONC injury. Our experimental and simulated results concluded that RGC axon bundle cross-sectional area is most sensitive to RGC damage. Our current study presented a new possibility to establish an in vivo quantifiable biomarker for RGC degeneration with glaucoma development and progression in future. 

Acknowledgments

The authors thank William Tucker for his technical support on bundle quantification and Ignacio Provencio for his insightful discussions. 

Supported in part by NIH grants R01EY029121, R01EY019949, U01EY033001, and R44EY026466, and Glaucoma Research Foundation. 

Disclosure: M. Grannonico, None; D.A. Miller, None; J. Gao, None; K.M. McHaney, None; M. Liu, None; M.A. Krause, None; P.A. Netland, None; H.F. Zhang, Opticent Health (F); X. Liu, None