Mingyi Zhou, John Bear, Paul A. Roberts, Filip K. Janiak, Julie Semmelhack, Takeshi Yoshimatsu, Tom Baden
Tissue preparation, immunolabeling, and imaging
For immunohistochemistry, larvae were euthanized by tricaine overdose (800 mg/l) and fixed in 4% paraformaldehyde in phosphate-buffered saline (PBS) for 30 minutes at room temperature before being washed in calcium-negative PBS. Retinae were then incubated in permeabilization/blocking buffer (PBS with 0.5% Triton X-100 and 5% normal donkey serum) at 4°C for 24 hours, and thereafter transferred to the appropriate labeling solution. For nuclear labeling, tissue was incubated at 4°C in blocking solution with Hoechst 33342 nuclear dye (Invitrogen, H21492, 1:2000) for 24 hours. For membrane staining, tissue was incubated at 4°C in blocking solution with BODIPY membrane dye (Invitrogen, C34556, 1:1000) for 24 hours. For immunostaining, tissue was incubated at 4°C for 72 hours in primary antibody solution (chicken anti-GFP (AbCam, 13970, 1:500), rabbit anti-cox iv (AbCam, 16056, 1:500), diluted in permeabilization/blocking solution). Samples were rinsed three times in PBS with 0.5% Triton X-100, then transferred to secondary antibody solution (donkey anti-chicken IgG CF488A conjugate (Sigma, SAB4600031, 1:500), donkey anti-rabbit IgG CF568 conjugate (Sigma, SAB4600076, 1:500)), diluted in permeabilization/blocking solution and incubated at 4°C for 24 hours. Finally, samples were rinsed three times in PBS with 0.5% Triton X-100 before being mounted in mounting media (VectaShield, Vector, H-1000) for confocal imaging.
]. Briefly, whole retinas were fixed in 2% PFA /2% glutaraldehyde for 24 hours at 4°C, rinsed in PBS, treated with 0.1% sodium borohydride (NaBH4) in 0.2% Triton X-100 in PBS for 10 minutes at room temperature, and rinsed again to remove excess NaBH4. For immunolabeling, all steps are as described above, with the following exceptions: blocking buffer consisted of 10% normal donkey serum, 0.1% Tween-20, and 0.5% Triton X-100 in PBS; primary and secondary antibodies were also diluted in this blocking buffer.
Confocal stacks and individual images were taken on Leica TCS SP8 using 40x water-immersion objective at xy resolution of 2,048×2,048 pixels (pixel width: 0.162 μm). Voxel depth of stacks was taken at z-step 0.3-0.5 μm. Contrast and brightness were adjusted in Fiji (NIH).
Cell density mapping
] using custom-written scripts in Igor Pro 6.37 (Wavemetrics). The density map of RGC somata was computed by subtracting the density map of dACs from that of GCL cells. Similarly, the density map of ACs was computed by summing the density maps of dACs and ACs from the inner nuclear layer. From here, RGC maps were also mapped into a sinusoidal projection of visual space [
The lipophilic tracer dye DiO (Invitrogen, D307) was used to trace RGC axons from the retina to their arborization fields in the pretectum and tectum. 1 mg/mL stock solution was prepared in dimethylformamide and stored at −20°C. For injection into Tg(Islet2b:nls-trpR, tUAS:MGCamp6f) retinas, the lenses of whole fixed larvae were removed and a sufficient amount of tracer dye injected into one of either the left or the right eye so as to completely cover the exposed surface of the GCL. Tissue was then incubated at 37°C for 3 days to allow the dye time to diffuse all the way up RGC axons to their terminals in the midbrain.
Prior to photoconversion, 6-8 dpf Islet2b:PA-GFP larvae were injected with BODIPY membrane dye (1 nL of 1 mg/mL; Sigma, D3821) into the space behind the right eye and underlying skin to demarcate retinal anatomy and facilitate subsequent targeting. Larvae were left for 10-20 minutes at 25°C to allow the dye to diffuse into the retina. After 20 minutes, the IPL was uniformly stained, and the individual somata of GCL neurons showed nuclear exclusions which were used for subsequent targeting.
Cells were photoconverted under the same 2-photon microscope as used for functional imaging (below). In each animal, we randomly photoconverted 2-5 cells per eye in the nasal retina and/or strike zone, with a minimum spacing of 30 μm between them. For photoactivation, the femtosecond laser was tuned to 760 nm and focused onto one single soma at a time for up to ∼2 minutes. After a typically > 40 minutes cells were visualized under 2-photon (927 nm) and imaged in a 512×512 pixel (1 μm z-steps) stack which encompassed each cell’s soma, axon initial segment, and the entirety of the dendritic structure. Throughout, the BODIPY signal was included as an anatomical reference.
Two-photon functional imaging and stimulation parameters
]; purchased through Sutter Instruments/Science Projects) equipped with the following: a mode-locked Ti:Sapphire laser (Chameleon Vision-S, Coherent) tuned to 927 nm for imaging GFP and 960 nm for imaging mCherry/BODIPY in combination with GFP; two fluorescent detection channels for GFP (F48x573, AHF/Chroma) and mCherry/BODIPY (F39x628, AHF/Chroma), and; a water-immersion objective (W Plan-Apochromat 20x/1,0 DIC M27, Zeiss). For image acquisition, we used custom-written software (ScanM, by M. Mueller, MPI, Martinsried and T Euler, CIN, Tübingen) running under Igor Pro 6.37 (Wavemetrics). Structural data was recorded at 512×512 pixels, while functional data was recorded at 64×32 pixel resolution (15.6 Hz, 2 ms line speed). For each functional scan, we first defined a curvature of the imaged IPL segment based on a structural scan, and thereafter “bent” the scan plane accordingly (“banana scan”). This ensured that the imaging laser spent a majority of time sampling from the curved IPL and INL, rather than adjacent dead space. The banana-scan function was custom-written under ScanM.
] to yield ‘natural white’: red, “100%” (34×105 photons /s /cone); green, “50%” (18 x105 photons /s /cone); blue, “13%” (4.7 x105 photons /s /cone); ultraviolet, “6%” (2.1×105 photons /s /cone). We did not compensate for cross-activation of other cones. Owing to 2-photon excitation of photopigments, an additional constant background illumination of ∼104 R∗ was present throughout [
]. For all experiments, larvae were kept at constant illumination for at least 2 s after the laser scanning started before light stimuli were presented. Two types of full-field stimuli were used: a binary dense “natural spectrum” white noise, in which the four LEDs were flickered independently in a known random binary sequence at 6.4 Hz for 258 s, and a natural-white chirp stimulus [
] where all four LEDs were driven together. To prevent interference of the stimulation light with the optical recording, LEDs were synchronized to the scanner’s retrace [
Quantification and Statistical Analysis
No statistical methods were used to predetermine sample size.
Data analysis was performed using IGOR Pro 6.3 (Wavemetrics), Fiji (NIH) and MATLAB R2018b (Mathworks).
ROI placements and quality criterion
]. To allocate ROIs to dendritic and somatic datasets a boundary between the GCL and IPL was drawn by hand in each scan – all ROIs with a center of mass above the boundary were considered as dendritic, and all ROIs below were considered as somatic. Since the lower part of the IPL tends to be dominated by On-circuits, it is possible that a small number of On-dendrites were incorrectly classed as somata which may go part-way to explaining the generally stronger On-bias among somatic compared to dendritic ROIs (e.g., Figure 2A). Moreover, due to the ring-like nature of mGCaMP6f expression profiles in somata when optically sectioned, it was possible that two ROIs could be inadvertently placed on different halves of the same soma. However, since whether or not a soma was split in this way was likely non-systematic over functional types, we did not attempt to correct for this possibility. Only ROIs where at least one of the four spectral kernels’ peak-to-peak amplitudes exceeded a minimum of ten standard deviations were kept for further analysis (n = 2,414/2,851 dendritic ROIs, 84.7%; 411/796 somatic ROIs, 51.6%). Equally, all individual color kernels that did not exceed 10 SDs were discarded (i.e set to NaN).
A note on ROI segmentation and identity
]. Specifically, the high density and overlap of dendritic processes across the IPL means that it is impossible to tell if groups of dendritic ROIs belong to the same RGC (Figures 1D and 1E). Nevertheless, functional dendritic data is indicative of the local computations that occur within RGC dendrites as they integrate signals from BCs and ACs in different layers of the IPL and in different positions of the eye [
]. Further, our single cell data (Figures S2A–S2E) suggests that dendritic signals in population recordings are probably also a reasonable proxy for somatic signals, with the added benefit that their signal-to-noise was generally higher (e.g., Figure 2A). To what extent the indicated close similarity of dendritic and somatic signals in zebrafish RGCs applies across all RGC types, and to what extent this can be linked to their generally small absolute size (e.g., compared to RGCs in larger eyes), will be important to address in the future.
]), ought to be considered with some caution. For example, determining a binary value for a kernel’s polarity (On or Off) can be conflicted with the fact that a neuron might exhibit both On and Off response aspects. Moreover, different possible measures of On or Off dominance in a kernel can generate different classification biases. Here, we defined On and Off based on a measure of a kernel’s dominant trajectory in time. For this, we determined the position in time of each kernel’s maximum and minimum. If the maximum preceded the minimum, the kernel was classified as Off, while vice versa if the minimum preceded the maximum, the kernel was defined as On. Examples On and Off kernels classified in this way can for example be seen in Figure 3B (cf. Figure 3A central horizontal column for a lookup of how each kernel was classified).
Digitizing photoactivated cells
Quantifying dendritic tilt
] for ϕ, since this variable is (2π-)periodic. In comparing SZ and nasal RGCs, the dendritic CoM positions, r, are predicted to be from different distributions (p = 0.0209, 3 s.f.); the dendritic tilt strengths, θ, are predicted to be from the same distribution (p = 0.894, 3 s.f.); and the dendritic tilt angles, ϕ, are predicted to be from different distributions (p = 0.001).
The morphological data consists of sets of points in three-dimensional Cartesian coordinates (x,y,z) describing the dendritic architecture for each of 131 RGCs, 67 from the nasal (N) region and 64 from the strike zone (SZ) region. The coordinate axes are orientated such that the y axis is perpendicular to the plane of the retina, spanning the width of the IPL, while the x and z axes are tangential to the plane of the retina. The coordinates in the y-dimension are scaled so as to lie in the interval [0,10] for any processes within the IPL, and > 10 or < 0 for INL and GCL processes (where applicable), respectively. The position of the soma, which always lay in the GCL, was not used for clustering.
Three summary statistics, each of which capture some aspect of the dendritic architecture, were defined for use in clustering: i) y_span: the width of the dendritic tree in the y-direction; ii) y_mean: the mean position of the points in the dendritic tree in the y-direction; and iii) num_pts: the number of points in the dendritic tree. While we experimented with other summary statistics, these three were found to be sufficient to differentiate the RGCs into their basic morphological groups.
We also defined one further summary statistic: iv) xz_area: the area spanned by the dendritic tree in the xz-plane, calculated as the convex hull using the MATLAB routine convhull. This statistic was not used for clustering since the information contained in xz_area is largely captured between y_span and num_pts. While not required for clustering, this summary statistic nonetheless captures important characteristics of the dendritic morphology and hence is represented in the results section alongside y_span, y_mean and num_pts.
Each of the summary statistics was standardized by subtracting the mean and dividing by the standard deviation. In this way, we ensured that each of the summary statistics was equally weighted by the clustering algorithm.
Clustering was performed in two stages, using agglomerative hierarchical clustering in both cases. The first stage of clustering used all three summary statistics (y_span, y_mean and num_pts), splitting the data into 18 clusters. Two of the resulting clusters were large and contained a variety of morphologies as discerned from visual inspection. These clusters were split further via a second round of clustering, using just the y_span summary statistic. The first cluster was split into 6 subclusters and the second into 3 subclusters, resulting in a total of 25 clusters, where the 13 clusters containing a minimum of 4 members were included for presentation.
Hierarchical clustering was performed using the MATLAB routines pdist, linkage and cluster. The function pdist calculates the distances between each RGC in (y_span,y_mean,num_pts)-space, while the function linkage operates on the output of the pdist routine to encode an agglomerative hierarchical cluster tree. There are a number of options for defining the distances between RGCs for pdist and the distances between clusters for linkage. We used the ‘city block’ distance metric for pdist and the ‘average’ distance metric for linkage as, in general, these were found to result in a larger cophenetic correlation coefficient (CCC) than any other combination of distance metrics. The CCC is a measure of the fidelity with which the cluster tree represents the dissimilarities between observations. It was calculated using the MATLAB routine cophenet and takes values between [-1,1], where values closer to positive unity represent a more faithful clustering. In the results presented here, the first stage of clustering had a CCC of 0.77 (2 d.p.), while the two subclusterings in the second stage had CCCs of 0.77 (2 d.p.) and 0.83 (2 d.p.).
Lastly, RGCs were assigned to clusters using the MATLAB routine cluster. The number of clusters was determined by specifying a cutoff distance which was chosen following visual inspection of the cluster tree dendrogram so as to respect a natural division in the data.
Functional data pre-processing and receptive field mapping
]. Next, the Ca2+ traces for each ROI were extracted and de-trended by high-pass filtering above ∼0.1 Hz and followed by z-normalization based on the time interval 1-6 s at the beginning of recordings using custom-written routines under IGOR Pro. A stimulus time marker embedded in the recording data served to align the Ca2+ traces relative to the visual stimulus with a temporal precision of 1 ms. Responses to the chirp stimulus were up-sampled to 1 KHz and averaged over 3-6 trials. For data from tetrachromatic noise stimulation we mapped linear receptive fields of each ROI by computing the Ca2+ transient-triggered-average. To this end, we resampled the time-derivative of each trace to match the stimulus-alignment rate of 500 Hz and used thresholding above 0.7 standard deviations relative to the baseline noise to the times ti at which Calcium transients occurred. We then computed the Ca2+ transient-triggered average stimulus, weighting each sample by the steepness of the transient:
Here, S(l,t) is the stimulus (“LED” and “time”), τ is the time lag (ranging from approx. −1,000 to 350 ms) and M is the number of Ca2+ events. RFs are shown in z-scores for each LED, normalized to the first 50 ms of the time-lag. To select ROIs with a non-random temporal kernel, we used all ROIs with a standard deviation of at least ten in at least one of the four spectral kernels. The precise choice of this quality criterion does not have a major effect on the results.
To summarize average functions of RGC processes across different positions in the eye and across IPL depths, we computed two-dimensional “Eye-IPL” maps. For this, we divided position in the eye (-π:π radians) into eight equal bins of width π/4. Similarly, we divided the IPL into 20 bins. All soma ROIs were allocated to bin 1 independent of their depth in the GCL. while all IPL ROIs were distributed to bins 3:20 based on their relative position between the IPL boundaries. As such, bin 2 is always empty, and serves as a visual barrier between IPL and GCL. From here, the responses of ROIs within each bin were averaged. All maps were in addition smoothed using a circular π/3 binomial (Gaussian) filter along eye-position, as well as for 5% of IPL depth across the y-dimension (dendritic bins 3:20 only).
On-Off index (OOi)
Where nOn and nOff correspond to the number of On and Off kernels in a bin, respectively. OOi ranged from 1 (all kernels On) to −1 (all kernels Off), with and OOi of zero denoting a bin where the number of On and Off kernels was equal.
Ternary response classification
Each ROI was allocated to one of 81 ternary response bins (three response states raised to the power of four spectral bands). One of three response-states was determined for each of four spectral kernels (red, green, blue, UV) belonging to the same ROI: On, Off or non-responding. All kernels with a peak-to-peak amplitude below ten standard deviations were considered non-responding, while the remainder was classified as either On or Off based on the sign of the largest transition in the kernel (upward: On, downward: Off).
Feature extraction and Clustering
Clustering was performed on four datasets, each containing the functional responses of RGCs to chirp stimuli and kernels derived from color noise stimuli: 1) pan retinal inner plexiform layer (PR-IPL) dataset (n = 2,851), sampling RGC dendritic responses at all eccentricities and across a range of depths in the IPL; 2) strike zone inner plexiform layer (SZ-IPL) dataset (n = 3,542), sampling RGCs at the SZ only and across the IPL; 3) pan retinal ganglion cell layer (PR-GCL) dataset (n = 796), sampling RGC responses at all eccentricities from the RGC somata in the GCL; and 4) strike zone ganglion cell layer (SZ-GCL) dataset (n = 1,694), sampling RGCs at the SZ only from the RGC somata. Mean responses to chirp stimuli were formatted as 2,499 time points (dt = 1 ms) while color kernels were formatted as 649 time points (dt = 2 ms, starting at t = −0.9735 s) per spectral channel (red, green, blue and UV).
For each dataset we clustered using only the kernels portion of the data since this was found to produce a cleaner clustering than when clustering chirp responses and kernels together, or chirp responses alone. ROIs with low quality kernels, determined as the maximum standard deviation across the four colors, were identified and removed from the dataset. For clustering, a kernel quality threshold of 5 was chosen, such that any ROI with a kernel quality below this threshold was eliminated from the data to be clustered.
Following quality control, the datasets had the following sizes: 1) PR-IPL: n = 2,414 (84.7% of original); 2) SZ-IPL: n = 2,435 (68.8% of original); 3) PR-GCL: n = 411 (51.6% of original); 4) SZ-GCL: n = 721 (42.6% of original).
We scaled the data corresponding to each kernel color by dividing each one by the standard deviation through time and across ROIs. In this way we ensured an even weighting for each color. This is important, since the red and green kernels tended to have larger amplitudes than the blue and UV kernels.
We used principal component analysis (PCA) to reduce the dimensions of the problem prior to clustering. PCA was performed using the MATLAB routine pca (default settings). We applied PCA to the portions of a dataset corresponding to each of the kernel colors separately, retaining the minimum number of principal components necessary to explain ≥ 99% of the variance. The resulting four ‘scores’ matrices were then concatenated into a single matrix ready for clustering. The following numbers of principal components were used for each of the four datasets: 1) PR-IPL: 8 red (R) components, 8 green (G) components, 13 blue (B) components, 33 ultraviolet (UV) components (62 in total); 2) SZ-IPL: 15 R, 17 G, 25 B, 18 UV (75 in total); 3) PR-GCL: 13 R, 11 G, 24 B, 36 UV (84 in total); and 4) SZ-GCL: 20 R, 21 G, 27 B, 34 UV (102 in total).
We clustered the combined ‘scores’ matrix using Gaussian Mixture Model (GMM) clustering, performed using the MATLAB routine fitgmdist. We clustered the data into clusters of sizes 1,2,…,100, using i) shared-diagonal, ii) unshared-diagonal, iii) shared-full and iv) unshared-full covariance matrices, such that (100∗4 = ) 400 different clustering options were explored in total. For each clustering option 20 replicates were calculated (each with a different set of initial values) and the replicate with the largest loglikelihood chosen. A regularization value of 10−5 was chosen to ensure that the estimated covariance matrices were positive definite, while the maximum number of iterations was set at 104. All other fitgmdist settings were set to their default values.
In datasets PR-IPL and SZ-IPL the optimum clustering was judged to be that which minimized the Bayesian information criterion (BIC), which balances the explanatory power of the model (loglikelihood) with model complexity (number of parameters), while clusters with < 10 members were removed. In datasets PR-GCL and SZ-GCL the BIC did not give a clean clustering; therefore, we specified 20 clusters for the PR-GCL and 10 clusters for the SZ-GCL, with unshared-diagonal covariance matrices, removing clusters with < 5 members.
Using the above procedure, we obtained the following optimum number of clusters for each dataset: 1. PR-IPL: 15 clusters (2 clusters with < 10 members removed); 2. SZ-IPL: 12 clusters (1 cluster with < 10 members removed); 3. PR-GCL: 13 clusters (7 clusters with < 5 members removed); 4. SZ-GCL: 9 clusters (1 cluster with < 5 members removed). Unshared-diagonal covariance matrices gave the optimal solution in all cases.