Stochastic single cell migration leads to robust horizontal cell layer formation in the vertebrate retina

During development, tissues and organs change their morphology to acquire the characteristic sizes, shapes and cellular patterns that ensure correct function. These processes are tightly regulated to yield reproducible outcomes and precise tissue patterning. The vertebrate retina is a fascinating example of such pattern formation. As an extension of the central nervous system, it develops from a pool of multipotent progenitors into a neural tissue composed of five neuronal cell types. Neurons of the same type, as well as the one glial cell type, Müller glia, are positioned in distinct layers that give the mature retina its laminated appearance (Fig. 1A). Interestingly, retinal development contains stochastic elements. For example, stochasticity in the apical-basal movement of neuroepithelial cell nuclei has been proposed to contribute to proliferative versus neurogenic division outcome in zebrafish (Baye and Link, 2007; Del Bene et al., 2008). In addition, stochasticity has been suggested to be involved in retinal neuronal fate determination (He et al., 2012). Studies so far have concentrated on exploring stochasticity of fate decisions; however, to what extent other phenomena, particularly retinal lamination, also show variability or stochastic elements remains unknown. To address this issue, we explored single cell behaviour of horizontal cells (HCs) before lamination. HCs are inhibitory interneurons that localize beneath the photoreceptor (PR) cell layer and form synaptic connections with PRs and bipolar cells (BCs) (Amini et al., 2018). BCs ultimately transmit visual information to the retinal ganglion cells (RGCs) (Fig. 1A). Committed HC precursors (HCpr) divide once more before integrating into their laminar array, producing two HCs (Godinho et al., 2007; Weber et al., 2014). The timing and dynamics of this last cell cycle have not yet been studied. Furthermore, HC migration patterns are not fully understood. Unlike RGCs, which keep an apical connection and reach their final destination via stereotypic somal translocation (Icha et al., 2016), HCs are not attached but instead undergo multipolar bidirectional migration to reach their final position (Chow et al., 2015; Edqvist and Hallböök, 2004). Although different phases of movement have been characterized (Chow et al., 2015), the stereotypicity of these movements remains unexplored.

Here, we use long-term in vivo light sheet imaging to quantify single cell behaviour over a large cohort of HCs following cell-cycle parameters and bidirectional migration (Chow et al., 2015; Godinho et al., 2007; Weber et al., 2014). Surprisingly, we find that cell cycle parameters substantially vary among single cells and that migration patterns of single cells are non-predictable before robust final lamination. Thus, we provide the first evidence that, as seen for fate decisions, retinal lamination also contains variable and stochastic elements.

Before final lamination, committed HCprs undergo one more division. The position of this division is not restricted to the final HC location but can occur along the inner nuclear layer (INL) (Godinho et al., 2007; Weber et al., 2014). As previous studies used immunostaining of fixed samples, no appreciation of HCpr division dynamics existed so far.

Following 428 cells from 21 embryos (i.e. n=428, N=21) (for an experimental example see Movie 1 and Fig. 1B), we found that, as postulated previously (Godinho et al., 2007; Weber et al., 2014), all HCprs underwent an additional division resulting in two HCs that later entered the HC layer (data from all pooled movies, Table S3). None of the tracked HCprs underwent more than one terminal division (data from all pooled movies Table S3). Single cell analysis of mitotic location (Fig. 1D,D′) (n=428, N=21) revealed variations in HCpr mitotic position. Although most cells entered and underwent mitosis near the outer plexiform layer (OPL) where they later resided, divisions were also seen along the whole apico-basal length of the INL, even close to the inner plexiform layer (IPL) (Fig. 1C-E). All markers used (see Material and Methods) showed a similar mitotic positioning pattern (Fig. S1E). Only in very rare cases did HCprs divide in the RGC layer (Fig. S2A). Frequency analysis confirmed a bias towards OPL positions as 41% of HCpr divisions occurred close to the OPL (distance=0-6 µm), whereas the remaining 59% spread along the INL (Fig. 1D′,E,H). This suggests that HCprs favour divisions closer to the prospective HC layer. Although this bias has not been shown previously, it was postulated that divisions close to the OPL facilitate rapid synaptic contacts between HCs and their partners (Godinho et al., 2007).

