An integrated single-cell RNA-seq atlas of the mouse hypothalamic paraventricular nucleus links transcriptional and functional types

The hypothalamic paraventricular nucleus (PVN) is a highly complex brain region that is crucial for homeostatic regulation through neuroendocrine signalling, outflow of the autonomic nervous system, and projections to other brain areas. The past years, single-cell datasets of the hypothalamus have contributed immensely to the current understanding of the diverse hypothalamic cellular composition. While the PVN has been adequately classified functionally, its molecular classification is currently still insufficient. To address this, we created a detailed atlas of PVN transcriptional cell types by integrating various PVN single-cell datasets into a recently published hypothalamus single-cell transcriptome atlas. Furthermore, we functionally profiled transcriptional cell types, based on relevant literature, existing retrograde tracing data and existing single-cell data of a PVN-projection target region. In our PVN atlas dataset, we identify the well-known different neuropeptide types, each composed of multiple novel subtypes. We identify Avp-Tac1, Avp-Th, Oxt-Foxp1, Crh-Nr3c1 and Trh-Nfib as the most important neuroendocrine subtypes based on markers described in literature. To characterize the pre-autonomic functional population, we integrated a single-cell retrograde tracing study of spinally-projecting pre-autonomic neurons into our PVN atlas. We identify these (pre-sympathetic) neurons to co-cluster with the Adarb2 + clusters in our dataset. Finally, we identify expression of receptors for Crh, Oxt, Penk, Sst, and Trh in the dorsal motor nucleus of the vagus, a key region that pre-parasympathetic PVN neurons project to. Concluding, our study present a detailed overview of the transcriptional cell types of the murine PVN, and provides a first attempt to resolve functionality for the identified populations.

parasympathetic PVN neurons project to.Finally, we identify Trh-Ucn3 and Brs3-Adarb2 as some centrally projecting populations.In conclusion, our study presents a detailed overview of the transcriptomic cell types of the murine PVN and provides a first attempt to resolve functionality for the identified populations.

| INTRODUCTION
The hypothalamus is a highly complex brain region that is crucial for homeostatic regulation.It consists of several nuclei that are functionally and molecularly diverse.Single-cell RNA sequencing has contributed immensely to the current understanding of the cellular composition of the hypothalamus.A recently published single-cell atlas based on 17 single-cell datasets currently provides the most high-resolution insight into the cellular makeup of the murine hypothalamus. 1Yet, not all hypothalamic nuclei are represented at equal resolution in this atlas.
One of the underrepresented nuclei is the paraventricular nucleus (PVN), a small nucleus of the anterior hypothalamus, adjacent to the third ventricle.[4][5] The neuroendocrine subset of neurons is usually described as being subdivided into the large magnocellular and smaller parvocellular neurons projecting to either the posterior pituitary or the median eminence.The magnocellular neurons release the neurohormones AVP and OXT into the systemic bloodstream via the posterior pituitary, whereas parvocellular neurons release CRH, TRH, and SST in the median eminence to be transported to the anterior pituitary via the hypophyseal portal system.In the anterior pituitary, CRH/AVP, TRH, and SST act on the different populations of trophic cells to stimulate or inhibit the hypothalamus-pituitary-adrenal (HPA), hypothalamuspituitary-thyroid, and hypothalamus-pituitary-somatic (HPS) axes, respectively.
The preautonomic neurons project to the sympathetic and parasympathetic preganglionic neurons in the intermediolateral nucleus (IML) of the spinal cord and dorsal motor nucleus of the vagus (DMV) in the brainstem, respectively.In addition, these preautonomic neurons project to other brainstem nuclei such as the nucleus of the solitary tract, the central gray, and the raphe nucleus. 2 Though the existence of these projections has been described, the molecular identity of these neurons has not been completely elucidated yet.[8] These studies have revealed the expression of several neuropeptides, like AVP, OXT, enkephalin (PENK), and dynorphin (PDYN) in these neurons.However, these studies were only able to characterize roughly half of retrogradely traced neurons, thus necessitating further characterization of the neuropeptides and other markers involved.
Neuropeptidergic projections to other central locations upstream of the brainstem exist as well and are often associated with behavioral modulation.For instance, magnocellular Oxt + neurons were reported to project to the amygdala, 9 and parvocellular Oxt + neurons were reported to project to the nucleus of the solitary tract and nucleus accumbens. 5,9Both magnocellular and parvocellular Oxt + neurons were also reported to project to the ventral tegmental area. 10Further literature on oxytocinergic central connections has been extensively reviewed by Grinevich and Neumann. 11Similarly, the central projections and function of Avp + neurons have been reported extensively. 12r other neuropeptides, like SST and TRH, the literature on central projections is unfortunately less extensive.
While the PVN has been adequately classified functionally (Figure 1A), the molecular classification is currently insufficient.Even less well-described is the relation between the molecular and functional classifications.For instance, it is not known whether the centrally projecting neurons are a separate population of PVN neurons or a subset of the neuroendocrine and/or preautonomic neurons.This lack of insight in the complexity of this region is impeding research progress, since targeting a single functional pathway is difficult without knowledge of the transcriptomic identity of the different functional types.To address this issue, we aimed to create a comprehensive atlas of PVN transcriptomic types from existing scRNA-seq datasets.Furthermore, we aimed to associate transcriptomic and functional types based on literature of neuroendocrine neuron-enriched genes, existing retrograde tracing data, and a singlecell dataset of a PVN-projection target region.

