The temporal structure of the inner retina at a single glance

The retina decomposes visual stimuli into parallel channels that encode different features of the visual environment. Central to this computation is the synaptic processing in a dense layer of neuropil, the so-called inner plexiform layer (IPL). Here, different types of bipolar cells stratifying at distinct depths relay the excitatory feedforward drive from photoreceptors to amacrine and ganglion cells. Current experimental techniques for studying processing in the IPL do not allow imaging the entire IPL simultaneously in the intact tissue. Here, we extend a two-photon microscope with an electrically tunable lens allowing us to obtain optical vertical slices of the IPL, which provide a complete picture of the response diversity of bipolar cells at a “single glance”. The nature of these axial recordings additionally allowed us to isolate and investigate batch effects, i.e. inter-experimental variations resulting in systematic differences in response speed. As a proof of principle, we developed a simple model that disentangles biological from experimental causes of variability and allowed us to recover the characteristic gradient of response speeds across the IPL with higher precision than before. Our new framework will make it possible to study the computations performed in the central synaptic layer of the retina more efficiently.

. Overview of the two-photon (2 P) microscope equipped with an electrically tunable lens (ETL). For simplicity, most lenses and silver mirrors (M) were omitted. For parts, see Table 1. (A) Schematic diagram of the microscope's main optical paths, with the EL-16-40-TC (Optotune) inserted before the x-y galvo scan mirrors (M Scan ). Inset: Cross section and working principle of the ETL; a voice-coil actuator generates pressure on a container, which in turn pushes optic fluid into the lens volume sealed by polymer membrane and, thereby, modulating the curvature of the lens surface. (B) Photograph of the excitation path before the scan mirrors. CM, cold mirror; DM, dichroic mirror; BP, band pass filter; PMT, GaAsP photomultiplier tube; RL, relay lens; T, telescope. Panel A adapted from Euler et al. 22 , inset adapted from Optotune website (https://www.optotune.com/).

Methods
Equipping the 2 P microscope with an ETL. To allow for axial scanning, we modified a movable objective microscope (MOM, designed by W. Denk, now MPI Martinsried; purchased from Sutter Instruments/Science Products, Novato, CA). Design and configuration of the MOM have been described elsewhere 7,22,23,28 . In brief, the microscope is driven by a mode-locked Ti:Sapphire laser (MaiTai-HP DeepSee, Newport Spectra-Physics, Darmstadt, Germany) and equipped with two fluorescence detection channels and a 20x water immersion objective (W Plan-Apochromat 20×/1.0 DIC M27, Zeiss, Oberkochen, Germany). An ETL with an open aperture of 16 mm (EL-16-40-TC, Optotune, Dietikon, Switzerland) was introduced into the laser path before the scanning unit (Fig. 1a); the optical path after the scanners was left unchanged. Before the ETL, the laser beam is expanded to a diameter of 15 mm using a telescope (T1, a 4f-system with relay lenses RL1,2; f R1 = 30 mm, f R2 = 250 mm; for complete parts list, see Table 1). The expanded beam is then reflected perpendicularly to the optical table by a silver mirror towards the horizontally placed ETL, which is housed in a 60 mm cage plate (Thorlabs, Dachau, Germany) using a custom-made adapter (Fig. 1b). After the laser beam passed the ETL, it is narrowed by a second telescope (T2, another 4f-system consisting of relay lenses RL3,4; f R3 = 80 mm, f R4 = 30 mm) to a diameter of approx. 5.6 mm, which approx. matches the size of the scanning mirrors. Telescope T2 relays the pupil of the ETL to a conjugate pupil on the scan mirrors. Here, RL4 needs to be positioned precisely at its focal distance and in the centre of the two scan mirrors for correct refocusing of the laser beam. By changing the current supplied to the ETL, it changes its optical power. According to the EL-16-40-TC's specifications (see link in Table 1), its optical power can be tuned from −2 to +3 dioptres (f ETL ranging from −500 to 333 mm, for currents of approx. ±250 mA), rendering the beam divergent or convergent, respectively. Adapting the calculation in Fahrbach et al. 29 , the shift in focal plane (Δz, in [µm]) under the objective lens can be roughly estimated using = ′ , as the objective tip is immersed in solution. To drive the ETL with our imaging software (ScanM, see below and Table 1), we used custom electronics (designed by R. Berndt, HIH, Tübingen) that translates a voltage signal from one of the analogue-out channels of an PCI 6110 card (NI, Austin, US) controlled by ScanM into a stable current signal. To obtain the relationship between ETL driver input and resulting shift in focal plane, we applied voltage steps of varying amplitude (n = 11 amplitudes, n = 5 trials), presented in a randomized order, and monitored the shift in focal plane by reading out the z position of the microscope's motorized scan head. For the used combination of lenses, this resulted in a measured Δz range of +80 and −120 µm (for ETL driving currents of −100 and +100 mA; cf. Fig. 2B). Due to technical limitations (i.e. size of the scan mirrors), the Δz range with largely constant laser power spanned approx. 50 µm (cf. Fig. 2C), which is sufficient to scan across the entire mouse IPL without adapting the laser power. To characterize the spatial resolution of our system, we measured the point spread function (PSF) of fluorescent beads (170 nm in diameter, λ Emission = 515 nm; P7220; Invitrogen) at different axial planes. Animals and tissue preparation. All animal procedures were approved by the governmental review board (Regierungspräsidium Tübingen, Baden-Württemberg, Konrad-Adenauer-Str. 20, 72072 Tübingen, Germany) and performed according to the laws governing animal experimentation issued by the German Government. For all experiments, we used adult mice of either sex from the following lines: B6;129S6-Chat tm2(cre)Lowl J (n = 3; ChAT:Cre, JAX 006410), and B6;129P2-Pvalb tm1(cre)Arbr /J (n = 6; PV:Cre, JAX 008069). All lines were purchased from The Jackson Laboratory (Bar Harbor, ME). The transgenic lines were crossbred with the Cre-dependent red fluorescence reporter line Gt(ROSA)26Sor tm9(CAG-tdTomato)Hze (Ai9 tdTomato , JAX 007905) for all experiments. Owing to the exploratory nature of our study, we did not use blinding and did not perform a power analysis to predetermine sample size. For details on the mouse lines and all reagents, see Table 2.