To understand whether these positional differences were connected to differences in HCpr mitotic behaviour, we analysed the duration of mitosis from rounding to anaphase onset (Fig. 1C,F). Single cell and frequency distribution analysis (n=404, N=20) revealed that while 62% of HCprs divided within the first 30 min after rounding, mitosis length varied from 10 min to 260 min (Fig. 1F,I). This is markedly different from what has been reported for retinal progenitor cells that have very reproducible mitotic timing lasting 25±3 min (Leung et al., 2011). To reveal potential differences between cells with short versus long mitotic length, we used a nuclear marker to differentiate between chromatin condensation, metaphase plate formation and DNA segregation. We found that the time between metaphase and anaphase was similarly diverse as overall mitotic length (Fig. S2B). Cells that took longer between metaphase plate formation and final division showed intense spindle rocking before anaphase onset (Fig. S2C, Movie 2). Moreover, final division angles showed a wide distribution variety (n=377, N=21) (Fig. 1D,G), suggesting that HCprs can divide with any directionality. Nevertheless, a bias towards 0° and 90° division angles was seen. No strong correlations were observed between division angle and timing of spindle positioning (Fig. S2D,E). In addition, the position of mitotic cells showed no strong connection to mitotic length (Fig. 1J), suggesting that duration of spindle positioning was the main factor in mitotic length distribution. Together, our data shows that division position, duration of final mitosis and spindle orientation are not stereotypical but vary widely among HCprs. However, independently of these variations, the resulting HCs consistently reach their correct layer (Fig. 1B, Movie 1).

The fact that divisions of HCprs can, in principle, occur at any position within the INL made us ask whether these cells follow an intrinsic cell cycle timer that is responsible for cell divisions occurring independently of HCpr position. To test this idea, we tracked HCprs from apical birth to final division (n=23, N=10) (Fig. 2A, Movie 3). Cell cycle duration ranged from 11 h 15 min to 16 h 41 min, with an average of 13 h (Fig. 2C). Therefore, HCpr cycles were much longer than those reported for progenitor cells, which, albeit varied, have an average cell cycle length of 7.3 h at 30 hpf and 5.3 h at 42 hpf (Matejčić et al., 2018). As we used light sheet microscopy to analyse cell cycle timing, it is unlikely that differences were due to increased phototoxicity (Icha et al., 2016, 2017). No clear correlation between cell cycle length and distance from apical surface during division was seen (Fig. 2D).

Studies of fixed chick retinae suggested that HCprs undergo S phase during apical-basal migration and arrest in G2 while stalling in the INL (Boije et al., 2009). We tested this notion using PCNA (proliferating cell nuclear antigen) as a marker to unambiguously follow HCpr cell cycle phases (Leung et al., 2011) (Fig. 2A,B, Movie 3). Interestingly, we found that HCprs passed G1 and S phase at diverse positions along the INL (Fig. 2B). Cell were seen in G1 in the INL when HCprs are already multipolar. In addition, HCprs could be in S phase during basal as well as apical migration (Fig. 2B). In addition, we found that G2 [defined as time between disappearance of PCNA dots and cell rounding (Leung et al., 2011)] covered a rather short span of the overall cell cycle length (Fig. 2E) often occurring close to the OPL (Fig. 2B).

Together, these data make it unlikely that HCprs arrest in G2 before apical migration in zebrafish. Furthermore, our results argue against an intrinsic cell cycle timer responsible for mitotic onset (Fig. 2).

It has previously been reported that maturing HCs undergo bidirectional migration before final lamination (Chow et al., 2015; Edqvist and Hallböök, 2004; Weber et al., 2014). However, so far only small cohorts of cells have been analysed using confocal microscopy with time resolution of 15 min or longer and focusing on different movement phases (Chow et al., 2015). To understand the stereotypicity of HC single cell migration, we followed cells from birth to final positioning using light sheet microscopy at 5 min intervals (from here on we will use the term HC for all migration behaviour before and after final division) (n=21 cells N=11) (Fig. 3A, Movie 4). The timing between birth and final positioning of HCs varied, broadly spanning between 9 h and almost 16 h (Fig. 3B). In addition, timing between birth, retraction of the apical process and final position varied (Fig. S3B). Single cell trajectory analysis revealed that, in contrast to the highly stereotypic migration behaviour of RGCs (Icha et al., 2016), each HC followed its own individual path and migration tracks took varying amounts of time (Fig. 3C, Fig. S3A). Trajectories are faster and more directed at early migration stages, as verified by mean square displacement and directionality analysis (Fig. 3C,D, Fig. S3C), most likely due to the maintenance of HC apical attachment (Chow et al., 2015). However, trajectories become more diverse and less directed after this stage (Fig. 3C,D, Fig. S3C). Tangential migration differed significantly among HCs, showing some cells are more constrained in nasal-temporal direction and others moving substantially along this axis (Fig. 3E). Despite these single cell differences, all HCs robustly reached their correct layer where they start polarizing (Fig. 3A, Movie 4).