| PVN neuropeptidergic neurons consist of multiple subpopulations
4][15][16] These datasets were integrated using scANVI, 17 an integration method derived from scVI, 18 the method used for the original HypoMap integration.From the integrated dataset, PVN neurons were selected, and the resulting subset was clustered.Finally, the identified clusters were annotated based on enriched genes (Figure 1B).
After integration and annotation of the dataset, we first validated the integration by visually inspecting the intermixing of datasets in the dimensionality reduction plot (Figure S1A).Then, we set out to characterize the neuronal populations and subpopulations.All expected major neuropeptidergic neuron populations were found (Figure 1C,D), as well as some nonneuropeptidergic populations.
Three populations in the dataset were excluded based on in situ expression data: Tcf7l2-Shox2, Slc32a1-Dlx1, and Trh-Zic1 neurons (Figure S2A).Anatomically, markers for these populations are found primarily in the thalamus and/or subparaventricular zone, especially compared with the rare immunoreactive neurons within the PVN (Figure S2B-D).Thus, the relevance of these neurons for PVN function is unclear, given the low frequency of occurrence in the PVN and their high abundance in bordering regions.All remaining clusters were positive for the vesicular glutamate transporter gene Slc17a6 and overall lacked the vesicular GABA transporter gene Slc32a1 (Figure S1B).A rare Slc32a1 + population has been described in the PVN, 19 raising the question if this population could be reproduced in the atlas.Although sparse Slc32a1 expression can be observed in the atlas, it is unclear if this is of biological relevance, as these neurons did not form distinct Slc32a1 + (sub)clusters within the atlas taxonomy.
Interestingly, despite the lack of Slc32a1, some of the neurons were positive for GABA biosynthesis enzymes, in particular for Gad2 (Figure S1B).
The identified clusters were annotated as neuropeptidergic subtypes, except the Asb4-Adarb2, Brs3-Adarb2, and Slc17a6-Adarb2 populations, which were not annotated as neuropeptidergic because of relatively low neuropeptide expression levels.Neuropeptidergic populations were further divided into subpopulations based on secondary markers.
The Avp + and Oxt + neurons were found to often present some level of coexpression of both genes (Figure 1D), despite this being rarely observed at protein level. 20These clusters were characterized as either Avp + or Oxt + , based on which of the two genes was most abundantly expressed.This yielded Avp-Adarb2, Avp-Tac1, Avp-Th, Oxt-Foxp1, and Oxt-Adarb2 as the subpopulations of Avp + and Oxt + neurons (Figure 1D).We further identified one population of Penk + neurons, Penk-Adarb2; two subpopulations of Crh + neurons, Crh-Nr3c1 and Crh-Adarb2; two subpopulations of Sst + neurons, Sst-Calb2 and Sst-Sfrp2; and two subpopulations of Ghrh + neurons, Ghrh-Adarb2 and Trh-Ghrh (Figure 1D).
Finally, the Trh + neurons were split in multiple categories: Trh-Ghrh, Trh-Nfib, Trh-Omp, Trh-Tac1, and Trh-Ucn3 (Figure 1D).Notably, while the Brs3-Adarb2 cluster was not annotated with a neuropeptide, it actually did express Trh (Figure S1C).This expression was present at relatively low levels compared with the other Trh + clusters and was found in a sequencing method-dependent manner (Figure S1C).As such, the Brs3-Adarb2 population was included in further analyses pertaining to Trh + populations.
We then interrogated the expression of neuropeptidergic receptors in the various populations, to identify possible autocrine and paracrine signaling in the PVN (Figure 1E).Most populations expressed neuropeptidergic receptors that may facilitate local signaling, and some populations expressed receptors for their own neuropeptides.

