Whole-brain connectivity atlas of glutamatergic and GABAergic neurons in mouse dorsal and median raphe nucleus

The dorsal raphe nucleus (DR) and median raphe nucleus (MR) contain populations of glutamatergic and GABAergic neurons regulating diverse behavioral functions. Their whole-brain input-output circuits remain incompletely understood. We used viral tracing combined with fluorescence micro-optical sectioning tomography to generate a comprehensive whole-brain atlas of inputs and outputs of glutamatergic and GABAergic neurons in the DR and MR. We discovered that these neurons receive inputs from similar upstream brain regions. The glutamatergic and GABAergic neurons in the same raphe nucleus have divergent projection patterns with differences in critical brain regions. Specifically, MR glutamatergic neurons project to the lateral habenula via multiple pathways. Correlation and cluster analysis indicated that glutamatergic and GABAergic neurons in the same raphe nucleus receive inputs from heterogeneous neurons in upstream brain regions and send different collateral projections. This connectivity atlas provides insights into the cell heterogeneity, anatomical connectivity and behavioral functions of the raphe nucleus.


Introduction 14
The dorsal raphe nucleus (DR) and median raphe nucleus (MR) are important modulatory 15 centers involved in a multitude of functions (Domonkos et al.,2016;Huang et al., 2019;Szőnyi et 16 al., 2019). They are implicated to have different and even antagonistic roles in the regulation of 17 specific functions, such as emotional behavior, social behavior and aggression (Balázsfi et al., 18 2018;Ohmura et al., 2020;Teissier et al., 2015). The diverse regulatory processes are related to 19 the connectivity of heterogeneous raphe neuron groups (Muzerelle et al., 2016;Nectow et al., 20 2017;Schneeberger et al., 2019). Deciphering precise neural circuits regarding the input and 21 output circuits of cell-type-specific neurons in the raphe nucleus is fundamental for understanding 22 their specific functions. 23 The DR and MR are heterogeneous and contain diverse types of neurons, including a large 24 proportion of glutamatergic and GABAergic neurons (Huang et al., 2019;Pinto et al., 2019;Sos et 25 al., 2017). Several studies have revealed that they are involved in specific functions. In the DR, 26 glutamatergic neurons play an important role in reward processing (McDevitt et al., 2014;Liu et 27 al., 2014), while GABAergic neurons are involved in regulating energy expenditure 28 (Schneeberger et al., 2019), and they have opposite effects on feeding (Nectow et al., 2017). In the 29 lines and virus tracing make it possible to label the whole-brain long-range connectivity of 48 specific neurons (Callaway et al., 2015;Huang et al., 2013;Wickersham et al., 2007). Several 49 studies have revealed a portion of the long-range connections of glutamatergic and GABAergic 50 neurons in the DR and MR through viral tracing techniques. For instance, DR GABAergic 51 neurons receive vast inputs and project to the dorsomedial nucleus of the hypothalamus (DMH) 52 and bed nuclei of the stria terminalis (BST) (Schneeberger et al., 2019;Weissbourd et al., 2014). 53 And MR glutamatergic neurons are innervated by certain aversion/fear or memory-related areas, 54 such as the LH, and they innervate the LH, medial ventral tegmental area, medial septum and the 55 vertical limbs of the diagonal bands of Broca (Szőnyi et al., 2019). Nevertheless, there is still a 56 lack of whole-brain quantitative results and comprehensive analysis of the input-output circuits of 57 glutamatergic and GABAergic neurons in the DR and MR. Furthermore, precise characterization 58 7 hypothalamus to DR and MR glutamatergic and GABAergic neurons, in which the LHA provided 145 the largest proportion, followed by the hypothalamic medial zone (MEZ), periventricular region 146 (PVR), zona incerta (ZI) and lateral preoptic area (LPO) (Figure 2C). The ZI provided more 147 inputs to MR glutamatergic neurons than to MR GABAergic neurons and DR glutamatergic 148 neurons ( Figure 3B,C). The LPO preferentially innervated MR GABAergic neurons in 149 comparison with DR GABAergic neurons and MR glutamatergic neurons (Figure 3B,D). contributed more inputs to MR glutamatergic neurons and DR GABAergic neurons than to MR 157 GABAergic neurons (Figure 3B,D). 158 The pons contributed dense inputs, in which the pons, motor related (P-mot) and pontine 159 reticular nucleus (PRNr) preferentially innervated MR glutamatergic neurons in comparison with 160 MR GABAergic neurons and DR glutamatergic neurons (Figure 3B,C). However, the pons, 161 sensory related (P-sen) provided more inputs to DR glutamatergic neurons than to MR 162 glutamatergic neurons ( Figure 3C). Moreover, the medulla, motor related (MY-mot) 163 preferentially provided inputs to MR glutamatergic neurons in comparison with MR GABAergic 164 neurons and DR glutamatergic neurons ( Figure 3B,C). These results indicated that the 165 glutamatergic and GABAergic neurons in the raphe nucleus receive inputs from similar upstream 166 brain regions with quantitative differences in certain brain regions. 167

