Heterogeneity of hypothalamic pro-opiomelanocortin-expressing neurons revealed by single-cell RNA sequencing

Objective Arcuate proopiomelanocortin (POMC) neurons are critical nodes in the control of body weight. Often characterized simply as direct targets for leptin, recent data suggest a more complex architecture. Methods Using single cell RNA sequencing, we have generated an atlas of gene expression in murine POMC neurons. Results Of 163 neurons, 118 expressed high levels of Pomc with little/no Agrp expression and were considered “canonical” POMC neurons (P+). The other 45/163 expressed low levels of Pomc and high levels of Agrp (A+P+). Unbiased clustering analysis of P+ neurons revealed four different classes, each with distinct cell surface receptor gene expression profiles. Further, only 12% (14/118) of P+ neurons expressed the leptin receptor (Lepr) compared with 58% (26/45) of A+P+ neurons. In contrast, the insulin receptor (Insr) was expressed at similar frequency on P+ and A+P+ neurons (64% and 55%, respectively). Conclusion These data reveal arcuate POMC neurons to be a highly heterogeneous population. Accession Numbers: GSE92707.


INTRODUCTION
The leptin-melanocortin pathway plays a key role in the control of food intake and body weight, with genetic disruption of multiple components of the pathway resulting in severe, early-onset obesity in both humans and mice (reviewed by [1]). In the brain, POMC expression is localized to neurons in the arcuate nucleus (ARC) of the hypothalamus and in the nucleus of the solitary tract (NTS) [2]. POMC is a pro-peptide that is processed to a range of bioactive peptide products including a and b-MSH [3,4], both of which act on melanocortin receptors 3 and 4, which are widely expressed in the central nervous system. The role of MC4R in mediating the effects of POMC on energy balance is well established, whereas the role of MC3R is less clear [2]. In the ARC, anorexigenic POMC neurons exist alongside orexigenic neurons that express the endogenous melanocortin antagonist agouti-related peptide (AgRP) [5,6]. The canonical view of the melanocortin pathway is one where leptin, which is produced from white adipose tissue to reflect fat mass, signals to leptin receptors on the POMC and AgRP neurons, increasing the activity of POMC neurons and reducing the activity of AgRP expressing neurons [1,2]. However, there is a body of evidence that indicates this view of the leptin melanocortin pathway to be an over-simplistic representation of the relevant circuitry. Woods and Stock, using FOS immunoreactivity of the immediate early gene c-fos as a marker of neuronal activation, reported that while neurons in the arcuate nucleus of leptin deficient ob/ob mice were activated by peripheral leptin administration, the area of the hypothalamus most strongly activated was the paraventricular nucleus [7]. Later, it became clear that not all POMC neurons express the receptors for leptin or insulin, with whole-cell patch-clamp electrophysiology indicating that a subpopulation is depolarized by leptin, while a separate and distinct population is hyperpolarized by insulin [8]. POMC neurons also vary in their expression of the neurotransmitters glutamate, Ƴ-aminobutyric acid (GABA) [9,10], and acetylcholine [11] and their downstream projections. For instance, 40% of POMC neurons appear to be GABAergic and hence inhibitory [9]. In addition, only about 40% of POMC neurons express the serotonin 2C receptor (5HT 2C -R), the target for the weight-loss drug Lorcaserin, through which it is thought to mediate much of its actions on satiety [12,13]. Finally, mice with a targeted deletion of the leptin receptor specifically from POMC neurons display a surprisingly mild phenotype [14], as compared to complete Pomc null mice [15,16] and db/db mice homozygous for a Lepr mutation [17]. A recent paper by Henry et al. reported the gene expression profiles of pooled AGRP and POMC neurons from fed and food deprived mice and found that the global pattern of gene expression in response to an acute energy deficit are far more dramatic in AGRP expressing than POMC expressing neurons [18]. However, because of its use of pooled cells, the transcriptomic profiles of individual neurons, and hence their variability in gene expression, was not examined. We now report the results of single cell RNA sequencing in isolated POMC expressing neurons of the arcuate nucleus, which has revealed a previously unappreciated degree of heterogeneity within this anatomically highly localised population.