| Neuroendocrine markers are enriched in distinct neuropeptide subpopulations
We next set out to identify which clusters may be linked to neuroendocrine function.In each neuropeptidergic type (Figure 2A), genes previously identified to be expressed in neuroendocrine cells were analyzed for presence in subpopulations.Furthermore, we compared similarities between subpopulations, revealing that not all neurons within a neuropeptide class were transcriptomically most similar to each other (Figure 2B).For instance, the Oxt-Foxp1 neurons were more similar to the Avp-Tac1 and Avp-Th neurons than the Oxt-Adarb2.In turn, the Oxt-Adarb2 neurons were most similar to other Adarb2 + neurons (Figure 2B).
To functionally classify the Oxt + neurons, we used the neuroendocrine Oxt + markers as described by Lewis et al. 5 In their study, a large overlap between magnocellular electrophysiology and neuroendocrine projection Oxt + neurons was found, and this magnocellular neuroendocrine population was reported to be enriched in Calb1, Kcnmb4, and Foxp1.Notably, a closer inspection of their single-cell data reveals that their magnocellular neuroendocrine cluster contains both Avp + and Oxt + neurons (Figure S3).Thus, these markers for the magnocellular neuroendocrine population are relevant to not only Oxt + neurons but also Avp + neurons.
In our data, these magnocellular neuroendocrine markers were found to be enriched in the Oxt-Foxp1 cells relative to the Oxt-Adarb2 cells, suggesting the former to be the neuroendocrine population (Figure 2C).Both Avp-Tac1 and Avp-Th clusters were found to be strongly correlated with the magnocellular neuroendocrine Oxt + neurons (Figure 2B), whereas the Avp-Adarb2 cluster was not.Furthermore, the Avp-Tac1 and Avp-Th populations were enriched in the magnocellular neuroendocrine markers relative to Avp-Adarb2 (Figure 2C).Between Avp-Tac1 and Avp-Th, the Avp-Th population was slightly more enriched in these markers, but the relative differences were smaller than the difference in the Avp-Adarb2 population or between Oxt + subpopulations.
Both the high degree of similarity between the Avp-Tac1 and Avp-Th populations, and the expression of neuroendocrine markers in both populations imply Avp-Tac1, and Avp-Th may be magnocellular neuroendocrine Avp + neurons.It is currently not known if the differences between Avp-Tac1 and Avp-Th neurons are functionally relevant, as the distinction has not been described before.It should be noted that, although positive for Th, the Avp-Th neurons do not express Ddc (Figure S4A).Thus, these neurons cannot perform the final dopamine biosynthesis conversion, and therefore are not dopaminergic.
For the Crh + neurons, multiple markers for neuroendocrine function have been described.First and foremost, the glucocorticoid receptor (Nr3c1; GR) is an essential mediator of regulatory feedback in Crh + neurons of the HPA axis. 21Another key component of these neurons is Avp, which strongly potentiates CRH-driven adrenocorticotropic hormone (ACTH) release from the anterior pituitary. 22Further, Scgn has been proposed as a marker of neuroendocrine Crh + cells, as Scgn silencing abolishes stress-driven ACTH release. 23Finally, Agtr1a has been found to be expressed in neuroendocrine Crh + and Trh + neurons. 24All these markers were enriched in the Crh-Nr3c1 population, suggesting these neurons to be the neuroendocrine Crh + population (Figure 2D).To validate this finding, we performed an immunofluorescence staining of SCGN in the rat PVN, combined with retrograde tracing of neuroendocrine and preautonomic neurons.
Expression of SCGN was found to be highly colocalized with Fluoro-Gold (FG)-labeled neuroendocrine neurons (Figures 2E and S4B).In contrast, no colocalization was found in SCGN with retrograde labeling of IML-projecting neurons (Figures 2F and S4B).
Like in Crh + neurons, neuroendocrine Trh + neurons are reported to express Agtr1a. 24Additionally, literature indicates that the thyroid hormone receptor beta is critical to the T 3 -mediated negative feedback on Trh expression. 25Besides these two, other markers for neuroendocrine Trh + cells have not been described in the literature.While several Trh + populations expressed only one of these markers, only the Trh-Nfib neurons were enriched for both (Figure 2G), suggesting those to be the neuroendocrine population.
Finally, no markers specific for neuroendocrine Sst + subpopulations have been described in the literature.The identified clusters Sst-Calb2 and Sst-Sfrp2 do seem distinct in their expression of various genes.In Sst-Sfrp2 neurons, Sfrp2, Ar, Npy1r, and Mc4r were enriched, whereas in the Sst-Calb2 neurons, Calb1, Calb2, Kcnip4, and Gpr101 were enriched (Figure S4C).

