ORIGINAL RESEARCH article

Front. Ecol. Evol., 12 June 2019

Sec. Evolutionary Developmental Biology

Volume 7 - 2019 | https://doi.org/10.3389/fevo.2019.00221

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

  • 1. GEM-DMC2 Unit, CABD (CSIC-Universidad Pablo de Olavide-Junta de Andalucía), Seville, Spain

  • 2. Department of Condensed Matter Physics, University of Sevilla, Sevilla, Spain

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 3ato-destabilized GFP (3ato-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

). When the segmented stacks are re-sliced (B–B″), nuclear boundaries (RED) appear correctly located in most nuclei. The insets show higher magnification views of the boxed regions. (C,D) 3D reconstruction of an early disc showing a sampling volume (white dashed box) oriented parallel to the AP axis and centered at the morphogenetic furrow, stained with Hth, Eya and DAPI. (x,y) (A) and (x,z) (B) views are shown. The double-headed arrow indicates the region anterior to the MF. In general, the 3D reconstruction may be sampled in specific directions (white dashed box in A,B), and the profiles for specific genes measured with high accuracy. The result is similar to a FACS analysis, but with the advantage of keeping the positional information.

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., ) 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 (Mangan and Alon, ):

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., ; Sanchez-Aragon and Casares, ). 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.

Figure 2

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 [3ato-dGFP and optix2/3-dGFP (Ostrin et al., )] (see Materials and Methods), which constitute a dense sampling of many key nodes of the gene network.

Figure 3

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., ). 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., ). 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., ; Singh et al., ), and Ey, Hth, and Tsh might form a protein complex (Bessa et al., ; Datta et al., ). 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., ). 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 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., ). 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., ; Ohto et al., ; Jemc and Rebay, ), itself a Pax6 target (Halder et al., ; Punzo et al., ). 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., , ). 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., ). 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., ).

Figure 4

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, ). 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, ; Heberlein and Treisman, ). 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, ; Kumar, ; Amore and Casares, ; Casares and Almudi, ), 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, ). 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 (dac3) 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. (), showing that in dac-mutant clones Ey expression is derepressed. Our interpretation of their results is that dac-clones cause a delay in MF movement (Mardon et al., ; Bras-Pereira et al., ) 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., ). 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, ; 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., ) and string [stg; (Lopes and Casares, )] is similar, being activated by both Ey and Eya:So. The mutual Hth/Eya:So repression is well established (Bessa et al., ; Lopes and Casares, ) 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., ; Lopes and Casares, ), and that Hth hinders Ey's capacity to activate stg (Lopes and Casares, ). 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).

Figure 5

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 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 (Mangan et al., ). 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, ). Dpp represses Hth at long distance (Bessa et al., ) and contributes to the repression of Ey at shorter range, just abutting the MF (Firth and Baker, ). 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, ).

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., ). As hth is repressed, a coherent feed-forward 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 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, ; Bessa et al., ; Neto et al., ), 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., ) 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., ). 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., ). 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., ).

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

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

Copyright

*Correspondence: M. Carmen Lemos Fernando Casares

‡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.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics