Abstract
To understand how the brain produces behavior, we must elucidate the relationships between neuronal connectivity and function. The medial prefrontal cortex (mPFC) is critical for complex functions including decision-making and mood. mPFC projection neurons collateralize extensively, but the relationships between mPFC neuronal activity and brain-wide connectivity are poorly understood. We performed whole-brain connectivity mapping and fiber photometry to better understand the mPFC circuits that control threat avoidance in male and female mice. Using tissue clearing and light sheet fluorescence microscopy (LSFM), we mapped the brain-wide axon collaterals of populations of mPFC neurons that project to nucleus accumbens (NAc), ventral tegmental area (VTA), or contralateral mPFC (cmPFC). We present DeepTraCE (deep learning-based tracing with combined enhancement), for quantifying bulk-labeled axonal projections in images of cleared tissue, and DeepCOUNT (deep-learning based counting of objects via 3D U-net pixel tagging), for quantifying cell bodies. Anatomical maps produced with DeepTraCE aligned with known axonal projection patterns and revealed class-specific topographic projections within regions. Using TRAP2 mice and DeepCOUNT, we analyzed whole-brain functional connectivity underlying threat avoidance. PL was the most highly connected node with functional connections to subsets of PL-cPL, PL-NAc, and PL-VTA target sites. Using fiber photometry, we found that during threat avoidance, cmPFC and NAc-projectors encoded conditioned stimuli, but only when action was required to avoid threats. mPFC-VTA neurons encoded learned but not innate avoidance behaviors. Together our results present new and optimized approaches for quantitative whole-brain analysis and indicate that anatomically defined classes of mPFC neurons have specialized roles in threat avoidance.
SIGNIFICANCE STATEMENT Understanding how the brain produces complex behaviors requires detailed knowledge of the relationships between neuronal connectivity and function. The medial prefrontal cortex (mPFC) plays a key role in learning, mood, and decision-making, including evaluating and responding to threats. mPFC dysfunction is strongly linked to fear, anxiety and mood disorders. Although mPFC circuits are clear therapeutic targets, gaps in our understanding of how they produce cognitive and emotional behaviors prevent us from designing effective interventions. To address this, we developed a high-throughput analysis pipeline for quantifying bulk-labeled fluorescent axons [DeepTraCE (deep learning-based tracing with combined enhancement)] or cell bodies [DeepCOUNT (deep-learning based counting of objects via 3D U-net pixel tagging)] in intact cleared brains. Using DeepTraCE, DeepCOUNT, and fiber photometry, we performed detailed anatomic and functional mapping of mPFC neuronal classes, identifying specialized roles in threat avoidance.
Introduction
The medial prefrontal cortex (mPFC) is a vastly interconnected brain region that controls complex functions including memory, decision-making, and mood (Euston et al., 2012; Giustino and Maren, 2015; Klune et al., 2021). Through its axonal projections, mPFC exerts top-down control over downstream regions, promoting adaptive responses in dynamic circumstances (Ye et al., 2016; C.K. Kim et al., 2017; Murugan et al., 2017; Otis et al., 2017; Jayachandran et al., 2019; Diehl et al., 2020). While studies have linked individual mPFC projections to discrete behavioral functions (Warden et al., 2012; Riga et al., 2014; DeNardo et al., 2015; Ye et al., 2016; C.K. Kim et al., 2017; Murugan et al., 2017; Otis et al., 2017; Vander Weele et al., 2018; Siciliano et al., 2019; Diehl et al., 2020; V.P. Mathis et al., 2021), mPFC neurons are anatomically heterogenous, with extensive axon collaterals (Gao et al., 2022). Many mPFC-dependent behaviors involve coordinated activity in multiple downstream brain regions (Y. Kim et al., 2015). Because of its role in survival behaviors and neuropsychiatric disorders (Euston et al., 2012; Giustino and Maren, 2015; Klune et al., 2021), there is an urgent need to integrate mPFC neuroanatomy with functional studies to understand how classes of mPFC projection neurons contribute to behavior. To address this, we developed new approaches for whole brain circuit mapping. By combining these approaches with fiber photometry, we examined the unique roles of three mPFC cell classes in learned threat avoidance, a survival behavior that is highly relevant to fear and anxiety disorders (Diehl et al., 2019).
Many studies use retrograde viral approaches to manipulate or observe the activity of projection-defined populations (Zhang et al., 2011; Roth, 2016; Dana et al., 2019). In many cases, activation of axon collaterals may inadvertently influence animal behaviors. Interpreting results from such studies therefore requires understanding the full projection patterns of neuronal classes of interest. Traditional studies of axon collaterals established foundational knowledge but have been limited by low-throughput approaches that require predefining collateral targets. These include measuring overlap of retrograde tracers following injections into multiple targets (Gabbott et al., 2005; Murugan et al., 2017) or recording from source neurons while antidromically stimulating target regions (Ishikawa and Nakamura, 2006; Ishikawa et al., 2008). Serial two-photon tomography and single neuron reconstructions have greatly improved our understanding of prefrontal projection patterns. But these techniques are expensive and computationally and time-intensive, making them impractical for many labs (Oh et al., 2014; Gao et al., 2022).
Conversely, tissue clearing (Chung et al., 2013; Renier et al., 2016) and light sheet fluorescence microscopy (LSFM) are fast, easy, and affordable and can preserve three dimensional (3D) circuit architecture without the need to align structures. These approaches are useful for comparing projection patterns of genetically, anatomically, or behaviorally defined neurons and generating new hypotheses that can be tested using functional methods. To enhance their throughput, we developed DeepTraCE (deep learning-based tracing with combined enhancement) for quantifying bulk-labeled fluorescent axons and DeepCOUNT for quantifying fluorescently labeled cell bodies. DeepTraCE and DeepCOUNT build on the machine learning-based image segmentation package TrailMap (Friedmann et al., 2020) to form a user-friendly whole-brain analysis pipeline based on tissue clearing and LSFM.
Advances in computer vision allow researchers to train classifiers to detect axons in tissue volumes (Berg et al., 2019; Friedmann et al., 2020). However, axon morphology varies across brain regions (Gao et al., 2022), causing classifiers to fail. To overcome this, DeepTraCE combines multiple models that have been separately trained using TrailMap to accurately identify axons with differing appearances. After segmenting the axons using a combination of TrailMap models, DeepTraCE registers the brains to a common coordinate framework and quantifies innervation density by brain region. Similarly, DeepCOUNT uses TrailMap to detect cells and then quantifies cell bodies by brain region.
We used DeepTraCE, DeepCOUNT and fiber photometry to reveal novel insights into mPFC threat avoidance circuits. The prelimbic (PL) subregion of mPFC is integral to both learned and innate threat avoidance (Stern et al., 2010; Bravo-Rivera et al., 2014; Capuzzo and Floresco, 2020). Yet our understanding of the underlying PL circuit mechanisms is incomplete. We focused on PL neurons that project to the contralateral PL (cPL), nucleus accumbens (NAc), or ventral tegmental area (VTA). These neuronal classes have established roles in cognition (Murugan et al., 2017; Lak et al., 2020; Lui et al., 2021) and reward-seeking (C.K. Kim et al., 2017; Otis et al., 2017), but the relationship between their connectivity and function remains poorly understood, especially in aversive learning. Here, we mapped their brain-wide collateral projections. DeepTraCE accurately detected axons, outperforming other model combination methods. Using DeepCOUNT, we mapped brain-wide neuronal activation patterns as mice avoided threats and then constructed functional networks, finding that PL was one of the most highly connected nodes. Using fiber photometry, we discovered that projection-defined PL classes have specialized roles in threat avoidance. Further, we demonstrated the utility of combining common techniques like fiber photometry with DeepTraCE and DeepCOUNT to understand the detailed structural and functional connectivity of the cells being studied.
Materials and Methods
Animals
Female and male C57B16/J mice (JAX stock #000664) or TRAP2;Ai14 mice (JAX stock #030323 and #007914) were group housed (two to five per cage) and kept on a 12/12 h light/dark cycle (lights on 7 A.M. to 7 P.M.). All animal procedures followed animal care guidelines approved by the University of California, Los Angeles Chancellor's Animal Research Committee.
Surgery
Mice were induced in 5% isoflurane in oxygen until loss of righting reflex and transferred to a stereotaxic apparatus where they were maintained under 2% isoflurane in oxygen. Mice were warmed with a circulating water heating pad throughout surgery and eye gel was applied to the eyes. The head was shaved and prepped with three scrubs of alternating betadine and then 70% ethanol. Following a small skin incision, a dental drill was used to drill through the skull above the injection targets. A syringe pump (Kopf, 693A) with a Hamilton syringe was used for injections. Injections were delivered at a rate 75 nl/min and the syringe was left in the brain for 7 min following injection. For collateralization mapping, 300 μl of AAVrg-Ef1a-mCherry-IRES-Cre-WPRE (Addgene 55632-AAVrg, 1.7 × 1013 vg/ml) was injected unilaterally into either left NAc (AP: 1.3, ML: −1.0, DV: 4.7), left VTA (AP: −3.3, ML: 0.4, DV: −4.5), or right PL (AP: 1.8, ML: −0.5, DV: −2.3). A total of 200 μl of AAV8-hSyn-DIO-hCHR2(H134R)-EYFP-WPRE (1.7 × 1013 vg/ml) was then injected into left PL (AP: 1.8, ML: −0.4, DV: −2.3). To control for Cre-independent expression of EYFP, control mice received a single injection of 200 μl of AAV8-hSyn-DIO-hCHR2(H134R)-EYFP-WPRE into left PL. For fiber photometry from cPL and NAc-projectors, we injected AAVrg-Ef1a-mCherry-IRES-Cre-WPRE (same as above) or AAVrg-EF1a-mCherry-IRES-Flpo-WPRE (Addgene 55634-AAVrg, 1.7 × 1013 vg/ml) into either cPL or NAc, and after two weeks of expression injected AAV8-EF1a-Con-Foff-2.0-GCamp6m-WPRE into PL (Addgene 137120-AAV8). For recordings from PL-VTA neurons, we injected AAVrg-Ef1a-mCherry-IRES-Cre-WPRE (same as above) into VTA and AAV5-CAG-FLExGCamp6m-WPRE (Addgene 100839-AAV5) into PL; 400-μm optic fibers (Doric) were implanted into left PL (AP: 1.8, ML: −0.4, DV: −2.3) and sealed in place using Metabond (Patterson Dental Company, 5533559, 5533492, S371). For topographical mapping of the origin of each projection type, mice were injected with 300 μl of cholera toxin subunit B (CTB; ThermoFisher, C34775, C34776, C34778) or AAVrg-hSyn-EGFP (Addgene 50465) at the same coordinates for cPL, NAc, and VTA listed above. For pain management mice received 5 mg/kg carprofen diluted in 0.9% saline subcutaneously. Mice received one injection during surgery and daily injections for 2 d following surgery. Samples with mistargeted injection sites were excluded from analysis. Samples with obviously poor antibody penetration or distribution following tissue clearing were also excluded.
Brain slice histology and immunostaining
Mice were transcardially perfused with PBS followed by 4% paraformaldehyde (PFA) in PBS. Brains were dissected, postfixed in 4% PFA for 12–24 h and placed in 30% sucrose for 24–48 h. They were then embedded in Optimum Cutting Temperature (OCT; Tissue Tek) and stored at −80°C until sectioning; 60-μm floating sections were collected into PBS. Sections were washed 3 × 10 min in PBS and then blocked in 0.3% PBST containing 10% normal donkey serum (Jackson ImmunoResearch, 17-000-121) for 2 h. Sections were then stained with chicken anti-GFP (AVES 1020 at 1:2000), rabbit anti-RFP (Rockland 600-401-379 at 1:2000) or rat anti-CTIP2 (Abcam ab18465, 1:200) in 0.3% PBST containing 3% donkey serum overnight at 4°C. The following day, sections were washed 3x5 min in PBS and then stained with secondary antibody (Jackson ImmunoResearch Cy2 donkey anti-chicken IgG(H+L) 703-225-155, 1:1000, Cy3 donkey anti-rabbit IgG(H+L) 711-005-152, 1:1000 or Alexa647 donkey anti-rat IgG(H+L) 702-605-150, 1:500) in 0.3% PBST containing 5% donkey serum for 2 h at room temperature. Sections were then washed 5 min with PBS, 15 min with PBS+DAPI (Thermofisher Scientific, D1306, 1:4000), and then 5 min with PBS. Sections were mounted on glass slides using FluoroMount-G (ThermoFisher, 00-4958-02) and then imaged at 10× with a Leica STELLARIS confocal microscope or at 5× on a Leica DM6 B scanning microscope.
Brain clearing
Mouse brains were collected and processed based on the published Adipo-Clear protocol (Chi et al., 2018a) with slight modifications. Mice were perfused intracardially with 20 ml of PBS (Invitrogen) followed by 4% paraformaldehyde (PFA; Electron Microscopy Sciences) on ice. Brains were hemisected ∼1 mm past midline and postfixed overnight in 4% PFA at 4°C. The following day, samples were dehydrated with a gradient of methanol (MeOH, Fisher Scientific):B1n buffer (1:1000 Triton X-100, 2% w/v glycine, 1:10,000 NaOH 10N, 0.02% sodium azide) for 1 h for each step (20%, 40%, 60%, 80%) on a nutator (VWR). Samples were then washed with 100% MeOH 2× for 1 h each and then incubated in a 2:1 dicholoromethane (DCM):MeOH solution overnight. The following day, two washes of 1 h in 100% DCM were performed followed by three washes of 100% MeOH for 30 min, 45 min then 1 h. Samples were bleached for 4 h in 5:1 H2O2/MeOH buffer. A cascade of MeOH/B1n washes (80%, 60%, 40%, 20%) for 30 min each rehydrated the samples followed by a 1 h wash in B1n buffer. 5%DMSO/0.3 M. Glycine/PTxWH permeabilized tissue for 1 h and then again for 2 h with fresh solution. Samples were washed with PTxwH for 30 min and then incubated in fresh PTxwH overnight. The following day, two more PTxwH washes lasted 1 h, then 2 h. For axon analysis, samples were incubated in primary GFP antibody (AVES Labs GFP 1020) at 1:2000 in PTxwH shaking at 37°C for 11 d, washed in PTxwH 2× 1 h and then 2× 2 h, then for 2 d with at least one PTxwH change per day while shaken at 37°C. For TRAP2 experiments, samples were incubated in primary RFP antibody (Rockland 600-401-379) at 1:300 in PTxwH shaking at 37°C for 11 d, washed in PTxwH 2× 1 h and then 2× 2 h, then for 2 d with at least one PTxwH change per day while shaken at 37°C. Samples were then incubated in secondary antibody (AlexaFluor 647, ThermoFisher Scientific) for 8 d shaken at 37°C. Samples were washed in PTxwH 2× 1 h, then 2× 2 h, then 2 d with at least one PTxwH change per day while shaken at 37°C. Samples were then washed in 1× PBS twice 1× 1 h, 2× 2 h and then overnight. To dehydrate samples, a gradient of washes in MeOH:H2O (20%, 40%, 60%, and 80%) were conducted for 30 min each, followed by 3 × 100% MeOH for 30 min, 1 h, then 1.5 h. Samples were incubated overnight in 2:1 DCM:MeOH on a nutator. The next day, samples were washed in 100% DCM 2× 1 h each. Samples were then cleared in 100% DBE. DBE was changed after 4 h. Samples were stored in DBE in a dark place at room temperature. Imaging took place at least 24 h after clearing.
Whole-brain imaging
Brain samples were imaged on a light-sheet microscope (Ultramicroscope II, LaVision Biotec) equipped with a sCMOS camera (Andor Neo) and a 2×/0.5 NA objective lens (MVPLAPO 2×) equipped with a 6-mm working distance dipping cap. Image stacks were acquired at 0.8× optical zoom using Imspector Microscope v285 controller software. For axons, we imaged using 488-nm (laser power 20%) and 640-nm (laser power 50%) lasers. The samples were scanned with a step-size of 3 µm using the continuous light-sheet scanning method with the included contrast adaptive algorithm for the 640-nm channel (20 acquisitions per plane), and without horizontal scanning for the 488-nm channel. For TRAP2-Ai14 brains, the 640-nm channel was imaged at 20% laser power.
Model training
Light sheet images of fluorescently labeled axons were segmented using the 3D U-net based machine learning pipeline TrailMap (Friedmann et al., 2020). We trained new models for segmentation of cortical axons using inference learning as described in the TrailMap pipeline (https://github.com/AlbertPun/TRAILMAP). Each training dataset included 6–20 hand-labeled image cubes (120 × 120 × 100 pixels). Training cubes were selected from a variety of brain regions across nine different mice. Each brain region was represented across multiple mice. A total of 10–20% of the image cubes were held out for use in model validation. Axons were hand-labeled in three or four planes from each cube. Where necessary, artifacts (pixels that could be mistakenly interpreted as axons, often bright clumps of antibody near the edge of the brain) were given a separate label. One-pixel-wide edges of axons were labeled using python to be given less weight in training of the 3D convolutional network, which accounts for slight variability in human annotation patterns.
We performed iterative training “sessions,” with each session having unique training data targeted to the weaknesses of the previous model. Each training session consisted of 5-20 epochs (20-120 steps per epoch). We selected the best model from each session by plotting loss in validation data across epochs. Loss was minimized to prevent overtraining. We used the weights from the best model as the starting point for the next session. Sessions were continued until all visible axons in the desired regions were reliably detected by the model.
We observed differences in axon appearance that likely stemmed from a combination of the speed-resolution trade-off of light sheet microscopy and the physical properties of the axons. Based on careful visual analysis of axons in different brain regions, we found that most axons were either bright and highly delineated, fuzzy and indistinct, or fell in the middle of these two categories. We therefore predicted that three models would best capture the axons in our dataset. We classified each brain region into one of three groups: delineated, moderately distinguished, or indistinct axons, and we trained Models 1, 2, and 3 to segment these types of axons, respectively. This empirical approach of classifying axons into three groups and then combining the models outperformed other model combination methods, as described below.
To train our models, we used the published TrailMap model weights as a starting point. First, seven training sessions were performed with data from regions assigned to Model 2. The result of the final session was used as Model 2. To generate Model 1, we used the weights of Model 2 as the starting point. Two sessions of training were performed using training data from regions assigned to Model 1. To generate Model 3, the weights of Model 2 were used as the starting point and four sessions of training were performed using training data from regions assigned to Model 3. In the DeepTraCE software, users can select the number of models they would like to use and assign brain regions accordingly. Or they can simply use a single model.
The TRAP2-Ai14 model was trained in the same manner. Initial weights were derived from a model that had been trained on fos+ cells. From this, six sessions of training with 20 steps per epoch and 4–150 epochs per session were performed to generate the TRAP2-Ai14 model based on tdTomato+ TRAPed cells.
Axon model validation
We validated our analysis pipeline including image segmentation, scaling and axon thinning by selecting 120 × 120 × 100-pixel cubes from regions assigned to the “delineated,” “moderately distinguished,” and “indistinct” groups in brains that were not included in initial training or validation datasets. Two human experts annotated two to four planes from each stack, marking pixels likely to contain axons. In stacks with difficult-to-distinguish axons, such as deep regions, experts used larger brush strokes to label broader groups of pixels likely to contain axons. One-pixel edges were then added to human annotations in python.
The raw image stacks were segmented as explained below. In brief, each stack was segmented with Models 1, 2, and 3 using TrailMap and scaled in ImageJ to a 10-μm space. Scaled images were skeletonized in python at 8 different thresholds for binarization, combined in MATLAB, and pixel values were adjusted in ImageJ. Using the same quantification threshold as used in all analyses (number of skeletonized pixels with intensity above 64 divided by total number of pixels in the region), labeling density was calculated in human-annotated “axon-positive” pixels (LP) and “axon-negative” pixels (LN), excluding edges to account for slight variability in human labels. Distinction score was calculated from these values as (LP − LN)/(LN + (1 − LP)). While direct true/false positive and true/false negative rates cannot be calculated because of the axon thinning process, which reduces the number of labeled pixels, LP correlates with the true positive rate, LN correlates with the false positive rate, and 1 − LP correlates with the false negative rate. Thus, LP − LN in the numerator will increase distinction score when better separation is obtained between true and false positives, and LN and 1 − LP in the denominator will decrease the distinction score if there is a higher rate of false positives or false negatives.
Alternate model combination methods were performed in FIJI. Maximum probability projections between all models were calculated by generating images with maximum pixel values across all models for each cube then processing the cubes for validation as above. Summation of probability maps between models was calculated by adding probability calculated by each model, then dividing by three before processing the cubes for validation as above. Distinction scores from alternate approaches were compared using repeated measures one-way ANOVA with Benjamini, Krieger, and Yekutieli posthoc test, comparing individual model selection (DeepTraCE) with each other method. Statistical comparisons of distinction scores were performed using GraphPad PRISM.
TRAP2-Ai14 model validation
To validate that cell counts produced by DeepCOUNT were accurate, we selected 120 × 120 × 100-pixel image cubes from several brain regions across multiple brains not included in the training dataset. Cubes were scaled and processed in the same manner as whole-brain data to produce a raw cell count from each cube. Two human experts then manually counted cells in each cube. Average counts from the two human experts were then compared with counts produced by DeepCOUNT.
DeepTraCE analysis pipeline
Image segmentation and registration
Whole-brain image stacks from the 640-nm channel from each brain were segmented three times in TrailMap, once each using Models 1, 2, and 3. Following segmentation, the 488 nm autofluorescence channel and axon segmentations from each model were converted to eight-bit and scaled to 10-μm resolution in FIJI with scaling values of 0.40625 in the x and y directions and 0.3 in the z direction using a bilinear interpolation algorithm. To improve image registration, each scaled 488 nm image was manually rotated in the x-, y-, and z-planes using the TransformJ ImageJ plugin (https://imagescience.org/meijering/software/transformj/) such that the midline blood vessels visible were all visible in the same z plane. The same manual rotation parameters were then applied to the three scaled model segmentation images from the corresponding brain.
The scaled and manually rotated 488 nm autofluorescence channel was registered using elastix to the Gubra Lab LSFM atlas average template, which has annotations based on the Allen CCF (Klein et al., 2010; Wang et al., 2020; Perens et al., 2021). The same transformation was applied to the scaled and manually rotated model segmentations using transformix. Image registration quality was manually verified by overlaying the atlas image and the registered 488-nm channels in ImageJ. Following segmentation and registration of each model, the transformix images were converted to 8-bit .tif format in ImageJ. When combining multiple probability maps, it is important that minimum and maximum pixel values are comparable between images. Model 1 produced slightly lower maximum probability values compared with the other two models. We corrected for this by brightening segmentations from Model 1 by 34% to match pixel values obtained by the other models before combination and thinning. We provide a unified python pipeline for automating these steps in the supporting software repository.
Model combination and thinning
Converted segmentations from the three models were then combined by generating a new image in which pixel values in each region were extracted from the model with best performance in that region (Fig. 1D). The same regional model assignments were used for all brains in the data set. Regional model assignments are provided in Extended Data Figure 3-2. Following model combination, axon segmentations were thinned (or “skeletonized”) as described in the TrailMap pipeline (Friedmann et al., 2020). In brief, images were binarized at eight different thresholds from 20–90% of the maximum intensity value using python. Skeletons were combined in MATLAB by summing values from each skeleton. Small objects unlikely to be axons were removed by calculating connected components within the combined skeleton and removing objects <90 voxels in size. Combined skeletons with small objects removed were optimized for visualization and quantification using an ImageJ macro that multiplied each pixel value by 17.
Axon quantification
Regional axon innervation was quantified in MATLAB (MathWorks) by counting the number of skeletonized pixels in each brain region above a threshold (64), then dividing this pixel count by the total number of pixels in a region. This regional pixel count was then divided by the total number of labeled pixels across the brain to normalize for differences in total fluorescence and viral expression. Regions were defined by a collapsed version of the LSFM atlas in which maximum granularity was balanced with the need to account for slight differences in registration which would lead to inaccurate quantification of small brain regions. This atlas was cropped on the anterior and posterior ends to match the amount of tissue visible in our data. Fiber tracts, ventricular systems, cerebellum, and olfactory bulb were excluded from analysis.
All statistical comparisons were performed in MATLAB and GraphPad Prism v9. For comparison of regional axon labeling between cell types, two-way ANOVA with Tukey's multiple comparison correction was performed on all brain regions. Anterior-posterior axon distributions within regions were calculated in MATLAB by binning the whole-brain image into 100-μm voxels and calculating the percentage of segmented pixels within each voxel. Voxels falling within a given region were summed across the medial-lateral and dorsal-ventral axis and normalized for total fluorescence as above. The averaged summation of axon counts from a given cell class was then averaged and plotted along with the SEM.
Axon visualization
To visualize axons, z-projections of raw light sheet data were created in FIJI by scaling images to a 4.0625-μm space, virtually reslicing images in the coronal plane, and performing maximum intensity z-projections of 100-μm depth followed by local contrast enhancement. Axon segmentations were created by overlaying skeletonized and registered axon segmentations from each sample of a cell type in a slightly different color and virtually reslicing in the coronal plane. 3D projections were created in Imaris using a representative sample from each cell type. Dotogram overlays (Figs. 2, 4) were created using MATLAB. Images were binned into 100-μm voxels and the percentage of segmented pixels within each voxel was calculated. Area of visible dot in the overlay corresponds with the averaged labeling intensity within a voxel across a condition. Outer dots represent the cell type with the highest labeling intensity within that voxel.
Computational analyses
Hierarchical clustering of regional axon quantifications was performed as described by Kebschull and colleagues (Kebschull et al., 2020).
DeepCOUNT analysis pipeline
Image segmentation and registration
Whole-brain image stacks from the 640-nm channel from each brain were segmented in TrailMap using the TRAP2-Ai14 trained model. Images were registered using the 488-nm autofluorescence channel as described above. Transformed images were converted to 8-bit in ImageJ.
3D maxima detection and single-pixel reduction
MATLAB was used to identify 3D maxima of the transformed probability map. Connected component analysis was then used to reduce any maxima that consisted of multiple pixels into a single pixel per cell.
Cell quantification
Regional TRAPed cell density was quantified in MATLAB by counting the number of labeled pixels (i.e., cells) in each brain region, then dividing this pixel count by the total number of pixels in a region. This regional pixel count was then divided by the total number of detected cells across the brain to normalize for differences in tamoxifen-induced recombination. Regions were defined in the same way as described for the DeepTraCE pipeline.
Functional network construction
In MATLAB, we calculated simple linear correlations between all pairs of brain regions based on the normalized, per area TRAPed cell counts. The nodes and connections in the networks represent brain regions and correlations that survived thresholding using Pearson's r ≥ 0.9 and p < 0.05 as cutoffs. While potentially interesting, we did not consider negative correlations in the current analysis. We used Cytoscape software to visualize and analyze networks. Degree was calculated by counting all the above-threshold connections for a given node. Node size is proportional to degree and lines represent above-threshold correlations.
4-Hydroxytamoxifen preparation
4-Hydroxytamoxifen (4-OHT; Sigma, catalog #H6278) was dissolved at 20 mg/ml in ethanol by shaking at 37°C for 15 min and was then aliquoted and stored at –20°C for up to several weeks. Before use, 4-OHT was redissolved in ethanol by shaking at 37°C for 15 min, a 1:4 mixture of castor oil:sunflower seed oil (Sigma, catalog #259853 and #S5007) was added to give a final concentration of 10 mg/ml 4-OHT, and the ethanol was evaporated by vacuum under centrifugation. The final 10 mg/ml 4-OHT solutions were always used on the day they were prepared. All injections were delivered intraperitoneally.
Behavioral assays
For PMA experiments, mice were placed in an operant chamber with a shock floor. The quarter of the floor furthest from the door of the chamber was covered in a white Plexiglas platform. Two odor pots with novel scents (vanilla, almond, coconut, or peanut butter) were placed underneath the part of the shock floor not covered by the platform to promote exploration.
For TRAP2 experiments, on training day, mice received three baseline tones (30 s, 4000 Hz, 75 dB) followed by nine tone-shock pairings (0.13-mA shock, 2 s, co-terminating), where mice could learn to avoid the shock by entering the safety platform. Tones were separated by a random interval between 80 and 150 s. The next day, mice received six tones with no shock. Mice were injected with 4-OHT solution immediately following the retrieval session. Nonshock control animals were placed in the operant chamber with no platform for 5 min on day 1, and on day 2 were placed in the operant chamber for 5 min and injected with 4-OHT solution immediately after this session. Brains were harvested two weeks after TRAPing for brain clearing, light sheet imaging, and analysis.
For fiber photometry recordings during PMA, mice received three baseline tones on day 1, followed by 12 tone-shock pairings on day 1 and 16 tone-shock pairings on days 2 and 3. For elevated zero maze (EZM) experiments, mice were placed on a custom-built elevated zero maze 24 inches in diameter for 15 min. Both assays were recorded using a Point Gray Chameleon3 USB camera (Teledyne FLIR).
Fiber photometry
Mice were habituated to the operant chamber and optic fiber for at least 2 d before recording. A Tucker-Davis Technologies (TDT) RZ10x processor in combination with the TDT Synapse software was used to simultaneously record the 405-nm isosbestic channel and the 465-nm signal channel during behavior. For each mouse, the light output was adjusted such that the 465-nm and 405-nm channel produced a signal of ∼80 mV as reported by the Synapse software.
Fiber photometry analysis
Point-tracking of PMA videos was performed in DeepLabCut (A. Mathis et al., 2018) and behavior was analyzed using BehaviorDEPOT (Gabriel et al., 2022). Elevated zero maze behavioral epochs were annotated manually. Fiber photometry analysis was performed in MATLAB in a modified version of an example provided by TDT. To align fiber photometry and behavioral data, Transistor-transistor logic (TTL) pulses marking the beginning and end of each tone were aligned between the fiber photometry signal and video frames. A lookup table was generated using linear interpolation between each TTL pulse to identify which behavior frame lines up with each photometry frame. For alignment of EZM fiber photometry data, a silent TTL pulse was generated every 30 s to be used for alignment in the same manner.
Fiber photometry signal was down-sampled by a factor of 10 before analysis. To account for potential movement artifacts and bleaching, the 405-nm isosbestic control channel was fit to the 465-nm signal using the polyfit function, and this curve was then subtracted from the 465 nm signal. Z scores of this signal were calculated using a baseline period of −10 to 0 s relative to the tone for tone-aligned responses (i.e., tone and shock responses) and −20 to −15 s relative to epoch onset for all other behaviors (i.e., platform entries and exits, closed arm entries and exits, head dips). The average of all traces for an individual animal was calculated and used for analysis. To generate plots, each animal's average trace was smoothed by averaging values from every 0.5 s (for time-locked tone and shock responses) or using a moving average of 0.5 s (for all other traces), and the mean ± SEM of smoothed traces across animals was displayed.
Basolateral amygdala (BLA) axon collaterals from photometry brains were quantified in 2D images acquired on a Leica DM6 B scanning fluorescence microscope. We took the ratio of fluorescence intensity in BLA to fluorescence intensity in a nearby background region of cortex that did not contain axons. 1 was subtracted from this value such that 0 represents the absence of axons.
Statistical analysis
Statistical analysis was performed using GraphPad Prism. For whole-brain analysis of axons, ANOVA examined the influence of neuronal class on projections patterns. We corrected for multiple comparisons using Tukey's post hoc tests. For analysis of TRAPed cells in key brain regions, we compared PMA animals to controls using Student's t tests. For fiber photometry experiments, ANOVA with Tukey's post hoc tests evaluated the influence of neuronal class on observed activity. All correlation analyses between photometry and behavior data or photometry and axon data were performed using a simple linear regression.
For whole-brain analysis of axons or TRAPed cells a single principal components analysis was performed that included all animals across groups. First, per area counts were normalized by region volume and total counts per brain (counts of axon-containing pixels or TRAPed cells, respectively). To visualize the whole-brain data in two-dimensional space, principal components analysis was applied to these normalized, per area counts to find the axes that captured the most variance across different brain areas. The per animal data across all brain regions was then projected onto the first two principal components. All computations were performed using MATLAB.
Code accessibility
Code, instructions, and sample data available at https://github.com/DeNardoLab/DeepTraCE and https://github.com/jcouto/DeepTraCE/tree/gui.
Results
DeepTraCE workflow
We developed DeepTraCE, an open-source, end-to-end analysis pipeline for quantifying bulk axonal projection patterns in cleared brains. DeepTraCE takes in raw images of fluorescently labeled axons and then applies TrailMap (Friedmann et al., 2020), a machine learning pipeline that trains a 3D U-net framework to automatically identify axons in cleared tissue (Fig. 1A). TrailMap provides a pretrained model that users can fine-tune to fit their samples. However, it was challenging to train a model that could sufficiently generalize across the diversity of mPFC axonal structures in all of their target regions; some axons were bright and delineated (e.g., in cortical areas) while others were dense and indistinct (e.g., in the basolateral amygdala and thalamus). To surmount this hurdle, DeepTraCE allows users to combine multiple trained models to optimize axon detection throughout the brain. After importing 3D image stacks, multiple models trained using TrailMap fully segment each brain. DeepTraCE then registers each brain to a standard atlas and generates a new image in which pixel values in each region are extracted from the assigned model. Axons in the resulting images are then thinned to single-pixel width (Friedmann et al., 2020) and the number of axon-containing pixels are quantified for each brain region (Fig. 1A).
DeepTraCE workflow. A, Overview of DeepTraCE workflow. B, Demonstration of segmentation using three different models followed by model combination. C, Example of human-labeled axons and artifacts for model training. D, Overlay of raw data and axon segmentation using single model or DeepTraCE concatenated models. E, Distinction score produced by models 1, 2, and 3 for images from regions assigned to axon groups 1, 2, and 3. Repeated-measures one-way ANOVA with Benjamini, Krieger, and Yekutieli post hoc test (FDR = 0.05; group 1: F(1.4,5.5) = 8.503, p = 0.0257; group 2: F(1.4,7.8) = 7.753, p = 0.0309; group 3: F(1.2,6.1) = 73.76, p = 0.0001). F, Comparison of DeepTraCE with alternate segmentation and model combination strategies. Repeated measures one-way ANOVA with Benjamini, Krieger, and Yekutieli post hoc test (FDR = 0.05; F(2,33) = 14.46, p < 0.0001). See Extended Data Figure 1-1 for detailed statistics. Scale bars, 200 μm. Descriptive statistics: mean ± SEM; *p < 0.05, **p < 0.01, ***p < 0.001. Distinction score: (LP − LN)/(LN + (1 − LP)). LP = axon+ human label. LN = axon− human label.
Extended Data Figure 1-1
ANOVA tables for regional axon density comparisons. Related to Figures 1, 3, and 4. Download Figure 1-1, DOCX file.
Validation of model combination method
To determine how many models to combine, we examined raw images of fluorescently labeled axons taken with LSFM. We grouped brain regions into three main categories based on the visual quality of PL axons (e.g., thick and delineated or dense and fuzzy). These categories largely segregated with the physical location of brain regions in tissue volumes (Fig. 1B). We used TrailMap (Friedmann et al., 2020) to train three separate models that we fine-tuned for each set of regions (Fig. 1B–D).
To validate the DeepTraCE model assignments, we performed quantitative comparisons of model accuracy in representative brain regions from each category, using human annotations as a reference. Human experts traced pixels containing axons in image stacks from each group of brain regions (Fig. 1C). We segmented images using each of the three models in addition to the original TrailMap model that was trained on serotonergic axons (Friedmann et al., 2020). We then calculated the resulting axonal labeling density with reference to human labels.
In LSFM of cleared tissue, the size of small fluorescent structures is amplified such that a thin object such as an axon occupies several more pixels in the image than it does in true biological space. We account for this overrepresentation by thinning segmented axons to a single pixel width (Friedmann et al., 2020). However, because rapid imaging of intact rodent brains requires lower optical resolution, it can be challenging in some regions to hand-label axons with single pixel precision. This caveat to whole-brain light sheet microscopy approach prevents the calculation of “true positive” and “true negative” rates based on human annotation. So, common metrics such as F1 score, precision, and recall do not provide interpretable measures of segmentation accuracy. Instead, we calculated a “distinction score” that measures the density of segmented axons in regions that humans estimated as containing axons versus not containing axons. If there is a large mismatch between human and computer labels, the distinction score will approach 0. Better overlap with human labels will produce higher distinction scores, but the upper limit is arbitrary because of the axon thinning process, which inherently reduces machine labeling even in human-annotated axon-positive pixels. We empirically observed that well-aligned labels produce values between 0.2 and 0.6 depending on the brain region. When applied to our validation set, Models 1, 2, and 3 produced the highest distinction scores in images from the brain regions assigned to each model (Fig. 1E).
Previous studies used alternative approaches to overcome the limitations of a single-model approach in whole-brain image segmentation. Some excluded deeper, more densely innervated areas from analysis (Kramer et al., 2021). Others combined probability maps from different segmentation methods by extracting the maximum probability for each pixel (Kebschull et al., 2020). To compare DeepTraCE with these alternative approaches, we calculated the distinction score across all human-annotated images in the validation set using the original TrailMap model that was trained on serotonin neurons, our best-trained single model, the maximum probability between our three models, the sum of probabilities, and the DeepTraCE combined models. DeepTraCE had the highest distinction score (Fig. 1F), indicating that it provided the most accurate segmentation of axons.
Viral circuit mapping and brain clearing strategy
We mapped the brain-wide axon collaterals of PL neurons that project to cPL, NAc, or VTA. Our goals were (1) to demonstrate the utility of DeepTraCE and (2) to provide an important anatomic resource for researchers studying mPFC function. Therefore, we used a commonly employed dual-virus approach to fluorescently label populations of PL projection neurons (Schwarz et al., 2015). In contrast to multisite retrograde tracer injections, this approach allows us to visualize the bulk brain-wide collateralization patterns of each projection-defined class of neurons. We first injected an axon-terminal-transducing adeno-associated virus (Tervo et al., 2016) expressing Cre recombinase (AAVrg-Ef1a-Cre) into cPL, NAc, or VTA. Next, into PL, we injected AAV expressing Cre-dependent EYFP-tagged Channelrhodopsin-2 (AAV8-hSyn-DIO-ChR2-EYFP), which traffics efficiently to axons (Fig. 2A). We used the iDISCO clearing variant, Adipo-Clear (Chi et al., 2018a,b) to immunostain intact brains for EYFP and render the tissue transparent. We then imaged the hemispheres ipsilateral to the PL injection with LSFM (Fig. 3).
Validation of Cre-dependent virus and injection site mapping. A, Viral injection strategy for targeting projection-defined PL neurons. B, Layer distributions of PL neurons retrogradely labeled via injections of AAVrg, AAVrg+AAV8, CTB-594, CTB-647, and AAVrg-EGFP in cPL, NAc, and VTA, respectively (cPL: Ftracer(2,6) = 4.76, p = 0.06; Fbin(2.5,15.3) = 20.04, p < 0.0001, n = 3/group; NAc: Ftracer(2,6) = 0.52, p = 0.62; Fbin(1.1,6.7) = 4.54, p = 0.07, n = 3/group; VTA: Ftracer(1,4) = 203.6, p = 0.0001; Fbin(1.7,6.8) = 4.28, p = 0.07, n = 3/group two-way ANOVA). Scale bars, 100 μm. C, Representative stain of layer marker CTIP2 in PL. Scale bars, 100 μm. D, Representative coronal sections of brains injected with AAVrg-Cre and AAV-DIO-ChR2-EYFP (left) or AAV-DIO-ChR2-EYFP in the absence of Cre and quantification of fluorescence in PL. E, Raw fluorescence images from LSFM showing mPFC injection sites for PL-cPL, PL-NAc, and PL-VTA collateral mapping experiments. Quantifications show fraction of fluorescently labeled cells detected in ACC, PL, and IL (each dot represent 1 mouse, n = 4 mice per group). F, Maps of fluorescently labeled cell bodies across the anterior-posterior axis of mPFC for each cell class. Error bars, mean ± SEM, ****p < 0.001, Student's t test. See Extended Data Figure 2-1 for quantification of layer limits.
Extended Data Figure 2-1
PL cortical layer estimates. Related to Figure 2. Download Figure 2-1, DOCX file.
Visualization and quantification of brain-wide projection patterns of mPFC. A–C, DeepTraCE segmentation of single 10-μm coronal optical sections of thinned axons registered to the standardized brain atlas. Each image represents four overlaid brains from each neuronal class. D, Dotogram overlay of three cell classes. Innvervation by cPL-projecting, NAc-projecting, and VTA-projecting PL neurons shown in yellow, blue, and purple, respectively. E, Left heatmap, Relative labeling density (normalized to region volume and gross label content per brain) across 140 regions defined by the Allen Brain Atlas. Middle heatmap, p values from multiple comparisons of axon innervation density. Right heatmap, Loadings for principal component (PC) 1 (arbitrary weight units). F, Dendrogram of hierarchical clustering of regional axon quantifications from 12 mice, colored by target region. G, Locations of individual mice projected in PC space defined by the first two PCs (arbitrary PC units, cPL, n = 4; NAc, n = 4; VTA, n = 4). See Extended Data Figures 1-1, 3-1, and 3-2 and Movie 1 for detailed statistics, abbreviations, and related data.
Extended Data Figure 3-1
List of abbreviations. Related to all figures. Download Figure 3-1, DOCX file.
Extended Data Figure 3-2
Quantification of axon density by brain region. Related to Figure 3. Download Figure 3-2, DOCX file.
3D visualization of whole-brain collateralization patterns of PL-cPL, PL-NAc, and PL-VTA neurons. Panels represent a single brain from each class.
In separate animals, we used confocal microscopy to examine the distribution of retrogradely labeled cell bodies in brain sections (Fig. 2B). Different retrograde tracers can have preferences for particular neuronal types. For instance, while AAVrg favors cortical layer (L)5 over L6, rabies virus has the opposite preference (Sun et al., 2019). Cholera toxin subunit B (CTB) is inefficient for labeling pyramidal tract (PT) neurons in Layer (L)5 (Leow et al., 2022). To assess potential biases in our retrograde labeling method, we analyzed the layer distribution for AAVrg-Cre-mCherry+ cells and the combination of AAVrg-hSyn-Cre-mCherry and AAV8-hSyn-DIO-ChR2-EYFP (assessed based on EYFP fluorescence). As a comparator, we measured the layer distribution of CTB+ cells in PL following injections into cPL and NAc. Consistent with previous reports observing other populations of PT neurons, we did not observe retrogradely labeled PL cells when we injected CTB into the VTA (Leow et al., 2022). We determined layer boundaries based on DAPI nuclear staining in combination with immunostaining for Ctip2, a marker of subcerebral projection neurons located in L5b–L6 (Arlotta et al., 2005). Layer boundaries and thicknesses aligned with previous literature (Extended Data Fig. 2-1; DeNardo et al., 2015; Anastasiades et al., 2019; Anastasiades and Carter, 2021). Consistent with previous findings (McGarry and Carter, 2016; Collins et al., 2018; Anastasiades et al., 2019), PL-cPL and PL-NAc neurons were preferentially located in L2–5a and PL-VTA neurons were restricted to L5b. We found no statistical differences in the distribution of cells labeled with AAVrg alone, AAVrg+AAV8 or CTB (Fig. 2B,C). As reported previously, we observed few AAVrg-labeled cell bodies in L6 (Sun et al., 2019; Leow et al., 2022). While L6 neurons were likely underrepresented in our datasets, we could still extract meaningful differences between classes. To ensure viral labeling was specific, we confirmed there was no cre-independent expression of our fluorescent constructs (Fig. 2D) and that injection sites were centered on PL. Fluorescent labeled neurons were mostly in PL, with some in the anterior cingulate cortex (ACC; Fig. 2E,F).
Whole-brain projection patterns of PL-cPL, PL-NAc, and PL-VTA neurons
We used DeepTraCE to determine the extent to which PL-cPL, PL-NAc, and PL-VTA neurons collateralize to other brain regions. Examination of thin optical sections revealed widespread collateralization for all three classes. Overall, innervation patterns were consistent with reports of general mPFC projection patterns (Gabbott et al., 2005; Anastasiades et al., 2019; DeNardo et al., 2019); mPFC axons were prominent in cortical association areas, striatum, midline thalamus, claustrum, amygdala, hypothalamus, and midbrain but largely absent from hippocampus, sensory thalamus, and relatively sparse in primary sensory areas (Fig. 3A–C; Movie 1).
We generated dotograms (Oh et al., 2014) to visually summarize the interclass distinctions. To create the dotogram, we averaged axon density across all samples for a given projection class and then overlaid the samples in brain space. Dot color indicates the neuronal class and dot size corresponds to the averaged axonal density within a voxel (Fig. 3D). Compared with PL-VTA neurons, PL-cPL and PL-NAc neurons more densely innervated cortical areas, especially the temporal association area (TEa), ectorhinal area, and entorhinal cortex. On the other hand, PL-VTA collaterals preferentially innervated thalamic (TH) and midbrain (MB) regions, a subset of which also received collateral input from PL-NAc neurons. All three subclasses innervated NAc and olfactory tubercle (OT), while PL-NAc neurons were the primary source of collaterals to piriform cortex (PIR).
We next quantified regional innervation densities and plotted them as a heatmap. For each brain, we normalized the number of axon-containing pixels in each region to the total number of axon-containing pixels in the brain. We then sorted regions according to their innervation density by projection classes; regions receiving more PL-cPL collaterals are on top and those receiving more PL-VTA collaterals are in the bottom half of the heatmap (Fig. 3E). While there was little overlap between projection patterns of PL-cPL and PL-VTA neurons, PL-NAc neurons shared several projection targets with both other classes. Brain-wide statistical comparisons confirmed the presence of 29 significantly differentially innervated brain regions, the most notable of which were in the cortex, subplate, and thalamus (Fig. 3E; Extended Data Figs. 1-1, 3-2).
Because PL-cPL, PL-NAc, and PL-VTA projection neurons all collateralize broadly, we wondered whether these neuronal classes could be separated based on their whole-brain projection patterns. To test this, we used two dimensionality reduction techniques. We first performed hierarchical clustering, in which Euclidean distance between individual brains based on axonal density in each subregion is marked by a higher branch point on the graph. PL-VTA neurons formed a distinct cluster, while PL-NAc and PL-cPL collaterals formed partially overlapping clusters (Fig. 3F).
We also used principal component analysis to assess differences in whole-brain projection patterns for each cell class. We plotted the location of each brain along the axes of the first two principal components. While PL-cPL and PL-VTA projection classes could be clearly separated along the first principal component, PL-NAc brains were positioned in between, with some overlap with PL-cPL brains (Fig. 3G). To assess which projection patterns distinguished the classes, we plotted the weights of each brain region in contribution to the first principal component. Positively weighted regions generally had the most innervation from PL-cPL collaterals including TEa, PERI, and ECT, and negatively weighted regions were those with most innervation from PL-VTA collaterals, including MTN and midbrain areas (Fig. 3E). These data show that while cPL and NAc-projectors have substantial overlap, PL-cPL and PL-VTA neurons represent separable classes, distinguished in large part by their collaterals to TEa/PERI/ECT, olfactory and limbic areas, and thalamus and midbrain, respectively. Together, the innervation patterns we observed were consistent with bulk projection mapping (Sesack et al., 1989; Vertes, 2004; Oh et al., 2014) and single neuron reconstructions (Winnubst et al., 2019; Gao et al., 2022), indicating that DeepTraCE accurately captured meaningful class-specific distinctions in collateral targeting.
Layer-specific and topographic innervation patterns in the cortex and subplate
Whole-brain analysis of bulk labeled neuronal classes is well-suited for understanding the geometric organization of axonal projections within brain regions, which can have important functional implications. The neocortex is organized into layers that contain distinct neuronal types with different morphology, physiology, and connectivity (Harris and Shepherd, 2015; Anastasiades and Carter, 2021). Long-range corticocortical axons may innervate superficial or deep layers depending on the hierarchical relationship with the target region (Harris and Shepherd, 2015; Anastasiades and Carter, 2021). Compared with sensory and motor cortices, much less is known about the layer organization of long-range mPFC connectivity, especially for specific projection classes. To investigate this, we plotted PL-cPL, PL-NAc, and PL-VTA collateral innervation density across layers in cortical target areas (Fig. 4A,B). In ORB, all classes preferentially targeted superficial layers (L1, L2/3). In TEa, ECT, and PERI, PL-cPL and PL-NAc collaterals preferentially targeted superficial layers (L1, L2/3) while PL-VTA collaterals targeted the deep layers (L5, L6; Fig. 4B; Extended Data Fig. 1-1). These findings suggest that PL-cPL, PL-NAc, and PL-VTA may all participate in feedforward or feedback connectivity depending on the target region.
Region-specific collateralization patterns of PL-cPL, PL-NAc, and PL-VTA neurons. A, Coronal view of 100-μm z-projections of raw 640-nm (axons, white) and 488-nm (autofluorescence, blue) channels from individual brains. Images show class-specific innervation patterns in anterior (left) and posterior (right) cortical areas. B, Layer distributions of axonal innervation in select cortical regions. C, Visualizations of axonal innervation density along the anterior-posterior axis in select regions. D, Raw images showing axons in the amygdalar complex. E, Quantification of axonal innervation density along the anterior-posterior axis of a given region. PL-cPL, PL-NAc, and PL-VTA neurons coded in yellow, blue, and purple, respectively. N = 4/group. Descriptive statistics are from two-way ANOVA with Tukey's post hoc multiple comparison test. See Extended Data Figures 1-1 and 3-1 for detailed statistics and abbreviations. Error bars: mean ± SEM; *p < 0.05, **p < 0.01. Scale bars, 200 µm.
Many PL target regions contain functional gradients that are orthogonal to layer-defined microcircuits. We therefore investigated whether PL classes project to topographically defined locations in target regions (Fig. 4C). Compared with PL-cPL and PL-NAc, PL-VTA collaterals had a tighter distribution in the anterior portion of agranular insula (AI) and claustrum (CLA). Compared with PL-VTA neurons, PL-cPL and PL-NAc collaterals were biased toward posterior visual cortex. PL-NAc collaterals were especially prominent in the posterior PIR. All three classes had similar distributions across entorhinal cortex, preferentially targeting the posterior part. Interestingly, PIR and AI both contain anatomic and functional gradients (Gehrlach et al., 2020; Poo et al., 2022). Our data suggest that class-specific PL projections may be a determining factor in these gradients.
While all projection classes sent collaterals to the amygdalar complex, amygdalar nuclei were differentially innervated. Notably, the anterior basolateral amygdala (BLAa) was robustly innervated by PL-cPL and PL-NAc neurons, but less so by PL-VTA neurons. Meanwhile, compared with PL-cPL neurons, PL-NAc neurons sent many more collaterals to olfactory amygdalar areas, including anterior cortical amygdalar area (COAa) and piriform amygdalar area (PAA; Fig. 4D,E). As PL-NAc collaterals were also more prominent in posterior PIR, which was recently found to play a role in spatial cognition (Poo et al., 2022), PL-NAc neurons may have privileged control over areas dedicated to olfactory processing during cognition and emotional learning.
Whole-brain cell counting pipeline
The application of computer vision to light sheet data extends beyond axon tracing and can be used to identify signatures of brain function. Expression of IEGs such as c-fos are frequently used as a proxy of neural activity (Sheng and Greenberg, 1990; Franceschini et al., 2020). For instance, recent studies have used tissue clearing, IEG-labeling, and LSFM to screen for changes in neuronal activation following experiences with drugs or fear learning (Wheeler et al., 2013; Renier et al., 2016; Vetere et al., 2017; DeNardo et al., 2019). Quantifying IEG-expressing cells on a brain-wide scale can serve as means to infer behaviorally relevant changes in functional connectivity. However, currently available open-source packages for whole-brain cellular quantification can be error prone, especially in regions of particularly dense labeling.
As a companion to DeepTraCE, we developed DeepCOUNT (deep-learning based counting of objects via 3D U-net pixel tagging; Fig. 5A). We first trained a TrailMap model to recognize fluorescently labeled cell bodies. This produces a probability map of which pixels are most likely to contain a cell body. After thresholding the images to a desired probability cutoff, they are registered to a standard brain, and DeepCOUNT uses a 3D maxima detection strategy in combination with a connected component analysis to ensure each neuron is represented by a single pixel. This single-pixel output can then be used to obtain regional cell counts across the brain.
Whole-brain IEG mapping of learned avoidance with DeepCOUNT. A, Overview of DeepCOUNT workflow. B, Overlay of raw data and cell segmentation using trained cell detection model. C, Comparison of residual error values for DeepCOUNT and ClearMap as compared with human annotations (p = 0.022, n = 6 brain regions; paired t test). D, Overview of platform-mediated avoidance assay. TRAP2;Ai14 mice were trained in PMA on day (D)1 and then TRAPed during a retrieval test on D2. E, Fraction of successful trials increased (F = 5.67, p = 0.012, n = 4, one-way ANOVA) and latency to enter the platform decreased (F = 7.21, p = 0.005, n = 4, one-way ANOVA) across training. On retrieval day, mice entered the platform with an average latency of 7.9 s after tone onset and were more likely to be on the platform at the end of the tone compared with pretraining (p = 0.0005, n = 4, paired t test). F, Left heatmap, Relative cell density (normalized to region volume and gross label content per brain) across 140 regions defined by the Allen Brain Atlas. Right heatmap, Loadings for PC1 (arbitrary PC weight units). G, Dotogram overlay of control (gray) and PMA (green) conditions. Dots represent the density of TRAPed cells in a given voxel. H, Locations of individual mice projected in principal component (PC) space defined by the first two PCs (arbitrary PC units, control, n = 5; PMA, n = 4). I, Comparison of TRAPed cell density in five brain regions of naive controls (gray) and PMA-trained mice (green; Student's t test; *p < 0.05, **p < 0.01; control: n = 5; PMA: n = 4). J, Network diagrams for control (left) and PMA-trained (right) mice based on brain-wide interregional correlations. Node size is proportional to degree and lines represent above-threshold correlations. K, Degree values for top 20 most connected regions for each condition (PL outlined in red). See Extended Data Figure 5-1 for detailed statistics and abbreviations. Error bars, mean ± SEM.
Extended Data Figure 5-1
Quantification of TRAPed cell density by brain region. Related to Figure 5. Download Figure 5-1, DOCX file.
We compared DeepCOUNT performance to that of ClearMap (Renier et al., 2016), a commonly used cell detection algorithm designed for tissue clearing and LSFM, and human observers. Two human observers analyzed image volumes from six different brain regions. Cell counts produced with DeepCOUNT were more accurate counts than ClearMap, using human counts as the ground truth (Fig. 5B,C).
Whole-brain threat avoidance networks
We used DeepCOUNT to map regions and functional connections involved in threat avoidance at a whole-brain level. PL is required for threat avoidance behavior (Bravo-Rivera et al., 2014; Diehl et al., 2018). Previous studies examining c-fos expression following active avoidance behavior analyzed only discrete regions of interest (Nikolaev et al., 1992; Duncan et al., 1996; Bravo-Rivera et al., 2015), some of which are PL target sites. How PL coordinates threat avoidance-related neural activity on a brain-wide scale is not understood. To begin to address this question, we used DeepCOUNT and activity-dependent genetic labeling to map brain-wide neuronal activation following platform-mediated avoidance (PMA; Bravo-Rivera et al., 2014), a form of learned threat avoidance.
We labeled activated neuronal populations using TRAP2 mice (Fosicre-2A-ERT2), in which the c-fos promoter drives expression of tamoxifen-inducible Cre recombinase. We crossed TRAP2 to the Ai14 Cre reporter line (Madisen et al., 2010; Allen et al., 2017; DeNardo et al., 2019). We then trained the TRAP2;Ai14 mice in PMA and TRAPed them after a threat avoidance retrieval session by pairing that experience with a tamoxifen injection (Fig. 5D). In PMA, we defined successful trials as those in which the animal remained on the platform for the last 2 s of the tone, when the shock normally occurs. During training, mice significantly increased their fraction of successful trials and decreased their latency to enter the safety platform. On the retrieval day, mice entered the platform with an average latency of 7.9 s after tone onset and were more likely to remain on the platform at the end of the tone compared with their pretraining baseline (Fig. 5E).
We used DeepCOUNT to quantify TRAPed cell density in each brain region and plotted these data as a heatmap in which regions with more TRAPed cells in PMA-trained mice are on top and regions with more TRAPed cells in naive controls are in the bottom half of the heatmap (Fig. 5F). We also plotted the data as a dotogram (Fig. 5G). Principal component analysis of whole-brain TRAPed cell counts separated PMA-trained and nonshocked control animals along the first principal component (PC1; Fig. 5H; Extended Data Fig. 5-1). PC1 weights revealed that PMA-activated brain regions were biased toward the cortex, hippocampal formation, and cortical subplate (Fig. 5F). We compared cell counts between PMA and control animals in PL and several key target regions, finding significant increases in the number of TRAPed cells in PL, BLA, and TEa but not in NAc or VTA (Fig. 5I). Overall, our DeepCOUNT screen identified regions known to be involved in threat avoidance (e.g., PL and BLA). We also identified several regions with unknown roles in threat avoidance [e.g., TEa and postsubiculum (POST) that will be of interest for future functional studies; Fig. 5F,G,I].
To relate the organization of PL circuits to their role in threat avoidance, we assessed functional connectivity based on whole-brain TRAPing. We computed inter-regional correlations for groups of PMA-trained and naive control mice. This allowed us to identify sets of regions where the numbers of TRAPed cells co-varied across mice. Regions that co-vary may constitute elements of a network engaged during learned threat avoidance. We then generated network graphs in which the nodes and connections represent brain regions and correlations that survived thresholding using r ≥ 0.9 and p < 0.05 as cutoffs. In the network plots, node size is proportional to degree, or the number of connections for that node (Fig. 5J).
To link these findings to our neuroanatomical studies, we examined PL functional connectivity while considering the collateral maps presented in Figures 3 and 4. Our network analysis revealed PL as one of the mostly highly connected nodes in the PMA group (Fig. 5K). PL was functionally connected to several cortical areas (TEa, ORBl, MOp, SSp, SSs), taenia tecta (TT), the lateral septal complex (LSx), claustrum (CLA), several amygdalar regions (LA, BLAa, BMAa), and the median eminence (ME). Of these regions, CLA, TEa, TT are preferentially innervated by PL-cPL collaterals and LSx is preferentially innervated by PL-VTA collaterals. Compared with PL-VTA, PL-NAc, and PL-cPL neurons send denser projections to BLAa. Compared with PL-cPL neurons, PL-NAc and PL-VTA neurons send denser projections to BMAa (Figs. 3, 4). To visualize this preferential anatomic innervation on the functional connectivity maps, we colored the PL-connected nodes based on their PCA weight from the analysis shown in Figure 3E (Fig. 5J). These data suggest that all three neuronal classes likely contribute to threat avoidance behavior and that they may act via particular collateral projections.
Fiber photometry during platform-mediated avoidance
As prelimbic cortex has been implicated in both innate and learned avoidance, we investigated how different PL cell classes contribute to threat avoidance by using fiber photometry to record population-level neural activity. Activity in PL-NAc neurons has previously been shown to reduce learned avoidance responses (Diehl et al., 2020). PL-VTA and PL-cPL neurons have yet to be studied in the context of learned avoidance; however, PL-VTA neurons encode aversive stimuli, and PL-cPL neurons are a strong candidate to mediate the prefrontal hemispheric synchrony observed in spatial avoidance (Vander Weele et al., 2018; Dickson et al., 2022). Further, all three of the populations we mapped have axon collaterals in BLA, a key region in fear and avoidance learning (Arruda-Carvalho and Clem, 2014; Diehl et al., 2020). In addition to revealing insights into the functions of understudied neuronal classes, our populations of choice allow us to dissect the heterogeneity of PL-BLA neurons and their role in avoidance.
Given the high overlap of collateral projection targets for PL-cPL and PL-NAc neurons (Fig. 3I,J), we further separated these classes using an intersectional viral strategy (Fenno et al., 2020). To record neuronal activity in PL-cPL neurons that do not project to NAc, we injected AAVrg-Cre into cPL and AAVrg-Flp into NAc. We then injected Cre-On Flp-Off GCaMP6M into PL. We switched the Cre and Flp injection sites to record from NAc-projecting neurons that do not project to cPL. To record from VTA-projecting neurons, we injected AAVrg-Cre into VTA and DIO-GCaMP6M into PL (Fig. 6A,B). Most GCaMP+ cell bodies from the intersectionally defined cPL-projecting and NAc-projecting neurons were in superficial layers of PL, while cell bodies from VTA-projecting neurons were restricted to the deeper layers (Fig. 6C), as expected from our initial layer analysis (Fig. 2B). We observed few GCaMP+ axon terminals in NAc from the cPL-projecting group, and few GCaMP+ axon terminals in cPL from the NAc-projecting group (Fig. 6D), suggesting our intersectional strategy was effective in separating these populations.
Neuronal class-specific activity during threat cues. A, AAV injection strategy. Left, To target cPL-projecting neurons, mice were injected with AAVrg-Cre in cPL, AAVrg-Flp in NAc, and Con-Foff-GCaMP6M in PL. Middle, To target NAc-projecting neurons, mice were injected with AAVrg-Flp in cPL, AAVrg-Cre in NAc, and Con-Foff-GCaMP6M in PL. Right, To target VTA-projecting neurons, mice were injected with AAVrg-Cre in VTA and DIO-GCaMP6M in PL. B, Representative images of GCaMP expression and fiber placement sites in PL. Scale bar, 1 mm. C, Distribution of cPL-projecting, NAc-projecting, and VTA-projecting GCaMP-expressing cells across cortical layers (Fbin(4,54) = 28.3, p < 0.0001, Fbinxclass(36,234) = 34.09, p < 0.0001; cPL: n = 6, NAc: n = 5, VTA: n = 6; two-way ANOVA). D, Representative images of axon terminals in cPL and NAc using the intersectional viral targeting strategy shown in A. Scale bars, 100 μm. E, Schematic of PMA assay and behavioral performance across sessions. Fraction successful trials: Fday(2,24) = 65.18, p < 0.0001, Fclass(2,14) = 0.58, p = 0.57, Fclass x day(6,36) = 0.43, p = 0.86, two-way ANOVA. Latency to enter platform: Fday(3,36) = 14.27, p < 0.0001, Fclass(2,14) = 0.09,038, p = 0.91, Fclass x day(6,36) = 0.48, p = 0.82, two-way ANOVA. F, GCaMP fluorescence in cPL-projecting, NAc-projecting, and VTA-projecting neurons. Signals are aligned to tone onset and separated by whether mouse is on (colored trace) or off (gray trace) the safety platform at the start of the tone. *p < 0.05 for Student's t test comparing on versus off platform activity in a given time window. G, AUC analysis of Ca2+ signal for tone periods (0–10 s) when mice were off the platform (F(2,14) = 4.106, p = 0.04, cPL: n = 6, NAc: n = 5, VTA: n = 6; one-way ANOVA with Tukey's post hoc test. H, Correlation between AUC during off-platform tone periods (0–10 s) and average latency for a mouse to enter the platform. Error bars, mean ± SEM.
We recorded from PL while mice performed PMA. We kept the shock active across 3 d of recording, which allowed us to average signals over many trials without extinction occurring. Mice learned quickly, successfully avoiding most shocks and entering the platform with an average latency of 8 s after tone onset on the first day of training, with minor performance improvements across days (Fig. 6E). We first analyzed responses to the conditioned tone. Interestingly, tone-evoked neural activity distinguished between epochs when mice were off versus on the safety platform (Fig. 6F). cPL-projecting neurons had significantly higher activity at the tone onset when the animal was off the platform. As the tone progressed, activity in these neurons gradually decreased. We observed a similar trend in NAc-projecting neurons. In contrast, PL-VTA cells lacked a tone onset response, but did have a gradual decrease in activity during the conditioned tone (Fig. 6G). These findings suggest neural activity in response to a conditioned stimulus is modulated by whether threat requires action, and that these patterns vary across projection classes.
To determine whether the observed tone onset responses related to avoidance behavior, we calculated correlations between the AUC for the first 10 s of the tone and the average latency to enter the platform. Activity in the PL-cPL neurons had a tight negative correlation with the latency to enter the platform. Activity in the PL-NAc neurons also had a negative correlation with latency to enter the platform, but it did not reach significance. Activity in PL-VTA neurons had no relationship with the latency to enter the platform (Fig. 6H). Together, these findings point to PL-cPL projectors as a key neuronal class that links threatening cues to avoidance behavior.
Activity in all classes sharply increased upon foot shock onset (Fig. 7A). NAc-projecting neurons had significantly higher activity than the other two classes in the 10-s preceding a shock (Fig. 7B), consistent with a larger role for those cells in risky exploration. We found no differences in peak shock response between the three classes (Fig. 7C). Further, while activity in PL-VTA neurons rapidly returned to baseline levels following the foot shock, activity in PL-cPL and PL-NAc neurons had prolonged deviations from baseline for up to a minute after the shock. Activity in NAc-projecting neurons remained elevated above baseline for up to 80 s following a foot shock, while activity in cPL-projecting neurons decreased steadily during the same period (Fig. 7D,E). These data suggest that while all PL classes encode aversive stimuli, prolonged postshock activity in PL-cPL and PL-NAc neurons may play specialized roles in action-outcome learning.
Neuronal class-specific activity during aversive stimuli. A, GCaMP fluorescence in cPL-projecting, NAc-projecting, and VTA-projecting PL cells. Signals are aligned to shock onset and separated by whether mouse is on (colored trace) or off (gray trace) the safety platform during the shock period. *p < 0.05, **p < 0.01 for Student's t test comparing on versus off platform activity in a given time window. Analysis of Ca2+ signals for circa-shock periods: (B) preshock AUC: F(2,14) = 6.97, p = 0.0079; (C) shock peak Z-score: F(2,14) = 1.87, p = 0.19; (D) postshock AUC (10–40 s): F(2,14) = 8.28, p = 0.0042; (E) postshock AUC (60–80 s): F(2,14) = 8.79, p = 0.0034. For B–E, cPL: n = 6, NAc: n = 5, VTA: n = 6, one-way ANOVA with Tukey's post hoc test, *p < 0.05, **p < 0.01. Error bars, mean ± SEM.
Prefrontal neuron classes distinguish between learned and innate avoidance
To further classify the effects of aversive learning on neural activity in PL cell classes, we compared their activity during learned versus innate threat avoidance. Using the same mice from the PMA experiments, we recorded neural activity in the elevated zero maze (EZM), in which two quarters of an elevated ring are protected by walls and the other two quarters are open. Mice innately avoid the open arms of this apparatus, where there is more perceived threat potential.
We compared activity during entries and exits from the safe zone in each assay. In PMA, population activity in all three neuronal classes increased before platform entries, an effect that grew stronger with training (Fig. 8A). In contrast, only NAc-projecting and cPL-projecting neurons increased their activity before entry into the closed arm of the EZM (Fig. 8B). During exits from the safe zone, a form of risky exploration, all three PL classes had similar increases in activity during PMA (Fig. 8C). In innate avoidance, however, exits from the closed arm were marked by a significantly larger response in NAc-projecting neurons compared with the other two cell populations (Fig. 8D). Together, these data indicate that aversive learning enhances encoding of avoidance behaviors and engages VTA-projecting neurons to behavioral circuits for threat avoidance. Further, PL-NAc neurons preferentially encode risky exploration.
Neuronal class-specific activity during learned versus innate threat avoidance behavior. A, B, GCaMP fluorescence aligned to safe-zone entry in learned (A) versus innate (B) avoidance for cPL-projecting (yellow), NAc-projecting (blue), and VTA-projecting (purple) neurons. Inset plots show AUC analysis of Ca2+ signal aligned to avoidance (−5–5 s). A, Platform entry AUC increases between training day (D)1 and D2/3 across cell types (PMA: Fday(1,14) = 17.59, p = 0.0009, Fclass(2,14) = 1.488, p = 0.26, Fclass x day(2,14) = 0.98, p = 0.40, two-way ANOVA). B, EZM: F(2,9) = 6.11, p = 0.02; one-way ANOVA with Tukey's post hoc test. C, D, GCaMP fluorescence aligned to onset of risky exploration in learned (C) versus innate (D) avoidance for cPL-projecting (yellow), NAc-projecting (blue), and VTA-projecting (purple) neurons. Inset plots show AUC analysis of Ca2+ signal aligned to exploration onset (0–20 s; PMA: F(2,14) = 2.43, p = 0.12; EPM: F(2,9) = 11.97, p = 0.0029, one-way ANOVA with Tukey's post hoc test). E, Analysis of Ca2+ signals during head dips (F(2,10) = 0.64, p = 0.55; one-way ANOVA with Tukey's post hoc test). **p < 0.01, ***p < 0.001. Error bars, mean ± SEM PMA: cPL n = 6, NAc n = 5, VTA n = 6; EZM: cPL n = 5, NAc n = 4, VTA n = 3.
A concern with comparing fiber photometry recordings across different classes of neurons is that differences in signal between the classes could lead to exaggerated results in one direction. While PL-VTA neurons had lower overall activity during risky exploration of the shock bars and open arm, all three PL classes had similar activity during head dips on the EZM, another form of exploratory behavior (Fig. 8E). This suggests that PL-VTA neurons did not simply have lower activity than the other classes during behavior, but instead that these differences are specific to situations when animals are navigating the environment.
PL-cPL and PL-NAc activity is correlated with BLA collateral axon density
PL-cPL and PL-NAc neuronal classes both send dense collateral projections to BLA (Figs. 3E, 4E), a key region in fear and avoidance learning and memory (Ledoux, 2000). However, recent single neuron reconstructions revealed that PL-cPL and PL-NAc classes are heterogeneous and contain neurons that do not send collaterals to the BLA (Gao et al., 2022). We took advantage of this heterogeneity to estimate the activity of BLA collaterals, calculating correlations between the activity measured using fiber photometry and the density of GCaMP+ axons in the BLA. To do this, after completing the fiber photometry experiments, we perfused the mice and quantified axonal GCaMP fluorescence in the BLA in 2D brain sections (Fig. 9A).
Correlation of neural activity and BLA axon collateralization in PL-cPL and PL-NAc populations. A, Representative images of GCaMP+ axons in BLA of animals with low and high off-platform tone responses. B, BLA axon fluorescence intensity is correlated with off-platform tone response (AUC: 0–4 s). cPL n = 6, NAc n = 5. Simple linear regression. C, BLA axon fluorescence intensity is not correlated with platform entry responses on day 2/3 (AUC: −5–5 s; PL-NAc: R2 = 0.38, p = 0.27; PL-cPL: R2 = 0.003, p = 0.92) or platform exit responses (AUC: 0–20 s; PL-NAc: R2 = 0.02, p = 0.82; PL-cPL: R2 = 0.35, p = 0.21). cPL n = 6, NAc n = 5. Simple linear regression.
In both PL-cPL and PL-NAc classes, we observed a significant positive correlation between tone-onset responses (AUC during first 4 s of the tone when mice were off the platform) and the density of GCaMP+ axon collaterals in BLA (Fig. 9B). On the other hand, we observed no significant relationships between activity and platform entries or exits (Fig. 9C). These data suggest that within the PL-NAc and PL-cPL populations, those neurons that collateralize to the BLA drive responses to threatening cues. Importantly, since we used intersectional approaches to separate PL-cPL and PL-NAc neurons, these data also indicate that multiple subclasses of PL-BLA projection neurons encode conditioned cues. Interestingly, for PL-NAc neurons, there was a trend-level correlation between activity during platform entries and BLA collateral axon density. This suggests there may be a subclass of PL neurons that projects to both NAc and BLA and encodes the behavioral action of threat avoidance, but more investigation is needed to confirm this.
Discussion
In this study, we integrated whole-brain mapping with the observation of neural activity to better understand how mPFC controls threat avoidance behavior. We introduce DeepTraCE and DeepCOUNT, open-source analysis pipelines for quantifying bulk axonal projection patterns and cells, respectively. We used DeepTraCE to produce whole brain projection maps of three populations of mPFC projection neurons: PL-cPL, PL-NAc, and PL-VTA. We combined activity-dependent genetic labeling with DeepCOUNT to relate the structural and functional organization of brain-wide networks for threat avoidance. With whole-brain projection maps as a foundation, we used intersectional viral targeting to separate the overlapping populations of PL-cPL and PL-NAc neurons. We then recorded the activity of PL-cPL, PL-NAc, and PL-VTA neuronal populations during learned and innate threat avoidance. Our study reveals PL class-specific roles in threat avoidance and demonstrates the utility of DeepTraCE and DeepCOUNT for linking high-throughput neuroanatomy with functional techniques to reveal mechanisms of complex behavior (Fig. 10).
Summary of findings. A, Visualization of collateralization density in key targets of PL-cPL, PL-NAc, and PL-VTA neurons. Dot radius correlates with average normalized labeling density within a region. B, Schematic of whole-brain collateralization patterns of PL-cPL, PL-NAc, and PL-VTA neurons. C, Summary of PL efferent connectivity patterns. D, Visualization of activity levels in PL-cPL, PL-NAc, and PL-VTA neurons during aspects of threat avoidance. Dot radius correlates with average normalized signal intensity. See Extended Data Figure 3-1 for abbreviations.
While we have gained invaluable insights into the complex anatomy of cortical neurons from single neuron reconstructions (Oh et al., 2014; Economo et al., 2018; Gao et al., 2022), the equipment, expertise, and time required to reconstruct individual neurons are out of reach for many labs. Other high-throughput techniques such as MAPseq (Kebschull et al., 2016) provide single cell resolution, but are costly and limited to predefined target regions. While bulk labeling approaches lose granularity, DeepTraCE's high-throughput workflow combined with its reliance on common viral techniques make it a useful tool to tackle many unanswered questions. Mouse Cre-driver lines and floxed alleles can be used with DeepTraCE to analyze the projection patterns of defined neuronal populations and how they can be altered by genetic manipulations. Or DeepTraCE and DeepCOUNT can be used to define differences in structural and functional connectivity resulting from environmental insults such as chronic stress. With DeepTraCE, in addition to simply quantifying relative densities of axonal projections in target regions, users can assess within-region changes in layer-specific or topographic targeting. This information can guide the design of precisely targeted functional studies.
There are caveats associated with retrograde viruses, brain clearing, and LSFM. AAVrg viruses have been shown to have preferential labeling of corticothalamic (CT) neurons in L5, and biases against CT neurons in L6 (Sun et al., 2019; Leow et al., 2022), while CTB showed the opposite pattern of layer bias (Leow et al., 2022). We observed few AAVrg+ neurons in L6, but some CTB+ neurons (Fig. 2B) and may have underestimated the contributions of L6 cells to the population data for PL-cPL and PL-NAc neurons (Gabbott et al., 2005; Anastasiades and Carter, 2021). Also, tissue clearing does not always achieve complete transparency. This issue, together with resolution limits of light sheet microscopes that can visualize intact mouse brains in a timely fashion, can produce blurring of structures deep in the tissue. While estimates of axonal density may be less accurate in deep brain structures, they are likely to be biased to a similar degree for different neuronal classes so are still well-suited for comparing interclass differences. Finally, like most mapping studies, we visualized axons but not synapses. Therefore, our DeepTraCE quantifications may also include axons of passage. To mitigate this, our models were trained to exclude major white matter tracts.
Despite these caveats, our brain-wide projection maps of PL-cPL, PL-VTA, and PL-NAc neurons aligned with known projection patterns (Sesack et al., 1989; Vertes, 2004; Gabbott et al., 2005; Oh et al., 2014; DeNardo et al., 2019; Winnubst et al., 2019; Gao et al., 2022), underscoring the accuracy of TrailMap (Friedmann et al., 2020) and DeepTraCE, and revealed unappreciated topographic innervation patterns. Consistent with bulk tracing studies (Sesack et al., 1989; Vertes, 2004; DeNardo et al., 2019), we observed the densest PL projections in association cortices, claustrum, striatum, amygdala, polymodal thalamic nuclei, hypothalamus, and midbrain. These bulk projection patterns are the sum of multiple classes.
mPFC contains intratelencephalic (IT) neurons that project within the cortex and PT neurons that project to subcortical areas. PL-cPL neurons fall in the IT class, PL-VTA neurons are from the PT class, and PL-NAc contain members of both (Anastasiades and Carter, 2021). In line with these classifications and single neuron reconstructions (Winnubst et al., 2019; Anastasiades and Carter, 2021; Gao et al., 2022), we observed that PL-cPL neurons were enriched in superficial layers (L2–L5a) and collateralize to other cortical association areas, subplate, and striatum. Consistent with the well-described features of PT neurons (Economo et al., 2018; Anastasiades and Carter, 2021; Gao et al., 2022), we observed retrogradely labeled PL-VTA neurons that were localized to L5b and collateralized to the thalamus, pons, striatum, and midbrain. PL-NAc neurons collateralized broadly to both IT and PT targets (Anastasiades and Carter, 2021). The diversity of the collateral targets of PL-NAc neurons is supported by single neuron reconstructions that revealed a variety of subclasses targeting only a few structures (Gao et al., 2022). Overall, this excellent alignment between our data and previous work supports the utility of DeepTraCE for capturing meaningful class-specific differences in connectivity in a high-throughput manner.
On a brain-wide scale, counts of IEG-expressing cells can reveal behaviorally relevant activity and functional connectivity (Wheeler et al., 2013; Vetere et al., 2017; DeNardo et al., 2019). Previous studies examining c-fos expression following active avoidance behavior analyzed only discrete regions of interest (Nikolaev et al., 1992; Duncan et al., 1996; Bravo-Rivera et al., 2015). Using DeepCOUNT, we recapitulated previous findings of activation in PL, BLA and OFC following PMA and found high numbers of TRAPed cells in understudied regions including POST and TEa. POST is important for spatial learning, and PL neurons that project there may facilitate learning of the location of the safety platform (Bett et al., 2012, 2013; Peckford et al., 2014). PL projections to TEa may facilitate long-term storage of the avoidance memory, as TEa is necessary for long-term retrieval of cued fear (Sacco and Sacchetti, 2010). Functional connectivity analysis revealed the major functional targets of PL. Of these, TEa, CLA, and TT are highly innervated by the collaterals of PL-cPL projection neurons, while LSx and BMAa are more densely innervated by PL-NAc and PL-VTA collaterals. These represent new targets to manipulate in future studies of threat avoidance and highlight how whole brain anatomic and functional maps can be used together to identify novel behaviorally relevant pathways.
We used fiber photometry to determine how each class of neuron encodes threat avoidance. We found that intersectionally defined PL-NAc and PL-cPL neurons increased their activity during conditioned tones, but only when mice were off the safety platform, suggesting those populations encode the predictability of threat. In contrast, PL-VTA neurons decreased their activity in response to the threat predictive cue. Single unit recordings have shown that distinct populations are excited and inhibited by a conditioned tone during PMA (Diehl et al., 2018). Our results indicate that the projection targets we studied may be a factor that separates these populations.
We found that PL-NAc neurons encoded risky exploration. This is consistent with previous work showing that excitation of PL-NAc projections decreases avoidance behaviors (Diehl et al., 2020). However, individual neurons within this projection class have heterogeneous functions. A subset of mPFC-NAc neurons restrain reward seeking under threat of punishment (C.K. Kim et al., 2017). Other work suggests mPFC-NAc neurons that suppress reward seeking preferentially target the NAc shell (Piantadosi et al., 2020). Our AAVrg injections targeted the NAc core and we likely captured a larger proportion of mPFC-NAc neurons involved in risk engagement rather than suppression. Our observation of PL-NAc neurons resembles the activity patterns seen in mPFC neurons projecting to dorsomedial striatum during avoidance (Loewke et al., 2021; Kajs et al., 2022). Indeed, we observed mPFC-NAc collaterals in the dorsomedial striatum, suggesting these populations overlap.
All three projection classes studied send axon collaterals to the BLA. Bulk stimulation of mPFC-BLA neurons promote avoidance in PMA by reducing the latency to enter the platform following tone onset (Diehl et al., 2020). Further, auditory fear conditioning induces synaptic plasticity in the mPFC-BLA pathway (Arruda-Carvalho and Clem, 2014). We found that within the PL-cPL and PL-NAc classes, the magnitude of the tone onset response correlates with the amount of axon collateral innervation within the BLA. This suggests neurons projecting to BLA are a major driver of the tone onset responses within each population. Single neuron reconstructions showed that mPFC-BLA projection neurons are heterogeneous, including a class with no collateralization (Gao et al., 2022). Our results demonstrate that PL-BLA neurons, which encode threat predictive cues, include neurons collateralizing to NAc and cPL, suggesting synchronized transmission of this signal to multiple downstream regions may facilitate the brain dynamics necessary for learning and/or behavioral responding. Collateralization has been suggested as an understudied mechanism through which synchronous or asynchronous patterns of brain activity emerge and contribute to cognitive states (Rockland, 2018).
PL-cPL connections are important for both hemispheric synchronization and lateralization of mPFC activity. Lateralization has been associated with anxiety-like behavior (Costa et al., 2016) and stress responses in rodents (Sullivan and Gratton, 1999; Stevenson et al., 2008; Lee et al., 2015), and increased mPFC hemispheric synchrony has been associated with spatial avoidance (Dickson et al., 2022). Our data show that PL-cPL neurons persistently tracked threat, threat predictive cues, and threat avoidance behaviors. In the future, studies manipulating these neurons will be required to understand whether these activity patterns are necessary for learned or innate avoidance.
PL-VTA populations encoded learned but not innate threat avoidance behavior. Associative learning depends on the mPFC to encode predictive cues and on VTA to encode prediction errors (Lak et al., 2020). PL-VTA neurons may regulate the separable actions of these two regions. To our knowledge, we are the first to directly study the role of mPFC-VTA neurons in threat avoidance, however, VTA-mPFC neurons have been shown to respond to aversive stimuli and facilitate associative learning (Gonzalez et al., 2014; Puig et al., 2014; C.K. Kim et al., 2016; Vander Weele et al., 2018). mPFC-VTA may form part of a feedback loop influencing aversive learning. Interestingly, PL-VTA neurons showed the greatest collateral innervation of LHb. In gerbils, VTA and LHb have been shown to have opposing effects on acquisition of learned avoidance, with stimulation of VTA enhancing avoidance and stimulation of LHb impairing it (Shumake et al., 2010). PL neurons that project to both VTA and LHb may balance the activity within each region to control avoidance responses.
Footnotes
This work was supported by the NIMH Grant K01MH116264 (to L.A.D.), the NIMH Grant 1R01MH127214-01A1 (to L.A.D.), a Whitehall Foundation Research Grant (L.A.D.), a Klingenstein-Simons Foundation Grant (L.A.D.), a National Alliance for Research on Schizophrenia and Depression Young Investigator Award (L.A.D.), a Fay/Frank Seed Grant (L.A.D.), an Natural Sciences and Engineering Research Council of Canada Postgraduate Fellowship (C.B.K.), the NIMH Grant T32MH073526-14 (to M.W.G.), and the University of California, Los Angeles Medical Scientist Training Program (M.W.G.). We thank Scott Wilke, Justus Kebschull, and Xiaoyin Chen for helpful comments on the manuscript.
The authors declare no competing financial interests.
- Correspondence should be addressed to Laura A. DeNardo at ldenardo{at}ucla.edu
