Introduction

A key milestone for the study of neural circuits is to delineate the composition and connectivity of each brain area that controls specific behaviors. Systematic identification of the individual components of the circuits i.e., neuronal cell types, has seen a dramatic breakthrough in recent years¹. This includes not only the transcriptomic and proteomic compositions of individual cells but also their morphology, biophysical properties, and connectivity with other brain areas. Understanding these basic parameters is essential to map information flow and to gain insights into potential computations within circuits². In addition, each circuit component should be characterized in the context of the patterns of activity during different behaviors.

The multimodal characterization of cell types, however, remains technically challenging. While in situ transcriptomic methods allow cell type assignment in the native tissue context and are amenable to combining with other neural readouts, existing methods including commercial platforms, suffer from two major limitations for wide adoption. First, in situ transcriptomics experiments are resource intensive. While they provide vital spatial information for transcriptomic cell types, they have not been used as routinely as one might need for multimodal experiments. Second, most platforms use decoding of barcoded signals, requiring a complex analysis pipeline. One way to address these issues is to limit in situ transcriptomics experiments to specific brain areas and use a small number of marker genes to perform the profiling. There are several existing platforms^3,4 one can use to perform these experiments efficiently for up to 24 genes with sequential gene readout. Higher number of genes, however, remains more challenging due to complex requirements on probe designs and extended experimental timelines. To address these points, we have developed hamFISH (highly amplified multiplexed fluorescence in situ hybridization), a method that uses branched DNA amplification on a well-established smFISH platform with gene-specific bridge readout probes to sequentially detect gene expression. Additionally, the use of an optimized tissue preparation protocol to detect mRNA transcripts in brain tissue sensitively and flexibly makes hamFISH more amenable for routine use. We then validated the robustness of hamFISH in the medial amygdala (MeA), a critical hub for processing social information that receives prominent pheromonal input^5,6. Transcriptomic cell types in the MeA were recently defined by single-cell sequencing^1,7,8. We therefore decided to use these datasets as a starting point to map these cells in a multimodal fashion, i.e., their spatial localization, projection targets, and activation profiles during distinct innate social behaviors.

Using a panel of 31 genes, we profiled 643,834 cells and classified MeA neurons into 16 inhibitory and 10 excitatory neuronal cell types, many of which have distinct spatial patterns. Furthermore, hamFISH has allowed us to identify a rare inhibitory cell type (Ndnf⁺) localized in the molecular layer of the MeA. We then combined hamFISH with MeA output mapping and found cell-type biases in projections to the medial preoptic area (MPOA), bed nucleus of the stria terminalis (BNST), and the ventrolateral portion of the ventral medial hypothalamus (VMHvl). Finally, by combining hamFISH with c-fos activity mapping upon interactions with conspecific males, females, or pups, we found that instead of specific MeA transcriptomic types being activated by unique social stimuli, most of these types were combinatorially activated, generating unique MeA population activities. Our in situ profiling of the MeA, therefore, provides a proof-of-principle example of overlaying transcriptomic types, projection, and activity in a behaviorally relevant manner and demonstrates the usefulness of hamFISH in multiplexed in situ gene expression profiling.

Results

Robust detection of RNAs in the brain using highly amplified multiplexed FISH

To achieve robust transcript detection under the high autofluorescence observed in mouse brain tissues, we first tested branched DNA amplification⁹ with single-molecule FISH. We designed preamplifier and amplifier oligonucleotides that contain repetitive sequences that together multimerize to form branched DNA. In the original scheme (highly amplified version 0), the preamplifiers were designed to bind common non-complimentary overhangs (20 nt) on RNA-specific encoding probes (Figure 1A). Preamplifiers also include 9x binding sites (19 nt) for amplifier oligonucleotides which themselves contain 10x binding sites (18 nt) for fluorophore-conjugated readout probes (all sequences are compiled in Table S1). We found that this method increased signal intensity by 2.7-fold (807 ± 20 A.U. vs 2157 ± 63 A.U.). However, we observed a high number of punctate signals (39.2 ± 4.4 puncta/10 µm²) in the extracellular space (Figure 1B,C and Figure 1 – figure supplement 1A). To test if these extra-somatic signals are true signals or represent non-specific probe binding events, we designed probe sets that recognize two different regions of Parvalbumin (Pvalb) mRNA (n=12 probes for each), each with a distinct fluorescence readout (Figure 1 – figure supplement 1B). Extra-somatic signals had an overlap rate of 6.2 ± 1.0%, significantly lower than perinuclear signals in Pvalb⁺ cells (39.1 ± 2.4 %), suggesting that the bulk of these spots may be non-specific signals (Figure 1 – figure supplement 1C, D). To alleviate this issue, we modified the probe design so that preamplifier binding sequences were split into two 14-nt half sites on adjacent encoding probes (highly amplified version 1.0) similar to splitFISH¹⁰. Indeed, we found that the extra-somatic signal decreased when using this design (17.4 ± 1.4 puncta/10 µm², Figure 1B and Figure 1 – figure supplement 1A). Furthermore, there was an increase in the overlap rate of perinuclear Pvalb⁺ spots (46.5 ± 2.5%) relative to our first design suggesting an improvement in the fidelity of signal detection (Figure 1 – figure supplement 1D). Compared with non-amplified smFISH, this branched DNA amplification scheme increased smFISH signal by 3.8-fold (807 ± 20 A.U. vs 3072 ± 50 A.U., Figure 1C).

Figure 1

Separately, we developed a system with preamplifiers and amplifiers containing 32x repeats (see materials and methods for details) to enable the reliable detection of low-abundance transcripts. These amplifiers, termed version 2.0, increased the single molecule signals by 12.3-fold compared with the non-amplified signal (807 ± 20 A.U. vs 9849 ± 172 A.U.), and 3.2-fold over our highly amplified version 1.0 smFISH method (Figure 1 – figure supplement 2A). Moreover, we performed rolling circle amplification (RCA) mediated detection of Pvalb transcripts¹¹. Although the probe chemistries of these methods are significantly different, we found that the signal intensity from the version 2.0 method and RCA-based signals are comparable (Figure 1 – Figure Supplement 2B). These data therefore confirm that branched DNA amplification enables tuneable amplification of single molecule FISH signals. To enable multiplexing with the branched DNA trees, we designed additional amplifiers and tested them using different encoding probes. Curiously, not all designed amplifiers produced specific signals with high dynamic range (n=13/32). We therefore replaced the amplifier sets by screening new sequences and selecting only those that produced a high signal-to-noise ratio in a more standardized fashion using Sst mRNA as a common target (see materials and methods for details).

To utilize encoding probes amplified from oligopools, which offer a significant cost advantage over individually synthesized probes, we modified our original scheme to include bridge-readout probes (Figure 1D). These gene-specific probes create an intermediate layer that links encoding and preamplifier probes. Importantly, while oligopool probes are synthesized in a pool and fixed in composition, bridge-readout probes are individually synthesized and can be customized for each experiment. This scheme therefore allowed us to flexibly select a subset of genes for readout without resynthesizing oligopool libraries.

The amplifiers were then tested for their orthogonality. We used an encoding probe library in which each gene is coupled to a unique amplifier set, confirming that the incorporation of the bridge-readout probe was successful at detecting marker genes in the MeA (Figure 1E, Figure 1 – figure supplement 3). Upon imaging signals from all the amplifiers sequentially, we performed single-molecule localization and examined signal overlaps. We then excluded those that showed cross-talk, resulting in an orthogonal set of 32 amplifiers (Figure 1F). Our amplified smFISH design therefore enables transcript detection in brain tissue that is compatible with multiplexed spatial profiling methods.

Application of hamFISH for profiling a subnucleus of the amygdala

Compared to in situ transcriptomics methods using barcoding at a single molecule level, non-barcoded, sequential readout of multiplexed FISH signals does not necessitate the precise alignments of single molecule signals from imaging rounds. In addition, amplified signals using hamFISH enabled us to speed up the imaging time by using a lower magnification objective (40x as opposed to 60x). Multiplexing by fluorescent readouts using AlexaFluor-488, Cyanine 3 (Cy3) and Cyanine 5 (Cy5) further reduced imaging time.

