Decoding the olfactory map: targeted transcriptomics link olfactory receptors to glomeruli


 Sensory processing in olfactory systems is organized across olfactory bulb glomeruli, wherein axons of peripheral sensory neurons expressing the same olfactory receptor co-terminate to transmit receptor-specific activity to central neurons. Understanding how receptors map to glomeruli is therefore critical to understanding olfaction. High-throughput spatial transcriptomics is a rapidly advancing field, but low-abundance olfactory receptor expression within glomeruli has previously precluded high-throughput mapping of receptors to glomeruli. Here we combined sequential sectioning along the anteroposterior, dorsoventral, and mediolateral axes with target capture enrichment sequencing to overcome low-abundance target expression. This strategy allowed us to spatially map 86% of olfactory receptors across the olfactory bulb and uncover a relationship between OR sequence and glomerular position.

To date, a systematic 3D model of OR glomerular positions has yet to be made for a 225 majority of the OR repertoire, but our high-throughput sequencing of spatial sections along 226 the 3 cardinal axes may be combined to achieve such a result. UMAP plots of the mean 227 position of ORs in each dimensional replicate indicated that ORs that are clustered on the OB 228 surface are also clustered within the UMAP positions ( Fig. 5A and fig. S10, A and B). 229 To account for the location of OR glomeruli on the outer surface of the OB, we 230 extracted coordinates from a diffusion tensor imaging (DTI) model of the mouse OB to 231 represent the approximate glomerular layer that would have been sampled by each section 232 along each dimension of our targeted transcriptomics data 61 . We based our 3D reconstruction 233 in a statistical approach to integrate our orthogonal sequencing data to determine which 234 voxels on the OB surface where the probability for a given OR is highest 62,63 . We then 235 were consistent with expectations based on OE zone, OR class, and our single-dimension 245 data (dorsal vs ventral P = 9.188 x 10 -98 , Class I vs Class II P = 5.96 x 10 -14 , TAAR vs OR P = 246 2.682 x 10 -6 , functional imaging surface enriched vs not-enriched P = 1.75 x 10 -24 (see 247 below), Mann-Whitney U-test) (Fig. 5, B, C and F and fig. S11, B and C). Additionally, the 248 distribution of predictions for the sets of Class I ORs and TAARs matched previously 249 published findings for target domains (Fig. 5C and fig. 11B). Our current model predicts 250 glomerular positions for 700 ORs and 9 TAARs, with predicted positions outperforming 251 randomly selected ORs from the same DV zone and medial/lateral side for the subset of ORs 252 with known positions with a median error of 141 μm (Fig. 5, D and E) 8 . The relative positions 253 of the Olfr1377 and Olfr881 glomeruli in the gene-targeted mouse lines (see below) were also 254 consistent with those predicted from the spatial transcriptomic data (Fig. 5D and figs. S11D, 255 S12D, S13 and S14), providing further support for our three-dimensional glomerular 256 predictions based on the spatial transcriptomic data. In summary, we found our three-257 dimensional reconstruction of OR glomerular positions to be both in agreement with 258 established OR spatial features and to show greater concordance with these established 259 features than the sets of single-dimension target capture sequencing data alone. Collectively, 260 our results thus provide the first large-scale, unified, and systematic database of OR 261 glomerular positions for the mouse OB. 262 263