RNA-sequencing of POMC neurons
Using a POMC-eGFP mouse [19], we successfully sequenced the transcriptome of 163 GFP positive (MBH Pos) and 7 GFP negative cells (MBH Neg) from hypothalamic tissue encompassing the arcuate nucleus ( Figure 1A), and 24 cells, as positive controls, from the pituitary (Pit) known to express high levels of POMC (GEO Acc# GSE92707). When we visualized the transcriptomic data using t-SNE dimension reduction analysis [20], the pituitary and hypothalamic cells segregated into two discrete clusters, with the hypothalamic neurons further segregating into GFP positive and GFP negative cells ( Figure 1B). All cells sequenced were of similar size, as measured by degree of forward scatter on the cell-sorter and expressed both Actb and Gapdh ( Figure 1C). All the hypothalamic cells expressed the neuronal marker Ncam, but, because the pituitary cells contained both neuronal derived posterior and the ectodermally derived anterior, there was a larger heterogeneity of Ncam expression in the pituitary cells ( Figure 1C). In contrast, the GFP signal, which reflects the degree of Pomc expression, was consistently high with very little heterogeneity in the pituitary POMC cells, whereas hypothalamic GFP neurons show a much larger variation in signal. Cartpt, which is canonically expressed together with Pomc in the same neurons, also showed variable expression in hypothalamic GFP neurons and was not detected in the pituitary ( Figure 1C). POMC is processed by the proprotein convertases, PC1 (PCSK1) and PC2 (PCSK2). Using the detection cut-off !5 linear RPM, most cells (156/163) expressed either or both Pcsk1 and Pcsk2 ( Figure 1D). On average, 6882 different transcripts were detected from each individual neuron (!5 linear reads per million or RPM). However, across all 163 POMC neurons, 12,187 different transcripts were expressed ( Figure 1E), indicating significant heterogeneity in gene expression repertoire among the neurons.

AgRP/NPY expressing POMC neurons
Unexpectedly, we found 27% (45/163) of Pomc-expressing neurons also expressed high levels of the endogenous melanocortin antagonist Agrp (Figure 2A, C). In addition, these neurons expressed a high level of Npy, while Pomc expression itself was a magnitude (>5 log 2 FC) lower than that seen in POMC neurons not expressing or expressing low levels of Agrp ( Figure 2B). Comparing our data with a published transcriptome of pooled AgRP neurons [18], we found a significant similarity in gene expression profiles (R 2 ¼ 0.71 p < 0.00001), with high expression, for instance, of the ghrelin receptor Ghsr ( Figure 2B). Analysis of the raw data from Henry et al., 2015 also indicated that Pomc and Cartpt were expressed in the AgRP neurons that they had analysed [18], albeit at low levels ( Figure 2D). Taken together, these 45 high Agrp and low Pomc expressing neurons (hereafter referred to as A þ P þ neurons) showed a striking similarity to canonical AgRP/NPY neurons, likely reflecting their developmental origins, with up to 25% of AgRP neurons sharing the same hypothalamic progenitors as POMC neurons [21]. We defined the rest of the 118 neurons, which expressed higher levels of Pomc, as canonical POMC expressing neurons (P þ ) and conducted our further studies on this population.

Expression of hypothalamic neuroendocrine genes in POMC neurons
We produced a matrix showing the expression distribution of 23 well characterized neuroendocrine genes, all known to play a role in the regulation of energy homoeostasis, on the individual POMC (P þ ) neurons ( Figure 3A). We found, as has been previously reported, that a large percentage of POMC neurons expressed Mc3r, Cnr1, Htr2c, and the Npy-y1, y2, and y5 receptors [1]. Of note, the percentage of Htr2c expressing neurons (34/118; 29%) was comparable to previous studies showing that 5HT 2C Rs (the protein product of Htr2c) are expressed on approximately 40% of POMC neurons [12]. We also found that a smaller percentage of neurons, 12/118, expressed Mc4r (co-localization validated as shown in Figure 3B), with 11/118 expressing Glp1r. ARC POMC neurons expressing GLP1R are thought to mediate, in part, the weight-loss effects of the GLP1 analogue liraglutide [22]. While we excluded the 45 A þ P þ neurons from this matrix, a large number of the remaining P þ neurons still expressed low levels of Agrp (66/118) and Npy (98/118). Our own data ( Figure 3A) and re-analysis of the raw data from previously published work indicated that the converse was also likely to be true, with a percentage 'bona fide' AgRP neurons also expressing low levels of Pomc ( Figure 2D).

Unbiased clustering reveals four 'classes' of POMC neurons
In order to determine whether P þ neurons segregate into distinct classes, we first adjusted for batch effects and then performed unbiased clustering of the 118 cells based on the consensus matrix of Euclidean, Pearson, and Spearman correlations e using the 'SC3' R package [23]. We obtained 4 main clusters ( Figure 4A), which could also be visualized using t-SNE [20] ( Figure 4B). Cluster 1 was the largest with 69 cells, followed by cluster 4 (27) and cluster 2 (15). Cluster 3 had the least number of cells (7). Cluster 2 showed resemblance to cluster 1 ( Figure 4A, B). We identified 2773 candidate genes that were differentially expressed among the 4 clusters, with a false discovery rate (FDR) cut-off of 5%, that drove the segregation of cells into the different groups ( Figure 4C). Interestingly, there was over-representation of membrane (530) and secretory genes (191) within the driver gene list ( Figure 4C, D). These membrane genes encode ion channels, G-protein coupled receptors (GPCRs) and other transmembrane receptors, and transporters ( Figure 4D). Figure 5A lists the top 25 genes, ranked by differential expression that drove the clustering of the four groups (with the exception of group 2, which was a small group and only has 3 significantly differentially expressed genes). Shown are the expression level for each of the genes and how unique each gene was to each group. Additionally, Table 1 contained genes that are not in the top 25, encoding selected extracellular secreted hormones, GPCRs and ion channels, indicating which groups had the highest expression of any given gene. For example, while the neurons in cluster 4 expressed the highest levels of Pomc and Cartpt, both genes, as expected, were also expressed in Clusters 1e3. Cluster 4 neurons, however, were characterized by expression of genes encoding transthyretin (Ttr), the cholesystekinin A receptor (Cckar), neuromedin U receptor 2 (Nmur2), and nascent helix loop helix 2 (Nhlh2) ( Figure 5A), as well as the GPCRs corticosterone releasing hormone receptor 1 (Crhr1), hypocretin receptor 1 (Hcrtr1), neuropeptide Y 2-receptor (Npy2r), and prolactin releasing hormone receptor (Prlhr) ( Table 1) at far higher levels than in the other clusters. Cluster 1, in contrast, was characterized by high expression of the serotonin-2a receptor (Htr2a), the melanocortin-3 receptor (Mc3r) and a large number of ion channels (Table 1). Thus, no single given 'driver' gene is responsible for these clusters, but rather the combined expression of a repertoire of genes drove the segregation. The entire repertoires of genes for each cluster are listed in Supplementary Table 1. These clusters are, of course, not fixed. For example, Figure 5B shows the t-SNE plot with the cluster 4 neurons highlighted and dividing clearly into two groups; one larger group of 22 and a smaller one of 5 neurons. Comparing the transcriptomes of the two groups, we could identify differentially expressed genes that defined the two groups ( Figure 5C, Supplementary Figure). The four clusters, as we report here, were achieved by using a specific statistical threshold. If we raised this threshold, and therefore statistical rigour, then all clusters, could be further segregated and, if taken to its extreme, could segregate into the original 163 individual neurons sequenced.

Glutamatergic and GABAergic POMC neurons
Previous studies have shown that POMC neurons can also be classified by the release of stimulatory neurotransmitter glutamate or inhibitory gamma-aminobutyric acid (GABA) [9,10,24]. We examined the two subtypes by looking at expression of vesicular glutamate transporter type 2 (Vglut2/Slc17a6) as a proxy for glutamatergic POMC neurons and glutamate decarboxylase (Gad67/Gad1) for their GABAergic counterparts. In Figure 5D, we found that the majority of cells (87%) express high levels of Gad1 (indicated in blue), and these included most of the A þ P þ neurons as well as 70% of the P þ neurons. On the other hand, only about half of the P þ neurons have high levels of Slc17a6 (indicated in red), and these tended to fall into clusters 1 and 4 ( Figures 4B and 5D). 24% of P þ neurons express both markers (indicated in purple) indicating the possibility of dual roles at different synaptic junctions.
2.6. Leptin receptor and insulin receptor expressing POMC neurons Next we examined the two most well studied types of POMC neurons, those expressing the leptin (Lepr) and insulin receptors (Insr). Looking at all 163 neurons, our sequencing data showed that 40/163 (25%) neurons expressed the leptin receptor and 94/163 (57%) expressed the insulin receptor. 11 (7%) expressed only Lepr and not Insr, 65 (40%) expressed only Insr and not Lepr, while 29 (18%) expressed both receptors. 58 (36%) did not express either of the receptors ( Figure 6A). When we split the 163 cells into the high Agrp expressing A þ P þ neurons and high Pomc expressing P þ neurons and looked at the distribution of Lepr in these two groups separately, we saw that 26/45 (58%) of A þ P þ neurons expressed Lepr as compared to only 14/118 (12%) of P þ neurons, almost five-fold more. In contrast, the Insr is expressed in 29/45 (64%) of A þ P þ neurons, which is comparable to its expression in 65/118 (55%) of P þ neurons. To determine if we had introduced any selection bias during the flow sorting, we sought to define the number of Lepr and Insr expressing POMC neurons in the intact mouse hypothalamus. Attempts to use immunocytochemistry to detect either the leptin or the insulin receptor on whole brain slices proved to be technically unreliable (data not shown). We therefore measured levels of Stat3 phosphorylation (pStat3) after leptin administration and levels of Akt phosphorylation (pAkt) after insulin administration using immunohistochemistry, as a proxy for expression of the respective receptors. We performed these experiments on POMC-GFP mice, allowing us to co-localize pStat3 or pAkt immunofluorescence in neurons expressing GFP ( Figure 6B, C). After leptin administration, pStat3 positive cells were largely concentrated to the base of the arcuate nucleus, surrounding the median eminence ( Figure 6B). Cell counting across multiple brain slices from three mice indicated that 23.8% of POMC-GFP neurons were positive for pStat3 in response to leptin, which was consistent with our RNAseq data ( Figure 6C). While the leptin responsiveness in POMC neurons was previously shown to be variable, dependent on the treatment protocols [25], the percentage we show here falls within the lower limit of what others have previously reported [26]. Insulin administration, however, resulted in a different pattern of activation, with pAkt  (D) All of the A þ P þ neurons express GABAergic marker Gad1 (Gad67, blue) and the same applies to 70% of the P þ neurons. We also found about 50% of the P þ neurons that express glutamatergic marker Slc17a6 (Vglut2, red), 24% of the P þ neurons express both markers indicating dual roles at synaptic junctions.   immunofluorescence evident more broadly through the arcuate nucleus ( Figure 6B). We found that 41.1% of POMC-GFP neurons were positive for pAkt in response to insulin, which was not significantly different (p ¼ 0.357) to the percentage of GFP neurons expressing Insr as determined by our sequencing data ( Figure 2C). These findings increased confidence that the FACS sorting has captured a representative snapshot of the POMC neuronal population.

DISCUSSION
We describe here the use of single cell RNA sequencing to determine the transcriptomes of 163 individual murine POMC neurons. Unexpectedly, 25% of POMC positive neurons also express high levels of Agrp. The transcriptomes of these neurons are comparable to those of previously published pooled AgRP neurons [18], indicating that these are not an artifact but likely to be a real subpopulation of canonical AgRP neurons. Additionally, when we examine the raw data from Henry et al., we find that AgRP/NPY GFP neurons also express low levels of POMC [18]. This phenomenon likely reflects the shared developmental origins of POMC and AgRP cells [21]. POMC has been shown to be broadly expressed in hypothalamic progenitors, with up to 25% of the mature AgRP/NPY population sharing a common progenitor with POMC neurons [21]. The A þ P þ population we describe here could be cells where POMC transcription is never fully silenced, even after the differentiation into AgRP/NPY neurons. Another possibility is that the high Agrp/Npy expressing neurons could be newly committed AgRP/NPY neurons from these progenitor cells, in which Pomc expression is yet to be fully suppressed. Regardless of the explanation, given that Pomc and AgRP genetic tools are in wide use, it is important for the field to know that a 'Pomc-specific' deletion of a given gene would also be removing the gene from a significant proportion of AgRP neurons. It is possible, for example, that the mild phenotype seen in the mice in which Lepr is 'specifically' deleted from POMC neurons [14] could be due in part to Lepr deletion from a proportion of AgRP neurons. The remaining 118 P þ neurons display striking heterogeneity and fall into four clusters. A few important points should be kept in mind while interpreting these data. First, these groups have emerged in an unbiased fashion after achieving a certain statistical threshold, which means that they are not invariant. Thus, increasing the threshold would invariably result in more clusters. Second, these groupings are not solely characterised by small numbers of unique genes; rather they are driven by a 'signature' expression profile of a whole repertoire of genes. Thus, while the driver genes (see Figure 5) are, on average, the highest expressed in their particular group, many transcripts are not unique to a single group and are not necessarily expressed in all of the neurons in any given group. Third, these different clusters are likely to be dynamic, depending on the nutritional state of the animal at the time the brain was obtained, with the mice in this particular experiment being ad libitum fed. The group with the highest expression of Pomc and Cartpt is cluster 4, and based on its expression profile of GPCRs, appears to be responsive to corticotrophin releasing hormone, hypocretin, prolactinreleasing hormone, Npy, CCK, and neuromedin U. The largest group is cluster 1, which, based on its repertoire of GPCR expression, responds to vasopressin, melanocortin peptides (via the MC3R), serotonin, and somatostatin. Cluster 1 is notable in that, compared to the other 3 groups, it expresses the highest levels of a wide variety of plasma membrane ion channels. We acknowledge that further experiments will be required to functionally characterize these neurons. However, we speculate that the enrichment of distinct classes of surface receptors on the different neuronal clusters possibly speaks to the temporal nature of the downstream effector function, with signalling through ion channels typically bringing about rapid changes, and circulating humoural factors providing longer term signals. Two drugs currently on the market that are used for weight-loss are the serotonin receptor agonist Lorcaserin [27] and the GLP1 analogue Liraglutide [28]. There is evidence that both mediate their effects on food intake via ARC POMC neurons expressing 5HT 2C R [12,13] and GLP1R [22], respectively. Our data provide two additional pieces of information ( Figure 3). First, less than 10% (11/118) of ARC POMC neurons express the GLP1R. Second, of the POMC neurons expressing either 5HT 2C R or GLP1R, only 3 neurons express both receptors. Thus, it is likely that different subgroups of POMC neurons mediate the satiety effects of each compound, which poses the question of whether a potential combinatorial therapy, without increasing individual drug-dose, might possibly result in increased weight-loss. As previously mentioned, the 45 A þ P þ neurons are very distinct from the 118 P þ neurons, clearly segregating into a separate group. This is particularly true when one examines the expression of the leptin and insulin receptors in these groups. The low percentage (12%) of Leprexpressing POMC (P þ ) neurons is consistent with previous evidence that the actions of leptin on these neurons are largely indirect [19,29]. Although we did not specifically or comprehensively evaluate an Agrp positive population of neurons, the much higher percentage of LepR positivity in the cells where we documented high levels of Agrp expression (57%) support the idea that AgRP neurons are more likely to be direct targets of leptin [19]. Our data also provide some evidence of autocrine/paracrine regulation by POMC neurons themselves, with 10% of P þ neurons expressing Mc4r and nearly 40% expressing the Mc3r. One could speculate that there may be amplification of the melanocortin system, whereby P þ neurons that respond directly to leptin go on to signal to Mc3r/Mc4r expressing neurons, but this remains to be determined. In contrast, the expression of insulin receptors is comparable between the A þ P þ and P þ neurons, indicating that insulins' actions on both these populations of neurons are likely to be direct [26,30].
In conclusion, single cell RNAseq has revealed previously unappreciated heterogeneity in POMC neurons. We have provided a comprehensive atlas, which we hope will be of utility to basic and translational researchers attempting to better understand the fundamental nature of regulation of energy balance by the hypothalamus and to manipulate these systems for therapeutic benefit.

Animals
7e9-week-old male Pomc-eGFP transgenic mice [19] (kind gift from Prof. Malcolm Low) were housed under a standard 12 h lightedark cycle and were fed ad libitum. All studies were approved by the local Ethics Committee and conducted according to the UK Home Office Animals (Scientific Procedures) Act 1986. 4.2. Isolation of single hypothalamic POMC-eGFP neurons using fluorescence-activated cell sorting (FACS) Ad libitum fed Pomc-eGFP mice were sacrificed by cervical dislocation, and tissue from the medial basal hypothalamus (MBH) region was micro-dissected into Hibernate A medium without calcium and magnesium (BrainBits, Springfield, IL, USA). The tissue was then digested with Papain (20 U/ml, Worthington, Lakewood, NJ, USA) for 30 min at 37 C with mild agitation, followed by trituration in Hibernate A medium containing DNase I (0.005%, Worthington). The cell suspension was then filtered through a 40 mm cell strainer into a new Falcon tube.
Fluorescence-activating cell sorting was performed using an Influx Cell Sorter (BD Biosciences, San Jose, CA, USA) utilizing a protocol that was previously used for single-cell isolation [31,32]. The cell gating was set according to cell size (FSC), cell granularity (SSC), FSC pulsewidth for singlets, fluorescence at 488 nm/532 nm for GFP and 647/ 670 nm for nuclear stain with DraQ5 (Biostatus, Shepshed, Leicester, UK). Cells were sorted direct into a 96 well plate containing lysis buffer with RNase inhibitor. 4.3. Whole-transcriptome amplification and single-cell RNA sequencing RNA from the single cell lysate was reverse transcribed and PCR amplified (19 cycles) into double stranded cDNA using Clontech SMART-Seq v4 Ultra Low Input RNA Kit (Takara Clontech, Mountain View, CA, USA). 150 pg of amplified cDNA was then used to generate barcoded sequencing libraries using Illumina Nextera XT library preparation kit (San Diego, CA, USA) according to manufacturer's manual. Samples failing to show amplification were removed from the experiment. The sequencing libraries were normalized to 10 nM concentration and combined into pools of 96. The pooled libraries were sequenced on 4 lanes of an Illumina HiSeq 4000 instrument at single-end 50bp (SE50), yielding an average of 15.7 million reads per cell. Library preparation was performed by the Genomics and Transcriptomic Core at the Institute of Metabolic Science. The sequencing was performed at the Genomics Core, Cancer Research UK Cambridge Institute. All the raw data in this manuscript has been submitted to GEO, accession number GSE92707.

Leptin and insulin administrations
For pSTAT3 immuno-labelling, mice were fasted overnight and leptin (3 mg/kg, PeproTech, Rocky Hill, NJ, USA) or saline was injected intraperitoneally. 45 min after leptin administration, the animals were sacrificed by a lethal overdose of pentobarbital (60 mg/kg) and perfused transcardially with PBS followed by 4% paraformaldehyde. For pAKT immunofluorescence, overnight fasted animals were administered Actrapid (2 mU/kg, Novo Nordisk, Bagsvaerd, Denmark) intravenously. Animals were sacrificed 15 min after administration and perfused as described before. 4.5. Immunohistochemistry for POMC-eGFP, pStat3 and pAkt Brains were collected immediately after perfusions and fixed for two additional hours in a 15% sucrose-4% paraformaldehyde solution. The brains were cryo-preserved in 30% sucrose PBS overnight, embedded in OCT (Bright, Huntingdon, UK), and stored at À80 C until further processing. Frozen coronal sections were sliced at 30 mm and processed for immunofluorescence. For pSTAT3 immunostaining, sections were pre-treated for 20 min in 1% NaOH and then placed overnight with a rabbit anti-pSTAT3 antibody (1:200, Tyr-705 Cell Signalling). For pAKT immunostaining, sections were incubated with anti-pAKT antibody (1:100, Ser-473 Cell Signalling) overnight. The primary antibodies were localized with Dylight 594 Goat anti-Rabbit IgGs (Vector; 1:500). Coronal images of the arcuate nucleus were taken using a Laser-scanning confocal microscope (LSM 510, Zeiss, Oberkochen, Germany), and eGFP positive cells and double labelled GFP-pSTAT3-or GFP-pAKT immune-reactive cells were counted using Image J. 4.6. Immunohistochemistry for Mc4r-MAPT, and Pomc B6.Cg-Tg(Mc4r-MAPT/Sapphire)21Rck/J male mice were perfused transcardially via a 23-gauge needle placed in the left ventricle with 100 ml of 0.1 M heparinized PBS, pH 7.4, followed by 100 ml of 4% paraformaldehyde in PBS, and the fixed brains were cryoprotected in 30% sucrose for 48 h at 4C. Coronal hypothalamic sections of 30 mm thickness were prepared on a freezing microtome. Free floating sections were incubated for 15 min in 1% hydrogen peroxide, washed 2 times in PBS, blocked 2 h in 0.3% Triton X-100 and 5% NGS in PBS, before being incubated overnight at 4 C in rabbit antiproopiomelanocortin precursor (1:200, Phoenix Pharmaceuticals) in 0.3% Triton X-100 and 5% NGS. Sections were then washed and incubated for 2 h with Alexa Fluor 594 goat anti-rabbit IgG. Sections were washed and incubated overnight at 4 C in chicken anti-GFP (1:500, Abcam). Sections were washed and incubated in Alex Fluor 488 goat anti-chicken IgG for 2 h. Sections were floated onto superfrost Plus microscope slides (Fisher), and coverslipped with Vectashield (Vector). Tissue sections were digitized using a Zeiss Axioplan fluorescent microscope, and areas of interest were outlined based on cellular morphology.

Bioinformatics
The quality of sequence reads was examined using FastQC, and Multiple Genome alignment (MGA) was used to rule out contamination from other DNA sources throughout the wet experimental procedures. Reads were aligned to the mouse genome (GRCm38, Ensembl v74) using Tophat version 2.0.11 [33], and the raw counts at the gene level were determined using Cufflinks version 2.2.1 [34]. The above analyses were performed at the High-performance computing cluster at the Research Computing Service, University of Cambridge. The raw counts from cells were normalized and batch correction was performed using edgeR [35] and limma [36] packages. Rtsne package was used to generate t-SNE plots [20] and SC3 package was used to perform the clustering analysis [23]. GeneSpring GX software (Agilent) was used to generate the heatmaps. The normalized gene abundance is represented in log 2 RPM, which stands for the base-2 logarithm of Reads per Million. The plots were generated using ggplot2 package, Microsoft Excel and GeneSpring GX. Gene annotations and pathway analysis were perform using Ingenuity Pathway Analysis (Qiagen, Redwood City, CA, USA).