A spatially-tracked single cell transcriptomics map of neuronal networks in the intrinsic cardiac nervous system

We developed a spatially-tracked single neuron transcriptomics map of an intrinsic cardiac ganglion - the right atrial ganglionic plexus (RAGP) that is a critical mediator of vagal control of the sinoatrial node (SAN) activity. We developed a 3D representation of RAGP with extensive mapping of neurons and used neuronal tracing to identify the spatial distribution of the subset of neurons that project to the SAN. RNAseq of laser capture microdissected neurons revealed a distinct composition of RAGP neurons compared to CNS neuronal subtypes. High-throughput qPCR of hundreds of laser capture microdissected single neurons led to a surprising finding that cholinergic and catecholaminergic neuronal markers Th and Chat were correlated, suggesting multipotential phenotypes that can drive neuroplasticity within RAGP. Interestingly, no single gene or module was an exclusive marker of RAGP neuronal connectivity to SAN. Neuropeptide-receptor coexpression analysis revealed a combinatorial paracrine neuromodulatory network within RAGP, informing follow-on studies on the vagal control of RAGP to regulate cardiac function in health and disease.


Introduction
In this study, we performed a spatially-tracked single cell transcriptomic analysis of an intrinsic cardiac ganglion in the pig heart to uncover the complex molecular landscape and putative paracrine neuromodulatory networks. The functional significance and complexity of the intrinsic cardiac nervous system (ICNS) has been studied for years with the majority of the focus on the physiological aspects. The ganglia at the heart are thought to constitute a "little brain" with afferent, parasympathetic and sympathetic components 1,2 .
Physiological studies using epicardial ablation have demonstrated that the intrinsic cardiac ganglia mediate central control of cardiac function through vagal and sympathetic circuits [3][4][5] . Yet, little is known about the distribution and organization of the molecular profiles of the neurons constituting the intrinsic cardiac ganglia. Recent single neuron gene expression profiling studies have uncovered the wide range of molecularly-defined subtypes 6,7 , gradient-based organization driven by inputs to neurons [8][9][10] , as well as molecular plasticity during homeostatic conditions and physiological perturbations 8,9,11 in the CNS as well as in the peripheral ganglia 7 . We set out to pursue such approaches and combine them with 3D positional information within the tissue to develop an extensive spatially-tracked molecular map of the intrinsic cardiac nervous system (ICNS).
In the present study, we focus on the right atrial ganglionic plexus (RAGP) as a key control point in the circuit that mediates vagal modulation of the SA node (SAN) activity. The routes of parasympathetic and sympathetic control of SAN were determined by several studies; surgical ablation of different areas of the heart and vagal stimulation showing a shift in pacemaker activity [12][13][14] . Stimulation of the RAGP fat pad was shown to cause a reduction in heart rate in dogs and human patients as well as cause a pacemaker shift [15][16][17] . Later studies ablating RAGP showed that this group of neurons are critical to the cardiac pacemaker response to vagal stimulation 18 . Molecular analysis of RAGP, and ICNS in general, has largely been targeted at the protein level using immunolabeling. The majority of the cardiac ganglia have been found to be cholinergic (ChAT+) [19][20][21][22] , while the proportion of observed catecholaminergic (TH+) neurons was widely variable across studies [21][22][23] . Additionally, some studies have reported co-expression of ChAT and TH in 10-20% of ICNS neurons 21,24 . The expression of other neurotransmitter/neuromodulator systems such as NPY, GAL, and SST have also been described within ICNS [25][26][27][28][29][30] . Here, we undertake a broad-based survey to characterize the gene expression of a wide range of neuronally relevant processes and their co-expression in single neurons with RAGP.
In a recent proof-of-principle study, we demonstrated a coordinated experimental approach that integrates imaging technologies with high throughput gene expression data (HT-qRTPCR) to develop a 3D anatomical and molecular map of rat ICNS 31 . Here, we build on that approach to incorporate single cell scale RNAseq and precisely integrate molecular data into a digitally reconstructed 3D RAGP with anatomical context of the pig heart. This approach contrasts with that of typical droplet-based single cell transcriptomics techniques in that the spatial and anatomical information of each sample is extensively tracked, which allows mapping of the molecular information into a 3D reconstructed anatomical organization of the tissue. Such an integrated anatomicalmolecular map permits analysis of relationships between spatial location and molecular profiles, within the tissue as well as reference to adjacent anatomical features 9,31 Our newly developed data demonstrate that the local cardiac ganglia harbor anatomical and molecular features necessary to function as complex signal processing units that critically mediate vagal control of heart function and health 32 .

