Single cell RNA sequencing of the developing Drosophila larval eye discs
We performed single cell RNA sequencing (scRNA-seq) from six time points during early larval Drosophila eye disc development using the 10x Genomics Chromium platform. To capture the gene expression dynamics from the initiation of the morphogenetic furrow (MF) to the differentiation of most major cell types, we first profiled eye disc cells from four time points spanning late second instar (69 hr after egg laying (AEL)) to third instar larval stages (72 hr, 75 hr and 78 hr AEL) (Fig. 1A). The majority of eye disc cells at these time points are undifferentiated and actively dividing and only a minority of differentiated cell types are present. We therefore combined the data from these time points (termed "early larval") for further analyses. Removal of dead or dying cells and multiplets, followed by dimension reduction using UMAP resulted in 3,624 eye cells. The UMAP plot shows clusters that correspond to all expected cell identities in the early larval eye disc (Fig. 1B). As expected, the anterior undifferentiated (AUnd) cell cluster comprises the most abundant cell type (1723 cells) in the combined early larval data set. The AUnd cluster was identified using Optix [18], eyeless (ey) [19] and teashirt (tsh) [20], which are predominantly expressed in AUnd cells in developing larval eye discs (Fig. 1C and 2B). The morphogenetic furrow (MF) cell cluster was identified using decapentaplegic (dpp) [21] and rotund (rn) [22]. These genes are expressed in the MF and dpp is required for the initiation and propagation of the MF across the eye disc [21] (Fig. 1C and 2C,D). We also identified two cell populations that express known photoreceptor markers. One of these cell clusters is identified as R8, as it expresses atonal (ato) [23, 24] and senseless (sens) [25]. ato and sens are specifically expressed in R8 and are required for the differentiation of R8. ato is also expressed in the MF and in a stripe immediately anterior to the MF (Fig. 1C and 2E,F). The second cell cluster expresses rough (ro) and seven up (svp) (Fig. 1C and 2G,H). While ro is expressed in R2/5 and R3/4 [26], svp expression is confined to R3/4 and R1/6 [27]. Since this second cluster shows expression of both genes without segregation into individual photoreceptor subtype clusters, we reasoned that this cluster comprises the R2/5, R3/4 and R1/6 subtypes and therefore named this the R cell cluster. Interestingly, the AUnd, MF and arrangement of photoreceptor clusters in the cluster plot appears as a progression that is reminiscent of the temporal nature of the larval eye disc in vivo. The posterior undifferentiated cell cluster (PUnd) was identified using Bar-H1 (B-H1) and lozenge (lz) which are known to be expressed in PUnd cells [28, 29] (Fig. 1C and 2L). Our data also shows cell clusters corresponding to the peripodium (PPD) and the posterior cuboidal epithelium (PC). Ance and ocelliless (oc) mark the PPD [14] (Fig. 1C), while the pair rule genes odd skipped (odd), drumstick (drm) and sister of odd and bowl (sob) are confined to the PC (Fig. 1C) and are required to initiate retinogenesis [30].
To complement our early larval scRNA-seq data, we performed scRNA-seq on two additional larval time points: 88 hr and 96 hr AEL. We combined these two data sets (termed "mid-larval"), filtered low quality cells and generated a UMAP cluster plot with 9,600 cells that shows all expected cell identities in the eye disc (Fig. 1D). Similar to the early larval data, undifferentiated cell clusters were identified using Optix (AUnd), hairy (h) (Preproneural (PPN)), dpp (MF), PCNA (second mitotic wave (SMW)), and lz (PUnd) (Fig. 1E). Photoreceptor cell clusters were identified using sens (R8), ro (R2/5 and R3/4), svp (R3/4 and R1/6), prospero (pros) [31] (R7) and spalt major (salm) (R3/4 and R7) [14] (Fig. 1E and 2I,J,K). Cone cells were identified using cut (ct) [32] and pros expression. pros is expressed in both R7 and cone cells (Fig. 1E and 2I). Importantly, the mid-larval eye scRNA-seq cluster plot closely resembles our previously published late larval cluster plot [14] with each cell cluster appearing with progressive patterns of gene expression suggesting a temporal component, reminiscent of the temporal nature of the larval eye disc. The photoreceptor clusters appear as strands of cells that are either connected to the MF (R8, R2/5 and R3/4) or the PUnd (R1/6 and R7) cell clusters. Though our mid-larval cluster plot shows R1/6, R7 and cone cells, they do not segregate distinctly and appear together as one cluster. Taken together, our early and mid-larval data show that we captured all major types from eye discs.
Gene expression profiles of anterior undifferentiated cells
Since our early larval data profiles only ~ 3,600 cells, we combined early larval and mid-larval data for downstream analyses. This merging of data enables good representation of each cell identity as well as profiles cell types at high resolution. This is possible because each larval eye disc is a spatiotemporal continuum with cells of different developmental ages in the same tissue. The combined data showed all expected cell identities that were observed in both early and mid-larval time points (Fig. 2A). We then performed differential gene expression analyses between cell clusters to understand their gene expression profiles. As expected, Optix, tsh, twin of eyeless (toy) and ey are among the top 15 differentially expressed markers in the AUnd cell cluster (Supplemental Data 1). It is well known that these genes are expressed in undifferentiated cells prior to the initiation of the MF as well as anterior to the MF after initiation [18–20]. We also identified CG17211 and Organic anion transporting polypeptide 74D (Oatp74D) (Fig. 3B,C) as two examples of novel AUnd markers. These are expressed predominantly in the AUnd cell cluster in which Optix is expressed and are likely AUnd markers. The complete list of AUnd markers is shown in Supplemental Data 1. We further performed Gene Ontology (GO) term enrichment using these marker genes and, as expected, terms related to imaginal disc development and proliferation showed enrichment (Supplemental Data 2). These include eye-antennal disc development (GO:0035214) and positive regulation of cell population proliferation (GO:0008284).
We also performed principal component analysis (PCA) on the AUnd cell cluster to identify genes that change as a function of time and may be involved in dynamic cellular processes such as cell division and differentiation. As expected, PCA on the AUnd cluster identified Optix with the second highest loading weight. The long non-coding RNA lncRNA:CR33938 is the top PC1 gene for the PUnd cluster. Odorant-binding protein 56a (Obp56a), Oatp74D, fatty acid binding protein (fabp) and narrow (nw) are the other top 6 PC1 genes. The top 1000 PC1 genes of the AUnd cluster with loading weights are shown in Supplemental Data 3. These top PC1 genes may be involved in dynamic processes of the AUnd cell cluster such as cell division, growth, and differentiation.
Gene expression profiles of the morphogenetic furrow
Differential gene expression analyses of the MF cluster show five Enhancer of split (E(spl)) genes in the top 10 marker gene list (Supplemental Data 1). E(spl) genes are the downstream effectors of Notch signaling pathway [33]. The role of Notch signaling in the initiation and propagation of MF is well known [33] and it is expected that the E(spl) genes appear as differentially expressed genes in the MF. In addition to other known MF genes such as dpp and rn, Dipeptidase B (Dip-b) and Phosphodiesterase 8 (Pde8) (Fig. 3D,E) are predominantly expressed in MF and may play a role in MF initiation and/or progression.
To identify putative maturation genes, we subjected the MF cluster to PCA. As expected, the E(spl) genes show the highest loading weights in PC1 along with other known MF genes such as rn. Similarly, scabrous (sca) and ato are in the top 30 PC1 genes. These genes are involved in the specification of R8 from the pool of undifferentiated progenitor cells [23, 34]. In addition, CG15282, gliolectin (glec), Protein kinase cAMP dependent regulatory subunit type 2 (Pka-R2), Bearded (Brd) and Dip-B are other top 30 PC1 genes whose function in the progression and maturation of cells in MF has not been previously reported. Taken together, these data suggest that we have identified several potential candidate genes that may drive specification of cell types in the MF.
Cell cycle and DNA replication genes are enriched in the second mitotic wave cluster
Differential gene expression analyses of the SMW cluster reveals known markers such as dac and the cell cycle-related genes Claspin and PCNA [14]. In addition, Odorant-binding protein 44a (Obp44a) and Tetraspanin 66E (Tsp66E) (Fig. 3F,G) and other markers were also identified and are shown in Supplemental Data 1. GO term analyses shows enrichment for terms involved in the cell cycle and DNA replication. Some of the GO terms include pyrimidine deoxyribonucleoside monophosphate biosynthetic process (GO:0009177), deoxyribonucleotide biosynthetic process (GO:0009263), leading strand elongation (GO:0006272) and cell cycle DNA replication initiation (GO:1902292). In particular, the DNA replication licensing factor and Minichromosome maintenance (MCM) complex genes [35, 36], which are replicative helicases for DNA replication, initiation and elongation, are associated with these GO terms. This is expected as the cells in the SMW are actively undergoing cell division and preparing to differentiate into R1/6, R7 or cone cells. We then subjected the SMW cluster to PCA to identify putative genes that may prime these cells to undergo division and differentiation. anterior open (aop/yan) shows the highest loading weight in PC1 of the SMW. The role of aop in regulating cell fate transitions in the eye is well documented [37]. target of wit (twit), Obp44a, CG42342 and Tsp66E are the other top 5 PC1 genes in the SMW, whose function in the SMW is currently unknown.
Gene expression profiles of posterior undifferentiated cells
Our cluster plot shows that the PUnd cells cluster separately from AUnd cells, reflecting transcriptomic differences between undifferentiated cells anterior and posterior to the MF. PUnd cells are developmentally more mature progenitors compared to AUnd cells and it is expected that the two undifferentiated cell populations will segregate on the UMAP plot [14]. In addition to the known PUnd markers lz and B-H1, differential gene expression analyses identify several PUnd markers such as Gasp and CG4374 (Fig. 3H,I). These genes are specifically expressed in the PUnd cell cluster and their function is currently not known. PCA on PUnd cell cluster identified CG13071, CG9691, Gasp, kekkon 1 (kek1) and scarface (scaf) as the top 5 genes with the most variation in the PUnd cluster. GO analyses of the PUnd cell cluster showed enrichment for negative regulation of compound eye photoreceptor cell differentiation (GO:0110118), positive regulation of compound eye retinal cell programmed cell death (GO:0046672) and negative regulation of cell fate specification (GO:0009996). These terms are expected because photoreceptors have undergone differentiation and some of the undifferentiated cells are eliminated through apoptosis.
Photoreceptor subtypes and cones cells appear as distinct clusters
Our data show five cell populations that appear as distinct strands, which we identified as photoreceptor subtypes (Fig. 2A). Three of these strands are connected to the MF and correspond to R8, R2/5 and R3/4, whereas R1/6 and R7 are connected to the PUnd cell cluster. This arrangement of clusters resembles the sequence of differentiation of cell types in the physical eye disc. R8, R2/5 and R3/4 differentiate first from cells in the MF, while R1/6, R7 and cones differentiate after the SMW and therefore are derived from posterior undifferentiated cells. Similar to our late larval data set [14], photoreceptor cells cluster as temporal strands with less mature cells located near the MF or PUnd clusters, while more mature cells are located at opposite ends of each strand. The expression of sens, ato and bride of sevenless (boss) reveals the temporal nature of the R8 strand (Fig. 4A-C). ato is expressed in the MF and in less mature R8 cells, sens is expressed in all R8 cells, and boss is expressed only in mature R8 cells. Our data show ato expression in the MF cluster and proximal tip of the R8 strand that is close to the MF but no ato expression is observed in the rest of the strand (Fig. 4A). sens mRNA is detected along the entire stream (Fig. 4B), while boss expression is observed only in mature R8 cells that are at the distal tip of the R8 strand (Fig. 4C). These data collectively indicate that the R8 photoreceptor strand is a continuum of cells arranged progressively according to their developmental age. Likewise, other photoreceptor subtype streams also exhibit such temporal dynamics. In addition, similar to our late larval eye scRNA-seq data, the R8, R2/5 and R3/4 strands appear to connect to a common cluster at their more mature ends and is named the 'Late R' cell cluster (Late R, Fig. 3A) [14]. Genes related to axon projection, guidance and synapse formation dominate the differentially expressed gene list of the Late R cell cluster. Cone cells appear as a separate cluster that is connected to the PUnd cluster. Our cluster plot shows that R1/6, R7 and cones are distinct but are part of one major cluster. Differential gene expression analyses identified many photoreceptor and cone subtype-specific markers and these are shown in Supplemental Data 1.
As an in vivo validation step we used Trojan-Gal4 (T2A-Gal4) lines [38] to test the expression of marker genes in the larval eye disc. T2A-Gal4 insertions express Gal4 under the control of endogenous promoters and recapitulate the endogenous pattern of expression of the gene in which the transgene is inserted. Our scRNA-seq data show that CG13532 is predominantly expressed in R2/5 and some expression is also observed in late R8 (Fig. 4D). We used CG13532-T2A-Gal4 to drive expression of a nuclear localized mCherry reporter (UAS-mCherry-nls) and costained larval eye discs with mCherry, Rough (Ro, an R2/5 marker) and Embryonic lethal abnormal vision (Elav, a pan-neuronal marker) antibodies. We observed that mCherry, Ro and Elav colocalize in two cells per ommatidium posterior to MF (Fig. 4E-E''') suggesting that CG13532 is expressed in R2/5 cells and is an R2/5 marker. Furthermore, these results suggest that our scRNA-seq data can predict the expression and distribution of unknown genes in the larval eye disc.
Gene expression profiles of peripodial cells are distinct from other cell types
The peripodial membrane is a thin squamous epithelium that covers the eye disc proper and is important for its patterning and growth. In addition, short narrow cells comprising the 'posterior cuboidal margin' (PC) are also present at the posterior border between the eye disc and peripodial layers. Our data show two clusters that correspond to the peripodial cells (PPD) and the PC (Fig. 3A). The PPD was identified using Ance and oc, which are expressed in the PPD. Differential gene expression analyses of the PPD cluster identified several PPD markers including Odorant-binding protein 99a (Obp99a) and CG14984, which are expressed in the PPD (Supplemental Fig. 1A,B). CG14984 is also expression in the AUnd. Several other larval eye PPD markers are shown in Supplemental Data 1. As expected, GO term enrichment analyses of the PPD shows terms related to tissue morphogenesis and appendage development as some of the head structures are derived from cells of the PPD. Examples of these terms include muscle tissue morphogenesis (GO:0060415), imaginal disc-derived appendage morphogenesis (GO:0035114), and imaginal disc development (GO:0007444). Surprisingly, we also see axon-related terms in our GO analyses such as axon guidance (GO:0007411), neuron projection guidance (GO:0097485) and neuron projection morphogenesis (GO:0048812).
The PC cluster shows expression of the pair-rule genes sob, odd and drm, which are specifically expressed in this cluster and trigger retinogenesis in the eye disc along with dpp and dac [30]. dac expression is also observed in this cluster. dac is required for proper cell fate determination in the eye imaginal disc and cells at the posterior margin of the eye disc adopt a cuticle fate in the absence of dac function [39]. Differential gene expression analyses identified CG17278 and crossveinless c (cv-c) (Supplemental Fig. 1C,D) as PC markers that are specifically expressed in this cluster. The function of these genes in PC cells and how they affect eye development is currently unknown. We also performed GO analyses using PC differentially expressed genes and, interestingly, terms related to plasma membrane extension show enrichment. These include cytoneme assembly (GO:0035231), cytoneme morphogenesis (GO:0003399) and filopodium assembly (GO:0046847). In addition, chitin-based cuticle attachment to epithelium (GO:0040005) and larval chitin-based cuticle development (GO:0008363) GO terms are also associated with genes in the PC cluster. This is expected as peripodial cells are known to form adult head cuticle [40].