A Toggle-Switch and a Feed-Forward Loop Engage in the Control of the Drosophila Retinal Determination Gene Network

Dipterans show a striking range of eye sizes, shapes and functional specializations. Their eye is of the compound type, the most frequent eye architecture in nature. The development of this compound eye has been most studied in Drosophila melanogaster. The early development of the Drosophila eye is under the control of a gene regulatory network of transcription factors and signaling molecules called the retinal determination gene network (RDGN). Nodes in this network have been found to be involved not only in the development of different eye types in invertebrates and vertebrates, but also of other organs. Here we have analyzed the network properties in detail. First, we have generated quantitative expression profiles for a number of the key RDGN transcription factors, at a single-cell resolution. With these profiles, and applying a correlation analysis, we revisited several of the links in the RDGN. Our study uncovers a new link, that we confirm experimentally, between the transcription factors Hth/Meis1 and Optix/Six3 and indicates that, at least during the period of eye differentiation, positive feedback regulation from Eya and Dac on the Pax6 gene Ey is not operating. From this revised RDGN we derive a simplified gene network that we model mathematically. This network integrates three basic motifs: a coherent feedforward loop, a toggle-switch and a positive autoregulation which, together with the input from the Dpp/BMP2 signaling molecule, recapitulate the gene expression profiles obtained experimentally, while ensuring a robust transition from progenitor cells into retinal precursors.


