Abstract
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, ; Smith et al., ).
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, ; Kumar, ; Amore and Casares, ; Casares and Almudi, ). 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., ; Zhang et al., 2002, 2006b; Bessarab et al., ; Purcell et al., ; Bumsted-O'Brien et al., ; Kaiser et al., ; Erickson et al., ; Pineiro et al., ; Spieler et al., ). The current Drosophila RDGN has been built compiling genetic (i.e., functional), expression and regulatory information. Genetic experiments include loss- and gain-of-function 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., ; Czerny et al., ; Baker et al., ). 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 spatio-temporal 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, ). 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., , ; Datta et al., ; 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, )]. 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, ; Basler and Hafen, ; Wartlick et al., ; Vollmer et al., ). 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, bi-stability and noise reduction properties to the RDGN.
Materials and Methods
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. () 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: ATCTCCATCCTCAATCAAAGCTAC and cloned into pCR8-TOPO. The cloned fragment was then transferred into the pBPUw-dGFP Gateway integration vector (Royo et al., ) according to the manufacturer's recommendations (Invitrogen). DNA constructs were microinjected into embryos from flies carrying the landing platformZH- attP- 22A (Bischof et al., ), using standard Drosophila transformation techniques. The hth-YFP protein trap strain [CPTI-001356; Flannotator (Ryder et al., )] was used in some experiments to follow hth expression and is described in (Choo et al., ).
dac3(Flybase) mutant clones were induced by flip-mediated mitotic recombination (Xu and Rubin, 1993) by subjecting yw, hs-flipase 122; FRT40A dac3/FRT40A Ubi-GFP larvae to a 30′ heat shock (37°C). dac3-mutant tissue was marked by the absence of GFP.
Flip-out clones: the Flip-Out method (Struhl and Basler, ) was used to knock down Drosophila Pax6 paralogues simultaneously. RNA-mediated interference (RNAi) of ey expression was achieved using UAS-eyRNAi(II) (VDRC 106628), while toy expression was knocked down using UAS-toyRNAi/TM6B,Tb(VDRC 15919). Both UAS-RNAi transgenes were combined in a single line using standard genetic techniques. Females of the genotype y,w,hsFLP,Act5C(FRT.hsCD2)Gal4;;UAS GFP/TM6b,Tb were crossed to males carrying both UAS-RNAi for ey and toy. Clones were induced 24–48 h after egg laying (AEL) by a 15′ heat shock at 35.5°C. Next, larvae were grown at 29°C to maximize UAS-RNAi expression on standard medium and dissected 48 h later. In order to recover clones of cells simultaneously expressing both ey and toy RNAi, larvae of the following genotype were selected: y,w,hsFLP,Act5C(FRT.hsCD2)Gal4/+;UAS-eyRNAi/+;UAS-toyRNAi/UAS-GFP. To induce hth expressing clones, y,w,hsFLP,Act5C(FRT.hsCD2)Gal4 females were crossed to UAS-131-GFPhth (Casares and Mann, ) males. Clones were induced 24–48 h AEL by a 15' heat shock at 35.5°C, and larvae raised at 25°C until dissection.
Immunofluorescence and Imaging
For immunofluorescence, discs were processed essentially as in Casares and Mann (). For the samples aimed at obtaining quantitative gene expression profiles, one extra step was introduced. In order to improve nuclei segmentation during image analysis, discs were briefly (10′) incubated in 0.75×PBS on ice just before fixation. This “hypotonic shock” induced a slight swelling of the cells resulting in larger spacing between nuclei that made the segmentation of the nuclei after confocal imaging easier. Primary antibodies used were guinea pig anti-Hth at 1:3000 (Casares and Mann, ), rat-anti-Ey (1/150; gift from P. Callaerts), mouse anti-Eya (10H6 at 1/100 from Developmental Studies Hybridoma Bank, DSHB), mouse anti-Dac (1/500 Mabdac1-1, DSHB), chicken anti-GFP (1/500 Abcamab 13970), rabbit anti-GFP (1/1000 Molecular Probes A11122), rabbit anti-Optix [1/500, gift from F. Pignoni, SUNY Upstate Medical University (Kenyon et al., )] and rat anti-Elav (1:10007EBA10, DSHB). Fluorescently labeled secondary antibodies were from Molecular Probes (1/500). Laser Confocal Microscopy (LCM) was carried out on Leica SPE (all imaging for quantitative analysis) or Leica SP2 (data in Figures 3C,D, 4E) confocal microscope set-ups. The gene combinations imaged and the number of samples per combination were: Hth;Eya (19), Hth;Ey;ato-dGFP (8), Hth:YFP;Dac (8), Ey;optix2/3-dGFP (5), Eya;optix2/3-dGFP (5), Dac;optix2/3-dGFP (3), Eya;ato-dGFP (2), Dac;ato-dGFP (4).
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 ().
Figure 1
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 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., half-maximal 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 Bk (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, αkK, 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 self-regulation 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.
Finally, we introduce the action of the Dpp morphogen signaling (“M”). Dpp is produced at the differentiation wave-front, 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.,
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).
Results
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.,
Figure 2