Similar to other parameters, the depth of migration also varied strongly among individual HCs (Fig. 3C,F, Fig. S3). To find tissue-wide parameters responsible for these differences, we analysed how migration depth correlates with the formation of the IPL, which starts around 54 hpf and contains processes connecting amacrine cells (ACs), BCs and RGCs (Godinho et al., 2005). When defining the distance between the final OPL and IPL as 100% of the migratory path, shallow cells reached less than 50% of this distance, while HCs moving more deeply could surpass the prospective IPL position by 20% (Fig. 3F,G). However, HC migration deeper than the prospective IPL occurred only before this structure emerged (Fig. 3G,H). Early HCs that pass the depth of the prospective IPL were able to return to apical positions (Fig. 3H). Once the IPL was formed, however, HCs did not migrate below it (Fig. 3G,H). Interestingly, regardless of migration depth, all HCs were in contact with the AC layer (Fig. 3H and Movie 1), suggesting that migration into the AC layer is necessary for returning apically. As HCs cannot go deeper than the IPL after its formation, this structure most likely acts as a steric hindrance for HCs.

Despite the variations found for HC cell cycle and migration parameters, HC layer formation occurred robustly in all embryos imaged (N=21). We therefore speculated that sister cell dispersal is important to fill the HC layer independently of single cell behaviour. We thus expected that this dispersal showed similar heterogeneity as other parameters to buffer the previous cell cycle variability and stochastic migration patterns observed. This was indeed the case (n=20, N=6). Although some sister HCs moved substantially in apicobasal direction after terminal divisions and ended up in close nasal-temporal vicinity (Fig. 4A,B and Movie 5), others moved very little in the apicobasal direction but dispersed much more laterally with final positions of sister cell bodies one-tenth of a micrometer apart (Fig. 4C,D). Lateral HC sister cell movement took between 1 h 40 min and 6 h 30 min (Fig. 4E), and lateral dispersion differed substantially among cohorts of sister cells upon final mitosis, ranging from 0.6 µm to 29 µm (Fig. 4F).

Thus, like other parameters, sister cell dispersion within the HC layer is highly heterogenous. Interestingly, sister HCs can end up far apart from each other, arguing against the idea that functional retinal columns exist, at least when it comes to HC clones (He et al., 2012), as suggested previously for mouse horizontal cells (Reese et al., 1995).

Overall, the use of long-term quantitative light sheet microscopy enabled us to present first evidence that, in addition to stochasticity in retinal fate decisions (Baye and Link, 2007; Del Bene et al., 2008; He et al., 2012), later stages of retinal development, in particular retinal laminations, also show high degrees of variability. The example presented here, the formation of the HC layer, demonstrates this to an extreme extent: all parameters analysed over hundreds of cells from many embryos showed heterogeneity. The term stochasticity applies, as single cell behaviour cannot be predicted. Instead, some degree of random chance seems to govern division timing and position and particularly the migratory path of the cells. Our findings are summarized in Fig. 4G.

At this point, we cannot completely exclude the possibility that common parameters could be found in different subsets of HCs due to the lack of subset-specific live imaging markers. However, the broad spread of our data argues against this possibility. To resolve this point, single cell sequencing of HCs would be an attractive possibility. Owing to the scarcity of HCs within the retinal cell population, however, this might be a demanding task.

Other cell types that laminate in the vertebrate retina show much more stereotypic behaviour than HCs. RGCs laminate via somal translocation, and their timing and migration trajectories are highly stereotypic (Icha et al., 2016). What could be the reason for these differences? We speculate that one cause could be the presence or absence of an apical cell attachment. RGCs are apically attached until they reach the RGC layer and only upon apical process detachment do they undergo random walk-like behaviour during fine positioning (Icha et al., 2016). HCs instead lose their apical attachment early during migration (Chow et al., 2015) (Fig. 4G) and from then on navigate freely through the tissue without any predisposed track information. Thus, timing and depth of migration can vary and this heterogeneity needs to be buffered by the heterogeneity of other parameters such as sister cell dispersal. Together, this ensures robust HC layer formation despite single parameter variability. Consequently, although the behaviour of single cells is not predictable, at the tissue level it still produces a robust pattern. The obvious next steps will be to understand how this emergence of a well-organized tissue from seemingly chaotic single cell behaviour is ensured. To achieve this, it is necessary to extract the tissue-wide parameters guiding it. Such parameters could be coded in chemical guidance cues; however, mechanical tissue properties could also play a role. Follow up studies will therefore raise an exciting new array of questions, such as which additional neuronal pattern formation phenomena involve variability of the cell cycle, the positioning parameters and stochastic migration patterns, and how this single cell variability forms reproducible lamination at the tissue level.

