Single-Cell Profiles of Retinal Ganglion Cells Differing in Resilience to Injury Reveal Neuroprotective Genes

Abstract

Neuronal types in the central nervous system differ dramatically in their resilience to injury or other insults. Here we studied the selective resilience of mouse retinal ganglion cells (RGCs) following optic nerve crush (ONC), which severs their axons and leads to death of ∼80% of RGCs within 2 weeks. To identify expression programs associated with differential resilience, we first used single-cell RNA-seq (scRNA-seq) to generate a comprehensive molecular atlas of 46 RGC types in adult retina. We then tracked their survival after ONC; characterized transcriptomic, physiological, and morphological changes that preceded degeneration; and identified genes selectively expressed by each type. Finally, using loss- and gain-of-function assays in vivo, we showed that manipulating some of these genes improved neuronal survival and axon regeneration following ONC. This study provides a systematic framework for parsing type-specific responses to injury and demonstrates that differential gene expression can be used to reveal molecular targets for intervention.

Figure 1.. scRNA-seq reveals 45 molecularly distinct RGC types in adult mice

A. RGCs (green) reside within the innermost layer of the retina, the ganglion cell layer (GCL). Their axons bundle together to form the optic nerve. IPL, inner plexiform layer, GCL, ganglion cell layer. B. Dendrites of different RGC types have distinct lamination patterns within sublaminae (S)1–5 of the IPL, which determines their choice of presynaptic partners. Stereotyped morphologies are illustrated here for several RGC subclasses and types. INL, inner nuclear layer. C. t-distributed stochastic neighbor embedding (t-SNE) visualization of the transcriptional heterogeneity of 35,699 adult mouse RGCs. Cells are colored by cluster assignments, determined using graph clustering. Clusters are numbered in order of decreasing frequency. D. Relative frequencies of RGC clusters C1–45 (mean ±SD, n=10 replicates). Clusters that matched to known types or subclasses are labeled. E. Dotplot showing the expression patterns of marker genes (rows) specific to different retinal classes across RGC and non-RGC clusters in the data (columns; see color bars, top and right). The size of each circle is proportional to the percentage of cells expressing the gene, and the color depicts the average normalized transcript count in expressing cells. GABA-AC and Gly-AC, GABAergic and glycinergic amacrine cells; HC, horizontal cells; BC, bipolar cells; PR, photoreceptors; Endo, endothelial cells. F. Dotplot showing gene combinations (rows) that uniquely mark RGC clusters (columns). Representation as in panel E for single genes, here normalized to 1. 2- or 3- marker codes, always involve the presence of a marker A, and the presence (e.g. A+B+ or A+B+C+) or absence (e.g. A+B-, or A+B-C+) of markers B and C. In such cases, the size of the circle indicates the percentage of cells satisfying the expression pattern, and the color depicts the average transcript count of positive markers in the cells, normalized to 1 for each combination. G. RGC type frequencies are highly similar between ActinomycinD (ActD)-treated (y-axis) and atlas (x-axis) retinas. H. Dotplot showing gene combinations that uniquely define each RGC type in nominal controls (as in F), are preserved in ActD-treated retinas. Row and column order as in F. I. Scatter plot showing tight correspondence (R_Pearson = 0.93) between relative frequencies of RGC groups found by scRNA-seq (y-axis) versus IHC (x-axis).

Figure 2.. Correspondence of scRNA-seq clusters to RGC types