| Preautonomic neurons are embedded in the Adarb2 + clusters
To identify preautonomic neurons in our data, we integrated a recent single-cell dataset that utilized retrograde tracing from the IML. 26om this dataset, all neurons identified by the authors to be from the PVN were integrated into our PVN atlas.After integration, the traced neurons were embedded within the populations using Adarb2 as secondary marker, here referred to as the Adarb2 + populations (Figure 3A).Similar to the Adarb2 + neurons, the traced neurons were high in their expression of Adarb2, Reln, Snca, and Ntng1 (Figures 3B   and S5A).Further similarity was found in the expression of Cacna1g (Figure S5B), a T-type calcium channel gene reported to mediate lowthreshold spiking in preautonomic neurons. 27,28 further validate the association between Adarb2 + neurons and preautonomic function, we did an immunofluorescence staining in the rat PVN for RELN, one of the Adarb2 + -enriched markers.We found strong colocalization of RELN with retrograde labeling of IMLprojecting neurons (Figures 3C and S5B).In contrast, no colocalization was found in RELN with FG-labeled neuroendocrine neurons (Figures 3D and S5B).
We then analyzed Adarb2 + subpopulations within the PVN atlas.
Among the Adarb2 + subpopulations, the traced neurons were embedded in particular within the Avp-Adarb2, Oxt-Adarb2, Penk-Adarb2, and Crh-Adarb2 clusters (Figure 3E).The traced neurons expressed all these neuropeptides to some extent, though the expression of Oxt was markedly the most pronounced (Figure 3E).
Besides the IML, the preautonomic neurons of the PVN also project to the DMV.Although there are currently no retrograde tracing sequencing datasets available, a recent single-cell sequencing dataset of the general DMV 29 may still give insight into possibly relevant circuitry, through its neuropeptide receptor repertoire.In their study, Tao et al. 29 microdissected the DMV based on Chat-driven fluorescence and identified seven distinct DMV populations.We analyzed their dataset and found neuropeptidergic receptors that relate to PVN-expressed neuropeptides.
The identified receptors were Crhr1, Trhr, Sstr2, Sstr5, Oprm1, Oprd1, and Oxtr (Figure S5C).In summary, our results indicate that sympathetic preautonomic neurons include Avp-Adarb2, Crh-Adarb2, Oxt-Adarb2, and Penk-Adarb2 neurons.Brainstem-projecting parasympathetic preautonomic neurons may be more diverse, with receptors for Crh, Trh, Sst, Penk, and Oxt present in the DMV.The Trh-Ucn3 population is synonymous with the perifornical population as described by Péterfi et al., 30 Autry et al., 31 and Horii-Hayashi et al. 32 Although generally described as perifornical, neurons of this Ucn3 + population exist within the PVN boundaries 33 (Figure S6), hence its inclusion in our atlas.This population has been described as being involved in various behaviors, ranging from feeding, 30 to infant-directed aggression, 31 to risk assessment behavior. 32 addition to the Trh-Ucn3 population, we investigated the Mc4r + population of the PVN, which has previously been described to  26 and the rest of the dataset, and columns are split by expression of Avp, Crh, Oxt, and Penk.be involved in feeding behavior. 34,35Expression of Mc4r is found in three populations, Sst-Sfrp2 and Brs3-Adarb2, and to lesser extent Trh-Nfib (Figure 4A,B).In literature, the expression of Mc4r has been described in neuroendocrine Trh + and Sst + neurons. 36,37Thus, the remaining undescribed Brs3-Adarb2 population is left as an interesting candidate for the feeding-related Mc4r + population.Besides Mc4r, the Brs3-Adarb2 population expresses Brs3, Npy1r, Npy5r, Vgf, and Asb4 (Figure 4B).[40][41][42] Interestingly, Mc4r, Brs3, Npy1r, and Vgf have also been associated with regulating energy expenditure through the sympathetic nervous system, 35,42,43 suggesting the Brs3-Adarb2 population to also be involved in sympathetic outflow.This indicates that at least part of the preautonomic neurons may also project centrally, which further highlights the difficulty in comprehensively classifying the central projections of the PVN.Nevertheless, we were able to identify a putative population involved in feeding behavior and energy expenditure, while harmonizing existing literature on the various effector genes expressed within this population.