Identification of ORs within the dorsal functional imaging window 264
We sought to validate and extend our findings by examining specific ORs that map to 265 glomeruli on the dorsal-central OB surface, which has been extensively characterized by 266 functional imaging in vivo 9,11,12,46,64 . To define the set of ORs accessible under standard 267 functional imaging approaches, we collected tissue samples from OBs from C57BL6 mice (2 268 male, 2 female). Each OB was dissected into two parts, one containing the approximate 269 dorsal-central imaging area and the other containing the remainder of the OB ( fig. S12A). 270 These 16 samples were processed for target capture sequencing and differential expression 271 analysis to identify ORs enriched in the functional imaging area. 272 A total of 121 ORs, including 27 Class I ORs and 94 Class II ORs, were consistently 273 enriched in the imaging surface (FDR ≤ 0.05 and LogFC > 0) (Fig. 5B, Fig. S12, B, E and F), 274 with 96% of these ORs known to localize in dorsal OE zones 38 . We also found nine of the 15 275 Olfr1377 and Olfr881 exhibited robust responses to low concentrations of 4-312 methylacetophenone, with Olfr1377 exhibiting a stronger response than Olfr881 (Fig. 6, B and 313 C). In addition to 4-methylacetophenone, we also tested a large odorant panel including 314 multiple cyclic ketones structurally related to 4-methylacetophenone, as well as more diverse 315 odorants, all at relatively low concentrations. From this panel, we identified multiple new, high-316 affinity ligands for each OR, including many cyclic ketones, and ultimately uncovered 317 overlapping but distinct response spectra for Olfr1377 and Olfr881 (Fig. 6C, fig. S15C). For 318 example, Olfr1377 showed strong responses to p-anisaldehyde, acetophenone, and the 319 aliphatic ketone 4-methyl-3-penten-2-one, while Olfr881 proved unresponsive to these 320 odorants. Here, we demonstrate a unique application of target capture sequencing, enriching low-350 abundance target transcripts found in the axon terminals of OSNs and comparing post-351 enrichment abundances across samples. Our high-throughput approach enabled us to 352 generate the first repertoire-scale dataset for OR glomeruli positions. This approach also 353 allowed us to create a three-dimensional positional estimate for the glomeruli of a majority of 354 the OR repertoire with precision within the ranges previously observed for positional variability 355 between animals 5-8, 23 . 356 With a systematic, anatomical mapping of OR identities to OB glomeruli, we can begin to 357 understand which features of sensory information impact organizational strategies utilized 358 within the mouse OB, the central circuit of olfaction. We identified functional odor-OR 359 relationships at the peripheral interface using the in vivo pS6-IP RNA-Seq deorphanization 360 approach. By identifying the set of ORs within functionally imaged glomeruli future studies will 361 be able to determine how OB circuits transform peripheral input into neural representations of 362 odors at the OR level. Additional gene targeted mouse lines amenable to functional imaging 363 will also enable us to better understand interglomerular connectivity and how patterns of 364 inhibition are established between nearby glomeruli. The combination of our newly generated 365 OR-OB positional data with ligand-receptor deorphanization data will therefore enable large-366 scale investigations into how the odor space relates to topographical features in the OB and 367 will facilitate the development of unified, systematic models for how odorants are processed, 368 interpreted, and transformed into behavior. 369 370

ACKNOWLEDGMENTS: 371
We thank members of the Matsunami lab for thoughtful discussions and feedback. Luis R.   Statistic is Mann-Whitney U-test. 691

DATA AVAILABILITY 709
All sequencing data will be made publicly available under NCBI BioProject 710 PRJNA773191 upon publication. Additional raw and processed data that support the findings 711 of this study are available from the corresponding authors upon request. 712

CODE AVAILABILITY 714
Key code for the project has been made publicly available at 715 https://github.com/kanazian/obmap. 716

Statistics 719
All statistical tests were performed in R. All Mann-Whitney U-tests were two-tailed and 720 performed with continuity corrections using the wilcox.test function from the stats package. before placing on ice. Once solidified, the mold was removed and the agarose-embedded 742 brain was prepared for vibratome sectioning by cutting away agarose to leave a triangular 743 shape with the tip forming at the bulb on the surface that will be cut first for that specific 744 direction. The triangular block was glued (Loctite 404) onto the vibratome stand (Leica 745 VT1000S with Feather carbon blades), submersed in cold 1X HBSS (Gibco) and 100 μm 746 sections serially cut. Sections were placed into 1.5mL tubes and stored at -80C. 400 µL of 747 Buffer RLT (Qiagen) with 10% BME was added to each frozen tissue section before 748 homogenization at 30,000RPM for 4 seconds using a mounted Biogen PRO200 with 5mm flat 749 tip generator probe. The homogenizer probe was rinsed three times with DI water between 750 samples. RNA was extracted using a Qiagen RNeasy kit according to the manufacturer's 751 instructions. cDNA synthesis was performed using a SMART-Seq v4 (Takara) kit with 10 ng 752 total RNA used for input in half-sized reactions. 10 ng cDNA was used in the KAPA Hyperplus 753 library construction kit following the manufacturer's instructions. All quantifications were 754 performed using Qubit.          Table S8: Three-dimensional model posterior medians 1133 Table S9: Differential expression analysis for functional imaging surface and pS6-IP 4-1134 methylacetophenone studies 1135