HCs have been identified in vertebrate retinas from fish to humans and are divided into two main types: axon bearing and axon less. In the adult zebrafish retina, three types of axon-bearing HCs (H1, H2 and H3) and one type of axon-less HC (HB) have been identified (Song et al., 2008). HCs can be labelled for different proteins and transcription factors (TFs) that are often expressed in a subtype-specific manner (Fig. S1A). In zebrafish, the gap junction protein connexin55.5 is expressed in HCs. Although some studies propose that this marker is predominantly expressed in HA types of HCs (Godinho et al., 2007; Shields et al., 2007), a more recent study suggests that it is expressed in more diverse subtypes (Klaassen et al., 2016). The pancreas-specific transcription factor 1a (Ptf1a) is expressed in all HC subtypes in the zebrafish retina (Jusuf et al., 2011), while the TF thyroid hormone receptor β2 (trβ2) labels 70% of Ptf1a-positive HCs in the zebrafish retina (Suzuki et al., 2013). In the chick retina, the Lim family TFs Islet1 (Isl1) and Lim1/Lhx1 are expressed in axon-less and axon-bearing HC subtypes, respectively (Edqvist et al., 2008). Whether Lim1/Lhx1 has a similar heterogenous expression pattern in HC subtypes in zebrafish is not currently known but we see that the majority of Ptf1a-positive HCs are labelled by this marker. Here, we use a combination of these HC markers (Fig. S1B,C), and show that no marker-specific HC behaviour is observed (Fig. S1E).

Zebrafish were maintained and bred at 26.5°C. Zebrafish embryos were raised at 28.5°C or 32°C and staged in hours post fertilization (hpf) according to Kimmel et al. (1995). Embryos were kept in E3 medium, which was supplemented by 0.003% 1-phenyl-2-thiourea (PTU) from 8±1 hpf onwards and renewed daily. PTU was added to prevent pigmentation. Embryos were anesthetized in 0.04% tricane methanesulfonate (MS-222; Sigma-Aldrich) prior to sorting for fluorescent signal (>35 hpf). Refer to Table S1 for a list of transgene lines used in this study. See Table S3 for a list of experiments, and numbers of cells and embryos. All animal work was performed in accordance with the European Union (EU) directive 2010/63/EU, as well as the German Animal Welfare act.

To mosaically label HCs, DNA constructs were injected into the cytoplasm of one-cell stage embryos. Constructs were diluted in double-distilled H₂O supplemented with 0.05% Phenol Red (to visualize the injection material) (Sigma-Aldrich). Injected volumes ranged from 1 to 1.5 nl. DNA concentrations were 20-30 ng/µl and did not exceed 30 ng/µl even when multiple constructs were injected. See Table S2 for a list of injected constructs.

To induce expression of heat shock promoter (hsp70)-driven constructs and transgenes, embryos were incubated in a water bath 1.5-3 h before imaging for 20 min at 39°C. See Tables S1 and S2 for a list of heat-shock induced transgenes and constructs, respectively.

Transplantation dishes were prepared by floating a plastic template in a Petri dish that was half-filled with 1% low-melting-point agarose in E3. When the agarose solidified, plastic templates were gently removed, leaving an agar mould that contained rows of wells to hold embryos. Embryos at stages high to sphere were dechorionated in pronase (Roche) and dissolved in Danieu's buffer. Dechorionated embryos were transferred to wells in agarose moulds using a wide-bore fire-polished glass pipet. At the 1000-cell stage, cells from the donor embryos were transplanted into the animal pole of the acceptor embryos using a Hamilton syringe. The animal-pole region is fated to become retina and forebrain tissue. Transplanted embryos were kept on agarose for about 3-5 h and then transferred onto glass dishes that contained E3 medium supplemented with 0.003% PTU and antibiotics (100 U of penicillin and streptomycin, Thermo Fisher Scientific). Transplanted embryos were identified via fluorescence and imaged from 42 hpf (see Table S3).

Imaging started between 35 hpf and 48 hpf (see Table S3). Embryos were manually dechorionated and mounted in glass capillaries in 1% low-melting-point agarose as previously described (Icha et al., 2016). The sample chamber was filled with E3 medium containing 0.01% MS-222 (Sigma-Aldrich) (to immobilize embryos) and 0.2 mM PTU (Sigma-Aldrich). Imaging was performed on a Zeiss Light sheet Z.1 microscope with a Zeiss Plan-Apochromat 20× water-dipping objective (NA 1.0) at 28.5°C. Z-stacks spanning the entire eye (100-120 μm) were recorded with 1 μm optical sectioning. Z-stacks were recorded every 5 min for 6-42 h, with double-sided illumination mode. The system was operated by the ZEN 2014 software (black edition).