| DISCUSSION
In this study, we integrated single-cell PVN datasets into a general hypothalamic atlas and subsequently selected PVN-specific neurons from the integrated dataset.We characterized the resulting dataset, revealing multiple transcriptomic subpopulations for neuropeptidergic neurons.For the Avp + , Oxt + , Crh + , and Trh + neurons, we identified distinct subpopulations that were enriched for neuroendocrine marker genes.We describe how spinally projecting preautonomic neurons cluster entirely with various Adarb2 + neurons, and we identify putative neuropeptidergic types of brainstem-projecting preautonomic neurons.
Finally, based on markers described in literature, we identify two putative populations involved in the central projections from the PVN.
Since we used a microdissected high-depth sequencing dataset as benchmark, 16 our final dataset is unlikely to contain any non-PVN neuron populations.Rare cell types for the PVN that are much more prominent in adjacent brain regions were excluded based on in situ hybridization data.Nonetheless, this approach also has its limitations, as the sample size of the benchmark dataset was relatively small.Rare transcriptomic populations, or populations restricted to undersampled PVN subregions, are less likely to be uncovered through this method.
Another limitation of our study is linked to the use of the retrogradely traced single-cell data from Beine et al. 26 Their study injected the retrograde tracer into the L1 spinal cord, while the IML spans from T2 through L2 in mice. 44Thus, the traced population is likely only a limited representation of the spinally projecting preautonomic neurons, as the axon terminal location could be a source of transcriptomic variation in preautonomic neurons.
Despite the limitations, we were able to identify multiple subpopulations for all neuropeptidergic neuron populations previously shown to reside in the PVN.For Avp + neurons, we identified three distinct subpopulations that have not been described before: Avp-Adarb2, Avp-Tac1, and Avp-Th.While the presence of Th in Avp + neurons is not novel, 45 the presence of Th as marker for a distinct subpopulation of Avp + neurons has not been described before.The Avp-Tac1 and Avp-Th subpopulations were highly correlated with the magnocellular neuroendocrine Oxt + neurons, and both were high in expression of the neuroendocrine markers as described by Lewis et al. 5 These results indicate that these Avp + populations may be magnocellular neuroendocrine neurons.On the other hand, the Avp-Adarb2 population was similar to the retrograde traced neuron, indicating a function as preautonomic neurons.
The identified Oxt + neuron subpopulations were transcriptomically the same as identified by Lewis et al. 5 -with the side note that the magnocellular population in their dataset was composed of both Avp + and Oxt + neuroendocrine neurons.We identify the Oxt-Foxp1 neurons to be the magnocellular neuroendocrine Oxt + population and identify the parvocellular Oxt-Adarb2 neurons to be preautonomic, based on the transcriptomic overlap with retrogradely traced neurons.
However, it is important to note that Oxt-Adarb2 neurons are not only functionally preautonomic but also centrally projecting.For instance, these neurons project to nucleus accumbens and nucleus of the solitary tract. 5,9r the Crh + neurons, two subpopulations were identified: Crh-Nr3c1 and Crh-Adarb2.These Crh + subpopulations are synonymous to the populations, as discussed by Romanov and Harkany 46 : Crh-Scgn and Crh-Fam150b, respectively.While these markers would have worked for our Crh + populations as well, clusters were annotated with different markers more fitting for the context of our dataset.We found that the Crh-Nr3c1 neurons were enriched for multiple neuroendocrine markers genes, while the Crh-Adarb2 neurons were transcriptomically similar to a subset of spinally projecting preautonomic neurons.In contrast, Romanov and Harkany 46 posit that both subpopulations are neuroendocrine and that other markers-like Npy1rdefine subsets of these populations that are preautonomic instead of neuroendocrine.We do not find any evidence for this theory.First, neuroendocrine markers like Agtr1a, Nr3c1, and Scgn are lacking entirely in the Crh-Adarb2 neurons, while Avp was also relatively low in expression in these cells.Second, no retrogradely traced preautonomic neurons were found coclustering with the Crh-Nr3c1 neurons.We identify this population to be synonymous with the perifornical Ucn3 + population.Although generally described as perifornical, neurons from this population reside in the PVN as well.[32] The identified Sst + subpopulations could not be linked to a functional identity due to lacking literature on functional Sst + subpopulations.Despite this, two of the identified differential genes do stand out.First, the Ar has been described in literature to define a subgroup of Sst + neurons that is present in both sexes, but in larger numbers in males. 47This population is postulated to be relevant for sexually dimorphic action of SST on the HPS axis.The other notable gene is Gpr101, which has been reported as the causative gene to X-linked acrogigantism-a disorder involving the growth hormone axis. 48Since both populations can be linked to the HPS axis, through either Ar or Gpr101, we hypothesize that both subpopulations are neuroendocrine in nature, with the two populations acting as separate substrates for HPS axis inhibition.Furthermore, we analyzed the neuropeptide receptor repertoire of the DMV.Receptors for CRH, OXT, PENK, and SST were present, as well as the receptor for TRH to a lesser extent.This receptor expression profile is consistent with literature on PVN-DMV neurons, which reports the expression of AVP, OXT, PENK, SST, and CRH in preparasympathetic neurons. 8,49The presence of TRH receptor might indicate relevance for this neuropeptide in PVN-DMV circuitry as well, though TRH projections could also originate from brain regions other than the PVN.As the DMV-projecting AVP and OXT neurons were reportedly parvocellular, we can infer that the Avp-Adarb2 and Oxt-Adarb2 populations are the most likely candidates for the vasopressinergic and oxytocinergic PVN-DMV connections.This would implicate these populations in both presympathetic and preparasympathetic functions, with currently no known differentiating factor to define these functional subpopulations.In summary, our findings reinforce and expand on current literature on PVN-DMV circuitry and introduce novel avenues to elucidate neuropeptidergic subpopulations involved in this circuit.
Finally, we characterized some central projections of the PVN.
We identified a likely candidate for the feeding-related Mc4r + population, Brs3-Adarb2.9][40][41][42] Like the Oxt-Adarb2 population, the Brs3-Adarb2 population seems to be involved in both preautonomic and central functions.It is unclear whether the dual-function populations are homogenous, in which these functions are inseparably linked, or if discrete subpopulations may be responsible for these different functions.Further research and data are needed to increase the resolution of these populations, and thereby resolve this question.
In conclusion, our study presents a detailed overview of the neuronal transcriptomic subpopulations in the PVN and attempts to resolve functional identities for the identified populations (Figure 5).
The study pinpoints interesting markers of neuropeptidergic subpopulations that may be used for research into functionalities that are specifically associated with these subpopulations.