This scalability makes sequential hamFISH an attractive method to be coupled with other physiological and connectivity readouts, where increased numbers of samples are inevitable due to the lower throughput nature of physiological and neural tracing methods. As an important case study, we applied hamFISH for cell type profiling linking gene expression, connectivity, and activity in the medial amygdala (MeA). The MeA is a critical hub for the neural circuits underlying innate social and defensive behaviors in mammals. While the behavioral roles of the MeA have been well established and specific MeA subpopulations have been linked to the execution of specific innate behaviors^12,13, circuit-level mechanisms of how each MeA cell drives specific physiological and behavioral effects remain incompletely characterized. Addressing cell type diversity in multimodal fashion is an important step in delineating the circuit basis of how these neurons control a repertoire of behavior that includes mating, aggression, and parenting.

To choose a panel of genes that captures cell type diversity in the MeA, we first re-clustered previously published scRNA-seq MeA data⁸ and found 26 clusters in cells that are Slc32a1 (inhibitory) or Slc17a6 (encoding Vglut2 and therefore excitatory) positive (Figure 2 – figure supplement 1A, Table S2). Eight of these clusters we suspected of either being located outside the MeA, based on Allen Brain Atlas in situ data (clusters 6, 11, 14), or expressed ubiquitously throughout the MeA and surrounding areas (clusters 4, 5, 20, 21, 26). In addition, we found three clusters that did not have specific marker genes that were differentially expressed or seemed to be a mix of excitatory and inhibitory cells (clusters 16,17,25), leaving 15 MeA cell clusters. Therefore, we assembled a gene panel consisting of marker genes for these scRNA-seq clusters and those used in previous studies. We finalized 31 genes for subsequent hamFISH experiments (Table S3).

We then used sequential hamFISH to profile cells along the full anterior-posterior axis of the MeA (at 100 µm interval over 1 mm). Initially, experiments were performed on 13 animals with a total of 163,925 cells profiled (mean of 12,607 ± 351 per animal). Of these, 40.3 ± 0.5% were inhibitory and 39.0 ± 1.4% excitatory, as defined by the expression of Slc32a1 and Slc17a6 respectively. The remaining cells were made up of either cell types that are localized outside of the MeA and did not express either Slc32a1 or Slc17a6 (8.6 ± 0.7%) or were MeA cells whose identity could not be established due to the absence of any markers other than Slc32a1 or Slc17a6. Nevertheless, cells in which the expression of three or more marker genes is detected made up ∼80% of the cells in the MeA. This is consistent with the scRNA-seq data showing that approximately 20% of the MeA cells are glia. Our 31-gene panel can therefore be used to define most of the neurons in the MeA.

We then clustered our hamFISH data for de novo discovery of cell types and identified 16 inhibitory (i1-i16) and 10 excitatory (e1-e10) neuronal classes, along with 4 classes of cells that expressed neither Slc32a1 nor Slc17a6 (n1-n4) (Figure 2A-C, Figure 2 – figure supplement 2A). The clusters were consistent between animals and sexes (Figure 2 – figure supplement 2B). To determine cluster reproducibility between different profiling methods, we correlated the mean expression of our gene panel markers in our hamFISH clusters with clusters from Chen et al. (2019) scRNA-seq data (Figure 2 – figure supplement 1B). We found reasonable correspondence between the datasets. For example, scRNA-seq clusters 1, 7 and 8 contained Slc17a6, Adcyap1, Nos1, Penk, and Tac1 positive cells, which correspond to hamFISH clusters e1-e6. The excitatory classes fall approximately into 3 classes: e1-6, e7-9 and e10 when compared with sequencing data. There were also hamFISH clusters that were not found in the scRNA-seq dataset. We attribute them to the clusters with very low number of cells (i3, i4) or low expression of marker genes that may not be reliably detected by scRNA-sequencing (i7-i9). In more recently published scRNA-seq data¹, however, we find these clusters except i3, perhaps due to the deeper sequencing coverage.

Figure 2

We then proceeded to analyze the spatial distributions of these clusters. Consistent with previous work¹⁴, most inhibitory classes are localized in dorsal and ventral subdivisions of the MeA, whereas excitatory neurons occupy primarily the ventral MeA (Figure 2D, Figure 2 – figure supplement 2C, Figure 1D). We found, for example, that the posterodorsal MeA (MeApd) contained six unique Ecel1⁺Prlr⁺ inhibitory clusters (i11-i16) with specific spatial patterns. For example, among the clusters marked by the expression of Cartpt and Cck, i15 (Esr2⁺) and i16 (Esr2^-) occupy the most posterior part of the MeApd, and these clusters are flanked laterally by i13 (Calcr⁺) neurons and intermixed with the sparser i14 (Dlk1⁺) population (Figure 2B). Rostral to these posterior MeApd populations are two subtypes that do not express either Cartpt or Cck but are positive or negative for Calb1 (i11 and i12, respectively). Excitatory classes also show specific spatial structures. For example, Cck⁺ excitatory cells can be split into a Dlk1⁺ population (e8) that solely occupies the most ventral part of the anterior MeA (MeAa), and a Dlk1^- population (e9), which sits more laterally. Other striking spatial classes of MeA neural subtypes include sparse inhibitory Sst-expressing populations (i5, i6) that span the entire MeA in a salt-and-pepper pattern and Ndnf⁺ inhibitory cells (i3) that line the medial border with the optic tract, or ventral surface of the MeA.

Our MeA hamFISH data also include cell populations outside the MeA based on the discontinuity of marker expression. For example, Penk⁺Pdyn⁺ (i2) cells primarily occupy a location lateral to the MeA, possibly the central amygdala (CeA), with a smaller proportion that appeared to localize within the MeA. There are four clusters located in the CeA (n2) and posterior cortico-medial amygdala (PMCo, n1, n3 and n4) which are Slc32a1^- and Slc17a6^- (n1-n4). These are putatively Slc17a7⁺ excitatory cells based on scRNA-seq and Allen Brain Atlas in situ data (mouse.brain-map.org). hamFISH therefore enabled the spatial mapping of at least 16 inhibitory and 10 excitatory neuron types.

Moreover, to compare differences in the spatial distribution of each cluster and quantify reproducibility between experiments, we coarsely aligned sections from different experiments and clustered cell classes using metrics that describe the spatial components of each cluster (Figure 2E). These were: mean x (medial-lateral axis) position, standard deviation (SD) of x, mean y (dorsal-ventral axis), SD of y, mean z (anterior-posterior axis), SD of z, total number of cells, distance to the molecular layer (the input layer of the MeA), and the distribution of clusters that harbor nearest neighboring cell to each neuron in consideration. The clustering revealed that inhibitory and excitatory classes generally have different spatial properties (Figure 1E, left), although the salt-and-pepper, sparse nature of e10 (Nts⁺) cells is more similar to inhibitory cells than other excitatory classes. Within the inhibitory classes, i3 (Ndnf⁺) and i4 (Gal⁺) are separated into distinct clusters, consistent with a more restricted spatial pattern compared with other inhibitory clusters (Figure 1E, middle). Within the excitatory classes, e6, e8, e9 and e10 were separable from other clusters, consistent with more ventral positions (e6 and e8), lateral position (e9) and sparse salt-and-pepper position (e10) compared with other excitatory cells (Figure 1E, right). These data therefore show that many of the transcriptomic cell clusters have distinct spatial patterns. By combining our amplified smFISH design with a panel of 31 genes imaged sequentially, our data provide a comprehensive in situ map of the major transcriptomic neural clusters in the MeA.

hamFISH uncovers a rare population of inhibitory neurons in the MeA