Computer-aided stretching of the tissue improves the registration of gene expression profiles. The folding of the disc epithelium is an important source of analytical error. We have developed an iFLIC plug-in (see Supplementary Material) that solves sufficiently this problem by estimating the actual z coordinate for each nucleus centroid of the sampled tissue and calculates the orthogonal projection of each cell on that surface (A). The disc shown is from an ato-dGFP individual, stained for GFP (ato), Hth (yellow), Ey (red), and DAPI (blue). The results of the computer-aided stretching are shown in (B,B′: top views) and (C,C′: lateral views). B′ and C′ are the stretched data. (D) Comparison of the Pearson's correlation coefficient (PMCC) of gene profiles in several independent samples by combining them in pairs, either non-stretched (gray) or after stretching (red) (D). Stretching leads to a general improvement of the PMCC. PMCC was computed by Fisher's transform, showing confidence intervals at α = 0.05.
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.,
Figure 3

Correlation analysis of gene expression profiles. (A) Local weighted correlation for Hth/Ey, Ey/Optix, and Eya/Hth showing three parameters: their mean (solid line), their Fano's factor (as a measure of signal variability, or “noise,” dashed line) and their weighted local correlation (gray ribbon). Ribbon width corresponds to a confidence interval with α = 0.05. (B) Potential regulatory interactions that are compatible with the correlation analysis (see main text for details). (C) Hth-expressing clones (marked with anti-Hth) in an optix2/3-dGFP eye disc. Most hth clones repress optix2/3 cell-autonomously (yellow arrowheads). A clone adjacent to the MF does not show this repression (white arrowhead). Merged and optix2/3-dGFP-only channels are shown. (D) Clones simultaneously expressing ey-RNAi and toy-RNAi (marked in green, yellow arrowheads). In these clones Hth expression is upregulated. Merged and Hth-only channels are shown; a white line delineates two confocal planes obtained from the same imaginal disc. Nuclei in (C) and (D) are counterstained with DAPI.
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.,
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.,
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.,
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.,
Figure 4