A. Characterization of novel RGC types by combining FISH (magenta) and IHC on sparsely labeled RGCs in the YFP-H line (green). Examples of S2/S4 laminating C10 and C24 RGCs expressing Gpr88 (left) and Fam19a4 (middle), respectively, and an S5 laminating C25 RGC expressing Slc17a7 (right). IPL sublaminae are drawn based on CALB1 or CALB2 staining (white dashed lines). YFP-H line. Examples of S2/S4 laminating C10 and C24 RGCs expressing Gpr88 (left) and Fam19a4 (middle), respectively, and an S5 laminating C25 RGC expressing Slc17a7 (right). IPL sublaminae are drawn based on CALB1 or CALB2 staining (white dashed lines). Merge in “C25” panel shows labeled cell at higher gain to reveal dendritic morphology. B-F Dotplots highlighting transcriptional distinctions among RGC types within subclasses. Dotted lines separate previously described markers (above) from novel markers identified in this study (below). B: αRGC types. C: T-RGC types. D: F-RGC types. E: ipRGC types. F: S2/S4 laminating RGC types. G. C16 comprising D/V-ooDSGCs can be partitioned into Calb1+ (putative D-ooDSGCs) and Calb1- (putative V-ooDSGCs) cells. H. Consistent with the interpretation in panel G, GFP+ cells in the Hb9 mouse line, which labels V-ooDSGCs, are CALB1- and CALB2+ (magenta). I. Dotplot showing consistent patterns of DE gene expression between W3 types (rows) detected in the droplet-based scRNA-seq atlas (red) and plate-based data from FACS-sorted W3 RGCs (green). Labeled by atlas cluster id. J. Transcriptional relatedness of RGC clusters visualized as a dendrogram reveals subclasses of RGC types (annotation bar, bottom). Dotplot shows expression of key subclass-enriched or -defining genes (rows) in clusters (columns).

Figure 3.. scRNA-seq profiling of RGCs following ONC

A. scRNA-seq was performed on RGCs collected before and at six times following ONC. 8,456–13,619 RGCs were collected at each time point. B. Illustration of a single step of the iGraphBoost procedure to classify RGCs collected at time t_n+1 based on an atlas of RGC types at the previous time point t_n. The procedure is initiated with Atlas RGCs at t₀. In Step 1, gradient boosted trees trained on t_n RGC types are used to classify t_n+1 RGCs. Only high-confidence assignments are applied, and a large number of RGCs remain unclassified at this stage. In Step 2, a Jaccard-weighted k-nearest neighbor graph built on all t_n+1 RGCs is used to propagate labels via nearest-neighbor voting to unassigned RGCs, using the classified RGCs in step 1 as anchors. Successfully classified t_n+1 RGCs are used to classify t_n+2 RGCs in the next iteration. C. Fraction of RGCs that can be confidently assigned to types (y-axis) at each time point following ONC (x-axis). The “one-step” approach (grey) using the atlas RGCs as training data results in a significantly lower proportion of assigned cells among late injured RGCs compared to iGraphBoost (black). D. Dotplot showing that gene combinations uniquely defining each RGC type (row and column order as in Figure 1F) are maintained in 14dpc assigned by iGraphBoost, though reduction in expression level of some markers was observed. E. RGC type-specific resilience at 14dpc relative to control (Ctrl) rank ordered based on decreasing values of the relative frequency ratio at 14dpc vs. Ctrl. RGC types exhibit a wide spectrum of survival at 14dpc ranging from 1–98%. F. 14d survival ranking (as in E) colored by RGC subclasses. Overlapping subclasses are denoted by two-tone color bars. G. 14d survival ranking (as in E) colored by relative abundance in control. H. Scatter plot showing correspondence between the 14dpc survival rates of RGC groups as determined by scRNAseq and IHC (R_Pearson =0.97). 26 combinations of antibodies and transgenic lines (Table S3) were used label groups of RGC types covering a broad frequency range. I. Loss of RGC somas as determined by IHC for RBPMS in this study (diamonds; see Figure S4F for example images) or by retrograde labeling from superior colliculus (triangles; redrawn from (Galindo-Romero et al., 2011)). J. -M. Each RGC type can be assigned to one of three survival groups based on the pattern of cell loss across time. Individual graphs of relative survival, defined as the fraction of cells surviving at each time point, shown for 7 resilient types (J), 11 intermediate types (K) and 27 susceptible types (L), (see also S3G). Fluctuations in sampling frequency resulted in relative survival values >1 through 2dpc (where there is little death) for rare RGC types (frequency < 0.5%). Error bars are not included for individual types in panel K-L for clarity of presentation. Grey lines, relative survival for each type within the survival group; colored lines, mean relative survival across types; shaded ribbons, standard deviation of relative survival values across types. Fluctuations observed through 2dpc were within expected error (colored ribbons), in contrast to later time points. Solid lines, mean relative survival across types within a survival group; shaded ribbons, standard deviation. Group means are superimposed in M.