Among the three classes of MeA GABAergic neurons that have been characterized in vitro, only one type, which is localized to the molecular layer, harbors short axons with fast-spiking characteristics, and has therefore been proposed to be a local interneuron^15,16. Strikingly, our in situ map of MeA neuron types identified a single class of cells preferentially localized to the MeA molecular layer (i3), and we decided to characterize this population further.. These cells express Ndnf and account for 95% of inhibitory neurons localized there. Notably, these cells are rare, comprising only 1.8% of neurons in the MeA, and we failed to identify this subtype in the reference scRNA-seq data. These cells also express Slc32a1 and Reln, but not Npy or Nos1, and therefore resemble a specific class of cortical L1 interneuron that has been linked to gain modulation¹⁷, plasticity, and cross-modal integration¹⁸. As the acetylcholine receptor a7 (encoded by Chrna7) is also expressed in GABAergic L1 neurons¹⁹ we tested if the remaining population of GABAergic cells in the molecular layer correspond to a different type. Unexpectedly, we found that 74.5% of Ndnf⁺ cells in the molecular layer of the MeA express also Chrna7 (Figure 2 – supplemental figure 3A), indicating that there is further heterogeneity within the Ndnf⁺ population.

To further characterize the connectivity of MeA^Ndnf neurons, we mapped their inputs and outputs. Using Ndnf-Cre animals in which Cre-dependent AAV harboring synaptophysin-EGFP and tdTomato was injected into the MeA, we identified the axonal projections of these neurons (Figure 2 – figure supplement 4A). We found GFP signal in four major brain areas: MeA, AOB, cortical amygdala (olfactory cortex) and BNST, indicating that, unlike L1 Ndnf⁺ neurons, MeA^Ndnf neurons harbor both local and long-range axonal arbors (Figure 2 – figure supplement 4B,C). MeA i3 neurons especially send a strong long-range input to the AOB, including mitral and granular cell layers, and to the BNST. Conversely, mapping of the inputs to i3 neurons by rabies tracing showed that they receive inputs from various brain areas including olfactory areas, cortex, and hippocampal formation (Figure 2 – Figure supplement 4D-F). Notably, a similar retrograde tracing experiment using Sst-expressing neurons in the MeA as starter cells showed these neurons have different input structures. In particular, MeA^Ndnf neurons receive more inputs from olfactory areas and less inputs from the dorsal striatum compared to MeA^Sst neurons. Altogether, these anatomical data establish that MeA^Ndnf neurons represent a unique neural subtype embedded in the circuit processing both main olfactory and vomeronasal inputs.

The projection specificity of transcriptomic subtypes of MeA neurons

Combining gene expression profiling of neural tissue with other biological readouts such as connectivity is vital to understand how neural circuits function. Our simple method to profile transcriptomic types in situ facilitates further multimodal interrogation of these cells. MeA neurons send projections to both proximal and distal targets^14,15,20. We sought to determine the projection targets of specific cell classes by combining hamFISH with multicolor retrograde tracing. We injected fluorescently labeled cholera toxin subunit b (CTB) into three limbic regions known to receive strong MeA inputs and play important roles in social behaviors: the medial preoptic area (MPOA), the posterior medial bed nucleus of stria terminalis (BNST), and the ventrolateral subdivision of ventromedial hypothalamus (VMHvl) (Figure 3A, Figure 3 – figure supplement 1A). Initial quantification of labeled cells in the MeA revealed differences in the distribution of projecting cells along the anterior-posterior axis (Figure 3B, Figure 3 – figure supplement 2A). Specifically, VMHvl-projecting neurons were enriched in the anterior over posterior, with a significant decrease at bregma -1.7 mm. In contrast, both BNST- and MPOA-projecting cells were distributed evenly through the anterior-posterior axis. We also compared the total number of CTB-positive cells per target area in the MeA and found that BNST-projecting cells were nearly 2-fold more abundant than MPOA- and VMHvl-projecting cells (Figure 3 – figure supplement 2B). Furthermore, triple color retrograde tracing in single animals allowed us to identify a minority of double- or triple-labeled cells, indicating that some MeA neurons target multiple brain regions (Figure 3 - figure supplement 2C, D). Approximately 75% of neurons were singly labeled, which is unlikely due to inefficient labeling as the rate of triple BNST-injected positive neurons was near 100% (Figure 3 – figure supplement 1B). These observations suggest specific connectivity of each MeA cell to downstream targets.

Figure 3

We then classified the labeled cells into the transcriptomic classes as defined above. Using hamFISH we assigned 263,474 CTB-labeled cells from 13 animals to neuronal cell types (Figure 3C, Figure 3 – figure supplement 3). These experiments uncovered two major insights into MeA projection neuron specificity. First, at the mesoscale level, we found a significant enrichment in VMHvl-projecting neurons that were excitatory (Figure 3D). In fact, of all VMHvl-projecting cells, 80 ± 9 % were excitatory in contrast to 43 ± 4 % of all neurons assigned. BNST and MPOA-projecting cells were instead made up of a similar proportion of excitatory cells as non- projecting cells (46.7 ± 9 % for BNST-projecting cells, and 49 ± 10 % for MPOA-projecting cells). Second, we found CTB labeling in most of the transcriptomic clusters, indicating most projection neurons project to all three areas examined. There was, however, enrichment of labeling in specific clusters for certain projection targets (Figure 3E, Figure 3 – figure supplement 2E). We found that VMHvl-single-projecting neurons, for example, were over-represented in the excitatory clusters defined by the marker genes Adcyap1, Nos1 and Penk (e.g., e2-e6), when compared with BNST and MPOA-single projecting cells. Broadly, MPOA and BNST-projecting neurons showed similar cell type enrichment to each other despite our findings that BNST/MPOA-double-projecting neurons were relatively low in abundance (6-13%) (Figure 3 – figure supplement 2D). Although MPOA and BNST-projecting neurons were assigned to excitatory cell types in the ventral MeAa, a greater enrichment could be seen in inhibitory Ecel1⁺Prlr⁺Lhx6⁺ MeApd clusters (clusters i8, i11-i16) when compared with VMHvl-projecting cells. Interestingly, BNST-projecting cells also showed enrichment in three excitatory clusters (e1, e4 and e8) when compared with MPOA-projecting cells, with e8 (Cck⁺Dlk⁺) cells also enriched compared with VMHvl-projecting cells (Figure 3 – figure supplement 2E). In addition, we identified several subtypes that are enriched in cells sending projections to multiple targets (Figure 3 – figure supplement 2F,G). For example, the cluster i16 was more abundant in double-projecting BNST/VMHvl and MPOA/VMHvl cells compared with single VMHvl- projecting cells. Conversely, clusters e3 and e6 were enriched in double-projecting BNST/VMHvl and BNST/MPOA cells compared with single BNST-projecting cells. These results therefore uncover a mixed but biased projection profile of individual transcriptomic MeA populations.

Mixed selectivity of the activation patterns of MeA transcriptomic types upon innate social behaviors

We next tested the response of MeA neurons to social stimuli to determine the input organization to each MeA cluster. We prepared MeA sections from the animals exposed to different conspecifics in a freely behaving fashion. For this experiment, we also included probes against the immediate early gene c-fos, which has been commonly used to survey activity in this brain area^14,21. Virgin male mice in their home cage were presented with either another male, a female, or pups, so that neuronal activation could be correlated with aggression, mating, and infanticide, respectively (Figure 3F). We found that while many of the transcriptomic classes can be commonly activated during social behaviors, the pattern of activation appears to be stereotyped and different among different behavioral contexts (Figure 3G and Figure 3 – figure supplement 3A). We found only a few clusters selectively activated during specific social behavior (Figure 3 – figure supplement 3B). The excitatory cluster e6 (Tac1⁺) was more highly activated in both aggressive (14.4 ± 1.2 % of e6 cells were c-fos⁺) and infanticidal (15.0 ± 2.2 %) behaviors compared with reproductive behavior (10.1 ± 0.7 %, p<0.01 vs male, p<0.001 vs pup, paired t-Test), whereas e8 (Dlk1⁺) cells showed the opposite preference, being more highly activated during mating behavior (25.4 ± 3.0 % of e8 cells were c-fos⁺ when exposed to a female conspecific compared with 2.8 ± 1.1 % in male exposed samples (p<0.0001), and 5.0 ± 2.3 % in pup exposed samples (p<0.05)). Activation of inhibitory clusters was similar during different behaviors, although i14 and i16 were more activated during aggressive behavior than infanticide. Taken together, these results suggest that different social stimuli activate the ensemble of transcriptomically distinct MeA neurons in unique patterns.