The raw LSFM data were deconvolved in ZEN 2014 software (black edition, release version 9.0) using the Nearest Neighbour algorithm. Minimal image pre-processing was implemented prior to image analysis, using open source ImageJ/Fiji software (fiji.sc). Processing consisted of extracting image subsets or maximum intensity projections of a few slices. Processed files were then analysed in Fiji and the data were plotted using MATLAB or Microsoft Excel.

HC rounding was defined as the onset of mitosis and the point of anaphase as the end of mitosis. The time required for this was defined as the mitosis length (Fig. 1C,F).

To visualize HC nuclei, Tg(lhx-1:eGFP) animals were crossed with Tg(hsp70:H2B-RFP) animals. Metaphase was detected by the condensed and elongated nature of the chromatin in the centre of the cell. Anaphase was assessed by chromatin separation (Fig. S2C).

To assess cell cycle phase of HCs, the Hsp70::EGFP-PCNA DNA plasmid was injected into one-cell stage embryos. Individual PCNA⁺ HCs were manually followed in Fiji from birth to terminal mitosis. Onset of G1 was defined as the time of HC birth at the apical side. The end of G1 was defined as last timepoint prior to PCNA foci appearance (Fig. 2A). S phase was defined between the appearance and disappearance of PCNA foci (Fig. 2A). G2 was measured from the disappearance of PCNA foci (end of S phase) to cell rounding (Leung et al., 2011).

Division orientation was analysed manually using the angle tool in Fiji RRID:SCR_002285 (Schindelin et al., 2012) from maximum intensity projections of five slices of LSFM images. For each observed anaphase cell, a line was drawn parallel to the division plane, and the angle was measured with respect to the temporal-nasal tissue axis (horizontal image axis) (Fig. S1D). Data were plotted as a polar histogram in MATLAB.

The Fiji Line tool was used to manually measure the distance from the centre of the mitotic HCpr (at rounding) and the OPL in the maximum projection of three to five consecutive z slices (Fig. S1D).

First, maximum projected sub-stacks (five z slices) of the raw live images were generated in Fiji. xy-drift of 2D stacks was then corrected using a manual drift correction Fiji plug-in (imagej.net/Manual_drift_correction_plugin).

Three approaches were used to label HCs: (1) double transgenic animals were generated by crossing transgene lines in Table S1 [Tg(lhx-1:eGFP);Tg(ath5:gap-RFP),Tg(lhx-1:eGFP);Tg(actb2:mCherry-Hsa.UTRN) and Tg(lhx-1:eGFP); Tg(ptf1a:DsRed) (Table S3)]; (2) blastomere transplantation was performed to mosaically label HCs; and (3) microinjection of plasmid DNA encoding trβ2:tdTomato or Connexin:55.5:rasGFP was used to mosaically label HCs in the retina (Tables S2 and S3).

The migrating HCs were manually tracked by following the centre of the cell body in 2D images using MTrackJ plug-in in Fiji (Meijering et al., 2012). The resulting trajectories were analysed as described previously (Icha et al., 2016). Mean square displacements and directionality ratios were calculated using the DiPer program (Gorelik and Gautreau, 2014) and executed as a macro in Excel (Microsoft).

The relative percentage of maximum HC depth was defined by the apex of HC basal position using individual HC trajectories. At this point, the distance of the centre of the cell body from the OPL was divided by the OPL to IPL thickness (Fig. S1F).

All statistical analysis was performed using Prism software (version 6.0c for Mac OS; GraphPad Software) to create graphs. In total, data were collected from 432 HCs in 24 independent time-lapse imaging experiments (Table S3). Pearson correlation coefficient was used to examine associations between length, position and angle of HCpr mitosis as well as HCpr cell cycle length and position. The correlation coefficient, r, ranges from −1 to +1, where 1 indicates perfect correlation, 0 to 1 indicates the two variables tend to increase or decrease together, 0 indicates that the two variables do not vary together at all, −1 to 0 indicates that one variable increases as the other decreases and −1 indicates perfect negative or inverse correlation.

We thank M. Matejcic, I. Yanakieva and the Norden lab for useful project discussions and helpful comments on the manuscript. C. Zechner is thanked for discussion on stochasticity. We are grateful to I. A. Deniz, H. Hollak, S. Kaufmann, the Light Microscopy Facility, the Scientific Computing Facility and the Fish Facility of the MPI-CBG for experimental help. We thank S. Maddu for help with analysis, and Wong and L. Godinho for generously sharing constructs.