Authors: David H. Brann (1Department of Neurobiology, Harvard Medical School, Boston, MA 02115, USA), Tatsuya Tsukahara (2Lunenfeld-Tanenbaum Research Institute, Sinai Health System, Toronto, ON M5G 1X5 Canada; 3Department of Molecular Genetics, University of Toronto, Toronto, ON M5S 1A8 Canada), Cyrus Tau (1Department of Neurobiology, Harvard Medical School, Boston, MA 02115, USA), Dennis Kalloor (1Department of Neurobiology, Harvard Medical School, Boston, MA 02115, USA), Rylin Lubash (1Department of Neurobiology, Harvard Medical School, Boston, MA 02115, USA), Lakshanyaa Thamarai Kannan (1Department of Neurobiology, Harvard Medical School, Boston, MA 02115, USA), Nell Klimpert (4Department of Neuroscience, Division of Biology and Medicine, Brown University, Providence, RI 02912; 5Robert J. and Nancy D. Carney Institute for Brain Science, Brown University, Providence, RI 02912), Mihaly Kollo (6Sensory Circuits and Neurotechnology Laboratory, The Francis Crick Institute, London, NW11AT UK), Martín Escamilla-Del-Arenal (7Department of Biochemistry and Molecular Biophysics, Vagelos College of Physicians and Surgeons, Columbia University, New York, NY 10032, USA; 8Mortimer B. Zuckerman Mind, Brain, and Behavior Institute, Columbia University, New York, NY 10027, USA), Bogdan Bintu (9Department of Cellular and Molecular Medicine, University of California San Diego, San Diego, CA, 92093, USA), Andreas Schaefer (6Sensory Circuits and Neurotechnology Laboratory, The Francis Crick Institute, London, NW11AT UK), Alexander Fleischmann (4Department of Neuroscience, Division of Biology and Medicine, Brown University, Providence, RI 02912; 5Robert J. and Nancy D. Carney Institute for Brain Science, Brown University, Providence, RI 02912), Thomas Bozza (10Department of Neurobiology, Northwestern University, Evanston, IL 60208, USA), Sandeep Robert Datta (1Department of Neurobiology, Harvard Medical School, Boston, MA 02115, USA; 11Lead contact)
Categories: Article
Source: Cell
Authors: David H. Brann, Tatsuya Tsukahara, Cyrus Tau, Dennis Kalloor, Rylin Lubash, Lakshanyaa Thamarai Kannan, Nell Klimpert, Mihaly Kollo, Martín Escamilla-Del-Arenal, Bogdan Bintu, Andreas Schaefer, Alexander Fleischmann, Thomas Bozza, Sandeep Robert Datta
Although topographical maps organize many peripheral sensory systems, mouse olfactory sensory neurons (OSNs) are thought to randomly choose which one of ~1100 possible olfactory receptors (ORs) to express, with spatial organization in the olfactory epithelium limited to a handful of broad anatomical “zones” that modestly restrict OR choice. Here we reveal a stereotyped receptor map in the olfactory epithelium, in which each receptor is expressed at a unique mean dorsoventral position. OSN dorsoventral identities are encoded by a coherent gene expression program, which includes key transcription factors and axon guidance molecules; use of this program reflects a dorsoventral gradient in retinoic acid signaling, translates each physical location into a spatially appropriate distribution of potential OR choices, and aligns receptor maps in the nose and brain. Spatial order in the olfactory system therefore arises from a continuously varying transcriptional code, one that precisely organizes the many discrete channels responsible for smell.
The nervous system uses space to map sensation^1,2^. The cochlea, for example, encodes different frequencies at distinct positions along its unrolled length; this tonotopic map is serially propagated through patterned projections to the brainstem, thalamus, and cortex^3^. Like these auditory maps, a similar cascade of maps — originating in the sensory periphery and hierarchically running through higher brain centers — supports touch and vision^4,5^. In the olfactory system, odors are detected by mature olfactory sensory neurons (OSNs) in the nasal main olfactory epithelium (MOE), each of which expresses only one of the 1,172 functional olfactory receptors (ORs) encoded in the mouse genome^6,7^. However, the extent to which OR expression is spatially organized in the MOE into a stereotyped peripheral sensory map remains unclear^8–10^. Addressing this foundational question is critical for our understanding of the sense of smell, as the answer has significant implications for how OSNs choose their receptors, encode information about odors, and find their targets in the higher brain.
Two challenges have vexed our understanding of whether OSNs are spatially patterned in the epithelium. First, characterizing the relationship between physical location and gene expression is difficult in the MOE, because the underlying cartilaginous turbinates transform an otherwise smooth and continuous epithelial sensory sheet into a complex labyrinth^11^, challenging analysis of spatial order. Second, most of what we know about spatial patterning of OR expression arises from in situ hybridization and related techniques, which can query only a few genes at once and are difficult to align across samples^12–15^; although microdissection has been recently used to assign ORs to particular epithelial locations, the resulting positional estimates are indirect and the degree of inter-individual stereotypy remains unclear^16,17^.
Despite these constraints, current models propose individual receptors are assigned to one of between four and nine dorsoventral “zones” (with some positing additional subzones) but are randomly chosen by the OSNs within each zone, leading to spatially distributed expression of each OR across a broad swath of epithelium (Figure 1A)^12–16,18^. Because the MOE continuously regenerates during adulthood from a population of basally-localized stem cells, this zonal model implies that OSN precursors acquire zone-specific cell identities that limit OR choice to an appropriate subset^19^. Consistent with this possibility, OSNs are known to be segregated into a small number of fate-restricted lineages that each express a different receptor gene class^20,21^. Late-stage OSN precursors also appear fated to choose from restricted subsets of ORs, as the OR “switches” occasionally observed during OSN differentiation are largely among ORs from the same epithelial zone^22–24^. However, only a handful of marker genes uniquely identify the dorsal-most zone, and no genes have been identified that demarcate the different ventral zones^25,26^.
It also remains unclear how OSN precursors might restrict OR choice to the spatially-appropriate set of ~150–250 ORs per zone^27^. Such restrictions must ultimately influence the mechanisms through which precursors choose a single OR, which are not yet fully understood but low-level co-expression of ORs in OSN precursors, which may license ORs for potential choice; heterochromatin-mediated silencing of ORs that are not ultimately chosen; inter-chromosomal interactions between OR genes, which facilitate both OR repression and choice; and, finally, high-level expression of a single OR allele^7,28–33^. These choice mechanisms may be sensitive to zonal position, as recent results suggest that zonal positioning may influence the set of ORs available for co-expression in precursor cells^28^.
Information from OSNs is projected to insular structures called glomeruli in the olfactory bulb (OB), the first waystation for olfactory processing in the brain. All OSNs of a given subtype (i.e., those expressing the same OR) innervate one of a small number of glomeruli whose spatial location is largely invariant across animals, thereby creating a stereotyped spatial map of odor identity^34–37^. This architecture demands a means to align the proposed zones in the nose and the precise, spatially-organized map in the OB^38,39^. This mechanism likely involves the axon guidance genes Robo2 and Nrp2, which are expressed in complementary gradients in the epithelium and causally regulate the targeting of OSN axons along the dorsoventral axis in the bulb^40–43^. However, it is not known whether these gradients systematically vary across OSN subtypes, or whether this variation is sufficient to ensure proper glomerular targeting^44,45^. Furthermore, the disconnect between the broad spatial distribution of ORs in the MOE and the precision of their axonal targeting in the bulb has led to the suggestion that the OR itself determines key aspects of OSN maturation by regulating the expression of axon guidance genes and by facilitating homophilic and heterophobic interactions among axons^38,46–49^; in this model, OSN cell identities are largely derived from their late-stage expression of ORs, rather than from spatial position within the epithelium.
Here we show that the nose harbors a receptor map, one built from ~1100 distinguishable peaks of OR expression whose mean positions are reliably stereotyped across mice. This patterning is the consequence of a coherent gene expression program comprised of ~250 genes, including key transcription factors, retinoic acid signaling components, and axon guidance molecules. The relative use of this program distinguishes all mature OSN subtypes, varies systematically as a function of dorsoventral spatial position in the epithelium, and reflects retinoic acid gradients in the mesenchyme underneath the OE proper. OSN stem cells and precursors take advantage of this gene expression program — which both influences many of the known mechanistic processes that underlie OR choice and predicts the spatial location of target glomeruli — to encode information about spatial position before receptors are chosen for expression and before the process of axon guidance has begun. Thus, the transcriptional identity of OSN precursors is defined by spatially-organized gradients of gene expression, which are then transformed into discrete and predictable OR choices and glomerular targets; this mechanism coordinately organizes the first two stages of the peripheral olfactory system by constructing a detailed receptor map in the nose aligned with the precise receptor map apparent in the olfactory bulb.
We hypothesized that olfactory sensory neurons (OSNs) contain transcriptional signatures that reflect their spatial positions within the epithelium. To explore this idea, we used single cell RNA sequencing (scRNA-seq) to characterize ~5 million MOE cells across hundreds of mice, which yielded a dataset of ~2.3 million OSNs that collectively expressed ~1,100 ORs (Figure S1A). Although we expected to identify specific marker genes for each zone, our analysis failed to reveal such genes (Figure S1B)^28,50^.
However, matrix factorization identified several OSN-specific gene expression programs (GEPs), each composed of 100–200 individual genes whose expression levels co-vary across different OSNs^51,52^. Two such GEPs (provisionally named GEPDorsal~ and GEPVentral) included the known dorsal marker Nqo1 and the ventral marker Nfia/b/x, respectively; these GEPs also included transcription factors, axon guidance molecules, and cytoskeletal regulatory genes (Figure S1A and Table S1). For simplicity, we distilled the expression of the genes comprising each GEP in each OSN into a single number, the GEP “usage”. Because individual OSNs used either GEPDorsal or GEPVentral but not both, we could further summarize the use of these two GEPs by computing a dorsoventral (DV) score (i.e., the usage of GEPDorsal subtracted from that of GEPVentral, Figure 1B). Unlike what would be expected if OSNs were organized into homogeneous and discrete “zones”, Nqo1 expression tonically and smoothly increased as a function of each OSN’s DV score, as did the 100 other genes participating in GEPDorsal~ (including e.g., Robo2); conversely, the expression of Nfia/b/x and the 150 other genes participating in GEPVentral~ (including e.g., Nrp2) tonically decreased as a function of each OSN’s DV score (Figure 1C–E)^40–43^. Thus, the expression of the 250 genes that make up GEPDorsal~ and GEPVentral systematically varies in a graded fashion across OSNs.
The use of GEPDorsal and GEPVentral by a given OSN reflected its expressed odor receptor, as each OSN subtype (defined as the collection of OSNs expressing the same OR) was associated with a unique average DV score; across subtypes, average DV scores varied smoothly without obvious clustering (Figures 1F–H, S1C–F, and Table S2). The relationship between OSN subtypes and DV scores reflected structured covariation across many genes, as DV scores were well-approximated by random subsets of 10–20 DV genes and changed little after excluding dorsal and ventral marker genes (Figures S1G–H). The individual OSNs belonging to each subtype exhibited closely related (but not identical) DV scores, such that OSNs belonging to different subtypes could be effectively distinguished based upon their DV scores alone (median auROC = 0.96, Figures 1I–J and S1I). The ability to distinguish OSN subtypes depended upon the difference in the DV scores of each pair, both within and across proposed zones (Figures 1J and S1J–K). Similarly, expression levels of nearly every DV gene were better predicted by DV scores than by each subtype’s putative zonal identity (Figure S1L).
The relationship between each OSN subtype and its associated mean DV score was remarkably stable, as the rank ordering of OSN subtypes based on their DV scores was nearly invariant (ρ = 0.987) across hundreds of independent samples (Figures 1K and S1M–P). Thus, despite variation in DV scores among individual OSNs expressing the same OR, the mean DV score associated with each OSN subtype — which summarizes the coordinated expression of several hundred genes — is unique and stereotyped across mice and is sufficient to differentiate each OSN subtype from all others.
We hypothesized that the precise dorsoventral location adopted by each OSN is encoded by its DV score. If so, the mean physical location of each OSN subtype — and therefore OR — along the DV axis should be unique, predictable and stereotyped across mice. To test this alternative hypothesis to the discrete zonal model, we used multiplexed error-robust fluorescence in situ hybridization (MERFISH) to measure the spatial positions of hundreds of OR and DV genes throughout the epithelium (Figures 2A–B and S2A–B)^51^.
As expected, individual ORs were expressed in a spatially restricted manner, with the most-dorsal receptors being more broadly expressed than more-ventral receptors (Figure 2B). However, simply coloring each OSN based upon its associated patterns of DV gene expression (the DVISH score, determined directly via MERFISH) revealed a previously unappreciated fine organization to the along each segment DVISH scores smoothly varied from dorsal to ventral (Figures 2C and S2C). Furthermore, identifying the OR expressed within each OSN and coloring it based upon its DV score (as measured via scRNA-seq) also yielded a similar smooth and continuous gradient (Figures 2D and S2C). Physical space also organized OSN DV scores when the epithelium was “unrolled” by aligning all segments to a common spatial axis (Figures 2E and S2D). Given the stereotyped association between each OSN subtype and its mean DV score, these observations reveal that each OR occupies a unique mean position along the DV axis of the epithelium.
The smoothly-varying relationships between epithelial position, DV scores, and OR expression were conserved across animals and enabled accurate predictions of the spatial location of each OSN subtype (and hence OR) in our MERFISH data, in a separate spatial transcriptomics dataset generated via Stereo-seq, and in five independently obtained datasets in which the position of OSN subtypes could be precisely or approximately assigned (Figures 2F–G and S2E–H)^16,17,28,52^. In contrast, we did not observe any discrete spatial clustering of ORs into clear zones, and DV scores outperformed zonal models at predicting OSN spatial positions (Figures S3A–B).
While ventral OSN subtypes had distinct and predictable mean locations, their constituent OSNs occupied a distribution of locations (with a sharp peak and long tails encompassing 5–15% of the epithelium); the degree of overlap in DV score distributions between OSN subtypes predicted their extent of spatial overlap (Figures 2H–I and S3C). The relationship between spatial positions and DV scores was less apparent in the dorsal-most region (Figures 2D–E and S3D), likely reflecting the broader spatial distribution observed for dorsal receptors. Importantly, DVISH~ scores accurately predicted spatial positions both within and across OSNs of a given subtype (Figures S3E–F), suggesting that DV score variation largely reflects differences in physical location rather than noise.
Together, these data reveal that the precise location of each OSN is well-captured by its DV score and that each OSN subtype occupies a unique and restricted spatial distribution in the olfactory epithelium, thereby creating a receptor map with ~1100 overlapping but distinguishable peaks. This pattern is the consequence of each OR having a preferred epithelial spatial position, and of those OSNs residing at each spatial position having a preferred OR.
Our factorization approach identified an additional pair of GEPs (together comprised of 300 total genes; see Table S3) which we provisionally named GEPanterior~ and GEPposterior, as they included genes (like Nrp1 and Plxna1) known to mediate axon targeting along the bulb anteroposterior (AP) axis^46,47,53,54^. Like the DV GEPs, OSNs used either GEPanterior or GEPposterior (Figure S3G), enabling the generation of an “AP score” that varied across subtypes. We wondered whether this variation was also driven by epithelial location, and via MERFISH observed that differences in AP scores across OSN subtypes corresponded to their apical-basal position within the OSN layer (Figures 2J–K, and S3H). This relationship was strongest for dorsal OSN subtypes (Figure 2L), which notably are the subset of OSNs in which DV scores are least informative about spatial position.
Prior work suggests that OR-specific spontaneous activity drives the expression of some (and possibly all) AP GEP genes^54^. We therefore wondered whether the DV score of each OSN subtype was similarly determined by the expressed OR. To test this idea, we performed scRNA-seq on OMP-P2 mice, in which the P2 receptor is ectopically expressed across the dorsoventral axis; this process frequently silences the previously-chosen OR in each OSN (Figure S4A)^55^. P2 OSNs from OMP-P2 mice adopted a much wider distribution of DV scores than did P2-expressing OSNs from wild-type mice (Figures S4B); this result demonstrates that the identity (and associated activity) of the expressed OR does not determine OSN DV scores, suggesting instead that DV scores are determined by OSN spatial position. Consistent with this alternative hypothesis, the ~15% of OSNs in which P2 was co-expressed with another OR (i.e., those switching from their original OR to P2) exhibited DV scores that matched the originally-chosen OR rather than P2 (Figure S4C–D). Odor-evoked activity downstream of the OR also appeared unlikely to determine DV scores, which were essentially unchanged after three distinct activity naris occlusion, exposure to novel olfactory environments, and knocking out CNGA2, the key ion channel that transforms OR activation into calcium entry and spiking activity (Figures S4E–F)^54,56^. Thus, the OR expressed by each OSN does not influence its DV score, and yet under normal circumstances the identity of the expressed OR correlates with each OSN’s spatial position and DV score.
These observations suggest a model in which spatially-organized patterns of DV gene expression are initiated prior to OR expression, in either the stem cells or immediate neuronal precursors (INPs) that reside in organized layers beneath their mature progeny (Figure 3A)^57,58^. To test this hypothesis, we analyzed DV gene expression in 400,000 differentiating cells extracted from our scRNA-seq dataset (Figures 3B–D and S5A). Many DV genes were expressed during OSN differentiation, which allowed us to compute DV scores for all precursor cells (Table S1). This revealed a gradient of DV scores across OSN precursors that was present throughout all stages of differentiation and which spanned all three OSN lineages (corresponding to expression of class I ORs, dorsal class II ORs, and ventral class II ORs) apparent in our dataset (Figures 3E and S5B–E).
Principal components analysis revealed that DV score variation was a major driver of transcriptional variation among precursors within each developmental stage, though expression of individual DV genes peaked at different stages (Figures S5F–G). Differences in precursor DV gene expression reflected differences in spatial positioning within the epithelium, while at each dorsoventral location precursors and overlying OSNs exhibited similar DV scores (Figure 3F). Spatial position therefore determines the expression levels of GEPDorsal and GEPVentral genes across all stages of OSN differentiation, thereby transcriptionally diversifying OSN precursors before they commit to expressing a single OR.
The observation that precursors and nearby OSNs have similar DV scores suggests that stem cell dorsoventral position might restrict OR choice. To test this idea, we performed in vivo clonal lineage tracing by chemically ablating OSNs with methimazole and infecting basal stem cells intranasally with lentiviral DNA-barcodes. After allowing the stem cells to regenerate the epithelium for several weeks^59^, we performed scRNA-seq to identify the set of OSN subtypes (and therefore ORs) generated by each marked stem cell (Figures 3G and S6A–D).
Most identified clones (~800 out of ~1,100) included OSNs, which collectively expressed more than 750 different ORs (Figure S6E). Individual OSNs derived from the same clone often expressed different ORs, demonstrating there is not a deterministic 1 relationship between stem cell positional identity and OR choice. However, clonally-related OSNs had much more similar DV scores than did OSNs from different clones, and as a result expressed the same OR at above-chance rates (Figures 3H–J, and S6F–I). Clonally-related precursors also had similar DV scores; importantly, precursor DV gene expression could be used to predict the ORs expressed by the OSNs found in the same clone (as determined by their associated DV scores, Figures 3K and S6J). Interestingly, the DV scores of clonally-related OSNs were even more restricted than expected based on their chosen ORs, suggesting a tighter mapping between epithelial space and DV identity than between epithelial space and OR choice (Figure 3L). These observations are consistent with the process of OR choice in OSN progenitors being biased towards a spatially optimal receptor — thereby generating a patterned map of OR identity in the nose — but also being sufficiently stochastic such that OSN progenitors often pick ORs associated with OSN subtypes nearby in space.
Our results demonstrate that OSN precursors somehow “know” where they are located physically along the dorsoventral axis, enabling location-aware OR choice. Consistent with this idea, OSNs regenerated following methimazole treatment exhibited the same relationships between DV scores and OR choice observed in untreated epithelia (Figure S6C). We searched for candidate signaling programs that might inform stem cells about their dorsoventral position, and observed several genes related to retinoic acid (RA) signaling whose expression varied spatially across different epithelial cell types (Figure 4A)^60–62^. Consistent with prior reports^40,62^, retinaldehyde dehydrogenase 2 (ALDH1A2, the enzyme that converts retinaldehyde into RA) is expressed in lamina propria mesenchymal cells — below the stem cell layer — in a gradient that negatively correlated with the DV scores of nearby OSN subtypes (Figures 4B–C and S7A). OSN precursors expressed RA receptors like Rara, Rarg, Rxrb and Rere, and both precursors and OSNs expressed additional genes capable of producing or responding to RA in a graded manner including Aldh1a3, Rarb, Dhrs3, and the cellular retinoic acid binding protein 1 (Crabp1, Figures 4A, 4C and S6B–C).
Mesenchymal RA might therefore signal dorsoventral position and influence OSN fates. To test this hypothesis, we systemically delivered either an inhibitor of ALDH1A2 (WIN 18,446) or an activator of RA signaling (all-trans retinoic acid, atRA) to adult mice during methimazole-induced regeneration and evaluated their effects via scRNA-seq (Figures 4D and S7D). Compared to control mice, we observed bidirectional shifts in the number of regenerated “dorsal” or “ventral” OSNs: mice given the RA inhibitor had fewer OSNs with more-ventral DV scores (and therefore fewer cells expressing ventral ORs), whereas mice given atRA exhibited the opposite phenotype (Figures 4E and S7E); these bidirectional, RA-dependent shifts in DV scores were observed in both precursors and OSNs (Figure 4F). Importantly, these manipulations neither altered the relationship between DV scores and the chosen ORs per se, nor changed the co-variation structure of the genes that make up the DV score (Figures 4H–J). Off-target drug effects are unlikely to explain our observations, as treating mice with a panretinoic acid receptor (RAR) inverse agonist (BMS493) also led to dorsal shifts in OR choices and DV scores (Figure 4K).
These results predict that inhibiting RA signaling should shift the relative positions of OSN subtypes along the dorsoventral axis of the epithelium, which we tested by performing in situ hybridizations for a pooled set of ORs expressed in the dorsal part of the ventral epithelium (see Methods). Consistent with our scRNA-seq results, treatment with an RAR inhibitor caused more cells to express relatively-dorsal ORs; furthermore, the spatial distribution of these cells was shifted ventrally (Figures 4L–M), consistent with the ventral epithelium being dorsalized under conditions of RA inhibition.
Our data suggest that a gradient of sub-epithelial RA informs overlying stem cells and precursors about their dorsoventral position, thereby translating spatial location into a particular pattern of DV gene expression and ultimately into the selection of spatially-appropriate ORs. We therefore wondered how precursor DV gene expression might influence the multiple mechanisms that support singular OR choice, which in broad strokes involves the low-level co-expression of multiple receptors, the repression of unchosen receptors through targeted heterochromatinization and genomic compartmentalization of OR genes, and the selective high-level expression of a single OR gene^28–31,63–66^.
Pseudotemporal ordering of our scRNA-seq data confirmed that OSN precursors often co-express 10–20 ORs at low levels before they ultimately choose a single OR to express at high levels. We wondered whether dorsal precursors preferentially co-expressed dorsal ORs and conversely whether ventral precursors co-expressed ventral ORs. Surprisingly, all precursors regardless of spatial position (as inferred by the DV score of each cell) initially express a limited subset of the dorsal-most receptors (Figures 5A–C). The dorsal-most precursors subsequently express dorsal ORs at higher levels, such that dorsal ORs like Olfr1513 are exclusively expressed at high levels in dorsal cells. In contrast, more ventral precursors instead silence the weakly co-expressed dorsal ORs and start to co-express slightly more ventral ORs like Olfr194 (Figure 5B). As differentiation proceeds these precursors then similarly resolve their co-expression into a singular choice, and the process of co-expression, selection (if spatially appropriate), and down-regulation (if not) progressively repeats itself in a dorsal-to-ventral pattern until the ventral-most ORs are finally co-expressed and subsequently chosen in the ventral-most cells.
As a result of this process, there is a sliding window of OR co-expression that systematically evolves as OSN precursors differentiate, with dorsal-most receptors being co-expressed and chosen first, and ventral-most receptors being co-expressed and chosen last (Figure 5C–F). The delay between dorsal and ventral choice likely reflects the fact that ventral precursor cells must first co-express and silence more-dorsal ORs before they can express a spatially-appropriate ventral OR. Note that class I receptors located in the dorsal-most zone violate this pattern of co-expression and repression, consistent with class I OR-expressing OSNs having a distinct developmental trajectory (Figure S5D).
Do the spatial identities of each OSN precursor predict which of the co-expressed ORs is chosen for singular expression? We find that at earlier stages there is only weak alignment between the weakly co-expressed ORs and cellular DV scores, but that once precursors start expressing ORs at higher levels the single highest-expressed OR becomes well aligned to cellular DV scores (Figures 5G and S5H–I). Thus, as OSN precursors transition from expressing multiple to single ORs, expression levels of “incorrect” ORs systematically decline, yielding an OSN that on average expresses the most-appropriate OR for its given spatial location (although such OSNs still frequently choose ORs associated with nearby spatial positions). Consistent with RA gradients informing OSN precursors about spatial location, RA manipulations systematically shifted the set of ORs that were co-expressed in pre-choice cells (Figure S7F). Dorsoventral positioning therefore not only influences which of the ~1100 different ORs a given OSN chooses to express but also shapes how this process unfolds during differentiation.
OR genes not chosen for singular expression are repressed by heterochromatin, which is deposited onto OR loci during differentiation. We wondered whether regulators of this process might participate in establishing the relationships between precursor DV scores, the timing of OR choice during differentiation, and the set of ORs available for singular choice. Indeed, recent work identified the NFI transcription factors as being causally responsible for tuning the levels of OR heterochromatinization in more-ventral “zones”, as knocking out Nfi genes in OSN precursors reduces heterochromatin levels on ventral OR loci and causes OSNs in putative ventral zones to express more-dorsal ORs^28^. The Nfi genes also plausibly couple spatial identity to OR choice, as they are part of GEPVentral, are expressed in a dorsoventral gradient in both precursors and OSNs, and are bidirectionally sensitive to RA manipulations (Figures 1C, S2B, S5B, and S7G). The expression of Nfi and potentially other regulatory DV genes — which reflect dorsoventral position — may therefore directly bias the choice process in OSN precursors by regulating heterochromatin deposition in a spatially-sensitive manner.
Given these observations, we hypothesized that OR heterochromatin deposition may not vary by “zones” but instead be organized into a continuous dorsoventral gradient, one capable of repressing suboptimal ORs to allow each precursor to express spatially-appropriate ORs. To explore this idea, we reanalyzed datasets in which OSNs from the entire epithelium were pooled and the levels of heterochromatin marks on each OR locus were quantified using chromatin immunoprecipitation followed by DNA sequencing (ChIP-seq)^28,31^. If OSN differentiation requires heterochromatinization of co-expressed ORs — and the most dorsal ORs are transiently expressed and subsequently silenced in all precursors while ventral ORs are expressed only in ventral precursors — then dorsal OR loci should have relatively higher levels of heterochromatin marks in this type of pooled analysis.
Indeed, we observe a strong monotonic relationship (ρ = 0.89) between the DV score associated with each OR and the level of heterochromatin observed at its locus (Figure 5H–I). In contrast, relatively uniform levels of heterochromatin were observed across all OR loci in ChIP-seq data from OSNs expressing the ventral Olfr1507 receptor (Figure S7H), consistent with ventral OSNs silencing all more-dorsal OR loci. We also observe a continuous dorsoventral gradient of inter-chromosomal OR–OR genomic contacts, in which dorsal OR loci systematically interact more frequently with other OR loci than do ventral OR loci when analyzed in bulk across all OSNs (Figures 5J and S7I). This result implies that more-ventral OSNs have nuclei with relatively more interchromosomal OR–OR genomic contacts than do more-dorsal OSNs, likely because they have silenced more OR loci. Thus — like OR expression patterns — both heterochromatin levels on OR loci and genomic contacts between them are not organized zonally but systematically vary along a dorsoventral gradient (Figure 5D).
If heterochromatin levels regulate which OR is chosen for expression at a given spatial position, then disrupting heterochromatin might alter the relationship between spatial position and OR choice. One strategy for doing so is suggested by the finding that heterochromatinized OR loci are bound by the heterochromatin protein 1β (HP1β), which helps to compact DNA into its repressed form^67,68^. Because HP1β loss leads to neonatal lethality, we performed scRNA-seq on “HP1 swap” mice, in which HP1β loss is rescued via HP1α, which also interacts with chromatin but plays a less prominent role in transcriptional silencing (Figures S7J)^69–71^.
OSNs from HP1 swap mice were dorsalized relative to control mice and had fewer OSNs expressing ventral receptors, similar to what we previously observed after treating mice with drugs that decrease RA signaling (Figures 5K and S7K–L). However, unlike what we observed after RA manipulations, OSNs and precursors derived from HP1 swap mice of a given DV score expressed even more dorsal receptors than would otherwise be expected (Figure 5L–M and S7M–N). The fact that these manipulations alter the typical relationship between DV gene expression and OR choice demonstrates that OR heterochromatinization — the result of DV genes like Nfi — causally translates spatial position into a spatially-appropriate distribution of ORs.
While variation in heterochromatin levels on OR loci provides a compelling means for suppressing the expression of spatially-inappropriate receptors, it remains unclear how each OR is linked to its specific spatial location. One possibility is that differences in the binding and/or activity of transcription factors upstream of OR genes organize OR expression in space. Motif enrichment analysis confirmed the presence of putative binding sites for the LHX and EBF transcription factors within most OR promoters^72,73^. Consistent with the early expression of dorsal ORs during OSN differentiation, ORs with increased number of binding sites for EBF (which is required for OR expression) were more likely to be dorsal (Figure S8A). Furthermore, motifs for NFI and RAR transcription factors were enriched in the promoters of more-ventral ORs, consistent with the graded expression of both NFI and RA-related genes observed in ventral OSNs (Figure S8A).
These observations suggest that local genomic context renders each OR competent for expression at a particular dorsoventral location in the epithelium. Indeed, analysis of “receptor swap” mice, in which the coding sequence of one receptor replaces that of another receptor at a distinct genomic location^22^, revealed that OSN DV scores (and expression levels of NFI and other DV genes) reflect the genomic locus in which the expressed OR resides (Figures S8B–E).
To further explore the sufficiency of cis-regulatory elements in determining OR dorsoventral positioning^74^, we characterized DV scores in OSNs from F1 hybrid crosses between C57BL/6J mice and wild-derived CAST/EiJ mice (Figure 6A). Consistent with OR genes being expressed monoallelically, each OSN exclusively expressed either the C57 or CAST allele of its chosen OR (Figure 6B). While these two mouse strains have largely homologous sets of ORs, their genomes also harbor many single nucleotide polymorphisms (SNPs) that in principle could impact cis-regulatory elements and thereby alter OR spatial positions^75^.
Consistent with this prediction, a subset of OSN subtypes (and therefore a subset of OR alleles) exhibited significant across-strain differences in their DV scores in both F1 hybrid mice and in the parental F0 C57 and CAST strains (Figures 6C and S9A–B). We wondered whether these strain-specific differences in DV scores reliably predicted differences in the epithelial positions of individual ORs. We therefore performed in situ hybridizations for Olfr938 and Olfr916, two ORs with higher DV scores and correspondingly higher levels of heterochromatin in C57 mice compared to those in CAST mice (Figures 6D–E). Consistent with the observed differences in their DV scores, OSNs expressing Olfr938 and Olfr916 were located exclusively in the Nqo1^+^ dorsal region of the epithelium in C57 mice and in the Nqo1^−^ ventral region in CAST mice (Figures 6F and S9C). Furthermore, in situs for Olfr938 and Olfr916 in the F1 hybrids revealed that the same OR was found both within and outside of the Nqo1^+^ dorsal region — and thus effectively resided in two “zones” simultaneously — matching the spatial positions of the respective F0 animals (Figures 6F and S9C).
How might OR genomic variation lead to spatial variation between strains? Genomic SNPs are located both in the upstream and coding regions of OR genes. However, neither strain-specific mutations in OR coding regions nor the number of SNPs in upstream regions predicted whether the position of a given OR varied across strains (Figure 6G–H), suggesting that the specific set of upstream SNPs determines their influence on spatial position. Indeed, comparing the spatial pattern of expression of Olfr938 in CAST, C57 and two additional wild-derived strains (MOLF/EiJ and PWK/PhJ) with the set of SNPs present in each strain^75,76^ nominated a single upstream SNP as a potentially causal determinant of ventral expression in CAST mice (Figure S9D–F). This analysis — like our genetic manipulation and F1 hybrid experiments — argues that the spatial position of each OR is governed by non-coding elements in its genomic locus.
ORs are the product of gene duplication events, which can generate OR genes with similar promoters and primary sequences at nearby genomic locations^21,77^. This phenomenon raises the possibility that OR genomic position, phylogenetic relationships or chemical tuning properties might correlate with DV scores and therefore spatial positioning. Consistent with prior work that suggested that homologous ORs within a genomic cluster may have similar zonal identities^16,78^, we find that genomically adjacent OR loci have, on average, more-similar protein sequences (as measured by their phylogenetic similarity) and more-similar epithelial spatial positions (as captured by the DV score) than do more distant loci (Figure S10A). However, the overall relationship between phylogenetic distance and DV scores was weak (Figures S10B–C), suggesting that these effects may only hold for a small subset of ORs, such as those that were the recent products of gene duplication.
Although OR phylogeny only modestly relates to spatial position, it remains possible that OR tuning properties themselves are spatially organized in the epithelium. To address this possibility, we used Act-seq to identify the set of OSNs activated by various monomolecular odorants and to reveal their associated ORs and DV scores^54^. Monomolecular odors delivered at concentrations that activate many receptors triggered patterns of neural activation that were spatially distributed across the entire epithelium (Figure S10D). In contrast, acids largely activated class I ORs in the dorsal region of the epithelium, and odors delivered at concentrations that activated fewer ORs tended to recruit ORs that had, on average, similar DV scores (Figure S10E). However, this spatial clustering of odor tuning properties was weak and appeared to be largely driven by close phylogenetic relationships among a subset of ORs. These results indicate that ORs with highly related sequences can sometimes be both similarly tuned and located near each other in the epithelium, but that in general odors activate OSNs that are spatially distributed across the entire epithelium.
The olfactory bulb contains a precise and stereotyped map in which each OSN subtype typically projects to one medial and one lateral glomerulus. In addition to transcription factor genes like Nfi and Rar that likely influence OR choice, the DV score includes known regulators of OSN axonal targeting like Robo2 and Nrp2 as well as over thirty additional cell surface and cytoskeletal regulatory molecules that could facilitate axon guidance (e.g., Ncam1, Plxna4, Plxnb3, Cdh3; see Table S1, Figures 1C, and S5B)^40–43^. Because the expression of DV score genes is independent of OR activity, often precedes OR choice, and is bidirectionally altered upon manipulations of RA signaling, DV GEPs represent a potentially parsimonious mechanism to couple epithelial space to both the process of OR choice and axonal guidance.
To address the possibility that DV genes might link epithelial spatial positions to axonal targets in the olfactory bulb, we took advantage of a recent spatial transcriptomics dataset that identified the glomeruli for many OSN subtypes^79^. Coloring identified glomeruli based on the DV scores of their associated receptors revealed that DV scores varied smoothly across the olfactory bulb and that neighboring glomeruli had similar DV scores (Figures 7A–C). As a result, the 3D position of each glomerulus predicted the DV score of the innervating OSN subtype with high accuracy (ρ = 0.95, Figure 7D), and conversely the DV score of each OSN subtype predicted the 3D location of its glomerulus (median error = ~400 μm, 3–4 glomerular widths, Figure 7E). Thus variation in DV gene expression (and by extension, epithelial location) across OSN subtypes can at least in part explain the spatial organization of their axonal projections to the bulb.
Notably, the AP GEPs also include genes (like Nrp1 and Plxna1) known to mediate axon targeting along the bulb anteroposterior (AP) axis^46,47,53^, as well as additional axon guidance and cytoskeletal genes (e.g., Plxna3, Sema6c, Pcdh11x, Cdh15, see Table S3). Like the DV GEPs, the use of the AP GEPs also smoothly varied across olfactory bulb glomeruli, although along an orthogonal spatial axis to that used by the DV GEPs and in a manner that was somewhat less predictive of glomerular position (Figures 7A–E). Glomerular 3D positions for each OSN subtype were best predicted by combining DV and AP GEP usage; the error of these predictions (~300 μm) was similar to that observed using models trained with all ~1,300 OSN variable genes (Figure 7E), indicating these two GEPs together effectively summarize the relevant transcriptional axes that map OSN subtypes onto bulb position. Thus, the two cardinal axes (DV and AP) that organize the olfactory bulb are topographically and transcriptionally represented in the epithelium itself.
The peripheral olfactory system has long been thought to balance determinism and randomness by limiting the expression of individual ORs to one of a handful of dorsoventral zones, while randomly distributing ORs within each zone^13,14^. Here we show instead that the olfactory epithelium harbors an organized map in which each of the ~1100 receptors predictably occupies a unique and distinguishable distribution along the dorsoventral axis of the epithelium. This pattern reflects a developmental program in which RA signaling transforms spatial position into transcriptional gradients across ~250 co-regulated genes. This coherent transcriptional axis — summarized as the DV score — allows us to demonstrate that precursors and OSNs harbor spatially specific transcriptional identities that restrict OR choice and coordinate the construction of spatial maps in the epithelium and bulb (Figure 7F).
These findings reveal that development can reliably organize ~1100 different sensory channels in space. Dorsoventral position influences many key steps that precede OR choice during OSN differentiation, including OR co-expression, heterochromatin deposition, and inter-chromosomal genomic interactions (Figure 5D). However, while the mean relationship between spatial position and OR choice is essentially invariant at the population level, precursors do not always pick their most-preferred OR, resulting in each OR occupying a tight distribution of positions. We speculate that variability in the timing of the steps that translate space into OR choice might allow noise to accumulate during differentiation; consistent with this possibility, DV score distributions are wider for ventral OSN subtypes (Figure S1P).
Our work is consistent with transplantation experiments that have suggested that spatial determinants influence OR choice^80^. RA-related genes like Aldh1a2 have also previously been proposed to signal information about physical space, but testing the role of RA on OSN differentiation has been challenging given the pleiotropic roles of RA on neuronal development^40,62,81–83^. Here we implicate RA in organizing the olfactory epithelial receptor map via pharmacological manipulations, which demonstrate that RA controls the coordinated expression of DV genes and spatially restricts OR choice. The ability of a continuous RA gradient to generate discrete variation in OR expression along the dorsoventral axis of the olfactory epithelium is reminiscent of its role in development where, for example, a continuous RA gradient ultimately supports the ontogeny of different somites along the rostrocaudal axis^84,85^.
What is the meaning or purpose of the receptor map in the nose? While at high concentrations monomolecular odors activate spatially distributed sets of ORs, at lower concentrations the sparser set of activated ORs have, on average, more similar associated DV scores. These activated ORs are nevertheless separated by ~100–200 intervening ORs, and as such the observed spatial clustering is modest at best. We hypothesize that the “weak” chemotopy observed here (and by others^17,86^) is driven by recently duplicated ORs that are more likely to respond to similar odorants; as different OR genes accumulate genetic variation, their spatial positions likely diverge due to variation in upstream promoter regions, which are less conserved than coding sequences.
If odor tuning properties are, to a first approximation, spatially distributed, why is receptor position organized in the nose? One obvious possibility is robustness, as spatially distributing tuning properties renders the epithelium resilient to local insults. In addition, recent results suggest that odor experience can influence the number of OSNs that respond to an experienced odor^87–90^. Our observation that OSN precursors of a given DV score and therefore spatial identity sit immediately underneath the OSNs to which they give rise might also allow for long-term adaptation in OSN frequency. We speculate that odor responses in OSNs generate signals that increase proliferation and differentiation of underlying stem cells by either diffusing locally or by stimulating adjacent sustentacular cells. However, because single OSN precursors give rise to more than one OSN subtype, a main prediction from our results is that the odor experience will recruit not just OSNs whose receptors respond to the experienced odors but also OSNs with similar DV scores; future work will be required to test this hypothesis.
Spatial variation may also help OSNs to simultaneously know both which OR they are to express and where in the brain they should send their axons. In this model, axonal targets are largely determined by the physical location of each OSN in the epithelium rather than by the OR chosen for expression^12,48^. Nevertheless, ORs likely still play causal roles in glomerular targeting and refinement^48,49^, as the mapping between spatial positions, DV scores, and OR choice is not fully deterministic. Furthermore, AP genes are receptor dependent, include key causal genes for guiding axons along the anteroposterior bulb axis, improve predictions of glomerular location, and likely play important roles in these processes^47,54,91^. Thus, while OSN identities reflect continuous variation in DV gene expression, this continuity is segmented by the expression of single ORs, which ultimately enables accurate targeting to discrete glomeruli.
First, much remains to be discovered about how RA patterns the MOE and translates space into coherent patterns of DV gene expression. Additional diffusible factors may also convey spatial information to epithelial cells^92,93^. This is likely true in the dorsal-most region of the epithelium, which appears less sensitive to RA manipulations than the ventral epithelium and harbors distinct precursor types that generate OSNs expressing Class I and TAAR-type receptors^94–96^; future work is required to understand how space relates to choice for these receptors, and for those few ORs that are organized into epithelial patches^97–99^. Second, although we find that upstream genomic sequences largely determine the epithelial positioning of each OR (extending past work showing OR promoters can confer spatially-restricted transgene expression^74,100–102^), it remains unclear how a transcriptional code might license individual ORs for choice. Third, while DV genes like Robo2 and Nrp2 are required for appropriate glomerular targeting, dissecting the causal role of additional DV genes in OSN differentiation, fate determination, and axon guidance may be challenging due to either redundancy or pleiotropy. Finally, how the OR-dependent AP GEPs influence OSN positioning along the epithelial apical-basal axis remains a mystery.
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Sandeep R. Datta (srdatta@hms.harvard.edu).
This study did not generate new unique reagents.
C57BL/6J mice were obtained from Jackson Laboratory (Stock No. 000664). Mice expressing OMP-IRES-tTA and tetO-P2-IRES-GFP (OMP-P2 mice) were obtained from the Lomvardas lab and maintained in the Datta Lab^55,104^. Cnga2 knockout animals were obtained from the Dulac lab and maintained by breeding heterozygous females with C57BL/6J males^56^. HP-1 swap mice (mice with loxP-HP1β-STOP-loxP-HP1α alleles at the HP1β locus crossed to Foxg1-Cre [Jackson Laboratory Stock No. 006084]^105^ such that knockout of HP1β is rescued via the expression of HP1α) were generated by Martín Escamilla-Del-Arenal using standard homologous recombination, as recently described^71^. Wild-derived CAST/EiJ (Stock No. 000928), MOLF/EiJ (Stock No. 000550), and PWK/PhJ (Stock No. 003715) mice were obtained from Jackson Laboratory, and CAST/C57 F1 hybrids were generated in lab by crossing CAST males with C57 females. As described, we found little effect of sex on the DV score (Figure S1N), and thus mice of either sex between 6–16 weeks old were used for experiments. Mouse husbandry and experiments were performed following institutional and federal guidelines and were approved by Harvard Medical School’s Institutional Animal Care and Use Committee.
Data from wild-type mice housed in home-cage environments, or from mice that either underwent transient naris occlusion for a week or were housed in novel olfactory environments, were all previously generated^54^. Act-seq was performed as previously described, by exposing mice to monomolecular odorants on filter paper for two hours prior to OSN dissociation and scRNA-seq^54^. Odorants were diluted in dipropylene glycol (DPG) based on their vapor pressures to give a nominal concentration of 500 ppm at their highest concentrations.
The main olfactory epithelium (MOE) was ablated and induced to regenerate via a single intraperitoneal injection of methimazole (50 mg/kg) in adult mice (>6 weeks old), as described elsewhere^118^. Intranasal lentiviral injections were performed 36–48 hours after methimazole injection, at timepoints in which the OSNs have sloughed off and the epithelium is a single monolayer of activated horizontal basal stem cells (HBCs). Mice were anaesthetized with Ketamine/Xylazine (100/10 mg/kg), and intranasal lentiviral injections were performed via a Hamilton syringe connected to thin PE-10 tubing; 15–25 μL of virus (at a titer of 1–10 × 10^8^ TU/mL) was injected over the course of 2–3 minutes. scRNA-seq was performed 4–6 weeks after methimazole injection, at timepoints under which much of the epithelium has regenerated and the infected progenitors have differentiated to give rise to clones of newly-born cells, many of which are OSNs.
Manipulations of retinoic acid signaling were performed in adult mice that received a single dose of methimazole (50 mg/kg). Starting 72 hours after methimazole injection, mice received either the ALDH1A2 inhibitor WIN 18,446 (5 mg/kg) all trans-retinoic acid (atRA, 5 mg/kg), or the RA receptor inverse agonist BMS493 (10 mg/kg). WIN 18,446 was dissolved in DMSO, atRA was dissolved in 5 % DMSO (dissolved in corn oil), BMS493 was dissolved in 5 % DMSO (dissolved in corn oil), and drugs were delivered subcutaneously. Littermates within each cohort received either drug or vehicle daily for 3–4 weeks after methimazole-induced regeneration, and scRNA-seq was performed after 4–5 weeks; animals receiving each drug or vehicle control were dissociated and sequenced or processed for in situs in parallel. The effects of the RA inhibitors WIN 18,446 and BMS493 were verified via observation of the testes^119^, which were typically 30–50% the weight of littermate controls.
Isolation of single cells from the MOE and fluorescence-activated cell sorting (FACS) was performed as described^54^. In brief, mice were euthanized via carbon dioxide inhalation, and the epithelial tissue of the MOE was dissected. Individual mice were used for each replicate. Cells were dissociated with Papain and DNase-I (Worthington), incubated at 37°C for 60 minutes, triturated and then washed and resuspended in Hibernate-A medium supplemented with fetal bovine serum; these solutions also contained transcriptional and translation inhibitors (5 mg/mL of Actinomycin D, 10 mg/mL of Anisomycin and 10 mM of Triptolide, all obtained from Sigma). FACS was performed to enrich for OSNs expressing fluorescent reporters in the OMP-P2 and lentiviral experiments. In all other experiments, FACS was performed to enrich for live singlets.
scRNA-seq was performed using the 10x genomics platform using the default protocols (CG000315 Rev E) for the Chromium Next GEM Single Cell 3ʹ Kit v3.1 and the default protocols (CG000731 Rev B) for the Chromium GEM-X Single Cell 3’ Reagent Kits v4. For each replicate, cells were loaded at a concentration predicted to maximize yield (10,000 cells for v3.1 kits and 20,000 for v4 kits), or to obtain the maximum number of cells given the measured concentrations when the total number of cells was less than 10,000. For the OMP-P2 experiments, GFP-positive and GFP-negative cells were mixed and loaded together in a single lane of the Next GEM chip. In the lentiviral experiments, Venus-positive cells and negative cells from each mouse were not mixed and were loaded into separate lanes. For experiments with genetic or experimental controls (like the CNGA2, RA, and HP1 swap experiments), control and experiment animals were processed together, loaded onto the same chips, and were sequenced together.
The resulting libraries were quantified via qPCR and via the Qubit and were sequenced with paired-end sequencing on the Illumina NovaSeq or NextSeq platforms (with Read1 = 28 cycle, Index (i7) = 10 cycle, Index (i5) = 10 cycle, Read2 = remaining number of cycles). Demultiplexed fastq files were processed in a manner similar to that previously described^54^, using a custom Nextflow pipeline that ran the Cell Ranger (version 6.1.2 for v3.1 and v9.0.1 with --include-introns= false for v4) count pipeline, as well as additional post-processing steps to identify and remove multi-mapped reads and regenerate the resulting cell by gene matrix. The Ensembl GRCm39 mouse genomic index (version 105) was used with a custom GTF file containing extended 3’ UTR annotations for some ORs and TAARs, to facilitate their identification in the 3’-biased scRNA-seq data. Cell doublets and low-quality cells were excluded, and OSNs were identified via iterative subclustering of the latent spaces of scVI models, as previously described^54,120^, based on their expression of known marker genes and ORs.
Clonal lineage tracing was performed using a lentiviral library with ~10M DNA barcodes at the 3’ UTR side of transcripts for Venus fluorescent protein and upstream of the polyA sequence (Cellecta CloneTracker XP). To amplify the lentiviral barcodes from the RNA transcripts of in each cell, a primer (5’-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTCCGACCACCGAACGCAACGCACGCA-3’) whose 5’ end overlapped with the Truseq Read 2 (underlined) and whose 3’ end was complementary to a constant region upstream of the barcode was spiked in (final concentration 0.5 μM) during the initial RT and cDNA amplification steps; this approach increased the yield of the lentiviral barcodes compared to targeted amplification from the already-amplified cDNA libraries. Amplified barcode fragments were purified from the cDNA library and barcode sequencing libraries were constructed by performing PCR on these fragments with 10x Library Index PCR primers (which bind to Truseq); the resulting barcode libraries were mixed and sequenced together with the final 10x libraries and were subsequently demultiplexed computationally.
MERFISH was performed using the Vizgen MERSCOPE platform. A custom 300-plex codebook, which contained nearly 200 ORs and 100 non-OR genes consisting of cell type marker genes, DV-related genes, and other highly-variable OSN genes, was designed. Bulk sequencing and scRNA-seq was used to identify suitable ORs and genes, and the total abundance of the gene panel was ~7600 FPKM to avoid optical crowding. ORs were chosen based on their DV scores and frequency of expression to obtain a set of ORs predicted to uniformly span the MOE at relatively-higher frequencies. Due to their high abundance, additional probes for OSNs (Omp, Calm2), and sustentacular cells (Cbr2, Cyp2g1) were imaged separately and sequentially (along with the staining for DAPI, polyT, and cell boundaries).
Fresh frozen 10 μm coronal samples from the MOE from young (4–5 week) mice were processed according to the manufacturer’s user guides (Vizgen 91600002 Rev E and 91600112 Rev C). In brief, tissue sectioning, fixation, permeabilization, and autofluorescence quenching were performed used the guidelines for fresh frozen tissue. However, the RNA anchoring steps from the formalin-fixed, paraffin-embedded (FFPE) tissue sample preparation guide were added to help the thin MOE tissue adhere to the MERSCOPE slides. Cell boundary staining, gel embedding, clearing, and probe hybridization steps were performed without modification from the fresh-frozen sample preparation guide. Samples were imaged on the MERSCOPE instrument using the 300-plex imaging kit, and transcripts were decoded using the MERLIN pipeline provided by Vizgen.
Stereo-seq was performed using the STOmics Stereo-seq Transcriptomics T v1.3 kits following the manufacturer’s user guide (STUM-TT001). In brief, 10 μm coronal sections from fresh frozen MOE samples (from young 4–5 week mice) were sectioned onto 1 mm^2^ Stereo-seq Chip T Slides, dried at 37°C for 5 minutes and methanol-fixed for 30 minutes at −20°C. The cells were stained with Qubit ssDNA reagent, and the slide was then rapidly imaged on an Olympus VS200 slide scanner with a 10x objective. Following imaging, slides were permeabilized at 37°C for 12–15 minutes, and the captured RNA was reverse transcribed for 2 hours at 45°C. The resulting cDNA was released at 55°C for 10 minutes and then amplified via PCR and size-selected with SPRIselect beads. 100ng cDNA was used for each sample and sequencing libraries were prepared according to the manufacturer’s user guide (STUM-LP002), using an equimolar mix of 4 index barcodes per sample. Two samples were pooled per sequencing run, and the resulting libraries were sequenced on DNBSEQ-T7 flowcells in the PE75 (50+100) configuration; demultiplexed fastq files were used for downstream analyses.
The Hi-C variant Micro-C, which uses MNase instead of restriction enzymes to digest cross-linked DNA^121^, was used to survey the chromatin conformation of mature olfactory sensory neurons. Micro-C was performed using the Dovetail Genomics Micro-C Kit (PN 20101E) following the instructions in its user guide (v2.1 for cells). In brief, ~1 million EYFP^+^ OSNs were dissociated and FACS-purified from OMP-IRES-Cre;Ai3(RCL-EYFP) animals for each sample, following the approach described above in the scRNA-seq experiments. Sorted cells were pelleted, frozen, and then thawed and cross-linked with DSG and 37% formaldehyde. Next, cells were digested with MNase for 15 minutes at 22°C. The digested lysate was proximity ligated, and DNA was purified. To increase library complexity, three technical replicates were prepared for each sample, using 150 ng of purified DNA per replicate. Fragments with ligation junctions were enriched with streptavidin beads, and index PCR and size selection were performed. Libraries were pooled and sequenced with NovaSeq X Plus 25B flowcells, targeting ~1B reads per sample. The resulting fastq files from each technical replicate were processed in parallel and were pooled for downstream analyses.
CUT&RUN was performed to measure the levels of H3K9me3 in C57/CAST F1 animals, using the reagents and instructions from a commercial kit (EpiCypher CUTANA CUT&RUN Kit, Version 4). In brief, the MOE was dissociated as in the scRNA-seq experiments, and 500k unsorted cells were used for each reaction, and two technical replicates were used for the cells from each animal. Cells were incubated overnight with ConA beads and 0.5 μg of a Rabbit polyclonal H3K9me3 antibody (Abcam ab8898), and on the following day permeabilized cells were incubated with pAG-MNase to target the digestion of DNA, which was purified with SPRIselect beads. E. coli spike-in DNA was used to normalize the yields across samples. Sequencing libraries were prepared from 5 ng CUT&RUN DNA, or the entire yield for replicates where the yield was <5 ng. Some experiments also included a H3K4me3 positive control (EpiCypher 13–0041), which yielded the expected enrichment for transcriptional start sites, IgG negative control (EpiCypher 13–0042) antibodies, which showed minimal enrichment, but these were omitted for subsequent experiments given the low background in the H3K9me3 CUT&RUN data. Sequencing libraries were generated using the NEBNext DNA library kit, following the instructions and formulations in the EpiCypher CUT&RUN library prep kit (14–1002). CUT&RUN libraries were pooled and sequenced on the Element AVITI machines with the Cloudbreak Freestyle sequencing kit (2×75 cycles).
In situ hybridization for genes in the MOE was performed via the Hybridization Chain Reaction (HCR) method^122^, using 16–20 μm cryosections from either fresh or fixed frozen MOE tissue, from 4 to 8-week-old animals. Paraformaldehyde-fixed tissue was post-fixed overnight, decalcified in 0.45M EDTA overnight, and equilibrated in 30% sucrose prior to freezing. In situs were performed according to manufacturer instructions for HCR v3.0 (Molecular Instruments, Rev. 4), except the proteinase K solution was changed to 1 μg/mL or omitted entirely and probe concentration was increased to 1–2 pmol per 100 μL.
Split DNA oligo probes for ORs (~15–20 probes per gene) were manually designed using regions in each OR’s cDNA sequence that were homologous between C57BL/6J and CAST/EiJ animals but lacked homology to other ORs. Probes for the dorsal marker Nqo1 and reference ORs whose positions remained constant across strains were also designed in a similar manner and used as reference landmarks to identify the positions of ORs of interest. For in situs in animals treated with the RAR inhibitor, a pooled OR probe was designed to target the P2 OR (Olfr17) along with eight additional ORs with similar DV scores (Olfr76, Olfr393, Olfr713, Olfr873, Olfr1302, Olfr1341, Olfr1370, and Olfr1490) whose frequency also increased in the RA-inhibited animals. The pooled OR probe was ordered as an oPool oligo pool from IDT (at the 50 pM scale) and was resuspended to a concentration of 1 μM per oligo and used at a concentration of 1 pmol per 100 μL. The resulting probe mix labeled OSNs in the same region of the epithelium as P2 alone, consistent with their similar DV scores. The sequences of all oligo probes can be found in Table S4.
To further validate the in situs, antibody staining was also performed for Nqo1 in a subset of animals. After the HCR protocol, sections were post-fixed with 4% PFA for 15 minutes, permeabilized and blocked for 30 minutes in PBS with 0.3% Triton X-100 (PBST) and 5% Normal Donkey Serum, and then incubated overnight with Goat anti-Nqo1 (1:100) in blocking buffer. The following day, sections were washed with 0.1% PBST and then incubated for 2 hours at room temperature with Donkey anti-Goat Alexa Fluor 488 secondary antibodies.
Sections were counterstained with DAPI and confocal imaging was performed using a Nikon Ti inverted spinning disk microscope equipped with a Yokogawa CSU-W1 Spinning disk (50 μm pinhole), a Plan Apo 20x/0.8 NA air objective and a Plan Fluor 40x/1.3 NA oil objective, and 405, 561, and 640 nm laser lines. Tiled images of the MOE were captured using an Andor Zyla 4.2 Plus sCMOS monochrome camera and Nikon Elements Acquisition Software (AR 5.02). For sections probed with individual ORs, the OSNs expressing each OR were manually counted and the fraction of OSNs within the dorsal zone was assessed by comparing to positions of each OSN relative to the Nqo1^+^ region of dorsal OSNs. Note, although imaging and quantification was performed across the entire section, images depict OSNs close to the edge of the Nqo1+ region due to the difficulty in visualizing the expression of single ORs across the entire MOE. Due to their increased density and to avoid experimenter biases, in the RAR inhibitor experiments ORs were automatically identified and counted using the cell detection algorithms in QuPath^117^. To compare across samples, endoturbinate III (using the nomenclature of Barrios et al.^123^) was chosen because it is readily identifiable, and the number of cells and their distance along an annotated path that spanned the turbinate from its tip were quantified.
Nearly 5 million cells from over 360 replicates were uniformly processed using a custom Nextflow pipeline to run Cell Ranger and custom postprocessing steps for each replicate. These replicates included the 150 replicates (and 780 thousand mature OSNs) described in a previously published dataset^54^, as well as over 200 additional replicates (containing another 1.4 million OSNs). The additional data consisted of the mice used as part of the experimental conditions described herein (e.g. the lentiviral experiments), as well as additional unpublished datasets that were generated for other work, including from mice that had been exposed to monomolecular odorants, as well as from smaller numbers of mice housed in different environmental conditions, and control and knockouts for various genetic manipulations. Mice of both sexes were used for experiments. Note, although some of these lines had not been backcrossed to C57BL/6J backgrounds, few differences in DV scores or gene expression were observed across replicates, even though the entire dataset consisted of 360 replicates from many separate experiments.
A combination of unsupervised and semi-supervised methods was used to map the cells from all replicates into a common latent space to perform cell type identification. The gene set used for this embedding were the top 3000 genes (excluding ORs, mitochondrial and ribosomal genes) identified using the Scanpy highly variable gene function, with the “seurat_v3” flavor. First, unsupervised learning was used to fit a scVI model (n_hidden: 128, n_latent: 30, n_layers: 2, gene_likelihood: “nb”, dropout_rate: 0.2, with a batch key to identify each replicate) on the raw counts of these genes in a smaller reference dataset of ~250 thousand cells from wild-type animals housed in standard home-cage conditions. The resulting scVI latent embeddings from this dataset were clustered via the Leiden algorithm and cluster labels were manually annotated based on known cell type markers. The trained scVI model and the manually-annotated cluster labels from this subset were then used to train a semi-supervised model using scANVI. The trained scANVI model was then applied to query the entire dataset, resulting in a common embedding of all cells into a 30-dimensional latent space. Cell types in this common latent space were identified by transferring the labels from the manually-annotated reference subset, using previous weighted nearest neighbor approaches for label transfer of scANVI models^124,125^. For cells from each replicate, their resulting neighbors (in the 30-dimensional latent scANVI space) from the annotated subset were found using PyNNDescent, distances were converted into affinities using a gaussian kernel, and cell type probabilities and matching uncertainty were calculated for each cell.
Non-neuronal cells, as well as low-quality doublets and dying cells were excluded, and mature OSNs were extracted from the cell type labels of the integrated dataset. Due to the dissociation and FACS approaches used, roughly 50% of cells from most replicates were OSNs. OSNs were filtered to retain only cells expressing a single OR, and the resulting dataset of 2.3 million OSNs was used for downstream analyses. Except where noted, only cells expressing single ORs were used for analyses, and OR expression was defined with a threshold of at least 3 transcripts (unique molecular identifiers (UMIs)) that mapped to that OR.
Zonal indices for each OR were obtained from prior work^15,16^, and were discretized into five zones to match prior conventions, with zone 1 being the most-dorsal zone and zone 5 the most-ventral. The mean for each subtype for each variable gene was calculated and genes whose expression was on average (log2 fold-change > 1) higher among the subtypes in its maximum versus second-highest zone were considered as zonally-restricted genes, and the number of genes for each zone were computed. Note, no genes were enriched in intermediate zone ORs (e.g. zone 2–4), consistent with gradients of gene expression that monotonically increase or decrease with the dorsoventral score. To account for the possibility that genes might enriched in alternate regions (like 2.5–3.5, or 2.0–2.5, or 2.1–2.6), bin widths of 0.5 and 1.0 were tested, and were slid across the entire zonal range (at 0.1 intervals) to quantify the number of genes with log2 fold-changes > 1 for each bin location. Such analyses also failed to identify any zonal genes for intermediate zones. Since no genes were identified, no correction for multiple comparisons was performed.
Consensus non-negative Matrix Factorization (cNMF) was used to decompose the transcriptomes of each OSN into a smaller set of interpretable factors^54,109^. The same NMF factors (referred to as gene expression programs (GEPs)) that were previously identified were used, and the gene loadings for these GEPs were applied to all datasets^54^. Importantly, cNMF was performed on a set of 1,300 highly variable genes (HVGs) that did not include any OR genes and thus the resulting gene loadings and GEP usages do not depend on the expressed ORs. These GEPs included two associated with ligand-dependent activity, (GEPHigh~ and GEPLow), two containing dorsal and ventral genes (GEPDorsal and GEPVentral), and two relating to anteroposterior positions (GEPAnterior and GEPPosterior). Therefore, even though we did not explicitly attempt to identify GEPs that might reflect spatial position, the fact that this approach yielded spatially-varying GEPs indicated that epithelial positions might be reflected in OSN transcriptomes. These GEPs were used in an orthogonal manner and their differences in usage was summarized into single the environmental state (ES) score = GEPHigh - GEPLow, the dorsoventral DV score = GEPDorsal - GEPVentral, and the anteroposterior (AP) score = GEPAnterior - GEPPosterior) to capture the continuous variation across these three independent axes (ES, DV, and AP). On the other hand, cNMF also identified GEPs distinguishes OSNs that express CD36, or are located in the “patch” region of the MOE, or that have high levels of genes related to the unfolded protein response. The usage of these GEPs was less graded, suggesting that the orthogonal nature of GEPDorsal - GEPVentral is not a necessary consequence of our cNMF approach but rather reflects the continuous transcriptional variation with respect to spatial position (which we validate with subsequence experiments).
The 248 genes (out of the set of highly variable genes) whose expression across OSNs was correlated with the DV scores (R^2^ > 0.5 for linear regression to predict the DV score from each gene and spearman’s ρ > 0.5 at the OSN subtype level) were considered as DV genes and were used for downstream analyses. DV genes were categorized into 10 different functional groups (transcription factors, axon guidance and cell adhesion, cellular metabolism, cell signaling, cell growth, cytoskeletal, GTP-activity, RNA-related, ion channels and transporters, and synaptic genes) via a combination of GO pathway analysis, the presence of protein features like DNA binding domains, or via manual annotation. The full list of DV genes and their associated categories (for the 162 belonging to the above categories) can be found in Table S1.
As described above, GEP usages and DV scores were calculated at the cell level using the gene loadings for each GEP. A small fraction (3%) of cells had DV scores that were zero, which can occur for cells that had little expression of either dorsal- or ventral-related genes and because NMF identifies usages for all GEPs at once. For a subset of analyses that required comparisons at the cell level, the DV scores of these cells were imputed using a multilayer perceptron (MLP) trained to predict the DV score of mature OSNs using the expression of the 248 DV genes in each cell. The MLP (with two hidden layers of size 50 and 10) was trained in PyTorch with the Adam optimizer (learning rate 0.01) and gave predictions that were highly correlated (ρ > 0.99) at both the cell and OR level for training, testing, and validation data.
Mature OSNs expressing a given OR (i.e. those from the same OSN subtype) had similar DV scores. Therefore, the DV score associated with each OSNs expressing each OR (also referred to as an OR DV score) was calculated by taking the mean DV score across all the OSNs expressing a given OR in the entire integrated dataset. Except where noted, for all downstream analyses the DV score of the expressed OR in each cell was considered as the OR DV scores from the entire dataset and thus could be compared to each cell DV’s score, especially in experiments that decoupled a given cell’s DV score from that of its expressed OR. In some experiments, DV scores were percentile normalized at the cell level via sklearn’s QuantileTransformer. This transformer was fit on data from all cells to maintain the distribution of DV scores for cells from each subtype.
Only ORs with functional protein activity (as assessed by their singular expression in mature OSNs). ORs with mean DV scores above 60 were considered as dorsal; this threshold was set such that ORs that had previously been classified as being part of the most-dorsal zone (zone 1) were considered as dorsal with respect to their DV scores. ORs were clustered based on their amino acid sequences and class I ORs were identified as the large cluster of ORs in the resulting phylogenetic tree that share homology with fish ORs and were located in a single genomic cluster on Chromosome 7. Phylogenetic distances were defined as the cophenetic distance of the resulting phylogenetic tree. Most rodent odorant receptors are class II ORs, which are located in multiple genomic clusters of genes that are scattered across the genome. ORs were considered part of the same genomic cluster if they were located on contiguous regions of the chromosome (separated by less than 3 Mb from the start of the previous OR gene). Neighboring ORs (i.e. those at neighboring genomic locations) were more likely to be phylogenetically similar. A small fraction of cells also express trace amine-associated receptors (TAARs). Except where specified, analyses used OSNs expressing either class I, class II, or TAARs though most results reflect properties of class II ORs, which make up the majority (~87%) of all ORs. The set of 1007 OSN subtypes with at least 150 cells (median 1,500 OSNs per subtype) was used for most analyses.
The DV genes all had strong positive/negative correlations with the DV score and with each other, suggesting they act as a coherent program. To verify these results, and to ensure that the DV scores weren’t being driven by small numbers of genes, we took a series of validation approaches. First, every single DV gene was individually removed and GEP usages were recalculated (by reweighting the gene loadings for each GEP to sum to 1 and applying the GEPs via non_negative_factorization with update_H=False). Second, to avoid any confounds in relating the DV score to the expression of known marker genes, the marker genes Nqo1, Nfia, Nfib, Nfix, Ncam2, Nrp2, Robo2, and Acsm4 were all removed and GEP usages were similarly refit. In both cases, the DV scores for each OSN subtype were perfectly correlated (ρ = 0.999) with those calculated using all genes, suggesting that individual genes have redundant roles in this larger program. Alternatively, the sufficiency of small numbers of genes at predicting the DV score was tested using regression models that predicted the DV score for each OSN subtype via the mean gene expression for that subtype. Elastic-net-regularized linear regression models were used (with l1_ratio=0.95) with cross-validation (training on 80% of OSN subtypes and testing on the held-out 20%). Random sets of 1–248 either random HVGs or DV genes were used across 1000 restarts, and the distribution of prediction accuracies (r2_score) for each type of gene set was summarized. These models were also compared to those that attempted to identify the best set of n DV genes for values of n from 1–248 (i.e. using the best DV gene or all DV genes). To do so, the best set of the DV genes in the training data on each restart for each value of n was identified using recursive feature elimination (RFE) with the same Elastic-net model.
Regression analysis was also used to evaluate how well the expression of each DV gene in each OSN subtype could be predicted by either the DV GEPs or a zonal model. Model comparison was performed between six regression models. First, gene expression across OSN subtypes was fit using the default zone indices for the canonical 5-zone model (with zone 1, 1–2, 2–3, 3–4, and 4–5 as categorical features). Five alternative models were 1) the zones, but offset by half (≤1.5, 1.5–2.5, 2.5–3.5, 3.5–4.5, >4.5); 2) a piecewise-linear model (PiecewiseLinFit in the pwlf library) fit using the DV score with 4 segments; 3) a LinearRegression model that used the usage of the two DV GEPs for each OSN subtype; 4) a KNN regression model (with uniform weights and 20 neighbors) fit on DV GEP usage in the training data; or 5) a LinearRegression model in which the two DV GEPs were expanded into a larger quadratic basis using the PolynomialFeatures(degree=2, include_bias=False) transformation. All models were fit via five-fold cross-validation, and the change in prediction accuracy (Δ r2_score) for each model for each gene was compared to that from the default zonal model.
Correlations of the DV scores across replicates were computed at the OSN subtype level, by comparing the mean for each subtype in a given replicate to that from all other data. To evaluate the consistency across data subsets, 150 cells were sampled for each OSN subtype, the mean for each OSN subtype was calculated separately for the first versus second half of cells for each subtype, the spearman’s correlation between data halves was evaluated, and this process was repeated 1,000 times. The Spearman’s correlation reflects changes in the rank order of OSN subtypes and the absolute change in ranks was also evaluated on each restart, as well as for data in which OR labels were either permuted across cells whose ORs shared the same zonal index (discretized into 5 zones as above) or permuted across all cells. The mean for each OSN subtype was also evaluated separately for data from male and female animals. In an alternative approach, the correlation in DV scores was computed for all pairs (363C2) of mice using, for each pair, the set of ORs with at least 4 cells in both samples. Even with only 4 cells per sample, the median spearman’s correlation was still 0.975.
The uniqueness of the DV score for each OSN subtype was determined using linear regression. Across 100 restarts, cross-validation was performed using five stratified folds for training and testing. The mean for each OSN subtype was computed using the training data and in the testing data the observed DV score for each OSN was compared to that of the either the mean of its respective OSN subtype in the training data or the mean of the n^th^ nearest OSN subtypes (as determined based on the means of each OSN subtype in the training data); the difference in median absolute errors was computed for models with n ranging from 0–200.
For computing deviations from the mean and area under the receiver operator curve (AUROC), equal numbers of OSNs (150) were sampled for each OSN subtype on each of 1,000 restarts, and AUROC values were calculated empirically for all pairs of OSN subtypes (1007C2 = 506,521 pairs). For comparisons to shuffled data, OSN labels were shuffled across the subsampled cells on each restart. The difference in DV scores for each cell from the respective mean of its corresponding OSN subtype was also summarized via the median absolute deviation (MAD). AUROC values were also evaluated for pairs of ORs found within the same zone, using either the default 5-zone model described above (and plotted as a function of their difference in DV scores), or the “Zolf” zone annotations described in a recent paper.
Data from different epithelial sections were aligned to a common template using a two-step approach. First, since the sections were from similar regions of the epithelium, they were all roughly aligned to a reference template with a simple rigid transformation. Second, these sections were then aligned STalign^126^, an approach that relies on diffeomorphic mapping (to help account for slight deformations in the thin epithelial layer). Overlaying the raw images and marker gene expression following alignment confirmed this approach successfully aligned the various tissue layers. MERFISH gene expression was also visualized using Baysor^127^, in its segmentation-free configuration. In brief, local neighborhoods were evaluated by finding the nearest transcripts (across all genes) for every gene, reducing the dimensionality of this data via principal component analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP)^128^, and converting the resulting low dimensional embedding into colors for visualization purposes. Because OSNs are small and densely packed and difficult to segment via DAPI or cell boundary staining alone, cell segmentation from the Vizgen MERlin pipeline (via cellpose) was further refined using Proseg^129^, a probabilistic approach that helps to better match transcripts to cells. The resulting log-normalized cell by gene matrix from Proseg was processed in Scanpy. PCA was performed and the top 30 PCs were kept. Leiden clustering was performed on a random subset of 100,000 cells (on a nearest neighbor graph with n=30 neighbors and leiden resolution=1.2) and the labels were transferred to all cells using the weighted nearest neighbor approach described above for the scRNA-seq data. Clusters were annotated and merged based on their expression of known marker genes (e.g for OSNs, immature precursors, and non-neuronal genes); the spatial distribution of these genes and their respective clusters were also evaluated, and matched prior knowledge. As indicated, downstream analyses that used a specific cell type (e.g. INPs or mesenchymal lamina propria cells) used only the cells and or transcripts found within those cells of a given cluster. Similar spatial segregation was also observed at the transcript level when the Baysor local neighbors for each transcript were clustered using KMeans clustering (with k=3–10 clusters).
OSNs expressing single ORs were identified in an iterative process, using the spatial locations of OR transcripts. First, only transcripts in the OSN layer (identified via the clustering described above) were considered. Next, to avoid false positives from background expression (or errors in spot decoding) OR transcripts were filtered for each gene to only keep those with at least 3 other transcripts within a 10 μm radius (the average size of an OSN), and the resulting transcripts were then clustered into individual cells using the DBSCAN algorithm (eps=5 and min_samples=3). The set of OSNs was further refined by removing 0.1% of cells in which the median pairwise distance between OR transcripts within the cell was greater than 10 μm. OSNs were identified separately for each sample, using the STaligned coordinates. Four OR genes out of ~200 showed nonspecific binding across the MOE (and were also decoded in damaged areas and areas with high autofluorescence) and were not considered for downstream analyses. Only ORs detected in at least 20 different OSNs were used for analyses.
The DVISH score was computed directly using the MERFISH-measured DV gene expression for each cell. To infer DVISH scores in OSNs, the total expression of each gene for each cell (all transcripts within a 10 μm radius of the cell centroid) were total-count normalized. Next, rather than applying cNMF directly (because cNMF was fit on scRNA-seq data and the MERFISH panel only contained a subset of all HVGs), PCA was performed on the total-count normalized expression of the DV genes that were part of the MERFISH panel and the first PC was taken to be the DVISH scores for each cell. DVISH scores smoothly varied across the MOE and matched the underlying gradients of individual DV genes. DVISH scores were calculated in a similar manner for INPs, using the set of INP cells that were located basally to OSNs that did not have detectable OR transcripts and the counts of the nearest 100 transcripts to each cell. The DV score for each cell was taken to be the DV score associated with its expressed OR, as measured via scRNA-seq.
Predictions of the DV score associated with the OR expressed in each cell were made in a cross-validated manner by holding out all the cells of a given OSN subtype at a time using either the DVISH score (linear regression with the top PC of DV expression, where PCA was fit on the training data) or the DV score associated with the ORs expressed in neighboring cells (with KNeighborsRegressor with n_neighbors=10 weighted by distance). Predictions were summarized at the OSN subtype level. Similar predictions were also observed using the nearest neighbor regression approach on a smaller sagittal dataset.
The complex 3D geometry of the nose made it difficult to compare data across sections, or across turbinates directly. An alignment procedure was developed to map data into common coordinate systems. Subsections of the epithelium, like the septum, had previously been shown to contain ORs from multiple zones. Using the expression of known dorsal and ventral ORs the boundaries of all segments of the epithelium that spanned from dorsal to ventral in the were manually annotated using the coordinates of the aligned sections post-alignment with STalign. Then, for each segment, the ranked distance of each cell from the most dorsal vertex of the polygon boundary of the segment was computed. These positions, expressed in terms of the fractional length along each segment, were considered as the “unaligned” coordinate system for downstream analyses.
If ORs are arranged continuously across the epithelium, then their positional ranks (along the dorsal to ventral axis of each segment) should be a good proxy for distance, assuming OR ranks are preserved across segments on average even if such ranks change with different rates across segments (e.g. due to differences in the total length of segments or their location within the MOE). The mean rank for each OR across all segments was calculated and used to develop a mapping between position (i.e. the fractional length along each segment) and the mean OR rank of the OR expressed at each position and these curves were fit with piecewise linear functions (with at most 5 knots and constrained to increase monotonically constraints) via the pwlfit library. The resulting piecewise linear regression thus gives a mapping, for each segment, between the positional coordinates of that segment and the common mean OR ranks; these predicted values were taken to be the new coordinates for each cell, thereby mapping each segment into a common coordinate space. Given the monotonicity constraint of this approach, the rank ordering of cells was relatively preserved by this approach. Accordingly, the overall mean OR ranks changed little by this approach, but the resulting positions in the new coordinate system were less confounded by differences in absolute positions across segments and better captured the underlying latent “position” of each cell, and in the aligned dataset the mean OR rank at a given position was more correlated. While this annotation and alignment procedure was done independently of the DV scores, both the mean unaligned and aligned positions were well-correlated with the DV score.
The resulting unaligned or aligned positions for the cells from each segment were used for downstream analyses. Common coordinate positions for each OR were averaged across all detected sections, and these mean positions strongly correlated (ρ = 0.969) with the DV score. Cross-validated regression models were used to quantitatively compare the predictability of position via the DV GEPs to that predicted by the zonal model. Three formulations of a five-zone model (i.e. linear regression with categorical features) were 1) with the standard zone indices (zones 1, 1–2, 2–3, 3–4, and 4–5); 2) the zones, but offset by half (≤1.5, 1.5–2.5, 2.5–3.5, 3.5–4.5, >4.5); or 3) the optimal bins from the training data, found using sklearn’s DecisionTreeRegressor with max_leaf_nodes=5. These models were compared to 1) a LinearRegression model that used the usage of the two DV GEPs and DV score for each OSN subtype; 2) via a KNN regression model (with uniform weights and 10 neighbors) fit on the same features in the training data, or 3) a LinearRegression model in which the two DV GEPs were expanded into a larger quadratic basis using the PolynomialFeatures(degree=2, include_bias=False) transformation. All models were fit via five-fold cross-validation, and the prediction accuracy (r2_score) for each model was summarized across 500 restarts.
AUROC analysis was performed between pairs of OSN subtypes, as described above for the scRNA-seq data. For analyses that looked at pairwise distances between cells, only pairs of cells within the same segment were considered, and pairs were accumulated and summarized across all segments. Dorsal ORs had broader and more intermingled spatial positions, which became compressed during our alignment procedure (since the mean ranks in the dorsal part of the epithelium were often flat). To avoid any confounding results, subsequent analyses were performed, where indicated, using only cells expressing ventral class II ORs.
Regression was also used to predict the position of each cell in the aligned coordinate system using each cell’s DVISH score. Predictions were performed using a Histogram-based Gradient Boosting Regression Tree (HistGradientBoostingRegressor with loss= “absolute_error”, max_leaf_nodes=15, max_depth=15). Equal numbers of cells were subsampled for each OSN subtype (for the set of OSN subtypes with at least 20 cells). The regressor was trained with five-fold cross-validation across 250 restarts and the median absolute error in the predicted positions on each restart was computed. The observed accuracy was compared to models trained on data in which DVISH scores were shuffled across cells within an OSN subtype or data in which the positions of each cell were shuffled across all cells.
Stereo-seq data were processed using the standard SAW pipeline^106^, which maps reads to the Ensembl v105 GRCm39 reference and decodes the spatial position of each barcode using the fiducial markings and tissue masks from the fluorescent ssDNA image. Downstream analyses were performed on the resulting gene expression file. OSNs were identified in data that was binned at 20 μm, using OR–bin pairs found within the epithelial layers that were detected in at least 3 UMIs. To assess how well DV scores mapped onto space, a similar spatial k-nearest neighbor approach was used, as described above for the MERFISH data. In brief, the DV score for each cell was predicted using the DV score associated with the ORs expressed in neighboring cells (the 10 nearest spatial neighbors, weighted uniformly by distance), using cross-validated KNN models that held out all cells for a given OR (for the 839 ORs detected in at least 6). Predictions were either averaged across all cells from all chips for each OSN subtype or averaged separated for each chip to calculate the correspondence in predictions across animals; predictions were also compared to those made using the MERFISH data, using the set of OSN subtypes identified with both technologies. For quantifying the AP genes via Stereo-seq, the expression of each gene associated with GEPAnterior and GEPPosterior in a 15 μm radius around each OSN was evaluated. OSNs were grouped into 10 deciles based on their associated AP scores, and the z-scored expression for anterior and posterior-associated genes was averaged for each decile. For visualization purposes, different sections from the same chip were also aligned to each other using affine transformations and STalign, as described above for the MERFISH data. In these aligned coordinates, the DV scores associated with each OR were then averaged in 50 μm bins.
DV scores were compared to two previous measurements, via bulk RNA-sequencing, of the spatial positions of ORs within the epithelium. One approach measured OR RNA levels in microdissected subregions of the epithelium^16^, whereas the other performed RNA tomography and measured OR RNA levels in cryosections across all three axes of the epithelium^17^. These spatial indices were compared to the mean DV score for OSNs expressing each OR in wild-type data (as measured via GEP usages, which are calculated using non-OR genes) for the sets of ORs detected both in these datasets and via scRNA-seq.
DV scores were also predicted based on the expression of ORs in manually-dissected “zonal” subregions of the MOE from a third dataset^28^. The raw fastq files for these experiments were obtained from the NIH sequence read archive (SRA, dataset SRP285789). Bulk RNA-seq experiments were uniformly processed and mapped to the same mouse index used for scRNA-seq experiments (Ensembl GRCm39 v105) using a Nextflow pipeline (nf-core/rnaseq) that ran STAR and RSEM to quantify the bulk RNA abundance of genes in each sample^112,130,131^. The DV score for each OSN subtype was predicted as a function of OR abundance in the zonal dissections using support vector machine regressors (C=500 and kernel=“rbf”), trained via 5-fold cross validation across 1,000 restarts.
MERFISH data from the MOE were generated in the Dulac lab and were publicly released (CC BY 4.0) as part of the NIH’s Brain Research through Advancing Innovative Neurotechnologies (BRAIN) Initiative - Cell Census Network (BICCN)^52^. Preprocessed data, which indicated the decoded olfactory receptor gene detected in each of ~3.3 million cells, were downloaded from the Brain Imaging Library (BIL, dataset ID: ace-bag-tin)^132^. This MERFISH dataset consisted of data from 10 mice in which coronal sections were probed with two OR gene panels that each expressed ~500 different ORs (and ~1,100 ORs in total). Analyses were performed on 49 sections from 8 mice, excluding 2 mice whose sections had lower quality. The epithelial sections were manually segmented as above, and the resulting 704 sections were also aligned to common coordinate system using the global OR ranks and the ranking of the cells in each segment. While the resulting ranks were well-correlated with the DV score overall, some ORs had appeared to be out of place compared to their neighbors, suggesting they might be the result of false-positive detections in this dataset. Lastly, to predict the DV scores of each OR in a segmentation-free manner, the same nearest-neighbor approach described above for the other MERFISH dataset was performed in which all the OSNs of a given subtype were held out and the DV scores associated with the held-out OR was predicted using the DV scores associated with the ORs expressed in the nearest neighbors spatially. Predictions were summarized at the OSN subtype level.
AP scores were also evaluated at the OSN subtype level in both MERFISH datasets using the AP scores associated with each OR in the integrated scRNA-seq dataset. To relate AP scores to epithelial positions, the apical and basal axes of each epithelial segment were manually annotated and the difference in the distance to the apical and basal axis (i.e. the side of the epithelium closer to the lumen/sustentacular cell layer and the side closer to the mesenchymal and basal stem cell layer, respectively) was computed for each OSN. Only OSNs whose distances to both axes were within the bottom 99^th^ percentile were included, to exclude outliers (likely false-positives outside of the epithelium or damaged areas of the tissue). The mean apical-basal distance for each OSN subtype was computed, which were then compared to the associated DV and AP scores for each subtype. The correlation in AP scores and apical-basal distances was also evaluated separately for dorsal and ventral ORs. The apical-basal positions of a subset of ORs were also manually verified via in situ hybridizations.
Data from wild-type mice housed in home-cage environments, or from mice that either underwent transient naris occlusion for a week or were housed in novel olfactory environments, were all previously generated^54^. DV scores were evaluated separately for each subtype for each condition (e.g. across open vs. closed nostrils or for mice housed in each environment), for OSN subtypes found in both conditions. In mice housed in home-cage environments, the DV score for each subtype was also compared to each subtype’s ES score, which summarizes the chronic activity induced by the odorants in the home-cage.
scRNA-seq was performed on female mice heterozygous for the CNGA2 allele. Because CNGA2 is X-linked its expression is mosaic in female heterozygotes because individual cells inactive either the wild-type or loss of function allele. Cells with the loss of function allele were identified either based on SNPs in the Cnga2 cDNA that distinguished it from the C57 allele, or via reads that mapped to lacZ. Cnga2 loss of function cells had decreased ES scores, consistent with the usage of that GEP being downstream of OR-driven activity.
Identification of activated OSN subtype in the Act-seq experiments was performed as previously described^54^. In brief, activated OSN subtypes were identified via immediate early gene induction, and the activation score for each subtype in the odor condition was calculated using a previous set of acutely-responsive “activation genes.” Activation scores were consistent across replicates. The DV scores for each activated OSN subtype were calculated using the entire set of cells in our integrated dataset.
P2 cells from OMP-P2 mice were identified via the expression of the P2 receptor, which was expressed in a large fraction of cells across the entire epithelium of mature mice. Two-OR cells were defined as those in which both ORs were expressed with at least 3 UMIs; many of these cells contained P2 in OMP-P2 mice, and P2 was often expressed at higher levels than the other non-P2 OR. DV scores in P2-expressing cells in the OMP-P2 mice were compared to those in our integrated dataset expressing P2 or to those expressing the non-P2 OR (for two-OR cells).
Cells expressing each OR in the OR-swap mice were identified based on the 3’ UTR associated with each OR transcript, which was unchanged in these mice and therefore identified the genomic locus of the expressed OR for each cell. The crossed design of the OR swap experiment (in which cells express either S50 or M72 receptors from each of the S50 and M72 genomic locus) facilitated the identification of genes that followed the locus or the expressed OR. Classification was performed using cross-validated SVM models (with kernel=“rbf”) trained on the expression of either the ~1,300 HVGs or the ~250 DV genes in each cell (with a feature selection step that kept the top 100 genes from the training data). Models were trained to predict the genomic locus of each cell using a subset of cells from each mouse and were tested on held-out cells.
Differentiating OSNs were also identified based on the transferred cell type labels of the scANVI latent space. Both mature and immature OSNs were excluded to obtain a dataset of differentiating cells in the process of choosing ORs. Iterative subclustering was performed on the latent space of new scVI models (n_hidden: 128, n_latent: 10, n_layers: 1, gene_likelihood: nb, dropout_rate: 0.1, with a batch key for each replicate and the number of genes as a continuous covariate) trained on the top variable genes from this dataset. In this process, any remaining doublets, low-quality cells, additional non-neuronal cells that were not excluded by the above filtering steps, and the small subset of cells from the non-canonical Emx1^+^ lineage that gives rise to Gucy1b2^+^ and Gucy2d^+^ sensory neurons were removed and the scVI models were retrained until the dataset consisted of only differentiating cells. The resulting differentiation dataset contained ~400k cells and densely captured cells from the earliest globose basal stem cells to late INP/early immature OSN cells that express single ORs. Note, although differentiating cells are labeled as INP to match prior literature^28^, cell cycle analysis suggests that many INPs are likely postmitotic. To remove the effects of the cell cycle on downstream analyses, cells scores for the G2/M and S phases were computed (scanpy’s score_genes_cell_cycle) using published marker gene sets^133^. Any gene whose expression was correlated with either G2/M or S score (pearson’s R > 0.2) was removed and the remaining 2773 were used to train a new scVI model with the same training parameters as listed above. Cells were clustered based on the nearest neighbor graph of the resulting 10-dimensional scVI latent space and this latent space was also used to reduce the dimensionality of the data to two dimensions for visualization purposes via UMAP (using a KNN-graph with 100 neighbors); these UMAPs are shown in the figures of the differentiating dataset.
Smooth differentiation trajectories were visible in the resulting UMAP projections, and thus an approach was designed to pseudotemporally order the cells. An approximate KNN graph was constructed based on the cosine similarity in the scVI latent space, using k=20 neighbors. A vector the size of the number of cells of the dataset was constructed. This vector was all zeros except for a single GBC cell, which was given a starting weight of 1. This vector was then multiplied by the kNN graph (weighting each cell by 0.1 and the influence of its nearest 20 neighbors by 0.9), renormalized from 0 to 1, and this procedure was repeated 1500 times until the values at each iteration remained stable. The resulting vector was rank-ordered, and these ranks were considered the pseudotime values for each cell; similar pseudotime values were also observed starting from other cells. A similar KNN-smoothing approach was used to smooth gene expression and GEP usages across neighboring cells, following procedures described previously^134^. Here the starting vector was the gene expression or GEP usages for each cell and the smoothing procedure was only run for 5 iterations, which results in locally-weighted values for each cell. The smoothed expression of example genes are shown in the UMAP representation, and the locally-smoothed version of each cell’s DV score was used for all downstream analyses, to help counteract the fact that fewer DV genes were expressed pre-choice and thus estimates at the cell level were noisier. Gene expression as a function of pseudotime was also fit via Generalized Additive Models using pyGAM. The alignment of DV score with axes of gene expression was assessed for cells of a given pseudotime bin (width = 0.1 and steps of 0.05) by taking the cosine distance of the DV score and the top principal components, where PCA was performed using the log-normalized expression of all OSN highly-variable genes in each pseudotime bin. DV gene expression as a function of pseudotime was also visualized for each gene, using pseudotime bins of width 0.05. Gene expression was z-scored across bins and then smoothed and sorted for visualization purposes using a rolling mean of 3 bins.
OR expression was considered with different thresholds. To identify multi-OR cells, a threshold of 1 UMI was used. Such a threshold likely induces some false positives (from e.g. background “soup” expression), given that both non-neuronal cells and GBCs sometimes expressed ORs at this 1 UMI threshold. Therefore, analyses of OR expression were restricted to cells expressing class II ORs with pseudotime values of at least 0.4, the value at which OR expression was consistently observed across cells. At subsequent pseudotime values, cells transiently expressed multiple ORs at the 1 UMI threshold. Cells with “competing” ORs were defined as those expressing multiple ORs at 3+ UMIs each (the threshold used for singular expression in mature OSNs). For cells expressing multiple ORs, the mean DV score of the co-expressed ORs was computed by weighting the DV scores associated with each OR (as measured in mature OSNs) by their expression levels (so an OR expressed at 2 UMIs would be weighted twice as much as one with 1 UMI). For individual ORs the median (and other percentiles) of the pseudotime values of all OSNs expressing that OR was computed, both for OSNs expressing each OR at any level, as well as for those expressing at either low (1–2 UMIs) or moderate levels (3+ UMIs).
To assess the correlation between each cell’s DV score and that associated with the OR, the correlation was evaluated both for cells at different differentiation stages as well as for OR expressed in those cells at either the highest, second-highest, or lowest (1-UMI) levels. Across 1,000 restarts, single ORs among the ORs expressed at each level in each cell (e.g. one of the 1-UMI ORs was picked for each cell) and the correlation was evaluated on each restart separately for cells at each differentiation stage (pre-choice with lots of low expression, cells with “competing” ORs, and cells expressing a single OR at high levels and others at lower levels).
To identify cells expressing barcodes and to extract the lentiviral barcode, raw sequencing reads that contained the constant start of the barcode (i.e. those that came from the amplified barcode libraries) were extracted from the BAM files. Only barcode reads with at least 10 reads for a given UMI were evaluated. The barcode library contained two variable regions (that were 14 and 30 bp long, respectively). Barcodes with hamming distances smaller than that observed between the whitelist of barcodes (3 for the 14 bp region and 5 for the 30 bp region) were collapsed and barcodes that were part of the whitelist were included in downstream analyses. Cells that expressed a single barcode (e.g. one barcode accounted for all of the barcode reads/UMIs for that cell) were considered for downstream analyses, and cells from the same replicate that expressed the same lentiviral barcode were considered as clonally-related. The resulting list of barcode-expressing cells from each experiment did not have any overlapping barcodes across mice, suggesting that the lentiviral libraries were sufficiently diverse to uniquely label each clone in each mouse. However, the overall UMIs per barcode were relatively low, and, even with the targeted amplification, many of the cells in the Venus^+^ libraries did not have any detectable barcodes. This might reflect false positives from FACS, or transcriptional silencing of barcode loci, especially given the extensive heterochromatinization of OSNs; the fluorescent signal from many infected cells was also quite weak, which is also consistent with silencing of the integrated cassette. Therefore, the modal number of cells per barcode was one, and the mode across clones with at least two cells was two. Nevertheless, results held true across clone sizes, further indicating that “small” clones are likely undersampled (due to the e.g. difficulties in barcode detection or in capturing all cells during the single-cell dissociated), but not fundamentally different than those with more cells. Additionally, the OSNs identified following methimazole-induced regeneration expressed single ORs at normal levels and had GEP usages that was similar to that of control OSNs that expressed the same OR. Similar results were observed for the lentiviral-infected cells, which expressed a wide array of ORs.
Clones expressing at least two OSNs, with each singularly expressing any OR or only class II ORs, as indicated, were considered for downstream analyses. DV scores for the expressed OR were calculated by using the DV score associated with that OR in the entire integrate dataset. The clonal mean OR DV was calculated by averaging the DV scores across clonally-related OSNs. Differences between the DV scores of clonally related ORs were calculated between all pairs of cells from the same clone and were either summarized across clones or across all pairs. Shuffles for all analyses were performed by shuffling clonal labels across cells. Clonal data was simulated based on the observed clonal restriction in DV scores by sampling, for each of the cells in the observed clones, cells from the wild-type data whose ORs had DV scores that matched the observed t-distribution of deltas from the clonal mean (e.g. ORs with DV scores similar to the clonal mean DV score had the highest probability of being chosen); the probabilities were also scaled by the observed frequency of each OR in the wild-type data. The percent of clones that had multiple cells that expressed the same OR (and therefore fewer number of ORs than cells) were then analyzed for both the observed and simulated data.
DV scores were also evaluated at the cell level, in which each cell’s DV scores was obtained by taking the difference in GEPDorsal and GEPVentral usage irrespective of the expressed OR. To test whether restrictions at the cell level were stronger than that at the OR level, the MAD in cell DV scores within a clone was compared the MAD obtained when DV scores were shuffled across cells of a given OSN subtype (irrespective of their clonal identity) or shuffled across all cells.
INP gene expression was related to OR DV scores, using the subset of clones that contained both INP cells and OSNs. For the cells in these clones support vector machine regression models (SVR with C=50 and kernel=“rbf”) were fit on the INP gene expression (where the input to the linear models was the top 250 genes in the training data, reduced to 30 dimensions via PCA) to predict the mean OR DV score of the clonally-related ORs expressed in the sister cells. The regression was performed via cross-validation, leaving out all the cells from a clone; additionally, the top genes and PCA transformation was identified using only the cells of the training data to avoid any data leakage.
Data from adult mice given either the ALDH1A2 inhibitor WIN 18,446 (RA inh.) or all trans-Retinoic Acid (atRA) during methimazole-induced regeneration were analyzed relative to their respective vehicle controls. Cell types were identified via clustering the latent space of scVI models and cNMF GEP loadings were applied to each cell, as described above for the integrated dataset. Analyses were either performed at the cell level, using the observed DV score for each cell, or at the OR level, using the DV score associated with each OR as measured using all cells in the integrated datasets. To assess distributional shifts, the median DV score across all cells or the number of cells expressing ORs with given associated DV scores were compared between the data from each condition and its respective control; given that fewer cells express the most ventral ORs to begin with, the changes in OSN frequency were also expressed as log2 fold-changes. Because RA manipulations changed the overall distributions of cells without altering the relationship between cell identity and OR choice, distributional shifts were observed when evaluating all OSN subtypes at the cell level, but minimal changes were observed across conditions for cells of a given OSN subtype. Changes in gene expression were evaluated by z-scoring the expression of each DV gene across all cells and the mean change in z-scored expression was evaluated for cells from drug versus control conditions. Analyses were performed on mature OSNs expressing single ORs, except where, as noted, INPs were used to evaluate changes in DV scores and OR expression in differentiating cells pre-choice.
ChIP-seq data were previously generated and were reanalyzed^28,64^. The raw fastq files for ChIP-seq data for H3K9me3 and H3k79me3 marks were obtained from the SRA (datasets SRP285789 and SRP096660) and were reprocessed with the Ensembl GRCm39 version 105 mm39 genome, using the default parameters of the nf-core/chipseq Nextflow pipeline (v2.0.0)^112^. As part of this pipeline, reads were mapped to the genome with BWA, and the normalized read density was summarized at base-pair resolution with UCSC-bedGraphToBigWig. The density of reads across the entire coding and non-coding regions of each OR for each sample were averaged, and the reads for each OR locus (for class II ORs) were correlated with the DV score of OSNs singularly expressing that same OR. The samples evaluated in these analyses were of sorted GBCs, INPs, iOSNs, or mOSNs that were sorted from the entire MOE (and thus contain a mix of cells from the dorsal and ventral regions). The average heterochromatin from these bulk dissections was also compared to those found in Olfr1507+ cells (sorted from Olfr1507-IRES-GFP animals).
Micro-C data was analyzed following standard approaches. In brief, paired end fastq files were aligned to the Ensembl v105 GRCm39 genome using bwa (bwa mem −5SP -T0). Valid ligation junction events were identified using pairtools parse (--min-mapq=40, --walks-policy=5unique, and --max-inter-align-gap=30), and then sorted (pairtools sort), deduplicated (pairtools dedup), and split (pairtools split) into pairs and pairsam files. Contacts in the pairs files were summarized at a 1kb resolution with cooler and converted into multi-resolution mcool files. To best analyze interchromosomal contacts for individual ORs, data from all biological and technical replicates were combined for downstream analyses, but results were consistent across samples. Using data binned at the 1kb resolution, the ± 5kb region from the start of each class II OR gene was considered as an OR locus, and the number of contacts between that locus and all other OR loci were summarized. To avoid distance-dependent effects, downstream analyses were restricted to interchromosomal (trans) contacts. Class II ORs were grouped into deciles based on their associated DV scores, and the mean pairwise (10 × 10 deciles) or one to all (10 deciles) contact frequency were evaluated.
Similar results were also observed via Hi-C from OMP-IRES-GFP mice, and from data from individual OR lines. Hi-C experiments were processed with the Nextflow nf-core/hic pipeline (with digestion== “dpnii”), and were aligned to the same reference genome using bowtie2. The resulting multi-resolution mcool files were then evaluated as described above.
scRNA-seq was performed as described above in the HP1 swap and control animals (i.e. with or without the Foxg1-Cre allele) in young (4–6 week) animals of either sex. Cell type identification, DV scores, and changes to the cell and DV scores in the HP1 swap animals were evaluated in a similar manner to those in the RA manipulations. To assess changes in OR choice for cells of a given DV scores, cells were binned based on their cell DV scores and the percent-normalized DV score associated with each OR was computed for cells in each bin. The cell vs OR DV score mapping was evaluated using robust linear regression (HuberRegressor), and equal numbers of cells were sampled for each OSN subtypes, for subtypes detected in both control and HP1 swap animals. To evaluate changes in the cell to OR DV score mapping, a nearest neighbor approach was used, and a nearest neighbor regression model (KNeighborsRegressor) was fit on all cells from the control mice to predict the DV score associated with the ORs expressed in cells of a given DV score. This model was then applied to the cell DV scores of the HP1 mice, and the residuals from this model, which capture the difference in the associated percentile-normalized DV score of the OR detected in each HP1 swap cell with that predicted based on each cell’s DV score, were evaluated and summarized for HP1 swap cells of a given DV score.
scRNA-seq was performed in F1 hybrid animals generated by crossing wild-derived CAST/EiJ female mice with C57BL/6J males. CAST/EiJ mice have on average single nucleotide polymorphisms (SNPs) every 150 bp. Because OR expression is monogenic and monoallelic, the reads that mapped to the chosen OR were used in each mature OSN to infer the strain of the chosen OR. The GRCm39 coordinates of SNPs that distinguished CAST/EiJ and C57BL/6J were obtained from the Mouse Genome Project^75,135^, and SNPs in ORs of interest were also confirmed by evaluating bulk and scRNA-seq data from homozygous CAST/EiJ animals, generated in house and in past work (Ibarria-Soria). For each cell, any reads (with MAPQ > 30) that overlapped with homozygous SNPs in CAST/EiJ mice for its chosen OR were evaluated to assess what fraction of reads expressed the CAST/EiJ or C57BL/6J variant. Cells that had reads that overlapped at least one SNP and in which at least 80% of such reads came from a single allele (mean 99.7% of reads) were used for downstream analyses. Similar results were also obtained using Demuxalot to demultiplex the strain of each cell via the set of SNPs across all OR genes^136^. The strain of the chosen OR was able to be inferred in 81% of cells. The remaining 19% of cells had too few reads that overlapped SNPs, likely because they expressed ORs with minimal variation between strains in their 3’ UTRs. SNPs that led to missense mutations were annotated via SnpEff^137^.
OSN subtypes with significant changes in GEP usages were evaluated empirically using permutation testing, and subtypes with a mean change between strains larger than 1% of shuffles for that GEP were considered as significant. Similar changes were observed across F0 and F1 animals, as well as across individual samples; however, some OSN subtypes had few cells for either CAST/C57 (as OR frequency varies across strains) and were likely underpowered.
Demultiplexed fastq files from CUT&RUN samples were mapped to the Ensembl GRCm39 v105 genome via Bowtie2 (-I 10 -X 700 --no-mixed --end-to-end --no-discordant --very-sensitive)^113^. Spike-in reads from E. coli DNA were mapped to the E. coli genome. However, because F1 animals have within-animal controls (the reads mapped to each allele), the number of reads that mapped to the E. coli genome were not used for normalization. Duplicates were not removed, and data were combined across replicates. As in the F1 scRNA-seq data, the allele of each read was inferred for reads overlapping SNPs that distinguished CAST/EiJ and C57BL/6J. The number of reads that mapped to each allele were then summarized for given ORs across the entire genomic coordinates for that gene (i.e. from TSS to TES).
Promoter motif enrichment was performed using HOMER^138^, using the 1kb region upstream of the TSS (i.e. the promoter region) for each OR (findMotifsGenome.pl with -size=given). This approach confirmed the enrichment of known OR promoter motifs (i.e. EBF, the homeobox TF LHX, and ATF); the de novo Motif results also identified additional variants of EBF and LHX motifs that were enriched in the OR promoters. EBF motifs were especially enriched in the 0–250 bp upstream of the TSS, as expected for promoter motifs. Next, differential motif enrichment between the top 200 dorsal and ventral class II OR promoters was evaluated using HOMER and MEME, but this approach merely returned the presence of EBF and LHX motifs, suggesting that it may have been data limited. Finally, supervised promoter annotation was performed using HOMER (annotatePeaks.pl) using the default motif files and threshold provided by HOMER for LHX2, EBF, ATF (given their general promoter enrichment identified above) and for NFI, and RARa (given their presence in the DV GEPs). This approach revealed dorsal OR promoters were enriched for EBF binding motifs in, and NFI/RAR motifs were more-likely to be found in ventral ORs. Comprehensive scanning for motifs for additional TFs and factors (including Mef2a, Emx2, Sox2, CTCF, AP-1, and other common promoter motifs, all via their default motif files provided by HOMER) failed to identify other TFs with dorsoventral biases, suggesting that the factors above may be sufficient, or such differences may not be captured by the differential presence of putative binding motifs. The odds ratio for the co-occurrence of EBF, NFI, and RAR motifs was evaluated empirically via bootstrapping.
The olfactory bulb glomerular map was generated via spatial transcriptomics using the 10x Visium platform, as recently described by Klimpert and colleagues^79^. In brief, putative glomeruli were identified using a Bayesian approach that used the expression of individual OR genes across the entire 3D volume to infer the location and number of glomeruli for each OR. The outputs of this model were manually curated to keep high-confidence identified glomeruli, and an axis of symmetry was identified that separated glomeruli from the medial and lateral domains. Importantly, the two hemibulbs are also symmetric, and once aligned to each other sister glomeruli expressing the same OR were found on average 230 μm away from each other (~1–2 glomerular lengths), reflecting a combination of both biological variability and spatial errors in the glomerular detection and alignment procedures; such deviations are also similar to those from previous estimates^34^.
Glomeruli were identified based on OR gene expression and could thus be linked to the scRNA-seq or MERFISH data for the associated OSN subtype expressing that OR. First, OB positions were used to predict the scRNA-seq measured GEP scores for each OSN subtypes, via a cross-validated support vector machine regression model (SVR with C=150 and “rbf” kernel). Predicting GEP usages from glomerular positions is easier than the reverse, due to the higher dimensionality of the 3D positions and due to the fact that the DV score mapped onto an axis in the OB that was not fully aligned with the cartesian D-V axis but rather correlated with both the D-V and A-P position of glomeruli. Similar results were also obtained through canonical correlation analysis, which identified optimal rotations to align the GEP and glomerular spaces, further demonstrating that the DV score was well-aligned (ρ > 0.9) with the top canonical correlation axis. Second, predictions of the 3D glomerular position for each OSN subtype was performed separately for glomeruli in each domain. Elastic-net-regularized linear regression models (alpha=1, l1_ratio=0.9) were trained on either GEP usage (e.g. GEPDorsal, GEPVentral, and the DV score) or all ~1,300 variable genes to predict the glomerular position along each axis and the 3D error of these predictions was evaluated for each glomerulus. For both models, five-fold cross validation was performed; this procedure was repeated 100 times, and the distribution of the accuracies and or errors across restart was reported.
Single-cell analyses were performed in python (versions 3.6–3.12) using the Scanpy and scvi-tools packages^107,108^, as well as custom-written scripts using the open-source python scientific stack (pysam, SciPy, NumPy, scikit-learn, umap-learn, pandas, statsmodels, matplotlib, seaborn, and numba). cNMF was performed using modified versions of the code from https://github.com/dylkot/cNMF.
Hypothesis testing was performed using parametric or non-parametric statistical tests, as indicated, or p-values were calculated empirically using resampling-based permutation tests. The precision of sample statistics and regression trend lines were evaluated using bootstrapping, and, except where noted, plots and error bars depict the mean and the 95% confidence intervals of the mean across 1,000–10,000 bootstraps. Throughout the paper, a nonparametric version of the box plot (also known as a letter-value plot) was used to represent multiple quantiles and tails of large distributions of data (e.g. to summarize across OSN subtypes or sets of ORs or OR pairs) in an agnostic way that does not require setting bandwidth parameters as in violin plots or kernel density estimates. Like a conventional box plot, the largest box represents the interquartile range (25–75 percentile) and the median (dotted-line). Subsequent boxes recursively represent exponentially-smaller quantiles (the 12.5–25 and 75–87.5 percentiles, then 6.25–12.5 and 87.5–93.75 percentiles, then 3.125–6.25 and 93.75–96.875 percentiles, and so forth. The results of additional statistical tests are listed in^57^S5.