Discussion

Development of highly amplified multiplexed FISH as a routine method for multiplexed transcript detection

Various in situ methods have been developed to capture the transcriptomic complexity of neurons while preserving spatial information, e.g., MERFISH^22–24, split-FISH¹⁰, seqFISH^25,26, STARmap¹¹, coppaFISH²⁷, Hyb-ISS²⁸, EASI- FISH³ and BARseq2²⁹, each offering unprecedented opportunities for understanding cell type heterogeneity within tissues. Furthermore, commercial platforms, e.g. MERSCOPE (Vizgen), Xenium (10x Genomics) and CosMx (NanoString) facilitate in situ transcriptomic studies using many (>100) genes, however, the high costs of these platforms preclude them from being used for routine multi-modal survey of cell types. In contrast, platforms such as coppaFISH²⁷, EASI-FISH³, SABER⁴, and RNAscope³⁰ allow cell type analysis with fewer genes (fewer than 100 but more than 5) and have been more amenable for multimodal analysis. These methods, however, pose challenges in implementing from scratch when commercial systems are not available.

Branched DNA amplification has been available for decades and is proven to be a robust platform for signal amplification in tissues⁹. Indeed, our initial tests using two complementary probe sets for a single gene confirm the detection of gene transcripts with rare dropouts. However, it remains challenging to extend this method beyond a dozen amplifier sets. We tested nearly 100 sequences to select specific amplifiers with high signal-to-noise and validated their orthogonality. This resulted in the development of 32 orthogonal amplifier modules. Our split encoding probe design is a feature that can be used in other methods such as splitFISH¹⁰ and HCRv3³¹ to ensure a low incidence of false positives. In addition, we designed and synthesized two different lengths of amplifiers to provide differential signal strength. Long amplifiers were used to read out c-fos to ensure the high dynamic range of expression of this gene. Short amplifiers (v1) produced sufficient signals for many of our target transcripts. In contrast, long amplifiers (v2) produced more signals, which some bled-through into unintended fluorescence channels under our filter setting.

Furthermore, the use of oligo pools to prepare probes rather than individually synthesized oligonucleotides that are commonly used for sequential detection of transcripts, helped reduce experimental costs. Different probe libraries can be independently amplified from single oligo pools that contain as many as 10,000 oligonucleotides. With the ‘bridge-readout’ scheme, one can further narrow down genes that should be visualized to optimize the marker gene sets. We also demonstrated that our probe scheme could be combined with other amplification systems such as hybridization chain reaction (HCR)^31,32 by incorporating HCR initiator sequences into bridge-readout probes. In addition, the imaging of each round was faster than HCR-based methods^3,33 which require stripping and polymerization steps. We estimate that the start-up reagent cost of hamFISH using 32x multiplexing is 7,000 GBP as of writing, which can accommodate approximately 140 experiments assuming that each sample covers 10 mm² tissue. The system is versatile so that one could choose to multiplex less, which further reduces the cost of the experiments. Taken together, these developments should democratize in situ cell type identification in individual laboratories.

Towards multimodal in situ map of MeA transcriptomic cell types

Previous behavior studies elucidated the functional specialization of unique MeA populations in specific social behaviors, e.g., MeApd GABAergic neurons are involved in aggression while glutamatergic neurons in MeApv are involved in self-grooming¹². In contrast, despite single-cell RNA sequencing of MeA neurons uncovering further transcriptomic heterogeneity in the MeA^7,8,34, the circuit role as well as behavior and physiological roles of most of these clusters remain uncharacterized. We therefore surmised that a constructive direction would be to overlay the anatomical and activity information on these transcriptomics datasets.

By using hamFISH as a workhorse method for in situ transcriptomics to survey over 600,000 cells, we learned that this focused approach complements unbiased sc/snRNA-sequencing screens. Our data revealed further subtypes that were missing from this dataset when compared to amygdala cell types found in scRNA-seq data available at the time of our data collection⁸, likely due to low abundance (i3) or those (i8, i9, i15) that are defined by genes with low expression levels (e.g., Lhx6, Esr2). Notably, we designed and collected the in situ data before new sets of data using 10x v3 platforms became available^1,7. Still, we have since found evidence of all these cell types in scRNA-sequencing data except for i3. Furthermore, our study failed to identify sexually dimorphic cell types that had been previously reported ⁷, indicating that our gene selection provides limited resolution.

The mapping of MeA transcriptomic types uncovered a striking spatial organization of transcriptomic types within MeA subregions. The neuronal types that occupy specific locations include e8 and e9, which are excitatory ventral Cck⁺ cells, as well as inhibitory i13, i15, and i16 cells which express Cartpt, Prlr, Cck and Ecel that mark the posterodorsal MeA. We also found that some clusters exhibit a salt-and-pepper distribution within a specific region, which might add further heterogeneity in a given population. Moreover, we also identified a rare cell type localized in the input layer of the MeA. Notably, Ndnf⁺ inhibitory neurons (i3) were not detected in the previous scRNA-seq data. Furthermore, we found expression of Chrna7 in a proportion of Ndnf⁺ cells, implying at least 2 different subtypes of cell make up the i3 population.

Recent studies identified two major developmental lineages of MeA neurons, Dbx1-derived and Foxp2- neurons^13,35,36. Although Dbx1 is not expressed in the adults and Foxp2 was not included in our gene panel, the Allen ABC atlas¹ shows that Foxp2 is broadly expressed across several of our hamFISH clusters, for example Sst⁺/Npy^- (i6), Htr2c (i7), Calb1 (i9), and seems to overlap with the Esr2^- MeApd population (i16). Therefore, given the roles of Foxp2 neurons in aggression¹³, it would be interesting to study the role of each of the Foxp2 subpopulations during conspecific male-male contact. Moreover, our projection mapping experiments uncovered that excitatory MeAa cells project to the VMHvl (Figure 3 – figure supplement 2G), an area that has been identified to play an important role in the experience-dependent escalation of aggressive behavior³⁷.

Whether mammalian pheromones activate dedicated neural circuits as observed in insect systems remains unresolved. Indeed, ESP1 and ESP22 have been shown to activate c-fos within MeA neurons in stereotyped patterns, some of which appear to mediate the action these pheromones in the control of lordosis and mating^38,39. In contrast, in vivo physiological recordings consistently found population coding within this area^40,41. To this end, we found no molecularly defined clusters uniquely activated during specific behavioral outputs, with a notable exception of e8. Importantly, the activation pattern of each cluster as a MeA ensemble appears unique in different social contexts, providing a potential molecular link to the neural representation of social information in the MeA^40,41. Our experiments also found the robust activation of glutamatergic MeA neurons during social behaviors, perhaps consistent with the activation of ventral MeA by single pheromonal compounds^38,39, but their circuit and behavioral roles need to be uncovered by specifically manipulating these populations rather than using pan-glutamatergic cell markers, the activation of which produces self-grooming¹². A prominent future direction therefore is to study a specific population of neurons based on molecular identity and projection targets, which would enable a systematic investigation of connections in this circuit.

Pitfalls of the study

A major pitfall of hamFISH is that this method is currently limited to thin section samples (<20 µm) due to the large molecular size of amplifiers, which creates a barrier to diffusion. HCR provides better scalability to thicker tissues in this respect. In addition, hamFISH currently allows only 32 genes to be sequentially read out, thus this method may prove inadequate for brain areas that might need more than 32 genes to capture transcriptomic diversity. While we chose to characterize the subclasses of MeA neurons, this may not be fine enough to fully capture the functional heterogeneity of the MeA. Future investigations therefore could use multiplexing strategies to profile the functional properties of MeA cells at their finest (i.e., “cluster”) level in the MeA.

