Nuclear Transcriptomes of the Seven Neuronal Cell Types That Constitute the Drosophila Mushroom Bodies

The insect mushroom body (MB) is a conserved brain structure that plays key roles in a diverse array of behaviors. The Drosophila melanogaster MB is the primary invertebrate model of neural circuits related to memory formation and storage, and its development, morphology, wiring, and function has been extensively studied. MBs consist of intrinsic Kenyon Cells that are divided into three major neuron classes (γ, α′/β′ and α/β) and 7 cell subtypes (γd, γm, α′/β′ap, α′/β′m, α/βp, α/βs and α/βc) based on their birth order, morphology, and connectivity. These subtypes play distinct roles in memory processing, however the underlying transcriptional differences are unknown. Here, we used RNA sequencing (RNA-seq) to profile the nuclear transcriptomes of each MB neuronal cell subtypes. We identified 350 MB class- or subtype-specific genes, including the widely used α/β class marker Fas2 and the α′/β′ class marker trio. Immunostaining corroborates the RNA-seq measurements at the protein level for several cases. Importantly, our data provide a full accounting of the neurotransmitter receptors, transporters, neurotransmitter biosynthetic enzymes, neuropeptides, and neuropeptide receptors expressed within each of these cell types. This high-quality, cell type-level transcriptome catalog for the Drosophila MB provides a valuable resource for the fly neuroscience community.

, and also play important roles in other behavioral contexts such as temperature preferences (Hong et al. 2008), sleep (Artiushin and Sehgal 2017) and responses to ethanol exposure (Kaun et al. 2011).
MB dependent plasticity is one of the most intensely studied aspects of invertebrate neurobiology. The morphology and developmental lineage of the neurons that populate the MB in Drosophila, as well as the identity and morphology of most of their neuronal inputs and outputs, have been fully characterized Jefferis et al. 2002;Aso et al. 2014a;2014b). Many functional manipulations of both neural activity and signaling pathways relevant to plasticity have been conducted within each of the identified neuronal cell types in this circuit (Connolly et al. 1996;Zars et al. 2000;Dubnau et al. 2001;McGuire et al. 2001;Isabel et al. 2004;Krashes et al. 2007;Blum et al. 2009;Trannoy et al. 2011;Qin et al. 2012;Huang et al. 2012;Cervantes-Sandoval et al. 2013;Perisse et al. 2013;Bouzaiane et al. 2015). Functional imaging studies have established neural activity correlates in behaving animals (Davis 2011). Together, these studies support the conclusion that the neurons of the MB play unique roles in memory acquisition, storage and retrieval. Moreover, memory storage over the course of minutes and hours after training relies on an evolving requirement for reverberating neural activity within a circuit that includes MB intrinsic neurons and the so-called extrinsic neurons that provide inputs and outputs (Dubnau and Chiang 2013;Cognigni et al. 2018). In contrast to the increasingly deep understanding of the development, connectivity and functional requirements of each cell type in this circuit, there is a surprising paucity of data on differences in their transcriptional profiles.
MB KCs can be divided into three major classes, g, a9/b9 and a/b, based on the projection patterns of the axons (Crittenden et al. 1998;Lee et al. 1999). Extensive anatomical and functional characterization corroborates this classification as biologically relevant Davis 2011;Dubnau and Chiang 2013). Gene expression differences between these three classes of MB KCs have been studied using microarray (Perrat et al. 2013), RNA sequencing (RNA-seq) (Crocker et al. 2016) and single-cell RNA-seq (Croset et al. 2018;Davie et al. 2018). However, it is known that the three KC classes can be further separated into seven subtypes: gd, gm, a9/b9ap, a9/b9m, a/bp, a/bs and a/bc KCs. The functional relevance of this further subdivision is supported by expression of split-GAL4 lines, analysis of the axonal projection patterns of individual neurons from each cell subtype (Luan et al. 2006;Aso et al. 2014a) and investigation of their functional roles in behavior (Perisse et al. 2013;Sitaraman et al. 2015a;Vogt et al. 2016;Yang et al. 2016). Until now, there have been no attempts to identify the unique transcriptional programs that control/establish the identity of each of these seven cell subtypes, while cell clustering using a resource of enormous single-cell transcriptional profiles showed only three MB KC clusters (Davie et al. 2018). But the availability of intersectional genetics approaches that make use of split-GAL4 provide the means to investigate each of the subtypes individually (Luan et al. 2006;Pfeiffer et al. 2008;Aso et al. 2014a).
Several methods have been developed to characterize the transcriptional programs of specific cell types in flies, such as TRAP (Heiman et al. 2008), TU-tagging (Miller et al. 2009), and EC-tagging (Hida et al. 2017). Here, we used an improved version of the INTACT method (Henry et al. 2012), called tandem affinity purification of intact nuclei and RNA-sequencing (TAPIN-seq; Davis et al. 2018), to profile the nuclear transcriptomes of all seven MB neuronal subtypes that constitute the MB (Aso et al. 2014a). These transcriptomes revealed 350 genes with either class-or subtype-specific expression, including several well-known and many new class or subtype markers. Moreover, our data provide a full accounting of the input-output signaling properties for each of these neuron subtypes including neurotransmitter biosynthetic machinery, neuropeptides and neurotransmitter and neuropeptide receptors.