INTRODUCTION
During organ development, specification of cell fates depends on gene regulatory networks (GRNs). These GRNs represent as directed graphs biochemical reactions that result in changes in gene expression (mRNA and protein production), which ultimately control cell function. GRNs comprise intracellular as well as extracellular components. Within a cell, nodes represent transcription factors (TFs) that regulate further tiers of genes, and the links connecting the nodes represent the activation/repression action of TFs on their target genes. Extracellular signals modify the activation/repression rates and thereby are key modulators of the dynamics of these GRNs. In general, GRNs operating during organ development must account for several biological phenomena: the generation of patterns of gene expression in space and time and the reliability in the generation of these patterns (i.e., robustness). Equally important, variations in the GRN of a particular organ underlie the evolutionary changes in the morphology and function of this organ (Levine and Davidson, 2005;Smith et al., 2018).
The Drosophila eye has been used as a paradigm to describe and study a GRN controlling organ specification and early development, the so-called retinal determination gene network (RDGN) (Silver and Rebay, 2005;Kumar, 2009;Amore and Casares, 2010;Casares and Almudi, 2016). Major genes in this network have been implicated in the development of other organs in Drosophila and, interestingly, in vertebrates, suggesting some degree of conservation in the processes controlled by those genes (Ikeda et al., 2002;Zhang et al., 2002Zhang et al., , 2006bBessarab et al., 2004;Purcell et al., 2005;Bumsted-O'Brien et al., 2007;Kaiser et al., 2007;Erickson et al., 2010;Pineiro et al., 2014;Spieler et al., 2014). The current Drosophila RDGN has been built compiling genetic (i.e., functional), expression and regulatory information. Genetic experiments include loss-and gain-offunction experiments. Some of the latter, performed through ectopic gene expression in other organs (mostly in the developing wing), showed that some of the RDGN genes were sufficient to drive eye development, with the Pax6 gene eyeless (ey) being a paradigmatic case of the capacity of a TF to re-specify tissues toward eye development (Halder et al., 1995;Czerny et al., 1999;Baker et al., 2018). Expression data included transcripts, protein products and transcriptional reporters. For some interactions, enhancer elements have been identified and direct biochemical proof of transcription factor binding obtained, establishing direct regulatory links. This GRN should then result in specific spatiotemporal patterns of gene expression. However, the RDGN has not been challenged against a comprehensive quantitative analysis of the expression of its key transcription factors yet. This analysis can potentially identify inconsistencies between the GRN and its actual output (gene expression patterns) or give support to network topologies, as has been shown, for example, for the Drosophila embryonic segmentation network (Jaeger and Manu, 2012). Further, if the quantitative analysis comes from space-resolved single-cell data, it allows measuring important parameters such as gene expression variability ("noise") and expression correlations among genes, from which potential regulatory relationships can be inferred.
The eye primordium is a monostratified epithelium. Before differentiation onset, the progenitor state is characterized by the expression of the Meis1 homolog homothorax (hth) and of two paralogue pairs: the Pax6 genes ey and twin of eyeless (toy) plus the paralogues teashirt (tsh) and tip-top (tio) (Bessa et al., 2002(Bessa et al., , 2009Datta et al., 2009;Weasner et al., 2009). The network is animated by two secreted signals: Hh (Hedgehog) and the BMP2 Dpp (Decapentaplegic). Hh and Dpp are initially expressed at the primordium's posterior margin and facilitate the repression of hth and the upregulation of a set of nuclear/transcription factors that include eyes absent (eya/Eya), sine oculis (so/Six2) and dachshund (dac/Dach), so that progenitors are converted into cell cycle quiescent precursors. The RDGN culminates with the activation of the proneural gene atonal (ato) which is required for the further differentiation of precursor cells into photoreceptors ("R"). R cells express Hh, and Hh, in turn, induces the expression of Dpp. In this way, Hh and Dpp set in motion a differentiation wave that sweeps across the primordium leaving on its wake differentiating retinal tissue [reviewed in (Treisman, 2013)]. The front of this wave (that marks the transition from precursors to R cells) is characterized by an indentation of the epithelium, the morphogenetic furrow (MF), which acts as landmark for the wave-front. The wave/MF advances across the primordium at about constant speed during most of the differentiation process (except at the beginning and ending) (Campos-Ortega and Hofbauer, 1977;Basler and Hafen, 1989;Wartlick et al., 2014;Vollmer et al., 2016). One important consequence of this fact is that gene expression patterns remain stationary relative to the MF throughout most of the process. This allows, in principle, to use data from different time points -as long as the initiation and ending of the process are not included-to generate gene expression curves, registering all data relative to the MF.
In this paper we have generated registered spatial expression curves for hth, optix, ey, eya, dac, and ato in the Drosophila eye primordium, extracting expression information as fluorescence intensity from single-nuclei using laser confocal microscopy data. By performing double-labeling experiments, correlated data for gene-pairs was obtained to analyze potential regulatory relationships. With these data at hand, we revise the current gene regulatory model and explore, using a computational model, a core network topology that might confer, simultaneously, bistability and noise reduction properties to the RDGN.

Genotypes and Genetic Manipulations
Eye imaginal discs from the Drosophila melanogaster wild type strain Oregon-R were used when no transgenes were present in the genotype. Transgenes used: The optix2/3-dGFP enhancer reporter transgene is described in Ostrin et al. (2006) and recapitulates the endogenous expression of Optix (see Supplementary Figure 1). To generate the 3 ′ ato-destabilized GFP (3 ′ ato-dGFP) enhancer reporter transgenic line, the 3 ′ ato enhancer sequence, which controls the onset of ato expression anterior to the morphogenetic furrow (Zhang et al., 2006a), was PCR-amplified from genomic DNA using the primer set: FW: ATCGGGAGCAGTAACAAACTTAAC and RW: ATCTCCAT CCTCAATCAAAGCTAC and cloned into pCR8-TOPO. The cloned fragment was then transferred into the pBPUw-dGFP Gateway integration vector (Royo et al., 2011) according to the manufacturer's recommendations (Invitrogen). DNA constructs were microinjected into embryos from flies carrying the landing platformZH-attP-22A (Bischof et al., 2007), using standard Drosophila transformation techniques. The hth-YFP protein trap strain [CPTI-001356;Flannotator (Ryder et al., 2009)] was used in some experiments to follow hth expression and is described in (Choo et al., 2014).

Gene Expression Profile Generation, Including Segmentation and Alignment
Although gene expression profiles were generated from discs of different strains, we expect all relative positions and correlations/anticorrelations of the gene expression profiles to be the same irrespective of strain. Samples imaged by LCM were processed in order to retrieve DAPI-stained nuclei in the three dimensions within a "sampling volume" (see Figure 1). This sampling volume spans the central region of the disc, avoiding the inclusion of the lateral folds of the disc, as these would complicate the analysis, contained the anterior-most region of the eye disc and straddled the morphogenetic furrow. An ad hoc software, iFLIC, was designed that scans for ellipsoid intersection patterns with the confocal planes in any position and orientation within a rank of semiaxis lengths (RSL). These RSL were determined by a pre-scanning in which blobs of variable radius were tested. The patterns were checked for voxels whose value is greater than an optimized threshold, i.e., Otsu. Once these patters were established, a limited flooding was performed around the detected ellipsoids in order to cover the nuclei volume. The software also allowed performing a manual user transformation of the coordinates (i.e., the centroids of nuclei segments) which results in the selection of anterior-posterior stripes of constant width centered at the MF, i.e., the transformed reference frame (TRF). Tissue folds (which generally would introduce unrelated variability due to differential allocation of cells) were later corrected by the estimation of the z-coordinate in the TRF using a Gaussian radial function with an optimal bandwidth (shape parameter) and the centroids as positional data, which rendered a surface model of the selected volume. Centroids were projected onto the closer position of this surface and their trajectories within (up to the TRF origin) computed for both the x and y directions, replacing the original coordinates. The full description of iFLIC, the image segmentation software, can be found in Sanchez-Aragon and Casares (2019).
Finally the lack of calibration for each transcription factor signal was further corrected by a method we call Linear Scaling Minimization, in which variability across samples was reduced by linear transformations (implemented separately in an R package, Riflic, described in Supplementary Material). This produced the final data set for the model validation. This software can be downloaded at http://www.pvcbacteria. org/maxf/. The mean gene expression profiles shown in Figure 3 for each gene are the assembled profiles using all data sets for that gene. However, the correlation and coefficient of variation metrics for gene-pairs were computed exclusively using data from samples co-stained for each gene pair. These co-stained data are shown individually within the Supplementary Material on Mathematical Methods for Variability Reduction in Spatially Distributed Samples of Cell Segmentation.

Model and Model Analysis
In order to reproduce qualitatively the gene expression profiles obtained experimentally (see Figure 5B), we modeled the dynamics of the concentrations of the key transcription factors FIGURE 1 | Quantification of gene expression using confocal microscopy data in eye discs. (A) Confocal image of a late third larval stage (L3) eye disc, stained for Hth, Eya, and the photoreceptor marker Elav. These genes mark three major domains in the primordium: progenitors ("G"), precursors ("P"), and differentiating retina ("R"). The dashed line marks the approximate location of the morphogenetic furrow (MF). The retinal determination gene network (RDGN) comprises the regulatory interactions that span the specification of progenitors and their transition to atonal-expressing retinal precursors. "a" and "p" indicate anterior and posterior, respectively. The confocal data used in this study included the nuclear marker DAPI. Nuclei are segmented in 3D using iFLIC, a custom made software of the RDGN using coupled differential equations. The model equations were built as follows: The regulatory links are modeled as Hill equations, of the type: for activation and for repression. We consider the production rate of protein A controlled by a single transcription factor B and, therefore, rate of production of A is equal to f(B). n is the Hill coefficient, β is the maximum A production rate, and K, which is named the activation/repression coefficient and has units of concentration, is the concentration of B necessary to obtain β/2-i.e., halfmaximal A production rate.
In the following equations, the subindex "X" makes reference to the TF X, and the subindex "XY" to a (positive or negative) interaction from X to Y.
The equations that describe the coherent feed-forward loop (cFFL) for Pax6 ("P"), Eya:So ("E"), and target ("T") are: in the form P activates both E and T, while E activates T and besides, E is self-regulated. The equations contain different parameter types: the first term B k (k = E, T) accounts for the basal production rate; β k (k = E, T) is the maximum production rate; α k (k = E, T) is for the degradation/dilution rates; and hence the last term, α k K, represents the decay of each TF (including degradation and dilution as cells grow). The profile of P distribution is an input in this case, and was modeled as a sum of two Gaussian functions in order to approximate qualitatively the experimental profile: Equation (3) describes the concentration variation with time of E. This variation is determined by P activation (second term of the equation) and E's autoregulation (third term). As both terms are independent of each other they appear as a sum. The selfregulation allows the concentration of E to be maintained even if P concentration falls to zero. Equation (4) describes the concentration variation with time of T. The product of the contribution of P and E represents the AND integration logic. This "AND" logic of the FFL imposes a delay in T because it makes necessary minimum concentrations of P and E to activate T. If the integration logic becomes "OR, " T concentration appears earlier, even before that of E. In addition, an OR logic would cause T be expressed even after P were no longer expressed.
Next, including the toggle-switch motif that represents the mutual repression between E (Eya/So) and H (Hth), the model equations are: with (8) describing the final dynamics of the target T.
Frontiers in Ecology and Evolution | www.frontiersin.org Finally, we introduce the action of the Dpp morphogen signaling ("M"). Dpp is produced at the differentiation wavefront, the morphogenetic furrow (MF), that separates anterior proliferative undifferentiated cells from posterior differentiating, cell cycle quiescent, retinal cells. Dpp signaling is required for Hth ("H") repression. As the Dpp signaling cascade is transduced in the nucleus by a transcription factor [Mad (Wiersdorff et al., 1996)], we model the action of Dpp signaling as a repressor transcriptional input on H, and approximate the distribution of Dpp signaling as a Gaussian representing the spatial distribution of the active form of Mad [see (Neto et al., 2016) and references therein]. The model is now described by the following equations, where the repression of M on H has been included. Note that the integration of the repressor links from E and Dpp ("M") on Hth ("H") (Equation 10) are additive, equivalent to an "OR" logic : In our model, the profile of M distribution has been approximated as a Gaussian function as this is a one-dimensional diffusion process centered around the MF. M is simulated as: The full model is represented by the differential equations (9), (10), and (11), being M and P the inputs to the system. In order to obtain a set of values for all the parameters capable of simulating computationally the expression profiles obtained from experiments, the model was implemented using Vensim software (Vensim PLE, Ventana Systems, https://vensim.com/ vensim-software/), a visual tool for solving differential equations that allows modifying parameter values in run-time. Note that the values of the parametric set have been manually adjusted to obtain qualitatively similar profiles to the experimental ones. Once this was done (see Supplementary Figure 2 for the Vensim model and parameter values used), the model was also implemented using MATLAB (http://www.mathworks.com/). To test whether the network is stable with respect to noise and whether the FFL acts as noise filter, we added white Gaussian noise to the inputs of the system-that is, to M and P. This implies that E, H, and T become stochastic variables. We computed a sufficiently large number of trajectories and calculated the standard deviation, which is a function of time. After averaging in time, we compared it to the noise in the input and obtained that the propagated noise constitutes only an ∼8.5% (See below).

Image Analysis Pipeline: Obtaining Gene Expression Profiles
Our first aim was to describe quantitatively changes in the expression of key RDGN genes as cells get closer to the differentiation wave. As a correlate of gene expression we used immunofluorescence intensity signal obtained from dense confocal imaging stacks. Rather than mean fluorescence profiles, we aimed at obtaining single cell expression data. This type of data is equivalent to performing cell cytometry but preserving spatial information, and allows a precise measurement of gene expression variation among cells. In addition, to obtain gene expression correlation profiles, we obtained data from eye discs co-stained with pairs of genes.
As a first step, we used an image analysis pipeline, based on watershed, to segment nuclei in 3D, using DAPI as nuclear marker, and recording the information on the position of the centroid of each nucleus. Next, each nucleus is assigned gene expression values as fluorescence intensity (Figure 1) (Naval-Sanchez et al., 2013;Sanchez-Aragon and Casares, 2019). Although the eye primordium (or "eye disc") is an epithelial monolayer, it becomes folded as development proceeds. In addition, even using spacers for mounting the samples, there is some experimental variation in the degree of folding of the discs. Therefore, to minimize the impact of this folding on the registration of gene expression profiles from different samples, we implemented a computational "flattening" procedure of the data. With it, nuclei are connected orthogonally with the mean plane going across the sample (Figures 2A-C) and then their position is reassigned as the plane is flattened. The correlation analysis of pairs of gene expression profiles before and after flattening shows that correlation is generally improved after flattening ( Figure 2D; see Supplementary Material for a full description of the procedure). The expression profiles were registered using the position of the MF as coordinate x = 0.

Link Inference Based on Expression Correlations
We generated expression profiles for key nodes of the gene network ( Figure 3A) using specific antibodies (Hth, Ey, Eya, and Dac), a hth protein-trap (hth-YFP) and transcriptional reporters for ato and optix [3 ′ ato-dGFP and optix2/3-dGFP (Ostrin et al., 2006)] (see Materials and Methods), which constitute a dense sampling of many key nodes of the gene network.
Samples were co-stained with pairs or triads of these genes (see Materials and Methods) and gene expression profiles relative to the MF, which is given a position x = 0, were quantified. The regulatory links may be inferred by analyzing the degree of correlation (or anti-correlation) between different nodes in the network, so that direct interactions show the highest degree of correlation (in the case of activating links) or anticorrelation (in case of inhibitory links) (Dunlop et al., 2008). As the distance in the connection graph increases, the signal (correlation/anticorrelation) blurs. For each correlation, we obtained three metrics: mean, weighted local correlation and Fano factor [variance (σ 2 )/mean], this latter a measure of expression variability or "noise" (Figure 3A). Noise may be particularly informative, as increasing noise is associated with fast expression changes (Munsky et al., 2012). In Figure 3A we represent all the spatial expression profiles and the noise profiles for each gene.
We first analyzed the correlation between the TFs Hth and Ey, which display high levels of expression during the progenitor (i.e., initial) state of the GRN. Both profiles keep a high correlation throughout the anterior region that drops slowly as cells approach the MF. This may point to a direct relation between the two TFs or their combined upregulation by a third gene that is directly linked to both Hth and Ey ( Figure 3B). In fact, the Tsh/Tio paralogues have been shown to regulate Hth (Bessa et al., 2002;Singh et al., 2002), and Ey, Hth, and Tsh might form a protein complex (Bessa et al., 2002;Datta et al., 2011). Therefore, this profile is compatible with a mutual interdependence between Ey, Hth and some other regulator, likely Tsh/Tio. We next checked the correlation between Ey and the activity of the optix enhancer, optix2/3-dGFP, which is a direct Ey's target (Ostrin et al., 2006). The destabilized GFP reporter ("dGFP") allows following the dynamics of GFP expression avoiding GFP perdurance. As expected, Ey and optix-dGFP profiles show peaks of positive correlation when both appear in the most anterior part of the disc and when they are shut down at the MF. This latter simultaneous loss of Ey and optix2/3 at the MF is consistent The function of ey and toy has been grouped as "Pax6." The mutual activating loop between Eya and So has been simplified as a positive feedback autoregulatory loop ("FB"). Coherent feedforward loop motif (FFL) and toggle-switch motif ("toggle") have been boxed. A potential repression from Hth to targets has been included (as light negative link), but is not required to retrieve the experimental pattern. "T" represents a generalized target gene, such as dac, stg, or ato. A gradient of Dpp, produced by the differentiation wavefront, animates the network. (B) Normalized experimental expression profiles. (C) Gene network model-generated expression profiles. Note that while the x axis in (B) is distance, x has time units in (C). This is because, as the differentiation speed is constant, the gene expression pattern (Continued) FIGURE 5 | observed in the disc is equivalent to the gene expression changes that a cell experiences as time passes (i.e., as the differentiation wave gets closer to the cell). This means that earlier times are equivalent to more anterior positions in the primordium. (D) Weakening the H to E repression within the toggle-switch (K HE from 0.1 to 0.9; indicated by the green repressive link) leads to a faster increase in E levels and, to a lesser extent, of T (orange and brown arrows, respectively). This results in the anterior shift of H. (E) More intense repression from M on H (K MH de 0.9 to 0.1; blue repressor link) also results in anterior shifts (arrows) of H (green) and E (orange) and T (brown). (F) Changing the regulatory integration logic in the FFL from "AND" to "OR" results in the premature expression of T. However, this expression does not extend to follow P, as a negative link from H to P prevents this further expansion. (G) Introducing noise in the Pax6 and Dpp profiles has minor effects on the expression profiles of the other nodes, including T (compare to C). See main text for details.
with Ey being a direct, linear optix activator -so that loss of Ey is closely followed by a loss of optix2/3 expression. The onset of optix2/3 expression is also positively correlated with Ey, with a peak of noise marking exactly the anterior correlation peak. However, optix activation is delayed relative to Ey. This delay is compatible with the need for a coactivator or the removal of a repressor. The anterior limit of Optix expression is almost mutually exclusive with Hth, which might point to Hth acting as the anterior repressor of optix. In order to test this hypothesis, we generated random clones of cells expressing Hth in an optix2/3-dGFP background. Indeed, Hth-expressing clones do repress optix2/3 activity in a cell autonomous manner (Figure 3C), in agreement with a repressor role for Hth. We note, though, that this repression ability is lost at the posterior border of optix2/3-dGFP expression, which coincides with the MF. Therefore, it is possible that signals from the MF, directly or through the regulation of gene intermediaries, prevent Hth from repressing this enhancer close to the MF. A third interaction that we analyzed in detail was that between Hth and Eya. It has been described that Hth and Eya repress each other (Bessa et al., 2002). This is reflected in the complementary pattern of their expression profiles ( Figure 3A). In addition, the two profiles show a strong anti-correlation, which suggests that the mutual repression is likely direct. This interaction seems critical in the GRN, as progenitors are maintained in their proliferative undifferentiated state as long as Hth expression is ON. Therefore, the transition to a precursor state requires the simultaneous loss of Hth and the upregulation of Eya. Eya, which lacks a DNA-binding domain, acts as co-transcriptional activator of So/Six2 (Pignoni et al., 1997;Ohto et al., 1999;Jemc and Rebay, 2007), itself a Pax6 target (Halder et al., 1998;Punzo et al., 2002). The Eya:So transcriptional complex then activates the transcription of their two genes, thereby establishing a positive feedback that maintains their expression and the precursor state. One of the targets of Eya:So is dac (Chen et al., 1997(Chen et al., , 1999. When all profiles are plotted together (Figure 4A), the onset of Dac expression follows that of Eya with a short delay, in agreement with it being an Eya:So target. Interestingly, the activation phase of Dac is associated to two noise peaks, instead of one ( Figure 4B). This might be indicative of two regulatory inputs separated in space acting as Dac activators. These two noise peaks coincide with the Ey and Eya expression maxima (Figures 4A,B). Recent work on the regulatory sequences of Dac and its upstream regulators has discovered two Dac enhancers, 5 ′ EE and 3 ′ EE. The 5 ′ EE is regulated by Ey, with a contribution of Eya:So, while the 3 ′ EE is activated by Eya:So and the Dpp signal produced by the MF (Pappu et al., 2005). Therefore, Dac regulation receives Ey input first and, only when Eya:So have been upregulated, also that of Eya:So. These two inputs, which are separate in time, might be responsible for the two noise peaks that we detect. However, the fact that Dac follows with a short delay the activation in Eya expression suggests that this regulatory architecture induces Dac expression only when both Ey and Eya:So are expressed (Pappu et al., 2005).
The endpoint of the RDGN is the activation of the proneural gene atonal (ato). Here we have monitored the transcriptional activity of the ato3 ′ enhancer, which is responsible for the initiation phase of ato expression (Zhang et al., 2006a;Tanaka-Matakatsu and Du, 2008). Ey and Eya:So have been shown to be direct regulators of this enhancer (Zhang et al., 2006a;Zhou et al., 2014), similar to what has been described for Dac. However, the expression of ato3 ′ -dGFP raises just anterior to the MF (Figures 2A-C, 4A). Hh, produced by the differentiating photoreceptors posterior to the furrow has been proposed as an ato activator (Dominguez, 1999;Heberlein and Treisman, 2000). Therefore, initiation of ato expression just anterior to the MF would result from the short range action of Hh, produced by cells just behind the MF, on Ey and Eya:So-expressing precursors.
Integrated with previous RDGN models (Silver and Rebay, 2005;Kumar, 2009;Amore and Casares, 2010;Casares and Almudi, 2016), these analyses result in a consensus Drosophila RDGN, shown in Figure 4C, comprising its major players. The GRN model includes activator feedback links from Dac to Eya:So and Ey/Toy, and from Eya:So to Ey/Toy. Experimental evidence for these links derives from ectopic gene expression. A suggested role for these positive feedbacks is to lock-in, once activated, a precursor state characterized by coexpression of Pax6, Eya, So, and Dac (Desplan, 1997). Our analysis of gene expression profiles does not allow us to clearly support or reject any such feedbacks, so we decided to test them experimentally. If there were a positive feedback between Dac and Ey, we would expect to lose Ey expression in Dac-mutant clones. However, this is not the case. In dac-null (dac 3 ) clones, Ey expression remains unchanged relative to the expression of control tissue surrounding the clones (Figures 4D, D'). This result contrasts with those obtained by Atkins et al. (2013), showing that in dacmutant clones Ey expression is derepressed. Our interpretation of their results is that dac-clones cause a delay in MF movement (Mardon et al., 1994;Bras-Pereira et al., 2016) and this results in a delay in the turning off of Ey, rather than a true derepression. To test if the feedback operates upstream of Dac, we generated clones expressing an Eya-specific RNAi and analyzed Ey levels, quantitatively, in individual nuclei. In eya-RNAi clones, labeled with GFP, Ey levels remained comparable to control ones, even where Ey expression should have been lost (i.e., behind the MF) ( Figure 4E). In fact, some eya-RNAi cells within the normal region of Ey/Eya co-expression express higher than normal Ey levels (Figure 4E'). These results confirm those obtained by the Mardon group (Atkins et al., 2013). Therefore, our results seem to indicate that during the development of the eye, positive feedbacks from Dac and Eya:So to Ey (and presumably Toy) are not necessary to maintain Ey levels.

A Core RDGN Integrates a Feed-Forward Loop and a Toggle Switch
The RDGN might have specific regulatory properties that would be functionally relevant beyond the biochemical properties of its individual TFs. To investigate what these properties might be, we focused on some central and well-established links to identify network motifs and study their dynamic properties. This "core" network is shown in Figure 5. In it, Ey and Toy are considered a single transcriptional function ("Pax6"), as these two genes have been shown to be partly redundant (Zhou et al., 2014;Lopes and Casares, 2015;Zhu et al., 2017). The Eya-So mutual positive feedback loop is represented by a single autoregulatory node ("Eya/So"). The output of the core network is a generalized target ("T"), as it seems that the regulatory logic that controls molecularly characterized enhancers of these targets, such as ato (Zhou et al., 2014), dac (Pappu et al., 2005) and string [stg; (Lopes and Casares, 2015)] is similar, being activated by both Ey and Eya:So. The mutual Hth/Eya:So repression is well established (Bessa et al., 2002;Lopes and Casares, 2010) and reinforced by the analysis presented here. In addition, it has been reported recently that Hth might also contribute to the repression of dac and stg (Bras-Pereira et al., 2015;Lopes and Casares, 2015), and that Hth hinders Ey's capacity to activate stg (Lopes and Casares, 2015). This activity was represented as an indirect repressive link from Hth onto the activation of the target ("T") by Pax6 ("P"). Next, we modeled this GRN, including a source of Dpp, the major signal produced at the MF as differentiation proceeds ( Figure 5A).
Ey ("P"), Eya:So ("E") and T are organized as a feed-forward loop (FFL) with Ey and Eya inputs having the same sign (in this case positive), making it a "coherent" FFL [cFFL (Alon, 2007)]. The integration logic of both inputs seems to be "AND." For example, the delayed expression of dac relative to Ey suggests that although targets typically require Ey, they only become expressed after Eya:So are also expressed . This cFFL loop has an addition: the autorregulation of E (Eya:So).
The GRN integrates another network motif: Hth ("H") and Eya:So ("E") repress each other. This regulatory structure is a "toggle switch, " a motif that allows the selection of either of two, mutually exclusive, states. In this case H:ON/ E:OFF or H:OFF/ E:ON, that correspond to proliferative progenitors or cell cycle quiescent precursors, respectively.
Finally, the eye primordium is polarized by the action of the moving MF, that produces Dpp, a BMP2 type morphogen ("M") (Gelbart, 1989). Dpp represses Hth at long distance (Bessa et al., 2002) and contributes to the repression of Ey at shorter range, just abutting the MF (Firth and Baker, 2009). The full model, including the intertwined toggle switch and the cFFL with a positive feedback (see Materials and Methods) recapitulates very well the measured experimental profiles (Figures 5B,C), indicating that, although not including all the TFs of the full RDGN, this core network captures the main regulatory processes occurring during normal eye specification. We tested specifically the role that the model predicts for the toggle switch between Hth ("H") and Eya/So ("E") by weakening the repression of H on E (K_HE from 0.1 to 0.9). In this case, E levels rise prematurely. This causes a reduction in the H domain (as E is a H repressor), a broadening of the E domain and an anterior shift of the target T, although less prominent (as T expression depends on the joint contribution of E AND P) (Figure 5D). Since the proliferative progenitor pool requires H for its maintenance, a weaker H-to-E link should result in a premature exhaustion of this pool and, as a consequence, in a smaller eye. Similarly, an increase in the intensity of repression of H by M also leads to a retraction of the H domain ( Figure 5E) and would result similarly in the reduction of the progenitor pool. Therefore, the regulatory links in the toggle-switch (including the action of Dpp signaling, that acts as a modulator of this motif) may have an impact on the final eye's size. When we change the logic operator within the FFL from "AND" to "OR" we observe that the delay in T expression relative to E disappears ( Figure 5F). However, T expression does not follow P as would be expected from a cFFL (in which P alone should be able to activate T even without E). This is because the model introduced a negative input from H to P, as shown for some Pax6 targets in the eye (see, for example, Lopes and Casares, 2015).

DISCUSSION
In this work, through a combination of quantitative singlecell imaging, correlation analysis, genetics and modeling, we have shown that a toggle-switch motif regulates the transition from progenitor to precursor cells. Hth exerts a general negative regulatory action on the establishment of the precursor state by regulating not only the expression of Eya, but also that of the Six3 TF Optix. The transition from progenitors to precursors is facilitated by Dpp. As Dpp is induced by differentiating photoreceptors through their production of Hh, ultimately it is the differentiating retina that controls this progenitor-toprecursor transition. As increasing the size of the eye requires the expansion of the progenitor pool, the Dpp-Hth link is key. Indeed, increasing the production of Dpp, itself an "eyepromoting" signal, is predicted to result in a smaller eye, as the progenitor cell pool would be prematurely exhausted. This control may be more complex than anticipated, as optix, itself controlled by Hth, also regulates Dpp production and signaling (Li et al., 2013). As hth is repressed, a coherent feedforward loop, accelerated by Eya/So mutual positive regulation, ensures a robust establishment of the precursor state and the onset of ato expression, the last step before the initiation of retinal differentiation.
The analysis of the network shows how Hth delays the engagement of the cFFL, through its participation in the toggleswitch. Here, a major role is played by Dpp which, acting as a repressor of Hth, tips the switch favoring Eya/So expression. The delay imposed by Hth seems especially important. As Hth maintains the Ey-expressing cells as progenitors (Pichaud and Casares, 2000;Bessa et al., 2002;Neto et al., 2017), the longer the delay in entering the FFL state, the larger the expansion of the progenitor pool, thereby having a potential impact on final eye size. The repressor link from Hth to T does not seem critical, as removing it just shortens the delay in T induction (not shown). In addition, our expression correlation analysis and genetic experiments further include Optix among the "pro-eye" transcription factors repressed by Hth. This repression may vary in intensity, as the region of overlap with Hth is larger for Optix than for Eya (and minimal for Dac). This fact reinforces the notion that Hth delays the full acquisition of a retinal precursor state. The fact that Optix has been recently shown to be necessary for Dpp signaling (Li et al., 2013) and that Dpp signaling contributes to Hth repression draws another feedback (Hth-Optix-Dpp-Hth) that requires further investigation. Finally, we noted that the establishment of the precursor cell state by the Pax6 genes should result in the repression of Hth. Indeed, the simultaneous attenuation of ey and toy through the combined expression of ey-RNAi and toy-RNAi in clones results in the maintenance of Hth (Figure 3D).
The coherent FFL imposes a delay in the activation of target genes downstream of Pax6 and Eya/So that ensures that these targets are activated only when the precursor cell state, characterized by the coexpression of Pax6 and Eya:So, has been established. The further addition of the "E" positive feedback loop on top of the FFL accelerates Eya and So expression and makes it independent of Ey, something that may allow the maintenance of this gene pair behind the MF, once the expression of Ey and Toy has been turned off. This is important in order to stay in the state H:OFF/ E:ON. In addition, the FFL may act as a noise filter (Gui et al., 2016). This noise may derive, for example, from fluctuations in biochemical rates or from genetic variation that could affect the regulatory constants that are included in our model. To test whether the network is stable with respect to noise, we carried out some computational experiments in which a random component was included in the network as white Gaussian noise, g(σ), on P [g(σ)·P], M [g(σ)·M] or both signals. This implies that E, H, and T become stochastic variables. The standard deviation (σ) was varied from low (5%) to high (20%) amplitudes with respect to the original value of P or M. We computed a sufficiently large number of simulations and, even with the highest noise, which would be considered greater than biological noise, added to both P and M, the behavior of the system is qualitatively similar to the case without noise, with a noise reduction of 91.5% (Figure 5G). This shows that the core network is capable of filtering noise efficiently. This is particularly interesting, because the network can absorb variations in two major inputs in the system: Ey/Toy ("P"), that initiates the FFL, and Dpp ("M"), which regulates the toggle-switch's output. One implication of this property is that the early RDGN can absorb the effects that might be caused by genetic variation affecting many of its parameter rates. This noise buffering property of the gene networks controlling early organ development has been observed experimentally analyzing the endomesoderm GRN in the sea urchin (Garfield et al., 2013). Thus, early organogenetic GRNs would be very robust, allowing the maintenance of genetic variation affecting these early stages, while morphological variation would rest more on terminal differentiation branches of the network (Garfield et al., 2013).

DATA AVAILABILITY
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The study involved wild type and transgenic Drosophila melanogaster strains.
The study did not entail ethical considerations.

FUNDING
This work was funded by MINECO and the Agencia Estatal de Investigacion (AEI) of Spain, co-financed by FEDER funds (EU) through grants BFU2012-34324 and BFU2015-66040-P to FC, MDM-2016-0687 in which FC is participant researcher, and TIN2017-89842 P in which MCL is participant researcher.