Materials and methods

Animals and behavioral assays

All animal experiments were performed under the UK Animals (Scientific Procedures) Act 1986 (P108E8CF8) following local ethical procedures. Adult wild-type CD1 and C57BL/6J (Charles River, UK) mice aged 8-12 weeks were used in this study and kept on a 12h:12h light/dark cycle with ad libitum access to food and water. Only CD1 mice had tissue collected to be used for hamFISH experiments. Prior to behavioral assays for c-fos measurements, animals were individually housed for 7 days and then habituated in a testing room under dim red lighting for 3 hours before behavioral assays were performed (at 1-2h after the onset of the dark cycle). Animals were sacrificed 30 min after exposure to a conspecific before tissue was harvested. Behavioral assays that involved exposure to conspecifics were video recorded to ensure the selection of tissue from mice that displayed the appropriate behavioral responses (see below) for subsequent hamFISH experiments. The specific details of the behavioral experiments are described below.

Infanticidal behavioral assays – Pups were generated by mating CD1 males with C57BL/6 females. Infanticidal males were exposed to two freshly euthanized pups (P2-4) and tissue was selected for subsequent hamFISH assays if evidence of pup-directed aggression was confirmed (e.g., biting and eating).

Inter-male aggression assays – One male mouse (intruder) was taken from his cage and placed in the cage of another male mouse (resident). Tissue was used from the resident mouse if evidence of fighting was confirmed. This included bouts of biting and tail-rattling.

Mating assays – Male mice were introduced to the cage of a female mouse. Tissue was used from these mice if evidence of sexual behavior such as anogenital sniffing, mounting, and intromission was observed. The female oestrus cycle timepoint was confirmed by vaginal swab.

Probe and bridge-readout design

Encoding probes were constructed by amplifying oligonucleotide sequences from complex oligo pools (synthesized by Twist Bioscience), as described below. To design the gene-specific sequences, we first used the BLAST program as 30 nt probe candidates against the mouse RefSeq dataset, then discarded sequences with hits with length 14 nt or more. Using R scripts, pairs of 30 nt were then selected in a manner that best optimized the number of consecutive probes along a transcript with a target of 90 sequences per gene, although some genes did not reach this due to sequence space. The reverse complement of the chosen sequences was then flanked with bridge-readout sites designed in the following way.

To design the flanking bridge-readout sites, a set of orthogonal 25mer DNA barcodes⁴² were expanded to 28mer sequences with the addition of random nucleotides, before a further BLAST search to exclude off-target hits. Only those where each half of the sequence had a similar T_m were chosen. Selected bridge-readout sequences were then split into 14mer half sites and placed on either side of the gene-specific sequences described above. Sequences were tested as described below. For the bridge-readout probes themselves, the complementary sequences of the 28mer bridge-readout sequences were conjugated to four preamplifier binding sites (two sites 5′ and two sites 3′) and synthesized as Ultramers by IDT. Flanking the encoding probe combined sequences, we placed sequences necessary for probe production (PCR primer amplification sites and T7 promoter), as described below. All sequences for encoding probes can be found in Table S3.

Encoding probe preparation