Fly stocks
Because quantitative traits such as gene expression profiles can be sensitive to genetic background, we created a newly isogenic subline of w 1118 (isoCJ1), which was itself derived from Canton-S wild type as an inbred line many years ago (Yin et al. 1994). To generate the isogenic strain, we used ten-generations of single male and female sibling crosses to generate 10 independent isogenic strains. MJ2 was selected based on its ability to form comparable olfactory short-term memory performance to the parental strain in the standard Pavlovian task (Supplemental Material, Figure S1). The nuclear envelope epitope tagged transgene P{5XUAS-unc84::2XGFP}attP40, the pJFRC28 strain P {10XUAS-IVS-GFP-p10}attP2 and each of the split-GAL4 inserts were backcrossed into this new MJ2 wild type strain for five generations to equilibrate each to this isogenic background. For each split-GAL4 combination, we separately backcrossed the GAL4 activating domain and DNA-binding domain components, and then combined the two hemidrivers as a split-GAL4 line in the MJ2 background thereafter using standard balancer chromosomes that had themselves been equilibrated to the MJ2 strain. The UAS-WM P{5XUAS-myr::GFP-V5-p2A-His2B:: mCherry-HA} reporter strain was generated using standard approaches (Chang et al. 2018). Flies were cultured on standard cornmeal food using the standard cornmeal recipe from Bloomington stock center. Food was supplemented with antibiotics.
For imaging to characterize expression with each split-GAL4 strain that had been reconstituted in the MJ2 background, we used 2 -5 day old male flies. For RNA-seq sample preparation, each split-GAL4 line was crossed to the P{5XUAS-unc84::2XGFP}attP40. 2 -5 day old adult progeny for each genotype were collected and frozen in liquid nitrogen between 10 AM and 7 PM.
Released chromatin and broken nuclei were adsorbed to carboxyl coated magnetic beads for 30 min at 4°with constant rotation. Unbound antibody was removed by incubating the sample on ice for 20 min with 100 mL of washed UNOsphere SUPra resin (Bio-Rad: 1560218). After the resin was removed on a 10 mm cup filter and the carboxyl beads on a magnet stand, the nuclei-containing supernatant was mixed with an equal volume of 500 mM sodium acetate pH 8.5, 250 mM sucrose, 6 mM EGTA, 6 mM EDTA, 0.6 mM spermidine, 0.2 mM spermidine, 1 mM DTT, 1· complete protease inhibitor, 0.25 mg/mL torula RNA and 30 mL Protein A Dynabeads (ThermoFisher: 10002D). A 2-hour incubation on ice with occasional agitation was used to recover tagged nuclei. Bead-bound nuclei were then recovered on a magnet stand and washed twice with 250 mM sodium acetate pH 8.5, 250 mM sucrose and 0.1% NP-40. Nuclei were then released at 37°for 1 hr by incubation in 50 mL of 10 mM Tris pH 7.5, 2.5 mM MgCl 2 , 0.5 mM CaCl 2 , 250 mM sucrose, 0.1% NP-40, 1 mg/mL torula RNA, 40 units RNAsin (Promega: N2515), 2 units DNAseI (NEB: M0303L), 320 units IdeZ protease (NEB: P0770S). The sample was diluted to 100 mL with 10 mM Tris pH 7.5, 2.5 mM MgCl 2 , 0.5 mM CaCl 2 , 250 mM sucrose and 0.1% NP-40, EGTA was added to 1 mM and the suspension was rapidly triturated 100 times. After returning the sample to a magnet stand, 90 mL of buffer containing released nuclei was removed and added to 1.5 mL of Protein G Dynabeads (ThermoFisher: 10004D) that were previously resuspended in 10 mL of 10 mM Tris pH 7.5, 2.5 mM MgCl 2 , 0.5 mM CaCl 2 , 250 mM sucrose and 0.1% NP-40. The second binding reaction was run for 1 -3 hr on ice with occasional agitation, followed by 2· 250 mL washes in 10 mM Tris pH 7.5, 2.5 mM MgCl 2 , 0.5 mM CaCl 2 , 250 mM sucrose and 0.1% NP-40. Prior to the last wash a 5 mL aliquot was removed for quantitation and the remainder of the sample was solubilized in Arcturus Picopure RNA extraction buffer (ThermoFisher: KIT0204).

