The wing disc is composed of multiple cell types, including the columnar cells of the disc proper, the squamous cells of the PE, and the mesoderm-derived AMPs (Figure 1A). In addition, the wing disc is in intimate contact with branches of the tracheal system and circulating blood cells called hemocytes. With the goal of generating a spatiotemporal atlas of the developing wing disc, we used single-cell RNA sequencing to collect transcriptional profiles of cells at mid and late third instar, which correspond to 96 and 120 hr after egg lay (AEL) (Figure 1A). Two biological replicates were obtained at each time point that, after filtering for low-quality cells, generated data from 6922 and 7091 cells in the 96 hr samples and 7453 and 5550 cells in the 120 hr samples. Harmonization of the different samples and dimensionality reduction was performed using scVI (Lopez et al., 2018). Clustering and differential expression analysis was done with the Seurat v3 R package (Stuart et al., 2019), and two-dimensional visualization of the data was performed with Uniform Manifold Approximation and Projection (UMAP) (McInnes et al., 2018) (see Materials and methods).

Our single-cell analysis identified four major cell types via known gene markers (Figure 1B–E, Figure 1—figure supplement 1A–G): the AMPs, the wing disc epithelial cells (disc proper and PE), and small numbers of tracheal cells and hemocytes. We observed expression of a similar set of marker genes in each of these cell types to those described by others (Deng et al., 2019; Bageritz et al., 2019; Zappia et al., 2020). Altogether, we recovered profiles for 19,885 AMPs, 7104 wing disc epithelial cells, 15 tracheal cells, and 12 hemocytes. Notably, our dataset shows an overrepresentation of AMPs, probably the result of our collagenase-based dissociation protocol (see Materials and methods) that dissociated AMPs far more effectively than epithelial cells. This enrichment of AMPs has enabled an especially detailed analysis of this cell type, which has previously received less attention. Unlike bulk RNA-sequencing approaches that would average changes in gene expression across multiple cell types, we were able to observe expression changes between 96 and 120 hr that occurred in both the epithelium and the AMPs together, as well as those changes that were confined to either the epithelial cells or the AMPs (Figure 1—figure supplement 1H and Supplementary file 1).

To search for the signals that might be exchanged between the epithelium and the AMPs, we first characterized the cell types within each of these populations separately. The wing disc epithelium is often divided into four broad domains – the notum, hinge, pouch, and PE – based both on morphology and gene expression patterns. The genes encoding the transcription factors nubbin (nub) and Zn finger homeodomain 2 (zfh2) were used to define the pouch and hinge, respectively (Zirin and Mann, 2007; Terriente et al., 2008; Ayala-Camargo et al., 2013). Expression of these proteins is shown for wing discs at 96 and 120 hr AEL (Figure 1F, G). To characterize the epithelial compartment, we clustered the cells and then classified these clusters as originating from the larger domains of notum, hinge, pouch, and PE based on the expression patterns of marker genes (Figure 1H–K, Figure 1—figure supplement 2A). Many of these genes were already expressed in a domain-specific manner by 96 hr, suggesting that cells in the epithelium had already been partitioned into these domains (Figure 1I–K, Figure 1—figure supplement 2B).

The proximodistal axis was a primary feature in stratifying cells within our analysis, while the anteroposterior axis separated the data to a lesser degree. This is somewhat surprising because anterior and posterior cells arise from distinct embryonic subpopulations (Garcia-Bellido et al., 1973; Mandaravally Madhavan and Schneiderman, 1977; Worley et al., 2013; Requena et al., 2017). We classified cells as anterior or posterior cells based on expression of compartment marker genes (e.g., cubitus interruptus, hedgehog, and engrailed) (Figure 1—figure supplement 3A–E) and found that the concise representation of the data in two dimensions (with UMAP) resulted in stratification of the cells based on the proximodistal axis first, with secondary stratification of cells based on anteroposterior identity (Figure 1—figure supplement 3E). Additional analysis indicated that there are more differentially expressed genes across the proximodistal axis than the anteroposterior axis (Figure 1—figure supplement 3F–H).