The protocol for probe production is described on a Protocols.io webpage (https://www.protocols.io/view/highly-amplified-multiplexed-fluorescence-in-situ-ditu4enw). Of note, encoding probes were constructed with an improved protocol. Briefly, the oligo pools were amplified by limited-cycle PCR (40 ng per reaction) using KAPA HiFi polymerase (Roche, KR0369) (10 mM dNTP, 1x KAPA HiFi fidelity buffer) with the following protocol: 98°C for 3 min, then 22 cycles of 98°C for 10 s and 72°C for 30 s, followed by 72°C for 5 min. In total, four sub-libraries were amplified with the following four primer pairs; lib1_fwd (5′- GCGTTGCGTGCTAACTCGGA-3′) and lib1_rev (5′-TTACTAACCCGGTCGTGCGG-3′), lib2_fwd (5′- GAGTGCGGATACGTCGTCGT-3′) and lib2_rev (5′-AGGGACGCTTAAGTCACGGC-3′), lib3_fwd (5′- GCAGTTCGTGCGACCGTGTA-3′) and lib3_rev (5′-ATCGCACGGTTCGTAATCCG-3′), lib4_fwd (5′- TCCGTTAACGTTACGCGGTG-3′) and lib4_rev (5′-AGTGGCGCCGAATAAGCAAT-3′).

Using the HiScribe T7 High Yield RNA Synthesis Kit (New England Biolabs, E2040), the sense strand was amplified further and we then produced antisense ssDNA by reverse transcription (RT) with Maxima H Minus Reverse transcriptase (Thermo Fisher, EP0751). To increase encoding probe tissue penetration, RT primers contained RNA bases to ensure degradation of the shared RT sequences in the final encoding probes upon hydrolysis. The RT primer sequences used are as follows: lib1_RT_fwd (5′- GCGTTGrCrGTGCTAArCTCGGArU-3′), lib2_RT_fwd (5′-GAGTGCGrGrATACGrUCGTCGTrU-3′), lib3_RT_fwd (5′-GCAGTTrCrGTGCGACrCGTGTArU-3′), lib4_fwd_RT (5′- TCCGTTrArACGTTACrGCGGTGrU-3′) where rX denotes an RNA base. Our final probe length following RT sequence degradation was therefore 58-62 bases rather than 78-82. The RNA was then hydrolyzed (0.5M NaOH, 0.25M EDTA, incubated at 95°C for 10 min) and the reaction was neutralized with 0.33x volume of 1M Tris-HCl pH 7.0. We purified the ssDNA probes by ethanol precipitation (0.5x volume of 7.5M ammonium acetate and 2.5x volume of 100% ethanol) and resuspended the pellet in ultrapure water.

Amplifier design and screening

Sequences used for amplifier and readout oligonucleotides were based on a collection of orthogonal 25mer DNA barcodes⁴². For the sequences used for preamplifier (19 nt) and amplifier (18 nt) repeats we shortened these sequences and selected those with a Tm range of 50-55°C. Fluorescent-readout probes were 25 nt, despite binding sites on amplifier oligonucleotides being only 18 nt. This allowed us to strip fluorescent-readout probes during hamFISH experiments using 25 nt ‘toehold’ oligonucleotides that are complementary to the fluorescent-readout probes. Additionally, amplifiers were synthesized with a 3′ amino modification for covalent linkage to the hydrogel. We developed 32 orthogonal sets of amplifiers over three rounds of tests. First, we designed 16 amplifier sets (consisting of preamplifier, amplifier, and fluorescent-readout probes), which we found to produce signals in varying intensities. We subsequently tested additional 18 sets to compare signals with the original 16 sets using Sst as the common target. Signal intensity and signal-to-noise ratio were evaluated using a Zeiss AxioImager 2 using a Plan Apo 20x objective (NA 0.8) and Orca Flash v2 (Hamamatsu Photonics) camera. In a third iteration of screening, we devised a method to screen 48 candidates in a more systematic fashion. First, we screened the bridge sequences using a fixed amplifier set. This was done by using Sst encoding probes with candidate bridge sequences and bridge probes that contain a common preamp binding sequence for the readout. After selecting bridge sequences that produced highly specific signals, we screened the preamplifier repeat sequences using a common Sst encoding probe and a bridge probe that contains a candidate preamplifier sequence. We then used another oligonucleotide that contains a sequence complementary to this candidate sequence plus 4x of the amplifier binding site. The signal was read out using a common amplifier and a fluorophore-conjugated probe. Using the screened sequences, we designed new amplifier sets consisting of preamplifier, amplifier, and label probes (label probes were not screened). Finally, we tested these new sets using Gad1 encoding probes and corresponding bridge probes containing binding sites for 4 different amplifier sets plus an existing set as a standard. We then chose the new amplifier sets that perform comparably to the reference.

For synthesis of longer amplifiers for use in highly amplified v2 smFISH, preamplifier and amplifier sequences containing 4 amplifier-binding and readout probe-binding repeats respectively, were synthesized by overlapping PCR amplification of three primers using KAPA HiFi polymerase (Roche, KR0369) (10 mM dNTP, 1x KAPA HiFi fidelity buffer) with the following protocol: 95°C for 3 min, then 30 cycles of 95°C for 30 s and 65°C for 30 s, followed by 72°C for 3 min). The repeated sequences were also flanked with inverted nicking endonuclease (Nt.BbvCI) sites (underlined) and contained a T7 promoter sequence (italicized). Taq Polymerase (Invitrogen, K461020) was used for this, according to the manufacturer’s protocol, for its compatibility with TOPO cloning. This resulted in the following sequences: preamplifier_4x = CATTCGGTTGGGTATCGGTAATCAGCGTCCTCAGC(GCCTAAATTGACAGATTCC)₄GCTGAGGCCCTATAGTGAGTCGTATTA, generated from the primers preamp_fwd1 (5′- CATTCGGTTGGGTATCGGTAATCAGCGTCCTCAGCGCCTAAATTGAC-3′), preamp_fwd2 (5′- CCTCAGCGCCTAAATTGACAGATTCCGCCTAAATTGACAGATTCCGCCTAAATTGACAGATTCCGCCTAAATTGACAGATTCCGCTGAGG-3′) and preamp_rev1 (5′-TAATACGACTCACTATAGGGCCTCAGCGGAATCTGTCAA-3′), amplifier_4x = GGAATCTGTCAATTTAGGCCCTCAGC(GATCTGCTATCCGTTAGT)₄GCTGAGGCCCTATAGTGAGTCGTATTA (generated from the primers amp_fwd1 (5′- CTCTCAATGCGCAATTAAGCCTCAGCGAGCCACGTTGA-3′), amp_fwd2 (5′- CCTCAGCGAGCCACGTTGAATTATCGAGCCACGTTGAATTATCGAGCCACGTTGAATTATCGAGC CACGTTGAATTATCGCTGAGG-3′) and amp_rev1 (5′- TAATACGACTCACTATAGGGCCTCAGCGATAATTCAACGT-3′). These were cloned into pTOPO using the TOPO cloning method (Invitrogen, K461020).

Following plasmid DNA preparation and confirmation of successful clones by sequencing, the supercoiled portion of the DNA was separated from any nicked or single-stranded DNA by purification after electrophoresis on a 1% agarose gel (50 µg plasmid DNA loaded). Purified supercoiled DNA was then incubated with Nt.BbvCI (New England Biolabs, R0632) according to the manufacturer’s protocol followed by column purification. To duplicate amplifier repeats within the nicked plasmid DNA, we incubated the samples with DNA polymerase I, Large (Klenow) fragment (New England Biolabs, M0210) which fills in 5′ overhangs (left by the nicking endonuclease) in the presence of single-stranded binding proteins (SSBP) to stabilize secondary structure (2 µg DNA, 1x NEBuffer 2, 1.6 mM dNTPs, 5 µg SSBP, 2.5 units Klenow, 37°C for 1h). The resulting single-stranded (filled-in) DNA was then gel purified and ligated with T4 DNA ligase (New England Biolabs, M0202) before transformation into Stbl3 competent cells (Thermo Fisher, C737303). Plasmid DNA that had been successfully nicked and filled-in therefore had amplifier 8 repeats, which was confirmed by sequencing. This process was repeated 2 more times, resulting in plasmid DNA with 32 repeats of either amplifier-binding or readout probe- binding sequences.

To synthesize single stranded amplifier sequences from these plasmids, DNA was digested with EcoRI to separate the amplifier template from the plasmid backbone. This template was gel purified and used as a template for in vitro transcription (New England Biolabs, E2040). The resulting RNA was reverse transcribed with Maxima H Minus Reverse transcriptase (Thermo Fisher, EP0751). Note that for the amplifier sample, dNTPs were supplemented with 100 nM aminoallyl-dUTP (Thermo Fisher, R0091) for incorporation into the tissue hydrogel during gel embedding. The RNA template was removed by alkaline hydrolysis, as described in the probe preparation methods above, and the remaining ssDNA was neutralized with 0.33x 1M Tris-HCl pH7.0, then column purified (Omega Biotek, D6943).

Tissue preparation for hamFISH experiments

Tissue was harvested by freezing freshly dissected brains of animals in OCT compound (Tissue-Tek, 4583). For animals injected with cholera toxin subunit b, the lobe containing the left MeA was cut away from the rest of the brain (which was drop-fixed in 4% v/v paraformaldehyde (Thermo Fisher, 28908) in 1x phosphate-buffered saline (PBS) (Severn Biotech, 20-74-10) and used for injection site confirmation, see below) before freezing in OCT. Embedded brains were stored at -80°C. Prior to sectioning, 40 mm diameter #1.5 coverslips (Bioptechs, 0420- 0323-2) were treated with 0.5% v/v (3-mercaptopropyl)trimethoxysilane (Sigma, 175617) and 0.34% glacial acetic acid in ethanol, and poly-L-lysine (Sigma, P8920) for silanization and promotion of tissue-adhesion respectively.

Frozen brains were manually sectioned onto these coverslips at -18°C on a cryostat (Leica, CM3050S). Surrounding tissue was trimmed to ensure approximately 18-22 sections (9-11 from each animal, 2 animals per coverslip) could fit onto the imaging area of each coverslip. Sections were cut at 10 µm thick with every 10^th section mounted, ensuring representative sections from the full anterior-posterior axis of the MeA was present (Bregma - 1.2 mm to -2.2 mm).

Samples were then fixed with 4% v/v paraformaldehyde, 0.4% v/v glyoxal (Sigma, 50649), 0.1% v/v methanol in 1xPBS⁴³ followed by treatment in pre-chilled 100% ethanol at -80°C for 30 min. Subsequent drying out of these samples was followed by rehydration in 2x saline-sodium citrate (SSC) (Severn Biotech, 20-6500-50) and clearing in 4% wt/v sodium dodecyl sulfate (SDS) (Sigma, 75746) in 1x PBS for 5 min each, on a rotating platform. After three quick washes in 2xSSC, sections were washed with 1x TE (Severn Biotech, 20-5310-05) (2 washes at 1 min each) followed by permeabilization using protease (0.34 µg/ml) (Sigma, P5380) in 1x PBS for 7.5 min. Samples were then washed again in 2x SSC followed by a final wash in 2x SSC containing 20% v/v formamide (Sigma, F7503).

Encoding probe staining, amplification, and gel embedding

The molarity per gene was estimated based on the number of genes in the probe library, despite some differences in the number of probes per gene. For example, a probe library containing 30 genes with a total concentration of 30 µM ssDNA was defined as 1 µM per gene and diluted accordingly. Encoding probes were diluted in hybridization buffer; 2x SSC, 22.5% v/v formamide (Sigma, 47671), 5% wt/v dextran sulfate (Millipore, S4031), 0.5% v/v TWEEN-20 (VWR, 437082Q), 10 mM ribonucleoside vanadyl complexes (New England Biolabs, S1402), 0.1% wt/v yeast RNA (ribonucleic acid from Torula utilis, Sigma, 83850). hamFISH encoding probe libraries were used at a final concentration of 25 nM per gene as defined above, and individually synthesized probes at 100 nM per gene. A polyA-anchor probe containing locked nucleic acids (LNA) as described in ref²³ was also included at a concentration of 1 µM (/5Acryd/TTGAGTGGATGGAGTGTAATT+TT+TT+TT+TT+TT+TT+TT+TT+TT+T). This probe permits RNA anchoring to a polyacrylamide gel as described below.

Coverslips were inverted onto a drop of hybridization buffer (approximately 160 µl) and placed onto a layer of parafilm sitting in a 10 mm petri dish, with a 0.75 mm thick spacer to ensure even hybridization. The sample was incubated at 40°C for 16-18 h before being washed in 2x SSC, 20% v/v formamide (wash buffer), initially at room temperature for 3 quick washes, followed by incubation for 1h at 40°C. For signal amplification samples were incubated in sequential hybridization of bridge-readout, preamplifier, and then amplifier oligonucleotides in amplifier hybridization buffer (2x SSC, 20% v/v formamide, 10% wt/v dextran sulfate) for 30 min at 40°C. This was interspersed with 2x 5 min washes in wash buffer on a rotating platform. Bridge-readout and blocking probes were used at a concentration of 1 nM per gene, whereas preamplifiers and amplifiers were used at 10 nM per set. After the amplification and wash steps, samples were fixed again for 10 min with 4% v/v paraformaldehyde in 1x PBS followed by a 2x SSC wash, then either frozen at -80°C after drying or processed for gel embedding and signal development (see below).

For gel embedding, samples were treated as previously described¹¹. Briefly, coverslips were washed in 20 mM acrylic acid NHS ester (Sigma, A8060) for 90 min, followed by a 30 min incubation in monomer buffer; 2x SSC, 4% acrylamide and 0.2% bis-acrylamide (Bio-Rad Laboratories). Aspirated samples were then inverted onto a 40 µl drop of polymerization mixture (monomer buffer with the addition of 0.2% ammonium persulfate and 0.2% tetramethylethylenediamine) upon a Gel Slick (Lonza, 50640) treated glass plate. After polymerization, coverslips were removed from the glass plate while submerged in 2x SSC, and we proceeded with digestion with 1 mg/ml proteinase K (Ambion, AM2548) in 2x SSC, 2% wt/v SDS, 0.5% triton-X (Sigma, X100) for 3 h at 40°C on a rotating platform.

Imaging platform

Samples were imaged on a custom-built imaging platform based on that described by^22,44 with 405, 488, 568 and 647 nm solid-state lasers (OBIS, Coherent). The base of the microscope was Nikon Ti-U (Nikon) equipped with a custom FocusLock system as described in²². The coverslip was assembled in a flow chamber (Bioptechs, FCS2) attached to a home-built fluidics system that controlled hybridization flow through the sample. The fluidics system comprised of a peristaltic pump (Gilson, MINIPULS 3) and three eight-way valves (Hamilton, MVP RS-232 valves), all computer-controlled by HAL/Dave and Kilroy software (https://github.com/ZhuangLab/storm-control/blob/master/storm_control/README.md) ensuring fully automated control of the whole system.

Imaging

The flow chamber was filled with imaging buffer during acquisition (2x SSC, 50 mM Tris-HCl pH8.0, 10% wt/v glucose (Sigma, G8270), 2 mM 6-Hydroxy-2,5,7,8-tetramethylchromane-2-carboxylic acid (Sigma, 238813), 250 µg/ml glucose oxidase (Sigma, 345386), 40 µg/ml catalase (Sigma, C30) and 40 U/ml RNase inhibitor (New England Biolabs, M0314)). Each hybridization readout buffer (2xSSC, 40% v/v formamide, 5% wt/v dextran sulfate) contained AlexaFluor-488, Cy3 or Cy5 fluorophore-conjugated readout probes (50 nM each) and toehold probes (1 µM each) corresponding to the previous imaging round’s fluorescent readout probes. The hybridization protocol for each round was as follows (at an approximate flow rate of 0.15 ml/min); 600s readout and toehold hybridization solution, 1800s stop (incubation time for hybridization), 350s 2x SSC, 100s stop, 180s 2x SSC with 8 ng/µl DAPI, 180s 2x SSC, 150s stop, 200s imaging buffer. Note that a DAPI step was included in our protocol to supplement DAPI signals that became weaker in each round. For barcoded hamFISH, a 60x Plan Apo 1.4NA oil objective (Nikon) was used for imaging. A Z-stack for each channel was acquired at each field of view, covering a thickness of 7 µm with 14 steps at 500 nm each. Typically, 30-40 FOVs covering 1.25-1.65mm² were needed to fully image each MeA section. The imaging area for each coverslip (containing the MeA from 1 animal, 12 sections from each) was therefore approximately 17mm². A 40x Plan Apo 1.3NA oil objective (Nikon) was used for imaging. A Z-stack for each channel was acquired at each field of view, covering a thickness of 5 µm with 6 steps at 1 µm each. Typically, 12 (4x3) FOVs covering 1.25mm² were needed to fully image each MeA section. The imaging area for each coverslip (containing the MeA from 2 animals, 10 sections from each) was therefore approximately 25mm².

Stereotaxic injection of retrograde tracers

Cholera toxin subunit b (CTB) -AlexaFluor488 (Thermo Fisher, C34775), CTB-AlexaFluor555 (C34776), and CTB- AlexaFluor647 (C34778) were used for retrograde tracing of MeA-projection neurons. Briefly, adult male CD1 mice were anesthetized in an induction chamber containing 4% isoflurane mixed with oxygen. After shaving the head, animals were then transferred to a stereotaxic frame and anesthesia was maintained with ∼1% isoflurane delivered through a nose cone. A small incision was made and a hole was drilled at each of the injection sites. CTB (0.5%, 75nl) was delivered through pulled glass micropipettes to the target location using Nanoject III (Drummond Scientific); MPOA (AP -0.1 mm, ML -0.5 mm, DV +5.4 mm), BNST (AP -0.2 mm, ML -1.1 mm, DV +4.3 mm), VMHvl (AP -1.4 mm, ML -0.7 mm, DV +5.9 mm), PMv (AP -2.5 mm, ML -0.5 mm, DV +5.8 mm). After injections, the scalp was sutured and animals were allowed to recover in their cage placed on a heated platform for 1h. Mice were euthanized by cervical dislocation 5-7 days after CTB injection. After dissection of the lobe containing the left MeA (which was frozen in OCT compound), the rest of the brain was drop-fixed in 4% v/v paraformaldehyde 1x PBS overnight at 4°C, before vibratome sectioning at 75 µm thickness and mounting with DAPI-containing mounting medium: 80 ng/µl DAPI (Sigma, D9542), 0.2% v/v n-propyl gallate, 90% v/v glycerol, 1x PBS. Injection sites were confirmed using Zeiss Axio Scan imaging.

Prehybridization imaging of CTB signal in the MeA

The left MeA lobe (previously frozen in OCT, to be used for FISH) was sectioned as described above. After fixation in 4% v/v paraformaldehyde, 0.4% v/v glyoxal, 0.1% v/v methanol in 1xPBS, and staining with 2x SSC 80 ng/µl DAPI, samples were washed in 2x SSC containing RNase inhibitor (40 U/ml, New England Biolabs, M0314) then mounted onto slides. Sections were then imaged using a Zeiss Axio Imager 2 using the 20x objective and Orca Flash v2 (Hamamatsu Photonics) camera, with a 4-step Z-stack covering 12 µm thickness. After imaging, slides were submerged into 2x SSC to loosen the coverslip without damaging the sample, and the tissue preparation protocol described above was followed starting with the ethanol incubation step.

Transcriptomic analysis pipeline

To analyze the previously published scRNA-Seq data⁸, clustering was performed using Seurat⁴⁵. We set the number of dimensions to 10 and a resolution of 2. After quality control (set as requiring 200-2500 RNA per cell), we selected the top 250 genes that showed significant variation for clustering.

For hamFISH experimental data, we built a smFISH analysis pipeline that creates maximum intensity projections from each imaging round (ImageJ). To correct the illumination bias of the excitation beams in the field of views, the images were corrected using a blurred average image compiled from all tiles across all hybridization rounds for each channel to calculate intensity differences. These images were then used to normalize intensity differences for DAPI images. This was then combined with DAPI-based Cellpose 2.0 cell segmentation⁴⁶. To account for slight shifts in the stage position across hybridization rounds, DAPI signal was used to calculate the translation in XY for every tile. This was done using the phase cross-correlation of scikit-image¹⁰ (https://github.com/khchenLab/split-fish). The modal value for each hybridization round was then taken and applied to all images for that round, as well as to Cellpose masks to define individual cell positions across the experiment.

A Fast Fourier transform (FFT) based frequency filter (200Hz) was applied to signal images (i.e., non - DAPI channels) to exclude background haze. These then underwent shading correction as described above for DAPI. The intensity of signal within each cell (chosen by a manually defined intensity cut-off) was then used to define positive and negative cells for each transcript. The intensity value below the threshold was subtracted, and the intensity was divided by the standard deviation of each bit to adjust for the differences in variance in each bit. This data is therefore comparable to normalized cell x gene data. To test for orthogonality, the centroids of each transcript signal were found using 3D-DAOSTORM⁴⁷ and overlap was determined to be positive if two or more bits were within 0.25 pixels. For CTB-injected smFISH experiments, the DAPI signal from hamFISH experiments was aligned to the prehybridization CTB-imaging DAPI signal using a custom Python script that performed affine transformation to account for distortion (cv2.warpAffine). We then manually verified the alignment of each section. All code for hamFISH data processing is available at https://github.com/mathew-d-edwards/hamFISH.

Clusters were determined using a combination of iterative Seurat-based clustering and manual confirmation with hamFISH signal overlays. Briefly, Seurat was used to define an initial set of clusters which were then checked by manual overlays of signal from differentiating genes confirm clusters. After further manual overlays with other genes, higher-resolution cluster definitions could be determined. After cells that fit into these determined clusters were removed from the analysis, the remaining cells were subject to further Seurat-based clustering to try and determine their type. This iterative process gave an initial list of clusters that we then tried to find in multiple datasets and allowed us to ensure clusters were consistent in number across animals and experiments. A final cluster definition was then used to assign all cells from every experiment into clusters. The data was then analyzed and presented using R, Sklearn (tSNE), Scanpy (violin plots)⁴⁸, Seurat (heatmaps), GraphPad Prism, and Python-based plotting scripts for cluster spatial plots. Comparison of scRNA-Seq clusters with hamFISH clusters was done using Pearson correlation between group averages of the expression of each transcript in our gene marker panel.

For spatial analysis of hamFISH clusters, we defined several parameters that were used to cluster cell types based on their spatial features. To do this we identified common sections across samples, first in z by manually defining anterior-posterior Bregma position based on cluster locations. Then in x and y by manually defining the ventral and medial edge of the MeA, as well as the most ventral-medial point that the MeA reaches, and by using the angle between all three to align sections between samples. Spatial features used to cluster cell types were: mean x (medial-lateral axis) position, standard deviation (SD) of x, mean y (dorsal-ventral axis), SD of y, mean z (anterior- posterior axis), SD of z, total number of cells, distance to molecular layer (manually defined in reference images), proportion of cells with nearest neighbor cells belonging to cluster i1, and proportion of cells with nearest neighbor cells belonging to cluster i2 etc. for all clusters.

Input/output tracing of Ndnf and Sst neurons

Following craniotomies as described above, the following AAV or rabies viruses were used for injection into the MeA (AP -1.6 mm, ML +/- 2 mm, DV +5.6 mm) of 2-3 month old NDNF-Cre or SOM-Cre male mice (C57BL/6J mice backcrossed to CD1 background). For monosynaptic retrograde rabies tracing experiments, 20-60 nl of AAV2/1.syn.flex.H2BG.N2CG and AAV2/1.flex.GFP-TVA (6.94E+13 vg/ml and 3.84E+13 vg/ml, respectively; Francis Crick Institute viral vector core) were injected unilaterally into the MeA followed two weeks later by 20- 60nl of N2C△GEnvA.Rabies.mCherry (1E+09 vg/ml; Sainsbury Wellcome Centre viral vector core (SWC-VVC)). To determine the projection targets of Ndnf cells, 100nl of AAV2/2.EF1a.FLEX.Synaptophysin-EGFP and AAV2/2.CAG.FLEX.tdTomato mix (6.93E+13 vg/ml and 5.77E+13 vg/ml, respectively (SWC-VCC)) were injected into the MeA of 8-12-week-old male CD1 mice. Tissue was harvested 4-5 weeks after injection. After perfusion, dissection, and post-fixation of the brain in 4% v/v paraformaldehyde 1x PBS overnight at 4°C, brains were embedded in agar and imaged using serial section two-photon microscopy^49,50. Imaging was controlled by ScanImage Basic (Vidrio Technologies, USA) with BakingTray, a custom software wrapper for setting up imaging parameters (https://github.com/SainsburyWellcomeCentre/BakingTray). Sections were cut at 50 µm and the signal was imaged with 780 nm excitation wavelength. Subsequent image assembly was performed using StitchIt (https://github.com/SainsburyWellcomeCentre/StitchIt). Sections were then registered to the 10 µm Allen Brain Atlas using BrainReg software^51,52. For the analysis of our input mapping, cellfinder software was used to automatically detect input cells⁵³. Although a low error rate of detection was seen for most areas, the main and accessory olfactory bulbs consistently gave an unacceptable error rate due to high background levels and inaccurate registration and therefore input cells were manually counted in these areas. For the analysis of our output mapping, the puncta-like signal from Synaptophysin-EGFP did not allow for automatic cell detection with cellfinder. The serial section two-photon microscopy images were therefore used to identify regions of interest, before re-mounting and imaging these areas using a Zeiss Axio Scan. These images were registered to the Allen CCFv3 atlas using AP_histology software (https://github.com/petersaj/AP_histology) and the signal was then manually counted as puncta using ImageJ to give a density value for each area of interest.

Hybridization Chain Reaction (HCR)

All HCR sequences can be found in Table S4. Briefly, Ndnf, Sst and Chrna7 encoding probes were designed with overhangs complementary to bridge sequences that served as initiators for hairpin amplification. HCR was performed using the same tissue preparation protocol used for hamFISH (see above) except that 20 µm thick sections were used instead of 10 µm, and samples were fixed for 10 min at room temperature in 1x PBS with 4% v/v paraformaldehyde. After tissue preparation, the manufacturer’s protocol (Molecular Instruments) was used to develop the HCR signal with the following alterations. We used encoding probes at a final concentration of 10 nM for our initial overnight incubation. On the following day a bridge probe hybridization step was performed in HCR Hybridization buffer for 30 min at 37°C with a bridge probe concentration of 10 nM. The section was washed with 0.2x SSC to remove the unbound bridge probe. To develop HCR signal, after denaturation (95°C for 90 sec) and cooling, B2-549 and B4-647 (30 nM each) were hybridized in an Amplification buffer (Molecular Instruments) overnight at room temperature before final washes, mounting and imaging.

Supplemental Figures

Supplemental Figure 1 - figure supplement 1

Supplemental Figure 1 - figure supplement 2

Supplemental Figure 1 - figure supplement 3

Supplemental Figure 2 - figure supplement 1

Supplemental Figure 2 - figure supplement 2

Supplemental Figure 2 - figure supplement 3

Supplemental Figure 2 - figure supplement 4

Supplemental Figure 3 - figure supplement 1

Supplemental Figure 3 - figure supplement 2

Supplemental Figure 3 - figure supplement 3

Acknowledgements

We thank the Sainsbury Wellcome Centre Neurobiological Research Facility for technical support; the Sainsbury Wellcome Centre Advanced Microscopy Facility for help with serial 2-photon microscopy; and Johannes Kohl for manuscript feedback. Fig 1 contains a panel that used Biorender. This work was supported by the Sainsbury Wellcome Centre core grant from the Gatsby Charitable Foundation and Wellcome Trust (090843/F/09/Z).

Additional information

Author contributions

Additional files

Amplifier sequences

scRNASeq clusters

Sequential hamFISH probe sequences

HCR sequences

Note

This reviewed preprint has been updated to use the correct text; previously, a version prior to the one reviewed was presented here.