QuPath Automated Analysis of Optic Nerve Degeneration in Brown Norway Rats

Abstract

Purpose: A novel application of QuPath open-source digital analysis software is used to provide in-depth morphological analysis of progressive optic nerve (ON) degeneration in rats.

Methods: QuPath software was adapted to assess axon and gliotic morphology in toluidine blue-stained, Brown Norway rat ON light micrographs. QuPath axon numbers, density, size distributions, and gliotic areas were obtained from test images and ON cross-sections separated by damage grade. QuPath results were compared with manual counting, AxonJ, and electron microscopy axon estimates.

Results: QuPath-derived axon number, density, and diameter decreased with increasing ON damage. Axon density negatively correlated with gliotic areas in test images (R² = 0.759; P < 0.0001; N = 40) and in ON cross-sections (R² = 0.803; P < 0.0004; N = 10). Although axon losses occurred across most axon diameters, large axons were more susceptible to degeneration. The exception was swollen axons > 2 µm, which increased in moderately but not severely damaged images. QuPath axon counts correlated strongly with manual counts of test images (R² = 0.956; P < 0.0001). QuPath outperformed AxonJ on test images and total ON axon counts. Compared to electron microscopy analysis, QuPath undercounted ON axons; however, correlation between the methods was robust (R² = 0.797; P < 0.001; N = 10).

Conclusions: QuPath analysis reliably identified axon loss, axon morphology changes, and gliotic expansion that occurred in degenerating ONs.

Translational Relevance: QuPath is a valuable tool for rapid, automated, analysis of healthy and degenerating ONs. Reproducible preclinical studies for new glaucoma treatments depend on unbiased in-depth analysis of ON pathology. This was provided by the QuPath approach.

Introduction

Glaucoma is the leading cause of irreversible blindness in the world.¹ Glaucoma pathology is characterized by optic nerve (ON) degeneration, including loss of retinal ganglion cell axons² and gliosis.^3–6 Signs of gliosis include activation of glial cells, glial hyperplasia and hypertrophy, astrocyte-mediated extracellular matrix remodeling, and changes in astrocyte gene expression and morphology.^4,5,7–11 Although qualitative grading schemes describing ON degeneration include gliosis in the damage criteria,^12,13 automated quantification of gliotic structures is rare. Accurate assessments of both ON axon loss and gliosis are required to improve rigor and reproducibility of preclinical animal trials for glaucoma treatments. 

Despite valuable contributions from numerous vision scientists, quantifying ON damage still poses significant challenges. The gold standard for postmortem quantification of axon numbers within the ON is axon counting from electron microscope (EM) images.¹⁴ Although manual axon counting across entire ON cross-sections is labor intensive and requires numerous EM images, electron microscopy easily detects small axons that are difficult to resolve by light microscope. Sampling schemes covering 3% to 20% of the total ON area reduce the axon counting effort^15–19; however, variability in total axon counts increases dramatically when sampled areas are less than 8% of the total ON area.²⁰ Given the regional variability of axon degeneration within the ON,^21–23 axon counts from sampled images may not accurately reflect the health of the entire ON. 

As computer technology advances, semi-automated and fully automated approaches to axon quantification are increasingly accessible to all investigators. Initial studies using computer algorithms for semi-automated axon analysis of sampled images often required proprietary software or in-house applications.^{19,21,24–28} Few investigators have attempted fully automated counts of all axons across the entire ON.^19,29 Development of ImageJ and Fiji open-source digital image analysis software^30,31 has increased access to software tools required for automated axon counting.^32–34 A major breakthrough was development of AxonJ, a publicly available ImageJ plug-in that quantifies rodent axons across entire ONs.²⁹ AxonJ is an important tool for axon counting; however, the AxonJ program does not assess gliosis. 

Here, we investigated the novel application of QuPath open-source image analysis software³⁵ for characterization of ON degeneration in Brown Norway rat retired breeders. QuPath software was originally designed to analyze tumor tissue pathology; however, the flexible design and general image analysis capabilities of QuPath allow for its adaptation to many other applications. In this study, specific QuPath parameters were selected for standard QuPath cell counting algorithms enabling recognition of both axons and gliotic elements. QuPath analysis was combined with the separation of images by damage grade to evaluate changing axon morphology and gliotic areas with progressive ON degeneration. QuPath results were compared to existing counting methods. 

Methods

Animals

A total of 10 ONs came from six Brown Norway rats, 6.5 to 13 months of age (Charles River Laboratories, Wilmington, MA) (Supplementary Table S1). Rats were maintained at an Augusta University animal facility in accordance with the Public Health Service Policy for Humane Care and Use of Animals and the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Vision Research. All work was performed under animal protocols approved by the Augusta University Institutional Animal Care and Use Committee. 

Microbead Injections