Whole-brain outputs of glutamatergic and GABAergic neurons in the DR and MR 168
To systematically map whole-brain outputs of glutamatergic and GABAergic neurons in the DR 169 and MR, we stereotaxically injected Cre-dependent AAV-DIO-EYFP into the DR or MR in 170 Vglut2-Cre and Gad2-Cre mice (n=4 per group). The virus-labeled and GMA resin-embedded 171 samples were imaged using fMOST system ( Figure 4A). To generate quantified outputs of 172 glutamatergic and GABAergic neurons in the DR and MR, we registered the high-resolution 173 8 whole-brain image datasets to Allen CCFv3 and segmented the projection signal to calculate the 174 proportion of projection signal across brain regions (Figure 4A 29.9% of the total projections from MR GABAergic neurons but only 4.0% from MR 204 glutamatergic neurons (Figure 4B,E). 205

Habenula-raphe circuits 206
The habenula, consisting of the medial habenula (MH) and lateral habenula (LH), appears to be 207 a node connecting the forebrain and midbrain regions that are related to emotional behaviors 208 (Hikosaka, 2010). The LH has been closely connected to the DR and MR both anatomically and 209 functionally, and their connections are involved in aversion-related behavior and depression (Hu 210 et al., 2020;Zhao et al., 2015). The LH provided dense inputs to glutamatergic and GABAergic 211 neurons in the DR and MR with a preference to MR neurons in comparison with corresponding 212 DR neurons (Figure 3C, D; Figure 5A). And the input neurons to MR glutamatergic and 213 GABAergic neurons were assembled more caudally ( Figure 5A). Specifically, as MR GABAergic 214 neurons received more inputs from the LH than MR glutamatergic neurons ( Figure 3D), we found 215 that the lateral part of LH sent dense inputs to MR GABAergic neurons but sparser inputs to MR 216 glutamatergic neurons and that the input neurons to MR GABAergic neurons were assembled 217 more laterally than the input neurons to MR glutamatergic neurons on the whole ( Figure 5A). to the LH. The LH has been revealed to play a critical role in aversion and depression (Cui et al., 229 2018;Hu et al., 2020;Yang et al., 2018), and previous studies discovered that MR Vglut2+ 230 neurons could activate the LH and that activation of MR Vglut2+ neurons could induce aversive 231 10 behaviors and depressive symptoms (Szőnyi et al., 2019). These results highlight the importance 232 of the structural characteristics of the MR-LH pathway for their function roles. 233 Although the MH sent few inputs to MR glutamatergic and GABAergic neurons, and scarce 234 inputs to the DR, it has been thought to strongly project to the IPN (Lima, et al., 2017;Qin et al., 235 2009). In our results, the IPN had remarkable inputs to MR glutamatergic and GABAergic 236 neurons with very sparse inputs to DR glutamatergic and GABAergic neurons ( Figure 3F) neurons receive massive inputs from the isocortex, striatum and medulla but sparsely project to 255 these regions ( Figure 2C; Figure 4B). 256 There were massive reciprocal connections of glutamatergic and GABAergic neurons in the DR 257 and MR, which implied the feedback regulation of specific functions. To assess the reciprocity, 258 we calculated the ratio of the proportion of outputs to the proportion of inputs for each brain 259 region (Figure 6-figure supplement 1B). Approximately 45% brain regions had ratio value 260 11 between 0.25 and 4, indicating relatively balanced reciprocal connectivity. In contrast, 261 approximately 45% brain region had input bias (ratio value<0.25), and only a few brain regions 262 showed output bias (ratio value>4). Particularly, for MR GABAergic neurons, the IPN contributed 263 2.5% of all inputs but received 29.9% of outputs. 264 The glutamatergic and GABAergic neurons in the DR and MR receive a vast range of inputs, 265 whereas the relationship between different upstream brain regions needs to be explored. And the 266 axons of neurons have collateral branches targeting different areas, but the relationship of these 267 targets also needs further exploration. Based on the fact that DR and MR were heterogeneous and 268 each injection only labeled a part of neurons that might have different connectivity, we performed 269 correlation analysis and hierarchical cluster analysis to explore the similarities and variances of 270 inputs and outputs of brain regions connected with glutamatergic and GABAergic neurons in the 271 DR and MR. We selected 14 brain regions that have close long-range connections with 272 glutamatergic and GABAergic neurons in the DR and MR for analysis. As a result, the clusters 273 were not completely consistent regarding the inputs or outputs of glutamatergic and GABAergic 274 neurons in a particular raphe nucleus (Figure 6A,B). Regarding inputs to DR glutamatergic and 275 GABAergic neurons, upstream brain regions formed clear clusters, but the clusters were not 276 consistent, where only the substantia nigra, reticular part (SNr), substantia nigra, compact part 277 (SNc) and VTA formed the same cluster. Regarding outputs of DR glutamatergic neurons, the 278 MEZ, LPO, LHA, ZI, SNr and IPN formed a cluster that had very high correlation. In contrast, for 279 outputs of DR GABAergic neurons, different clusters were presented. These results implied that 280 DR glutamatergic and GABAergic neurons might have different collateral projection patterns. In 281 addition, a pair of brain regions might display opposing correlations regarding the inputs or 282 outputs of glutamatergic and GABAergic neurons in the DR. For instance, the substantia 283 innominate (SI) and NDB were positively correlated for inputs to DR glutamatergic neurons but 284 negatively correlated for inputs to DR GABAergic neurons (Figure 6A), which suggested that DR 285 glutamatergic and GABAergic neurons might receive distinct inputs from basal forebrain. 286 Regarding inputs and outputs of MR glutamatergic neurons, clear clusters formed. Specifically, as 287 for the clusters of inputs to MR glutamatergic neurons, the ZI was separate from other regions, 288 which indicated the ZI-MR glutamatergic neurons pathway might execute special functions 289 12 different from other upstream brain regions. Regarding inputs and outputs of MR GABAergic 290 neurons, no obvious clusters formed, which implied that MR GABAergic neurons might have few 291 collateral projections to these brain regions. And there were also pairs of brain regions displaying 292 opposing correlations regarding the inputs or outputs of MR glutamatergic and GABAergic 293 neurons, such as the SNc and IPN were positively correlated for inputs to MR glutamatergic 294 neurons but negatively correlated for inputs to MR GABAergic neurons ( Figure 6B). Overall,295 these results implied that glutamatergic and GABAergic neurons in the DR and MR might receive 296 inputs from and project to various unions of brain regions. 297

Discussion 298
In this study, we used virus tracing and whole-brain high-resolution imaging to generate a 299 comprehensive whole-brain atlas of inputs and outputs of glutamatergic and GABAergic neurons 300 in the DR and MR (Figure 7A,B). And systematic quantitative analysis has been performed to 301 elucidate the convergences and divergences in input and output patterns. 302 We found that glutamatergic and GABAergic neurons in the DR and MR receive inputs from 303 similar upstream brain regions but send distinctive outputs. There were essential differences 304 between the connections of DR and MR neurons. DR glutamatergic and GABAergic neurons had 305 close connections with the CEA, while MR glutamatergic and GABAergic neurons received few 306 inputs from CEA and did not project to the CEA (Figure 6-figure supplement 1A). And the 307 CEA has been revealed to regulate reward and food intake (Carter, et al., 2013;Janak et al., 2015;308 Zséli et al., 2018) that are also regulated by DR glutamatergic and GABAergic neurons (Nectow 309 et al., 2017). The IPN had dense reciprocal connections with MR glutamatergic and GABAergic 310 neurons, but almost no direct connections with DR glutamatergic and GABAergic neurons 311 (Figure 6-figure supplement 1A). The IPN and MR are both considered as important parts of 312 the midline network involved in regulating hippocampal theta rhythm (Lima et al., 2017). The 313 consistency between anatomical connectivity and behavioral function indicated the significance of 314 dissecting whole-brain connectivity for elucidating functions. 315 The differences between connections of glutamatergic and GABAergic neurons in the same 316 raphe nucleus might provide insight into their functions. DR GABAergic neurons preferentially 317 projected to the CEA, unlike DR glutamatergic neurons (Figure 4-figure supplement 1A,D). 318 13 Optogenetic activation of the CEA inhibits food intake (Carter et al., 2015), while activation of 319 DR GABAergic neurons increases food intake (Nectow et al., 2017). Furthermore, DR 320 GABAergic neurons uniquely innervate the PVT (Figure 4-figure supplement 1C,D), which is 321 connected with the CEA and involved in inhibiting food intake (Kirouac, 2015). Thus, we 322 speculated that activation of DR GABAergic neurons might inhibit CEA and PVT neurons to 323 increase food intake. DR GABAergic neurons send considerable projections to DMH, while DR 324 glutamatergic neurons scarcely project to the DMH (Figure 4-figure supplement 1B,D). It should be noted that a pair of brain regions might display different or even opposing 339 correlations regarding the inputs to glutamatergic and GABAergic neurons in the same raphe 340 nucleus (Figure 6 A,B), which implied that different neuron groups in upstream brain regions 341 might individually target heterogeneous raphe neuron groups. More advanced labeling methods 342 that could label upstream inputs to different cell types in a single brain sample is needed to 343 explore this problem further. The target regions of projections of glutamatergic and GABAergic 344 neurons in the same raphe nucleus also showed different correlations (Figure 6 A,B). This result 345 indicated that glutamatergic and GABAergic neurons in the same raphe nucleus might have 346 different collateral projection patterns, which are worth illustrating through complete single 347 14 neuron reconstruction. 348 Our results were in accordance with the input and output circuits of neurons in DR and MR 349 identified using classic tracing techniques (Marcinkiewicz et al., 1989;Oh, et al., 2014;Peyron et 350 al., 1997;Vertes et al., 2008). However, the DR and MR are heterogeneous and contain diverse 351 types of neurons, including glutamatergic, GABAergic, serotonergic and dopaminergic neurons. In summary, we constructed a comprehensive whole-brain atlas of inputs and outputs of 383 glutamatergic and GABAergic neurons in the DR and MR, revealing similar input patterns but 384 divergent projection patterns. The differences in connectivity patterns are related to specific 385 regulatory processes of specific functions. As the whole-brain connections of genetically targeted 386 neurons are key factors in characterizing cell types, our result would contribute to generating 387 whole-brain cell atlases that are under ongoing effort. Our work could form the foundation for 388 exploring the relationship of cell heterogeneity, anatomical connectivity and behavior function of 389 the raphe nucleus. 390

Materials and Methods 391
Animals 392 In this study, adult Vglut2-Cre mice (stock number: 016963) and Gad2-Cre mice (stock 393 number: 010802) purchased from the Jackson Laboratory were used. All mice were housed in an 394 experiment environment with 12-h light/dark cycle, 22 ± 1 °C temperature, 55 ± 5% humidity and 395 food and water ad libitum. All animal experiments were approved by the Institutional Animal 396 Ethics Committee of Huazhong University of Science and Technology and were conducted in 397 accordance with relevant guidelines. 398

Imaging and image preprocessing 421
For whole-brain high-resolution imaging, the virus-labeled and GMA resin embedded samples 422 were imaged with propidium iodide (PI) simultaneously staining cytoarchitecture landmarks using 423 our home-made fMOST system at a resolution of 0.32 μm× 0.32 μm × 2 μm. The acquired two-424 channel raw data are processed by mosaic stitching and illumination correction to piece together 425 into entire coronal sections using the a previously described workflow (Gong et al., 2016). Each 426 channel dataset of one brain sample contains approximately 5,500 coronal slices. For starter cells, 427 the samples were sectioned in 50 μm coronal slices using the vibrating slicer (VT1200S, Leica) 428 and imaged using the automated slide scanner (VS120 Virtual Slide, Olympus). 429

Data processing 430
Registration. To quantify and integrate the whole-brain connections, the coordinates of the soma 431 of input neurons and high-resolution image stack of labeled outputs were registered to Allen 432 CCFv3 using the transformation parameters acquired by the previously described methods (Ni et 433 al., 2020). In brief, we segmented several brain regions as landmarks through cytoarchitecture 434 17 references, such as the outline, caudoputamen, medial habenula, lateral ventricle, third ventricle. 435 Based on these landmarks, we performed affine transformation and symmetric image 436 normalization in Advanced Normalization Tools (ANTS) to acquire transformation parameters. 437 Nomenclature of brain regions. Demarcation and annotation of brain regions were based on 438 Allen CCFv3. The superior central nucleus raphe (CS) corresponds to the median raphe nucleus 439 (MR) with the consultation of the mouse brain atlas by Paxinos and Franklin (Paxinos and 440 Franklin, 2012). Based on Allen CCFv3's hierarchy of brain regions, as there are no or few input 441 neurons and projections in many brain regions, we collapsed some brain regions to their "parent" 442 region as needed, thereby divided the whole-brain into 117 brain regions (see Supplementary File 443 1) and identified 71 brain regions for analysis (areas that have small proportion of connections are 444 merged into "Others"). The STR-NA, PAL-NA, TH-NA, HY-NA, MB-NA, P-NA and MY-NA 445 respectively refer to the non-annotated area in the striatum, pallidum, thalamus, hypothalamus, 446 midbrain, pons and medulla. 447

Detection and quantification of whole-brain inputs.
Regarding the input circuits, we 448 automatically identified and localized the soma of input neurons using NeuroGPS (Quan et al., 449 2013) and manually checked the results to eliminate some indiscernible mistakes, then warped the 450 soma coordinates to Allen CCFv3 using the transformation parameters from registration described 451 above. We calculated the number and proportion of input neurons in each brain region of interest 452 (excluding target area) to generate the quantified whole-brain inputs. 453 Detection and quantification of whole-brain outputs. Regarding the output circuits, we generated 454 quantified whole-brain outputs by taking following steps: 455 We resampled the image stack of labeled neural structures to isotropic 1 μm, segmented the 456 outline of brain and set the intensity of pixels outside the outline to 0, then used the transformation 457 parameters of registration described above to warp them to Allen CCFv3 at 1 μm scaling. Then we 458 manually segmented injection sites on registered coronal sections. 459 To detect projection signal from background, each registered coronal section was background 460 subtracted, Gaussian filtered, and threshold segmented to binary image. The background image 461 * I was calculated as * min( , ) I I backgrou n d  followed by ten convolutions with the 462 averaging template of 9×9 size, where I is the gray level of coronal section and the backgroud is 463 18 an approximate estimated background intensity (Quan et al., 2013). The size of gaussian filter was 464 5×5. The filtered image was binarized by * max(4 , ) I threshold where threshold was the 465 result calculated by the Yen method that clipped to the predetermined threshold range (Yen et al., 466 1995). 467 The whole-brain images were divided into 10×10×10 μm 3 grids. In each division, we calculated 468 signal density by the definition of the sum of detected pixels divided by the sum of all pixels in a 469 three-dimensional grid, therefore generated a three-dimensional signal density matrix of 10 μm 470 voxel resolution. Then we calculated the computational path based on the signal density matrix 471 using multistencils fast marching algorithm and removed the voxels that could not back-track to 472 injection site or back-track to injection site with low confidence (Oh et al., 2014;Hassouna and 473 Farag, 2007;Liu, et al., 2018). The confidence of the path was defined as the proportion of back-474 tracking points that were located in the foreground voxel in the path, and the foreground voxel 475 refers to the voxel whose signal density was greater than a threshold. Finally, we manually 476 inspected the results and removed the remaining confusing noise voxels. 477 The outputs were quantified as projection signal volume in each brain region normalized by 478 signal volume across whole brain (excluding the injection site and target area). As soma and 479 dendrites of labeled neurons contributed a lot of signals in the injection site, we excluded the 480 injection site for more accuracy. 481

Visualization and statistical analysis 482
The Amira software (v6.1.1, FEI) and Imaris software (v9.5.0, Bitplane) were used to visualize 483 the inputs and outputs of glutamatergic and GABAergic neurons in the DR and MR. To compare 484 the inputs to glutamatergic and GABAergic neurons in the DR and MR across brain regions, we 485 performed one-way ANOVA followed by multiple comparisons with Tukey's test. To compare 486 the outputs of glutamatergic and GABAergic neurons in the same nucleus across brain regions, we 487 performed one-way ANOVA. To quantify the similarities of input and output patterns, we 488 calculated Pearson's correlation coefficients. To explore the similarities and variances of 489 inputs/outputs of brain regions connected with glutamatergic and GABAergic neurons in the DR 490 and MR, we performed correlation analysis and hierarchical cluster analysis. These processes 491 were performed using MATLAB (v2017a, MathWorks) and Python 3.6.4. To compare the whole-492 19 brain inputs across all samples, hierarchical clustering and bootstrapping were performed using 493 pvclust that is a package of R (Suzuki and Shimodaira, 2006). All histograms were generated 494 using GraphPad Prism (v.6.0, GraphPad). 495