| scRNA-seq data collection and preprocessing
PVN scRNA-seq datasets were acquired from the relevant Gene Expression Omnibus (GEO) repositories (Table 1).The hypothalamic atlas data (HypoMap) were acquired from CELLxGENE (Table 1).Each dataset was preprocessed separately to filter low-quality cells based on transcript count, gene count, and percentage of mitochondrial gene expression (Table S1).Datasets using Ensembl IDs were converted to MGI symbols.Furthermore, dataset metadata was standardized to include assay, disease, organism, sex, SRA_ID, Sample_ID, and Dataset, using the HypoMap atlas metadata as template.Within the atlas, cells were annotated at various levels of granularity, ranging from C7_named to C465_named.The number in these annotation levels denotes the number of clusters at that annotation level.Nonatlas datasets were annotated at the C7_named level by comparison of cluster marker genes with the respective expression patterns in the atlas, assigning clusters in the nonatlas dataset based on these comparisons.Finally, the datasets were merged on the raw counts and  metadata.For further analysis, we selected neuronal cells only, based on the C7_named annotation key.Additionally, the first principal component of the immediate early genes (IEGs) Fos, Fosb, Jun, Junb, Egr1, Egr2, Npas4, Nr4a1, Nr4a2, and Nr4a3 was calculated on normalized gene expression, to regress out the effect of IEGs during integration.

