No statistical methods were used to predetermine sample size.
All procedures were performed in accordance with the law on animal protection issued by the German Federal Government (Tierschutzgesetz) and approved by the institutional animal welfare committee of the University of Tübingen. For all experiments, we used 4- to 8-week-old mice of either sex. In addition to C57Bl6 (wild-type) mice, we used the transgenic lines PvalbCre (‘PV’, JAX 008069, The Jackson Laboratory; ref. 43), Pcp2Cre (‘Pcp2’, JAX 006207; ref. 44) and ChATCre (JAX 006410; ref. 51), cross-bred with the red fluorescence Cre-dependent reporter line Ai9:tdTomato (JAX 007905). Owing to the exploratory nature of our study, we did not use randomization and blinding.
Animals were housed under a standard 12 h day/night rhythm. For activity recordings, animals were dark-adapted for ≥1 h, then anaesthetized with isoflurane (Baxter) and killed by cervical dislocation. The eyes were enucleated and hemisected in carboxygenated (95% O , 5% CO ) artificial cerebral spinal fluid (ACSF) solution containing (in mM): 125 NaCl, 2.5 KCl, 2 CaCl , 1 MgCl , 1.25 NaH PO , 26 NaHCO , 20 glucose, and 0.5 l-glutamine (pH 7.4). Bulk electroporation of the fluorescent Ca2+ indicator Oregon-Green BAPTA-1 (OGB-1) into the ganglion cell layer (GCL) was carried out as described before19, 47. In brief, the retina was dissected from the eyecup, flat-mounted onto an Anodisc (#13, 0.2 μm pore size, GE Healthcare) with the GCL facing up, and placed between two 4-mm horizontal plate electrodes (CUY700P4E/L, Nepagene/Xceltis). A 10 μl drop of 5 mM OGB-1 (hexapotassium salt; Life Technologies) in ACSF was suspended from the upper electrode and lowered onto the retina. After application of 10–12 pulses (+9 V, 100 ms pulse width, at 1 Hz) from a pulse generator/wide-band amplifier combination (TGP110 and WA301, Thurlby handar/Farnell), the tissue was moved to the recording chamber of the microscope, where it was continuously perfused with carboxygenated ACSF at ~37 °C and left to recover for ~60 min before the recordings started. In all experiments with wild-type mice, ACSF contained ~0.1 μM Sulforhodamine-101 (SR101, Invitrogen) to reveal blood vessels and any damaged cells in the red fluorescence channel20. All procedures were carried out under very dim red (>650 nm) light.
We used a MOM-type two-photon microscope (designed by W. Denk, MPI, Martinsried; purchased from Sutter Instruments/Science Products). Design and procedures were described previously20. In brief, the system was equipped with a mode-locked Ti:Sapphire laser (MaiTai-HP DeepSee, Newport Spectra-Physics) tuned to 927 nm, two fluorescence detection channels for OGB-1 (HQ 510/84, AHF/Chroma) and SR101 (HQ 630/60, AHF), and a water immersion objective (W Plan-Apochromat 20x/1,0 DIC M27, Zeiss). For image acquisition, we used custom-made software (ScanM, by M. Müller, MPI, Martinsried, and T.E.) running under IGOR Pro 6.3 for Windows (Wavemetrics), taking 64 × 64 pixel image sequences (7.8 frames per s) for activity scans or 512 × 512 pixel images for high-resolution morphology scans.
For light stimulation, we focused a DLP projector (K11, Acer) through the objective, fitted with band-pass-filtered light-emitting diodes (LEDs) (‘green’, 578 BP 10; and ‘blue’, HC 405 BP 10, AHF/Croma) that roughly match the spectral sensitivity of moose M- and S-opsins. LEDs were synchronized with the microscope’s scan retrace. Stimulator intensity (as photoisomerisation rate, 103 P* per s per cone) was calibrated as described previously52 to range from 0.6 and 0.7 (black image) to 18.8 and 20.3 for M- and S-opsins, respectively. Owing to two-photon excitation of photopigments, an additional, steady illumination component of ~104 P* per s per cone was present during the recordings (for detailed discussion, see ref. 22). For all experiments, the tissue was kept at a constant intensity level (see stimuli below) for at least 30 s after the laser scanning started before light stimuli were presented. Four types of light stimulus were used (Fig. 1b, top): (i) full-field ‘chirp’ stimulus consisting of a bright step and two sinusoidal intensity modulations, one with increasing frequency and one with increasing contrast; (ii) 0.3 × 1 mm bright bar moving at 1 mm s−1 in eight directions19; (iii) alternating blue and green 3-s flashes; and (iv) binary dense noise, a 20 × 15 matrix with 40 μm pixel-side length; each pixel displayed an independent, perfectly balanced random sequence at 5 Hz yielding a total running time of 5 min for receptive field (RF) mapping. In some experiments, we used in addition dark moving bars (like (ii) but contrast-inversed), and stationary bright or dark 0.2 × 0.8 mm bars flashed for 1 s in six orientations (see Extended Data Fig. 7h–s). Except for (iii), stimuli were achromatic, with matched photo-isomerization rates for M- and S-opsins.
TdTomato- or OGB-1-labelled RGCs were targeted using two-photon imaging for juxtacellular recordings using borosilicate electrodes (4–6 MΩ) filled with ACSF with added SR101 (250 μM). Data were acquired using an Axoclamp-900A or Axopatch 200A amplifier in combination with a Digidata 1440 (all: Molecular Devices), digitized at 10 kHz and analysed offline using IGOR Pro. We presented the same light stimuli as for the Ca2+ imaging. In some experiments, electrical recordings and Ca2+ imaging was performed simultaneously. After the recording, the membrane under the electrode was opened using a voltage ‘buzz’ to let the cell fill with dye by diffusion for approximately 30 min; then two-photon image stacks were acquired to document the cell’s morphology.
Filled cells were traced semi-automatically offline using the Simple Neurite Tracer plugin implemented in Fiji (http://fiji.sc/Fiji), yielding cell skeletons. If necessary, we used the original image stack to correct the skeletons for any warping of the IPL using custom-written scripts in IGOR Pro. To this end, we employed the SR101-stained blood vessel plexuses on either side of the IPL as landmarks to define the IPL borders (see below). Only cells where the filling quality allowed full anatomical reconstruction were used for analysis (see below).
Following Ca2+ imaging, retinas were mounted onto filter paper (0.8 μm pore size, Millipore) and fixed in 4% paraformaldehyde (in PBS) for 15 min at 4 °C. Immunolabelling was performed using antibodies against ChAT (choline-acetyltransferase; goat anti-ChAT, 1:100, AB144P, Millipore), GAD67 (glutamate decarboxylase; mouse anti-GAD67, 1:100, MAB5046, Millipore), SMI-32 (mouse anti-SMI32, 1:100, SMI-32R-100, Covance), and melanopsin (rabbit anti-melanopsin, 1:4000, AB-N38, Advanced Targeting Systems) for 4 days. Secondary antibodies were Alexa Fluor conjugates (1:750, 16 h, Life Technologies). For each retina, the recorded region was identified by the local blood vessel pattern and confirmed by comparing size and position of individual somata in the GCL. Image stacks were acquired on a confocal microscope (Nikon Eclipse C1) equipped with a ×60 oil objective (1.4 NA). The degree of immunolabelling of GCL cells was evaluate and rated (from 0, negative, to 4 positive) using z-stacks. Attribution of labelled somata to recorded ones was performed manually using ImageJ (http://imagej.nih.gov/ij) and IGOR Pro.
Data analysis was performed using Matlab 2012 and 2014a (The Mathworks Inc.), and IGOR Pro. Data were organized in a custom written schema using the DataJoint for Matlab framework (http://datajoint.github.io/; D. Yatsenko, Tolias lab, Baylor College of Medicine).
Regions of interest (ROIs), corresponding to somata in the GCL, were defined semi-automatically by custom software (D. Velychko, CIN) based on a high-resolution (512 × 512 pixels) image stack of the recorded field. Then, the Ca2+ traces for each ROI were extracted (as ΔF/F) using the IGOR Pro-based image analysis toolbox SARFIA (http://www.igorexchange.com/project/SARFIA). 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 2 ms. Stimulus-aligned Ca2+ traces for each ROI were imported into Matlab for further analysis.
First, we de-trended the Ca2+ traces by high-pass filtering above ~0.1 Hz. For all stimuli except the dense noise (for RF mapping), we then subtracted the baseline (median of first eight samples), computed the median activity r(t) across stimulus repetitions (typically three to five repetitions) and normalized it such that .
We mapped the linear RFs of the neurons by computing the Ca2+ transient-triggered average. To this end, we resampled the temporal derivative of the Ca2+ response ċ(t) at 10-times the stimulus frequency and used Matlab’s findpeaks function to detect the times t at which Ca2+ transients occurred. We set the minimum peak height to 1 s.d., where the s.d. was robustly estimated using:
We then computed the Ca2+ transient-triggered average stimulus, weighting each sample by the steepness of the transient:
Here, S(x, y, t) is the stimulus, τ is the time lag (ranging from approximately −320 to 1,380 ms) and M is the number of Ca2+ events. We smoothed this raw RF estimate using a 5 × 5 pixel Gaussian window for each time lag separately. RF maps shown correspond to a s.d. map, where the s.d. is calculated over time lags τ:
To extract the RF’s position and scale, we fitted it with a 2D Gaussian function using Matlab’s lsqcurvefit. The time course of the receptive field F (τ) was estimated by the average of the eight pixels closest to the fitted RF centre (according to the Mahalanobis distance) weighted by a Gaussian profile. RF quality (QI ) was measured as one minus the fraction of variance explained by the Gaussian fit
To extract time course and directional tuning of the Ca2+ response to the moving bar stimulus, we performed a singular value decomposition (SVD) on the T by D normalized mean response matrix M (times samples by number of directions; T = 32; D = 8; Extended Data Fig. 7a, b):
This procedure decomposes the response into a temporal component in the first column of U and a direction dependent component or tuning curve in the first column of V, such that the response matrix can be approximated as an outer product of the two:
An advantage of this procedure is that it does not require manual selection of time bins for computing direction tuning, but extracts the direction tuning curve given the varying temporal dynamics of different neurons.
To measure direction selectivity (DS) and its significance, we projected the tuning curve V on a complex exponential , where α is the direction in the kth condition:
This is mathematically equivalent to computing the vector sum in the 2D plane or computing the power in the first Fourier component. We computed a DS index as the resulting vector length
correcting for the direction spacing. We additionally assessed the statistical significance of direction tuning using a permutation test53. To this end, we created surrogate trials (that is, stimulus repetitions) by shuffling the trial labels (that is, destroying any relationship between condition and response), computed the tuning curve for each surrogate trial and projected it on the complex exponential ϕ. Carrying out the procedure 1,000 times generated a null distribution for K, assuming no direction tuning. We used the percentile of the true K as the P value for direction tuning (Extended Data Fig. 7c). Importantly, a large DSi does not necessarily result in a small P value, for example, in the case of large trial to trial variability. As a result, the DSi distributions of significantly and not significantly direction tuned neurons show substantial overlap (Extended Data Fig. 7d, e). Therefore, a simple threshold as a DS criterion (for example, DSi > 0.4) does not provide a good separation into direction selective cell types and others.
Orientation selectivity (OS) was assessed in an analogous way. However, we used the complex exponential , corresponding to the second Fourier component.
To measure how well a cell responded to a stimulus (chirp, moving bar, colour), we computed the signal-to-noise ratio
where C is the T by R response matrix (time samples by stimulus repetitions) and 〈 〉 and Var[ ] denote the mean and variance across the indicated dimension, respectively. If all trials are identical such that the mean response is a perfect representative of the response, QI is equal to 1. If all trials are completely random with fixed variance (so that the mean response is not informative about the individual trial responses at all), QI is proportional to 1/R.
For further analysis, we used only cells that responded well to the chirp and/or the moving bar stimulus (QI > 0.45 or QI > 0.6).
The full-field index was computed as
comparing the response quality to a local stimulus (moving bar) and a global stimulus (chirp).
where r and r are defined as the activity during the response to the leading edge of the moving bar (the first 400 ms of the ON response) and the trailing edge of the moving bar (the first 400 ms of the OFF response).
Colour selectivity was measured for the ON response using
and for the OFF response using an analogous definition. Here, r and r are the responses in a time window of 1,280 ms after onset of the green and blue stimulus, respectively.
We used sparse principle component analysis54, as implemented in the SpaSM toolbox by K. Sjöstrang et al. (http://www2.imm.dtu.dk/projects/spasm/), to extract sparse response features from the responses to the chirp, colour, and moving bar stimulus, resulting in features which use only a small number of time bins. The extracted features are localized in time and therefore readily interpretable (for example, ‘high-frequency feature’), although this constraint was not explicitly enforced by the algorithm (Extended Data Fig. 2e). We also explored alternative feature extraction techniques such as regular PCA, but these resulted in inferior cluster quality. In addition, they required manually defining regions corresponding to specific parts of the stimulus (for example, frequency chirp) to yield localized and interpretable features.
We extracted 20 features with 10 non-zero time bins from the mean response to the chirp (averaging across trials) and 6 features with 10 non-zero time bins from the mean response to the colour stimulus. For the moving bar stimulus, we extracted 8 features with 5 non-zero time bins from the response time course (see above) and 4 features with 6 non-zero time bins from its temporal derivative. All features were in the temporal domain, ensuring spatial invariance. In addition, we used two features from the time course of the RF, extracted with regular PCA. Overall, this procedure resulted in a 40 dimensional feature vector for each cell. Before clustering, we standardized each feature separately across the population of cells.
DS and non-DS cells were clustered independently, classifying cells as DS if the permutation test resulted in P < 0.05 (see above). We fit each data set with a Mixture of Gaussians model using the expectation-maximization algorithm (Matlab’s gmdistribution object). We constrained the covariance matrix for each component to be diagonal, resulting in 80 parameters per component (40 for the mean, 40 for the variances). We further regularized the covariance matrix by adding a constant (10−5) to the diagonal. To find the optimal number of clusters, we evaluated the Bayesian information criterion55
where L is the log-likelihood of the model, N is the number of cells and M is the number of parameters in the model, that is, M = 81C − 1 where C is the number of clusters and the contributions that arose from means, variances and mixture proportions (which have to add to 1). Although other choices such as the Aikaike information criterion (AIC) would have been possible, we found the BIC to yield a good compromise between model complexity and quality, since the AIC is known to find too many clusters for large sample sizes. We also computed log Bayes factors as 2ΔBIC for each candidate cluster number to test how strong the evidence for further splitting is. Values >6 were treated as strong evidence in favour of further splitting. The minimum of the BIC coincided well with the number of clusters after which there was no strong evidence for further adding more clusters. To avoid local minima, we restarted the EM algorithm 20 times per candidate number of clusters and used the solution with the largest likelihood. This procedure resulted in 24 and 48 clusters for DS and non-DS cells, respectively (Extended Data Fig. 2a).
To evaluate cluster quality, we rank-ordered the posterior probabilities for cluster assignment for each cluster, normalized for cluster size and averaged across clusters for non-DS and DS cells separately (Extended Data Fig. 2b). The steep decays of the sigmoidal functions indicate good cluster separability. To check how consistent the clustering was against subsampling of the data, we created 20 surrogate data sets containing random selections of 90% of the cells. We fit these surrogate data sets with a Mixture of Gaussians model with the optimal number of clusters determined on the original data set. For each cluster mean in these models, we computed the correlations with the most similar cluster for the model fit on the original data set. To summarize the similarity of clusterings, we computed the median correlation across clusters (Extended Data Fig. 2c). On average, the clusterings obtained on the surrogate and the original data set were very similar (mean median correlation: 0.96 ± 0.19 and 0.97 ± 0.01; mean ± s.d.; for DS and non-DS cells separately).
In addition, we performed an alternative clustering version, where we did not split the data in DS and non-DS cells but added DSi, OSi, soma and receptive field size as features. The identified clusters were very similar, but this strategy failed to identify most DS types as separate clusters, except for the ON–OFF DS cell. Therefore, we decided to first isolate significant DS cells and cluster them separately, before merging similarly responding DS and non-DS clusters (see below), if we did not find a reason to keep the DS group as a separate RGC type. Nevertheless, a strategy equally justified as ours could start with the alternative clustering and then split those clusters containing large fractions of DS cells.
A subset of cells was stained against GAD67 to identify dACs (see above). The intensity of this staining was manually rated as follows: −2 (absent), −1 (probably absent), 0 (uncertain), 1 (probably present), and 2 (present). For each cluster, we computed the average staining from the labelled cells (average number of cells with GAD67 information per cluster: 16.8 ± 10.0, mean ± s.d.). Clusters with an average staining <-0.2 were labelled RGCs (n = 30 clusters), those with average staining >0.2 were labelled ACs (n = 26). Clusters with average staining in-between those values (n = 5), or those that contained 6 or less cells with GAD67 information (n = 8) were labelled as uncertain, unless other clear criteria such as soma size or genetic labels indicated that they are ACs or RGCs. In this case they were manually allocated to RGC or AC (n = 3 and n = 2, respectively). Two clusters automatically classified as AC were included in the uncertain group due to their functional similarity with the OFF-suppressed types (G ). This procedure resulted in 33 RGC clusters, 10 uncertain clusters and 26 AC clusters.
We extracted all cells with large cell bodies (>136 μm2; mean + 1 s.d. of total soma size distribution; Extended Data Fig. 2i, j) from RGC and uncertain clusters. Predominantly, these cells had been assigned to nine of the clusters. We re-clustered those cells using a Mixture of Gaussians model as described above, resulting in 16 clusters (Extended Data Fig. 2j). Receptive field size was not used in this process. Five of these clusters could be clearly associated with the three known alpha-RGC types and their response profiles31 (trans. OFF alpha, 2; sust. OFF alpha, 2; ON alpha, 1). Cells in these clusters were SMI-32-positive, as expected from alpha RGCs (Fig. 3i, k). Probably, this procedure missed some alpha cells, as somata close to the edge of scan fields were cut and we thus underestimate the soma size of these cells (for example, G c, see Fig. 2a–c). Remaining cells were kept in their original cluster. Logistic regression was used to assess the effect cell type (alpha vs. mini) on SMI-32 staining (absent vs. present). We used the Matlab implementation fitglm with a binomial nonlinearity. 95%-confidence intervals on the proportion of SMI32-positive cells were computed using bootstrapping with 1,000 samples.
We used a standard linkage algorithm on the means of the RGC groups in the standardized feature space with correlation distance and average unweighted distance and plotted the result as a dendrogram (using Matlab functions linkage and dendrogram). The leaf order was optimized using the Matlab function optimalleaforder and modified for clarity of presentation.
with the number of cells in a group (n ), the median RF size (A ) within a group counting only cells that surpassed a RF quality criterion of 0.3, and the total scan area across all experiments (A ). We corrected n for 29% cells discarded by our quality criterion as well as an empirically estimated 8% of cells that did not yield a ROI in the first place due to weak or absent labelling. In addition, A was corrected for an empirically estimated 34.8% RF overhang (that is, where a cell’s RF exceeds the scan field edge). This procedure yielded a CF of 2.0 ± 0.7 for most RGC groups (Gaussian fit; see Fig. 2e, right). However, differences between studies in approaches to measure RFs (for example, checkerboards vs. bars), in the assumptions used for RF fitting (for example, homogenous RFs best fitted by Gaussians), or in the methods to estimate dendritic arbor area can easily yield different absolute estimates of CF (see also Supplementary Table 1).
To determine a cell’s IPL stratification profile, we calculated dendritic density as described previously10 with spatial smoothing of 1 μm3. The resultant 3D density cloud was projected on the z axis to estimate the mean IPL depth profile. The relationship between the depth profiles and the two ChAT bands was estimated in independent experiments using mice that express tdTomato in cholinergic ACs (ChATCre × Ai9:tdTomato). We compared the IPL depth of the tdTomato-labelled dendritic plexi to the two SR101-labelled blood vessel planes that line the inner retina. We estimated the error to be ~1.5 μm (s.d.), corresponding to 3–4% IPL depth (n = 13 measurements in 2 mice).
To relate each cell’s IPL profile to functional groups we calculated the mean correlation coefficient between a cell’s response to the chirp and moving bar stimuli and each group’s mean response. The correlation coefficient (−1…1) for each pair was then multiplied with the cell’s depth profile and a correlation-rank based weighting factor W = 0.9rank −1. Thus, each individual recording yielded a complete two-dimensional map, with IPL depth on one axis and functional group on the other. Next, we averaged across the maps for those cells that passed our response quality criterion (n = 31/51; n = 24/33; see above). The resultant matrices were normalized in two steps: First, we divided each group’s IPL depth profile by the mean depth profile of all included cells to eliminate any bias in sampling depth. Second, we divided each depth profile by its own maximum. This resulted in an automatic and unbiased estimate of dendritic stratification depth for all RGC groups (Fig. 5). Note that this automated approach is based on a relatively small sample of reconstructed cells and therefore can only provide an approximate prediction of stratification levels. This approach is invariant to differences in lateral dendritic field dimensions that may be associated with retinal position (for example, refs 10,23).