A dynamic cell recruitment process drives growth of the Drosophila wing by overscaling the Vestigial expression pattern

Organs mainly attain their size by cell growth and proliferation, but sometimes also grow through recruitment of undifferentiated cells. Here we investigate the participation of cell recruitment in establishing the pattern of Vestigial (Vg), the product of the wing selector gene in Drosophila. We find that the Vg pattern overscales along the dorsal-ventral (DV) axis of the wing imaginal disc, i.e., it expands faster than the DV length of the pouch. The overscaling of the Vg pattern cannot be explained by differential proliferation, apoptosis, or oriented-cell divisions, but can be recapitulated by a mathematical model that explicitly considers cell recruitment. By impairing cell recruitment genetically, we find that the Vg pattern almost perfectly scales and adult wings are approximately 20% smaller. Furthermore, using fluorescent reporter tools, we provide direct evidence that cell recruitment takes place in a specific time between early and mid third-instar larval development. Altogether, our work quantitatively shows when, how, and by how much cell recruitment shapes the Vg pattern and drives growth of the Drosophila wing.


Introduction 1 2
Organ growth during development is orchestrated by morphogens, which are signaling molecules 3 that determine gene expression in a non-cell autonomous manner and also act as mitogens (Day and 4 Lawrence, 2000;Schwank and Basler, 2010;Dekanty and Milán, 2011;Lander, 2011;Wartlick et 5 al., 2011;Bryant and Gardiner, 2016;Vollmer et al., 2017). However, growth may also be driven 6 independently of morphogen-induced cell proliferation, by a growth-by-patterning mechanism in 7 which a differentiation pattern expands at the expense of the incorporation of undifferentiated cells. 8 This mechanism, also known as inductive assimilation or cell recruitment, has been reported to 9 work in both vertebrate and invertebrate development such as in the eye (Heberlein et al., 1995;10 Strutt et al., 1995) and wing discs (Baena-López and García-Bellido, 2003;Zecca and Struhl, 11 2007a) of the fruit fly, Drosophila melanogaster, as well as in the mammalian thyroid (Fagman et 12 al., 2006) and kidney (Lindström et al., 2018). Little is known, however, about the dynamic 13 properties of patterning by cell recruitment and how much organ growth can be gained through cell 14 recruitment relative to cell growth and proliferation. These questions are generally difficult to 15 address during normal development because the dynamics of patterning and cell proliferation occur 16 simultaneously and it is challenging to isolate the contribution of each of these mechanisms to organ 17 size. The Drosophila wing imaginal disc is a useful system to tackle this problem since previous 18 studies have partly identified the molecular players of recruitment (Zecca and Struhl, 2007a;Zecca 19 and Struhl, 2010) and genetic tools allow manipulation and quantitative assessment of patterning 20 and growth (Hariharan and Bilder, 2006;Beira and Paro, 2016). 21 In Drosophila, wing fate is determined by the expression of the wing selector gene, vestigial 22 (vg), in the pouch of the wing imaginal disc (Williams et al., 1991). The absence of vg results in 23 loss of wing blade identity (Williams et al., 1991;Williams et al., 1993), whereas ectopic 24 expression of vg in other imaginal discs could induce transformation into wing-like tissue (Williams 25 et al., 1994;Kim et al., 1996;Halder et al., 1998;Baena-López and García-Bellido, 2003). 26 Therefore, adult wing size depends on the size and shape of the Vg pattern established during the 27 larval stage. 28 vg expression is controlled at least by two enhancers, the vg boundary enhancer (vg BE ) that 29 responds to Notch signaling and results in high expression of Vg at the dorsal-ventral (DV) 30 boundary (Irvine and Vogt, 1997), and the vg quadrant enhancer (vg QE ) that requires Wingless (Wg) 31 signaling (Kim et al., 1996;Zecca and Struhl, 2007b) resulting in a gradient of Vg along the DV 32 example, cells near the edge of the pouch may be experiencing higher apoptosis rates or Vg-23 expressing cells may proliferate faster than cells outside of the Vg domain, explaining the 24 overscaling of the Vg pattern. To test these possibilities, we first examined the expression of the 25 pro-apototic marker Caspase 3 (Cas 3) throughout development and found that apoptosis occurs at 26 low frequency and seems to be homogeneous in the wing pouch (Fig. S5); this is consistent with a 27 previous study (Milán et al., 1997). Second, we examined cell proliferation and found that during 28 most of the third instar cell proliferation within the pouch is approximately uniform, except for cells 29 at the DV boundary, confirming the results of previous studies (Fig. S6A;Schwank et al., 2011;30 Wartlick et al., 2011;Mao et al., 2013). We conclude that the overscaling of the Vg pattern cannot 31 be accounted by differences in cell proliferation or apoptosis between cells within and outside the 32 Vg domain.