| scRNA-seq data integration
The highly variable genes of the merged dataset were selected using the function highly_variable_genes from scanpy 50  To isolate PVN-specific cells from the integrated dataset, the cells from Xu et al. 16 were used as benchmark.In that study, PVN boundaries were identified based on fluorescent agouti-related peptide (AgRP) axon tracts, and PVN neurons were microdissected and sequenced using a plate-based sequencing method (SMART-SCRB).
The resulting dataset was thus deemed the most accurate representation of the PVN because of the combination of high sequencing depth and precise cell isolation method used.The other PVN datasets were either selective for a single neuropeptide neuron type, or not fully specific for the PVN, due to contamination with neurons from bordering regions.
As such, a subset was taken of predicted C66_named clusters that contained five or more cells originating from Xu et al. 16 This subset was clustered using the Leiden algorithm 51 with n_neighbors = 30 and resolution = 1.A subset was taken of the resulting clusters, again retaining only clusters with at least five cells from Xu et al. 16 These subsetting steps retained 98.58% and 97.56% of the cells from Xu et al., 16 respectively.The resulting subset was then reintegrated using the same steps and parameters as before, except for the highly variable gene selection with n_top_genes = 1000.
The scANVI latent embedding was used to cluster the data with the Leiden algorithm 51 with n_neighbors = 30 and resolution = 1.
Each of the resultant 20 clusters was analyzed for the presence of Xu et al. 16 neurons.Only clusters with the presence of minimally five of these neurons were annotated.An exception was made for Trh-Ucn3, as Allen Brain Atlas in situ hybridization data showed Ucn3 expression in the PVN (Figure S6).For annotation, these clusters were tested for potential subclusters using the FindSubClusters function, and subclusters containing no cells from Xu et al. 16 were excluded.The resultant clusters and subclusters were tested for enriched genes using the FindMarkers functions and were annotated accordingly.Finally, a subset was taken with all annotated clusters, containing only the annotated populations described in the Results section.
To analyze the retrogradely traced neurons from Beine et al., 26 the PVN neurons within that dataset were selected based on Sim1 expression.The selected neurons were integrated into the PVN atlas, using the same steps and parameters for scANVI integration as described before.Altered parameters were n_top_genes = 1000 and early_stopping = False.The latent representation of this final integration was used for visualization of the data, using the preserve_neighbors function from the minimally distorting embedding (MDE) framework. 52Parameters used for MDE visualization were attracti-ve_penalty = Huber, repulsive_penalty = Log, n_neighbors = 30, and device = "cuda." Finally, for the analysis of the neuropeptide receptors in the DMV, preprocessing and clustering steps as described by Tao et al. 29 were replicated.

| Production of virus
The AAV2-retro-CAG-nls-GFP virus was prepared following the protocol as outlined by Verhaagen et al. 53 Briefly, HEK293T cells were cultured in 15 cm dishes with DMEM containing 10% fetal calf serum and penicillin/streptomycin.The next day, the cells were cotransfected with AAV2-retro serotype helper plasmid, 54 helper plasmid pAdδ F6, and pAAV-CAG-nls-GFP 55   Because of the leaky blood-brain barrier at the median eminence and posterior pituitary, FG is taken up by the axonal terminals of PVN neuroendocrine neurons.

| Tissue preparation
The animals were anesthetized with an overdose of Pentobarbital and then perfused transcardially with 150 mL of 0.  Images were acquired using a Leica SP8 confocal microscope.

| Statistics
Single-cell transcriptomics pipelines and statistical computations were performed using both Python 3.7.10 and R 4.1.2. 56Packages used for the analysis were tidyverse 2.0.0,Seurat 4.

F I G U R E 1
Neuropeptidergic paraventricular nucleus (PVN) neurons consist of multiple subpopulations.(A) Schematic showing the various projections from the PVN.Arrows point to the dorsal motor nucleus of the vagus (DMV), intermediolateral nucleus (IML), anterior pituitary (AP), and posterior pituitary (PP).Other central projections are denoted with the upward pink arrows.(B) Schematic of the methods used for this section.Five PVN datasets were integrated with the hypothalamus atlas (HypoMap) using scANVI.The integrated low-dimensional embedding was used for clustering and visualization.Clusters were annotated based on differential gene expression.(C) Minimally distorting embedding (MDE) representation of the datasets, with the cluster annotation based on neuropeptide markers and secondary gene markers.(D) Violin plots of the expression of PVN-specific transcription factors (TFs), neuropeptides, and secondary cluster markers.(E) Neuropeptidergic receptors in the various populations of the PVN.Dot size represents the percentage of cells expressing a gene, and dot color represents the average expression level.F I G U R E 2 Legend on next page.BERKHOUT ET AL.

Finally, we set
out to investigate the other central projections, upstream of the brainstem.While we assume that the Trh-Omp, Trh-Ghrh, and Trh-Tac1 populations are centrally projecting, we unfortunately could not find any literature regarding the functional identities of these populations.F I G U R E 2 Expression of neuroendocrine markers in subpopulations of Avp + , Oxt + , Crh + , and Trh + neurons.(A) Minimally distorting embedding (MDE) representation of the data, colored by neuropeptide class.Coloring identical to violin plot colors in panels (C,D,E).(B) Heatmap of correlations between clusters, calculated on pseudobulk expression values of highly variable genes.(C) Violin plots of the expression of neuroendocrine Avp + and Oxt + markers.(D) Violin plots of the expression of neuroendocrine Crh + markers.(E) Immunofluorescent staining of SCGN in the rat paraventricular nucleus (PVN), showing colocalization with FluoroGold (FG).(F) Immunofluorescent staining of SCGN in the rat PVN, showing no colocalization with retrograde tracing from the intermediolateral nucleus.(G) Violin plots of the expression of neuroendocrine Trh + markers.