Mapping spatially-tracked single cell transcriptomics onto an imaging-based 3D tissue reconstruction of pig right atrial ganglionic plexus
We developed a 3D map of the single neuron scale gene expression within the pig RAGP.
We used our recently developed method pipeline that combines single neuron anatomical position using 3D mapping with gene expression data of the mapped neurons obtained from single cell scale RNA-seq and High-Throughput qRT-PCR (HT-qPCR) 31 (Fig. 1a).
The RAGP neurons that project to SAN were labeled by a tracer that was injected into the SAN 3 weeks prior to sacrificing the animal for tissue harvesting. The heart tissue corresponding to the location of RAGP was sectioned serially from superior to inferior end for laser capture microdissection of single neurons and neuron pools (n=4 animals).
During sectioning, blockface images were obtained, which were contoured and organized into a 3D image stack (2,698 images across n=4 animals). After sectioning and staining, tissue was subjected to laser capture microdissection where both FastBlue labeled (SANprojecting) and unlabeled (considered as Non-SAN-projecting in the present analysis) single neurons were collected for microfluidic HT-qPCR (405 single neurons X 241 genes per neuron) and regional neuronal lifts were collected for single cell scale RNAseq (90 neuron pools). Image tracking was used to digitally annotate the spatial locations for single cell neuronal samples on a digitally reconstructed RAGP using TissueMapper software (Fig. 1b, Supplementary Movie 1). Expression levels of Uchl1 (PGP9.5), a pan neuronal marker, showed a wide range that is persistent throughout the RAGP with no spatial bias for enhancement or depletion (Fig. 1b,c). We assessed all single neuron samples for expression of typical neuronal markers, cholinergic and catecholaminergic markers, and key neuropeptides (Fig. 1c-e). Nearly 100% of all collected samples not only showed detectable expression, but also abundant expression of NeuN, a common neuronal marker 33 . PGP9.5 and Map2, another neuron-specific gene 34,35 , were also abundantly expressed in a large majority of the sampled single neurons (Fig. 1c). A high percentage of neurons showed abundant expression of Chat, Th, and to some extent, Dbh, suggesting a high degree of co-expression between cholinergic and catecholaminergic markers across single RAGP neurons (Fig. 1d). Neuropeptides such as Neuropeptide Y (Npy), Galanin (Gal), and Somatostatin (Sst), also showed abundant expression in a high proportion of neurons in the RAGP (Fig. 1e).
Transcriptomic landscape of pig RAGP from a single-cell scale RNAseq analysis 142 regional neuronal samples from the RAGP of a female Yucatan minipig were collected through LCM and spatially tracked as described above and subjected to singlecell scale RNAseq using the Smart-3SEQ protocol 36 (described in methods). After quality check we were left with 90 samples with detectable expression in 15,000 genes. Using publicly available data in the GTEx database 37 , we retrieved a list of genes that are statistically enriched in neuronal tissues compared to other tissue types (p < 0.01). A total of 1,639 neuronally enriched genes showed detectable expression in our RNA-seq data ( Fig. 2a). Of note, genes for ion channels associated with calcium signaling as well as glutamatergic receptors were expressed at high levels throughout the RAGP (Fig. 2b).
The transcriptomic landscape of the neuronally enriched sample set was assessed through a nonlinear embedding approach, t-stochastic neighbor embedding (tSNE) 38 . The results show that the samples were not separated into groups corresponding to distinguishable phenotypes, but were instead organized as a single cloud suggesting a gradient of gene expression based organization of underlying neuronal molecular states 8 .
Coloring the tSNE map for relative expression of Chat, an important marker for cholinergic expression, reveals a distinct gradient across the transcriptomic landscape (Fig. 2c). We then visualized the expression pattern of Chat in the context of 3D anatomical space.
Similar to the spatially distributed expression of PGP9.5 (Fig. 1b), Chat expression was widely distributed throughout the RAGP with no spatial gradients along any axis (Fig. 2d).
We compared the transcriptomic profiles of the RAGP neurons with the molecular phenotypes identified from single neuron transcriptomics analysis in the CNS 6 . We used the available mouse CNS data, due to lack of comparable pig CNS single neuron transcriptomic data set. We reconstructed the tSNE map of single neurons in the mouse CNS based on the 1,968 expression of the highly variable genes 6 (Fig. 2e). To compare how the neuronal phenotypes identified in the CNS compared to those in the RAGP, we first extracted genes contributing most to the variability in the RAGP through Principal Component Analysis (PCA). Of the 1,814 most variable genes in the RAGP, a subset of 1,412 genes were found in the mouse CNS single neuron data 6 . Visualization by tSNE map shows that the well defined neuronal clusters based on variable genes in the CNS are largely lost when the neuronal heterogeneity was analyzed using the 1,412 most variable genes in the RAGP (Fig. 2e). These results demonstrate that the highly variable genes within RAGP may not be as variable in the CNS and vice versa, suggesting a different organization of neuronal heterogeneity between RAGP and CNS structures.
Considering that a large proportion of neurons showed expression of Th and Chat (Fig.   1d), we examined the co-expression patterns of these genes in the transcriptomics data.
Neurons in the mouse 6 and human CNS 39 show mutually exclusive expression of Th and Chat. In stark contrast to the lack of Th and Chat co-expression in the CNS, we found a high degree of co-expression between Th and Chat in the RAGP neurons (Fig. 2f). This finding further augments the results seen in the tSNE map ( Fig. 2e) suggesting that the neuronal molecular states within the RAGP may be organized in a manner that is not similar to the neuronal phenotypes observed in the CNS.