Mathematical modeling of cell recruitment predicts the overscaling of the Vg pattern 2
In order to investigate if the dynamics of the Vg gradient could be explained by a cell recruitment 3 mechanism, we modeled the distribution of Vg in the wing pouch by means of a multi-scale model 4 (see Supplementary Material for full description). The model combines an ordinary differential 5 equation to account for the rate of change of Vg concentration in each cell ( ) with a 2D Cellular 6 Potts Model (CPM; Graner and Glazier, 1992;Glazier and Graner, 1993) describing the cellular 7 dynamics ( Fig. 2A and Fig. S7A, B). We assumed that cells produce Vg by two mechanisms, while 8 there is only one sink modeled as a linear degradation. The first production term assumes that Vg 9 expression responds to the concentration of a given non-scaling morphogen M ( Fig. 2A and  10 Supplementary Material). The second cellular source of Vg is our quantitative formulation of the 11 recruitment mechanism. In particular, we proposed that this production term depends on the 12 difference in Vg concentration between the actual cell and the average of the concentration of its 13 neighbors by a second order Hill function ( Fig. 2A and Fig. S7C). Hence, this production term 14 becomes relevant when the concentration of the cell is different from the one of its neighbors and 15 negligible when they are similar. As a control, we also considered a model without the recruitment 16 term. To account for disc growth, we explicitly included homogeneous cell proliferation in which 17 the cell cycle parameters were fitted to the average number of cells in each of the four groups 18  Video S1). 2 Finally, since cell divisions in the wing pouch are oriented along the proximal-distal axis 3 (González-Gaitán et al. 1994;Resino et al. 2002;Baena-López et al. 2005;Mao et al. 2013), we 4 used our mathematical model to ask whether the overscaling effect could be enhanced by assuming 5 that cell divisions are oriented radially. However, essentially the same fit was achieved in a model 6 that assumes radially-oriented cell divisions vs. a model where cell divisions are randomly oriented 7 (Fig. S12). We conclude that a mathematical model encoding a mechanism of cell recruitment may 8 explain the overscaling dynamics of the Vg gradient regardless of the orientation of cell divisions. vivo, we predict that overscaling phenotype will be lost in ds>vg RNAi discs. We sorted the discs into 18 four groups according to their DV length (Fig. S13B) and quantified the patterns of Vg as we did 19 for the wild-type discs ( gradient. However, we found that cell proliferation, apoptosis, and the patterns of cell divisions are 1 very similar in control and ds>vg RNAi discs ( Fig. S6 and S14). We conclude that preventing Vg 2 expression in the ds expression domain blocks the recruitment-dependent overexpansion of the Vg 3

pattern. 4 5
A time-dependent cell recruitment process contributes to the Vg pattern during normal development 6 So far, the evidence of a cell recruitment mechanism that patterns the Vg gradient has been indirect. 7 In order to investigate more directly the dynamics of the cell recruitment process during normal 8 development, we used a recently-developed tool of dual-color fluorescent reporters, known as 9 TransTimer (He et al. 2019), to examine the dynamics of the vg QE . The TransTimer expresses a 10 rapid yet unstable GFP and a more stable, but slower RFP under vg QE Gal4 control. Hence, recently-11 recruited cells will express GFP but not RFP (or have a RFP to GFP ratio below a threshold), 12 whereas cells that have been activating the reporters for longer than 12 hours will express RFP as 13 well (Fig. 4A). We aimed to use this system, both to provide direct evidence of the recruitment 14 process and to show that it takes place at a particular time during normal development, as revealed 15 by our time course of Vg expression (Fig. 1G). Again, we used DV length as defined in Fig. 1

Impairment of cell recruitment results in smaller, but well-proportioned adult wings 26
We then asked whether cell recruitment could have an effect on adult wing size. Therefore, we 27 compared ds>vg RNAi vs. control adult wings (Fig. 5). Most adult wings of ds>vg RNAi animals show 28 the normal vein patterns (Fig. 5A-B). However, ds>vg RNAi wings are on average smaller than control 29 wings (Fig. 5C). We assume that cell recruitment takes place in dorsal-and ventral-most regions of 30 the wing pouch, which correspond to proximal regions of the adult wing. Therefore, we predicted 31 that ds>vg RNAi adult wings would be smaller in proximal regions, but unaffected in distal regions. 1 However, we found that representative areas of both proximal (Fig. 5D, inset) and distal (Fig. 5E, 2 inset) regions of the adult wing are similarly reduced in ds>vg RNAi animals (17 and 21 %, 3 respectively; Fig. 5D-E). Furthermore, we found that ds>vg RNAi wings maintain their proximal-distal 4 (longitudinal) to anterior-posterior (transversal) proportions, as the ratio of longitudinal to 5 transversal dimensions is not statistically different between control and mutant adult wings (Fig.  6   5F). Taken together, we conclude that the impairment of cell recruitment in ds>vg RNAi animals 7 significantly affects adult wing size but not its proportions. is uniform, this results in patterns that scale with size throughout development. In both of these 20 models, the patterning process adjust to organ growth and the final patterns are invariant to organ 21 size. Our data on patterning of the wing selector gene, vg, in the Drosophila wing disc support an 22 additional conceptual model. After being established by Notch and Wg signaling, Vg expression is 23 maintained through cell divisions by a mechanism involving Polycomb/Trithorax Responsive 24 Elements. In a tissue that grows uniformly like the wing disc, this would result in a pattern that 25 scales with tissue size. However, we find that the Vg gradient overscales relative to DV axis length 26 ( Fig. 1). This cannot be explained by non-uniform cell proliferation, apoptosis, or a substantial 27 change in pattern of oriented cell divisions (Figs. S5, S6, S14). Instead, our modeling results 28 suggest that Vg overscaling could be explained by a cell recruitment mechanism (Fig. 2). 29 Experimentally, when cell recruitment is genetically blocked, the overscaling phenotype is mostly 30 lost (Fig. 3). We provide direct evidence of the cell recruitment dynamics during Vg patterning 31 (Fig. 4). Finally, we show that cell recruitment contributes to about 20 % of adult wing size (Fig. 5). 32 Based on qualitative evidence on Vg expansion in genetic mosaics, Zecca and Struhl 1 identified the molecular mechanisms of recruitment and proposed that Vg-dependent cell 2 recruitment could contribute to growth of the Drosophila wing disc (Zecca and Struhl, 2007a; Zecca 3 and Struhl 2010). However, it remained unclear when and to what extent cell recruitment actually 4 contributes to normal patterning and growth of the disc. Our study directly shows that cell 5 recruitment works by overexpanding the endogenous pattern of a selector gene, thus contributing to 6 both patterning and organ growth. To our knowledge, this is the first study where recruited cells are 7 directly detected during wild-type development and where its contributions to pattern formation and 8 organ growth are quantitatively determined. 9 Wg signaling is absolutely required for Vg expression and propagation along the DV axis 10 In this study, we focused on Vg patterning along the DV axis, but other signaling pathways 21 participate in Vg expression. For example, Decapentaplegic (Dpp), a member of the BMP family, 22 patterns the disc along the anterior-posterior (AP) axis and has been proposed an input signaling for 23 vg expression (Kim et al., 1996). Importantly, our analysis using the TransTimer system does shows 24 new activation of the GFP reporter in cells located in the anterior and posterior ends of the disc 25 (arrowheads in Fig. 4B''). Therefore, how Dpp signaling contributes to cell recruitment along the 26 AP axis and participates in the 2-dimensional Vg pattern is a promising perspective left to future 27 work. 28 Our TransTimer experiment shows that new cells are being incorporated into the vg QE 29 pattern, i.e. new cells get recruited, and reveal that this recruitment process is temporally controlled 30 (Fig. 4). This is consistent with the major overscaling of the Vg width that occurs in discs of similar 31 sizes in Fig. 1G. However, it is likely that the TransTimer reporters only show the earliest 32 manifestation of the recruitment process and that this continues over the rest of the third larval 33 instar by increasing the levels of Vg until they reach an appropriate level for wing-fate 1 differentiation. Thus, we interpret the dynamics of Vg pattern formation as a two-step recruitment 2 process: The first step takes place between early and mid third-instar, and results in the overscaling 3 of the Vg gradient width (Fig. 1G,H). The fact that this step occurs during a narrow window of time 4 (Fig. 4), suggests that the Ft-Ds polarization signal can propagate several cells without the need of a 5 cell-by-cell expansion of Vg expression, i.e. several layers of cells could begin to be recruited at 6 once (Wortman et al. 2017). The second step extends from mid to late third-larval instar and results 7 in increasing levels of Vg at the tails of the gradient. This is captured by a decrease in the slopes of 8 the Vg gradient in our time-course examination of wild-type Vg patterns ( Fig. 1G and Fig. S4). The 9 fact that the slopes at the tails of the Vg gradient never completely flatten suggests that the rates of 10 Vg upregulation are position-dependent, possibly dictated by the strength of the Ft-Ds polarization 11 which decays from its source near the DV boundary towards ventral-and dorsal-most positions. 12 This two-step model of cell recruitment suggests that the spatial range of the process is limited by 13 the length-scale of the Ft-Ds polarization signal or by the availability of recruitable cells within the 14 wing pouch. 15 We show that cell recruitment contributes to approximately one fifth of the total adult wing 16 size (Fig. 5C). An interesting question that remains is to understand the purpose of cell recruitment 17 as a growth mechanism. In other words, what is the advantage of this developmental design over 18 simply adjusting cell sizes or cell proliferation rates to achieve a target final size? It is possible that 19 cell recruitment plays a role in conferring some sort of robustness to developmental growth control. 20 For example, perhaps recruitment compensates for variations in cell proliferation and could explain 21 why final wing disc size is robust to perturbations against cell proliferation rates or cell size (Day 22 and Lawrence, 2000). 23 Recent studies have showed that a cell recruitment mechanism is present in early thyroid New York, USA); vg QE -Gal4 (III) (BDSC #8229). Imaginal discs were dissected from third-instar 7 larvae of both sexes. For normal endogenous Vg patterns quantification, the y, w stock was used as 8 wild-type flies (Fig. 1) and kept at 25°C during egg laid and development. For recruitment-9 impairment experiments (Figs. 3 and 5), dsGal4, UAS-GFP flies were crossed to UAS-vg RNAi (for 10 recruitment inhibition) or the y, w stock (control) and were kept at 25°C during egg laid and then 11 changed to 29°C during development to increase Gal4 system efficiency. For the evaluation of the 12 vg QE dynamics (Fig. 4), UAS-TransTimer reporter, a destabilized GFP (dGFP) protein combined 13 with a stable RFP (He et al. 2019), was placed under vg QE Gal4 driver to quantify cells that recently 14 activated this enhancer. Taking advantage of the different maturation rates of both proteins (dGFP 15 about 0.1 h, and RFP about 1.5 h), we selected all the GFP pixels with intensity values above a 16 threshold in the wing pouch for each disc; this threshold was set by ranking all the pixels within by 17 their GFP intensity in three different groups using the machine learning K-Means function from 18 scikit-learn python package (https://scikit-learn.org/stable/index.html) and all the pixels in the 19 lower group were taken away. Next, we looked for the corresponding RFP channel and computed a 20 ratio of RFP over GFP intensity. If this ratio was lower than 0.4 at a specific pixel (this value 21 corresponds to the half maximum for GFP taken from was used to stain nuclei. 5-ethynyl-2'-deoxyuridine (EdU) labeling was performed using the Click-32 iT EdU Alexa Fluor 647 Imaging Kit (Invitrogen Cat#C10340) following manufacturer instructions. 1 Primary antibodies were detected with Alexa Fluor 488 anti-mouse and 647 anti-rabbit (1:1000). 2 Imaging was done with a confocal microscope (Leica TCS SP8 Confocal Microscope) using a 63X 3 oil-immersion objective. 4 5

Wing mounting 6
Adult flies were dehydrated overnight in 70% ethanol and then separated by gender. Wings were 7 dissected in 50% ethanol. The isolated wings were mounted and dried in a plate at 60 °C. Imaging 8 of adult wings was done using a bright-field microscope (Nikon eclipse Ci-L/Ci-S) using a 4X 9 objective. Wings were dissected from female and male flies and were independently analyzed. 10 11

Wing disc and adult wing quantification 21
Vg patterns in the wing disc were quantified as explained in Fig. S1. Wing representative distances 22 and areas were quantified using the straight line and polygon selection tools respectively in 23 ImageJ/Fiji software (https://imagej.net/). 24

25
Simulations 26 All simulations were performed using Morpheus 1.9.3 software (https://imc.zih.tu-27 dresden.de/wiki/morpheus). The time dependence of the Vg concentration was modeled by means 28 of an ordinary differential equation for each cell ( Fig. 2A), and the tissue was modeled by a Cellular 29 Pots Model (CPM). Examples of the .xml files used to perform simulations are provided as 1 supplementary files. 2 3 Image processing and data visualization 4 All the images were processed and analyzed using ImageJ/Fiji software, the matplotlib 5 (https://matplotlib.org/), pandas (https://pandas.pydata.org/), and NumPy (http://www.numpy.org/) 6 python packages. The data, after being analyzed, were visualized with the seaborn 7 (https://seaborn.pydata.org/) and Matplotlib python packages. 8 9 Acknowledgments 10 11 We thank Li He from Norbert Perrimon's Lab for kindly providing us the UAS-TransTimer flies 12 and Fanis Missirlis for critical comments on the manuscript. We also thank the members of the 13 Nahmad and Chara laboratories for interesting discussions.  pouch where the quantification of the Vg pattern was performed (yellow rectangle). This 3 rectangular region is centered in the anterior posterior border of the disc and is delimited by the 4 ventral and dorsal folds that separate the pouch from the hinge (throughout the paper, we use the 5 distance between these folds as a measure of disc size). The intensity of Vg in each vertical line of 6 pixels within the rectangle is averaged to obtain a spatial Vg expression along the DV axis. (B) 7 Third-instar (ages 85-120 h After Egg Laying (AEL) at 25ºC) yellow-white (y, w) discs (throughout 8 the study, we refer to y, w discs as wild type) are classified into four groups according to the 9 distance between dorsal and ventral folds (DV length). The groups are defined by dividing the 10 shortest to the longest DV length into four intervals of equal length. (C-F) Representative images of 11 discs within each group immunostained with a Vg antibody (C-F) and DAPI (C'-F') and the 12 quantification of the Vg patterns in the region defined in A in absolute (C''-F'') and relative (C'''-13 the Pearson Correlation Coefficient, r, of the regression is displayed. The p-value corresponds to a 23 Student-t statistical test assuming a zero-slope of the regression line as a null hypothesis. (Note that 24 a positive slope corresponds to overscaling; a zero-slope, i.e., when the null hypothesis cannot be 25 rejected, corresponds to perfect scaling; whereas a negative slope corresponds to underscaling). (I) 26 Comparison of the location of the Vg gradient peak in RU for the discs in each group; error bars 27 correspond to the SEM. (Statistical significance is analyzed using a Kruskal-Wallis non-parametric 28 test; ***, p-value<10 -5 . Pairwise statistical comparisons between groups 1 and 2, and groups 3 and 4 29 using a Mann-Whitney non-parametric test reveal non-statistical differences between these groups; 30 ns, p-value<0.01).