All surgical and experimental procedures were in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and approved by the Animal Care and Use Committees of the Allen Institute or Johns Hopkins University. Allocations of mice to experiments were randomized and experimenter blind. Mice were housed at 20–22 °C and 38–42% humidity.

Definition of LC-NE neurons

We used Dbh-Cre mice⁸⁰ (Dbh^{tm3.2(cre)Pjen}, The Jackson Laboratory, 033951; RRID:IMSR_JAX:033951) backcrossed with C57BL/6J mice (The Jackson Laboratory, 000664; RRID:IMSR_JAX:000664). Dbh-Cre mice were crossed to RCL-H2B-GFP mice (The Jackson Laboratory, 036761; RRID:IMSR_JAX:036761), and the native fluorescence of the labelled nuclei was imaged on a SmartSPIM microscope after tissue clearing using LifeCanvas active delipidation, agar embedding, and refractive-index matching using EasyIndex (SmartSPIM, LifeCanvas Technologies). Images from 8 brains were acquired at 1.8 × 1.8 × 2.0 μm per voxel. Raw image data were passed through a customized 3D UNet (Trailmap)⁸¹ with a model trained for immunolabelled nuclei⁸². The resulting probability maps were convolved with a spherical kernel approximating a nucleus diameter and then centroids were estimated with 3D local maxima detection. Stitching and registration to the CCFv3 (RRID:SCR_020999) followed the SmartSPIM Pipeline (https://github.com/AllenNeuralDynamics/aind-smartspim-pipeline).

Points were subsequently transferred to the Allen atlas reference space via an ANTS transform (RRID:SCR_004757)⁸¹. After restricting analysis to caudal midbrain and dorsal pons, and after reflecting all points into a single hemisphere, local densities were estimated via nearest-neighbour analysis and points at uniform density levels were used to define meshes. In brief, normals for points were estimated and used to generate surfels which were wrapped in a watertight mesh⁸³ by isosurface extraction (https://www.github.com/fwilliams/point-cloud-utils). A mesh representing LC-core at the 67th percentile represents the density threshold at which sub-coeruleus cleanly separates from the main portion of LC. Initial mesh estimates were made using dynamic radius and resolution parameters to compensate for point density and mesh size across the displayed percentiles.

Single-neuron reconstructions

Whole mouse brain specimens were prepared for single-neuron reconstructions using previously published methods^33,84. In brief, 2 male and 2 female adult Dbh-Cre mice received systemic injections, via the retro-orbital sinus, of a 100 μl mixture of Cre-dependent Tet transactivator (AAV-PHP-eB_Syn-FlexTRE-2tTA, Addgene plasmid id: 191210; RRID:Addgene_191210; dosage range 1.0 × 10⁸ genome copies (gc) ml⁻¹ to 3.0 × 10⁹ gc ml⁻¹) and a reporter virus (either AAV-PHP-eB_7x-TRE-3xeGFP or AAV-PHP-eB_7x-TRE-tdTomato, typical dose 1.8 × 10¹¹ gc ml⁻¹; RRID:Addgene_191206 and RRID:Addgene_191207). Viruses were obtained from either the Allen Institute for Brain Science viral vector core, the University of North Carolina, or the BICCN-Neurotools core and were prepared in an adeno-associated virus (AAV) buffer consisting of 1× phosphate-buffered saline (PBS), 5% sorbitol, and 350 mM NaCl. Six-to-eight weeks after viral transfection, mice were anaesthetized with an overdose of isoflurane and transcardially perfused with 10 ml 0.9% saline at a flow rate of 9 ml min⁻¹ followed by 50 ml 4% paraformaldehyde in PBS at a flow rate of 9 ml min⁻¹. Brains were extracted and post-fixed in 4% paraformaldehyde at room temperature for 3–6 h and then left at 4 °C overnight (12–14 h). The following day, brains were washed in 1× PBS to remove all traces of excess fixative. Subsequent tissue processing including clearing, immunolabeling, and whole-brain expansion steps were carried out as previously described^33,84. For immunolabeling, delipidated brains were equilibrated in a detergent buffer (PTxw) followed by incubation in PTxw containing 20 µg per brain in 4.5 ml of the primary antibody, either rabbit anti-GFP (ab290, Abcam) or goat anti-tdTomato (AB8181, SICGEN). After thorough washing in PTxw, the corresponding secondary antibody was applied at a dose of 40 µg per brain in 4.5 ml; donkey anti-rabbit Alexa Fluor 488 (A-21206, Invitrogen) for GFP-labelled brains or donkey anti-goat Alexa Fluor 568 (A-11057) for brains labelled with tdTomato. All other steps were carried out exactly as described previously^33,84. Gelled brains were soaked in 0.05× saline sodium citrate (SSC) to achieve approximately 3× expansion, and 24 h prior to imaging the expanded brains were equilibrated in a solution of 0.05× SSC that contained 10 mM ascorbic acid included as an antifade.

Expanded brains were imaged on a custom SPIM microscope (ExA-SPIM)^33,85. Brains were imaged with both 488 nm and 561 nm excitation and an 8× binned autofluorescence image volume was collected for the ‘non-signal’ channel for registration to the CCF. All subsequent data processing steps including image illumination correction, stitching and fusion into a coherent image volume were carried out via automated cloud pipelines³³.

ExA-SPIM whole mouse brain image volumes were registered to CCFv3. A custom ExA-SPIM template was generated by aligning and averaging whole-brain images (8× binned, autofluorescence channel) from 7 brains, including flips (14 specimens total). Registration proceeded in two steps. First, an automated two-stage registration was performed using Advanced Normalization Tools (ANTs; https://github.com/ANTsX/ANTs), consisting of a per-sample affine and SyN deformable registration to the ExA-SPIM template, followed by a fixed ExA-SPIM-template-to-CCF registration. Second, local misalignments remaining after the automated step were manually refined via landmark-based registration, matching anatomical landmarks across the automatically aligned sample and the CCF reference, using 3D Slicer³⁴. The resulting displacement field, computed by the automated sample-to-template and template-to-CCF transforms, was applied to the node coordinates of each neuronal reconstruction (SWC) via ANTs point transformations, thereby mapping reconstructed neurons from sample space into CCF space.

To generate whole-brain single-neuron reconstructions, we used HortaCloud⁸⁶, an open-source streaming 3D annotation platform enabling fast visualization and collaborative proofreading of terabyte-scale image volumes. Using this platform, human annotators proofread individual neurons via a web browser on a personal workstation. Proofreading entailed starting from the soma and tracing out all axonal and dendritic segments using a depth first search tree traversal approach^34,87. Neuronal trees generated by manual placement of control points were refined offline to generate dense (one node per voxel) and sparse (uniform sub-sampling of the dense tree) representations of the reconstruction data before use for subsequent quantitative analysis. Quantitative comparisons were made with neurons in the MouseLight database (RRID:SCR_016668).

Models of single-neuron axonal distributions

We developed two models to generate whole-brain distributions of axons to compare to the distribution of LC-NE axons (Fig. 2 and Extended Data Fig. 2d). We first create an axonal projection data matrix where each row is a reconstructed neuron, and the columns are coarse brain regions. Each entry in the matrix is the fraction of that cell’s reconstructed axon that lies within that coarse brain region. The coarse brain regions we used were: olfactory areas (OLF), isocortex (CTX), hippocampal formation (HPF), cortical subplate (CTXsp), cerebral nuclei (CNU), thalamus (TH), hypothalamus (HY), midbrain (MB), cerebellum (CB), pons (P), medulla (MY), cerebellum related fibre tracts (cbf), lateral forebrain bundle system (lfbs), medial forebrain bundle system (mfbs) and ventricular systems (VS).

The first model (‘random projections’; Extended Data Fig. 2d, middle) posits that each LC-NE neuron projects randomly across the brain. We assess this model by shuffling the projection matrix across brain regions and comparing the shuffled data with the original. As expected, the shuffled data loses the structured correlation matrix of the original data. This model ignores important anatomical constraints of the brain, namely that an axon cannot go straight from the pons to isocortex.

The second model (‘Markov with renewal’; Extended Data Fig. 2d, right) posits that LC-NE neurons distribute their axons by a Markov process walk across the brain, constrained by connections between brain regions and the minimum and maximum axonal lengths. We first created a transition matrix between brain regions. Using the axonal reconstructions, we iterate across cells and axonal segments (approximately 10 μm in length). We tally the brain region where that segment starts and stops. After aggregating across all cells, we normalize the matrix into a transition matrix where each entry represents the probability of transitioning from one brain region to another. An axonal branch that meanders near a region boundary might artificially inflate the rate of transitions between regions. We thus only counted transitions between two brain regions if the axon stayed in the new region for a minimum of 2 mm.

Each simulated neuron starts in the LC (pons). The axonal arbor is generated in 10-μm steps. On each step, the axon transitions to a new brain region according to the transition matrix. The axons can also branch or terminate, based on rates learned from the data, which depended on axon length from the soma but not the identity of the current brain region. Over the first 10 cm from the soma, both rates increase with the branching rates larger than the termination rate. After 10 cm from the soma, both rates are roughly constant and equal. Importantly, each branch grows independently according to the same Markov process. We limit the minimum axon length (summed across branches) to 10 mm to match the shortest neuron in our dataset. A Weibull renewal process stops all active axon processes. The Weibull process has a hazard function that is dependent on the overall axon length of the cell (summed across branches) and then results in a Weibull distribution of axonal length. We fit a Weibull distribution to the LC-NE neurons in our dataset, using neurons with complete reconstructions as uncensored observations, and neurons with incomplete spinal cord projections as censored observations.

This simple growth model produces axonal arbors that are qualitatively similar to the data. Model axons under-project to the isocortex and over-project to the thalamus and cerebellum. The correlation matrix matches the anterior/posterior block structure observed in the data. The grouping of axons into coarse brain regions is imperfect, especially around narrow fibre tracts. This explains why the model creates a stronger correlation between the cerebellum and brainstem than is observed in the data.

Somatodendritic properties

We quantified somatic bipolarity by fitting a 3D Gaussian distribution to a local patch of the image volume around each cell’s soma centre. The square root of the largest eigenvalue of the covariance matrix (normalized length of the principal axis) was taken as the somatic bipolarity coefficient. A similar approach was applied to dendrites, by first calculating a unit vector in the direction of each dendritic ‘stem’ point (primary branches within 50 μm or points where each dendrite crosses a 50 μm radius from the soma centre), then finding the normalized length of the principal axis of the set of dendritic stem direction vectors. We also calculated the angular separation between the somatic and dendritic primary axes, to show that each independently captured the neurons’ bipolarity. The harmonic mean of these bipolarity metrics was used as a combined bipolarity coefficient.

Retrograde tracing

Heterozygous Dbh-Cre mice (19 female, 16 male) crossed with homozygous Ai65 mice (The Jackson Laboratory, 021875; P40–P50; RRID:IMSR_JAX:021875) were anaesthetized with 2.0% isoflurane in O₂ and maintained at 0.7–1.0% isoflurane throughout the surgery. Mice were mounted in a stereotaxic frame (David Kopf Instruments) for injections of AAV-pEF1a-DIO-FLPo-WPRE-hGHpA (Addgene 87306-AAVrg, lot v188462; RRID:Addgene_87306; Supplementary Table 1). Small craniotomies were made in the skull above the injection target for intracerebral injections (using iontophoretic injections⁸⁸ and pressure injections⁸⁹ as described). Spinal cord injections were made between C4–C5 vertebrae without drilling⁹⁰. Intracerebral injections were in the right hemisphere. Spinal cord injections were bilateral.

Mice were euthanized with an isoflurane overdose 4–6 weeks after AAV injection and transcardially perfused with PBS followed by 4% paraformaldehyde (PFA; EMS 50-980-495 or equivalent⁹¹). Brains were dissected, post-fixed in 4% PFA overnight, and transferred to PBS. Tissue was cleared using LifeCanvas active clearing, followed by agar embedding, refractive-index matching using EasyIndex solution, and imaging (SmartSPIM, LifeCanvas Technologies).

Retrogradely labelled LC-NE neurons were identified by tdTomato fluorescence in the soma and primary processes. Positions of retrogradely labelled cells in whole-mount stitched and CCF-registered samples were identified either by manual labelling or automated segmentation followed by manual verification. Labelling accuracy was examined in at least two distinct planes (coronal, horizontal).

For automated segmentation, brain volumes were divided into 3D chunks (512 × 512 × 512) optimized for OME-Zarr image processing. Each block underwent background subtraction using the photutils photometry package⁹². Processed images were passed to a modified version of the CellFinder algorithm⁹³, where a median filter followed by a Laplacian of Gaussian filter were applied to each slice along the imaging plane. A threshold 6× s.d. above the mean was used to identify potential cell locations and the image was binarized. 2D regions were merged via an ellipsoid filter and large regions, resulting from the merging of cells, were split using iterative ellipsoid filters. To minimize false positives from fluctuations in tissue autofluorescence, cell proposals from the detection phase were passed to an 18-layer ResNet for classification. For each proposal a small block (28 × 28 × 50) centred on the proposed cell location was taken from the signal, as well as the autofluorescence channel. The input image to the network was downsampled by a factor of two across each dimension prior to classification to minimize the contributions of fluorescent signal from axonal and dendritic processes. Interactive Neuroglancer images were created using proposals classified as true cells, and a final round of manual refinement was done to remove any remaining false positives and add false negatives.

Annotations were registered to CCFv3 using the affine and warp transforms calculated during image registration. Annotations related to the LC were identified using a signed distance function comparing point locations within the CCF relative to the exterior mesh of the pons. All annotations were categorized as existing inside, outside, or on the boundary of the mesh. Retained points inside and on the pons mesh boundary were then manually examined to remove false-positive cells and add back in false negatives.

BARseq and MAPseq

The sindbis virus barcode library HZ120⁹⁴ was generated by the MAPseq core facility at Cold Spring Harbor Laboratory (CSHL) and used for MAPseq and BARseq experiments as described previously⁹⁵. The HZ120 library exhibits a diversity of approximately 8 million barcodes and was not fully sequenced in vitro.

Two 8-week-old C57BL/6J male mice were anaesthetized using oxygenated 4% isoflurane and maintained with oxygenated 0.7–1.0% isoflurane throughout the surgery. We injected 300 nl of HZ120 sindbis virus in LC in each hemisphere at anterior–posterior (AP) −5.2; medial–lateral (ML) 0.85; dorsal–ventral (DV) −3.20, −2.90 using pressure injection (Nanoject III⁸⁹). Mice recovered and were euthanized 22–28 h post-injection.

Brains were dissected⁹⁶, embedded in OCT and snap-frozen in an ethanol dry ice bath. Extruded spinal cords were flash frozen, straightened out on a razor blade secured in a conical tube. Samples were stored at −80 °C until cryosectioning (Leica 3050S) in the coronal plane⁹⁷. We cut 300-μm slices for regions outside LC (MAPseq) and 20 μm for tissue containing LC (BARseq). To avoid cross-contamination, we used a fresh unused part of a blade to cut each slice and cleaned the brush and the holding platform with 100% ethanol between slices. Cryostat chamber temperature was allowed to equilibrate for 30 min to account for temperature changes between cutting thick and thin slices.

For MAPseq experiments, 3× 300-μm coronal sections were mounted in a row onto Superfrost Plus Gold slides. For BARseq experiments, 4× 20-μm sections were mounted onto Superfrost Plus Gold slides according to microfluidic chamber template. Once all cut sections were melted onto the slide, slides were rapidly frozen on dry ice and stored at –80 °C in slide boxes in vacuum sealed bags until microdissection (MAPseq) or library prep (BARseq).

Regions of interest were dissected from each chilled brain slice with microscalpels pre-chilled (Fine Science Tools 10316-14). Glass slides were kept on dry ice and placed on a pre-chilled metal platform mounted above a dry ice/ethanol bath. Sections were first bisected down the midline. Then major anatomical regions (for example, isocortex, hippocampus or midbrain) were cut along their respective boundaries. Each cut was performed with a cleaned blade. Blades were rinsed in consecutive 4× 50 ml conical tubes containing ethanol stored on ice. To keep blade temperature consistently low, 3–4 scalpels were rotated and stored on the metal platform. After right hemisphere dissection was completed, a reference image for dissection boundaries was taken again for each of the 3 sections. After that, the left side was dissected along the same boundaries. Individual tissue pieces were then collected into pre-chilled microtubes (Qiagen 19560 tubes, 19566 caps) on dry ice, capped and stored at −80 °C. Left and right hemisphere samples were kept separate, and corresponding regions of interest were combined across 3 consecutive 300-μm sections mounted on the same slide. To assess cross-sample contamination, we collected tissue from uninjected brains as negative controls. After all regions of interest were collected, with samples kept on dry ice we added pre-chilled Qiagen bead (69989) to each tube. Spinal cord samples then received 800 μl Trizol, while the rest used 400 μl Trizol (Thermofisher 15596026). Samples were shipped overnight on dry ice to the CSHL MAPseq facility.

RNA extraction and sequencing library preparation were performed following the protocol in ref. ⁹⁸, as described previously^37,95,99. In brief, total RNA for each region was extracted using Trizol (Thermo Fisher Scientific, 15596018) and eluted in 13 μl of H₂O. Before library preparation, a few randomly selected samples were run on a Bioanalyzer (Agilent, 5067-1513) to ensure good RNA quality. To prepare the library for sequencing, we mixed 4 μl total RNA of each sample with 1 μl spike-in RNA (GTCATGATCATAATACGACTCACTATAGGGGACGAGCTGTACAAGTAAACGCGTAATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNNNNNNNNNNNNNNNNNNNATCAGTCATCGGAGCGGCCGCTA CCTAATTGCCGTCGTGAGGTACGACCACCGCTAGCTGTACA), where ATCAGTCA is the barcode tag of the spike-in. Spike-in RNA was transcribed in vitro by T7 RNA polymerase and diluted into 10⁴ molecules per μl. Barcode mRNA was reverse transcribed into the first-strand cDNA using gene-specific barcoded reverse transcription primers and SuperScript IV (Thermo Fisher Scientific, 18090010). The barcoded reverse transcription primer sequence was 5′-CTTGGCACCCGAGAATTCCAXXXXXXXXXXXXZZZZZZZZTGTACAGCTAGCGGTGGTCG-3′, where X12 are the barcoded unique molecular identifiers (UMI) and Z8 are barcoded sample specific identifiers (SSIs). The synthesized first-strand cDNA of each sample was labelled with a 12-nt UMI (unique for each RNA molecule) and 8-nt SSI (unique for each sample). After synthesis of the first-strand cDNA, every 7–8 samples of target brain regions were pooled together for clean-up by 1.8× AMPure XP beads (Beckman Coulter, A63881) and synthesis of the second strand cDNA (Thermo Fisher, A48571). The double-stranded cDNA samples were cleaned up by 1.8 x AMPure XP beads and treated with Exonuclease I (New England Biolabs, M0293S) to remove the single-stranded reverse transcription primers before PCR amplification.

The double-stranded cDNA samples were amplified by two rounds of polymerase chain reaction (PCR) using standard Accuprime Pfx protocol (Thermo Fisher Scientific, 12344032) with 2 min extension for each cycle. Primers 5′-CTGTACAAGTAAACGCGTAATG-3′ and 5′-CAAG CAGAAGACGGCATACGAGATCGTGATGTGACTGGAGTTCCTTGGCACCCGAGAATTCCA-3′ were used for the first PCR reaction, and primers 5′-AATGATACGGCGACCACCGA-3′ and 5′-CAAGCAGAAGACGGCATACGA-3′ were used for the second PCR. In the first PCR, the pooled cDNA samples were amplified for 15 cycles in 250 μl of reaction volume per 7–8 samples. After treatment with Exonuclease I (New England Biolabs, M0293S) to remove excess primers, all the PCR1 products of the target areas were pooled together. A quarter of the pooled PCR1 products of target areas were amplified by 10 cycles in 12 ml of PCR II reaction volume. The PCR2 products were cleaned up and concentrated by SV wizard PCR clean-up kit (Promega, A9282) and loaded into a 2% agarose gel for electrophoresis. The 233 bp PCR product was cut from the agarose gel, cleaned up by Qiagen MinElute Gel Extraction Kit (Qiagen 28606) and quantified by Bioanalyzer (Agilent, 5067-4626) and quantitative PCR.

The pooled library was sequenced on two lanes of Illumina Nextseq550 high output at paired end 36 with Illumina compatible Read 1 Sequencing Primer and Illumina compatible small RNA read 2 primer as described previously³⁷.

MAPseq data processing

The raw MAPseq data consist of FASTQ files, where paired end 1 covers the 30-nt barcode sequence and paired end 2 covers 12-nt UMI and the 8-nt SSI³⁷. The raw sequencing data were processed with code available at https://github.com/ZadorLaboratory/mapseq-processing. We first assembled the reads from the paired-end input, trimming each to the length of the sequence of interest. These sequences were then aggregated by exact duplicate and the read count for each unique sequence kept. We removed any data with ambiguous bases or long runs (>7) of the same base. We labelled and categorized all the reads by SSI sequence and expected spike-in sub-sequence. To account for expected replication, amplification, and sequencing errors, we identified all barcodes that differed from each other by an edit distance of 3 or less. This parameter is informed by the barcode sequence length, and the known diversity of the virus library. All groups of barcodes were then collapsed to the most common variant, as measured by read count. We then aggregated by barcode, SSI and UMI sequence, retaining the number of UMIs in each. To exclude low-confidence projections, we required a barcode to have at least 5 molecule counts in at least one projection region. We normalized the variation in library preparation and sequencing across different samples within each brain based on recovered spike-in RNA as described previously⁹⁹.

BARseq library preparation

BARseq samples for in situ sequencing were prepared using a protocol described previously^{95,99,100,101} with modifications¹⁰². In brief, brain sections were fixed in 4% PFA in PBS for 1 h, dehydrated through 70%, 85%, and 100% ethanol, and incubated in 100% ethanol for 1.5 h at 4 °C. After rehydration in PBST (PBS and 0.5% Tween-20), slices were incubated in the reverse transcription mix (RT primers, 20 U μl⁻¹ RevertAid H Minus M-MuLV reverse transcriptase, 500 μM dNTPs, 0.2 μg μl⁻¹ BSA, 1 U μl⁻¹ RiboLock RNase Inhibitor, 1× RevertAid RT buffer) at 37 °C overnight.

On day two, cDNA was cross-linked using BS(PEG)9 (40 μl in 160 μl PBST) for 1 h, the remaining crosslinker neutralized with 1 M Tris pH 8.0 for 30 min, followed by a wash with PBST. Sections were then incubated with non-gap-filling ligation mix (1× Ampligase buffer, padlock probe mix for endogenous genes, 0.5 U μl⁻¹ Ampligase, 0.4 U μl⁻¹ RNase H, 1 U μl⁻¹ RiboLock RNase Inhibitor, additional 50 mM KCl, 20% formamide) for 30 min at 37 °C and 45 min at 45 °C. Then sections were incubated in gap-filling ligation mix (same as the non-gap-filling mix with the Sindbis barcode padlock probe [XCAI5] as the only padlock probe, and with 50 μM dNTPs, 0.2 U μl⁻¹ Phusion DNA polymerase, and 5% glycerol) for 5 min at 37 °C and 45 min at 45 °C. After the second round of ligation, samples were washed with PBST, hybridized with 1 μM RCA primer (XC1417) for 10 min in 2× SSC with 10% formamide, followed by two washes in 2× SSC with 10% formamide and twice in PBST, then incubated in the RCA mix (1 U μl⁻¹ phi29 DNA polymerase, 1× phi29 polymerase buffer, 0.25 mM dNTP, 0.2 μg μl⁻¹ BSA, 5% glycerol (extra of those from the enzymes), 125 μM aminoallyl dUTP) overnight at room temperature to complete rolling circle amplification. On day three, rolonies were cross-linked using BS(PEG)9 (40 μl in 160 μl PBST) for 1 h, the remaining crosslinker neutralized with 1 M Tris pH 8.0 for 30 min, followed by a wash with PBST.

BARseq sequencing and imaging

Sequencing was performed using Illumina MiSeq Reagent Nano Kit v2 (300-cycles) and imaged on a Nikon Ti2-E microscope with Crest xlight v3 spinning disk confocal, photometrics Kinetix camera, and Lumencor Celesta laser. All images were taken with a Nikon CFI S Plan Fluor LWD 20×0.7NA non-immersion objective. For each sample, first 7 sequencing cycles probe the selected gene panel, followed by 15 sequencing cycles for barcodes, followed by hybridization cycle. In each imaging round, a z-stack of 10 images centred around the region of interest was captured with a 1.5-μm step size at each field of view (FOV). Adjacent FOVs had an overlap of 24%.

A detailed protocol is available at ref. ¹⁰³. In brief, for sequencing the first gene cycles, the sequencing primer (YS220) was hybridized in 2× SSC with 10% formamide for 10 min at room temperature. Subsequently, the sample underwent two washes in 2× SSC with 10% formamide and twice in PBS2T (PBS with 2% Tween-20). Following this, the sample was incubated in MiSeq Incorporation buffer at 60 °C for 3 min, followed by a wash in PBS2T. The sample was then exposed to Iodoacetamide (9.3 mg vial in 2.5 ml PBS2T) at 60 °C for 3 min, followed by a single wash in PBS2T and twice in Incorporation buffer. It was then subjected to two incubations in IMS at 60 °C for 3 min, followed by four washes with PBS2T at 60 °C for 3 min. Finally, the sample was washed and imaged in SRE. For sequencing the first barcode cycle, the sequenced products were initially stripped with three 10-min incubations at 60 °C in 2× SSC with 60% formamide. Subsequently, the sample was washed with 2× SSC with 10% formamide, and the barcode sequencing primer (XCAI5) was hybridized. The sequencing process then proceeded similarly to the first gene cycle.

For subsequent sequencing cycles for both genes and barcodes, the sample underwent two washes in Incorporation buffer, followed by two incubations in CMS at 60 °C for 3 min. This was succeeded by two more washes in Incorporation buffer, after which the sample was treated with iodoacetamide as in the first sequencing cycle. During hybridization cycles, the sequenced products were stripped with three 10-min incubations at 60 °C in 2× SSC with 60% formamide. The sample was then washed with 2× SSC with 10% formamide. Following this, the hybridization probes (XC2758, XC2759, XC2760, YS221) were hybridized for 10 min at room temperature. The sample was subsequently washed twice with 2× SSC with 10% formamide and incubated in PBST with 0.002 mg ml⁻¹ DAPI for 5 min. Finally, the sample was transferred to SRE for imaging. For one brain all the above steps were completed manually. For the subsequent samples in situ sequencing was performed on an automated microfluidics set up mounted on the microscope using the same reagents, temperatures and incubation durations.

Analysis of MAPseq and BARseq data

BARseq data processing

BARseq data were processed according to previous methods¹⁰⁰ with slight adjustments. Max-projections were generated from the image stacks, followed by noise reduction using Noise2Void¹⁰⁴. Subsequently, background subtraction, correction for channel shift and bleed-through, and registration of images across all sequencing cycles were applied. Cell segmentation was performed using Cellpose (RRID:SCR_021716)¹⁰⁵, while rolonies were decoded using BarDensr¹⁰⁶, and then assigned to cells. During BarDensr decoding, negative control gene identification indexes (GIIs; GIIs not associated with any padlock probe) were utilized to estimate the false discovery rate (FDR), with the decoding threshold automatically adjusted to target approximately 5% FDR. Whole-slice images were generated by stitching individual imaging FOVs. Processing steps were applied separately to each FOV to prevent stitching errors, minimize alignment artefacts, and facilitate parallel processing. Stitched images were solely used to generate a transformation, applied to each rolony and cell after completing all other steps.

Overlapping cells from neighbouring FOVs were identified using a custom implementation of sort and sweep collision detection algorithm¹⁰⁷. In cases where two or more cells overlapped, cells with higher read counts (indicating higher read quality) were retained. Each stitched slice was registered to Allen CCFv3 using manual registration. QuickNii (v.3 2017; RRID:SCR_016854) was used to manually select the CCF plane for each slice, followed by alignment of area borders within each slice using Visualign (v. 0.9) using nonlinear adjustments to the coronal plates.

BARseq transcriptomic analyses

MATLAB (RRID:SCR_001622) outputs from the BARseq data processing pipeline (.mat, v7.3/HDF5) were imported into R by reconstructing the sparse gene-by-cell count matrix from the stored CSC components and assembling per-cell metadata into a SingleCellExperiment object. Cells were assigned stable unique identifiers by concatenating batch number, slice, and cell ID. Negative control genes were removed from the data matrix. When duplicated gene symbols were present, the duplicate entry corresponding to a hybridization cycle was removed.

To identify transcriptomic types, an iterative clustering pipeline adapted from single-cell RNA sequencing studies¹⁰⁸ was used. Cells were QC filtered prior to normalization and clustering to ≥20 total counts and ≥5 detected genes to ensure correct transcriptomic identity assignment. Counts were normalized to CP10 (counts per 10) values and then subjected to log1p transformation. As the panel consisted of pre-selected genes, highly variable gene selection was skipped, and PCA was conducted on all probed genes. Unsupervised clustering was performed on log-normalized expression matrix using PCA (30 components), UMAP on the PCA space (100 neighbours), and shared-nearest-neighbour graph clustering (k = 50 neighbours; Jaccard weighting) followed by Louvain community detection to assign cluster labels.

LC neuron isolation by iterative clustering and anatomical gating

LC-NE neurons were isolated by iterative clustering. First, all neurons were clustered using the PCA/UMAP/SNN-Louvain workflow described above. For each sample, the cluster enriched for canonical noradrenergic markers (Dbh, Th, Ddc and Slc18a2) was identified and subset. The subset was then renormalized to CP10, log1p-transformed and reclustered to increase resolution within the LC-enriched compartment. At each iteration, candidate LC and LC-NE clusters were evaluated using (i) marker enrichment and (ii) anatomical localization across slices using CCF coordinates. Clusters whose cells localized clearly outside the LC region were excluded as entire clusters.

As the final clean-up step to obtain high-confidence LC-NE neurons, spatial coherence and density filtering was applied after the last round of LC-NE refinement. Using CCF coordinates (ML and DV converted to μm by ×25) and a depth coordinate defined as slice index × 20 μm, k = 10 nearest neighbours were computed per cell in 3D space. Spatial coherence was defined as distance-weighted cluster purity: the fraction of neighbours sharing the same Louvain label, weighted by 1/(distance + ϵ) (ϵ = 10⁻⁶). Spatial density was defined as 1/(mean kNN distance + ϵ). Cells were retained if spatial coherence ≥0.05 and spatial density was at least the 5th percentile of the dataset-wide density distribution. Cells failing either criterion were excluded, yielding the final high-confidence LC-NE cell set used for all subsequent downstream transcriptomic and projectome analyses.

BARseq-MAPseq barcode matching

BARseq 15-cycle barcodes from high-confidence, barcoded LC-NE neurons were converted from numeric base calls to nucleotide strings using a fixed mapping (1 = G, 2 = T, 3 = A, 4 = C). MAPseq 32-nt barcodes were truncated to 15 nt to match BARseq length. For each BARseq cell, candidate MAPseq matches were identified by Hamming distance at thresholds of 0–3 mismatches. Matches were resolved conservatively by retaining only the minimum-distance hit per cell and accepting a match only if (1) a single best hit remained and (2) the matched MAPseq 15-mer was unique in the MAPseq dataset (15-mer collisions excluded). MAPseq projection tables were exported for each mismatch threshold. False-positive match rates were estimated by applying the same matching and resolver rules to (1) random 15-mer barcodes and (2) shuffled BARseq barcodes (with the 9th nucleotide fixed and the remaining positions permuted), averaging across 20 runs. For each mismatch threshold, the fraction of simulated barcodes that were selected as unique best matches was used as the estimated false-positive rate under the full resolver. All projection analyses used data from one Hamming distance mismatch as false positive rate was estimated under 1%.

After matching, we removed duplicate entries arising from segmentation/FOV overlaps or barcode reuse. Duplicate structure was detected in two ways: (1) identical MAPseq profiles across all ROIs and (2) repeated BARseq barcodes. All identical MAPseq profiles were the result of repeat BARseq barcodes. For these duplicate BARseq groups, segmented cell images were exported and manually reviewed. If the same barcode appeared in different structures (barcode reuse/mixed sources), all entries for that barcode were discarded. If entries were the same biological cell segmented more than once, the better-segmented instance (operationalized by manual review, with transcriptomic QC metrics such as total UMI counts) was retained and the others discarded.

Projection matrix construction and harmonization

Per-brain matched MAPseq projection tables were converted into analysis-ready projection matrices using brain-specific loader functions. First, matched barcodes were filtered using the curated blacklists derived from manual review of duplicate or ambiguous barcoded cells. Next, MAPseq ROIs were renamed using sample metadata to encode target identity, hemisphere it was dissected from and slide number. Then we matched uniquely barcoded MAPseq cells to BARseq metadata, and soma hemisphere was assigned from CCF coordinates.

To enable cross-brain data pooling, per-brain matrices were harmonized to a shared ROI set. This included resolving minor ROI naming inconsistencies, combining predefined subdivisions, and collapsing slide-specific replicates by summing columns that shared the same base ROI and hemisphere. As an internal integrity check, hemisphere assignment for selected targets was evaluated by comparing mean projection strengths to left versus right targets stratified by soma hemisphere. In one instance when a consistent inversion was detected, the corresponding right hemisphere (RH) and left hemisphere (LH) target columns were swapped.

Per brain ipsilateral and contralateral matrices were then created by relabelling projection columns relative to soma hemisphere and retaining only common ROIs across RH and LH samples. Projection strengths were row-normalized by each cell’s total projection count, log1p-transformed, and scaled by the global maximum value for the given sample. Normalized matrices were concatenated across brains by row-binding. Metadata tables were combined and reordered to match the final matrix row order.

Projection-pattern analyses

The harmonized ipsi/contra-normalized projection matrix was visualized as heat maps. Cells were ordered by their top-projection region, defined as the region with the maximum normalized value per cell.

Regional module aggregation and coarse projection-pattern visualization

To align single-neuron MAPseq projection profiles with the coarse anatomical summaries used in ExA-SPIM, we mapped fine-grained target ROIs to 12 regions (OLF, CTX, HPF, CTXsp, CNU, TH, HY, MB, CB, P, MY and SP) and summed ipsi/contra-normalized projection strengths within each group to generate a cell × group matrix. Values were then row-normalized to fractions of each neuron’s total grouped projection strength (rows sum to 1) and visualized as a heat map after ordering neurons by their top (maximum-fraction) group to emphasize coarse projection motifs. To quantify coarse co-variation between region groups, we computed a group-by-group Spearman rank correlation matrix across neurons using the row-normalized grouped profiles (pairwise complete observations) and visualized the resulting 12 × 12 correlation matrix. For anatomical plots, per-cell summary annotations were generated from the grouped, row-normalized matrix: (1) the top-projection group label and (2) its corresponding maximum normalized value (‘top-projection strength’). The distribution of top-group assignments across cells was tabulated. Finally, top-group annotations were merged with BARseq-derived CCF coordinates from the combined metadata, and the merged table was exported for plotting.

snRNA-seq

LC dissection

Twenty male and 20 female C57BL/6J mice (age 57–69 days) were used for two rounds of LC dissection. Half the mice received viral injections that aimed to label LC-NE neurons with YFP. There was no significant labelling, and analyses were aggregated across batches of mice. Each mouse was anaesthetized with 2.5–3% isoflurane and transcardially perfused with oxygenated artificial cerebrospinal fluid (ACSF) ((in mM) calcium chloride, 0.5; d-glucose, 25; HCl, 98; HEPES, 20; magnesium sulfate, 10; monosodium phosphate, 1.25; myo-inositol, 3; N-acetylcysteine, 12; N-methyl-d-glucamine, 96; potassium chloride, 2.5; sodium bicarbonate, 25; sodium l-ascorbate, 5; sodium pyruvate, 3; taurine, 0.01; thiourea, 2; 300–310 mOsm). Brains were immediately dissected, mounted and cut into 350-μm coronal slabs in oxygenated ACSF with 5 μM anisomycin (A9789, Sigma) and 22.7 μM actinomycin-D (A1410, Sigma) on a Compresstome (Precisionary Instruments). LC with surrounding areas was dissected manually under a dissection microscope using local landmarks under bright field illumination. The dissected tissues were transferred to a 2-ml tube, immediately frozen in dry ice, and stored at –80 °C until all the samples in a cohort were collected.

Single-nucleus isolation

Nuclei were isolated using the RAISINs method with a few modifications as already described in a nuclei isolation protocol developed at Allen Institute for Brain Science. All reagents and methods are described in ref. ¹⁰⁹. In brief, dissectates were thawed and pooled separately in a 12-well plate containing CHAPS with salts and Tris (CST) extraction buffer (146 mM NaCl, 21 mM MgCl₂, 1 mM CaCl₂, 10 mM Tris-HCl pH 7.5, 0.49% (w/v) CHAPS, 0.1 mg ml^–1 BSA (Thermo Fisher), 1 U µl^–1 RNase inhibitor (2313B, Takara)). Tissues were chopped using spring scissors in ice-cold CST buffer for 10 min. The suspension was then transferred to a 50-ml conical tube while passing through a 100-μm filter and the walls of the tube were washed using salts and Tris (ST) buffer (146 mM NaCl, 21 mM MgCl₂, 10 mM Tris-HCl pH 7.5). Next, the suspension was gently transferred to a 15-ml conical tube and centrifuged in a swinging-bucket centrifuge for 5 min at 500 rcf and 4 °C. The supernatant was removed, and pellets were resuspended in 100 μl 0.1× lysis buffer and incubated for 2 min on ice. After addition of 1 ml wash buffer, samples were gently filtered using a 20-μm filter and centrifuged as above. After supernatant was discarded, pellets were resuspended in chilled nuclei buffer.

Sequencing library preparation

For 10x snRNA-seq, the cell suspensions were processed using Chromium GEM-X Single Cell 3′ Kit v4 (1000691, 10x Genomics). We followed the manufacturer’s instructions for cell capture, barcoding, reverse transcription, cDNA amplification and library construction¹¹⁰.

Sequencing data processing and QC

Processing of 10x Genomics snRNA-seq libraries was performed as described previously¹¹¹. We targeted a sequencing depth of 120,000 reads per cell. In brief, libraries were sequenced on an Illumina NovaSeq system at the Northwest Genomics Center at the University of Washington, and sequencing reads were aligned to the mouse reference transcriptome (M21, GRCm38.p6) using the 10x Genomics CellRanger pipeline (v.8.0.1) with the default parameters. The cell by gene matrix along with metadata was exported as R Data (.rda).

Data import, preprocessing and quality control

snRNA-seq data were processed in R using Seurat (Seurat v5.4.0; SeuratObject v5.3.0; RRID:SCR_016341) with batch correction performed using Harmony (RRID:SCR_022206). Raw UMI count matrices and accompanying metadata were used to generate a Seurat object. Nuclei with fewer than 2,000 or more than 15,000 detected genes were excluded, as were nuclei with >4% mitochondrial or >4% ribosomal transcript content. Genes with <3 cells detected were removed and cells with <200 detected genes were also excluded. Cells annotated with \(\ge \)0.4 doublet scores were removed. After quality control, 231,500 nuclei (from 398,912 initial profiles) were retained for downstream analyses. Data were count-per-10k normalized followed by log-normalizations. Then highly variable genes were identified, and scaled expression values were used for PCA. Batch effects were corrected by sequentially applying Harmony over vendor, batch, and RNA-amplification covariates, with each run initialized on the previous Harmony embedding (30 principal components). The integrated low-dimensional space was used to construct a shared-nearest-neighbour graph for community detection. Clusters were identified using graph-based clustering and visualized using UMAP computed on the Harmony-corrected embeddings.

Identification of noradrenergic cells

To classify noradrenergic cells, clusters with significant noradrenergic marker gene expression (including Dbh, Th, Slc6a2 and Slc18a2; Extended Data Fig. 6) were considered. Specifically, we calculated the average expression of each gene for each cluster, and only the clusters with larger than 1 unit of log CP10k expression for all the marker genes were kept. The remaining clusters were subset for further analysis. The LC subset was reprocessed independently, including normalization, identification of 2,000 variable features, scaling, PCA, and Harmony-based integration restricted to LC nuclei. A new neighbour graph, UMAP embedding, and graph-based clustering were computed to resolve LC subpopulations at higher resolution. Low-quality clusters were found by computing the median of the doublet score and total UMI counts, and clusters with the lowest median UMI (corresponding to the lowest detected genes) as well as the cluster with the highest median doublet score were excluded, and the final LC subcluster object was retained for downstream analyses, yielding 4,769 LC-NE cells. This subcluster was further examined for its quality, and we found 26 cells were outliers when reclustered (detected in UMAP), corresponding to a population with high doublet scores (average 0.173 versus 0.068 for the rest of the populations). We removed these 26 cells, yielding a final number of 4,743 LC-NE cells.

LC-NE transcriptomic analysis

Analyses were done in Python using the Scanpy package (RRID:SCR_018139). The selected cluster was saved as a separate AnnData object while retaining all genes. This dataset was used for downstream analyses. Gene chromosomal locations were queried using the mygene Python package (https://pypi.org/project/mygene/; RRID:SCR_018660). Genes located on chromosomes X or Y were excluded. From the remaining genes, 1,500 highly variable genes (HVGs) were identified (Seurat v3 implementation in Scanpy).

Batch correction and normalization with scVI

We ran batch correction first using scVI (v1.3.3, https://scvi-tools.org/; RRID:SCR_026673). Batches were defined as combinations of experimental date (two categories) and sex (two categories), resulting in four batches. Gene expression was modelled using a zero-inflated negative binomial (ZINB) distribution. A single-layer variational autoencoder with four latent dimensions was used to account for the relatively homogeneous cell population. All other parameters were set to default values. Normalized, batch-corrected expression values (ρ) were extracted from the trained model, representing reconstructed gene expression under an assumed library size of 1 (or converted to counts per million (CPM) when specified), with batch effects averaged across all batches. To calculate relative expression for genes in Extended Data Fig. 6g,h, many of which were not in the HVG set, we used the library-normalized counts of the dominant batch (3785 cells, 79.8% of the total snRNA-seq data).

Scaling and dimensionality reduction

Prior to PCA, gene expression values were z-score normalized by scaling each gene to zero mean and unit variance. Extreme standardized values (z-score normalization with scanpy function pp.scale) were clipped to ±10. PCA was performed on the 1,500 HVGs, and the top 50 principal components were retained. UMAP embeddings were computed from the scVI latent space using default Scanpy parameters for visualization.

Pseudospace analysis

To capture continuous transcriptional variation that we observed in PCA space, we developed a metric termed pseudospace, representing an ordering of cells along the primary transcriptional axis.

1. 1.

   Initial cell ordering. Cells were embedded in 2D UMAP space derived from PCA. The UMAP coordinates were centred, and singular value decomposition (SVD) was applied. The first principal component of the centred UMAP coordinates defines the primary axis of variation. Cells were projected onto this axis to obtain an initial ordering.

2. 2.

   Trajectory fitting in PCA space. Using the initial ordering, a smooth trajectory was fitted in 10-dimensional PCA space. Cells were ordered according to the initial ordering, and the first 10 principal components (explaining approximately 79.5% of variance) were extracted. For each principal component dimension, a cubic polynomial was fitted as a function of the ordering parameter using least squares. This yielded a parametric trajectory f(t) = [f₁(t), f₂(t),…, f₁₀(t)]. Here, subscripts indicate the principal component dimensions and t denotes the continuous pseudotime (ordering) parameter. The fitted function fk(t)(k = 1,…, 10) defines a smooth curve embedded in 10-dimensional space that approximates the underlying cell distribution.

3. 3.

   Pseudospace score calculations. The parametric trajectory was sampled at 10,000 evenly spaced points. Each cell was projected onto this trajectory by finding the value of t that minimized the squared Euclidean distance between the cell’s principal component coordinates x_i and the sampled points (on the fitted trajectory), \(_^=_\mathrm_Keep following us for the latest insights.-f(t)x^\). The closest trajectory point was identified, and the corresponding parameter value t* was assigned to the cell. These values were normalized to the range [0, 1] to produce final pseudospace scores.

Validating the robustness of pseudospace scores

Robustness was assessed using bootstrap analysis with 10,000 resamples. For each bootstrap sample, the trajectory was re-fitted and 95% confidence intervals were computed.

Identification of U-shaped genes

To identify genes exhibiting U-shaped (non-monotonic) expression patterns along the principal transcriptomic gradient we first fit the 1-dimensional pseudospace scores as described above, where each cell was assigned a pseudospace score. The gene expression profiles were reordered according to the pseudospace score, and the first and second derivatives of each gene’s smoothed expression with respect to pseudospace were computed. U-shaped genes were further selected based on the criterion that the first derivative at the beginning of the trajectory was negative (<−8) and at the end was positive (>8), indicating a decrease-then-increase (namely U-shaped) expression profiles. Candidate U-genes were ordered by their centres of mass along the pseudospace axis.

MMIDAS for identifying the potential clusters

MMIDAS, an unsupervised clustering framework based on coupled autoencoders with GAN-based augmentations³⁸, was applied as an independent validation approach. A two-arm variational autoencoder architecture was used with 10 latent dimensions, including 4 continuous dimensions and 15 categorical variables. All other parameters followed the implementation provided in the code repository.

Cluster quality analysis

Clustering quality was evaluated using silhouette scores across K = 1,…, 15 clusters. Analyses were performed on 10-dimensional PCA and scVI representations. Three types of data are used: (1) synthetic data with known ground truth (K = 2); (2) synthetic data generated from extreme pseudospace endpoints using a mixture of Gaussian distributions (with K = 2); (3) the LC-NE data. For the LC-NE data, we used two methods for the clustering: (a) Gaussian mixture models applied to PCA space; and (b) MMIDAS clustering evaluated in PCA space. For (1) and (2) data we used the same approach as (a). For each method and K, five independent runs were performed with different random seeds, and mean ± s.d. silhouette scores were reported. Higher silhouette scores (ranging from −1 to 1) indicate well-separated clusters, and they were always highest with a single cluster. As shown in Extended Data Fig. 6e, we observed a large decrease when k increased from 1 to 2 regardless of which method was used, indicating that one cluster best described the population.

Analysis of MERFISH data

Data source and image registrations

MERFISH datasets were downloaded from the original publication repositories³⁹. We analysed data from 4 mice (2 female, 2 male) from the original dataset, which gave 24 slices in total. The original data contain the two-dimensional coordinates and the gene counts for each cell. For each image slice, we ran Leiden clustering and coloured each cell based on the assigned clusters, before performing image registration manually using QuickNII (https://www.nitrc.org/projects/quicknii) for affine transformations and VisuAlign (https://www.nitrc.org/projects/visualign/; RRID:SCR_017978) for deformable transformations, to map to CCF coordinates. Transformation matrices were inverted and applied to images, which then yielded the three-dimensional spatial coordinates for each cell, which were aligned to anatomical reference frames.

Data preprocessing and LC-NE cell selection

MERFISH gene expression matrices and CCF-aligned cellular coordinates were processed independently for each tissue section. For each section, we first removed blank probes and then computed the product of raw counts for the canonical LC-NE markers Dbh, Th and Slc6a2 for every cell, retaining only cells with a marker product value >10. The filtered cells from all sections were concatenated into a single AnnData object with raw counts preserved and associated metadata (slicename and CCF coordinates). To further restrict the dataset to anatomically consistent LC neurons, cells with dorsal–ventral (CCF y axis) coordinates greater than 220 pixels (5.5 mm) were excluded. The resulting high-confidence LC-NE population was visualized in three-dimensional CCF space to confirm spatial consistency with the LC template.

Batch correction, dimensionality reduction and clustering

The LC-NE MERFISH AnnData object was used as input to scVI, with tissue slice section specified as the batch covariate. A variational autoencoder (one hidden layer, four latent dimensions) was trained under a ZINB likelihood to jointly perform library size normalization, denoising and batch correction. The resulting latent representations were used to construct a 30-nearest-neighbour graph, followed by UMAP embedding for visualization and Leiden community detection (resolution = 1) to identify transcriptionally distinct subpopulations. Leiden clusters were retained as the final LC-NE population when the mean expression of all three markers (Dbh, Th and Slc6a2) exceeded the cross-cluster average (z-score > 0), yielding 2,308 cells. Batch-corrected normalized expression values were obtained from the posterior mean for downstream interpretation. For all downstream analysis, library-normalized raw counts were used unless specified to use the batch-corrected counts.

CCA for space–gene relationships

Canonical correlation analysis (CCA) was performed between gene expression matrices and 3D spatial coordinates to identify space–gene relationships within the LC. Because the LC is approximately bilaterally symmetric, spatial coordinates were folded across the midline to remove left–right redundancy: coordinates from the right hemisphere were reflected across the midline and mapped onto the left hemisphere. Both gene expression and spatial features were then standardized (zero mean, unit variance) to remove scale differences across features. The top two canonical components were retained. The canonical variables (scalar) were then used as the cell weight for the visualization, termed weighted gene score.

Gene–space correlation analysis

Rank correlations were computed between individual gene expression levels and projection scores along identified spatial directions. Genes were ranked by absolute correlation values. The top 20 genes were listed in the table and visualized in CCF space.

Gene expression imputation

To expand gene coverage beyond the MERFISH panel and enable spatial visualization of genes not directly measured, we performed kNN imputation using matched snRNA-seq data as reference. Highly variable genes were identified in the snRNA-seq dataset, and the union of these genes with the MERFISH panel defined the target gene space for imputation. For similarity computation, both datasets were restricted to their shared genes. Expression profiles were normalized using a row-wise rank transformation followed by standardization to improve cross-platform comparability. For each MERFISH cell, we identified its k nearest neighbours (k = 200) in the snRNA-seq dataset using Euclidean distance in the normalized expression space. Gene expression for unmeasured genes was imputed as a weighted average of the corresponding snRNA-seq neighbour profiles across the full union gene set. Neighbour weights were derived from a distance-based softmax kernel, w_ij ∝ exp(−d_ij/τ), and normalized to sum to one for each MERFISH cell.

Choosing the number of neighbours for gene imputation

To determine the number of neighbours for the imputations, we took five genes that were included in the original MERFISH gene panel and we know a priori to have large variance, and imputed their expressions from the rest of the genes. We varied the number of neighbours and tested how many neighbours would give the optimal correlations between imputed and the ground truth expression values. Based on these results, we chose 200 neighbours for the following gene imputations.

Pseudospace score imputations

Similar to gene expression imputation, the pseudospace scores derived from snRNA-seq analysis were transferred to MERFISH cells (with k = 50). For each MERFISH cell, we computed the weighted average pseudospace score from k nearest snRNA-seq neighbours using the shared genes. To quantify uncertainty, we measured the weighted s.d. of the neighbours’ pseudospace scores, as well as the confidence intervals. For the confidence intervals, we used weighted sampling of 500 samples among the neighbours of each cell, based on the pre-calculated weights (which sum to one across all the neighbours) to get the resampled means. We then took the 2.5th and 97.5th percentiles. The width of the confidence interval (which contains 95% of the means, and the value is constrained between 0 and 1) was computed. The larger width indicates larger uncertainty.

To quantify mapping reliability, we estimated a baseline cross-modal distance distribution from randomly paired cells and compared each MERFISH cell’s mean neighbour distance to this baseline. We derived per-cell confidence metrics, including a normalized distance score and an empirical p-score (fraction of random distances exceeding the observed mean neighbour distance), which were retained for downstream quality control of imputed profiles.

Retro-seq

Experimental mice were generated by crossing homozygous RCFL-H2B-GFP mice¹¹² (Jackson Laboratory, 028581; RRID:IMSR_JAX:028581) to heterozygous Dbh-Cre mice (9 female, 10 male). The resulting offspring heterozygous for both Dbh-Cre and RCFL-H2B-GFP alleles were used for retrograde injections and subsequent SSv4 snRNA sequencing experiments following FACS sorting for GFP.

Dbh-Cre:RCFL-H2B-GFP mice (40–50 days old) were anaesthetized using oxygenated 4% isoflurane and maintained with oxygenated 0.7–1.0% isoflurane throughout the surgery. Animals were mounted in a stereotaxic frame (David Kopf Instruments), and small craniotomies were performed in the skull above the injection target (pressure injections using Nanoject III⁸⁸). Spinal cord injections were performed between C4–C5 vertebrae without drilling (Nanoject III⁹⁰). Intracerebral injections were always targeted to the right hemisphere. Spinal cord injections were performed bilaterally. Supplementary Table 2 lists injection targets, types, volumes and AP/ML/DV coordinates from Bregma. Virus AAV pEF1a-DIO-FLPo-WPRE-hGHpA was procured from Addgene (87306-AAVrg, lot v188462).

Tissue collection and processing

Mice were euthanized with an isoflurane overdose 4–5 weeks after AAV injection and transcardially perfused with cold, pH 7.4 HEPES buffer containing 110 mM NaCl, 10 mM HEPES, 25 mM glucose, 75 mM sucrose, 7.5 mM MgCl₂ and 2.5 mM KCl to remove blood from brain. Brains were quickly dissected out and flash frozen in liquid nitrogen. Frozen brains were placed into a 5-ml centrifuge tube containing 1 ml frozen OCT cushion at the bottom and capped.

Data preprocessing and quality control

Raw retro-seq counts were filtered to remove genes detected in fewer than 3 cells and cells with fewer than 200 genes. Cells from 20 donors with a known injection target (frontal cortex, cerebellum, spinal cord, or thalamus) were kept, and quality control retained cells with 2,000–15,000 detected genes, <10% mitochondrial counts, and <3% ribosomal counts. Non-neuronal contamination was removed by low-resolution clustering (Leiden, resolution = 0.1), keeping the cluster with the highest mean expression of the marker genes (Dbh, Th, Slc6a2 and Slc18a2). A second-pass threshold (>4,000 detected genes) yielded 1,076 high-quality LC-NE neurons. Retro-seq data were subset to genes shared with the snRNA-seq reference. Thalamic projecting neurons were kept but the projection information was excluded from downstream analysis to avoid potential confusion with the frontal cortex-projecting neurons.

Gene length normalization by TPM normalization

SmartSeq2 technology is known to exhibit gene length-dependent bias; therefore, expression values for all samples were normalized by the corresponding gene lengths. Gene lengths were derived from a reference GTF annotation file from GENCODE (vM38, https://www.gencodegenes.org/mouse/; RRID:SCR_014966) by summing exon lengths per gene. Exons lacking a ’gene_name’ annotation in the GTF file were excluded. For each gene, exon lengths were summed across all annotated exons to obtain the total exonic gene length in base pairs.

Gene length-normalized expression values were computed by dividing raw transcript counts by the corresponding gene lengths. Counts for each sample were scaled to CPM, then rounded to integers, and used for all downstream analyses.

Assessment of batch effects using supervised classification

Residual batch effects were assessed using supervised classification. Following TPM normalization, log transformation, scaling, and PCA, random forest classifiers were trained to predict donor identity, sex, or projection target from gene expression profiles. Model performance was evaluated using fivefold stratified cross-validation and summarized as mean accuracy ± s.d. Classification accuracy was compared to the theoretical chance level (1/K where K is the number of classes), with above-chance performance indicating detectable batch or covariate structure in the data.

Projection site predictions from gene expression

Projection sites were predicted using supervised classification. The dataset was split into training (60%), validation (20%) and test (20%) sets using stratified sampling to preserve class proportions. Five algorithms were evaluated: logistic regression, random forest, support vector machine, kNN and gradient boosting. Hyperparameter optimization was performed using fivefold cross-validation on the training set. Models were first evaluated on the validation set, then retrained on the combined training and validation set, and finally assessed on the test set. Classification performance was quantified using accuracy.

Projection site imputations for snRNA-seq data

As in the MERFISH analysis, projection site labels were transferred from retro-seq data to snRNA-seq data using kNN based on rank-based correlations. For each snRNA-seq cell, the 10 nearest neighbours in the retro-seq dataset were identified. A projection target was assigned if at least 8 of the 10 neighbours agreed on the projection site. The quality of imputation was computed as disagreement score, which is 1 − (number of agreement)/(number of valid projections). (Invalid projections indicate thalamus-projecting neurons, which is not considered in this case).

Patch-seq

Transgenic mice and labelling projection neurons

Heterozygous Dbh-Cre mice crossed with homozygous Ai65 mice (P28–P30) were anaesthetized with 2.0% isoflurane in O₂ and maintained at 0.7–1.0% isoflurane throughout surgery. Mice were mounted in a stereotaxic frame (David Kopf Instruments) for injections of AAV pEF1a-DIO-FLPo-WPRE-hGHpA (Addgene 87306-AAVrg, lot v188462). Small craniotomies were made in the skull above the injection target for intracerebral injections⁸⁹. Spinal cord injections were made between C4–C5 vertebrae without drilling⁹⁰. Intracerebral injections were in the right hemisphere. Spinal cord injections were bilateral. Mice recovered in their home cage for approximately 21 days before Patch-seq experiments.

Tissue processing

For preparation of acute brain slices, adult male and female mice (50–74 days old) were first fully anaesthetized by 5% isoflurane inhalation. A transcardial perfusion was then performed with ~10 ml of ice-cold cutting artificial cerebrospinal fluid (ACSF; 0.5 mM calcium chloride (dihydrate), 25 mM d-glucose, 20 mM HEPES buffer, 10 mM magnesium sulfate, 1.25 mM sodium phosphate monobasic monohydrate, 3 mM myo-inositol, 12 mM N-acetyl-l-cysteine, 96 mM N-methyl-d-glucamine chloride (NMDG-Cl), 2.5 mM potassium chloride, 25 mM sodium bicarbonate, 5 mM sodium l-ascorbate, 3 mM sodium pyruvate, 0.01 mM taurine, and 2 mM thiourea (pH 7.3), which had been continuously bubbling with a mixture of 95% O₂/5% CO₂). Coronal or sagittal sections (thickness 350 μm) containing the LC were sliced on a vibrating microtome (VT1200S Vibratome, Leica Biosystems). Immediately after slicing, brain slices were placed in warm (34 °C) oxygenated cutting ACSF for 10 min, then allowed to further recover in holding ACSF (2 mM calcium chloride (dihydrate), 25 mM d-glucose, 20 mM HEPES buffer, 2 mM magnesium sulfate, 1.25 mM sodium phosphate monobasic monohydrate, 3 mM myo-inositol, 12.3 mM N-acetyl-l-cysteine, 84 mM sodium chloride, 2.5 mM potassium chloride, 25 mM sodium bicarbonate, 5 mM sodium l-ascorbate, 3 mM sodium pyruvate, 0.01 mM taurine, and 2 mM thiourea (pH 7.3)), bubbling with a mixture of 95% O₂/5% CO₂ at room temperature until transferred to the microscope for recordings.

Patch-seq recordings

Slices were bathed in warm (34 °C) recording ACSF (2 mM calcium chloride (dihydrate), 12.5 mM d-glucose, 1 mM magnesium sulfate, 1.25 mM sodium phosphate monobasic monohydrate, 2.5 mM potassium chloride, 26 mM sodium bicarbonate, and 126 mM sodium chloride (pH 7.3)) and continuously bubbled with 95% O₂/5% CO₂. The bath solution contained blockers of fast glutamatergic (1 mM kynurenic acid) and GABAergic (γ-aminobutyric acid-dependent) synaptic transmission (0.1 mM picrotoxin). Thick-walled borosilicate glass (Warner Instruments, G150F-3) electrodes were manufactured (Narishige PC-10) with a resistance of 4–5 MΩ. Before recording, the electrodes were filled with approximately 1.0 to 1.5 μl of internal solution with biocytin (110 mM potassium gluconate, 10.0 mM HEPES, 0.2 mM ethylene glycol-bis (2-aminoethylether)-N,N,N,N-tetraacetic acid, 4 mM potassium chloride, 0.3 mM guanosine 5-triphosphate sodium salt hydrate, 10 mM phosphocreatine disodium salt hydrate, 1 mM adenosine 5-triphosphate magnesium salt, 20 µg ml⁻¹ glycogen, 0.5 U RNAse inhibitor (Takara, 2313 A) and 0.5% biocytin (Sigma B4261), pH 7.3). The pipette was mounted on a Multiclamp 700B amplifier headstage (Molecular Devices) fixed to a micromanipulator (PatchStar, Scientifica). Electrophysiology signals were recorded using an ITC-18 Data Acquisition Interface (HEKA). Commands were generated, signals processed, and amplifier metadata were acquired using MIES (https://github.com/AllenInstitute/MIES; RRID:SCR_016443), written in Igor Pro (Wavemetrics; RRID:SCR_000325). Data were filtered (Bessel) at 10 kHz and digitized at 50 kHz. Data were reported uncorrected for the measured liquid junction potential of −14 mV between the electrode and bath solutions. Prior to data collection, all surfaces, equipment and materials were thoroughly cleaned in the following manner: a wipe down with DNA away (Thermo Scientific), RNAse Zap (Sigma-Aldrich), and finally with nuclease-free water. After formation of a stable seal and break-in, the resting membrane potential of the neuron was recorded (typically within the first minute). A bias current was injected, either manually or automatically using algorithms within the MIES data acquisition package, for the remainder of the experiment to maintain that initial resting membrane potential. Bias currents remained stable for a minimum of 1 s before each stimulus current injection. To be included in the analysis, neurons needed to have a >1 GΩ seal recorded before break-in and an initial access resistance <20 MΩ and <15% of the R_input. For an individual sweep to be included, the following criteria were applied: (1) the bridge balance was <20 MΩ and <15% of R_input; (2) bias (leak) current within ±100 pA; and (3) root mean square noise measurements in a short window (1.5 ms, to gauge high frequency noise) and longer window (500 ms, to measure patch instability) <0.07 mV and <0.5 mV, respectively. After electrophysiological recording, the pipette was centred on the soma or placed near the nucleus (if visible). A small amount of negative pressure was applied (approximately −30 mbar) to begin cytosol extraction and to attract the nucleus to the tip of pipette. After approximately 1 min, the soma visibly shrank and/or the nucleus was near the tip of the pipette. While maintaining negative pressure, the pipette was slowly retracted; slow, continuous movement was maintained while monitoring the pipette seal. Once the pipette seal reached >1 GΩ and the nucleus was visible on the tip of the pipette, the speed was increased to remove the pipette from the slice. The pipette containing internal solution, cytosol, and the nucleus was removed from pipette holder, and its contents were expelled into a PCR tube containing the lysis buffer (Takara, 634894). Electrophysiological features were computed using the Intrinsic Physiology Feature Extractor (IPFX) Python package.

cDNA and sequencing

For Patch-seq experiments, we reverse transcribed the collected nuclear and cytosolic mRNA, and sequenced the resulting cDNA using a SMART-Seq v4 method^41,48. We used the SMART-Seq v4 Ultra Low Input RNA Kit for Sequencing (Takara, 634894) to reverse transcribe poly(A) RNA and amplify full-length cDNA according to the manufacturer’s instructions. We performed reverse transcription and cDNA amplification for 21 PCR cycles in 0.65 ml tubes, in sets of 88 tubes at a time. At least 1 control 8-strip was used per amplification set, which contained 4 wells without cells and 4 wells with 10 pg control RNA. Control RNA was either Mouse Whole Brain Total RNA (Zyagen, MR−201) or control RNA provided in the SMART-Seq v4 kit. All samples proceeded through Nextera XT DNA Library Preparation (Illumina FC-131-1096) using either Nextera XT Index Kit V2 Set A–D (FC-131-2001,2002,2003,2004) or custom dual-indexes provided by IDT (Integrated DNA Technologies). Nextera XT DNA Library prep was performed according to manufacturer’s instructions except that the volumes of all reagents including cDNA input were decreased either to 0.4× or to 0.2× by volume. Each sample was sequenced to approximately 500,000 – 1 million reads.

Fifty-base-pair paired-end reads were aligned to mm10 GENCODE vM23/Ensembl 98 reference genome, downloaded from 10X cell ranger (refdata-cellranger-arcmm10-2020-A-2.0.0). Sequence alignment was performed using STAR aligner (v2.7.1a; RRID:SCR_004463) with default settings. PCR duplicates were masked and removed using STAR option bamRemoveDuplicates. Only uniquely aligned reads were used for gene quantification. Gene counts were computed using the R Genomic Alignments package¹¹³ summarizeOverlaps function using IntersectionNotEmpty mode for exonic and intronic regions separately. Exonic and intronic reads were added together to calculate total gene counts; this was done for both the reference dissociated cell dataset and the Patch-seq dataset. Data were analysed as CPM.

Slices from Patch-seq experiments were mounted on slides and imaged on an upright AxioImager Z2 microscope (Zeiss, Germany) with an Axiocam 506 monochrome camera and 0.63× optivar lens. Two-dimensional tiled overview images were also captured (Zeiss Plan-NEOFLUAR 20×/0.5) in brightfield transmission and fluorescence channels. Light was transmitted using an oil-immersion condenser (1.4 NA). High-resolution, multi-tile image stacks were captured (Zeiss Plan-Apochromat 63×/1.4 Oil or Zeiss LD LCI Plan-Apochromat 63×/1.2 Imm Corr) at an interval of 0.28 μm (1.4 NA objective) or 0.44 μm (1.2 NA objective) along the Z axis. Image tiles were stitched in ZEN software and exported as single-plane TIFF files.

Physiology and behaviour

Mice and surgery

We used heterozygous Dbh-Cre mice for all in vivo experiments. Forty-two mice (6 female) were used for behavioural experiments, 28 mice (5 female) were used for electrophysiological recordings during the behavioural task, and 14 mice (1 female) for calcium indicator recordings. Surgery was performed on mice between the ages of P56–P112, under isoflurane anaesthesia (1.5–2.0% in O₂) and under aseptic conditions. During all surgeries, titanium headplates were surgically attached to the skull using dental adhesive (C&B-Metabond, Parkell). After the surgeries, analgesia (ketoprofen, 5 mg kg⁻¹ and buprenorphine, 0.05–0.1 mg kg⁻¹) was administered to minimize pain and aid recovery. Surgeries performed at the Allen Institute followed the described protocol¹¹⁴.

For all experiments, mice were given at least one week to recover prior to water restriction. During water restriction, mice had free access to food and were monitored daily in order to maintain 80% of their baseline body weight. Mice were housed in a reverse light cycle (12 h dark:12 h light, dark from 08:00–20:00) and all experiments were conducted during the dark cycle between 12:00 and 20:00.

Behavioural task

Before training on the tasks, water-restricted mice were habituated to head fixation for 1–3 days with free access to water from the two spouts (21 gauge stainless steel tubes separated by 4 mm) placed in front of the 38.1-mm diameter acrylic tube in which the mice rested. For LC-NE recordings with tetrodes and calcium indicator recordings in LC axons, the spouts were mounted on a micromanipulator (DT12XYZ, Thorlabs) with a custom digital rotary encoder to measure the position of the lick spouts with 5–10 μm resolution. Each spout was attached to a solenoid (ROB-11015, Sparkfun) to enable movement in the anterior–posterior axis of the mouse. The tones used for the cues (randomly assigned to go and no-go cues per mouse) were 7.5 and 15 kHz square waves generated by microcontrollers (ATmega16U2 or ATmega328), amplified and delivered through speakers (CUI Devices, GF0401M).

Licks were detected using custom circuits (Janelia Research Campus 2019-053). Task events were controlled and recorded using custom code (Arduino or Bonsai: https://github.com/AllenNeuralDynamics/dynamic-foraging-task) written for microcontrollers (ATmega16U2, ATmega328, or Harp boards). Water rewards were 2–3 μl, adjusted for each mouse to maximize the number of trials completed per session and to keep sessions around 60 min. Solenoids (LHDA1233115H or LHDB1233518H, The Lee Co) were calibrated to release the desired volume of water and were mounted on the outside of the sound-attenuated chamber used for behaviour. For mice with tetrodes, white noise (2–60 kHz, Sweetwater Lynx L22 sound card, Rotel RB-930AX two-channel power amplifier, and Pettersson L60 Ultrasound Speaker), was played inside the chamber to mask ambient noise.

For LC-NE recordings with Neuropixels probes, behaviour was controlled by Harp devices using Python and Bonsai (https://github.com/AllenNeuralDynamics/dynamic-foraging-task). Lick spouts were mounted on a New Scale microcontroller (M3-LS-3.4-15) to track and control lick spout locations across days. The tone used for the cue was 7.5 kHz at 70 dB, generated by a Harp Soundcard. Licks were detected using a custom Harp device (Lickety-split, https://github.com/AllenNeuralDynamics/harp.device.lickety-split).

During the 1–3 days of habituation, mice were trained to lick both spouts to receive water. Water was delivered following a lick to the correct spout at any time. Reward probabilities were chosen from the set and reversed every 20 trials. In the second stage of training (5–12 days), the trial structure with tone presentation was introduced. Each trial began with the 0.5 s delivery of either an auditory go cue (P = 0.95) or a no-go cue (P = 0.05) in tetrode and GCaMP experiments. Following the go cue, mice could lick either the left or the right spout. If a lick was made during a 1.5 or 1.8 s response window (in separate experiments), reward was delivered probabilistically from the chosen spout. In tetrode and GCaMP experiments, the unchosen spout was retracted when the tongue contacted the chosen spout to prevent mice from sampling both spouts within a trial. The unchosen spout was replaced 3.1 s after cue onset (or 3.0 s in experiments with axonal GCaMP). Following a no-go cue, lick responses were neither rewarded nor punished. Reward probabilities during this stage were chosen from the set and reversed every 20–35 trials. During this period of training only, water was occasionally manually delivered to encourage learning of the response window and appropriate switching behaviour. During this second stage of training, we introduced a 100-ms delay between choice and outcome. This delay was gradually increased to 300 ms for electrophysiological experiments and 200 ms for axonal GCaMP. If a directional lick bias was observed in one session, the lick spouts were moved horizontally 50–300 μm in the opposite direction prior to the following session.

After the 3.0 or 3.1 s trial duration, inter-trial intervals (the times between consecutive cue onsets) were generated as draws from an exponential distribution with a rate parameter of 0.3 and a maximum of 20 s. This distribution results in a flat hazard rate for inter-trial intervals such that the probability of the next trial did not increase over the duration of the inter-trial interval. Inter-trial intervals were 4.83 s on average (range, 1.1–20 s, including a no-lick window; see below). In more than 90% of all well-trained sessions, mice made a leftward or rightward choice in greater than 90% of trials.

In the final stage of the task, the reward probabilities assigned to each lick spout were drawn pseudorandomly from the set . The probabilities were assigned to each spout individually with block lengths drawn from a uniform distribution of 20–35 trials. To stagger the blocks of probability assignment for each spout, the block length for one spout in the first block of each session was drawn from a uniform distribution of 6–21 trials. For each spout, probability assignments could not be repeated across consecutive blocks. To maintain task engagement, reward probabilities of 0.1 could not be simultaneously assigned to both spouts. If one spout was assigned a reward probability greater than or equal to the reward probability of the other spout for 3 consecutive blocks, the probability of that spout was set to 0.1 to encourage switching behaviour and limit the creation of a direction bias. If a mouse perseverated on a spout with reward probability of 0.1 for 4 consecutive trials, 4 trials were added to the length of both blocks. This prevented mice from choosing one spout until the reward probability became high again.

To minimize spontaneous licking, we enforced a 1 s no-lick window prior to tone onset. Licks within this window were punished with a new randomly generated inter-trial interval, followed by a 1 s no-lick window. Implementing this window significantly reduced spontaneous licking throughout the entirety of behavioural experiments.

Electrophysiology with viral injections

To express ChR2¹¹⁵, Chrimson¹¹⁶, or ChRmine¹¹⁷ in LC-NE neurons, we pressure-injected rAAV5-EF1a-DIO-hChR2(H134R)-EYFP (3 × 10¹³ GC ml⁻¹), AAV5-Syn-FLEX-rc[ChrimsonR-tdTomato] (2.2 × 10¹³ GC ml⁻¹) or pAAV-Ef1a-DIO-ChRmine-eYFP-WPRE (9.69 × 10¹² GC ml⁻¹, Addgene; packaged in house, lot number VT8453G) into the LC of Dbh-Cre mice at a rate of 1 nl s⁻¹ (MMO-220A, Narishige). pAAV-EF1a-double floxed-hChR2(H134R)-EYFP-WPRE-HGHpA was a gift from K. Deisseroth (Addgene viral prep 20298-AAV5; RRID:Addgene_20298). pAAV-Syn-FLEX-rc[ChrimsonR-tdTomato] was a gift from E. Boyden (Addgene plasmid 62723; http://n2t.net/addgene:62723; RRID:Addgene_62723). pAAV-Ef1a-DIO ChRmine-eYFP-WPRE was a gift from K. Deisseroth (Addgene plasmid 130996; http://n2t.net/addgene:130996; RRID:Addgene 130996).

For ChR2, we made three injections of 200 nl at the following coordinates: 0.28–0.35 mm anterior to the junction of the inferior colliculus and cerebellum, 0.85–0.90 mm right from the midline, and Check back often for more exciting news! mm ventral from the pial surface. Before the first injection, the pipette was left at the most ventral coordinate for 5 min. Before each injection, the pipette was withdrawn 50 μm and left in place for 5 min after injection. For electrophysiology experiments with rAAV5-EF1a-DIO-hChR2(H134R)-EYFP injections, the microdrive was implanted through the same craniotomy.

For Chrimson and ChRmine, we injected 300 nl at the following coordinates: from bregma, 5.4 mm posterior for male mice, 5.2 mm posterior for female mice; 0.85 mm lateral; and (2.9, 3.1) ventral from surface as described¹¹⁴.

Electrode recording

We used two types of electrodes for extracellular recordings: tetrodes and Neuropixels probes. For electrophysiological experiments with tetrodes, we implanted a custom microdrive targeting right LC, entering through a craniotomy at 0.28–0.35 mm anterior to the junction of the inferior colliculus and cerebellum and 0.85–0.90 mm lateral from the midline (identified using the vasculature as a landmark). For recordings with tetrodes, we sampled at 32 kHz (Digital Lynx 4SX, Neuralynx, Inc.). The recording system was connected to 8 implanted tetrodes (32 channels, nichrome wire, PX000004, Sandvik) fed through guide tubes that could be advanced with the turn of a screw on a custom, 3D-printed microdrive. The impedances of each wire in the tetrodes were reduced to 180–250 kΩ by gold plating. The tetrodes were adhered to a 200-μm optic fibre used for optogenetic identification. After each recording session, the tetrode–optic fibre bundle was driven down 30–75 μm. Channels were bandpass-filtered between 0.3–6 kHz. The bandpass-filtered signal, x, was thresholded at 2.5 σ_n, where \(_=\mathrm\left(\frac\right)\). Detected peaks were sorted into individual unit clusters offline (Spikesort 3D, Neuralynx) using peak waveform amplitude, minimum waveform trough, and waveform PCA. We used two metrics of isolation quality as inclusion criteria: L-ratio (<0.05) and interspike interval (ISI) violation ratio <0.1%.

For electrophysiological experiments with Neuropixels probes, we made acute craniotomies at 0.28–0.35 mm anterior to the junction of the inferior colliculus and cerebellum and 0.85–0.90 mm lateral from the midline prior to recording as described¹¹⁸ and inserted the probe either vertically or with a 6-degree anterior angle. A Kilosort (RRID:SCR_016422) pipeline¹¹⁹ was used for spike detection and sorting into single units. Signals adjacent to laser stimulation (0–4 ms) were removed and smoothed. Signals were filtered between 0.3–6 kHz, with phase shift correction and common reference removal. Probe displacement was estimated using the DREDge algorithm¹²⁰. For some sessions, we also made a second recording at the brain surface without moving the probe to estimate the brain’s surface location along the probe.

The following were applied as initial inclusion criteria for single units unless further specified: ISI violation ratio <0.1 (unless identified to project to cortex, where we used 0.2 as a threshold); presence ratio >0.8, amplitude cut-off <0.1. Unit drift detection was performed by checking for timepoints where abrupt spike rate changes could be predicted by changes of probe location (estimated by DREDge) and changes of waveform peaks. Only periods without abrupt changes were used for future analysis. See details in ‘Drift control’.

Spike waveform extraction

To minimize distortion of extracellular spike waveforms by filtering, raw traces were Butterworth filtered between 50–8,000 Hz before computing spike waveforms. In waveform analysis, only neurons recorded with tetrodes or Neuropixels 2.0 were included. To optimize the data for unit waveform characterization, unit inclusion criteria were modified as follows: ISI violation ratio bar was raised to 0.5, and a peak threshold of lower than −50 μV was applied. Each spike waveform was aligned to the time of its peak and normalized so that baseline started at zero and had peak values of 1. Six features were used to characterize each waveform: post-peak trough, post-peak time, post-slope, pre-slope, slope symmetry, trough distance, and trough symmetry (Extended Data Fig. 11a). All features were z-scored before further analysis.

Optical identification of LC-NE neurons

In tetrode recordings, light was delivered through a fibre attached to the tetrode at 5 or 10 Hz with pulse width of 10 or 20 ms. Laser irradiance was 0.2 mW, measured from the patch cord. In Neuropixels 2.0 recordings, light was delivered from brain surface above the inferior colliculus at 5 Hz with pulse width 4 or 5 ms. Laser irradiance was 10, 20, 30, 40 or 50 mW, measured at the collimator that focused the laser onto the brain surface.

Individual neurons were determined to be identified LC-NE neurons using the following criteria: (1) responses to brief pulses (5, 10 or 20 ms) of laser stimulation (473, 560 or 640 nm wavelength, depending on the opsin) with significant increase of firing rate compared to baseline on both early and late pulses in a train, with short latency (<20 ms but >2 ms) and with higher probability than sham (>0.55; Extended Data Fig. 10f). (2) Light-evoked spike waveforms were similar to spontaneous waveforms, with correlation coefficient greater than 0.95 and Euclidean distance (normalized to spike amplitude) less than 0.35.

In addition to observing responses to light stimuli, in tetrode experiments, LC targeting was confirmed by performing electrolytic lesions of the tissue (25 s of 20 μA direct current across two wires of the same tetrode) and examining the tissue after perfusion; in Neuropixels recordings, probe tracks were reconstructed (details below).

Antidromic stimulation

To identify if an LC-NE neuron projected to cortex, light was delivered above multiple sites in frontal cortex where the skull was thinned and cleared using cyanoacrylate, covered by Kwik-Cast or Kwik-Sil to avoid damage. Light (640 nm) was delivered at 5 Hz with a pulse width of 4 ms.

For each neuron, we defined antidromic response time (τ_anti) as the peak time of the peristimulus time histogram (PSTH) aligned to antidromic stimulation; response jitter as the half-peak-width of the PSTH around τ_anti; collision window as the time window τ_anti before and τ_anti after antidromic stimulation (Extended Data Fig. 10j). The collision window is the time during which if a spontaneous spike occurred, it would collide with the backpropagating antidromic spike.

We fit a regression to the spike rate λ_resp around τ_anti after real or sham (equal number of random timepoints as real) antidromic stimulations. Regressors were as follows: sham or real laser stimulation \(L\in (0,1)\) (0 for sham), spike rate in the collision window λ_col, interaction between the presence of laser stimulation and presence of spikes in the collision window, such that λ_resp ∼ 1 + L + λ_col + L × λ_col. A neuron was identified as projecting to cortex if it responded to laser stimulation consistently (t-statistic of laser stimulation L > 0 and P < 0.005) unless there were spontaneous spikes in the collision window (t-statistics of interaction between laser stimulation L and spike rate in the collision window λ_col less than −3.5 and P < 0.005) with low response jitter (<20 ms). Jitter threshold was applied here to ensure accuracy of τ_anti estimation.

Histology

After experiments were completed, mice were euthanized with an overdose of isoflurane, exsanguinated with saline, and perfused with 4% paraformaldehyde. For mice with tetrodes, the brains were cut in 100-μm-thick coronal sections and mounted on glass slides. We validated expression of rAAV5-EF1a-DIO-hChR2(H134R)-EYFP with epifluorescence images of LC (Zeiss Axio Zoom.V16) and confirmed targeting of the optic fibre–tetrode bundle to LC by colocalization of the electrolytic lesion with immunostaining against TH (rabbit anti-TH, ab152, Millipore, RRID:AB_390204, 1:500, followed by donkey anti-rabbit Alexa 488, Life Technologies, RRID:AB_2556546, 1:1,000, or goat anti-rabbit Alexa 594, Life Technologies, RRID:AB_2534079, 1:1,000). In electrophysiological experiments with Neuropixels probes, probe tracks were reconstructed as described below.

Probe reconstruction

In Neuropixels recordings, probes were dipped into fluorescent dye before penetrating the brain. After 1 week of recording (3–5 recording sessions), mice were perfused with PFA, brains were extracted and cleared with a LifeCanvas protocol, and imaged using SPIM. Images were processed through a SmartSPIM pipeline (https://github.com/AllenNeuralDynamics/aind-smartspim-pipeline), where images were stitched and registered to CCF.

In stitched images, probe tracks were manually labelled and transformed into CCF coordinates. To further refine the location of the probe¹²¹, electrophysiological features along the probe were aligned to brain regions along the probe track using a modified version of International Brain Laboratory software (https://github.com/AllenNeuralDynamics/ibl-ephys-alignment-gui). Key features for identifying the locations of recorded cells dorsal to LC include the discontinuity of spatial cross-correlation of local field potentials at the dorsal edge of the pons, change of signal amplitude among cerebellum layers, and transient signal change at the brain surface.

Drift control

In recordings with Neuropixels probes, spike detection and clustering are subject to movement artefacts when the brain moves relative to the probe. We controlled for both abrupt movement and slow drift of the probe.

To remove the effects of abrupt displacement, we restricted our analyses to time windows without such events. To identify when abrupt changes occurred, we tested whether sudden changes in spike rate could be predicted by concurrent changes in spike waveform amplitude and probe position. Specifically, we estimated the first derivative of spike rate λ′(t), probe displacement d′(t) and spike waveform peak amplitude p′(t) from one channel for each neuron, at each time point t. For a generic signal v(t) at time t₀, the first derivative was estimated using a windowed contrast measure defined as follows:

where m(v, t₀, [a, b]) denoted the windowed mean of signal v in the interval [t₀ + a, t₀ + b), defined as

where \(N=|\{i:{t}_{i}\in [{t}_{0}+a,{t}_{0}+b)\}|\). The first derivatives were estimated every 100 s. A regression model was fitted to predict |λ′(t)| using |p′(t)| and |d′(t)| as, |λ′(t)| ∼ 1 + |d′| + |p′|. The estimated absolute value of the first derivative of firing rate is \(|\hat{\lambda \text{‘}}|\). This analysis was performed at two different timescales, short and long, with pre = −300 s for both timescales, post = 100 s for short timescales, and post = 300 s for long timescale calculations. Abrupt drift events were detected at time t when both λ′(t) and \(|\hat{\lambda \text{‘}}|\) exceeded 0.5, with either short or long timescale calculations.

We also estimated slow drifts of the probe over long timescales. A neuron’s overall relative spike rate change over a whole session was calculated as the ratio of standard deviation of spike rate divided by mean spike rate, calculated by binning the session into non-overlapping 300-s bins. Noradrenergic neurons with this ratio higher than 0.3 (98th percentile) after being restricted to periods without abrupt drift events were excluded from analysis.

To avoid potential effects of slow electrode motion on spike sorting, one drift regressor was added to all single-neuron analyses with linear regression. The drift regressor was the deviation of mean waveform amplitude on its peak channel for each trial from the mode of waveform amplitudes.

Axonal calcium indicator measurements

To express GCaMP6s, or jGCaMP8s⁴⁴ in LC-NE neurons, we pressure-injected within LC AAV-hSyn-FLEX-axon-GCaMP6s (1.0 × 10¹³ GC ml⁻¹) or AAV-hSyn-FLEX-Axon-jGCaMP8s (1.0 × 10¹³ GC ml⁻¹) into the LC of Dbh-Cre mice at a rate of 1 nl s⁻¹. pAAV-hSynapsin1-FLEx-axon-GCaMP6s was a gift from L. Tian (Addgene viral prep 112010-AAV5; RRID:Addgene_112010)¹²². For Axon-jGCaMP8s, based on the Axon-GCaMP6s sequence (Addgene: 112010), we PCR-amplified AvrII-Axon-jGCAMP8s-NheI. The construct was then ligated into the hSynapsin1-FLEX pAAV backbone cut from pAAV-Syn-FLEX-rc[ChrimsonR-tdTomato] (Addgene: 62723) with AvrII/NheI.

Viruses were injected at following locations: from bregma, 5.4 mm posterior for male mice (400 nl), 5.2 mm posterior for female mice (300 nl); 0.85 mm bilateral; and (2.9, 3.1) mm from surface.

Fibre implantation

We implanted 200-μm-diameter optic fibres (0.39 NA) above prelimbic cortex (2.0 mm anterior to bregma, 0.5 mm lateral, 1.5 mm ventral from the pial surface) as described in detail¹²³. We measured GCaMP6s and GCaMP8s fluorescence using 60 μW irradiance of 488-nm (to excite GCaMP6s and GCaMP8s) and 40 μW 405-nm (as an isosbestic control) light using a custom photometry system¹²⁴. Signals were sampled at 60 Hz, alternative sampling from 405 nm, 488 nm and 560 nm (0 irradiance from 560 nm in these experiments) channels, resulting in 20 Hz for each channel. Example raw traces are in Extended Data Fig. 14c.

Data preprocessing

Raw photometry traces from both channels (488 nm and 405 nm) were first de-trended by fitting and removing a tri-exponential curve, each modelling signal decay with different timescales. ΔF/F was computed as the ratio of the de-trended signal and the trend. After this, we generated a low-pass filtered baseline at 0.05 Hz for both channels’ ΔF/F and fitted the isosbestic channel’s (405 nm) baseline to the 488 nm channel’s baseline. The fitted baseline was treated as a movement artefact estimated by the isosbestic channel. This baseline was removed from ΔF/F from the 488 nm channel, generating a motion-corrected de-trended signal. This signal was then low-pass filtered at 3 Hz before further analysis. Code can be found at https://github.com/AllenNeuralDynamics/aind-fip-dff.

Comparison between single-neuron activity and axonal GCaMP

To test whether LC-NE neurons projecting to cortex differed from the whole distribution of neurons, we compared task responses of single LC-NE neurons identified as projecting to cortex versus those that were not. To compare t-statistics from regression models between these two groups of neurons, we used Welch’s two-sample t-test (two-sided), which allows unequal variances. In addition to the asymptotic P value, significance was also assessed using a label-permutation test in which pooled observations were randomly reassigned to groups while keeping group sizes. Effect size was quantified using Cohen’s d, and 95% confidence interval was estimated from bootstrapping by resampling each group with replacement across bootstrap iterations.

Comparisons between single-neuron spiking and photometry-derived measures, are sensitive to modality-dependent differences in effect magnitudes or signal-to-noise ratios. Thus, we used a sign-based comparison to compare single-neuron spiking to photometry. To compare the sign of effects (t-statistics from regression models), we computed the proportion of positive values in each group and tested for differences between groups. Group differences were quantified using the risk difference, risk ratio, and odds ratio. Statistical significance reported in the main text was evaluated using Fisher’s exact test. A two-proportion z-test and a label-permutation test on the difference in positive rates obtained by shuffling group labels while keeping group sizes were additionally performed and are reported in the accompanying analysis code.

Circular distribution comparison

When comparing differences between outcome-related activity of neurons projecting to cortex versus not and neuronal spiking versus calcium indicator dynamics in PL, differences between circular distributions were assessed using a permutation-based implementation of the Mardia–Watson–Wheeler test. Angular observations were first converted to floating-point values. All angles were then wrapped to the interval [0, 2π) using a modulo 2π transformation to ensure a common circular reference frame.

For two samples containing n₁ and n₂ observations, respectively, angles from both groups were pooled and assigned ranks across the combined dataset using average ranks for tied values. The pooled ranks r_i were then mapped onto the circle as ϕ_i = 2πr_i/n, where n = n₁ + n₂. For each group g, the summed cosine and sine components of the transformed angles were computed as C_g = Σ_i ∈ g cos(ϕ_i), S_g = Σ_i ∈ g sin(ϕ_i).

The Mardia–Watson–Wheeler statistic was calculated as

which tests the null hypothesis that the two samples are drawn from the same circular distribution.

Statistical significance was evaluated using a permutation procedure. Group labels were randomly permuted across the pooled angles while preserving the original group sizes, and the test statistic was recomputed for each permutation. This procedure was repeated 5,000 times using a pseudo-random number generator. The P value was estimated as P = (k + 1)/(N_perm + 1), where k denotes the number of permutations yielding a statistic greater than or equal to the observed value and N_perm is the total number of permutations. Analyses were performed in Python using custom code based on the rank-based formulation of the Mardia–Watson–Wheeler statistic. Note that because angles compared here are independent of the magnitude of t-statistics, they are not subject to differences of signal-to-noise ratios from different data modalities.

Pupil tracking

Pupil video was acquired at 22–24 frames per second with a telecentric lens (Edmund Optics 58-430). Time synchronization was performed by adding an LED reflection into the mouse’s eye, next to the pupil, which turns on at the start of each trial. Synchronization signal and pupil edges were detected using a DeepLabCut model (v2.1.6.2; RRID:SCR_021391). Two models were trained using 314 and 208 frames from 9 and 13 mice, respectively, and were further refined based on tracking performance. Pupil boundaries were labelled on the left-most and right-most edges. Edge detection confidence lower than 0.9 was excluded from analysis. Signals were low-pass filtered with a 2nd-order Butterworth filter (cutoff, 5 Hz) using zero-phase forward–backward filtering.

Pose estimation

We collected video (SpinView 1.29.0.5) during task performance at a framerate of 500 Hz and resolution of 720 × 540 pixels using FLIR (forward-looking infrared) cameras (Blackfly S BFS-U3-04S2M). We trained a Lightning Pose model (RRID:SCR_024480)¹²⁵ to detect the distal tip of the tongue at the midline (mean error of 4.47 pixels). The model was trained on 1,200 labelled frames across 12 behavioural sessions and 8 mice. Tongue tip location was predicted across all analysed sessions, and subsequent analysis was limited to sessions with high model performance as measured by >90% mutual agreement with the capacitive lick sensor and >60 ms median tongue movement duration, which excluded sessions with poor model generalization.

For signal processing, kinematic data (tongue tip position over time) was masked for model confidence at 0.90. Data were then filtered using a 50 Hz, low-pass, symmetric, fourth-order Butterworth filter. Reaction time was defined per trial as the time elapsed from the go cue until the first timepoint in which the tongue was detected, with maximum cutoff set at 1.0 s.

Correlations between spike counts and reaction time were nonparametric. Baseline window was defined as 1.0 s prior to the go cue. Response window was defined as 0.2 s following the go cue, prior to the median reaction time across sessions. t-statistics were calculated from linear regression of reaction time to spike counts in specified windows on each trial.

Lick bout detection

Licks recorded from lick sensor contact were first processed to remove lick detection noise either from signal rebound or crosstalk between two lick sensors. Licks detected from pose tracking in videos were constrained to licks with (1) existing time longer than 0.02 s and shorter than 95th percentiles and (2) end of lick movement location, peak velocity, and total distance all smaller than 95th percentiles. Curated licks detected by lick sensors or by video were segmented into lick bouts where inter-lick-interval was longer than 0.5 s.

Data analysis

All data are presented as mean ± s.e.m. unless reported otherwise. All statistical tests were two-sided. For all analyses, no-go cues (when presented) were treated as part of the inter-trial interval. Code is available at https://github.com/AllenNeuralDynamics/aind_stan_fit_sim, https://github.com/AllenNeuralDynamics/LC-beh-physiology-analysis and https://github.com/JeremiahYCohenLab/sueAnalysis/tree/master/python/pupillometry.

Analyses and models of behaviour

We fit logistic regression models to predict choices and engagement as a function of outcome history for each mouse. To predict choices, we used this model:

where c_r(t) = 1 was a rightward choice and 0 was a leftward choice. R = 1 was a rewarded choice and 0 was an unrewarded choice, and C = 1 was a rightward choice and −1 was a leftward choice.

To predict engagement, we used this model:

where c(t) = 1 was a response to the go cue and 0 for ignoring a go cue, and R = 1 was a rewarded choice regardless of direction.

To predict switch versus stay choices, we used this model:

where s(t) is a switch choice, N_ipsi = 1 was an unrewarded choice ipsilateral to the choice at t − 1, −1 for an unrewarded choice contralateral to the choice at t − 1, and 0 for a rewarded choice, c_ipsi = 1 was a choice made ipsilateral to the choice at t − 1, −1 for an unrewarded choice contralateral to the choice at t − 1.

We applied a family of reinforcement-learning models of behaviour used previously^40,43. This model estimates action values (Q_l(t) and Q_r(t)) on each trial to generate choices. Choices are described by a random variable, c(t), corresponding to left or right choice, c(t) ∈ {l r}. The value of a choice is updated as a function of the RPE, and the rate at which this learning occurs is controlled by the learning rate parameter α. To account for asymmetric learning from rewards and no rewards, we used separate learning rates for each outcome. For example, if the left spout was chosen, then

where δ(t) = R(t) − Q_l(t) and ζ represents the forgetting rate parameter. The forgetting rate captures the increasing uncertainty about the value of the unchosen spout.

The Q-values were used to generate choice probabilities through a softmax decision function:

where β, the inverse temperature parameter, controls the steepness of the sigmoidal choice function. In other words, β controls the stochasticity of choice. This model was used for experiments with electrophysiology and GCaMP measurements.

To compare change in choice dynamics against RPE, we estimated the change of choice likelihood in consecutive trials. This quantity, defined as

determines how future behaviour is altered following RPE. \({\mathcal{L}}(c(t+1)=c(t))\) was the likelihood of a mouse choosing current choice in the future, estimated from the sequence of mouse choices using logistic regression model. \({\mathcal{L}}(c(t))\) represents animals’ choice policy before outcome, it was the likelihood estimated from the Q-learning model. Change of likelihood was centred in each session to account for choice autocorrelations.

Behavioural model fitting

We fit and assessed models using python and the probabilistic programming language, Stan (https://mc-stan.org/; RRID:SCR_018459) with the PyStan interface. Stan was used to construct hierarchical models with mouse-level hyperparameters to govern session-level parameters. This hierarchical construction uses partial pooling to mitigate overfitting to noise in individual sessions (often seen in the point estimates for session-level parameters that result from other methods of estimation) without ignoring meaningful session-to-session variability. For each session, each parameter in the model (for example, α) was modelled as a draw from a mouse-level distribution with mean μ and variance σ. Models were fit using noninformative (uniform distribution) priors for mouse-level hyperparameters (Supplementary Table 3).

Weakly informative priors were used for session-level parameters. Mouse-level hyperparameters were chosen to achieve model convergence under the assumption that individual mice behave similarly across days. The parameters were sampled in an unconstrained space and transformed into bounded values (for those parameters that were bounded) by a standard normal inverse cumulative density function. Stan uses full Bayesian statistical inference to generate posterior distributions of parameter estimates using Hamiltonian Markov chain Monte Carlo sampling¹²⁶. The default no-U-turn sampler was used. The Metropolis acceptance rate was set to 0.85–0.9 to force smaller step sizes and improve sampler efficiency. The models were fit with 5,000 iterations and 2,500 warmup draws run on each of 16 chains in parallel. Default configuration settings were used otherwise (https://github.com/AllenNeuralDynamics/aind_stan_fit_sim).

Extracting model parameters and variables, behaviour simulation

For extracting model variables (like RPE), we took at least 1,000 draws from the Hamiltonian Markov chain Monte Carlo samples of session-level parameters, ran the model agent through the task with the actual choices and outcomes, and averaged each model variable across runs. For comparisons of individual parameters across models, we estimated maximum a posteriori parameter values by approximating the mode of the distribution: binning the values in 50 bins and taking the mean value of the most populated bin.

Linear regression models of spike rates

To characterize neurons’ correlation with task components, we performed linear regression in different time windows. To account for the potential contribution of slow drift of the probe, mean spike amplitude around each trial was included as a regressor and contributed to prediction of spike counts in on average 25.4% of all neurons across different regression analysis.

To determine how neurons responded to different outcomes (presence or absence of reward) from different choices (ipsilateral or contralateral to the recording site), we regressed firing rates on outcome (R), choice side c(t), 1 for the direction ipsilateral to the recording site, −1 for contralateral to the recording site, reward prediction (Q_c), and their interaction, using the Python package scikit-learn (RRID:SCR_002577). Because there was no external sensory information to indicate that the mouse successfully made a choice and that water would or would not arrive, we defined separate windows for analysis of reward and no reward. Regressions were performed for each neuron in both a reward response window and a no-reward response window and the one with larger absolute t-statistics were used for each neuron.

To quantify how each neuron responds to outcome, the area under the receiver operating characteristic curve (AUROC) was computed in a sliding window of 1 s, from response time to 2 s after the response time (Extended Data Fig. 12). Each neuron’s peak outcome response time was when the AUROC values were furthest from 0.5. A neuron’s outcome response sign was defined as positive if this AUROC value was greater than 0.5 and negative if it was less than 0.5. Across the population, we estimated the distribution of outcome response time in neurons with positive and negative response signs. The modes of the distributions were taken as the positive and negative outcome response windows.

We quantified regressions using both t-statistics of regressors and the angle in Cartesian coordinates of the estimates of coefficients for R and Q_c for each neuron. The latter was plotted as a polar histogram (Fig. 6).

To quantify each neuron’s relationship to task engagement, we performed linear regression on the baseline firing rate 2 s before each go cue with presence or absence of a lick response to the go cue as a regressor.

To quantify each neuron’s relationship to choice (switch versus stay), we performed linear regression on the go cue period’s spike rate (0.5 s after go cue) with choice as a regressor.

Linear regression models of axonal GCaMP

To compare axonal activity correlation with behaviour, we performed linear regression in different time windows. To address movement artefacts, the isosbestic channel’s signal from same time windows were included as additional regressor.

We applied a quality control threshold on the LC-NE axon dynamics, including data with signal increases (z-scored across the whole recording) higher than 0.5 z-score ΔF/F.

To determine how LC-NE axons responded to different outcomes, we calculated regressions as for the single-neuron analysis (presence or absence of reward, ipsilateral or contralateral choice relative to the recording site, their interactions, and chosen value). Because axonal ΔF/F has lower temporal resolution than single-neuron spiking, a single larger time window (2 s) after outcome was used to compute the regression model for outcome encoding.

To quantify how LC-NE axons correlated with task engagement, we calculated linear regression on the baseline firing rate 2 s before each go cue with presence or absence of a lick response to the go cue as a regressor.

To quantify how LC-NE axons correlated with choice (switch versus stay), we performed linear regression on the go cue period’s firing rate (0.75 s after go cue) with choice as regressor.

We quantified regressions using both t-statistics of regressors and the angle in Cartesian coordinates of the estimates of coefficients for R and Q_c for each neuron. The latter was plotted as a polar histogram (Fig. 6).

Pupil diameter analysis with behaviour

To test how pupil diameter changes with behaviour events. A linear regression model was fitted to a sliding window of 1 s across the trial time, aligned to the start of go cues. Regressors included choice outcome, choice side (ipsilateral or contralateral to pupil recording side), and whether mice changed choice or repeated the previous choice. Distribution of t-statistics are reported in Extended Data Fig. 13f.

Pupil auto-correlogram

To capture temporal statistics of pupil diameter, auto-correlogram were computed across the whole session of pupil diameter recording with a sliding window of 0.5 s, with period of poor time alignment or poor pupil edge detection removed. An exponential curve was fitted to the autocorrelation curve: pupil = Ae^−t/τ + C. τ had a mean of 3.37 and s.d. of 1.03.

Pupil correlation with neuron activity

To analyse how LC-NE activity correlated with pupil diameter on different timescales, we performed three types of correlation analysis focusing on different time windows. First, to have a general estimate of task-independent correlation between pupil diameter and neuron activity on long timescales, we calculated cross-correlations between the two. Mean spike rate and mean pupil diameter were computed using a sliding window of 5 s. Cross-correlations were estimated with a maximum lag of 20 s, advancing in steps of 0.1 seconds. Pupil diameter was first de-trended using a linear fit of time to remove the effects of slow pupil diameter changes on the timescale of the whole session. This correlation was computed both including and excluding the trial period (0–3 s from go cues).

Second, to quantify spike-pupil relationships during inter-trial-intervals before go cues, we computed the trial-wise correlation between baseline spike rate and baseline pupil diameter 1 s before the onset of the go cue.

Third, to compare baseline neuronal activity to pupil dilation, we computed trial-wise correlations of spike rates 1 s before go cues with pupil dilation (maximum pupil increase from baseline pupil diameter, 1 s before go cue).

Comparison of phasic and tonic responses

To quantify how tonic activity (pre-trial spike rate) correlated with phasic responses to the go cue, we calculated correlation coefficients 0.5 s before the go cue 0–0.3 s after the go cue. We compared against a null model where we computed this correlation with ‘sham’ windows, randomly sampled from a uniform distribution between the first and last go cue.

Analysing relationships with space

We analysed the relationship between space and both one-dimensional and high-dimensional data. To test whether a particular feature (or group of features) was randomly distributed in space or instead shows spatial structure, we quantified two complementary forms of spatial dependence using the corresponding CCF coordinates X = (x, y, z). The analysis assessed (1) the presence of a global linear trend across space and (2) the extent to which values were more similar to their neighbours than randomly assigned neurons. All significance testing was performed using permutation tests that preserve the spatial sampling geometry by keeping coordinates fixed while randomly reassigning values across locations. When assessing a significant linear trend, we also estimated the primary spatial vector and its confidence interval for each feature (or group of features). Any neuron containing non-finite entries in either its location or feature values were removed. Analyses were only performed when at least 100 valid spatial samples remained after filtering. Medial–lateral coordinates were mirrored to the left hemisphere by taking the absolute value of ML, so that left and right locations contributed to the same axis estimate.

Global linear trend test for single features using linear regression

We first asked whether y exhibits a global linear gradient across space by fitting an ordinary least squares (OLS) regression model, y = β₀ + βX + ϵ, where β₀ is an intercept and β ∈ R^D are slope coefficients for each spatial dimension. The model fit was summarized using the coefficient of determination \({R}_{\mathrm{obs}}^{2}\), which quantifies the fraction of variance in y explained by a linear function of space. β was normalized and used as the primary axis of the spatial trend.

To determine whether the observed linear trend exceeded what would be expected by chance given the same spatial sampling, we constructed a null distribution by repeatedly permuting y while leaving X fixed. For each permutation b = 1,…, B (B = 5,000), we computed \({R}_{b}^{2}={R}^{2}({OLS}({\pi }_{b}(y) \sim X))\), where πb(⋅) denotes a random permutation. A one-sided permutation P value was then computed as \({p}_{\mathrm{trend}}=\frac{\left(1+\mathop{\sum }\limits_{b=1}^{B}[{R}_{b}^{2}\ge {R}_{\mathrm{obs}}^{2}]\right)}{(B+1)}\), testing whether the observed explained variance is larger than expected under spatial randomness. Reported outputs include unit primary spatial axis, \({R}_{{\rm{obs}}}^{2}\), F-statistic, permutation P value, and the mean and s.d. of the permuted R² distribution.

Local spatial predictability test for single features using kNN test

Global linear trends can miss nonlinear or locally structured spatial patterns (for example, patches or clusters). Therefore, we additionally quantified local spatial dependence by testing whether values at a location were predictable from nearby locations using distance-weighted kNN regression. We performed K-fold cross-validation (default K = 5), shuffling samples with a fixed random seed for reproducibility. For each fold, a kNN regressor was trained on the training set using the spatial coordinates X, and predictions were generated for held-out points. The number of neighbours was set to k = min(k_neighbors, n_train), (default k_neighbors = 20 for analyses of spike rates during behaviour and 30 for spike waveform analysis, unless noted otherwise) with inverse-distance weighting so that closer neighbours contribute more strongly.

Predictive performance was summarized as the cross-validated coefficient of determination, \({R}_{\mathrm{cv},\mathrm{obs}}^{2}={R}^{2}(y,{\hat{y}}_{\mathrm{cv}})\), where \({\hat{y}}_{\mathrm{cv}}\) denotes concatenated out-of-fold predictions across all folds (default, 5). To evaluate whether local predictability exceeded chance, we generated a permutation null distribution by permuting \(y\) across locations and recomputing the full cross-validated kNN R² for each permuted dataset: \({R}_{\mathrm{cv},b}^{2}={R}_{\mathrm{cv}}^{2}({\pi }_{b}(y))\).

A one-sided P value was computed analogously: \({P}_{\mathrm{cv}}=(1\,+\)\({\Sigma }_{b=1}^{B}I[{R}_{\mathrm{cv},b}^{2}\ge {R}_{\mathrm{cv},\mathrm{obs}}^{2}])/(B+1)\). Reported outputs include \({R}_{\mathrm{cv},\mathrm{obs}}^{2}\), permutation P value, and the mean and s.d. of the null distribution.

Global linear trend test for single features using LDA

To identify whether categorical labels exhibited a structured spatial organization, we estimated a primary spatial axis that best separated labels using linear discriminant analysis (LDA). Specifically, we modelled the relationship between class labels y and spatial coordinates X ∈ ℝ^n × 3 by fitting an LDA model, which finds a linear projection of the coordinates that maximizes between-class variance relative to within-class variance.

Given coordinates X and categorical labels y, LDA estimates a set of discriminant directions β ∈ ℝ³ such that the projected values Xβ optimally separate the label classes. For multiclass settings, the number of discriminant components is bounded by min(D, K − 1), where D = 3 is the number of spatial dimensions and K is the number of unique classes. In all analyses, we retained only the first discriminant component, corresponding to the direction of maximal class separability. The resulting discriminant vector β was extracted from the fitted model (using the feature-space projection matrix when available) and normalized to unit length to define the primary spatial axis as \(\hat{v}=\beta /{||}\beta {||}\).

The unit vector \(\hat{v}\) represents the dominant spatial direction along which class labels were most strongly differentiated. For consistency with linear regression-based analyses, we also report the raw coefficients β, while the intercept term was set to zero as it is not meaningful in this geometric interpretation. Reported outputs include the unit primary spatial axis \(\hat{v}\), the unnormalized discriminant vector β, and model validity checks based on sample counts and class balance.

Global linear trend test for groups of features

To test whether a group of related variables (for example, different waveform features or expression of multiple genes) had a linear trend in space, we applied CCA. Values were z-scored prior to analysis. CCA was performed using the scikit-learn implementation (sklearn.cross_decomposition.CCA), extracting canonical variates that maximize correlation between linear combinations of functional and anatomical variables. The strength of association was quantified as the linear correlation between paired canonical projections. The weight from the spatial location was then normalized and used as the primary axis of spatial variation.

To assess statistical significance, we performed a permutation test in which anatomical labels were randomly shuffled across units (10,000 permutations). For each shuffle, CCA was recomputed, and the maximum canonical correlation was recorded to generate a null distribution. Empirical P values were computed as the proportion of shuffled correlations exceeding the observed correlation.

Bootstrap estimation of spatial vector confidence intervals

To estimate confidence intervals for primary spatial vectors inferred from the methods mentioned above, we performed bootstrap resampling with replacement across units (5,000 iterations). In each iteration, rows of the original dataset (one row for each unit) were resampled with replacement, and spatial vectors were re-estimated using corresponding methods. Because the sign of axis direction was arbitrary, each bootstrapped axis was aligned to the observed axis by flipping its sign whenever its inner product with the observed axis was negative. This ensured that all bootstrapped axes were in the same hemisphere and prevented artificial bimodality in the bootstrap distribution.

To visualize uncertainty, bootstrapped axis distributions were plotted in azimuth-elevation space and summarized as a 95% confidence cone around the observed axis (Extended Data Fig. 15). The cone half-angle was defined as the 95th percentile of the angular deviation between the bootstrapped axes and the observed axis. For plotting in 2D anatomical planes, the boundary of this 3D cone on the unit sphere was projected into each plane and displayed as a shaded confidence region around the projected axis arrow (Fig. 3 and Extended Data Fig. 15).

Comparison of two primary spatial axes

To compare whether two inferred primary spatial axes differed in direction, each pair of unit vectors from bootstrapping was represented in a local tangent plane coordinate system. Given a reference axis b_x (one of the primary spatial axes to compare), we constructed an orthonormal basis (e₁, e₂, u₀), where u₀ = b_x/||b_x|| is the normalized reference axis. The first tangent direction e₁ was defined as the normalized projection of the second axis onto the plane perpendicular to u₀,

and the second tangent direction was defined as e₂ = u₀ × e₁. This basis defines a two-dimensional plane orthogonal to the reference axis. The observed directional deviation of b_y from b_x was represented by projecting the second axis into this plane, d_obs = Ab_y, where A is a 2 × 3 projection matrix,

Bootstrap uncertainty for the difference between axes was estimated from the corresponding bootstrap axis estimates \({b}_{x}^{(i)}\) and \({b}_{y}^{(j)}\). Each bootstrap vector was first aligned to its observed axis by flipping its sign when the inner product with the observed axis was negative. Bootstrap vectors were then projected into the same tangent plane, \({p}_{x}^{(i)}=A{b}_{x}^{(i)},{p}_{y}^{(j)}=A{b}_{y}^{(j)}\).

Since bootstrap samples were not paired across datasets, the sampling distribution of directional differences was approximated by randomly sampling combinations of bootstrap projections from the two datasets and computing \({d}_{\mathrm{boot}}={p}_{y}^{(j)}-{p}_{x}^{(i)}\). This procedure generates a bootstrap cloud that approximates the sampling distribution of the two-dimensional directional difference. In the case where there are 2,000 bootstraps, 10,000 pairs of bootstrapped vectors were randomly sampled to approximate the distribution.

The magnitude of the observed directional difference was quantified using the Mahalanobis distance of the observed deviation from zero, \(W={d}_{\mathrm{obs}}^{\top }\,{\Sigma }^{-1}{d}_{\mathrm{obs}}\), where \(\Sigma =\mathrm{Cov}({d}_{\mathrm{boot}})\) is the covariance matrix of the bootstrap difference cloud.

To compute a bootstrap-based P value, the bootstrap difference distribution was first recentered to estimate the null distribution, \({d}_{\mathrm{boot},\mathrm{null}}={d}_{\mathrm{boot}}-{\bar{d}}_{\mathrm{boot}}\). For each recentered bootstrap sample, a quadratic form was computed \({W}_{\mathrm{boot}}={d}_{\mathrm{boot},\mathrm{null}}^{\top }{\Sigma }^{-1}{d}_{\mathrm{boot},\mathrm{null}}\). The bootstrap P value was then defined as

where N_boot is the number of bootstrap difference samples and I(⋅) is the indicator function.

The bootstrap difference cloud was also used directly to visualize uncertainty and to assess whether the 95% confidence ellipse included the origin, corresponding to no directional difference. Finally, the angular separation between axes was reported as \(\theta ={\cos }^{-1}\,({b}_{x}^{\top }{b}_{y})\), which provides an intuitive summary of directional difference between the two spatial axes.

Relating behavioural features to spatial axes

To assess how behavioural encoding varied along each primary spatial axis inferred from other data modalities, each neuron’s anatomical coordinate was projected onto each axis by taking the dot product between its 3D coordinate and the unit spatial vector: \({p}_{i}={x}_{i}^{\top }b\). This yielded a one-dimensional coordinate along the waveform axis, the MERFISH axis, and the retrograde label axis.

For each behavioural feature, the projected coordinate was related to the feature value using Pearson correlation. Correlation coefficients and corresponding P values were calculated using only samples with non-missing values. For visualization, a simple linear regression line was fit to the relationship between projected coordinate and behavioural feature, and a 95% confidence band for the fitted line was computed from the standard linear regression uncertainty estimate. Each behavioural feature was therefore compared across multiple spatial embeddings by repeating the same projection-and-correlation analysis for each inferred axis (Fig. 5 and Extended Data Fig. 15).

Spike count cross-correlations

To estimate spike rate correlations among simultaneously recorded LC-NE neurons, we counted spikes in non-overlapping windows of 50 ms and calculated cross-correlations between them. To extract the main features of cross-correlograms, PCA was applied.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.