To induce ON degeneration in a subset of ONs, intraocular pressure (IOP) was raised by injecting polystyrene microbeads (Polysciences, Inc., Warrington, PA) into the anterior chamber.³⁶ Briefly, microbeads (10 µm and 6 µm; 2.7% solid suspension) were sterilized in 100% ethanol, washed three times in sterile phosphate-buffered saline (PBS), and resuspended in sterile PBS at one to four times the manufacturer's supplied concentration, resulting in a final concentration 2.7% to 10.8% by weight. Rats were anesthetized with a cocktail of 60 mg/kg ketamine (Zoetis, Kalamazoo, MI) and 7.5 mg/kg xylazine (Akorn Pharmaceuticals, Lake Forest, IL). The eye was numbed with one drop of proparacaine HCl, 0.5% (Akorn Pharmaceuticals). A 100-µl SGE Gas Tight Syringe, Luer Lock 100F-LL-GT (World Precision Instruments, Sarasota, FL), was coupled to a pre-pulled 30-µm inner diameter glass micropipette tip (#TIP30TW1, World Precision Instruments). The syringe was loaded from behind with 10 µl of 6-µm beads, 10 µl of 10-µm beads, and 10 to 20 µl of HEALON5 viscoelastic (Johnson & Johnson Vision Care, Inc., Jacksonville, FL). The beads and viscoelastic were rapidly discharged over a period of approximately 5 seconds into the anterior chamber. Acute IOP elevation (50–60 mm Hg) sometimes accompanied bead injections but dropped to 20 to 40 mmHg within 24 hours. Contralateral eyes were not injected. Naïve eyes were from rats receiving no injection in either eye. IOP was monitored two or three times per week with a handheld iCare tonometer (iCare Finland Oy, Vantaa, Finland) in awake, non-anesthetized rats. Average IOPs for bead and contralateral eyes were calculated as area under the IOP versus time curve divided by total time of experiment. 

Removal of ONs, Tissue Processing, EM Analysis

Five weeks after the microbead injection, the rats were euthanized by CO₂. A small piece of ON (1–2 mm in length) was removed approximately 1.5 mm posterior to the globe and placed in fixative (4% paraformaldehyde, 2% glutaraldehyde in 0.1-M sodium cacodylate buffer; pH 7.4); it was then postfixed in 2% osmium tetroxide, stained en bloc with uranyl acetate, dehydrated with a graded ethanol series, and embedded in Epon–Araldite resin. Semithin (0.5–1 µm) sections were stained with 1% toluidine blue for subsequent analysis by light microscopy. Reagents were obtained from Electron Microscopy Sciences (Hatfield, PA). 

For EM analysis, thin sections were cut, collected on copper grids, and stained with uranyl acetate and lead citrate. ON tissue was observed in a JEM 1230 transmission electron microscope (JEOL USA, Inc., Peabody, MA) at 110 kV and imaged with an UltraScan 4000 charge-coupled device camera (Gatan, Inc., Pleasanton, CA) and First Light Digital Camera Controller (First Light Imaging Corp., Cambridge, MA). EM images (2040 × 2040 pixels2 covering 1752 µm²) were taken at 1000×. EM images were taken randomly across the entire ON cross-section, with the sum of EM image areas covering at least 10% of the ON area. Average axon density was computed from manual counts of axons with intact myelin sheaths. Average axon density was multiplied by ON area to compute total ON axon number. 

Light Microscope Imaging

Of the 10 to 15 toluidine blue-stained ON cross-sections on each slide, one was chosen for light microscope imaging based on the presence of even staining and absence of artifacts. Light microscope imaging was performed using a Zeiss Axio Imager D2 microscope equipped with a high-resolution camera and Zeiss ZEN 2.3 imaging software (Carl Zeiss AG, Oberkochen, Germany). Overlapping 63× oil magnification images covering the entire ON were aligned and stitched together using Adobe Photoshop Elements 2.0 (Adobe, Inc., San Jose, CA) to form high-resolution composite ON cross-sections saved in the .jpeg format. ON composite images are available in Supplementary Data S1 (key) and Supplementary Data S2- to S11 (images). Forty test images, each 438 × 438 pixels covering 1000 µm², were selected from composite ON cross-sections (stitched montages). No particular image quality metric was used except to select test images from a range of healthy and mildly, moderately, or severely damaged ON regions. A listing of test sample grade and source is included in Supplementary Data Table S1. 

Image Grading and Manual Axon Counts

Test images were graded by three independent investigators on a scale of 1 to 5 according to Jia et al.¹³ and then further classified as healthy and mildly, moderately, or severely damaged. Healthy ON images (grade 1) had few degenerating axons. Mild damage (grade 2) had focal degeneration of axons with some swollen axons. Moderate damage (grades 3 and 4) had numerous degenerating axons and axon swelling with some normal axons. Severe damage (grade 5) had gliosis with degenerating axons across the entire ON. Manual counts of axons with intact myelin sheaths were averaged from three independent investigators. 