The RDGN gene expression profiles. Normalized gene expression (A) and noise (Fano factor; B) profiles relative to the MF (x = 0). The two bars spanning A and B mark the two “noise” peaks for the dac gene. (C) Revised schematic RDGN. Genes are named as “Drosophila gene name/vertebrate gene homolog.” Colored nodes are genes studied in this work. Boxes group genes with partly redundant functions (tsh/tio; ey/toy) or working jointly (eya/so). Green links are regulatory relationships supported by our study. Blue links are suggested relationships, including a potential negative regulation from ey to toy ensuring constant levels of Pax6 function. Red links are not supported (see D,E), at least, for eye development during the third larval stage. (D) Disc containing dac3-mutant clones (marked by the absence of GFP, and outlined in D'), stained for Ey and Elav (to mark the differentiating retina). The expression of Ey within the dac-mutant tissue is indistinguishable from its expression in the surrounding control tissue (D′). (E) Clones attenuating eya function (with an eya-RNAi), marked in green. The disc is stained for GFP (eya-RNAi), Ey and nuclei (DAPI). Quantitative analysis of Ey expression with single cell resolution (E′) shows that Ey expression is derepressed in eya-RNAi cells.
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,
Integrated with previous RDGN models (Silver and Rebay,
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,
Figure 5

A gene regulatory network with intertwined feed-forward loop and toggle-switch motifs recapitulates the RDGN expression profiles. (A) The basic RDGN. 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 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 (KHE 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 (KMH 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.
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,
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,
Discussion
In this work, through a combination of quantitative single-cell 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-to-precursor 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 “eye-promoting” 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.,
The analysis of the network shows how Hth delays the engagement of the cFFL, through its participation in the toggle-switch. 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,
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.,
Statements
Data availability statement
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.
Author contributions
FC: Conceptualization. MS-A, JC-G, and MCL: Formal analysis. MS-A, CML, CSL, and CB-P: Investigation. MS-A: Software. MS-A, JC-G, CML, CSL, and CB-P: Visualization. FC: Writing–original draft. FC, MCL, MS-A, JC-G, CML, CSL, and CB-P: Writing–review and final version. FC and MCL: Supervision. FC: Funding acquisition.
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.
Acknowledgments
We thank ALMI (Avanced Light Imaging and Analysis Platform, CABD) for imaging support. CML would like to thank Prof. Dr. Isabel Correas (Universidad Autónoma de Madrid, Spain) for her constant support and generosity.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2019.00221/full#supplementary-material
References
1
AlonU. (2007). Network motifs: theory and experimental approaches. Nat. Rev. Genet.8, 450–461. 10.1038/nrg2102
2
AmoreG.CasaresF. (2010). Size matters: the contribution of cell proliferation to the progression of the specification Drosophila eye gene regulatory network. Dev. Biol.344, 569–577. 10.1016/j.ydbio.2010.06.015
3
AtkinsM.JiangY.Sansores-GarciaL.JusiakB.HalderG.MardonG. (2013). Dynamic rewiring of the Drosophila retinal determination network switches its function from selector to differentiation. PLoS Genet.9:e1003731. 10.1371/journal.pgen.1003731
4
BakerL. R.WeasnerB. M.NagelA.NeumanS. D.BashirullahA.KumarJ. P. (2018). Eyeless/Pax6 initiates eye formation non-autonomously from the peripodial epithelium. Development 145. 10.1242/dev.163329
5
BaslerK.HafenE. (1989). Dynamics of Drosophila eye development and temporal requirements of sevenless expression. Development107, 723–731.
6
BessaJ.CarmonaL.CasaresF. (2009). Zinc-finger paralogues tsh and tio are functionally equivalent during imaginal development in Drosophila and maintain their expression levels through auto- and cross-negative feedback loops. Dev. Dyn.238, 19–28. 10.1002/dvdy.21808
7
BessaJ.GebeleinB.PichaudF.CasaresF.MannR. S. (2002). Combinatorial control of Drosophila eye development by eyeless, homothorax, and teashirt. Genes Dev.16, 2415–2427. 10.1101/gad.1009002
8
BessarabD. A.ChongS. W.KorzhV. (2004). Expression of zebrafish six1 during sensory organ development and myogenesis. Dev. Dyn.230, 781–786. 10.1002/dvdy.20093
9
BischofJ.MaedaR. K.HedigerM.KarchF.BaslerK. (2007). An optimized transgenesis system for Drosophila using germ-line-specific phiC31 integrases. Proc. Natl. Acad. Sci. U.S.A.104, 3312–3317. 10.1073/pnas.0611511104
10
Bras-PereiraC.CasaresF.JanodyF. (2015). The retinal determination gene Dachshund restricts cell proliferation by limiting the activity of the Homothorax-Yorkie complex. Development142, 1470–1479. 10.1242/dev.113340
11
Bras-PereiraC.PotierD.JacobsJ.AertsS.CasaresF.JanodyF. (2016). dachshund potentiates hedgehog signaling during Drosophila retinogenesis. PLoS Genet.12:e1006204. 10.1371/journal.pgen.1006204
12
Bumsted-O'BrienK. M.HendricksonA.HaverkampS.Ashery-PadanR.SchulteD. (2007). Expression of the homeodomain transcription factor Meis2 in the embryonic and postnatal retina. J. Comp. Neurol.505, 58–72. 10.1002/cne.21458
13
Campos-OrtegaJ. A.HofbauerA. (1977). Cell clones and pattern formation: on the lineage of photoreceptor cells in the compound eye of Drosophila. Wilehm Roux Arch. Dev. Biol.181, 227–245. 10.1007/BF00848423
14
CasaresF.AlmudiI. (2016). Fast and Furious 800. The retinal determination gene network in drosophila, in Organogenetic Gene Networks, eds Castelli-Gair HombríaJ.BovolentaP. (Cham: Springer), 95–124. 10.1007/978-3-319-42767-6_4
15
CasaresF.MannR. S. (1998). Control of antennal versus leg development in Drosophila. Nature392, 723–726. 10.1038/33706
16
CasaresF.MannR. S. (2000). A dual role for homothorax in inhibiting wing blade development and specifying proximal wing identities in Drosophila. Development127, 1499–1508. Available online at: http://dev.biologists.org/content/127/7/1499.long
17
ChenR.AmouiM.ZhangZ.MardonG. (1997). Dachshund and eyes absent proteins form a complex and function synergistically to induce ectopic eye development in Drosophila. Cell91, 893–903. 10.1016/S0092-8674(00)80481-X
18
ChenR.HalderG.ZhangZ.MardonG. (1999). Signaling by the TGF-beta homolog decapentaplegic functions reiteratively within the network of genes controlling retinal cell fate determination in Drosophila. Development126, 935–943.
19
ChooS. W.BehC. Y.RussellS.WhiteR. (2014). Characterisation of Drosophila Ubx CPTI000601 and hth CPTI000378 protein trap lines. ScientificWorldJournal2014:191535. 10.1155/2014/191535
20
CzernyT.HalderG.KloterU.SouabniA.GehringW. J.BusslingerM. (1999). twin of eyeless, a second Pax-6 gene of Drosophila, acts upstream of eyeless in the control of eye development. Mol. Cell3, 297–307. 10.1016/S1097-2765(00)80457-8
21
DattaR. R.LuryeJ. M.KumarJ. P. (2009). Restriction of ectopic eye formation by Drosophila teashirt and tiptop to the developing antenna. Dev. Dyn. 238, 2202–2210. 10.1002/dvdy.21927
22
DattaR. R.WeasnerB. P.KumarJ. P. (2011). A dissection of the teashirt and tiptop genes reveals a novel mechanism for regulating transcription factor activity. Dev. Biol.360, 391–402. 10.1016/j.ydbio.2011.09.030
23
DesplanC. (1997). Eye development: governed by a dictator or a junta?Cell91, 861–864. 10.1016/S0092-8674(00)80475-4
24
DominguezM. (1999). Dual role for Hedgehog in the regulation of the proneural gene atonal during ommatidia development. Development126, 2345–2353.
25
DunlopM. J.CoxR. S.3rdLevineJ. H.MurrayR. M.ElowitzM. B. (2008). Regulatory activity revealed by dynamic correlations in gene expression noise. Nat. Genet.40, 1493–1498. 10.1038/ng.281
26
EricksonT.FrenchC. R.WaskiewiczA. J. (2010). Meis1 specifies positional information in the retina and tectum to organize the zebrafish visual system. Neural. Dev.5:22. 10.1186/1749-8104-5-22
27
FirthL. C.BakerN. E. (2009). Retinal determination genes as targets and possible effectors of extracellular signals. Dev. Biol.327, 366–375. 10.1016/j.ydbio.2008.12.021
28
GarfieldD. A.RuncieD. E.BabbittC. C.HaygoodR.NielsenW. J.WrayG. A. (2013). The impact of gene expression variation on the robustness and evolvability of a developmental gene regulatory network. PLoS Biol.11:e1001696. 10.1371/journal.pbio.1001696
29
GelbartW. M. (1989). The decapentaplegic gene: a TGF-beta homologue controlling pattern formation in Drosophila. Development107 (Suppl), 65–74.
30
GuiR.LiuQ.YaoY.DengH.MaC.JiaY.et al. (2016). Noise decomposition principle in a coherent feed-forward transcriptional regulatory loop. Front. Physiol.7:600. 10.3389/fphys.2016.00600
31
HalderG.CallaertsP.FlisterS.WalldorfU.KloterU.GehringW. J. (1998). Eyeless initiates the expression of both sine oculis and eyes absent during Drosophila compound eye development. Development125, 2181–2191.
32
HalderG.CallaertsP.GehringW. J. (1995). Induction of ectopic eyes by targeted expression of the eyeless gene in Drosophila. Science267, 1788–1792. 10.1126/science.7892602
33
HeberleinU.TreismanJ. E. (2000). Early retinal development in Drosophila. Results Probl. Cell. Differ.31, 37–50. 10.1007/978-3-540-46826-4_3
34
IkedaK.WatanabeY.OhtoH.KawakamiK. (2002). Molecular interaction and synergistic activation of a promoter by Six, Eya, and Dach proteins mediated through CREB binding protein. Mol. Cell Biol.22, 6759–6766. 10.1128/MCB.22.19.6759-6766.2002
35
JaegerJ.Manu ReinitzJ. (2012). Drosophila blastoderm patterning. Curr. Opin. Genet. Dev.22, 533–541. 10.1016/j.gde.2012.10.005
36
JemcJ.RebayI. (2007). The eyes absent family of phosphotyrosine phosphatases: properties and roles in developmental regulation of transcription. Annu. Rev. Biochem.76, 513–538. 10.1146/annurev.biochem.76.052705.164916
37
KaiserR.PosteguilloE. G.MullerD.JustW. (2007). Exclusion of genes from the EYA-DACH-SIX-PAX pathway as candidates for Branchio-Oculo-Facial syndrome (BOFS). Am. J. Med. Genet. A143A, 2185–2188. 10.1002/ajmg.a.31875
38
KenyonK. L.Yang-ZhouD.CaiC. Q.TranS.ClouserC.DeceneG.et al. (2005). Partner specificity is essential for proper function of the SIX-type homeodomain proteins Sine oculis and Optix during fly eye development. Dev. Biol.286, 158–168. 10.1016/j.ydbio.2005.07.017
39
KumarJ. P. (2009). The molecular circuitry governing retinal determination. Biochim. Biophys. Acta. 1789, 306–31410.1016/j.bbagrm.2008.10.001
40
LevineM.DavidsonE. H. (2005). Gene regulatory networks for development. Proc. Natl. Acad. Sci. U.S.A.102, 4936–4942. 10.1073/pnas.0408031102
41
LiY.JiangY.ChenY.KarandikarU.HoffmanK.ChattopadhyayA.et al. (2013). optix functions as a link between the retinal determination network and the dpp pathway to control morphogenetic furrow progression in Drosophila. Dev. Biol.381, 50–61. 10.1016/j.ydbio.2013.06.015
42
LopesC. S.CasaresF. (2010). hth maintains the pool of eye progenitors and its downregulation by Dpp and Hh couples retinal fate acquisition with cell cycle exit. Dev. Biol.339, 78–88. 10.1016/j.ydbio.2009.12.020
43
LopesC. S.CasaresF. (2015). Eye selector logic for a coordinated cell cycle exit. PLoS Genet.11:e1004981. 10.1371/journal.pgen.1004981
44
ManganS.AlonU. (2003). Structure and function of the feed-forward loop network motif. Proc. Natl. Acad. Sci. U.S.A.100, 11980–11985. 10.1073/pnas.2133841100
45
ManganS.ZaslaverA.AlonU. (2003). The coherent feedforward loop serves as a sign-sensitive delay element in transcription networks. J. Mol. Biol.334, 197–204. 10.1016/j.jmb.2003.09.049
46
MardonG.SolomonN. M.RubinG. M. (1994). Dachshund encodes a nuclear protein required for normal eye and leg development in Drosophila. Development120, 3473–3486.
47
MunskyB.NeuertG.van OudenaardenA. (2012). Using gene expression noise to understand gene regulation. Science336, 183–187. 10.1126/science.1216379
48
Naval-SanchezM.PotierD.HaagenL.SanchezM.MunckS.Van de SandeB.et al. (2013). Comparative motif discovery combined with comparative transcriptomics yields accurate targetome and enhancer predictions. Genome Res.23, 74–88. 10.1101/gr.140426.112
49
NetoM.Aguilar-HidalgoD.CasaresF. (2016). Increased avidity for Dpp/BMP2 maintains the proliferation of progenitors-like cells in the Drosophila eye. Dev. Biol.418, 98–107. 10.1016/j.ydbio.2016.08.004
50
NetoM.Naval-SanchezM.PotierD.PereiraP. S.GeertsD.AertsS.et al. (2017). Nuclear receptors connect progenitor transcription factors to cell cycle control. Sci. Rep.7:4845. 10.1038/s41598-017-04936-7
51
OhtoH.KamadaS.TagoK.TominagaS. I.OzakiH.SatoS.et al. (1999). Cooperation of six and eya in activation of their target genes through nuclear translocation of Eya. Mol. Cell Biol.19, 6815–6824. 10.1128/MCB.19.10.6815
52
OstrinE. J.LiY.HoffmanK.LiuJ.WangK.ZhangL.et al. (2006). Genome-wide identification of direct targets of the Drosophila retinal determination protein Eyeless. Genome Res.16, 466–476. 10.1101/gr.4673006
53
PappuK. S.OstrinE. J.MiddlebrooksB. W.SiliB. T.ChenR.AtkinsM. R.et al. (2005). Dual regulation and redundant function of two eye-specific enhancers of the Drosophila retinal determination gene dachshund. Development132, 2895–2905. 10.1242/dev.01869
54
PichaudF.CasaresF. (2000). homothorax and iroquois-C genes are required for the establishment of territories within the developing eye disc. Mech. Dev.96, 15–25. 10.1016/S0925-4773(00)00372-5
55
PignoniF.HuB.ZavitzK. H.XiaoJ.GarrityP. A.ZipurskyS. L. (1997). The eye-specification proteins So and Eya form a complex and regulate multiple steps in Drosophila eye development. Cell91, 881–891. 10.1016/S0092-8674(00)80480-8
56
PineiroC.LopesC. S.CasaresF. (2014). A conserved transcriptional network regulates lamina development in the Drosophila visual system. Development141, 2838–2847. 10.1242/dev.108670
57
PunzoC.SeimiyaM.FlisterS.GehringW. J.PlazaS. (2002). Differential interactions of eyeless and twin of eyeless with the sine oculis enhancer. Development129, 625–634. Available online at: http://dev.biologists.org/content/129/3/625.long
58
PurcellP.OliverG.MardonG.DonnerA. L.MaasR. L. (2005). Pax6-dependence of Six3, Eya1 and Dach1 expression during lens and nasal placode induction. Gene Expr. Patterns6, 110–118. 10.1016/j.modgep.2005.04.010
59
RoyoJ. L.MaesoI.IrimiaM.GaoF.PeterI. S.LopesC. S.et al. (2011). Transphyletic conservation of developmental regulatory state in animal evolution. Proc. Natl. Acad. Sci. U.S.A.108, 14186–14191. 10.1073/pnas.1109037108
60
RyderE.SpriggsH.DrummondE.St JohnstonD.RussellS. (2009). The Flannotator–a gene and protein expression annotation tool for Drosophila melanogaster. Bioinformatics25, 548–549. 10.1093/bioinformatics/btp012
61
Sanchez-AragonM.CasaresF. (2019). A new image segmentation algorithm with applications in confocal microscopy analysis. BioRxiv. 10.1101/524389
62
SilverS. J.RebayI. (2005). Signaling circuitries in development: insights from the retinal determination gene network. Development132, 3–13. 10.1242/dev.01539
63
SinghA.Kango-SinghM.SunY. H. (2002). Eye suppression, a novel function of teashirt, requires Wingless signaling. Development129, 4271–4280. Available online at: http://dev.biologists.org/content/129/18/4271.long
64
SmithS. J.RebeizM.DavidsonL. (2018). From pattern to process: studies at the interface of gene regulatory networks, morphogenesis, and evolution. Curr. Opin. Genet. Dev.51, 103–110. 10.1016/j.gde.2018.08.004
65
SpielerD.KaffeM.KnaufF.BessaJ.TenaJ. J.GiesertF.et al. (2014). Restless legs syndrome-associated intronic common variant in Meis1 alters enhancer function in the developing telencephalon. Genome Res.24, 592–603. 10.1101/gr.166751.113
66
StruhlG.BaslerK. (1993). Organizing activity of wingless protein in Drosophila. Cell72, 527–540. 10.1016/0092-8674(93)90072-X
67
Tanaka-MatakatsuM.DuW. (2008). Direct control of the proneural gene atonal by retinal determination factors during Drosophila eye development. Dev. Biol.313, 787–801. 10.1016/j.ydbio.2007.11.017
68
TreismanJ. E. (2013). Retinal differentiation in drosophila. Wiley Interdisc. Rev. Dev. Boil.2, 545–557. 10.1002/wdev.100
69
VollmerJ.FriedP.Sanchez-AragonM.LopesC. S.CasaresF.IberD. (2016). A quantitative analysis of growth control in the Drosophila eye disc. Development143, 1482–1490. 10.1242/dev.129775
70
WartlickO.JulicherF.Gonzalez-GaitanM. (2014). Growth control by a moving morphogen gradient during Drosophila eye development. Development141, 1884–1893. 10.1242/dev.105650
71
WeasnerB. M.WeasnerB.DeyoungS. M.MichaelsS. D.KumarJ. P. (2009). Transcriptional activities of the Pax6 gene eyeless regulate tissue specificity of ectopic eye formation in Drosophila. Dev. Biol.334, 492–502. 10.1016/j.ydbio.2009.04.027
72
WiersdorffV.LecuitT.CohenS. M.MlodzikM. (1996). Mad acts downstream of Dpp receptors, revealing a differential requirement for dpp signaling in initiation and propagation of morphogenesis in the Drosophila eye. Development122, 2153–2162.
73
XuT.RubinG. M. (1993). Analysis of genetic mosaics in developing and adult Drosophila tissues. Development117, 1223–1237.
74
ZhangT.RanadeS.CaiC. Q.ClouserC.PignoniF. (2006a). Direct control of neurogenesis by selector factors in the fly eye: regulation of atonal by Ey and So. Development133, 4881–4889. 10.1242/dev.02669
75
ZhangX.FriedmanA.HeaneyS.PurcellP.MaasR. L. (2002). Meis homeoproteins directly regulate Pax6 during vertebrate lens morphogenesis. Genes Dev.16, 2097–2107. 10.1101/gad.1007602
76
ZhangX.RowanS.YueY.HeaneyS.PanY.BrendolanA.et al. (2006b). Pax6 is regulated by Meis and Pbx homeoproteins during pancreatic development. Dev. Biol. 300, 748–757. 10.1016/j.ydbio.2006.06.030
77
ZhouQ.ZhangT.JemcJ. C.ChenY.ChenR.RebayI.et al. (2014). Onset of atonal expression in Drosophila retinal progenitors involves redundant and synergistic contributions of Ey/Pax6 and So binding sites within two distant enhancers. Dev. Biol.386, 152–164. 10.1016/j.ydbio.2013.11.012
78
ZhuJ.PalliyilS.RanC.KumarJ. P. (2017). Drosophila Pax6 promotes development of the entire eye-antennal disc, thereby ensuring proper adult head formation. Proc. Natl. Acad. Sci. U.S.A.114, 5846–5853. 10.1073/pnas.1610614114
Summary
Keywords
eye development, Drosophila, gene regulatory networks, quantitative gene expression, modeling, cell specification, noise
Citation
Sánchez-Aragón M, Cantisán-Gómez J, Luque CM, Brás-Pereira C, Lopes CS, Lemos MC and Casares F (2019) A Toggle-Switch and a Feed-Forward Loop Engage in the Control of the Drosophila Retinal Determination Gene Network. Front. Ecol. Evol. 7:221. doi: 10.3389/fevo.2019.00221
Received
28 March 2019
Accepted
28 May 2019
Published
12 June 2019
Volume
7 - 2019
Edited by
Alistair Peter McGregor, Oxford Brookes University, United Kingdom
Reviewed by
Sebastian Kittelmann, University of Oxford, United Kingdom; Berta Verd, University of Cambridge, United Kingdom
Updates

Check for updates
Copyright
© 2019 Sánchez-Aragón, Cantisán-Gómez, Luque, Brás-Pereira, Lopes, Lemos and Casares.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: M. Carmen Lemos lemos@us.esFernando Casares fcasfer@upo.es
‡Present Address: Julia Cantisán-Gómez, Department of Physics, Universidad Rey Juan Carlos, Madrid, Spain
This article was submitted to Evolutionary Developmental Biology, a section of the journal Frontiers in Ecology and Evolution
†These authors have contributed equally to this work
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.