Author response: Temporal evolution of single-cell transcriptomes of Drosophila olfactory projection neurons (2020)

Article Figures and data Abstract Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Neurons undergo substantial morphological and functional changes during development to form precise synaptic connections and acquire specific physiological properties. What are the underlying transcriptomic bases? Here, we obtained the single-cell transcriptomes of Drosophila olfactory projection neurons (PNs) at four developmental stages. We decoded the identity of 21 transcriptomic clusters corresponding to 20 PN types and developed methods to match transcriptomic clusters representing the same PN type across development. We discovered that PN transcriptomes reflect unique biological processes unfolding at each stage—neurite growth and pruning during metamorphosis at an early pupal stage; peaked transcriptomic diversity during olfactory circuit assembly at mid-pupal stages; and neuronal signaling in adults. At early developmental stages, PN types with adjacent birth order share similar transcriptomes. Together, our work reveals principles of cellular diversity during brain development and provides a resource for future studies of neural development in PNs and other neuronal types. Introduction Cell-type diversity and connection specificity between neurons are the basis of information processing underlying all nervous system functions. The precise assembly of neural circuits involves multiple highly regulated steps. First, neurons are born from their progenitors and acquire unique fates through a combination of (1) intrinsic mechanisms, such as lineage, birth order, and birth timing; (2) extrinsic mechanisms, such as lateral inhibition and extracellular induction; and (3) developmental stochasticity in some cases (Jan and Jan, 1994; Johnston and Desplan, 2010; Kohwi and Doe, 2013; Holguera and Desplan, 2018; Li et al., 2018). During wiring, neurons extend their neurites to a coarse targeting region, elaborate their terminal structures, select pre- and post-synaptic partners, and finally form synaptic connections (Sanes and Yamagata, 2009; Jan and Jan, 2010; Kolodkin and Tessier-Lavigne, 2011; Luo, 2020; Sanes and Zipursky, 2020). Studies from the past few decades have uncovered many molecules and mechanisms that regulate each of these developmental processes. The development of Drosophila olfactory projection neurons (PNs) has been extensively studied (Jefferis et al., 2004; Hong and Luo, 2014). PNs are the second-order olfactory neurons that receive organized input from olfactory receptor neurons (ORNs) at ~50 stereotyped and individually identifiable glomeruli in the antennal lobe, and carry olfactory information to higher brain centers (Vosshall and Stocker, 2007; Wilson, 2013; Figure 1A). Different types of PNs send their dendrites to a single glomerulus or multiple glomeruli (Marin et al., 2002; Lai et al., 2008; Yu et al., 2010; Tanaka et al., 2012; Bates et al., 2020). PNs are derived from three separate neuroblast lineages—anterodorsal, lateral, and ventral lineages, corresponding to their cell bodies' positions relative to the antennal lobe (Jefferis et al., 2001). PNs produced from the anterodorsal and lateral lineages (adPNs and lPNs) are cholinergic excitatory neurons. The fate of uniglomerular excitatory PN types, defined by their glomerular targets, is predetermined by their lineage and birth order (Jefferis et al., 2001; Marin et al., 2005; Yu et al., 2010; Lin et al., 2012). PNs produced from the ventral lineage (vPNs), on the other hand, are GABAergic inhibitory neurons (Jefferis et al., 2007; Liang et al., 2013; Parnas et al., 2013). The connectivity and physiology of PNs have also been systematically studied (Bhandawat et al., 2007; Jeanne et al., 2018; Bates et al., 2020). Figure 1 with 1 supplement see all Download asset Open asset Overview of single-cell transcriptomic profiling of Drosophila olfactory projection neurons (PNs). (A) Schematic of the adult Drosophila olfactory system. Approximately 50 types of olfactory receptor neurons (ORNs) form one-to-one synaptic connections with 50 types of excitatory PNs at 50 glomeruli in the antennal lobe. Illustrated are two types each of ORNs (brown) and PNs (green), as well as two glomeruli to which their axons and dendrites target. (B) Schematic of the developmental process of the adult Drosophila olfactory system. The ~50 types of uniglomerular excitatory PNs are from either anterodorsal (adPN) or lateral (lPN) neuroblast lineages. PNs with cell body on the ventral side are inhibitory ventral PNs (vPNs). (C) Representative confocal images of PNs from four different developmental stages, 0 hr APF, 24 hr APF, 48 hr APF, and adult. APF: after puparium formation. Images are shown as maximum z-projections of confocal stacks. Antennal lobes are outlined. Scale bars, 40 μm. (D) Workflow of the single-cell RNA sequencing using plate-based SMART-seq2. FACS: fluorescence-activated cell sorting. (E) Summary of the number of high-quality PNs sequenced at each timepoint and driver lines used. Most PNs refer to PNs sequenced using either GH146-GAL4 or VT033006-GAL4. (F) Visualization of all sequenced PNs from four different developmental stages using tSNE plot. Dimensionality reduction was performed using the top 500 overdispersed genes identified from all sequenced PNs. Despite the fact that PNs are among the most well-characterized cell types in all nervous systems, their transcriptome-wide gene expression changes across different developmental stages with cell-type specificity are still unknown. This information can help us obtain a more complete picture of both known and unexplored pathways underlying neural development and function. Recently, the advent of single-cell RNA sequencing (scRNA-seq) has paved the way toward obtaining such data (Li et al., 2017; Kalish et al., 2018; Zhong et al., 2018; Li, 2020). Here, we profiled and analyzed the single-cell transcriptomes of most uniglomerular excitatory PNs. We identified the correspondence between two-thirds of transcriptomes and PN types at one stage and developed methods to reliably match transcriptomic clusters corresponding to the same types of PNs across different stages. We discovered that PN transcriptomes exhibit unique characteristics at different stages, including birth-order, neurite pruning, wiring specificity, and neuronal signaling. The identification of many differentially expressed genes among different PN types, such as transcription factors, cell-surface molecules, ion channels, and neurotransmitter receptors, provides a rich resource for further investigations of the development and function of the olfactory system. Results Single-cell transcriptomic profiling of Drosophila PNs at four developmental stages The development of PNs follows the coordinated steps as previously described (Hong and Luo, 2014). Eighteen out of 40 types of adPNs are born embryonically and participate in the larval olfactory system. Then, during the larval stage, the rest of adPNs and all lPNs are born (Jefferis et al., 2001; Marin et al., 2005; Yu et al., 2010; Lin et al., 2012). During early metamorphosis following puparium formation, embryonically born PNs first prune terminal branches of dendrites and axons, and then re-extend their dendrites into the future adult antennal lobe, and axons into the mushroom body and lateral horn following the neurites of larval-born PNs (Marin et al., 2005). From 0 to 24 hr after puparium formation (APF), PNs extend their dendrites into the developing antennal lobe and occupy restricted regions. ORN axons begin to invade antennal lobe at ~24 hr APF. PN dendrites and ORN axons then match with their respective partners beginning at ~30 hr APF and establish discrete glomerular compartments at ~48 hr APF. Thereafter, they expand their terminal branches, build synaptic connections, and finally form mature adult olfactory circuits (Jefferis et al., 2004; Figure 1B). To better understand the molecular mechanisms that control these dynamic developmental processes underlying neural circuit assembly, we performed scRNA-seq of PNs from four different developmental stages: 0–6 hr APF, 24–30 hr APF, 48–54 hr APF, and 1–5 days adult (hereafter 0, 24, 48 hr APF and adult) (Figure 1C). We used GH146-GAL4 (Stocker et al., 1997) to drive UAS-mCD8-GFP (Lee and Luo, 1999) expression in most PNs at 24 hr, 48 hr, and adult, which labels ~90 of the estimated 150 PNs in each hemisphere, covering ~40 of the 50 PN types. At 0 hr APF, GH146-GAL4 also labels cells in the optic lobes (Figure 1—figure supplement 1A), which are inseparable from the central brain by dissection. Therefore, we used VT033006-GAL4 to label PNs at 0 hr APF (Figure 1C and Figure 1—figure supplement 1B; Tirian and Dickson, 2017). VT033006-GAL4 labels most PNs from the anterodrosal and lateral lineage, but not PNs from the ventral lineage or anterior paired lateral (APL) neurons like GH146-GAL4. It is expressed in ~95 cells that innervate ~44 glomeruli which largely overlap with PNs labeled by GH146-GAL4 (Inada et al., 2017; Elkahlah et al., 2020). In addition to PNs labeled by GH146-GAL4 and VT033006-GAL4 (we will refer to them as 'most PNs' hereafter), we have collected single-cell transcriptomic data using drivers that only label a small number of PN types for mapping the transcriptomic clusters to anatomically defined PN types. For scRNA-seq, fly brains with a unique set of PN types labeled using different drivers at each developmental stage were dissected and dissociated into single-cell suspensions. GFP+ cells were sorted into 384-well plates by fluorescence-activated cell sorting (FACS), and sequenced using SMART-seq2 (Picelli et al., 2014; Figure 1D) to a depth of ~1 million reads per cell (Figure 1—figure supplement 1C). On average ~3000 genes were detected per cell (Figure 1—figure supplement 1D), and after quality filtering (see Materials and methods), we obtained ~3700 high quality PNs in addition to the previously sequenced ~1200 PNs (Li et al., 2017), yielding ~4900 PNs for analysis in this study (Figure 1E). All analyzed PNs express high levels of neuronal markers but not glial markers, confirming the specificity of sequenced cells (Figure 1—figure supplement 1E). Unbiased clustering using overdispersed genes from all PNs readily separates them into different groups according to their stage (Figure 1F), suggesting that gene expression changes across these four developmental stages represent a principal difference in their single-cell transcriptomes. Decoding the glomerular identity of transcriptomic clusters by sequencing subsets of PNs at 24 hr APF PNs labeled by GH146-GAL4 at 24 hr APF form ~30 distinct transcriptomic clusters. We previously matched six of these transcriptomic clusters to specific anatomically and functionally defined PN types (Li et al., 2017), hereafter referred to as 'decoding transcriptomic identity'. Unlike ORNs, whose identities can be decoded using uniquely expressed olfactory receptors (Li et al., 2020a), PNs lack known type-specific markers. Instead, PN types are mostly specified by combinatorial expression of several genes (Li et al., 2017), making it more challenging to decode their transcriptomic identities. To circumvent these challenges and decode the transcriptomic identities of more types of PNs, we took advantage of the extensive driver line collection in Drosophila (Luan et al., 2006; Jenett et al., 2012; Dionne et al., 2018). We searched for split-GAL4 lines that only labeled a small proportion of all PNs (Yoshi Aso, unpublished data). Using such drivers, we could sequence a few types of PNs at a time, plot those cells with most PNs, and then use differentially expressed markers among them to decode their identities one-by-one. split#28 GAL4 labeled two types of PNs—those that project their dendrites to the DC3 and DA4l glomeruli in developing and adult animals (Figure 2A,B; note that PN types are named after the glomeruli they project their dendrites to). We sequenced those PNs (split#28+ PNs hereafter) at 24 hr APF. We chose this stage because this is when different PN types exhibit the highest transcriptome diversity as hinted by the number of clusters seen in Figure 1F (see following sections for more detailed analysis). To visualize sequenced split#28+ PNs, we performed dimensionality reduction using 561 genes identified from most 24 hr PNs using Iterative Clustering for Identifying Markers (ICIM), an unsupervised machine learning algorithm (Li et al., 2017), followed by embedding in the tSNE space. Split#28+ PNs (orange dots) fell into two distinct clusters and intermingled with GH146+ PNs (gray dots) (Figure 2C). One cluster mapped to previously decoded DC3 PNs (Li et al., 2017), and the other cluster expressed zfh2 (Figure 2—figure supplement 1A). We validated that this cluster indeed represents DA4l PNs by visualizing the expression of zfh2 in PNs utilizing an intersectional strategy by combining zfh2-GAL4, GH146-Flp, and UAS-FRT-STOP-FRT-mCD8-GFP (hereafter referred to as 'intersecting with GH146-Flp') (Figure 2—figure supplement 1B). Figure 2 with 3 supplements see all Download asset Open asset Matching 15 transcriptomic clusters to specific projection neuron (PN) types at 24 hr APF. (A) Representative maximum z-projection of confocal stacks of split#28 GAL4 in adults. Dendrites of split#28 GAL4+ PNs target the DC3 and DA4l glomeruli. (B) Diagram of split#28-GAL4+ PNs. (C) tSNE plot showing newly sequenced split#28-GAL4+ PNs, which form two clusters that can be assigned to DC3 and DA4l PNs (see also Figure 2—figure supplement 1). (D) Representative confocal images of split#7 GAL4 labeled PNs using permanent labeling strategy. One anterior section and one posterior section of the antennal lobe are shown. Using permanent labeling, we found that this driver is expressed in eight PN types. Genotype: split#7-GAL4, UAS-Flp, Actin promoter-FRT-STOP-FRT-GAL4, UAS-mCD8-GFP. (E) Diagram of split#7-GAL4+ PNs. split#7 GAL4 labels eight types of PNs. four from the adPN lineage (green letters) and four from the lPN lineage (red letters). (F) tSNE plot of split#7 GAL4 PNs with GH146+ PNs (see Figure 2—figure supplement 2 for details on the decoding procedure). (G) Representative maximum z-projection of confocal stacks of kn+ PNs in the adult. kn-GAL4 was intersected with GH146-Flp to restrict the expression of GAL4 in only PNs. (H) Representative confocal images of split#15 GAL4 in adults, which labels two kn+ PN types. (I) Diagram showing that kn+ PNs include six types of adPNs and two vPNs. (J) tSNE plot of kn-GAL4 PNs with GH146+ PNs (see Figure 2—figure supplement 3 for details on the decoding procedure). (K) Dot plot summarizing drivers and marker genes we used to map 21 transcriptomic clusters to 20 PN types [14 adPNs, 5 lPNs—DA1 PNs form two clusters, one fru+ and one fru– (Li et al., 2017)—and 1 vPNs] and the anterior paired lateral (APL) neurons at 24 hr APF. Gene expression level [log2(CPM+1)] is shown by the dot color, and percentages of cells expressing a marker are shown by dot size. (L) tSNE plot showing 24 hr APF PNs colored by PN types (GH146+ PNs with split#7+/ split#28 PNs to increase cell number in some less abundant PN types). Scale bars, 20 μm. Axes, D (dorsal), L (lateral). In panel B, E, and I, orange glomeruli represent PN types of unknown transcriptomic identity prior to this study. Green glomeruli represent PN types whose transcriptomic identity were previously decoded. Note that the positions of cells on a tSNE plot are dependent on the random initialization of the program as well as every cell present in the dataset, therefore the position of GH146+ PNs clusters are different when we plot them with different set of newly sequenced PNs (gray in panels C, F, and J). split#7 GAL4 labeled three types of PNs in the adult stage (Figure 2—figure supplement 2A). However, when we sequenced cells labeled by this GAL4 line at 24 hr APF and visualized them using tSNE, we found eight distinct clusters (Figure 2F). We reasoned that this could be due to loss of driver expression in adult stage for some PN types. To test this hypothesis and reveal PNs that are labeled by this driver transiently during development, we used a permanent labeling strategy to label all cells that express split#7 GAL4 at any time of development (split#7+ PNs hereafter) by combining it with UAS-mCD8-GFP, Actin promoter-FRT-STOP-FRT-GAL4, and UAS-Flp. Using this strategy, we observed labeling of eight types of PNs (Figure 2D), consistent with number of clusters we observed by sequencing. Among split#7+ PNs, four types belong to the adPN lineage (acj6+) and the other four types belong to the lPN lineage (vvl+) (Figure 2E). Only one lPN type, DA1 (CG31676+), has previously been decoded (Figure 2—figure supplement 2B). We identified differentially expressed genes among split#7+ PNs and obtained existing GAL4 lines mimicking their expression. By intersecting those GAL4 lines with GH146-Flp, we mapped all seven previously unknown transcriptomic clusters to seven PN types (Figure 2—figure supplement 2C-H; see legends for detailed description). In addition to screening through collections of existing driver lines, we also utilized scRNA-seq data to find drivers that label a subpopulation of PNs. One such marker was the gene knot (kn), which was expressed in seven transcriptomic clusters among all GH146+ PNs (Figure 2—figure supplement 3A). One of the kn+ clusters expressing trol has been previously mapped to VM2 PNs (Li et al., 2017). When kn-GAL4 was intersected with GH146-Flp, six types of adPNs (acj6+) and several vPNs (Lim1+) were labeled (Figure 2G,J). Among the six adPN types, VM7 and VM5v PNs were also labeled by split#15 GAL4 (Figure 2H). Although it has been previously reported that GH146-GAL4 is not expressed in VM5v PNs (Yu et al., 2010), labeling of these PNs when GH146-Flp was intersected with either kn-GAL4 or split#15 GAL4 indicates that GH146-Flp must be expressed in VM5v PNs at some point during development. Using split#15 GAL4, we were able to decode the two clusters to be either VM7 or VM5v PNs (Figure 2—figure supplement 3B). Due to the lack of existing GAL4 drivers for differentially expressed genes between these two clusters, we could not further distinguish them so far; their identities can be decoded by creating new GAL4 drivers in future studies. Other than these two clusters, we were able to match transcriptomic clusters and glomerular types for the rest of kn+ adPNs one-to-one (Figure 2—figure supplement 3C–E). In addition to excitatory PNs, one kn+ vPN type innervated DA1 glomerulus (because DA1 glomerulus is innervated only by lPNs and vPNs, not adPNs). We found that DIP-beta was expressed in one kn+ vPN cluster but not in lPNs innervating DA1 glomerulus (Figure 2—figure supplement 3F,G). Intersecting DIP-beta-GAL4 with GH146-Flp confirmed that DIP-beta+ vPN indeed targeted their dendrites to DA1 glomerulus, illustrating the DIP-beta+ vPN cluster to be DA1 vPNs (Figure 2—figure supplement 3H). In summary, by sequencing a small number of known PN types at a time and analyzing the expression pattern of differentially expressed genes, we have now mapped a total of 21 transcriptomic clusters corresponding to anatomically defined PN types at 24 hr APF (Figure 2K,L). Ultimately, we aimed to match the transcriptomes of the same PN types across development. As an intermediate step, we carried out global analysis of gene expression changes across development, which could help us reliably identify transcriptomic clusters representing different PN types at different developmental stages. Global gene expression dynamics across four developmental stages All sequenced PNs segregated into different clusters according to their developmental stages using unbiased, over-dispersed genes for clustering regardless of PN types (Figure 1F). Even when we used the genes identified by ICIM for clustering, which emphasizes the differences between different PN types (Li et al., 2017), we still observed that individual PNs were separated principally by developmental stages (Figure 3A). Together, these observations illustrate global transcriptome changes of PNs from pupa to adult. Figure 3 Download asset Open asset Global gene expression dynamics of projection neurons (PNs). (A) Visualization of PNs from four different developmental stages: 0 hr APF, 24 hr APF, 48 hr APF, and adult sequenced using either VT033006-GAL4 or GH146-GAL4. tSNE dimensionality reduction was performed using 1216 genes identified by iterative clustering for identifying markers (ICIM) among them. (B) Hierarchical heatmap showing the expression of the top 52 out of 474 differentially expressed genes identified among PNs of different developmental stages. (C) Examples of the expression of the dynamic genes. Cells are colored according to the expression level of each gene. Akap200 (A kinase anchor protein 200, encodes a scaffolding protein that contributes to the maintenance and regulation of cytoskeletal structure), cib (ciboulot, encodes an actin binding protein), and fax (failed axon connections, a gene involved in axon development) have the highest expression in early pupal stage and are downregulated gradually. Rdl (Resistant to dieldrin, encodes a chloride channel), slo (slowpoke, encodes a subunit of calcium-activated channel), and are as PNs 474 differentially expressed genes can be into eight groups on their dynamic groups developmental groups of genes and two groups of genes lines represent individual genes and the line expression of genes in each The highest expression is as one for all genes. The top Gene for and downregulated genes are shown on To understand types of genes drive this we searched for genes that were differentially expressed in different developmental stages (Figure We the genes into different groups on their expression pattern development. groups of genes developmental groups were from pupa to adult and two groups were (Figure with our neural genes, including those with in and were in developing PNs (Figure genes to synaptic ion and on the other hand, were in mature PNs (Figure Li et al., 2017; Li et al., Single-cell transcriptomes of PNs reveal biological processes at different stages of development PN transcriptomes global we to find a to reliably and transcriptomic clusters representing different PN types at all stages. We first identified genes for clustering from each stage using ICIM and used them for further dimensionality However, using this we obtained different of clusters at each stage (Figure of each stage unique biological of PN development. Figure with 2 supplements see all Download asset Open asset PN transcriptomes distinct at different stages of development. (A) Visualization of most projection neurons (PNs) from 0 hr APF, 24 hr APF, 48 hr APF, and using tSNE on genes identified by ICIM at each clusters are identified using (B) Clustering of 0 hr APF PNs using identified two clusters. (C) of the molecular pathways for neurite pruning in (D) whose function have been in neurite pruning have higher expression in cluster and PN transcriptomic at the cell level using 1216 ICIM genes identified from PNs of all four and the cluster level using differentially expressed genes identified from 24 hr PN for adPNs (E) hr APF: 15 clusters, 24 hr APF: 15 clusters, 48 hr APF: clusters, adult 15 clusters, and lPNs (F) hr APF: clusters, 24 hr APF: clusters, 48 hr APF: clusters, adult eight clusters, (G) Schematic of PN transcriptome changes from early pupal stage to PN diversity during circuit assembly 24 hr APF and as they into mature neurons. (H) of and in adult PNs. At 0 hr APF, PNs two distinct cluster of both adPNs and and a one with only adPNs (Figure Figure supplement 2A). As all lPNs and many adPNs are born during the larval stage, some adPNs are born during the We that the cluster could represent embryonically born PNs, which undergo axon and pruning during early metamorphosis (Marin et al., 2005). pruning in Drosophila on cell (Lee et al., of the receptor et al., et al., et al., binding of the and form a to a of targets, including a transcription which in expression of the cytoskeletal and (Figure et al., et al., 2009; et al., 2011; et al., 2013). To test our we the expression of genes which are known to participate in neurite pruning and genes that expression in the mushroom body neurons during pruning et al., 2018). We found that and two sorting for and indeed higher expression levels in the cluster (Figure We also confirmed our hypothesis by mapping two types of embryonically born PNs, DA4l and PNs, to this cluster (Figure supplement see mapping details in Figure At 24 hr APF, we observed the highest number of clusters different PN types. dimensionality reduction using the top overdispersed genes also more distinct clusters at this timepoint to the (Figure supplement 1). of transcriptomic among PNs at each stage indeed confirmed the highest diversity among PNs at 24 hr APF (Figure This is by the fact that at this stage, PNs their dendrites to specific and begin to as for their ORN In PN axons at the lateral horn begin to establish their (Jefferis et al., All these processes high level of molecular diversity among different PN types to precise wiring, more between their transcriptomes at this In to the high transcriptomic diversity in 24 hr APF PNs, adult PNs only three clusters (Figure by The three clusters represent cholinergic excitatory PNs by and two GABAergic inhibitory cell and neurons (Figure This is because after wiring specificity is all excitatory PNs similar but distinct from inhibitory neuronal types. at different developmental stages, the differentially expressed genes we identified all the most biological processes those neurons are observations that PN transcriptomes reflect the pruning process of embryonically born PNs at 0 hr APF, PN type and wiring at 24 hr APF, and neurotransmitter type in adults. Identifying PN types at all developmental stages the of the 24 hr APF PNs, gene identified from each of the other stages could not distinct clusters PN type diversity (Figure Therefore, we to use the genes identified by ICIM from 24 hr APF PNs to cluster PNs of the other stages. We found that this gene set all other gene in different PN types at all (Figure In most gene found by different methods at 24 hr APF, including