RNA-seq library construction
Nuclear RNA was DNAseI treated and purified using the Arcturus PicoPure system exactly as instructed by the supplier. Purified RNA was mixed with a 1:100,000 dilution of ERCC standard RNA (ThermoFisher: 4456740) and amplified using the Nugen Ovation v2 system (Nugen: 7102-32). cDNA was then blunted, ligated to barcoded linkers (Nugen: 0319-32, 0320-32) and sequenced in paired-end mode on an Illumina HiSeq2500 to 125 nt read lengths.

RNA-seq data analysis
We trimmed RNA-seq reads (5nt from the 59 end of the forward read, using seqtk option "trim -b 5"; https://github.com/lh3/seqtk) to remove non-transcript sequences introduced by the NuGen Ovation kit and then pseudo-aligned these to the Drosophila transcriptome (ENSEMBL release 91, BDGP6) using kallisto (Bray et al. 2016) to estimate transcript abundances. We also included sequences for the synthetic ERCC spike-in species and TAPIN reporter in the transcriptome index. After pseudo-alignment, we removed ERCC, TAPIN reporter, and ribosomal RNA entries and renormalized the transcript abundance matrix to units of transcripts per million (TPM). To visualize TAPIN-seq signal across the genome, we also aligned trimmed reads to the whole genome (BDGP6, dm6) using STAR (Dobin et al. 2013), created bigWig genome tracks (deeptools; Ramírez et al. 2016), and visualized them in the IGV genome browser (Robinson et al. 2011).
To identify class-and subtype-enriched genes, we performed differential expression analysis using the estimated counts from kallisto as input to limma (Ritchie et al. 2015), voom (Law et al. 2014), and quantile normalizing the expression levels to account for differences in the number of genes detected in each sample (Table S1). We used criteria of at least 10 TPM abundance in one sample, at least twofold difference in expression, and 5% false discovery rate to identify differentially expressed genes.

Immunohistochemistry
Immunohistochemistry was performed essentially as in a previous report . Brains were dissected in isotonic PBS and immediately transferred to 4% paraformaldehyde in PBS for a 30-min fixation at room temperature. Fixed brain samples, were rinsed with isotonic PBS and incubated in PBS containing 2% Triton X-100, 10% normal goat serum (NGS; Penetration & Blocking Buffer) while being subjected to a degassing procedure (Chiang et al. 2011). Brain samples were agitated in the same buffer at 4°overnight. Brains were then transferred to primary antibodies diluted with PBS containing 0.25% Triton X-100, 1% NGS (Dilution Buffer) and agitated at 4°for 1-3 day. After primary antibody incubation, the brain samples were washed in PBS containing 1% Triton X-100, 3% NaCl (Washing Buffer) three times before they were moved to the 1:250 diluted secondary antibodies for one day agitation at 4°. For the biotin amplification staining (

Confocal imaging and post-imaging processing
Brain samples were imaged under a Zeiss LSM 800 confocal microscope with a 40X C-Apochromat water-immersion objective lens (N.A. value 1.2). The settings for scanning were manually adjusted. To overcome the limited field of view when imaging the GAL4 expression patterns, we scanned each brain in two parallel stacks of confocal images with some overlap between the two brain hemispheres, with a voxel size of 0.31 X 0.31 X 1.25 mm. We then stitched the two image stacks into a single data set with 'Pairwise stitching' function in Fiji (Preibisch et al. 2009;Schindelin et al. 2012), segmented the brain region based on the dlg1 staining channels with 3D Slicer (https://www.slicer.org; Kikinis et al. 2013), and made a 'Z projection' with Fiji (Schindelin et al. 2012). The MB subtype models were constructed from the GFP channel of the confocal images used for projections, by using the 3D Slicer to segment, show 3D, and conduct smoothing.

Cell counting
The confocal images for cell counting were acquired with a voxel size of 0.16 X 0.16 X 1.00 mm. The GFP channel was first used to identify KCs, and then the 'Cell Counter' plugin in Fiji was used to count all detectable nuclei in the mCherry channel (Preibisch et al. 2009). For each line, we counted three hemispheres from three different animals.

Data availability
All Drosophila strains are available upon request. Supplemental files are available at FigShare. Files S1 -S7 contain the expression pattern for each of the split-GAL4 drivers, while Files S8 -S14 contain one example of the cognate WM images used for cell counting. Figure S1 contains the olfactory short-term memory performance for 10 newly generated isogenic strains. Figure S2 contains overview of the workflow of TAPIN purification of nuclei and the molecular size distribution of amplified cDNA obtained from the 14 TAPIN-seq libraries. Figure S3 and S4 contains the view of the cell bodies of MB KCs in the confocal images that have double labeling for immunostaining and GAL4 expression. Table S1 contains the TAPIN-seq library statistics, including the numbers of raw reads, pseudoaligned reads and uniquely aligned reads, as well as the detected gene numbers. Table S2 contains the full list of genes enriched and depleted in individual MB classes and subtypes. The sequencing reads and processed data files, including the tables of estimated abundances, are available in NCBI GEO (GSE119629). All code used to analyze RNA-seq results and create the figures and tables in this manuscript are available in the GitHub repository (http://github.com/fredpdavis/ mushroombody). Supplemental material available at Figshare: https:// doi.org/10.25387/g3.7267481.

TAPIN-seq profiling of MB neuronal cell subtypes
To label MB subtypes, we used seven split-GAL4 lines: MB607B (gd), MB131B (gm+d), MB370B (a9/b9ap+m), MB418B (a9/b9m), MB371B (a/bp), MB185B (a/bs) and MB594B (a/bc) (Aso et al. 2014a). We first backcrossed all the split GAL4 hemi-drivers and the nuclear-tag reporter (Henry et al. 2012) into MJ2, an isogenic Canton-S derivative (Methods). Because of the change in genetic background, we re-characterized the expression pattern of each split-GAL4 combination to confirm that the expected cell subtype-specific pattern had not been altered. We used a novel membrane-GFP-P2A-nuclear-mCherry dual label reporter (Watermelon or WM for short; see Methods) (Chang et al. 2018), which expresses membrane-tethered GFP and nuclear mCherry from a single transcript via a viral ribosome skip peptide (P2A; Daniels et al. 2014). The GFP label revealed neuronal morphology, thereby confirming MB subtype specificity, and the nuclear-mCherry marked each nucleus, which permitted an accurate cell count. Using WM labeling, we found that each split-GAL4 combination yielded limited expression in only small numbers of neurons outside the MB, but exhibited strong expression in the annotated MB cell subtype ( Figures 1A -1G -S14), we were able to count the total number of MB KCs of each subtype that is labeled by a given spilt-GAL4 combination ( Table 1). The numbers of labeled cells for each MB KC subtype are in general agreement with a previous report (Aso et al. 2014a) and it appears that most if not all of the neurons of each subtype are labeled with these Split-GAL4 combinations. In fact the total number of KCs labeled by five non-overlapping split-GAL4 lines -MB131B, MB370B, MB371B, MB185B and MB594Bis very close to the estimated total number of KCs (1,855 vs. 2,000; Aso et al. 2014a), making it unlikely that any major KC subtype is missed. Together, these results confirm the previously reported specificity and comprehensiveness of these split-GAL4 lines to label each of the MB KC neuronal subtypes.
To profile the nuclear transcriptome of each MB subtype, we used TAPIN-seq (Davis et al. 2018), a modification of INTACT (Henry et al. 2012) that yields improved selectivity by use of a two-step purification ( Figure S2A). The 7 characterized split-GAL4 combinations were used to express the nuclear membrane protein UNC84 fused with 2 copies of GFP ( Figure 2A). For each split-GAL4 combination, tagged nuclei of a given MB KC subtype were purified from 400 fly heads. Nuclear RNA was then extracted and used to generate RNA-seq libraries ( Figure  S2B). We generated libraries from two independent biological replicates for each MB KC subtype, and paired-end sequenced them (Table  S1). We estimated transcript abundances in each library using kallisto (Bray et al. 2016;see Methods). The sequencing reads and estimated abundances are available in NCBI GEO (GSE119629). Transcript abundances were well correlated between replicate libraries ( Figure 2B). We observed strong expression of the neuronspecific marker elav in all subtypes (1,118 -1,784 TPM), contrasting with low levels of the glial-specific gene repo (0.1 -1.3 TPM; Figure 2C), consistent with a high-fidelity purification of TAPIN labeled nuclei. We also detected strong and broad expression of genes expected in all MB neuron subtypes, such as the transcription factor ey (919 -2,490 TPM) and components of the cAMP signaling pathway such as rutabaga (735 -1,629 TPM) which encodes the Ca 2+ /calmodulin-activated adenylyl cyclase ( Figure 2C; Crittenden et al. 1998). The TAPIN-seq profiles also recovered the expected pattern of genes known to be enriched in individual classes, including trio, Fas2, and sNPF.

Genes enriched in MB neuronal classes and subtypes
We next identified transcripts that are differentially expressed across MB neuron classes or subtypes using three criteria (Methods). We found 341 transcripts that are enriched or depleted in one of the three MB neuron classes ( Figure 3A) and 57 that are enriched or depleted in one of the seven MB cell subtypes ( Figure 3B; Table 2 for summary and  Table S2 for the full gene list). To evaluate the accuracy of our differential expression analysis, we examined several genes encoding proteins reported to differentially label MB cell classes. Antibodies against trio, for example, are reported to label a9/b9 and g classes (Awasaki et al. 2000). Indeed, the trio gene is identified in our analysis as an MB classspecific gene depleted in a/b KCs (average 179 TPM in a/b vs. 1,345 TPM a9/b9 and 674 TPM g), and we confirm that anti-trio immunoreactivity is localized in the MB a9/b9 and g lobes, and their cell bodies ( Figure 4A). Fas2 has been reported to exhibit strong immunoreactive signal in the MB a/b lobe class, weak signal in the g lobe class and no signal in the a9/b9 lobe class (Crittenden et al. 1998). Consistent with this protein distribution, we identify Fas2 transcripts as an MB classspecific gene depleted in a9/b9 KCs (124 TPM a9/b9 vs. 959 TPM a/b and 734 TPM g). At the cell type level, Fas2 was enriched in in the gd subtype relative to gm+d (1,052 vs. 417 TPM, respectively). Using immunolabeling, we confirmed this pattern of anti-Fas2 immunoreactivity. As previously reported, we detect strong immunoreactive signal of Fas2 in the a/b lobe class of neurons, somewhat weaker signal in the g lobe class and no signal in the a9/b9 lobe class neurons. We further demonstrate that Fas2 immunoreactivity appears weaker in the gm subtype compared to gd subtype KCs ( Figure 4B), which mirrors the prediction from the TAPIN-seq described above.
We also examined the expression of the neuropeptide gene sNPF which we identified as an MB class-specific gene depleted in a9/b9 KCs (20 TPM in a9/b9 vs. 414 TPM a/b and 725 TPM g). Immunolabeling with an antibody against the sNPF precursor confirmed a previous report and is consistent with the TAPIN-seq results. We observe no detectable signal in a9/b9 KC class neurons and strong signal in a/b and g classes ( Figure 4C; Johard et al. 2008). No noticeable signal was detected in the cell bodies of KCs ( Figure S3), consistent with a previous description (Johard et al. 2008). At the protein level, we further noted elevated immunoreactivity in the MB a/bp subtype, which was not reflected by the TAPIN-seq results ( Figure 4C). This discrepancy points either to limitations of TAPIN mediated profiling with these split-GAL4 lines and/or the importance of post-transcriptional regulatory mechanisms in determining the accumulation of the neuropeptide. Taken together, the identification of known differentially expressed transcripts and the immunostaining data for three identified examples broadly corroborate the fidelity of the TAPIN-seq results.
The differentially expressed genes also include several transcriptional regulators with class-and subtype-specific patterns (Figure 3). For example, we identify transcriptional regulators enriched in the gm+d (zfh2), a9/b9m (Ets96B, otp), a/bp (br, bru3, dl, disco-r, ind, rn), and a/bc (Tet) subtypes. Although the abundance and significance of DNA methylation in Drosophila is unclear, enrichment of the Tet DNA methyltransferase is intriguing given the role of DNA methylation in memory formation in other insects (Biergans et al. 2012) and mammals (Day and Sweatt 2010). The subtype-enriched expression could reflect a remnant of a functional expression pattern from an ancestral species with DNA methylation.
Several genes encoding cell-surface molecules were also differentially expressed across the subtypes. These genes included members of gene families previously implicated in specifying synaptic connectivity, including the defective proboscis extension response (Dpr) as well as the Dpr-interacting protein (DIP). Interactions between proteins from these families have been previously documented and shown to underlie synaptic connectivity in circuits including the visual system (Özkan et al. 2013;Carrillo et al. 2015). The expression patterns we observed for several cell-surface molecules suggest they might be involved in specifying class-and subtype-specific connectivity.

Neurotransmitter output
To explore the functional utility of our TAPIN-seq measurements, we next focused on genes related to the input and output properties of the MB cell subtypes. Specifically, we examined genes encoding neurotransmitter biosynthetic enzymes, neurotransmitter transporters ( Figure  5A), neurotransmitter receptors ( Figure 5B), neuropeptides ( Figure  5C), neuropeptide and protein hormone receptors ( Figure 5D), and gap junction components ( Figure 5E). A recent report has established that the primary neurotransmitter in the MB is acetylcholine (Barnstedt et al. 2016). Because our method profiled each of the seven MB neuronal subtypes, we re-assessed this conclusion at a higher resolution. We confirmed that all seven MB cell subtypes express high levels of both choline acetyltransferase (ChAT; 75 -198 TPM) and vesicular acetylcholine transporter (VAChT; 257 -517 TPM).
We next asked whether one or more KC subtypes might co-release other small molecule neurotransmitters, but found no strong argument to support such a conclusion. In a few cases, we see moderate expression of biosynthetic enzymes for other neurotransmitters including GABA, serotonin and dopamine. But in each case, there are findings that undermine the conclusion that these neurotransmitters are consistently produced/released by any of the 7 KC subtype. For example, we do see strong expression in all MB cell subtypes of vesicular GABA transporter (VGAT; 152 -279 TPM), which is an essential transporter that is responsible for packaging the neurotransmitter GABA into synaptic vesicles (Fei et al. 2010). On the other hand, the gene for GABA biosynthetic enzyme, Gad1, is moderately expressed in only the a/bp (121 TPM) and a/bc subtypes of KCs (54 TPM) and at lower levels in the remaining cell subtypes (4 -14 TPM). In principle, moderate Gad1 expression might be consistent with the hypothesis that GABA is released from a fraction of a/bp and/or a/bc subtype neurons. To examine this possibility, we conducted immunofluorescence experiments n with antibodies against GABA, Gad1, and VGAT but found no marked immunoreactivity in any a/bp KC subtype cell bodies ( Figure S4). This result suggests that the most likely explanation for our observed Gad1 pattern is that the split-GAL4 lines we used for a/bp and a/bc KC subtypes drive low levels of expression in some subset of GABAergic neurons outside the MB ( Figure S4A). A similar set of findings are apparent with Dopa decarboxylase (Ddc), a commonly used marker for dopaminergic or serotonergic neurons. Ddc catalyzes the decarboxylation of dopa to dopamine and 5-hydroxytryptophan to serotonin but not tyrosine to tyramine (Gramates et al. 2017). We see fairly strong Ddc expression in all MB cell types (55 -574 TPM), especially the gm+d (574 TPM) and gd (339 TPM) KCs. However, we do not see expression of Tryptophan hydroxylase (Trh; 0 -3 TPM), which provides the first and rate-limiting step in the synthesis of serotonin. Nor do we detect serotonin transporter (SerT; 0 -2 TPM; Giang et al. 2011) in any of the 7 cell subtypes ( Figure 5A). Pale (ple), which encodes tyrosine hydroxylase (TH) for dopamine synthesis, also is not highly expressed in any MB cell subtype (2 -16 TPM; Figure 5A). Thus, we conclude that all MB KC subtypes likely release acetylcholine (Barnstedt et al. 2016;Crocker et al. 2016) and no other small molecule neurotransmitters.

Neurotransmitter receptors
We next examined expression profiles of small molecule neurotransmitter receptors. Dopamine has been established as a key input to MB for many behaviors including aversive and appetitive olfactory learning (Kim et al. 2007;Claridge-Chang et al. 2009;Qin et al. 2012;Liu et al. 2012;Burke et al. 2012), modulation of motivational state (Krashes et al. 2009), regulated forgetting (Berry et al. 2012;Shuai et al. 2015), sleep (Sitaraman et al. 2015b), courtship behaviors (Kuo et al. 2015;Lim et al. 2018), and temperature preference (Bang et al. 2011). A network of dopaminergic neurons innervates all MB cell subtypes (Aso et al. 2014a), and all four dopamine receptors, Dop1R1, Dop1R2, Dop2R and DopEcR, are expressed in all 7 MB KC subtypes (594 -1152, Figure 2 TAPIN-seq profiling of MB subtypes. (A) Driver lines expressing in the seven Kenyon cell subtypes were crossed with the TAPIN-seq reporter, which labels the nuclei in each subtype. Nuclear RNA from each subtype was used to generate RNA-seq libraries, which were sequenced in paired-end mode. (B) We estimated reproducibility of the TAPIN-seq measurements by calculating the Pearson correlation between estimated transcript abundances (log2 transformed Transcripts Per Million + 1). (C) The TAPIN-seq transcriptomes recover the neuronal marker elav, while not detecting the glial marker repo. The transcriptomes also recover the expected expression patterns of the known pan-Kenyon cell markers ey and rut as well as class-enriched genes trio, Fas2, and sNPF.
Serotonergic signaling in the MB also is involved in olfactory memory formation, sleep regulation and stress response modulation (Yuan et al. 2006;Lee et al. 2011;Haynes et al. 2015;Ries et al. 2017). We found that among five serotonin receptors, only 5-HT1A is expressed in the MB, with an enrichment in the a/b lobe KC class ( Figure 5B). This is consistent with the finding that dorsal paired medial (DPM) neurons release serotonin onto 5-HT1A receptors expressed in a/b KCs to support anesthesia-resistant memory formation (Lee et al. 2011). A previous report observed 5-HT1B-GAL4 expression and 5-HT1B immunoreactivity in the g KC class (Yuan et al. 2005). We did not detect nuclear 5-HT1B transcripts by TAPIN-seq (0.1 -1.0 TPM), but we cannot rule out the possibility that higher levels of transcripts and protein are present in the cytoplasm, beyond detection in this nuclear transcriptome.
MB KCs receive nicotinic acetylcholine receptor (nAChR)-mediated synaptic transmission from antennal lobe projection neurons (Gu and O'Dowd 2006). We found the nAChR subunits nAChRalpha1, nACh-Ralpha4, nAChRalpha5, nAChRalpha6, nAChRalpha7, nAChRbeta1 and nAChRbeta2 are strongly expressed in all seven MB cell subtypes, but nAChRalpha3 and nAChRbeta3 transcripts are absent or undetectable ( Figure 5B). Two muscarinic acetylcholine receptors (mAChRs), mAChR-A and mAChR-B, are also expressed in the MB with an enrichment in a/b lobe KC class ( Figure 5B & Table S2).
We also examined the expression of the six known octopamine receptors (Gramates et al. 2017). We found that Oamb and CG18208 (recently characterized as Octa2R; Qi et al. 2017), which are a-adrenergic-like receptors, are class-and subtype-specific genes, respectively (Table S2). Oamb TAPIN-seq signal is strongly detected in a/b class KCs and far lower levels are seen in a9/b9 class KCs ( Figure 5B) (Crittenden et al. 1998;Kim et al. 2013). The three b-adrenergic-like receptors -Octb1R, Octb2R and Octb3Rare all detected in the MB with variable levels across cell subtypes (cf Wu et al. 2013 for Octb2R immunolabeling).
Although transient glutamate immunoreactivity has been shown in a/bc lobe KC class of young adult, VGlut expression has never been observed in the MB KCs (Daniels et al. 2008;Sinakevitch et al. 2010). Consistent with this conclusion, we observe low VGlut transcript abundance (5 -50 TPM; Figure 5A). We also confirmed that the NMDA receptors, both Nmdar1 and Nmdar2, are broadly expressed in the MB ( Figure 5B) (Xia et al. 2005;Ueno et al. 2017); ionotropic receptors GluRIA and GluRIB are also expressed ( Figure 5B). Flies also have a unique metabotropic glutamate receptor called mGluR, which has been  previously observed by immunolabeling throughout the adult brain, but minimally in the MB lobes (Devaud et al. 2008). Here, with the cell type resolution of our dataset, we identified mGluR as an MB class-specific gene that is expressed in a subset of KCs -mGluR is depleted in a/b class KCs (40 TPM; Table S2) relative to both a9/b9 and g classes of KCs (504 and 612 TPM, respectively; Figure 5B).

Neuropeptides
In addition to small neurotransmitter systems, we also examined expression of neuropeptides and neuropeptide receptors. We consistently detect expression of a small group of well-characterized and putative neuropeptides. Unexpectedly, this includes amnesiac (amn), which is detected at fairly robust levels in all 7 MB KC subtypes (52 -175 TPM)although this small gene resides in the intron of another gene, Hers, which is highly expressed (218 -239 TPM), thus complicating the accurate estimation of amn abundance. Previous work has established a requirement during memory formation for amn expression in a single pair of DPM neurons outside MB (Waddell et al. 2000). sNPF expression appears as a class-specific transcript, with high levels in a/b and g KC classes and very little expression in the a9/b9 KC class (414, 725, 20 TPM, respectively; Figure 4C). This is consistent with previous reports (Johard et al. 2008) and we further confirm this expression pattern by immunostaining ( Figure 4C). Another notable class-specific neuropeptide is Drosophila Insulin-like peptide 1 (Ilp1; Liu et al. 2016b), expressed at high levels in a9/b9 class of KCs (284, 382 TPM in a9/b9, a9/b9m, respectively) and a/bc subtype (220 TPM) but not in most other MB KC subtypes (0.8 -55 TPM; Figure 5C). As with neuropeptides, we also detect a panel of neuropeptide and protein hormone receptors, some of which are robustly expressed in all MB KC subtypes (e.g., Ecdysone receptor), some of which are class specific (e.g., Dh44-R1 and hector), and some of which are enriched or depleted in one or more subtypes of neurons (e.g., AstC-R1 and ETHR; Figure 5D). Finally, because of a report that gap junctions may form between MB KCs of different classes and play a role in visual learning (Liu et al. 2016a), we examined the expression levels of the eight gap junction genes in MB cell subtypes. shakB is strongly expressed (433 -751 TPM), and Inx3 moderately (29 -107 TPM), in all MB cell subtypes ( Figure 5E). Although Inx5 and Inx6 are reportedly required in the a/b and a9/b9 KC classes for visual learning and memory (Liu et al. 2016a), we do not detect either gene in any of the 7 MB KC subtypes (0 -0.2 TPM, 0 -0.5 TPM, respectively; Figure 5E). We cannot rule out the possibility that functionally relevant levels of expression are below our detection limit or that cytoplasmic RNA levels are higher.

DISCUSSION
Our results establish a high-quality, neuronal cell type-level transcriptome for Drosophila MB. We identified 350 differentially expressed genes, which includes most of the previously reported MB lobe (class specific) markers and many novel class-specific or cell subtype-specific profiles of expression. In addition to the subtype level resolution of our experimental design, the TAPIN approach that we used also offers several advantages and technical differences with these prior approaches. First, because TAPIN is compatible with flash frozen tissue as the input, the method introduces minimal disturbance to the endogenous transcriptome as compared to more lengthy procedures for purification of neurons for expression profiling. Second, it may be relevant that TAPIN explicitly profiles nuclear RNAs, likely enriching for actively transcribed/nascent transcripts vs. abundant ones that are stably maintained in the cytoplasm. Thus, it would be attractive to apply this method to profile transcriptional response to behavioral perturbations.
Several previous studies have used genome-wide methods to profile expression in the Drosophila MB (Perrat et al. 2013;Crocker et al. 2016;Croset et al. 2018;Davie et al. 2018;Jones et al. 2018). Perrat et al. used a microarray-based approach to profile expression of each of the three major classes of MB KCs (purified by flow cytometry from dissociated brains) and compared these profiles with expression in the rest of the brain. They first focused on the expression of transposons (Perrat et al. 2013), and subsequently used the same transcriptome dataset to discover that MB KCs are cholinergic (Barnstedt et al. 2016) based on expression of biosynthetic enzymes. Crocker et al. used an RNA-seqbased approach to profile expression in relatively small pools of physically isolated a/b and g class neurons to search for memory-related changes in gene expression (Crocker et al. 2016). Most recently, Croset et al. and Davies et al. used droplet-based single cell sequencing to profile the Drosophila brain, and by clustering the single cells they were able to identify the three MB classes, but not the further sub-division into neuronal subtypes (Croset et al. 2018;Davie et al. 2018).
Although we used a different profiling method and resolved transcriptomes at the cell subtype rather than class level, our findings are broadly compatible with prior reports (Perrat et al. 2013;Barnstedt et al. 2016;Crocker et al. 2016;Croset et al. 2018;Davie et al. 2018). Our dataset reveals strong expression of both ChAT and VAChT, consistent with the conclusion that MBs are cholinergic (Barnstedt et al. 2016). Our findings further support the conclusion that all of the individual MB KC subtypes are cholinergic. The datasets also are consistent in the expression of known class-specific markers. One notable difference is that Crocker et al. reported high levels of expression of the 5-HT1B receptor in both a/b and g classes of KCs (Crocker et al. 2016), and Davie et al. also observed 5-HT1B expression in a/b and g KCs single cell clusters (Davie et al. 2018). In contrast, we see no evidence for expression of this receptor in our TAPIN-seq profiles (0.1 -1 TPM). This difference could reflect methodology: Crocker and Davie both measured 5-HT1B receptor transcripts in the cytoplasm while we measured the levels that are actively transcribed or present in the nucleus. This technical difference could be especially relevant for neurotransmitter receptors, some of which can be translated locally at dendrites (Steward and Banker 1992).
Our dataset is the first to profile expression in this brain region at neuronal subtype resolution. This level of resolution is critical given the wealth of data on the functional differences of each MB KC subtype in Drosophila behaviors. Our dataset provides a full accounting within each of the MB KC subtypes of the profiles of expression of the cellular machinery to produce and receive neurotransmission, including small molecule transmitters and their receptors, neuropeptides and neuropeptide receptors and subunits of gap junctions. It is noteworthy that the TAPIN expression dataset supports the conclusions that all the adult MB KC subtypes are cholinergic, and that none of the subtypes express genes that would suggest the co-release of GABA, dopamine, glutamate, or serotonin. On the other hand, we detect expression of a spectrum of neuropeptides and their receptors. This observation is consistent with the hypothesis that MB KCs may co-release both acetylcholine and several neuropeptides (Takemura et al. 2017).
In addition to these findings with regards to the inputs and outputs, we identified 350 differentially expressed genes including many that distinguish MB KC classes or even individual cell subtypes. MB a/bp subtype showed 21 enriched genes and 11 depleted genes, contrasting with two other subtypes in the a/b class and two other classes. This uniqueness is supported by its unique odor responses (Perisse et al. 2013) and connectivity (Chen et al. 2012). Despite the limitation in the methodology that we profile the MB gm subtype using the spilt-GAL4 line MB131B that has minor expression in gd, and profile the a9/b9ap subtype using MB370B that has minor expression in a9/b9m ( Figure 1 & Table 1), we still identified distinct sets of enriched/depleted genes ( Figure 3 & Table S2), indicating the differences between two subtypes in MB g or a9/b9 classes.
Our dataset provides a valuable resource for the fly neuroscience community to conduct functional studies. For example, our data provide a list of previously unknown class specific and sub-type specific transcripts, whose impact on the functional differences between these neurons are not known (Figure 3; Table S2). An arsenal of genetic tools to manipulate any gene's function within each of these cell subtypes already exists (Caygill and Brand 2016;Kaya-Çopur and Schnorrer 2016). In addition to olfactory associative memory, MBs also play fundamental roles in other forms of memory including visual and gustatory (Vogt et al. 2014;Masek and Keene 2016), temperature preference (Hong et al. 2008), courtship behaviors (Kuo et al. 2015;Lim et al. 2018), stress response (Ries et al. 2017), food-seeking (Tsao et al. 2018), sleep (Artiushin and Sehgal 2017) and responses to ethanol (Kaun et al. 2011). This dataset will facilitate the discovery of neural mechanisms for each of these conserved behaviors.