QuPath Method

QuPath software³⁵ is available for Windows, Mac, and Linux platforms. QuPath (v0.1.2) for Windows was downloaded from https://qupath.github.io/. Minimum Windows requirements were a 64-bit processor and operating system, 4 GB RAM, and 500 MB hard disk space. This study used a Dell computer with Intel Core i5-2400 CPU with 3.10-GHz processor, 8 GB RAM, and ample hard disk space. The QuPath online manual is available at https://github.com/qupath/qupath/wiki. 

Basic steps for QuPath analysis are summarized in the flowchart shown in Fig. 1A and in Supplementary Data S12 (QuPath instructions). First, axons (green) and glia (purple) were added to the annotations menu and image type set to fluorescence. For test images, the entire image was selected for analysis. For ON cross-sections, the ON area excluding outer connective tissue was selected and subsequently subdivided into 200 × 200 pixel2 tiles to facilitate axon identification. Watershed cell detection routine (Fig. 1B) was used to identify objects. The choice of watershed cell detection parameters and optimization of threshold value to 10 reduced the effects of image exposure on axon counts. QuPath returned 12 different parameters for each detected object (Fig. 1C). Circularity, defined as 4πA/P² where A is area and P is perimeter, was used to classify objects as axons or gliotic elements.^19,34 Rat axons had circularity ranging from 0.5 to 1.0 and were colorized green.³⁴ Gliotic objects had circularity < 0.5 and were colorized purple. Size data were saved in the Excel spreadsheet format. 

Figure 1.

 

QuPath method and parameters. (A) QuPath flow chart. (B) Watershed cell detection parameters. (C) QuPath measures 12 different parameters for each object; area and circularity were used to investigate axon morphology and gliotic areas.

Percent Gliotic Areas and Axon Distributions

QuPath-generated percent gliotic area and axon diameter distributions were calculated from QuPath size data. Percent gliotic area was 100*(Sum of gliotic area)/(Total area). Axon diameter (D) was calculated assuming a circular axon area (Area = πD²/4). The Excel histogram function computed axon diameter frequencies. Relative axon frequency was calculated for each bin range as 100*(Axon frequency)/(Total axon number). Bin size was 0.25 µm and ranged from 0 to 2 µm, plus there was an additional category for axons with diameter > 2 µm. For manual axon diameter distributions, axon areas excluding the myelin sheath were outlined for all axons identified by manual counting of the healthy test images. Axon areas measured by digital analysis software were recorded in an Excel spreadsheet and axon diameter distributions calculated as described for the QuPath axon size data. 

AxonJ Software