To investigate how cells changed over developmental time, we analyzed transcriptional changes that occurred within epithelium cell clusters between 96 and 120 hr AEL and found that 337 and 408 genes were significantly downregulated and upregulated, respectively, within at least one cluster (see Materials and methods; Figure 1—figure supplement 2C; Supplementary file 2). One example is string (stg) (Edgar and O'Farrell, 1990), which encodes a regulator of the cell cycle and is upregulated within the wing margin while being downregulated in other regions of the disc. Thus, even though the major cell types are established by 96 hr of development (Figure 1I–K, Figure 1—figure supplement 2B), we still find evidence of further pattern refinement via highly localized gene expression changes (Figure 1—figure supplement 2C).

Initial analysis of the AMPs showed a clear partition of the cells with respect to two primary features: cell sex (discs were collected from both male and female larvae) and cell cycle phase (Figure 1—figure supplement 4A–G, Figure 1—figure supplement 5A–G). We utilized scVI to suppress the effects of these covariates, which enabled us to obtain a clearer view of other aspects of AMP cell biology (see Materials and methods).

The AMPs are known to differentiate into either DFMs or IFMs of the adult fly (Bate, 1993; Roy and VijayRaghavan, 1999; Sudarsan et al., 2001). The precursors of these two populations can be identified by their location within the tissue and are canonically classified by their relative expression of two transcription factors, Vg and Ct, at the late third instar larval (L3) stage (corresponding to our 120 hr time point) (Sudarsan et al., 2001). The precursor cells of the IFMs are localized more dorsally (closer to the wing disc stalk) and display relatively high Vg and low Ct protein expression (Sudarsan et al., 2001). The DFM cell precursors are localized more ventrally (closer to the wing hinge) and are identifiable by little or no Vg and high Ct protein expression (Sudarsan et al., 2001). Ct expression at both time points is shown in Figure 1L, M. After unsupervised clustering (Figure 1N), we found that one cluster was characterized by high levels of ct and low levels of vg, and the remaining clusters displayed relatively elevated levels of vg and low levels of ct (Figure 1O–Q, Figure 1—figure supplement 5E, F). Based on this distinction, we classified cells as representing the direct or indirect AMPs (Figure 1N, Figure 1—figure supplement 5F), obtaining 17,604 indirect AMPs and 2281 direct AMPs.

Although the expression of ct and vg was less distinct at 96 hr, our analysis still classified cells from 96 hr into both direct and indirect populations. This prompted us to investigate if there were more subtle differences besides ct and vg at 96 hr that distinguished the direct and indirect AMPs. We identified a small number of genes that at both time points showed differential expression between direct and indirect AMPs (Figure 1—figure supplement 6A; Supplementary file 3), possibly indicative of pathways that are necessary both for the initial establishment and subsequent maintenance of AMP cell types. This included several predicted targets of Wg signaling, specifically naked cuticle (nkd) and vg, and the modulator of Notch signaling fringe (fng) (Figure 1—figure supplement 6A, B). However, the majority of differentially expressed genes that we identified were developmentally regulated. For these genes, differential expression between the direct and indirect AMPs was only observed at 120 hr. Thus, unlike the epithelium, the subdivision of the AMPs into direct and indirect pools (as assessed by differential expression of genes) mostly occurs only by 120 hr.

From 96 to 120 hr, we noticed the activation and refinement of expression of many genes previously implicated in axon guidance (Figure 1—figure supplement 6A). These include the receptor–ligand pairs roundabout2/slit and Neurotactin/Amalgam (Kidd et al., 1999; Frémion et al., 2000), as well as the synaptic partner-matching genes Ten-a and Ten-m (Hong et al., 2012). We confirmed that Ten-m protein increases in expression from 96 to 120 hr by staining discs with anti-Ten-m antibody (Figure 1—figure supplement 6C). These observations suggest that pathways known to function in axon guidance could also function in myoblasts.

Since patterning of the epithelium precedes that of the AMPs, we looked for signals from the epithelium to the AMPs. First, to predict and visualize the location of gene expression within the wing disc, we generated a three-layered model, with gene expression levels inferred at different spatial positions within the AMP, disc proper, and PE layers. This enabled us to discover the expression of genes in regions of the disc epithelium that are closest to the AMPs as well as in the AMPs themselves. Second, we examined the expression of ligand–receptor pairs to identify those with complementary expression patterns between the disc epithelium and the AMPs.

To generate a three-layered transcriptomic map of the wing disc, we mapped our single-cell data to a reference model. The cells were mapped with the R package DistMap (Karaiskos et al., 2017), and the reference model was assembled from a manually curated set of gene expression patterns (Figure 2A and Materials and methods). This spatial mapping of cells was largely consistent with our manual cluster annotations based on marker genes (Figure 2—figure supplement 1). This included the four broad epithelial domains (notum, hinge, pouch, PE; Figure 2B) and more refined cell groups, such as the wing margin, posterior notum, outer pouch, and anterior hinge (Figure 2C). To test how well our virtual wing disc would predict novel gene expression patterns, we predicted the virtual in situ expression patterns of grain (grn) and pou domain motif 3 (pdm3), neither of which were included in generating the disc model. The predicted patterns largely matched the expressions of transcriptional reporters for grn and pdm3 within wing discs (Figure 2D–G), indicating that our virtual wing disc can successfully predict novel gene expression patterns.

Figure 2 with 1 supplement see all

Download asset Open asset

To look for potential cell communication between the different cell layers of the wing disc, we examined the expression of ligand–receptor pairs within the major domains of the disc epithelium and AMPs (Figure 2H) (for details on receptor–ligand pairs examined, see Materials and methods). Interestingly, we observed high levels of expression of the genes encoding two FGF-family ligands, thisbe (ths) and pyramus (pyr) (Stathopoulos et al., 2004), in the notum region of the epithelium, whereas the gene encoding their receptor, heartless (htl) (Beiman et al., 1996), was specifically expressed in the AMPs. Similarly, the ligand hedgehog (hh) appears to be expressed only in the disc epithelium, while its receptor patched (ptc) and signal transducer smoothened (smo) are both expressed in the epithelium as expected, but also unexpectedly in the AMPs. Our detailed investigations of the FGF and Hedgehog pathways are presented in this study.

Ths and Pyr are the ligands for one of the Drosophila FGF signaling pathways, and both interact with the receptor Htl (Stathopoulos et al., 2004; Figure 3A). While htl is detected in nearly all of the AMPs, the three-layered disc map predicts that ths and pyr are expressed primarily in the epithelial cells of the notum, with considerable overlap (Figure 3B–D). Specifically, ths is localized to the most proximal region of the notum, while pyr has a broader expression pattern extending into the posterior hinge (compare Figure 3C and D). The expression patterns of ths and htl reporters were consistent with the expression predicted using our virtual wing disc (Figure 3E, I).

Figure 3

Download asset Open asset

During embryogenesis, htl, ths, and pyr are known to influence mesoderm spreading along the embryonic ectoderm and formation of cardiac progenitor cells (Beiman et al., 1996; Stathopoulos et al., 2004; Kadam et al., 2009). Notably, in htl mutants and in ths and pyr double mutants, mesoderm cells are still present within the embryo, but they accumulate in multilayered arrangements instead of a monolayer along the ectoderm (as observed in wild-type embryos) (Beiman et al., 1996; Stathopoulos et al., 2004; Kadam et al., 2009). These observations suggest that FGF signaling might be primarily needed for proper mesoderm spreading rather than for cell proliferation and survival. By analogy, Ths and Pyr may regulate the localization of AMPs relative to the epithelium.

To examine the consequences of interfering with FGF signaling within the larval wing disc, we perturbed the expression of Pyr and Htl. To disrupt Pyr expression within the epithelial tissue, we utilized an apterous (ap) Gal4 driver (ap-Gal4) that expresses in the entire dorsal compartment of the disc proper, including all of the epithelial cells that overlie the AMPs (Figure 3F, G). Expressing an RNAi that targets pyr with this ap-Gal4 driver resulted in a reduction in AMPs and an obvious increase in apoptosis (visualized using anti-Dcp1), primarily observed in the more ventral- and posterior-localized AMPs (Figure 3H, compare with Figure 3F). A likely explanation for this result is that pyr knockdown within the wing disc restricts AMP survival to the Ths-expressing region of the dorsal notum as Pyr and Ths have been noted to have partially redundant functions (Stathopoulos et al., 2004; Kadam et al., 2009). Moreover, this result suggests that sufficient Ths cannot reach this region of the notum to compensate for the knockdown of Pyr at physiological levels of Ths expression. We next tested if the Ths and Pyr receptor Htl is required within the AMPs. Using the AMP-specific driver 15B03-Gal4 (Figure 3J), we expressed an RNAi for htl and observed an obvious decrease in the number of AMPs following knockdown of FGF signal transduction (Figure 3K). Altogether, we conclude that FGF signaling between AMPs and the disc epithelium is necessary for proper AMP survival.

To determine if the location and level of FGF signaling controls the position and number of the AMPs, we assessed the effects of increasing the levels of FGF ligands and also expressing them ectopically. First, we used a dpp-Gal4 driver that is expressed in a stripe of cells just anterior to the anterior-posterior compartment boundary of the epithelium, including in the notum (Figure 4A). Expression of either pyr or ths in this domain caused a massive increase in the number of AMPs, not only beneath the notum epithelium where AMPs are normally present, but throughout the entire dpp-Gal4 expression domain, including underneath the wing pouch (Figure 4B, Figure 4—figure supplement 1A). AMPs were observed along the entire stripe of ectopic FGF expression, all the way to the ventral hinge and even extending on the ventral side to the peripodial epithelium. From this, we conclude that FGF signaling does not just increase AMP number, but can also induce AMP spreading beyond the epithelial notum.

Figure 4 with 1 supplement see all

Download asset Open asset

We further investigated the role of FGF signaling in AMP migration by generating a separate patch of pyr expression, discontinuous from the domain of endogenous expression. To this end, we ectopically expressed pyr within the wing pouch using the TRiP-Overexpression VPR toolkit (Lin et al., 2015) with a nub driver. At the onset of nub expression, the pouch and the notum are separated by multiple cell diameters. We observed a large number of AMPs basal to the epithelium of the wing pouch (Figure 4C, D), suggesting that ectopic Pyr expression can even recruit AMPs from a distance of several cell diameters.

To ascertain whether the AMPs were capable of proliferation at ectopic locations, we assessed whether the cells were progressing through the cell cycle. AMPs throughout the A-P axis of the disc were found to be in both S- and M-phase of the cell cycle, as indicated by the incorporation of the thymidine analog EdU and staining for phospho-histone-H3 (PHH3) (Figure 4E, F, Figure 4—figure supplement 1D, E). This shows that when an FGF source is provided in the epithelium AMPs are capable of cell cycle entry and progression even in portions of the disc that are distant from the notum, consistent with the large increase in the number of AMPs. Thus, expression of Pyr and Ths in epithelial cells appears sufficient to both attract AMPs to that location and sustain their local proliferation.

To determine whether AMPs can be induced to migrate at different stages of larval development, we restricted the ectopic expression of FGF ligands to the later stages of larval development. The expression of dpp>pyr or dpp>ths was controlled by a temperature-sensitive Gal80, and when expression was initiated in mid-L3 (48 hr prior to dissection), AMPs were observed throughout the dorsal portion of the wing pouch (Figure 4G, Figure 4—figure supplement 1B). In contrast, when expression was initiated later in L3 (24 hr prior to dissection) AMPs were observed in the dorsal hinge but not the pouch (Figure 4H, Figure 4—figure supplement 1C). Taken together, these experiments indicate that the myoblasts found to be associated with other regions of the disc are derived from those underlying the notum. Additionally, these experiments show that AMPs can be induced to emigrate from the notum region even in the later stages of L3, suggesting that the localized expression of the FGF ligands in the notum both localizes AMPs to that region of the disc and also sustains their survival and proliferation (Figure 4I). Thus, FGF signaling effectively defines the AMP niche.

Our analysis indicated that patched (ptc), which encodes the transmembrane receptor for the ligand Hh and is also a transcriptional target for Hh signaling, is expressed at low levels in most AMPs and at a much higher level in a subset of the direct AMPs (Figure 5A). Moreover, Hh-signaling pathway components smoothened (smo) and cubitus interruptus (ci) are expressed in most AMPs at uniform levels (Figure 5B, C). However, we detected negligible levels of hh transcripts within our AMP data, and furthermore did not detect expression of hh-Gal4 within the AMPs (Figure 5D, Figure 5—figure supplement 1B). In contrast, hh transcripts were detected in approximately 32% of cells in the epithelium, roughly the size of the Hh-producing posterior compartment (Figure 1—figure supplement 3B). Together, these observations support the possibility that Hh from posterior cells of the disc epithelium activates Hh signaling in a subset of myoblasts, which respond by expressing higher levels of ptc.

Figure 5 with 2 supplements see all

Download asset Open asset

We detected Ptc protein in a subpopulation of AMPs, localized primarily beneath the posterior compartment of the disc epithelium and extending approximately 20–40 µm into the region underlying the anterior compartment (Figure 5E, Figure 5—figure supplement 1). Consistent with our scRNAseq data, Ptc was observed mostly in a group of direct AMPs, but also in two additional smaller groups of cells that are located more dorsally among indirect AMPs. The proximity of these Ptc-expressing AMPs to Hh-secreting epithelial cells suggests that they are responding to the Hh secreted by these cells rather than a circulating pool of Hh that should be available to all AMPs. Indeed, a recent study suggests that Hh ligand from the epithelium is transported via cytonemes to nearby AMPs (Hatori and Kornberg, 2020).

To determine if Ptc expression in the AMPs closest to the Hh-producing epithelial cells is a result of Hh pathway activation, we reduced smo expression in all AMPs and observed that Ptc expression in AMPs was abolished (Figure 5F). This indicated that, as in the disc epithelium, Hh signal transduction within the AMPs is required to establish high Ptc expression in the posterior-localized AMPs. To address whether all AMPs are capable of this response, we expressed an activated form of the transcription factor Ci (Ci^3m), which is resistant to proteolytic cleavage (Price and Kalderon, 1999), in all AMPs and observed that this leads to an uniformly high level of Ptc protein expression (Figure 5G). Thus, all AMPs appear capable of responding to Hh, but during normal development, only AMPs with close proximity to the posterior compartment of the disc epithelium receive the signal. Neither smo^RNAi-knockdown nor ci^3m overexpression caused obvious changes in AMP numbers (compare Figure 5E with Figure 5F, G), suggesting that Hh signaling is not controlling AMP proliferation but instead is likely important for patterning.

To investigate a possible role of Hh signaling in AMP cell fate specification, we examined adult flies after genetic perturbations for flight muscle defects. During the pupal phase, the AMPs give rise to three distinct muscle fiber types within the adult thorax: dorsal longitudinal muscles (DLM), dorsoventral muscles (DVM), and DFMs (Figure 5H). While DLMs and DVMs are IFMs that generate the mechanical movement required for flight by compressing the thorax, the DFMs are responsible for flight steering by fine-tuning the position of the wing blades (reviewed by Bate, 1993). Both DLMs and DVMs are formed from indirect AMPs, whereas the direct AMPs develop into the DFMs. Control adults displayed wild-type posture (Figure 5I), while after Hh signaling was downregulated in AMPs with smo^RNAi, we found that a majority of adults displayed an ‘outstretched’ wing posture phenotype (Figure 5J, L). When ci^3m was expressed in all AMPs, we observed that many adults displayed a ‘downtilted’ wing posture (Figure 5K, L). These wing posture phenotypes were reproducible with multiple smo^RNAi and in both sexes (Figure 5—figure supplement 2A). Adults with either the outstretched or downtilted phenotypes were incapable of flight. These observations suggest a crucial role of Hh signaling within the AMPs for the formation of functional adult flight muscles.

To examine if Hh-signaling perturbations affected the structure of adult muscle fibers, we dissected adult thoraxes (Figure 5M–R). When we reduced smo expression in the AMPs, we observed a misalignment of the DFM fibers (Figure 5U). In particular, the more posterior DFMs 52–57 (Miller, 1950; Bate, 1993; Ghazi et al., 2000) displayed improper position and overall disorganization (Figure 5—figure supplement 2B–M). Muscle 53, for example, inappropriately projects to the dorsal attachment site of muscle 54. In contrast, the DLM and DVM muscle fibers appeared relatively normal (Figure 5S, T). This indicates that the loss of Hh signaling within the AMPs causes defects in the adult muscles, specifically a subset of muscles formed by the direct AMPs.

Conversely, when the Hh-signaling pathway was constitutively active via the expression of ci^3m within the AMPs, we observed elimination of the DVMs (Figure 5W) and the DFMs were often severely disorganized and malformed (Figure 5X). Importantly, muscle 51, which is derived from a separate group of AMPs not associated with the wing disc (Lawrence, 1982), is unaffected by these manipulations. The DLMs had no noticeable defects (Figure 5V), likely because the DLMs do not arise in de novo, unlike the other adult flight muscles, but rather by the fusion of AMPs with histolyzing larval muscles that act as templates (Fernandes et al., 1991). Overall, our data suggest that Hh signaling is important for proper specification of a subset of the direct AMPs and that excessive Ci activity causes inappropriate patterning that perturbs the development of the IFMs.

What are the downstream targets that are activated by Hh signaling in the AMPs? Since the canonical target of Hh signaling in the disc epithelium, dpp, is not expressed in the AMPs, we searched for candidate genes that were specifically expressed within the subpopulation of AMPs that express high levels of ptc. We found that midline (mid) and Neurotactin (Nrt) displayed relatively high correlation with that of ptc (Pearson correlation of 0.44 and 0.34 with mid and Nrt, respectively) (Figure 6A–C). Mid, also known as Neuromancer 2, is a T-box transcription factor most related to mouse Tbx-20 (Buescher et al., 2004) that regulates cell fate in the developing nervous system (Leal et al., 2009). Nrt encodes a single-pass transmembrane protein expressed on the cell surface (Hortsch et al., 1990). Dimers of the secreted protein Amalgam (Ama), which are expressed in the direct AMPs, are able to bind to Nrt on two different cells and promote their adhesion (Frémion et al., 2000; Zeev-Ben-Mordehai et al., 2009).

Figure 6 with 1 supplement see all

Download asset Open asset

Nrt and mid expression increased dramatically in direct AMPs from 96 to 120 hr based on our single-cell data (Figure 6A, B) and antibody staining (Figure 6D–G), whereas ptc is expressed at similar levels at both time points (Figure 6C, H, I). Surprisingly, while Nrt and Mid expression patterns included all posterior-localized direct AMPs, the expression of both genes extended into anterior-localized direct AMPs, which do not express Ptc. Due to their high expression levels in the posterior-localized AMPs, we hypothesized that expression of Nrt and Mid, at least in the posterior AMPs, was influenced by Hh signaling.

To determine if Nrt and mid are downstream Hh-signaling targets, we examined if perturbing the Hh pathway within the AMPs would alter their expression. The knockdown of smo and consequent reduction in Hh signaling resulted in a dramatic decrease in both Nrt and Mid expression in the direct AMPs at 120 hr AEL (Figure 6J, K, N). Remarkably, this was observed in both the posterior- and anterior-localized AMPs alike. We tested if increased Hh signaling would be sufficient to induce expression of Nrt and Mid by driving ci^3m in the AMPs. This resulted in the ectopic expression of both Nrt and Mid in all of the AMPs, although expression was higher within direct AMPs (Figure 6L–N). Thus, all AMPs are capable of inducing Nrt and Mid expression in response to Hh signaling, but only the AMPs closest to Hh-secreting epithelial cells do so under physiological conditions. This relationship between ptc, mid, and Nrt is specific to the AMPs as we did not observe correlation between ptc and either Nrt or mid within the epithelium. These experiments indicate that Hh signaling is important for proper patterning of the AMPs. However, Hh signaling is not required to regulate all aspects of the direct AMP cell fate as high Ct expression was unaffected following the manipulation of Hh signaling by smo knockdown (Figure 6—figure supplement 1A, B). This suggests that Hh signaling is acting independent from ct to specify the fate of a subpopulation of the AMPs.

To evaluate the functional consequences of reducing Nrt and mid expression, we used RNAi to reduce their expression. The RNAi line directed against mid, however, failed to significantly reduce Mid protein levels. In contrast, Nrt knockdown reduced Nrt levels (Figure 6O) and resulted in adults with defects in the posterior DFMs (Figure 6P; compare to Figure 5R), albeit not as severely as knockdown of smo. This result shows that the Hh-signaling downstream target Nrt is critical for proper DFM development.

The AMPs with high Ptc expression are in the vicinity of Hh-expressing epithelial cells, which is consistent with the notion that Hh is a short-range morphogen capable of acting on cells within 40 µm from the source. However, we observed expression of the two target genes, Nrt and Mid, in anterior-localized AMPs far beyond 40 µm from the Hh source. Knockdown of Hh signaling via smo^RNAi expression abolished the expression of both genes in anterior- and posterior-localized AMPs alike (Figure 6J, K), ruling out the possibility of their activation in a Hh-independent manner. This raised the possibility that these cells were initially closer to the Hh source earlier in development and had subsequently moved anteriorly either as a result of passive displacement or active migration. In this scenario, the expression of Nrt and Mid would be expected to perdure for longer than the expression of Ptc. To test whether anterior-localized AMPs had past activation of the Hh-signaling pathway, we used a lineage-tracing method to identify cells that had previously expressed ptc (ptc^ts>FLP, GFP^FLP-Out; see Materials and methods). We labeled all cells that displayed activity of ptc-Gal4 during late second and early third instar, and observed lineage labeling (indicated by GFP expression) of all direct AMPs underlying the posterior compartment of the epithelium as well as a trail of anterior AMPs that recapitulate the domain of anti-Nrt staining (Figure 6Q). These results indicate that these anterior GFP-positive AMPs likely descended from cells that previously expressed high levels of ptc, and that past activation of the Hh pathway is likely to be responsible for the expression of both anterior- and posterior-localized Nrt and Mid (Figure 6R).

Data obtained from scRNAseq experiments can provide both spatial and temporal information that give us a better understanding of how cells diversify and then stabilize their transcriptomes during development, and point to ways in which they interact with each other. Moreover, we have been able to visualize the expression of genes from our dataset in a three-layered virtual wing disc – a format that will be useful for developmental biologists.

One interesting observation from our data is that spatial positioning within the wing disc is highly informative of the transcriptional state of cells. In particular, the proximodistal axis of the disc epithelium is one of the primary stratifying features within our single-cell data. Epithelial cell clusters were easily mapped back to sub-regions within the notum, hinge, pouch, and PE. In contrast, although the cells of the anterior and posterior compartments have been separated by lineage since early in embryogenesis, we observe less differential expression between the two compartments. Thus, position along the proximodistal axis has a far greater influence on the transcriptome of a cell than its ancestry.

In the epithelium, we observe that most of the major cell types observed at late L3 (120 hr AEL) are already present at mid L3 (96 hr AEL). However, the transcriptomes of the two major populations of AMPs, those that give rise to the DFMs and IFMs, diverge significantly during this time interval. At 96 hr, AMPs appear to be in a relatively naive state; canonical markers for the direct and indirect flight precursors, ct and vg, both show relatively uniform expression at the mid L3 stage. At the late L3 stage, we observe more distinguishable differences between the transcriptomes of direct and indirect cell types. Both ct and vg have greater differential expression in the AMPs at this time point. The earlier stabilization of epithelial cell fates then provides a stable platform for the generation of spatially localized signals that regulate myoblast numbers and provide instructive signals for fate specification.

Ths and Pyr have previously shown to regulate the spreading of mesodermal cells during embryogenesis, but a role for these ligands in regulating myoblast numbers was not previously appreciated. We show these ligands are necessary for AMP survival and that increased levels of these ligands can promote AMP proliferation even at sites distant from the notum. Indeed, increased expression can induce a dramatic overproliferation of the AMPs, even as late as third instar. Ectopic FGF signaling is sufficient to support AMP viability and cell proliferation at other locations in the wing disc, including underneath the pouch. Thus, the localized expression and level of Ths and Pyr secreted by epithelial cells in the notum could provide sufficient trophic support to generate the appropriate number of AMPs during normal development. While this work was in preparation, another group independently showed that the ths-Gal4 line is expressed in the notum epithelium and that reducing ths function reduces AMP numbers (Vishal et al., 2020).

We have also demonstrated that ectopic and elevated levels of expression of Ths or Pyr can draw AMPs out of the notum region, all the way to the ventral hinge and around the ventral edge of the disc proper onto the peripodial epithelium. With constitutive FGF ligand expression driven by dpp-Gal4, we did not observe AMP migration to the lateral regions of the disc, which are far from the ectopic FGF source. However, when FGF expression is initiated in L3 by controlled Gal80^ts repression, we observed a large number of laterally located AMPs. We attribute the difference in these two scenarios to ectopic AMPs along the dpp stripe possibly serving as an FGF sink. When FGF is turned on later in development, there are no ectopic AMPs. In their absence, FGF could reach more lateral portions of the disc, thus allowing AMP emigration to those regions. Furthermore, the expression of Pyr in the wing pouch, which is separated from the notum by the dorsal hinge, was sufficient to promote colonization of the pouch region by AMPs. The AMPs can be induced to colonize new domains late into larval development by the expression of either FGF ligand. These results suggest that under physiological conditions AMPs remain beneath the notum epithelium because there is insufficient FGF outside of this region. Altogether, our work illustrates that a source of epithelial FGF is critical for forming the AMP niche and that levels of FGF regulate both the location and number of the AMP cells.

We have shown that the anteroposterior identities of the disc epithelium are important for proper specification of gene expression within the underlying AMPs. One powerful advantage of scRNAseq as opposed to its bulk sample predecessor is the ability to measure gene co-expression within subpopulations of a tissue. We leveraged this advantage to identify two novel Hh-signaling targets, Nrt and mid, within AMPs. While we currently do not know whether Nrt and mid are direct targets of Ci, both genes do have consensus Ci-binding sites within potential regulatory regions. Nrt is a single-pass transmembrane protein. Its extracellular ligand, Amalgam, has more widespread expression in the direct myoblasts and is expressed at comparable levels at both 96 and 120 hr AEL. Two molecules of Amalgam can form homodimers, and each is capable of binding to Nrt on different cells (Frémion et al., 2000; Zeev-Ben-Mordehai et al., 2009). Thus, an effect of Hh-induced expression of Nrt in a subset of the direct adult flight muscles might be to promote aggregation of Nrt-expressing cells at a later stage of development.

An unexpected observation was that AMPs beneath the anterior compartment, distant from the epithelial source of Hh, express both identified Hh targets Nrt and Mid. However, this expression is dependent upon Hh signaling since knockdown of smo blocks gene expression. Although we cannot completely exclude the possibility that a second signal from posterior AMPs activates Nrt and mid expression in these cells, our lineage-tracing experiments favor a model where a subset of the direct AMPs are generated posteriorly and move anteriorly during the course of development. Such movement could be due to a process of active migration in response to hitherto unknown external cues or to displacement as a result of oriented cell division. Understanding the mechanistic basis of AMP migration would represent an exciting avenue of future research.

Our work has provided a base for the study of heterotypic interactions in the developing wing disc during conditions of normal growth and demonstrate that such interactions can have a major effect on cell number, cell migration, and cell fate in the wing disc. By examining receptor–ligand expression patterns in conjunction with spatial mapping of our single-cell data, our analysis provides many hints of signaling pathways that may function between the disc epithelium and the AMPs and also within subsets of cells with each of these populations that provide multiple avenues for future investigations.

For each sample, approximately 250 staged Drosophila wing-imaginal discs were dissected within 1 hr. The collected tissue was then transferred to a microcentrifuge tube and incubated within a dissociation cocktail consisting of 2.5 mg/mL collagenase (Sigma #C9891) and 1× TrypLE (Thermo Fisher #A1217701) in Rinaldini solution (modified from Ariss et al., 2018). The sample tube was placed horizontally on a shaker machine operating at 225 rpm for 25 min at room temperature (method modified from Ariss et al., 2018). At the 10, 20, and 25 min marks, the tube was flicked 20 times for additional mechanical dissociation. Dissociation was halted by centrifuging the sample at 5000 rpm for 3 min, aspirating the dissociation cocktail, and then adding in 1 mL of cold PBS-10% FBS. The cell pellet was mixed by pipetting up-and-down approximately 25 times with a 1 mL pipette for additional mechanical dissociation, and then centrifuged again at 5000 rpm for 3 min. The media was replaced with cold PBS-1% FBS, and the cell pellet was resuspended in preparation for FACS.

FACS of the sample was performed on a BD FACSAria Fusion flow cytometer. Dead cells were identified and removed via the addition of propidium iodide to the sample, and high-quality single cells were sorted into cold PBS-10% FBS. Cell concentration of the post-FACS sample was assessed by a hemocytometer and adjusted 1000 cells per µL.

Single-cell suspensions were barcoded for single-cell RNA sequencing with the 10X Chromium Single Cell platform (v2 chemistry). Barcoded samples were sequenced on an Illumina NovaSeq (S2 flow cell) to over 60% saturation.

The 10X Genomics Cell Ranger pipeline (v2.2.0) was used to align the sequencing reads to the Drosophila melanogaster transcriptome (version 6.24). The data was analyzed using the R and Python programming languages, primarily utilizing the packages Seurat v3 (Stuart et al., 2019) and scVI v0.4.1 (Lopez et al., 2018).

Our standard analysis pipeline is as follows: first, each dataset was analyzed separated using the standard Seurat pipeline, with no cells filtered, 30 principal components calculated, and clustering resolution set to 2.0 (all other parameters remained default). We then removed cell clusters with an abundance of low-quality cells (defined as clusters with mean number of genes detected per cell [nGene] was less than one standard deviation below the mean nGene of all cells in the dataset). Additionally, we found that each dataset had a cluster with markers for both AMP and epithelial cell types (e.g., SPARC and Fas3) and unusually high mean nGene; this cluster was suspected to be AMP-epithelial doublets and was also removed. Clusters were then split into AMP and epithelial cell subsets based on the expression of known marker genes. Cells within each subset were subsequently filtered if either (1) their nGene count that was outside the mean nGene of the subset ±1.5 standard deviations or (2) their percentage of reads for mitochondrial genes that was greater than 1.5 standard deviations above the mean mitochondrial read percentage of the subset.

Data subsets were harmonized into collective AMP or epithelium datasets using scVI. The scVI VAE model consisted of 2 layers (n_layers = 2) and 20 latent dimensions (n_latent = 20), with a negative-binomial reconstruction loss (reconstruction_loss=‘nb’). The model was trained on variable genes selected by Seurat’s variance-stabilizing transformation method; 1000 (for epithelial subsets) or 2000 (for AMP subsets) variable genes were calculated for each inputted batch, and then the union of these genes was supplied to scVI. The following parameters were used for model training: train_size = 0.75, n_epochs = 400, and lr = 1e-3 (other parameters were left as default). Cell clustering and UMAP was performed using Seurat on the latent space derived from the scVI model. After harmonization, clusters were re-examined for doublet characteristics; clusters with a mean nGene count greater than one standard deviation above the mean nGene count of all cells were removed, as were clusters that displayed markers for both AMP and epithelial cell types. Identified hemocyte and tracheal cells were also separated out. scVI and Seurat were both re-run on the datasets to generate our final AMP and epithelial cell atlases (Figure 1H, N).

To generate our full cell atlas consisting of all cell types (Figure 1B, C), we merged and harmonized the cells in the AMP and epithelium cells atlases along with the separated hemocyte and tracheal cells. scVI and Seurat were run as previously described, with the scVI model trained on the union of the top 2000 variable genes for each batch as calculated by Seurat. No additional cell filtering was performed after harmonization.

For visualizing data on UMAPs and dot plots, we calculated normalized and scaled expression counts using Seurat’s NormalizeData and ScaleData functions, respectively, with default parameters. For the normalized data, raw counts were normalized by total unique molecular identifiers (UMIs) per cell, multiplied by 10,000. Natural-log normalized data is used for expression levels visualized with UMAP. For the scaled data, the natural-log normalized data is scaled for each gene, such that the mean expression is 0 with a standard deviation of 1. Scaled data is used for expression visualization on the dot plots.

Cells were classified as male or female by their expression levels of the dosage compensation complex genes lncRNA:roX1 and lncRNA:roX2 (Franke and Baker, 1999; Meller and Rattner, 2002), which are both expressed almost exclusively in male cells. For both genes, we examined the natural-log normalized expression counts (calculated by Seurat’s NormalizeData function), computed the density over the data, and identified the first local minima as a threshold (see Figure 1—figure supplement 4C, D). Cells that were above the threshold for either lncRNA:roX1 or lncRNA:roX2 were classified as male; otherwise, they were classified as female. From this, we assigned 8097 cells as male and 11,788 cells as female, which roughly matches the size ratio between male and female wing discs given that female discs are larger. We removed cell sex stratification by processing male and female AMPs as separate batches (for each actual batch) within scVI (see Figure 1—figure supplement 4E, F for comparison of data before and after cell sex correction).

We observed significant data stratification that correlated with a number of cell cycle-related genes, such as proliferating cell nuclear antigen (PCNA) and cyclin B (CycB) (Yamaguchi et al., 1990; Lehner and O'Farrell, 1990), indicating that our data was split between S phase and non-S phase (Figure 1—figure supplement 5A–C). Definitive classification of cells into cell cycle stages is difficult because expression of these genes is not typically demarcated sharply into specific cell cycle stages. To remove cell cycle stratification from our data, we examined the correlation of each scVI latent dimension with the expression levels of highly variable cell cycle genes and found that one latent dimension was strongly related (Figure 1—figure supplement 5D). By masking this latent dimension from our downstream analysis (e.g., clustering and UMAP), we effectively diminished cell cycle stratification. Figure 1N shows the UMAP of our AMP data after subtraction of cell sex and cell cycle stratification, which allowed us to focus our analysis on different cell types within the AMPs (to see how each correction affected the AMP data, see Figure 1—figure supplement 5A–C).

When examining clusters, genes were considered to have significant differential expression if they had (1) a false discovery rate (FDR) <0.05 (as calculated via Wilcoxon test), (2) a natural-log fold-change of 0.15 or more, and (3) a percent expression of at least 15% in one of the two populations in the comparison. This test was performed with Seurat’s FindMarkers function, using a Wilcoxon test (test.use = ‘wilcox’).

When evaluating differential expression between clusters of epithelial cells (e.g., Figure 1—figure supplement 2A), we used a one cluster vs. all analysis. When evaluating differential expression between direct and indirect AMPs (e.g., Figure 1—figure supplement 5E), we compared cells of the two groups as classified in Figure 1—figure supplement 5F. In these cases, differential expression statistics (i.e., FDR and fold-change) are obtained by combining the cells (across batches) in each group and conducting a single comparison. When evaluating differential expression between time points (e.g., Figure 1—figure supplement 1H, Figure 1—figure supplement 2C, Figure 1—figure supplement 6A) (which would be inherently confounded with batch effects since time points were collected across separate sequencing experiments), we took a conservative approach and only considered genes that were consistently significant (by the criteria defined above) in each temporal pairwise comparison (i.e., DE analysis was conducted between all temporal pairs: 96h1 vs. 120h1, 96h2 vs. 120h1, 96h2 vs.120h1, 96h2 vs.120h2). We report the natural-log of the average value for these pairwise comparisons, and the maximum FDR calculated (see Supplementary files 1–3).

We assembled reference gene expression patterns from a number of sources (Held, 2002; Butler et al., 2003) and based our starting geometry on the disc proper from images in Bageritz et al., 2019. The images were processed in Adobe Photoshop and assembled in R with EBimage (Pau et al., 2010) to generate binarized gene expression reference for the AMPs, disc proper, and peripodial epithelium. The geometry of the three-layered model is provided in Supplementary file 4, and the binarized reference gene expression patterns are provided in Supplementary file 5. We used DistMap (Karaiskos et al., 2017) to statistically map single cells back to the reference. With this virtual wing disc model, we used DistMap to calculate a ‘virtual in situ’ or a prediction of gene expression patterns. This is based on the detected gene expression with the single-cell data and the mapping location to calculate relative expression values for our model. We mapped the AMP and epithelial cells separately as this improved how well the model predicted genes with known expression patterns. In addition, we used the scVI imputed gene expression values when mapping the cells to the reference.

From FlyBase, we assembled a list of genes encoding receptors and ligands from the following 19 pathways of interest: Wnt/Wingless, FGF, Hedgehog, PDGF/VEGF, JAK-STAT, activin, BMP, Fat-Ds, Slit-Robo, ephrin, toll/toll-like, semaphorin, Notch, insulin-like, fog, torso, miple, EGFR, and TNF. For our analysis, we only examined pathways in which at least one receptor or ligand was either (1) differentially expressed within one of the major domains of the epithelium (notum, hinge, pouch, or PE) when compared to all other epithelial cells, (2) differentially expressed within one of the major domains of the AMPs (direct or indirect cells) when compared to each other, or (3) differentially expressed between all epithelial cells vs. all AMP cells. These pathways (and their receptors and ligands) are shown in Figure 4H.

The stocks used in this study include the following lines from the Bloomington Stock Center: R15B03-GAL4 (BL49261); G-TRACE (BL28280, 28281) (Evans et al., 2009); UAS-FLP, Ubi-FRT-stop-FRT-GFP^nls (BL28282); pdm3-GFP (BL60560); grn-GFP (BL58483); ptc-GAL4 (BL2017); tub-GAL80^ts (BL7108); dpp-GAL4 (BL1553); dpp-GAL4, tub-GAL80^ts, UAS-dCas9.VPR (BL67066); nub-GAL4 (BL25754); hh-GAL4 (Tanimoto et al., 2000); ap-GAL4 (BL3041); fng-Gal4 (BL9891); htl-GAL4 (GMR93H07-GAL4, BL40669) is an enhancer within the first intron of the htl gene; UAS-smo^RNAi (primarily BL43134, but also BL27037, 62987 in Figure 5—figure supplement 2); UAS-Nrt^RNAi (BL28742). Drosophila stocks from other labs: UAS-ths and UAS-pyr (A Stathopoulos); UAS-ci^3m (D Kalderon). TRiP-CRISPR driven overexpression of pyr was conducted with a guide RNA that targets the upstream transcriptional start site, P{TOE.GS00085}attP40 (BL67537), and works together with a nuclease-dead Cas9 fused with a transcriptional activator domain, UAS-dCas9.VPR to cause gene activation (BL67055) (Lin et al., 2015).

Imaginal discs, unless otherwise noted, were fixed in 4% paraformaldehyde (PFA) for 15 min, permeabilized in PBS plus 0.1% Triton X-100, and blocked in 10% normal goat serum. For anti-Nrt antibody staining, we substituted the Triton X-100 for 0.05% saponin. The following antibodies were used from the Developmental Studies Hybridoma Bank (DSHB): mouse anti-Cut (1:200, 2B10); mouse anti-Ptc (1:50, Apa-1); mouse anti-Nrt (BP 106 anti-Neurotactin); and mouse anti-Wg (1:100, 4D4). The following antibodies were gifted: rat anti-Twist (1:1000, Eric Wieschaus), rabbit anti-Midline (1:500, James Skeath), and rat anti-Zfh2 (1:100, Chris Doe; Tran et al., 2010). The following antibodies are from commercial sources: rabbit anti-Dcp1 (1:250, Cell Signaling); rabbit anti-GFP (1:500, Torrey Pines Laboratories, Secaucus, NJ); chicken anti-GFP (1:500, ab13970 Abcam, Cambridge, UK); rabbit anti-beta-galactosidase (1:1000, #559762; MP Biomedicals, Santa Ana, CA); and rabbit anti-PHH3 (1:500, Millipore-Sigma). Secondary antibodies were from Cell Signaling. Nuclear staining with DAPI (1:1000).

For EdU staining, we incubated live discs in fluorescent EdU incorporation solution for 1 hr, following the protocol for the Click-iT EdU Cell Proliferation Kit, Alexa Fluor 555 (ThermoFisher C10338). After the incubation, discs were fixed in 4% PFA for 15 min, before proceeding with standard antibody stainings as detailed above.

To ectopically express FGF ligands starting at mid or late third instar development (e.g., Figure 4G, H), we used a temperature-sensitive dpp-Gal4 stock (dpp-GAL4, tub-GAL80^ts, UAS-dCas9.VPR) (BL67066) with UAS-pyr and UAS-ths lines. In-vial egg lays were collected over 8 hr, and larvae were initially raised at 18°C. Larvae were shifted to 30°C (to relieve Gal80^ts repression) at either 48 or 24 hr prior to dissection, corresponding to mid and late third instar FGF activation, respectively.

To lineage-trace cells that expressed ptc-Gal4 during the second and early third instars, we used a temperature-sensitive, Ptc-dependent FLP-Out system (ptc-Gal4, tub-Gal80^ts, UAS-FLP, Ubi-FRT-stop-FRT-GFP^nls). Under this system, inactivation of Gal80^ts repression permits cells that express ptc-Gal4 (and their descendant cells) to become permanently labeled with GFP by an FLP–FRT recombination event. In the experiment shown in Figure 6Q, first instar larvae were collected and reared at 18°C. At 5 days post-egg lay, larvae were shifted to 30°C for 24 hr, then shifted back to 18°C for 48 hr prior to dissection.

Wing discs were imaged on a Zeiss Axioplan microscope with Apotome attachment, using ×10 and ×20 objectives. Image files were processed with ImageJ software. For each of the genotypes examined, we examined at least eight discs and have reported representative results in this paper.

To image adult flight muscles, male flies aged a minimum of 2 days after eclosion were anesthetized and submerged in 70% ethanol with dry ice. The thorax was isolated by removing the head, wings, legs, and abdomen. Thoraces were bisected sagittally with a 11-blade scalpel blade. For DVMs, the DLMs, leg muscles, and excess cuticle were removed from hemithoraces. For DFMs, the DLMs, DVMs, leg muscles, and excess cuticle were removed from hemithoraces. The DLMs, DVMs, and DFMs were fixed in 4% PFA for 2 hr. Muscles were then rinsed three times and permeabilized in 0.3% PBST for three cycles, 15 min each on a nutator. Hemithoraces were incubated in rhodamine phalloidin (1:200) and DAPI (1:500) in 0.3% PBST, then rinsed three times and washed in 0.3% PBST for three cycles, 15 min each on a nutator. Hemithoraces were mounted in a depression slide using antifade mountant. DLMs and DVMs were imaged with a ×10 objective using a Zeiss Axioplan microscope. DFMs were imaged with the ×20 and ×63 objectives using confocal microscopy.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank A Stathopoulos, J Skeath, C Doe, and E Wieschaus for stocks and reagents and M Frolov for discussions. We thank current and former members of the Yosef lab for advice and guidance. We thank T Ashuach and D DeTomaso for discussions regarding computational analyses. We thank current and former members of the Hariharan lab for advice, guidance, and help with collection of samples. We thank C Wingert for technical assistance. We thank the Bloomington Stock Center, DRSC/TRiP Functional Genomics Resources, and Developmental Studies Hybridoma Bank for stocks and reagents. This work was funded by NIH grant R35 GM122490.

1. Received: July 21, 2020
2. Accepted: March 21, 2021
3. Accepted Manuscript published: March 22, 2021 (version 1)
4. Version of Record published: April 4, 2021 (version 2)

© 2021, Everetts et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.