Figure 4.. Physiological characteristics of resilient and susceptible RGCs

A. Representative recordings of two out of 32 channels in 1dpc and 14dpc mesh-implanted retinas. B. Sorted spike waveforms for two individual RGCs per channel (rows) represented in A recorded over multiple days. Ch1 shows spike waveforms of two sorted RGCs (purple and green lines) on 1dpc and 3dpc; cells have died by 8dpc. Ch2 shows waveforms of two sorted RGCs (blue and red) on 1dpc and 3dpc, but only one RGC was still detectable at 8dpc. C. Polar plots of responses of direction-selective (DS), orientation-selective (OS) and neither orientation- nor directions-selective (NS) RGCs to gratings moving in each of 8 directions. Each plot shows measurements from the same cell on different days. D. Proportion of RGCs by response type within each response category (columns) at 1dpc. S, sustained, T, transient, ON, OFF and ON/OFF, responds to light increments, decrements or both. E. RGC survival as a function of time in physiological recordings following ONC (black line) compared to uncrushed control (dotted line shows data replotted from (Hong et al., 2018)) F. Sustained RGCs survive better than transient RGCs as assessed by physiology (* = p<0.03 by Fisher’s Exact Test). G. OSGCs are more susceptible than DSGCs or NSGCs (* = p<0.04 by Fisher’s Exact Test). H. Among RGCs that are either OS or NS, ON-OFF cells are more susceptible than ON or OFF cells (* = p<0.03 by Fisher’s Exact Test). I. Among DSGCs, ON-OFF cells (ooDSGCs) are susceptible than ON or OFF cells (p=0.06 at 14dpc by Fisher’s Exact Test). J. Average firing rates for RGCs that survive until 14dpc or die by 8dpc. K. RGCs that are dead by 5dpc exhibit little changes in firing rate between 1–3dpc. L. RGCs that are dead by 5dpc exhibit little change in direction/orientation selectivity index (DSI/OSI) between 1–3dpc. M. En face morphology of resilient RGCs (αOFF-S, C42) at Ctrl, 4, 7, and 14dpc. N. Quantification of C42 morphological complexity (total branch points) and size (dendritic area) shows no significant difference between time points for either measure (one-way ANOVA with post-hoc Tukey HSD test). Data are shown as mean±SD. O. En face morphology of susceptible RGCs (αOFF-T, C45) at Ctrl, 3, and 4dpc. O’ showing zoomed in views of dendrites at Ctrl and 4dpc. P. Quantification of C45 morphological complexity as in N. * = p<0.04; one-way ANOVA with post-hoc Tukey HSD test. Data are shown as mean ±SD.

Figure 5.. Global changes in gene expression following injury

A. Heatmap of genes showing temporal variation following ONC. Expression values of each gene (row) is averaged across all RGCs at a given time point (columns), and then z-scored across times prior to plotting. Black bars separate genes into 8 modules (Mod) based on temporal dynamics. B. Mean temporal dynamics of individual genes (lines) from Module 1 that were associated with gene ontology (GO) biological processes related to axon and neuronal functions. Genes and the GO processes from which they were selected are listed in Table S4. C. As in B, for Modules 5 and 6 for genes associated GO biological processes related to apoptosis or various stress pathways. D. Expression dynamics of genes from B plotted for each RGC type (lines). Blue lines correspond to ipRGC types (C31, 22, 40, 33). Expression values for each type were z-scored to track relative changes. E. Same as D, for genes from C F. Expression patterns of DE genes (rows) distinguishing the 7 resRGC types and the 10 most susceptible RGC types (columns), based on 14dpc survival in the uninjured retina (Figure 3F). Values were z-scored along each row prior to plotting. G. -J. Averaged temporal dynamics of candidate genes selectively upregulated in resRGC or susRGC types (lines). Blue lines correspond the 7 resRGC types, including types that upregulate Ucn (C42, 43) or Nppb (C22, 31, 33, 40, 43) (left panels), which were not enriched for Tac1 or Cidea (right panels).

