Stimulus display was performed with a conventional DLP projector using custom optics to reduce and focus the image onto the photoreceptor layer, with size of a pixel of \(\approx\)\(4 \mu m\), maintaining an average irradiance of \(70nW/mm^2\). Stimuli were generated using the code available at https://github.com/NeuralEnsemble/MotionClouds. Timing of the images was controlled via a custom-built software based on Psychtoolbox for MATLAB \citep{Brainard:1997we}. Spike sorting analysis was performed using Spyking-Circus \citep{Yger:2016bi} and all data fitting procedures were performed with LMFIT \citep{newville_2014_11813}. Data was analyzed via redistributable Jupyter notebooks using SciPy \citep{Jones01} and running Python 3 kernels. Figures were generated with Matplotlib \citep{Hunter:2007ih}.

Characterization of spatiotemporal tuning of RGCs

To produce a standard characterization of RGC for different stimulations, we measured the response to white noise checkerboard pattern (block size = \(0.05mm\)) for \(1200\sec\) at  \(60\) Hz and sinusoidal drifting gratings at maximum contrast (minimum and maximum Weber contrast of \(\approx-94\) and \(\approx 129\) respectively), with spatial frequencies and speeds detailed in in sequences of three seconds with 10 repetitions of each. Due to the constraints in recording length stemming from tissue viability, we focused only on changes in response to the speed of the moving stimulus and not on direction, so the present protocol uses a single direction for all sets of stimuli.
Receptive fields were estimated from the response to the checkerboard stimulus by reverse correlation, yielding the Spike-Triggered Average (STA)  \citep{Chichilnisky:2001ua} as a three-dimensional spatiotemporal impulse response \FIGreportE. The spatial characterization was performed by fitting a bi-dimensional Gaussian function at the time point of maximum amplitude, then drawing an ellipse at one standard deviation. The size of each RF was calculated as the radius of the circle with the same area as the ellipse fit \citep{Petrusca:2007je}. The shape was quantified by the eccentricity \(\epsilon\) of the ellipse as \(\epsilon=\sqrt{1-\left(b/a\right)^{2}}\) where \(a\) and \(b\) are respectively the radius of the major and minor axis. The temporal profile was computed as the time course of the intensity at the point with the largest variance, and then fitted to a difference of two cascades of low-pass filters \citep{Chichilnisky:2002wu}. Basal activity was set to zero and amplitude normalized, such that it follows the form
\begin{equation} \label{eq:temp_fit} \label{eq:temp_fit}R(t)=p_{1}\left(\frac{1}{\tau_{1}}\right)^{n}e^{-n\frac{t}{\tau_{1}-1}}-p_{2}\left(\frac{1}{\tau_{2}}\right)^{n}e^{-n\frac{t}{\tau_{2}-1}}\par \\ \end{equation}
where \(t\) represents time in number of image frames before the spike, \(\tau_{1}\) and \(\tau_{2}\) describe the temporal decay of the response of each filter; scalars \(p_{1}\) and \(p_{2}\) are amplitude responses of each filter, and \(n\) constitutes a free parameter. Tuning to motion was evaluated from the response to the drifting gratings; the response to each set of speeds and spatial frequencies was measured by Peristimulus Time Histogram (PSTH) for \(10\) repetitions of each sequence.
Tuning to speed \(v\) was fitted to a skewed Gaussian in which the response to speed \(R(v)\) takes the form
\begin{equation} \label{eq:v_fit} \label{eq:v_fit}R(v)=A\left[\exp{\left({\frac{-(\log_{2}v-\log_{2}V)^{2}}{2\times(\sigma_{v}+\zeta\times(\log_{2}v-\log_{2}V))^{2}}}\right)}-\exp{\left({\frac{-1}{\zeta^{2}}}\right)}\right]\par \\ \end{equation}
where \(V\) is preferred speed, \(\sigma_{v}\) curve width, \(\zeta\) represents the skew of the curve and \(A\) a scaling factor \citep{Priebe:2006jk}. Since \(\sigma_{v}\) is highly correlated with tuning bandwidth (width of the function at half-height) \citep{Pinto:2010dt}, we used \(\sigma_{v}\) for all bandwidth related functions. Changes in tuning bandwidth were measured as the normalized difference in \(\sigma_{v}\) as
\begin{equation} \label{eq:delta_sigma} \label{eq:delta_sigma}\Delta\sigma_{v}=\frac{\sigma_{v1}-\sigma_{v2}}{\sigma_{v1}+\sigma_{v2}}\par \\ \end{equation}

Natural-like stimuli

To measure the response of the retina to motion with different levels of naturalness, in terms of spatio-temporal correlations and spatio-temporal spectrum, we used synthetic random textures (details in \BOXmotionclouds) parametrized to obtain the same contrast, average spatial frequencies and speeds as the drifting grating stimuli.