AxonJ²⁹ is a freely available plug-in for counting rodent ON axons in either ImageJ³⁰ or Fiji³¹ open-source image analysis software. Fiji for Windows 64-bit (https://imagej.net/Fiji/Downloads) was used in this work. AxonJ (https://imagej.nih.gov/ij/plugins/axonj/index.html) and CLAHE (https://imagej.nih.gov/ij/plugins/clahe/index.html) plug-ins were downloaded and installed using the Plugins>Install plug-ins command. FeatureJ_.jar (version 2.0.0) and imagescience.jar (version 3.0.0) were downloaded from https://imagescience.org/meijering/software/featurej/ to the Fiji plug-ins folder. The command prompt Plugins>AxonJ was used to run AxonJ. Results for axon numbers and sizes were saved in an Excel spreadsheet format. 

Statistics

Statistics were performed with GraphPad Prism 8 (GraphPad Software Inc., La Jolla, CA). For multiple comparisons with normal distribution one-way analysis of variance (ANOVA) was used with Welch's correction for unequal variances. For comparison of test image manual counts to QuPath and AxonJ, ANOVA with repeated measures was used. Post hoc testing was by the Tukey–Kramer method. Statistical analysis of manual axon diameter frequency distributions was by two-way ANOVA of manual data normalized to QuPath data within each bin range followed by Sidak's multiple comparison test (Fig. 4C). Analysis of test image and ON axon frequency data was by two-way ANOVA of axon frequency data normalized to healthy test image or healthy ON data within each bin range followed by Dunnett's multiple comparisons with healthy controls (Figs. 4E, 4H, 7B, and 7E). Significance was set at P < 0.05. 

Results

A novel application of QuPath³⁵ open-source digital image analysis software was used to assess axon morphology and gliotic areas in 10 ON cross-sections and a subset of 40 test images. The six Brown Norway rats used in this study were retired breeders, male and female, 6.5 to 13 months of age. Retired Brown Norway breeders are an aged rat model consistent with increased risk of glaucoma with age in humans. Charles River did not supply retired breeders at a specific age, resulting in a wide variance in rat age. The 10 ONs in this study came from six rats. Four ONs were from microbead-injected eyes (four rats) to induce increased ON degeneration, two ONs were from eyes contralateral to the bead-injected eyes, and four ONs were from two naïve rats not injected in either eye (Supplementary data Table S1). In our hands, the microbead model of glaucoma produced variable elevation of IOP depending on the effectiveness of bead delivery and responses of individual animals. Average IOP over 5 weeks after microbead injection ranged from 10.8 to 30.1 mm Hg in bead-injected eyes and from 13.2 to 18.0 mm Hg in contralateral eyes (Supplementary Table S1). Visual inspection of composite toluidine blue-stained ON cross-sections by 63× oil light microscope showed that ONs ranged from healthy to severely damaged. The 40 test images were chosen from healthy and mildly, moderately, or severely damaged regions of interest derived from the 10 high-resolution ON cross-sections as listed in Supplementary Table S1. Test images were separated by three independent investigators into healthy (n = 8) and mildly (n = 8), moderately (n = 14), or severely (n = 10) damaged categories. 

Test image results showed progressive changes in ON morphology with damage severity. Axon losses were accompanied by expansion of gliotic areas with increasing damage grade. Figure 2A shows representative 63× oil light micrographs of healthy and mildly, moderately, and severely damaged ON test images derived from the semithin (0.5–1 µm) ON cross-sections stained with 1% toluidine blue. Objects detected by QuPath are outlined in red (Fig. 2B). Objects were then classified based on circularity as axons and colorized green or as gliotic areas and colorized purple by the QuPath software (Fig. 2C). QuPath quantification showed that axon numbers decreased significantly in moderately and severely damaged ON test images compared to healthy and mildly damaged images and between moderately and severely damaged ON test images (Fig. 3A). Mean axon diameter decreased significantly between mildly damaged and healthy test images and in the severely damaged test ON images compared to all other image categories (Fig. 3B). 

Figure 2.

 

QuPath identification of axons and gliotic elements in representative test images. Test images were classified by damage grade, illustrating QuPath identification of axons and gliotic elements. (A) Representative healthy and mildly, moderately, and severely damaged light microscope images (63× oil magnification). (B) Corresponding images with QuPath detected objects outlined in red. (C) Corresponding images showing QuPath identification of axons (green) and glial elements (purple). Scale bar for all images is 5 µm.

Figure 3.

 

Quantification of axons and gliotic areas in test images. QuPath results showed progressive changes in axon morphology and gliotic areas with damage severity. (A) Axon number per test image and (B) mean axon diameter decreased significantly with damage severity. (C) Percent gliotic area increased significantly with damage severity. (D) Linear regression of axon density (axon number/1000 µm² image area) with percent gliotic area produced a line with negative slope –7.91 ± 0.72 and goodness of fit R² = 0.76 (P < 0.0001; N = 40). ^**P < 0.01, ^***P < 0.001, ^****P < 0.0001 compared to healthy images; ^##P < 0.01, ^###P < 0.001, ^####P < 0.0001 compared to mildly damaged ON images; ^P < 0.05, ^^P < 0.01, ^^^^P < 0.0001 between moderately damaged and severely damaged ON images; n = 8 to 14 in each image category, for a total of 40 test images.

In contrast to axon numbers, area occupied by gliotic objects increased with damage severity. Percent gliotic area increased significantly in mildly, moderately, and severely damaged image categories compared to the healthy images and between severely damaged images and all other image categories (Fig. 3C). Linear regression of axon density (axons/1000 µm²) with percent gliotic area (Fig. 3D) produced a line with negative slope –7.91 ± 0.72 and goodness of fit R² = 0.76 (P < 0.0001; N = 40). 

QuPath automated axon analysis allowed detailed examination of axon size distributions which were validated against manual axon size distributions for the healthy test images. Small, medium, and large axons were observed throughout healthy ON test images (Fig. 4A). For axon size distributions, bin size was 0.25 µm with a bin range of 0 to 2 µm, plus one size category with axons > 2 µm. QuPath relative axon diameter frequencies for healthy test images (Fig. 4B) peaked at 0.75 µm with a value of 35.2 ± 1.3%. The distribution was unimodal with a tail toward large axon sizes. Few axon diameters were less than 0.25 µm or greater than 2 µm. Manually derived axon diameter distributions were similar to those for QuPath (Fig. 4B) except for the lowest and highest bin ranges (Fig. 4C). Comparison of manual axon frequencies to QuPath within each bin range (Fig. 4C) showed that the fold changes in manual axon diameter frequencies were significantly lower in the 0.25-µm bin range and significantly greater in the >2-µm bin ranges compared to QuPath. 

Figure 4.

 

Axon diameter distributions in test images by damage category. (A) Representative 63× oil magnification image of small (0.4 µm), medium (0.75 µm), and large (3 µm) diameter axons (single-headed arrows). Double-headed arrow illustrates inner axon diameter. (B) Graph showing relative axon diameter frequency distributions derived from QuPath and manual axon size data from healthy test images. (C) Fold change of manual axon size distributions compared to QuPath within each bin range showing similar relative axon diameter frequencies except for the 0.25-µm and >2-µm axon diameter bin ranges. Analysis by two-way ANOVA of normalized data followed by Sidak's multiple comparison test (^***P < 0.001, ^****P < 0.0001; n = 8). (D) Graph of QuPath-generated axon diameter distributions showing axon frequencies of healthy and mildly, moderately, and severely damaged ON test image categories by bin range. (E) Graph showing fold change of mildly, moderately, and severely damaged ON axon frequencies relative to healthy test images within each bin range. (F) Chart showing that significant fold-change decreases in axon frequencies occurred among axon diameters 0.75 µm and above in moderately and/or severely damaged ON compared to healthy test images. Exceptions to this trend occurred within the axon diameter bin ranges of 0.25 µm and >2 µm. (G) Graph of QuPath-generated relative axon diameter distributions showing axon frequency percent of healthy and mildly, moderately, and severely damaged ON test images with bin range. (H) Graph showing fold change of mildly, moderately, and severely damaged ON axon frequencies relative to healthy test images within each bin range. (I) Chart showing significant fold-change increases in axon frequency percent for 0.25-µm and 0.5-µm diameter bin ranges and significant fold change decreases among 1.5-, 1.75-, 2.0-, and >2-µm diameter axons in severely damaged compared to healthy test images. Moderately damaged ON test images had significant fold-change increases in 0.25-µm and >2-µm bin ranges compared to healthy test images. Analysis of fold-change axon frequencies and axon frequency percent was by two-way ANOVA of data normalized to healthy test images within bin ranges followed by Dunnett's multiple comparisons with the healthy test images. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001, ^****P < 0.0001; n = 8 to 14 for each test image category; N = 40.

QuPath automated axon analysis was next used to compare axon size distributions in healthy versus damaged ON test images. Axon losses were observed across a range of axon diameter bin ranges in moderately and severely damaged ON test images (Fig. 4D). Normalization of damaged axon diameter frequencies to healthy values within each bin range (Fig. 4E) showed significant fold-change decreases in axon frequencies among axons with diameters of 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, and >2 µm in the moderately and/or severely damaged test ON image categories (Figs. 4E and 4F). Exceptions to this trend occurred within the 0.25-µm and >2-µm bin ranges of the moderately damaged ON test images. which had significant fold-change increases in axon frequencies compared to healthy test images (Figs. 4E and 4F). 

Next, changes in relative axon diameter frequencies were examined. The relative axon diameter frequency distribution is the fraction or percent of axons with diameters within each bin range. Relative axon frequency is calculated by dividing the axon frequency within each bin range by the total number of axons. Axon frequency percent gives important information about relative numbers of axons within each bin range. In Figure 4G, relative axon diameter distributions (axon frequency percent) shifted toward small-diameter axons at the expense of large axons in moderately and severely damaged ON test images. The fold change in axon frequency percent with damage grade compared to healthy test images (Fig. 4H) decreased within the 1.5-, 1.75-, 2.0-, and >2-µm bin ranges for the severely damaged compared to healthy test image category (Fig. 4I). In contrast, increased fold change in the percentage of very small diameter axons (0.25 and 0.5 µm) were observed in the moderately and/or severely damaged groups and in the very large (>2 µm) moderately damaged test ON images compared to healthy test images (Figs. 4H and 4I). 

After analyzing the test images, we examined whether the QuPath approach could be extended to entire ON cross-sections. For analysis purposes, 10 ONs were separated into three groups (healthy, moderately damaged, or severely damaged) based on the percent gliotic area. The healthy ON category contained ONs with gliotic areas < 20%, the maximum gliotic area observed in the mildly damaged category for the test images. None of the 10 ONs analyzed had as low a percent gliotic area across the entire ON as was observed in the healthy test images (see Supplementary Table S1). This finding was understandable, as healthy test images represented only the healthiest regions of interest, not the average health of all regions across the entire ON. For moderately damaged ONs, gliotic areas ranged from 20% to 26% (n = 3). Severely damaged ONs had gliotic areas ranging from 31% to 37% (n = 4). Representative images of healthy, moderately damaged, and severely degenerated ONs (Fig. 5A) illustrate the loss of axons (green) and increased gliotic areas (purple) with damage severity (Fig. 5B). 

Figure 5.

 

QuPath identification of axons and gliotic areas in entire ON cross sections. (A) Representative 63× oil stitched composite images of ON cross-sections from healthy ONs compared to ONs with moderate and severe ON damage. (B) Corresponding colorized images showing QuPath identification of axons (green) and gliotic elements (purple). Scale bar for each ON cross section is 100 µm.

ON cross-sections exhibited axon morphology and gliotic area changes with damage severity similar to those observed in the test images (Fig. 6). Total axon number, axon density (number/1000 µm²), and mean axon diameter significantly decreased in severely damaged ONs (Figs. 6A, 6B, and 6C, respectively). In contrast, percent gliotic area (Fig. 6D) significantly increased in severely damaged ONs compared to healthy and moderately damaged ONs. Linear regression of axon density (Fig. 6E) produced a line with negative slope –8.30 ± 1.30 and goodness of fit R² = 0.803 (P < 0.0004; N = 10), similar to results from test images (Fig. 3D). 

Figure 6.

 

Analysis of axons and gliotic areas in entire ON cross-sections. QuPath results revealed changes in axon morphology and gliotic area with ON damage severity. (A) ON axon number, (B) ON axon density (axon number per 1000 µm²), and (C) mean axon diameter all decreased significantly with ON damage severity. (D) In contrast, ON percent gliotic area increased significantly with damage severity. (E) Linear regression of axon density with percent gliotic area produced a line with slope –8.34 ± 1.30 and a goodness of fit R² = 0.84 (P < 0.0002; N = 10). ^*P < 0.05, ^**P < 0.01, ^***P < 0.001, ^****P < 0.0001 compared to healthy ONs; ^P < 0.05, ^^P < 0.01 between moderately and severely damaged ONs; n = 3 or 4 in each category for a total of 10 ONs.

Subsequent analysis of entire ON cross-sections showed significant differences in axon diameter distributions with increasing damage severity. The axon diameter frequency distribution (Fig. 7A) showed decreased axon frequencies with damage severity. Comparison of moderately and severely damaged ON axon distributions with the healthy ON axon distribution within each bin range (Fig. 7B) showed that, for severely damaged ONs, significant fold-change reductions in axon frequencies occurred for axon diameter bin ranges 0.75 µm and above (Fig. 7C). In contrast, relative axon diameter distributions shifted toward small-diameter axons at the expense of larger axons (Fig. 7D). Comparison of moderately and severely damaged ON axon distributions with the healthy ONs (Fig. 7E) showed that the fold-change percentages of 0.25-µm and 0.5-µm diameter axons increased significantly, whereas the fold-change percentage of axons with diameters of 1.5, 1.75, 2.0, and >2 µm decreased significantly in the severely damaged ONs. In contrast, within the moderately damaged ONs, the fold change in swollen axons > 2 µm increased significantly relative to the healthy ONs (Figs. 7E and 7F). 

Figure 7.

 

ON axon diameter distributions. (A) Comparison of axon diameter distributions for healthy and moderately and severely damaged ONs showed decreased axon diameter frequencies in moderately and severely damaged ONs. (B) Graph showing fold change in axon frequencies of moderately and severely damaged ONs relative to healthy ONs within each bin range. (C) Chart showing significant decreases in fold-change frequencies of axon diameters 0.75 µm and greater in the severely damaged compared to healthy ONs. (D) Graph of relative axon diameter distributions showed a shift toward the percentage of small-diameter axons at the expense of the percentage of large-diameter axons in damaged ONs. (E) Graph showing fold change in axon frequencies of moderately and severely damaged ONs relative to healthy ONs within each bin range. (F) Chart summarizing comparison of axon frequency percent among ON categories. Severely damaged ONs had significant fold-change increases in axon frequency percent for 0.25-µm and 0.5-µm axon diameter bin ranges and significant fold-change decreases among 1.5-, 1.75-, 2.0-, and >2-µm diameter axon bin ranges compared to healthy ONs. Moderately damaged ONs had significant fold-change increases in axon frequency percent of swollen axons > 2 µm compared to healthy ONs. Analysis of fold-change axon frequencies and axon frequency percent was by two-way ANOVA of normalized data to healthy ONs within each bin range followed by Dunnett's multiple comparisons with the healthy test images. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001, ^****P < 0.0001; n = 3 or 4 for each ON category for a total of 10 ONs.

To validate the QuPath analysis of ON degeneration, QuPath results were compared to existing axon counting methods. QuPath correlation with manual axon counts for test images (Fig. 8A) produced R² = 0.956 (P < 0.0001; N = 40). QuPath axon counts were not significantly different than manual counts in the healthy and mildly damaged ON test image categories (Figs. 8C and 8D); however, QuPath significantly overcounted axons compared to manual counts in the moderately damaged and severely damaged ON test image categories (Fig. 8A; boxed points 8E and 8F). 

Figure 8.

 

QuPath comparison with manual counts and AxonJ in test images. (A) Correlation of QuPath versus manual axon counts per test image produced R² = 0.956 (P < 0.0001; N = 40). Comparison to the dotted line, x = y, showed good agreement between QuPath and manual counts except for the severely damage test images (boxed points). (B) Correlation of AxonJ versus manual axon counts produced R² = 0.806 (P < 0.0001; N = 40). Comparison to the dotted line, x = y, showed slightly higher AxonJ counts compared to manual axon counts which was exacerbated in the severely damaged ON test images (boxed points). Breakout of QuPath axon counts per test image compared to manual counting and AxonJ for (C) healthy, (D) mildly damaged, (E) moderately damaged, and (F) severely damaged ON test image categories. ^***P < 0.001, ^****P < 0.0001 compared to manual counts; ^P < 0.05, ^^^^P < 0.0001 between AxonJ and QuPath axon counts; n = 8 to 14 in each category; N = 40; ns, not significant.

QuPath performance was also compared to AxonJ, an alternative open-source automated axon counting program. Correlation of AxonJ counts to manual counts for test images was lower but equally significant (R² = 0.806; P < 0.0001; N = 40) compared to QuPath (Fig. 8B). AxonJ overcounted compared to manual counts in all image damage categories, but this tendency was exacerbated compared to QuPath in the moderately and severely damaged ON test image categories (Fig. 8B; boxed points 8C–8F). 

QuPath and AxonJ automated light microscope counting methods were compared to the gold standard, EM-based axon estimates. EM axon counts were derived from average axon densities determined from random EM images covering at least 10% of the ON area (Fig. 9A). Correlation of QuPath axon counts with EM counts (Fig. 9B) produced R² of 0.797 (P < 0.001; N = 10). Correlation of AxonJ axon counts with EM-based axon counts (Fig. 9B) produced an R² of 0.50 (P < 0.02; N = 10), not as strong as for QuPath. Axon counts from all three methods were significantly different from each other (Fig. 9C). 

Figure 9.

 

QuPath and AxonJ ON axon counts compared to EM axon counts. EM analysis of rat ON axon numbers was performed via manual counting of EM micrographs obtained by random sampling at least 10% of the ON cross-sectional area. (A) Illustration of EM manual counting method and representative EM micrograph with manual counts; scale bar = 5 µm. (B) Correlation of total ON axon numbers calculated by QuPath and AxonJ with EM-based estimates produced R² = 0.797 (P < 0.001) for QuPath and R² = 0.498 (P < 0.03) for AxonJ. (C) Compared to EM axon estimates, QuPath significantly undercounted and AxonJ significantly overcounted total ON axon numbers. (D) Relative axon diameter distribution for one healthy ON. Unlike QuPath, AxonJ failed to account for small-diameter axons < 0.5 µm. ^*P < 0.05 compared to EM counts; ^^^^P < 0.0001 compared to QuPath; N = 10.

Interestingly, QuPath and AxonJ axon diameter distributions were dramatically different (Fig. 9D). Using AxonJ axon size data, we performed analogous calculations as with QuPath to derive AxonJ relative axon diameter distributions for healthy rat ON. AxonJ showed a peak axon diameter frequency of 50.4% at 0.75 µm, which was higher than for QuPath. AxonJ did not detect axons less than 0.5 µm in any of the ON cross-sections. 

Discussion

Application of QuPath software provided rapid automated analysis of ON pathology that quantified both axon and gliotic morphology changes during Brown Norway rat ON degeneration. QuPath takes approximately 4 to 5 minutes to analyze a high-resolution ON cross-section for both axon numbers and gliotic areas compared to 4 to 5 minutes for AxonJ axon quantification or 2 to 3 days of manual axon counting for a single ON. QuPath produced detailed quantification of axon losses and size data that allowed calculation of QuPath-derived axon diameter distributions. Axon distribution data supported a paradigm where axons swell before the preferential loss of large axons. QuPath estimates of gliotic areas confirmed that, in the degenerating ONs, gliotic areas expand to replace lost or shrinking axons. The favorable performance of QuPath compared to existing methods validated the value of QuPath for studying ON degeneration. 

QuPath analysis quantified the axon loss that is a hallmark of ON degeneration in humans and animals.^{3,11,23,24,37,38} Given the regional variability of axon degeneration within the ON,^21–23 healthy and damaged tissues often appear within the same ON. Thus, initial analyses used test images with homogeneous damage levels. Not surprisingly, test images had more significant declines in axon number, axon density, and axon mean diameter and more robust expansion of gliotic areas with damage grade than in entire ON cross-sections. 

This study extended our understanding of axon degeneration by quantifying axon size distribution changes in damaged ONs. QuPath-derived axon diameter distributions in healthy test images were strikingly similar to those obtained by EM analysis of axon morphology in the hooded rat.³⁹ In addition, QuPath axon diameter distributions for the healthy test images corresponded with the manually derived axon diameter distributions, except for the smallest and largest axon diameters, which were over- and underestimated by QuPath, respectively. The small- and large-diameter axons represented very small but important axon diameter fractions and, despite differences with manual analysis, QuPath assessment of axon diameter distributions confirmed that axon shrinkage and swelling were important features of ON degeneration. Although not perfect, QuPath axon diameter distributions provided an automated approach for studying axon morphology changes that would be prohibitively labor intensive to study via manual measurements. In addition, QuPath analysis in whole ONs showed that damaged ONs lost axons across a range of axon diameters (0.75 to >2 µm); however, the relative axon distribution shifted toward small-diameter axons (0.25 µm and 0.5 µm), suggesting that large-diameter axons were more prone to degeneration. This QuPath-generated observation was in accord with studies in primate models showing that large axons are particularly susceptible to glaucomatous damage.^24,38 

One exception to large axon vulnerability was a significant rise in the percentage of swollen axons > 2 µm. This increase occurred in moderately damaged ON test images and whole ON cross-sections but not in severely damaged ON test images or severely damaged ON cross-sections. Axon swelling is part of qualitative grading scales for ON degeneration and occurs in the DBA/2J mouse model of glaucoma before axon loss.^12,13,40 Our results also suggest that ON degeneration is a dynamic process where initial axon swelling may be followed by axon loss. 

This study quantitatively confirmed the qualitative observation that gliotic areas expand to replace lost or shrinking axons in damaged ONs. Here, the term gliotic object refers to all non-axonal objects excluding the myelin sheath. In damaged ON, degenerating axons are replaced by reactive glial cells, invading cells, and extracellular matrix to form glial scars.^4,5,7–11 Admittedly, QuPath identification of objects as gliotic does not distinguish between different elements contributing to gliosis or the vasculature, nor does QuPath account for darker staining striations within gliotic areas. Detecting specific contributions of glial cells to these regions would require specific labeling.⁶ Despite these shortcomings, QuPath identified a significant negative correlation between axon density and percent gliotic areas that prevailed throughout the test images and ON cross-sections. This finding was similar to the DBA/2J model of glaucoma where axonal loss is proportional to glial expansion.³ 

QuPath's value for ON analysis is evidenced by favorable comparisons with existing axon counting methods. QuPath axon counts differed from manual counts only when damage grades reached moderate and severe levels. Whether these differences reflected deficits in QuPath counting is questionable. As damage increased, axons became overly myelinated, making acquisition of high-quality images and identification of viable axons difficult. The QuPath automated method offered the advantage of applying consistent object-recognition criteria to reduce bias from individual investigators. Unfortunately, automated light microscope methods are also limited by their dependence on the availability of high-quality images. One caveat to the present investigation was the relatively small number of ONs and ON cross-sections analyzed. Examining ON cross-sections from ON regions more proximal or more distal to the globe could provide additional interpretations of the data in this paper. Automated methods such as QuPath make analysis of multiple ON regions feasible, something that was previously not practical with manual counting methods. 

Although QuPath undercounted compared to EM total axon estimates, differences between the two methods were consistent and understandable. EM easily detects small axons which are difficult to resolve by light microscope. Thus, QuPath axon counts were expected to be 20% to 30% lower than EM axon counts.^15,32 Despite undercounting, QuPath axon counts still showed significant and robust correlation with EM counts. 

QuPath outperformed AxonJ, an alternative open-source automated axon counting program. The tendency to overcount was exacerbated in AxonJ compared to QuPath in test images. AxonJ also overcounted compared to EM-based axon estimates, a finding that was unexpected, as light microscope-based methods typically undercount compared to EM approaches.^15,32 These differences may be because AxonJ was developed and tested on paraphenylenediamine-stained mouse ONs, not toluidine blue-stained rat ONs. AxonJ axon size distributions were also different than QuPath. AxonJ peak diameter frequency occurred at 0.75 µm, but the magnitude was outside the range reported by Reese et al.³⁹ In addition, AxonJ did not detect axons with diameters less than 0.5 µm, possibly due to software design that rejects axon candidates smaller than 0.144 µm.²⁹ 

Conclusions

In summary, QuPath analysis of ON morphology in Brown Norway rats identified differences in axons and gliotic elements within selected regions of interest (test images) and across entire ONs with variable degrees of ON degeneration. QuPath offered distinct advantages over existing publicly available ON analysis tools, particularly QuPath's ability to quantify gliotic expansion. QuPath's efficient automated analysis of ON morphology is needed to reproducibly evaluate therapeutic efficacy in preclinical glaucoma treatment trials. Whether the QuPath approach can be extended to other types of images or animal models is an intriguing line of future work. 

Acknowledgments

We thank the Augusta University Histology core for their help, including Libby Perry for processing of ONs, toluidine blue staining and sectioning, and preparation of samples for EM; Brendan Marshall for taking EM images; and Penny Roon and Donna Kuminski for additional help staining and processing samples. 

This project was supported by National Institutes of Health Award No. R01EY027406. 

Disclosure: B.A. Mysona, None; S. Segar, None; C. Hernandez, None; C. Kim, None; J. Zhao, None; D. Mysona, None; K.E. Bollinger, None