Figure 6.. Genes that affect RGC survival

A. Ucn is selectively upregulated in sustained αRGC’s (α-RGC-S; C42, 43) and Crhbp is selectively expressed in a subset of susRGC types (C14, 15, 17, 24, 26, 28, 39). Violin plots show merged expression for indicated clusters at 0 and 7dpc. The number above the violins indicates the percentage of cells expressing the marker within each subset. Box plots depict the median and interquartile range. B. FISH of retinal sections shows Ucn upregulation at 7dpc in Spp1+ RGCs (α-RGCs marker): white circles. Crhbp is expressed in a set of Spp1- RGCs (non-α-RGCs) before and after ONC: green circles. C. IHC in retinal whole mounts for RBPMS shows increased survival of RGCs at 14dpc following OE-Ucn, KO-Crhbp, or injection of UCN protein. D. Timp2 is selectively expressed in the resilient ipRGCs (C22, 31, 33, 40, 43) before and after ONC. Mmp12 is upregulated in a broad subset of susRGCs (C7, 8, 11, 12, 14, 17, 18, 23, 24, 27, 28, 43, 37, 39, 41) after crush but is low in ipRGCs in scRNAseq data. Violin plots as in A. E. FISH of retinal sections as in B. F. IHC in retinal whole mounts as in C. G. Expression in resRGC subsets at 0 and 7dpc of Ndnf (C22, 31, 43) and Prph (C31, 43). Violin plots as in panel A. H. FISH of retinal sections as in B. I. IHC of retinal whole mounts as in C. J. Total RGC survival (RBPMS+ cells; mean ± SEM) in whole mounts following interventions shown in C, F, and I. Red line and ribbon, mean RBPMS density ±SEM averaged from four sets of controls, which did not differ significantly from each other: no injection, PBS, UCN vehicle, and MMP12 inhibitor vehicle. n=18; details in STAR methods. *adjusted p-value <.05 (Bonferroni). K. IHC showing increased survival of CARTPT+ RGCs (circles) at 14dpc following OE-Ucn and OE-Timp2 compared to vehicle. Top row, CARTPT⁺ RGCs at 0dpc. L. IHC quantification showing selective survival of CARTPT⁺ RGCs (C12, 14, 16, 36) compared to NEUROD2⁺ RGCs (C12, 19, 20, 25, 26, 29, 35, 39) at 14dpc following indicated treatments. y-axis, #positive per section RGCs at 14dpc/control. Performed on retinal sagittal sections through the optic nerve. * p-value <.05 (FDR adjusted). Scale bar: 25μm for B,D,H,K; 100μm for C,F,I

Figure 7.. Genes that promote RGC axon regeneration

A. In vivo OE and KO. An AAV2 carrying the OE gene or KO sgRNA is injected intravitreally 14 days before the crush. At 12dpc regenerating axons are anterogradely labeled via CTB647 injection. UCN protein was injected at 2dpc. B. Maximum projections of cleared optic nerves showing anterograde-labeled RGC axons at 14dpc following vehicle injection or indicated treatment of Ucn (OE or protein) and KO-Crhbp (g1 and g2). C. Same as B, following OE-Timp2 and KO-Mmp9 (g1 and g2). D. Same as B, following OE-Ndnf and OE-Prph. E. -G. Quantification of axon regeneration. Control line represents mean±SEM from three groups, which did not differ significantly from each other: PBS only, UCN vehicle and AAV-Cre with no sgRNA; n=14; details in STAR methods. * p < 0.05 two-tailed Student’s t-test of area under the curve evaluated using numerical integration. Mixed effects analysis with Bonferroni correction for individual distances are shown in Table S6. In B-D, Scale bar, 250μm; X = crush site; red lines - 500, 1000, 1500μm distances from crush site.