Light stimulation.
A modified LightCrafter (DLPLCR4500, Texas instruments; modification by EKB Technology) was focused through the objective lens of the microscope 22,30 . Instead of standard RGB light-emitting diodes (LEDs), it was fitted with a green (576 nm) and a UV (390 nm) LED for matching the spectral sensitivity of mouse M-and S-opsins 31,32 . To prevent the LEDs from interfering with the fluorescence detection, the light from the projector was band-pass-filtered (ET dual band exciter, 380-407/562-589, AHF) and the LEDs were synchronised with the microscope's scan retrace 23 . Stimulator intensity was calibrated to range from 0.1 × 10 3 ("black" background image) to 20 × 10 3 ("white" full field) photoisomerisations P*/s/cone. The light stimulus was centred before every experiment, such that its centre corresponded to the centre of the recording field. Light stimuli were generated using QDSpy, a custom visual stimulation software written in Python 3 (see Table 1). To probe BC function, we presented After jumping back to the beginning of a frame, the ETL requires a few milliseconds to settle; this "settling" generates an artefact at the bottom of the frame and makes the film appear wider in frame 7 (for details, see Results). Inset: Frame 3 with intensity distribution along z-axis; for this scan configuration, the fluorescent band width was 4.4 pixels ± 0.1 (mean ± s.d. for width at half maximum, n = 5 measurements), corresponding to a pixel "height" of 1.1 µm. (E) Illustration of point spread function (PSF) measurements at three positions (−18 (red), 0 (black), 18 μm (blue)) along the z-axis (right); example images of fluorescent beads (170-nm beads, λ Em, Peak = 515 nm; 256 × 256 pixels, n = 60 z-planes, Δz = 0.2 µm, zoom XY = 8) at 0 μm, with mean Gaussian fits (n = 3 measurements/plane). PSF x and PSF z indicate the mean ± s.d. across the three axial planes (n = 9 measurements; see Table 3).
3-4 repeats of a "chirp" stimulus in two sizes, local (100 µm in diameter) and global (800 µm); it consisted of a light-On step followed by sinusoidal intensity modulations of increasing frequency and contrasts 7 .
Data analysis. Data preprocessing was performed in IGOR Pro 6 (Wavemetrics), DataJoint 33 and Python 3. Regions of Interest (ROIs) were defined automatically by custom correlation-based algorithms in IGOR Pro 7 . First, we estimated a correlation image by correlating the trace of every pixel with the trace of its eight neighbouring pixels and calculating the mean local correlation (ρ local ). In contrast to previous x-y recordings 7 , local correlation of neighbouring pixels varied with IPL depth in x-z scans (Fig. 4B) due to differences in iGluSnFR labelling ( Fig. 3B) and laser intensity (Fig. 2C). To account for that, an IPL depth-specific correlation threshold (ρ threshold ) was defined as the 70 th percentile of all local correlation values in each z-axis scan line. For every pixel  Table 3. Point spread functions (PSF) from n = 3 independent measurements at each z position.   www.nature.com/scientificreports www.nature.com/scientificreports/ with ρ local > ρ threshold ("seed pixel"), we grouped neighbouring pixels with ρ local > ρ threshold into one ROI. To match ROI sizes with the sizes of BC axon terminals, we restricted ROI diameters (estimated as effective diameter of area-equivalent circle) to range between 1 and 4 μm.
To register each ROI's depth in the IPL, we determined for each x-z scan the position of the On and the Off ChAT band ( Fig. 3; see also Results). In the case of ChAT:Cre x Ai9 mice, ChAT bands could be imaged directly due to their tdTomato fluorescence (Fig. 3B). We found the ChAT band positions can also be reliably determined from the IPL borders, which were determined from the location of iGluSnFR-labelled (or, in PV:Cre x Ai9, tdTomato-labelled) somata in GCL and INL. Here, we first determined the position of the IPL borders relative to tdTomato labelled ChAT bands in a subset of experiments (Fig. 3C). Following conventions 34 , we defined the On and the Off ChAT band position as 0 and 1, respectively. For every ROI, we then estimated the shortest distance to On and Off ChAT bands or IPL borders and expressed ROI IPL depth as a relative value between approx. −1 (GCL border) and +2 (INL border). Next, for every scan field (=batch), ROI depth estimates were corrected using the IPL depth at which the response polarity switched between On and Off BCs. The IPL position of this polarity switch was determined using the first principal component (PC) of the local chirp responses (0.24 ± 0.14, mean ± s.d., n = 37 fields) and subtracted from each scan field's depth estimates. We then added 0.5 to align the IPL centre, i.e. the separation between On and Off BC terminals to previous definitions 34 .
The glutamate traces for each ROI were extracted using the image analysis toolbox SARFIA for IGOR Pro 35 . Then, the traces were synchronised to the light stimuli using time markers that were generated by the stimulation software and acquired during imaging. Finally, we up-sampled the traces to 64 Hz temporal resolution and de-trended them by applying a high-pass Butterworth filter with a cut-off frequency of 0.1 Hz. www.nature.com/scientificreports www.nature.com/scientificreports/ Linear mixed effects model. All modelling was performed using DataJoint and Python 3. Batch effects were first estimated using a series of simple linear mixed effects models, which predicted the expected glutamate release of all ROIs across time (Y T Hz s , 64 32 2048 ) as a linear function of different predictors that were all encoded as dummy variables: (1) A model that used only the polarity (X polarity Nx2 ∈  ) to predict ˆ= Y X w polarity polarity . Thus, this model computed simply the average of all On and Off ROIs as the weight vector.
(2) A model that used only the IPL depth (X depth Nx20 ∈  ), where the response is estimated non-parametrically with 10 depth bins across the IPL for each polarity (cf. Fig. 5A), to predict Ŷ X w depth depth = .
(3) A model that estimated the response of each ROI as a function of the batch ( ∈  X batch Nx B 2 ; B 36 = , the first batch was left out as reference and to avoid a singular design matrix) where each batch had one predictor per polarity to predict ˆ= Y X w batch batch . (4) And finally, a combined model (2. and 3.) that predicted the response of each ROI as a function of both batch and depth (Ŷ X w X w batch batch depth depth In the last model, we wanted to make sure that any shared variance between the two predictors would be ascribed to the depth predictor -conservatively mistaking nuisance or noise from the batches rather as biological signal than vice versa. To achieve this, we orthogonalized X batch with respect to X depth , effectively by using X X P X _ . Then we define a temporal filter for this ROI as where we set w 0 i,0,2 = to represent the DC component with the sine, = T 64 the temporal extent of our filter reaching back 1 s in time (64 Hz), K 21 = the highest frequency of the Fourier basis for our kernels (i.e. well below the Nyquist, and reasonably smooth), and α i a temporal stretching factor. The sigmoidal factor before the sum is a soft-thresholding mask that sets the last part of our filter to 0 to avoid entering the next cycle when 1 i α > . The response is then predicted asˆβ Where β i is an offset of this ROI and * denotes convolution. We fit three models of decreasing flexibility: (1) A model with a separate kernel for each ROI, i.e. fixing α = 1 i and fitting w i k , and β i for all ROIs. This model is the most flexible, however, at the price of not yielding any abstract insight into our data.
(2) A model with one learned shared kernel, i.e.
, and one temporal stretch per ROI, i.e. fitting α i for all ROIs. Additionally, each ROI learned a scalar a f i i to scale and flip (for On and Off BCs) the shared kernel.
(3) Finally, a model with a shared kernel (like model 2.) and a stretch that is a function of the ROI's depth d i and Firstly, using effectively only the depth and fitting a weight j ξ for each depth bin c j (cf. Fig. 5A) to model the speed of each ROI: . Secondly, the same model but with an additive shift b i ψ for each batch: . And finally, with a nonparametric interaction between batch and depth, i.e. separate depth bin weights b j

Estimating explained variance.
For an observed y and a predicted ŷ response, we estimated the explained variance asÊ Estimating BC synapse number in x-z slices. We estimated the average number of synapses (~ROIs) per BC type to be expected in an x-z scan, which can be considered as a ~0.5 µm-thick optical section (cf. PSF measurements, Table 3). To this end, we used a published EM dataset (e2006, ref. 10 36 . Finally, we estimated the number of ribbons (synapses) per slice (x-z scan) and BC type by dividing each BC's axon terminal volume in a slice by its average axon terminal volume per ribbon (Fig. 4F).

Results
Setting up axial scanning. To allow axial scanning of the retina, we inserted an ETL into the optical pathway of the 2 P laser of our microscope (Fig. 1). By electrically modulating the optical power of the ETL, the beam of the 2 P laser entering the microscope´s objective converges or diverges, resulting in a focus shift along the z-axis of the recording. For the ETL used here, a change in electrical current is transformed into a pressure change, which in turn regulates the lens volume and, thereby, the curvature of the lens surface (Fig. 1A, inset).
When the lens is positioned vertically (w/optical axis parallel to the table), its fluidic core may be slightly deformed by gravitational forces, resulting in a deterioration of its optical properties. The simplest horizontal arrangement would be to place the ETL directly above the objective 25,37 . However, we decided against this possibility for two reasons: First, in this position, the ETL would introduce a focal plane-dependent change in image magnification 25 . Second, if the visual stimulus is coupled into the laser pathway after the scan mirrors -like in our www.nature.com/scientificreports www.nature.com/scientificreports/ setup (Fig. 1A) 22,30 -the ETL would also modulate the stimulation plane and spectrally filter the stimulus. The latter is critical if UV stimuli are used: Since prolonged exposure to UV light can degrade the optomechanical properties of the lens, ETLs are typically equipped with a UV-reflecting glass window (for specifications, see Table 1). Therefore, placing the ETL into the stimulus path hampers UV stimulation, but UV stimulation is required to properly drive the mouse retina with its UV-sensitive cone photoreceptors (λ Peak = 360 nm; as discussed in ref. 30 ). To avoid these issues, we positioned the ETL horizontally in the pathway that reflects the excitation laser onto the scan mirrors (Fig. 1A,B) (F. Voigt, F. Helmchen, personal communication).
To make use of the full aperture of the ETL, we pre-expanded the laser beam using a telescope system (T1). After the ETL, we used a second telescope system (T2) to refocus the laser beam onto the scanning mirrors. As T2 determines the beam diameter that enters the back aperture of the objective and the z-range of focus shift (together with the objective lens' magnification), it needs to be defined for the specific microscope and the desired z-range.
First, we evaluated the ETL's performance using a solution containing the red fluorescent dye SR101 ( Fig. 2A). Our combination of custom ETL driver and optical configuration allowed for a practical focus shifting range of approx. 200 µm (Methods). Importantly, the same voltage signal reliably resulted in the same z position of the focal plane (s.d. <1 µm) and when staying within a limited voltage range (e.g. ±0.4 V), the relationship between voltage input and z position was sufficiently linear (Fig. 2B). However, we consistently observed intensity fluctuations in the first few lines of each x-z scan frame (Fig. 2C,D). This "artefact" can be explained by the fact that www.nature.com/scientificreports www.nature.com/scientificreports/ for large, rapid changes in voltage input -as it happens, for instance, when jumping back from the last line of a frame to the beginning (retrace) -the ETL´s new refractive state requires a few milliseconds to settle 25,29 . These brief oscillations are visible at the bottom of the image (where the frame starts; see particularly Fig. 2D-7). We dealt with this problem simply by excluding the first couple of lines (~10 ms) of each frame from our analysis. In addition, we found that the fluorescence smoothly varied along the z axis (Fig. 2C, left scan). The decrease in fluorescence reflects a loss in laser power, which happens when the ETL changes the beam's collimation to an extent that the beam becomes too large for the scan mirrors or is partially blocked by down-stream optics. This can be improved by applying an offset voltage to the ETL driver signal, such that the laser intensity peak covers the required z scan range. As the intensity peak is relatively shallow, the IPL of the mouse (~40 µm) can be imaged with almost constant laser intensity.
Next, we tested whether pixel size and spatial resolution remained constant along the z axis of an x-z scan suitable for capturing the complete mouse IPL (e.g. 64 × 40 pixels corresponding to 50 × 40 µm). By moving the microscope's objective lens, we placed a thin fluorescent film (Methods) at different z positions within an x-z scan field (Fig. 2D) and measured the film's thickness. Apart from the artefact (see above), the recorded thickness of the film remained constant, suggesting that the pixel size is constant along the z axis. Finally, we quantified the spatial resolution of our system by measuring the point spread function (PSF) of fluorescent beads both in the x-y plane and at different z positions (Methods). This was done by first setting one of three z planes (Fig. 2E) using the ETL and then taking image stacks (using the microscope's motorized stage). The measured PSFs were around 0.4 and 1.8 µm along the x and z axis, respectively, and varied very little for the different ETL planes (Table 3).
In summary, our ETL configuration allows for spatially (nearly) linear, fast axial imaging without detectable loss in spatial resolution.
Axial scanning in the IPL of the mouse retina. Axial scans were evaluated by imaging light-evoked glutamate release from BC axon terminals (Fig. 3A). After AAV transduction, the glutamate biosensor iGluS-nFR 26 was ubiquitously expressed across the whole retina, including the IPL (Fig. 3B) 7 . As the axon terminals of different BC types stratify at distinct levels of the IPL 10,13,38 , registering IPL depth within the x-z scans is critical. Important landmarks in the IPL are the so-called ChAT bands, which are formed by the dendritic plexi of the cholinergic starburst amacrine cells (SACs) 39 . Accordingly, a commonly used metrics for IPL depth is to define the inner ("On") band as the origin (=0) and the distance to the outer ("Off ") band as 1 (Fig. 3C) (see also 7,34 ). To relate these positions to IPL borders, we used transgenic mice in which SACs were fluorescently labelled ( Fig. 3B; for details, see Methods). We found that the relative distance between ChAT bands and IPL borders was highly consistent across scans and mice (Fig. 3C). Hence, for mice lacking fluorescently labelled ChAT bands, IPL depth can be reliably estimated from the IPL borders.
The ubiquitous expression of iGluSnFR combined with axial scanning allowed sampling of glutamate release at all IPL depths (Fig. 4). To achieve scan rates >10 Hz, we used x-z scans with 64 × 56 pixels (1.6 ms/line) at a zoom factor that yielded a pixel size of ~1 µm (Fig. 4A). Regions of interest (ROIs) were based on local image correlations with an IPL depth-dependent threshold (Methods). This ensured ROI placement across the entire IPL (Fig. 4B,C). Subsequently, ROIs were quality-filtered using the reliability of their responses 7 .
To evaluate if the anatomical and functional properties of single ROIs are consistent across x-y and x-z scans, we compared the distribution of ROI sizes and response quality indices ( Supplementary Fig. S1). We found that ROI sizes were only slightly larger in x-z than x-y scans, which is likely due to the lower resolution in the z-compared to y-dimension (cf. Fig. 2E), but were still within the expected range of BC axon terminals (2-8 µm 2 , cf. Extended Data Fig. 1 in ref. 7 ). In addition, response quality was overall comparable between ROIs recorded in x-z and x-y scans. The higher number of low-quality ROIs in x-z scans (quality index <0.3) originates from ROIs located at the IPL borders, a region that was only rarely sampled in x-y scans. Together, this suggests that ROIs estimated from x-z scans are qualitatively comparable to ROIs estimated from x-y scans.
Recently, mouse BCs have been functionally characterized by recording BC glutamate release in x-y scans in response to local and global chirp stimuli activating the centre and both centre and surround of the cells´ receptive fields, respectively 7 . Here, we used the same stimuli to record glutamate traces from BC axon terminals in axial IPL x-z scans (Fig. 4D). As expected from the subdivision of the IPL into On and Off sublamina 19,40 , ROIs located in the outer half displayed increases in glutamate release upon light decrements (Off cells; see ROIs 1-3 in Fig. 4C,D). In contrast, ROIs in the inner part of the IPL showed increased glutamate release upon light increments (On cells; ROIs 4-6). In addition, response transience varied with IPL depth, with more transient (e.g. ROIs 3 and 4) and sustained BCs (e.g. ROIs 1 and 6) stratifying in the IPL centre and borders, respectively. This observation is in agreement with previous studies demonstrating a spatial segregation of the IPL in transient and sustained BC channels 15,20,21 .
Next, we wanted to evaluate if a single x-z scan -in principle -could capture the full complement of BC types. Overall, ROIs with chirp responses that met our quality criterion 7 were almost evenly distributed across the IPL (Fig. 4E), suggesting that generally we did not under-or oversample certain IPL levels. However, it has been long known that BC types come at different densities and vary in the size of their axon terminal system 9 . For example, S-cone selective type 9 BCs 11,41 are low-density and have sparse terminals, whereas rod BCs have small terminals but make up ~38% of all BCs in the mouse 9,42 . To estimate the number of ROIs (~output synapses) the different BC types are expected to contribute to an x-z scan (as in Fig. 4), we calculated the average volume each BC type's axon terminals occupy in an 0.5 µm (~PSF X ) thick slice of the IPL using electron microscopy data 10 and then determined the mean number of synapses in this volume using published ribbon densities 36 (Methods; Fig. 4F). This analysis showed that some BC types (i.e. X, 8 and 9) are expected to contribute only few synaptic terminals to an x-z scan.
Taken together, axial x-z scans allow recording glutamate signals virtually simultaneously across the IPL with similarly high fidelity compared to sequential, hence time-consuming, "traditional" x-y scans. While the (2020) 10:4399 | https://doi.org/10.1038/s41598-020-60214-z www.nature.com/scientificreports www.nature.com/scientificreports/ functional diversity that can be recovered from individual x-z scans qualitatively resembles that described in an earlier study 7 , reliably capturing signals from BC types with low terminal densities requires integrating data from multiple scans.

Identification of batch effects. We recorded local chirp light-evoked BC glutamate release from 5,379
ROIs (37 scan fields, 6 mice) across the entire IPL (cf. Fig. 4E). Of those, 3,893 ROIs passed our quality criterion as previously defined in 7 and were selected for further analysis. When visually inspecting the data obtained from different recordings, we noticed that the timing of recorded glutamate traces varied systematically across recordings (Fig. 6A). We refer to this variation as "batch effects", in accordance with similar inter-experimental variability in the single-cell genetics literature (e.g. refs. [43][44][45][46]. In our experiments, the variations between scans may have been caused by experimental factors such as slight temperature fluctuations, as well as differences in light adaptation and/or fluorescence biosensor expression. To investigate the variability in our data, we performed principal component analysis (PCA) on the recorded time-series and inspected the projection onto the first two principal components (PCs; Fig. 6B,C). While On and Off BCs could be distinguished clearly based on the first PC (Fig. 6B), a substantial portion of the variability observed within On and Off BCs seemed to stem from variability across recordings (Fig. 6C). Qualitatively, these batch effects were large enough to be a challenge for recovering the biological differences between cell types within the On and the Off BCs.
We quantified the relative contribution of three factors to the total variance of the observed signal: (1) polarity, i.e. whether the ROI was located in the On or Off sublamina, (2) IPL depth bin in which the ROI was recorded, and (3) the batch (scan field) from which the ROI originated. To this end, we fit a series of linear models (Fig. 5A), each of which included one or more of the three factors (On/Off, IPL depth, batch), and estimated the fraction of variance explained by the models (Fig. 5B,C). The first model, which captured only the polarity, accounted for 34% of the response variance. The second model, which used only IPL depth as a predictor, accounted for 48% of the variance. Note that the first model is a special case of the second one, obtained by splitting the IPL into On and Off sublaminae. As a third model, we used polarity × batch (scan field) ID as a predictor. This model, which amounts to estimating the average trace of On and Off cells in each scan field, accounted for 66% of the variance, substantially outperforming the previous models that only considered the biological source of variation. Finally, adding IPL depth bin as a predictor improved the explained variance only marginally (to 69%).
To summarize, we found that batch effects alone accounted for a larger fraction of the variance than IPL depth (Fig. 5B), which suggests that accounting for such variation can greatly facilitate any analysis of functional differences between BC types beyond On vs. Off.
What is the nature of these batch effects? The most salient difference across the three example batches was a shift in response speed (Fig. 6A). This is especially striking in the response to the chirp's frequency modulation, where the batch-averaged responses are almost entirely out of phase (Fig. 6A, bottom right). We found the same temporal misalignment in the predictions of our model that considered only IPL depth but ignored the batch effects (Fig. 5C). Comparing the predicted and the recorded traces, we observed that the model was too fast for the first (slow) batch, approximately aligned for the second (medium) batch and too slow for the third (fast) batch. This observation is in line with a previous study that reported systematic differences in the response speed of RGCs recorded from different macaque retinae 47 .
A possible explanation for shifts in response speed between experiments may be differences in recording temperature. While we used a closed-loop system to keep the temperature of the tissue at 36 °C, small temperature fluctuations in the order of ±1 °C cannot be excluded. The temperature coefficient (Q 10 ) of biological reactions, including neural activity, is typically between 2 and 4 (cf. ref. 48 ), therefore, following = ∆ ∆ Q Q T T 10 /10 , a temperature increase of 1 °C may result in a 7 to 15% increase in response speed, which is in the range that we estimated for the batch effects (cf. Fig. 7D). It is well known that temperature affects neural processing. For instance, swordfish actively raise their retinal temperature by 10 to 15 °C, thereby increasing temporal resolution up to ten times to gain a predatory advantage 49 .

Building batch and IPL depth variations into a shared BC encoding model. To investigate the
idea that batch effects effectively result in changes in response kinetics more directly, we fit a linear encoding model and estimated the temporal receptive field kernels of the ROIs in the three example scans shown before. As expected, the temporal kernels showed systematic differences between the three scans that seem to be largely explained by rescaling them in time (Fig. 8A,B). Moreover, within a single batch we could still discern the underlying IPL gradient: ROIs closer to the IPL centre (=lighter colours in Fig. 8C) had their leading edge closer to zero and, hence, responded faster. In addition, central ROIs displayed more biphasic kernels and, hence, responded more transiently.
The data presented above suggest that functional differences between individual ROIs can, to a large extent, be accounted for by modelling response speed, and this speed depends on two main factors: (1) laminar location within the IPL and (2) batch effects due to variability between scan fields. We therefore developed a very simple joint encoding model that reduces the functional differences between neurons (here ROIs) to a temporal rescaling of their response kernels (Fig. 7A). The model learned exactly one response kernel that is shared among all neurons and all scan fields. In addition, it allowed for a temporal rescaling of this kernel and a neuron-specific scaling of the response magnitude and polarity-flip (On/Off).
First, we observed that sharing the same kernel shape across all ROIs and only adjusting speed, scale and polarity yielded a predictive performance of 49% explained variance compared to 56% by the model with an individual temporal kernel for each ROI (Fig. 7B). While a 7% difference is not negligible, the difference in complexity between the models is considerable: The latter model used more than 40 parameters for each ROI to specify a response kernel, whereas the former (simplified) model had effectively just 3 parameters to model speed, Scientific RepoRtS | (2020) 10:4399 | https://doi.org/10.1038/s41598-020-60214-z www.nature.com/scientificreports www.nature.com/scientificreports/ scale and polarity for each ROI. Moreover, the simplified model allows us to assess how well we can predict the response speed of individual ROIs based on their IPL depth and a batch-specific speed adjustment. The complex model -while more accurate -does not allow us to make such abstraction, because it models each ROI separately www.nature.com/scientificreports www.nature.com/scientificreports/ and does not include any simplifying assumption of shared variability between ROIs. For us, it serves as a benchmark to judge the interneural variability in the shape of the response kernel, which the simpler model cannot capture by design. Therefore, for the remainder of the paper, we focussed on the simplified model.
Next, we tested the effect of additionally constraining this simple model: First, we assumed that the speed of each ROI is only a function of the ROI's depth in the IPL. This assumption, which meant that all ROIs with the same laminar location share the same response kernel, decreased the predictive performance to 40% explained variance (Fig. 7B). Alternatively, allowing for batch effects by adding a scan field-specific global shift to the speed estimates for all the ROIs in the same scan field, enabled us to capture 46% explained variance. This was similar to the model that allowed each ROI its own speed (49%). Additionally, including an interaction term between IPL depth and batch improved performances only slightly (47%), suggesting that batch variations had similar, approximately additive effects onto the speed across the IPL.
In summary, to the extent that variability between BCs can be modelled by differences in response speed, these speed differences can almost entirely be accounted for (46% vs. 49% explained variance) by laminar location within the IPL and a batch-specific global shift in response speed.
Moreover, with this model we can give a quantitative estimate of the speed gradient across the IPL (Fig. 7C). In line with earlier reports 7,15,20,21 , the BC response speed fell off from the centre of the IPL towards its borders. Notably, between the ChAT bands this speed gradient was nearly symmetrical between On and Off BCs. However, after that the On BCs levelled off and exhibited the same, relatively slow speed until the ganglion cell layer. The Off BCs, by contrast, continued decreasing in their speed almost all the way until the inner nuclear layer.

Discussion
We implemented fast axial x-z scanning of the whole-mounted mouse retina by equipping a 2 P microscope 24 with an ETL to rapidly shift the focal plane of the laser. We showed that this experimental setup is suitable to record the light-evoked glutamatergic output of BCs almost simultaneously across the complete IPL. Axial scans enabled comparing temporal response properties between IPL strata more directly than "traditional", time-consuming series of horizontal (x-y) scans. At the same time, x-z scans allowed identifying batch effects, characterized mostly by inter-experimental differences in signal speed. We showed that already batch correction with a simple linear model can improve recovering the characteristic response speed profile across the IPL. Our results indicate that careful consideration of inter-experimental variance is key for extracting functional differences between neurons.
Techniques for axial scanning. Several technical solutions that enable fast axial scanning in 2 P microscopy have been published. Here, two main approaches can be distinguished: In the first group, a focus change within the sample is realized by moving the objective lens relative to the sample (or vice versa). This has, for instance, been implemented using a piezo to move the objective along the z axis 50 . By coordinating the trajectories of galvo scanners (x-y) and piezo (z), fast volume scanning (i.e. 10 Hz frame rate for a 0.25 mm cube) can be achieved using 3D spiral patterns. While this solution is relatively easy to implement, the inertia of the objective limits the speed with which a focal plane can be selected. In addition, the objective's movements may introduce vibrations to the sample.
In the second group, focus shifting is achieved by changing the collimation of the laser beam ("remote focussing") while the objective-to-sample distance remains constant. This strategy eliminates movements close to the sample. A classical solution for remote focussing is to add a reference objective to the laser path to axially displace the focal plane in the sample 51 . Because only a lightweight mirror under this reference objective is moved, this arrangement can reach high axial velocities. However, its superb optical performance is offset by high complexity (i.e. optical alignment; ref. 52 ) and the costs of a second high-quality objective. In a different approach, an arrangement of inertia-free acousto-optical deflectors (AODs) replaces the galvo scanner and, thus, allows for very fast random-access scanning 25,[53][54][55][56] . Using a clever AOD arrangement, random-access in 3D is possible 54,57 . Due to the absence of mechanical parts, such a solution enables extremely fast focus changes (<1 ms) at large z ranges (>1 mm). However, AOD-based solutions typically are complex and expensive systems, which usually require substantially more laser power than galvo scanner-based systems.
Comparably fast remote focussing can also be achieved with an ETL, in which the curvature -and hence the focal length -of a liquid lens core is controlled electrically (reviewed in ref. 58 ). With high-quality www.nature.com/scientificreports www.nature.com/scientificreports/ inexpensive ETLs becoming available, they offer a cost-efficient and relatively simple way for equipping mechanical scanner-based (or "simple" 2D AOD-based) fluorescent microscopes, including confocal 59 , 2P 25,37,60,61 , and light-sheet systems 29,62 , with fast focussing. In the simplest configuration, the ETL is positioned directly on top of the objective lens 25,37 . However, in this position, shifts in focal plane are accompanied by image scaling in x-y 25 . Also, the ETL's transmission in the relevant spectral bands -like here, the UV transmission for light stimulation (cf. Results) -may need to be considered. By integrating the ETL into the laser path upstream of the scan mirrors 59 , image scaling and (potential) ETL transmission issues are avoided. The additional telescope needed to couple the ETL into the laser path slightly increases complexity but at the same time allows adapting the available z range to the experimental needs.
In the current study, we applied ETL-enabled axial scanning to the isolated but intact whole-mounted mouse retina. We found axial scans to work across the whole mouse retina 63 -except for areas very close to the optic disc (where the nerve fibre layer becomes very dense) or close to the edge (damage from excising the retina). Our approach should in principle be applicable to a wide range of experiments, including in vivo recordings in the brain. As with other optical techniques, a main limitation of scanning depth is scattering within the tissue.
Further improving ETL-based axial scanning in the ex-vivo retina. For x-z scans across the mouse IPL, we used a unidirectional scan mode, where at the end of a frame, the focal position of the excitation laser is shifted back to the first scan line in one ~50 µm "jump". This results in the aforementioned artefact in the first few lines of each frame and is caused by fast oscillations in the ETL's focal power after a rapid change in driving current (see link to specifications, Table 1). The stabilization time of ~10 ms (for travel distances of ~50 µm) we observed is consistent with earlier reports 25,64 . For simplicity, we here used scans with more lines, such that the artefact was outside the IPL. Alternatively, optimizing the current trajectory driving the ETL -e.g. by using a steep ramp and "overdriving" the current instead of just a step (cf. ref. 25 ) -may dampen the ETL's oscillations and, thus, decrease settling time. Also, bidirectional scans that do not require large and rapid changes in z position may reduce such z travel distance-dependent artefacts.
Another potential caveat of an ETL is thermal drift, because driving the lens may slowly heat it up. Since the resistance of the coil that shapes the lens' core is temperature-dependent, also the current-to-focal power relationship depends on the ETL's temperature (for details, see ETL specs). We embedded our ETL into a solid adapter ring made from aluminium, which seemed to have kept the ETL's temperature stable enough, as we did not observe any relevant thermal drifts during the course of a recording. In any case, the ETL model we used features a build-in temperature sensor that can be read out via an I²C (Inter-Integrated Circuit) bus connection to monitor the ETL's temperature.
One consequence of the ETL being positioned upstream of the scan mirrors is a focus dependent change in laser power. Because we needed a relatively small z focus range to scan across the mouse IPL, we applied an offset voltage to the ETL, shifting the shallow peak in laser power to the imaging range. For larger z ranges, one could automatically adjust the laser power with a sufficiently fast modulator (i.e. a Pockels cell) as a function of the ETL's control signal.
Identification and removal of batch effects in axial scans of the mouse IPL. Batch effects -referring to inter-experimental variability -recently became a prominent topic in the single-cell genetics community (e.g. refs. [43][44][45][46], which motivated us to look for such effects in our data. Note that the batch effects we observed in the presented data are not specific to axial scans (or optical recordings). Instead, batch effects reflect experimental variability that can result from small differences in recording conditions (i.e. temperature; see below) but also method-specific variations: For imaging, this could be labelling intensity; for electrical single-cell recordings, electrode tip geometry.
For the functional characterization of retinal cell types, our previous studies used data obtained from sequential x-y recordings 7,28 . In the IPL, one disadvantage of this approach is that the sample of BC types in each individual scan greatly varies between scans, and the cells recorded within any one scan will typically share similar response properties. In the middle of the IPL, a third of all BC types can theoretically be recorded in an individual x-y plane (see stratification profiles in refs. 10,13,38 ), but typical scans will mostly sample 2 or 3 BC types. In contrast, axial x-z scans established here allow less biased recordings of the glutamate output across the entire IPL, with BC types with very different response properties present in the same scan field. As a result, each scan exhibited a highly stereotypic functional organization of response kinetics across the IPL, as described before 7,15,20,21 . This property greatly facilitated comparison of data obtained by different recordings, which allowed us to detect, quantify and correct for batch-specific variability in BC responses (see below).
In our recent studies characterizing BC and RGC types 7,28 , we did not explicitly correct for batch effects. However, we used other measures to minimize the effect of such inter-experimental variations on clustering. For BC recordings, we estimated a prior probability for cluster allocation for each scan field taken at a specific IPL depth, which was based on the relative axon terminal volume of all BC types at the respective depth (cf. Figure 2c in ref. 7 ). This is similar to clustering separately within bins of IPL depth, which helps to identify -on average -BC type-specific response signatures. The fact that ROIs from single scan fields were consistently assigned to different functional types suggests that type-specific differences could be resolved despite the variability induced by batch effects. For GCL recordings, we judged response quality based on alpha RGCs 65 , which are easily recognized by their large somata. Only if alpha RGCs displayed their characteristic temporal response profile 65,66 , we included the data. By doing so, we implicitly minimized the variability induced by batch effects. Again, we found that most functional groups were present across experiments, suggesting that experiment-specific speed differences did not induce additional functional clusters. Thus, if a readily identified and well-calibrated reference exists, batch effects can be reduced by excluding recordings that deviate strongly. However, this approach can only reduce batch (2020) 10:4399 | https://doi.org/10.1038/s41598-020-60214-z www.nature.com/scientificreports www.nature.com/scientificreports/ effects at the cost of experimental yield. It is therefore not making optimal use of the available data and resolving more subtle differences between cell types will be difficult as they are likely to be masked by batch effects.
Recent evidence indicates that batch effects or inter-retina variability are also present in electrophysiological recordings of primate RGCs 47 , where the response profile of well-characterized parasol cells was used to normalize (by temporal rescaling) and then combine data from different experiments. This adds further support to the hypothesis of temporal rescaling, but it also highlights the presence of batch effects in other experimental modalities. In addition, subtler batch effects, such as changes in receptive field size, spiking nonlinearity or autocorrelation of cell response have been reported 67 . In our data, we observed that temporal effects made the largest contribution to the batch variability. Other contributions could come from the glutamate-sensing kinetics of iGluSnFR, which depending on its expression and distribution can alter the postsynaptic response profile 68 . In this study, we focused on modelling explicitly temporal rescaling and found that this recovered most of the predictive performance of much more flexible models.