F I G U R E 3
Preautonomic neurons cocluster with the Adarb2 + neurons.(A) Violin plots showing expression of Adarb2, Reln, Snca, and Ntng1 in the paraventricular nucleus (PVN) atlas and retrogradely traced PVN neuron from Beine et al. 26 (B) The minimally distorting embedding (MDE) representation, with cells from Beine et al. 26 highlighted in red.(C) Immunofluorescent staining of RELN in the rat PVN, showing colocalization with retrograde tracing from the intermediolateral nucleus.(D) Immunofluorescent staining of RELN in the rat PVN, showing no colocalization with FluoroGold (FG).(E) Zoomed-in MDE representation, showing the Adarb2 + clusters.Rows show cells from Beine et al.

F
I G U R E 4 The Brs3-Adarb2 population expresses various feeding-related and energy expenditure-related genes, including Mc4r.(A) Violin plots showing the expression of Mc4r in the paraventricular nucleus (PVN) atlas.(B) Dot plot showing the expression of Mc4r, Brs3, Npy1r, Npy5r, Vgf, and Asb4 in the PVN atlas.Dot size represents the percentage of cells expressing a gene, and dot color represents the average expression level.

Finally, our
immunofluorescence validations do not support a preautonomic function of the Scgn-expressing Crh-Nr3c1 subpopulation nor a neuroendocrine function of Reln-expressing Crh-Adarb2 neurons.Together, this suggests that these two Crh + subpopulations are functionally distinct.Multiple subpopulations of Trh + neurons were also identified in our analysis, but the functional interpretation [for these Trh subpopulations] was more difficult than [the functional interpretation] for other [non-Trh] populations.We observed an enrichment of two neuroendocrine markers in Trh-Nfib neurons.Aside from neuroendocrine function, only the Trh-Ucn3 population could be functionally resolved.

F I G U R E 5
Summarized findings from our PVN atlas.(A) Overview of projections from the paraventricular nucleus (PVN) to the intermediolateral nucleus (IML), dorsal motor nucleus of the vagus (DMV), and other central projections with relevant cell populations indicated.Neuropeptides and neuropeptide receptors are shown, where exact clusters could not be identified.(B) Overview of projections from the PVN to the pituitary, with relevant cell populations shown.T A B L E 1 Datasets used in this article.
Adult male Wistar rats (Charles River, Germany) weighing 200-250 g were group-housed in a 12/12-h light/dark cycle with ad libitum access to food and water.All animal experiments were carried out in accordance with the guidelines of the Animal Ethics Committee of the Royal Dutch Academy of Arts and Sciences (KNAW, Amsterdam,The Netherlands) and approved by the Netherlands Institute for Neuroscience (NIN, Amsterdam, The Netherlands).
using PEI (polyethylenimine 25 kDA, linear).After 72 h of transfection, the cells were harvested, lysed, and centrifuged.Viral particles were purified from the supernatant using iodixanol density gradient centrifugation with an Amicon Ultra-15 centrifugal filter.The final viral titers were determined by quantitative polymerase chain reaction for the WPRE gene in the viral genome (vg) using the following primers: 5 0 -GTGTTGCCACCTG-GATT-3 0 and 5 0 -CGAAGGGACGTAGCAGAA-3 0 .The titer of the viral batch used was 4E+12 vg/mL.

4. 5 |
Surgical proceduresRats were anesthetized by an intramuscular injection of Medetomidine (0.26 mg/kg) and ketamine (64 mg/kg) mix.For the retrograde tracing of the preautonomic motor neurons in the PVN, the spinal cord was carefully exposed via a small laminectomy.Two to three microliters of AAV2-retro-CAG-nls-GFP virus was delivered to the IML column of the spinal cord at T5-T8 using a Nanoject III, and coordinates 0.05 mm lateral and 0.1 mm ventral from the surface of the spinal cord.The glass capillary needle was left in place for 2 min after the injection to minimize leakage of the virus.Animals were sacrificed for further analysis after 4 weeks.To target the neuroendocrine neurons in the PVN, an intraperitoneal injection of 2% (w/v) FG (Santa Cruz Biotechnology) was carried out with a survival time of 1 week.
9% saline, followed by 200 mL of cold 4% paraformaldehyde in a phosphate-buffered (0.1 M, pH 7.4) solution.The brains were collected and left in 4% paraformaldehyde for postfixation overnight and then transferred to 30% sucrose until the tissue sank to the bottom of the container.The tissue was embedded in Tissue-Tek O.C.T. Compound and sectioned at 35 m using a Leica Biosystems Cryostat.