Landscape of neuronal transcriptional states in the pig RAGP
Using the 3D map of a representative RAGP, examination at the single cell scale allowed us to visualize the distribution of neurons within their 3D anatomical framework, revealing that neurons, while distributed throughout the RAGP, are more densely packed closer to the endocardium (Fig. 3a). Annotation of neurons based on their projection to the SAN indicates that while both projecting and non-projecting neurons are present throughout the RAGP, SAN-projecting neurons appear to be more concentrated closer to the SAN and less concentrated towards the epicardium (Fig. 3a-c, Supplementary Movie 1). We assayed 405 spatially-tracked single neuron samples from RAGP (n=4 animals) for expression of 211 genes within each neuron using HT-qPCR, representing both SANprojecting and non SAN-projecting neurons ( Supplementary Fig. 1a). At the molecular level, a set of only 6 genes (Cck, Gal, Grp, Hcrtr1, Ntrk1, Ret) showed significant differences in the distribution of expression between SAN-projecting and SAN nonprojecting neurons (K-S statistic, FDR-adjusted p < 0.01, fold change > 2, Supplementary   Fig. 1b). Comparing neurons across all 4 RAGP, the expression distribution of select neuronal markers NeuN, PGP9.5, Chat, Th, Dbh, and Npy were relatively consistent across animals and between SAN-projecting and non SAN-projecting neurons ( Supplementary Fig. 1c).
In order to identify gene expression modules that can better characterize the neurons based on the connectivity to SAN, we used a combination of clustering and template matching. Our approach yielded six transcriptional states within the RAGP neurons (  (Fig. 4a). Visualization using a tSNE map revealed that within the context of the transcriptional landscape, these neuronal states are distributed across parts of a single cloud, suggesting a gradient-based organization of the neuronal states ( Fig. 4b). Analysis of the potential patterns and gradients of these states within the 3D anatomical framework of a representative RAGP revealed that the neuronal states were evenly distributed throughout the RAGP (Fig. 4c, Supplementary Movie 2).
We further explored the heterogeneity of gene expression distribution across these neuronal states ( Fig. 4d-f, Supplementary Fig. 1a,3a). State A is largely defined by high expression of genes in module 2 which include neuron-specific genes such as Map2, Eno2, and PGP9.5 as well as genes involved in neurotransmitter processes such as Chat, Th, and Dbh (Fig. 4d). States B and C were both characterized by high expression of genes in modules 5 and 7, which include several adrenergic receptors and potassium channels (  Supplementary Fig. 3b,c).

Correlated cholinergic and catecholaminergic gene expression in RAGP neurons
Examination of genes involved in acetylcholine and catecholamine biosynthesis and transport processes shows consistent expression across all RAGP and between SANprojecting and non SAN-projecting neurons and also reveal a surprising level of coexpression ( Fig. 5a,b, Supplementary Fig. 1a, 4a-d). A closer look at the genes involved in the catecholamine biosynthesis showed that Th, the rate limiting enzyme in the production of all catecholamines, and Dbh, responsible for the conversion of dopamine to norepinephrine are expressed over a wide range and show a high degree of coexpression ( Figure 5b, Supplementary Fig. 4b). Ddc, the enzyme converting L-DOPA to dopamine, and Pnmt, responsible for the conversion of norepinephrine to epinephrine, were both stably expressed in the majority of samples as seen by the high abundance and narrow range of expression for both genes (Fig. 5b). This further underscores the regulation of catecholamine biosynthesis process at the level of Th and Dbh as the ratelimiting enzymatic steps whose gene expression levels are highly variable across single neurons.
Visualization of key genes such as Chat and Th within the 3D anatomical framework of a representative RAGP revealed that their wide range of expression is distributed spatially throughout the RAGP ( Visualization of the expression distribution of Chat and Th through a tSNE map revealed distinct and overlapping gradients within the transcriptional landscape (Fig. 5d,e). Single cell scale RNAseq data suggested regional co-expression of Th and Chat (Fig. 2f).
Building on these results, single neuron gene expression analysis showed that Chat and Th were highly correlated across single neurons in the RAGP (Fig. 5f). Interestingly, gene co-expression of these cholinergic and catecholaminergic markers was in stark contrast to protein expression patterns that showed much reduced overlap of expression between TH and VAChT, another cholinergic marker, in individual neurons ( Supplementary Fig.   4c). Immunohistochemistry for TH and VAChT revealed that a majority of neurons showed robust protein expression of VAChT with a subset co-staining for TH (Fig. 5g,h,i).
This finding suggests that post-transcriptional regulation plays a key role in shaping the neurotransmitter patterns within RAGP. In particular, the high correlation between Th and Chat at the mRNA level suggests that the RAGP neurons are poised to use both cholinergic and catecholaminergic processes, along the lines of multipotential neuronal phenotypes observed in CNS 8 .

Neuropeptidergic interaction networks in the pig RAGP
We examined the co-expression patterns of neuropeptides and their cognate receptors to identify putative local paracrine networks within RAGP. Identifying neuronal subsets based on co-expression patterns of a neuropeptide with its receptors allows us to classify groups of neurons that exhibit autocrine or paracrine signaling. We examined local paracrine networks for three important neuropeptides Sst, Gal, and Npy in detail ( Fig. 6a-f, Supplementary Fig. 5b-d,6).
In the case of somatostatin and its receptors, Sstr1 and Sstr2, neurons cluster into six groups based on the presence or absence of each gene (Fig. 6a). This revealed one neuronal subset that expresses somatostatin, but not its receptors, and therefore can  Fig. 6b,d,f). However, Npy (and not its receptors) showed a gradient of expression that largely mirrored the Th and Chat gradients, consistent with the co-expression of these genes across single neurons ( Supplementary Fig. 6f). Taken together, the results indicate that while Sst, Gal, and Npy each have distinct co-expression patterns that outline local paracrine networks, the signalling of any one neuropeptide alone is unlikely to be the main driver of transcriptional states of RAGP neurons.
Integrating all three paracrine networks at once revealed that the overall combination is more complex than any individual neuropeptide driven network (Fig. 6g,h, Supplementary  6,7 . Transcriptomics and proteomics with 3D anatomical location tracking have been established for the brain (Allen Brain Atlas) 10, [43][44][45] .To our knowledge, we are the first group to attempt this at the mammalian heart 31 . This work enabled creation of a 3D model of pig RAGP that precisely integrates molecular data of single neurons into the 3D anatomical framework. Unlike the CNS, where molecular profiles correspond to specific anatomically located nuclei, we have found no discernable connection between the molecular profiles and anatomical locations within the RAGP 10,[43][44][45] . This suggests that the physiological functions attributed to the ICNS are likely to not be restricted to a particular anatomical location within any ganglionic plexus, and that the RAGP, and possibly other ganglionic plexuses in the ICNS, is composed of a combination of neuronal types whose integrative control of the heart enables the complex response patterns observed in physiological studies 46 .
Available physiological data suggest the notion that RAGP may consist of 80% locally connecting neurons, with the remainder distributed between those receiving afferent and specifically targeted our single cell sampling to PNs (Fig. 3c). While it is possible that some sif cells were collected along with the PNs, we believe that the data and therefore the coexpression patterns between Chat and Th can be attributed mainly to the PNs in the RAGP. Testing the protein level expression of these enzymes in the RAGP showed consistent results with previous studies in the ICNS. Confocal microscopy studies of RAGP showed that 99% neurons were VAChT positive and 10-20% neurons also positive for TH (Fig. 5 g-i). We interpret these differences between gene and protein expression levels as representative of multipotential phenotypes that are observed at the transcriptional scale that then are shaped further post-transcriptionally to yield a specific distribution in a given physiological context. We and others have shown that cells retain their plasticity beyond development and are able to adapt to perturbation based on the inputs they receive in normal and diseased states by transcriptional, as well as posttranscriptional regulation 8,11 . For example, we have previously shown that TH While Gal has been shown to attenuate cardiac vagal activity with no effect on blood pressure, Npy does just the opposite, increasing blood pressure but having no inhibitory effects on vagal activity 51

Neural tracing experiments
For initial surgery, following sedation (induction: ketamine (10mg Kg-1 IM)/midazolam (1mg Kg-1 IM), maintenance: isoflurane 1-2% inhalation) and intubation, a right unilateral thoracotomy was performed by dividing the pectoral muscle, making a small incision in the pericardium, and exposing the right atrial-superior vena cava junction. Then 5mg of Fast Blue (Polysciences), a neuronal tracer that is retained and diffuses within the lipid bilayer, in 250uL of sterile water (2% weight/volume) was injected using a 27-gauge needle into the SAN region. A chest tube was placed, and the incision was closed.
Immediately prior to removal, the chest tube was aspirated. Tissues were harvested in a terminal procedure at least 3 weeks later as described below.

Porcine tissue collection
Following sedation (induction: tiletamine-zolazepam 6mg Kg -1 IM, maintenance: isoflurane 1-2% inhalation) and intubation, we performed a midline sternotomy and exposed the heart. A heparin bolus of 5000U IV was administered, and the pig was then placed in ventricular fibrillation with application of a 9V battery to the surface of the heart.
The heart was explanted and syringe-flushed with heparinized normal saline (5U mL -1 ) via the transected aorta. The area of interest (RAGP-SAN region) was then excised and rinsed in heparinized saline. RAGPs were separated from the SANs, immersed in 1x PBS at RT for 30 seconds and transferred to 25% Optimal Cutting Temperature compound (OCT, TissueTek; VWR 25608-930), followed by 50% OCT and then embedded in 100% OCT and placed in a cryomold. The cryomold was placed in a methanol dry ice bath for flash freezing. Samples with non-zero gene counts <6000 were considered as outliers. Additionally, 10800 genes that are present in very low quantities (<30 non-zero gene counts) were filtered out. A regularized log transformation was carried out using DESeq2 58 . We normalized the filtered data using a quantile regression method SCnorm 59 . Our fInal matrix consisted of 90 samples and 15000 genes.

Extraction of significant genes from GTEx
We generated a list of neuronal genes enriched in brain and neuronal tissues using data from the Genotype Tissue Expression (GTEx) project. We downloaded the median gene expression values for all tissue types from GTEx Analysis V8 (dbGaP Accession phs000424.v8.p2). Pavlidis Template Matching was used to find genes that were specifically enriched in brain and neuronal tissues (template-maximum for neuronal tissues, minimum for non-neuronal tissues) with p value < 0.01.

High-throughput real-time PCR
Single RAGP neurons in lysis buffer (Cells Direct Lysis Buffer, Invitrogen) were directly processed for reverse transcriptase reaction using SuperScript VILO Master Mix (Thermo Fisher Scientific, Waltham, MA), followed by real-time PCR for targeted amplification and detection using the Evagreen intercalated dye-based approach to detect the PCRamplified product. Intron-spanning PCR primers were designed for every assay using

HT-qPCR data analysis
Individual qRT-PCR results were examined to determine the quality of the qRT-PCR based on melt-curve analysis. Following this initial quality control, samples with >30% for at least two other states (Figure 3d).

Plots for showing abundance and range of gene expression in Ct space
In order to accurately display the range of the normalized data while also keeping in mind the overall abundance of each gene throughout the 4 RAGP, the -ΔΔCt data was

Data Availability
The authors declare that all the data supporting the findings of this study are available within the article and its supplementary information files or from the corresponding author (RRID:SCR_017001) project and is available at https://portal.brain-map.org/atlases-anddata/rnaseq/human-multiple-cortical-areas-smart-seq.       Samples had two replicates (a,c) and genes had three replicates (b,d). To robustly assess the three gene replicates, we compared each of the three sets to each other, shown in three different colors, where blue represents the comparison between sets 1 and 2, green represents the comparison between sets 1 and 3, and red represents the comparison between sets 2 and 3.

Supplementary Movie 1:
Acquisition of single neurons from the RAGP and 3D visualization of SAN-projecting and non-SAN-projecting neurons within a representative RAGP.

Supplementary Movie 2:
Visualization of identified transcriptional states within the 3D anatomical framework of a representative RAGP.

Supplementary Movie 3:
3D spatial distribution of Chat and Th expression across a representative RAGP.

Supplementary Movie 4:
Neuropeptidergic interaction networks within the 3D anatomical framework of a representative RAGP.

Supplementary Cytoscape Network File:
Cytoscape file containing all the network information related to Figure 6 and