Concentration-Dependent Recruitment of Mammalian Odorant Receptors

Abstract A fundamental challenge in studying principles of organization used by the olfactory system to encode odor concentration information has been identifying comprehensive sets of activated odorant receptors (ORs) across a broad concentration range inside freely behaving animals. In mammals, this has recently become feasible with high-throughput sequencing-based methods that identify populations of activated ORs in vivo. In this study, we characterized the mouse OR repertoires activated by the two odorants, acetophenone (ACT) and 2,5-dihydro-2,4,5-trimethylthiazoline (TMT), from 0.01% to 100% (v/v) as starting concentrations using phosphorylated ribosomal protein S6 capture followed by RNA-Seq. We found Olfr923 to be one of the most sensitive ORs that is enriched by ACT. Using a mouse line that genetically labels Olfr923-positive axons, we provided evidence that ACT activates the Olfr923 glomeruli in the olfactory bulb. Through molecular dynamics stimulations, we identified amino acid residues in the Olfr923 binding cavity that facilitate ACT binding. This study sheds light on the active process by which unique OR repertoires may collectively facilitate the discrimination of odorant concentrations.


Introduction
In nature, odors coming from foods, mates and predators are dynamic in concentration. The ability to discriminate odorants at different concentrations while recognizing the same odorant across different concentrations serves nutritional, reproductive and protective purposes, allowing animals to make appropriate behavioral decisions to maximize their fitness. The mammalian olfactory system, capable of detecting and discriminating many odorous volatile molecules, begins at the nasal passage, where millions of olfactory sensory neurons (OSNs) line the main olfactory epithelium. Odor detection is initiated by the activation of a specific set within the large repertoire of G-protein-coupled odorant receptors (ORs), encoded by ;400 and ;1100 intact OR genes in humans and mice, respectively (Buck and Axel, 1991;Healy et al., 1997;Zhang and Firestein, 2002;Godfrey et al., 2004;Malnic et al., 2004;Ibarra-Soria et al., 2014;Niimura et al., 2014). Mature OSNs express only one type of the available ORs at a high level (Chess et al., 1994;Malnic et al., 1999;Hanchate et al., 2015;Saraiva et al., 2015;Tan et al., 2015;Scholz et al., 2016). As OSNs expressing the same OR send convergent axonal projections to form glomeruli in the olfactory bulb, OR activation is translated into glomerular activation patterns which are further processed within the olfactory bulb and the olfactory cortical areas (Ressler et al., 1994;Vassar et al., 1994;Mombaerts et al., 1996;Rubin and Katz, 1999;Wachowiak and Cohen, 2001;Oka et al., 2006;Miyamichi et al., 2011;Sosulski et al., 2011;Storace and Cohen, 2017;Wilson et al., 2017). Odorants are represented by combinatorial codes whereby a given odorant, at a specific concentration activates a specific combination of ORs, which in turn activate a specific combination of glomeruli (Malnic et al., 1999;Rubin and Katz, 1999;Bozza et al., 2002;Oka et al., 2006;Saito et al., 2009). Previous studies found an increase in the recruitment of active ORs (in the OSNs) and glomeruli (in the olfactory bulb) in response to higher odorant concentrations (Rubin and Katz, 1999;Fried et al., 2002;Khan et al., 2010;Jiang et al., 2015;Wilson et al., 2017). There is, however, insufficient research determining the identity of these ORs and their collective responses at different odorant concentrations in vivo. This lack of research limited our ability to draw conclusions on how odor identity and intensity information is encoded by specific ORs.
To address this limitation, we identified ORs activated by acetophenone (ACT) and 2,5-dihydro-trimethylthiazoline (TMT) at varying odorant concentrations in vivo. This was made possible through a recently developed in vivo-based method which identifies ORs expressed in activated OSNs via mRNA profiling (Jiang et al., 2015). This method makes use of the finding that the phosphorylation of the ribosomal protein S6 is an indicator of neuronal activation (Knight et al., 2012;Biever et al., 2015). Immunoprecipitation of ribosome-mRNA complex containing phospho S6 (pS6-IP also known as phosphoTRAP) followed by next-generation sequencing (pS6-IP-Seq) identifies OR mRNAs present in activated OSNs (Fig. 1).
Besides pS6-IP-Seq, there are two additional approaches to study OR responses in vivo: one is the Kentucky method based on activation-induced S100a5 expression (McClintock et al., 2014) and another is the DREAM technique (deorphanization of receptors based on expression alterations of mRNA levels) which is based on OR-specific transcriptional downregulation after prolonged odor stimulation (von der Weid et al., 2015). It has been previously reported that DREAM and pS6-IP-Seq identified overlapping, but not identical sets of ORs (Jiang et al., 2015;von der Weid et al., 2015).
In this study, we used a combination of in vivo, in situ, and in silico approaches to investigate ORs with distinct sensitivities to the tested odorants. In addition, we compared data obtained by pS6-IP, DREAM and heterologous expression methods to clarify the degree of consistency in the ACT responsive ORs identified by these methods. We examined Olfr923, which we identified to be one of the most sensitive ACT ORs based on our pS6-IP-Seq data. To this end, we created a gene knock-in mouse line to test activation of the Olfr923 glomeruli in the olfactory bulb. Additionally, we performed molecular modeling to evaluate ligand binding and amino acid residues involved in ACT's activation of Olfr923.

Animals
All animal procedures were performed in accordance with the Duke University animal care committee's regulations. All experiments in this paper were conducted on both male and female mice. C57BL/6J and B6.Cg-Gt (ROSA)26Sor tm9(CAG-tdTomato)Hze/J were obtained from the Jackson Laboratory and crossed in the mouse facility. C57BL/6J mice of 20-22 d of age were used for pS6-IP-Seq, RNA-Seq and staining experiments. Littermates of the same sex were used for each pS6-IP-Seq replicate (n = 3) and data from both sexes were pooled to avoid sexspecific genes being enriched in a given odor exposure. Although we did not formally test between the two sexes for differences, we did not observe any obvious differences between males and females in the experiments.

Odor exposure
Odorant exposures of ACT or TMT were performed on C57BL/6J mice as described (Jiang et al., 2015). Mice were habituated in a disposable paper container for 1 h This work was supported by National Institutes of Health Grants DC014423 and DC016224 the National Science Foundation Grant 1556207 and the Kahn Neurotechnology Development Grant. X.S.H. was supported by the Holland-Trice Graduate Student Fellowship. K.W.Z. is supported by 5F31DC017394. and exposed to an odor cassette for 1 h. The starting concentrations of odorants used for stimulations contained 10 ml of tested odorant diluted with dW (v/v). Diluted odorants were spotted individually onto a 2 Â 2-cm blotting paper and placed into a Uni-cassette (Sakura) for odor stimulation. The mice stimulation environment was inside a fresh paper container with lid (International Paper) within the fume hood.

pS6-IP-Seq
Following odor stimulation, the mouse was sacrificed, and the dissected olfactory epithelium was processed for RNA library in reference to the method described by Jiang et al. (2015), with some modifications. Pieces of freshly dissected olfactory epithelium were homogenized and centrifuged. 6 ml of a monoclonal phospho-S6 ribosomal protein (Ser240/244; D68F8) XP rabbit mAb (Cell Signaling Technologies) was added to the supernatant, and incubated with rotation for 1.5 h at 4°C. Protein A Dynabeads (Invitrogen) were added to the antibody mixture and incubated for another hour with rotation at 4°C. After incubation and washing, the bead-captured RNA was eluted and purified with RNeasy Micro kit (QIAGEN). 10 ng of starting RNA was used to make cDNA libraries with the SMART-Seq v4 kit (Clontech) with 11 cycles of PCR. cDNA was purified with AMPure XP beads (Beckman Coulter) and 0.5 ng was used to make cDNA libraries using the Nextera XT DNA Library Preparation kit (Illumina). 50 base sequencing was performed on the HiSeq 2500 System (Illumina), obtaining 15-38 million reads per sample. Sequence reads were mapped using STAR alignment (Dobin et al., 2013) against the entire transcriptome and newly annotated ORs [untranslated regions (UTRs) annotations included for all genes]. Reads mapped onto 1088 intact OR genes (297 ORs with pseudogene labels removed) were used for the purpose identifying odorant activated OSNs based on differential expression analysis with edgeR (Anders and Huber, 2010;Robinson et al., 2010;McCarthy et al., 2012; for R code, see Extended Data 1). ORs with false discovery rate (FDR) corrected p value ,0.05 and a positive fold change were considered significantly enriched. Raw reads and quantification results can be accessed at GEO: data submitted and processing.

OR mapping and annotation
Based on the Jiang et al., publication, mRNAs captured by pS6-IP were mapped based on the Mus musculus (house mouse) genome assembly GRCm38 (mm10) coding exons using Bowtie (Jiang et al., 2015). OR transcript definitions were re-annotated to include 59 and 39 UTRs with the coding sequence. This is based on Ibarra-Soria et al's work which included new annotations for 1248 mouse ORs (Ibarra-Soria et al., 2014). University of California Santa Cruz genome assembly mRNA sequences were used to visualize transcripts that mapped onto coding and non-coding parts of the exons for all genes  al., 2016). In total, we included the gene definitions of 1385 ORs (pseudogenes included) in our new pS6-IP-Seq pipeline. To analyze the Illumina sequencing read files, we created a customized GTF annotation file. BLAT v. 35 (Kent, 2002) which was used to align intact OR sequences to the mouse genome represented by the Mus_musculus.GRCm38.dna_sm.primary_assembly.fa.gz file from Ensembl release 96 (April 2019). A custom R script was used to remove hits that did not overlap with the coding region of the queried OR and format the results in GTF format. OR pseudogenes and non-OR genes were defined using the ensembl standard file Mus_musculus.GRCm38.96.gtf.gz file. Snakemake v3.5.5 (Köster and Rahmann, 2012) was used to process read files through alignment and quantification. STAR v2.7.0d (Dobin et al., 2013) was used to generate a genome index using the Ensembl Mus musculus primary assembly genome sequence and custom GTF annotation file. Reads were aligned to this genome index using STAR with default options except for -quantMode TranscriptomeSAM which maps genome alignments to transcript coordinates. STAR output transcriptome SAM files were quantified using RSEM v1.3.1 (Li and Dewey, 2011) using default options to generate gene and transcript level counts for differential expression.

RNA-Seq from dorsal and ventral olfactory epithelium
After the dorsal and ventral olfactory epithelium were dissected, RNA was extracted using the TRIzol reagent (Invitrogen) followed by clean up with RNeasy Mini (QIAGEN). A total of 1000 ng of RNA were used for cDNA libraries using the SMART-Seq v4 kit (Clontech), with two cycles of PCR. cDNA was purified with AMPure XP beads (Beckman Coulter) and 0.5 ng was used to make cDNA libraries using the Nextera XT DNA Library Preparation kit (Illumina). The sequenced reads were then mapped and analyzed using the same bioinformatics pipeline as described above for pS6-IP-Seq. Out of the 117 unclassified ORs, 45 ORs were unclassified due to low expression. The vast majority of the remaining ORs (70/73) are assigned to zonal index positions between 1.5 and 2.5 (1 being most dorsal and 5 being most ventral) in Tan and Xie (2018), suggesting that the these unclassified ORs are expressed in areas close to the dorsal/ventral boundary (Extended Data Fig. 2-1).

In situ hybridization and immunohistochemical staining
Fluorescent in situ hybridization (FISH) and immunohistochemical staining were performed in reference to Jiang et al. (2015) with some modifications. Briefly, 18-mm frozen sections of the mouse snout containing the olfactory epithelium were fixed in paraformaldehyde and pretreated. Digoxigenin (DIG)-labeled complementary RNA probes were hybridized to the tissue sections overnight at 58°C. Unless specified, we used open reading frame of ORs for the RNA probes. We designed a 3'UTR specific probe for the labeling of Olfr376 because of open reading frame sequence similarity to other ORs, and the primers used for cloning into pCI vector were as follows: AAACGCGTCAGTAATATTTTAA CACTGA (forward) and AAGCGGCCGCCTGAGTGTACAGTTTTTGAG (reverse). Following washes, the sections were incubated for 45 min with a horseradish peroxidase (HRP)-conjugated antibody (Roche, 1:1000 in blocking solution) against DIG. Hybridization signals were labeled with a 10-min incubation of tyramide signal amplification (TSA) solution. For immunostaining, the FISH-labeled sections were incubated with a polyclonal rabbit anti-phospho-S6 (244/247; Thermo Fisher Scientific, 1:300 dilution in blocking solution) overnight at 4°C. Following washes, a 45-min incubation with the donkey Cy3-conjugated anti-rabbit IgG (Jackson ImmunoResearch Laboratories, Inc., 1:200) was used to visualize pS6 proteins. Olfactory epithelium images were taken using the Zeiss Axiocam MRm and upright inverted fluorescent microscope with ApoTome functionality at 200Â magnification. The filter sets used were as follows: Zeiss filter set #38 for fluorescein, #43 for Cy3, and #49 for bisbenzimide. ImageJ was used in the quantification of colocalization signals, where pS6 pixel intensity of DIG-positive OSNs was subtracted from the average pixel pS6 intensity of no odor exposed OSNs labeled with the same OR probe. The resulting normalized intensities were quantified and compared using one-way ANOVA with Dunnett's test.

Olfr923-CRE transgenic mouse
Targeting of the Olfr923 locus A targeting vector expressing Cre under control of the Olfr923 locus was generated by BAC recombineering (Liu et al., 2003) in SW102 cells (Warming et al., 2005). Briefly, an IRES-Cre-frt-neomycin-frt cassette was introduced by homologous recombination into BAC clone bMQ186d17 (SystemBioSciences) containing 147 kb of mouse chromosome 9 immediately downstream of the TGA stop codon within the Olfr923 coding region. Next, the targeting vector containing the IRES-Cre-Neo flanked by homology arms was retrieved by gap repair from the BAC into plasmid pL253. The resulting targeting vector was linearized with NotI restriction endonuclease and purified by phenol:chloroform extraction and ethanol precipitation. Murine ES cell line was electroporated with a linearized targeting vector containing 7 kb of 59 homology and 2 kb of 39 homology to the endogenous Olfr923 locus (G4 mES Cells were acquired from Samuel Lunenfeld Research Institute Mount Sinai Hospital).

ES targeting
A total of 15 million G4 ES cells (Nagy et al., 1993) passage 13, were electroporated with 35 mg of targeting construct in a 4 mm gap cuvette using a Bio-Rad Genepulser II with Bio-Rad Capacitance Extender, set at 0.25 KV, 0.5 UF Â 1000 for two pulses. Transfected cells were plated onto four 10 cm plates with neomycin resistant feeder cells. Targeted cells were fed selection medium with 250 mg/ml active G418 Geneticin at 24 h after transfection, followed by the addition of two Um ganciclovir to the selection medium at 48 h. Selection medium with G418/Ganc was changed once a day until clones were picked on day 7 after transfection. A total of 142 clones were picked and expanded on 96-well plates. DNA was analyzed by PCR for homologous recombination on the 39 end. 24/142 clones were identified as positive. Eighteen clones were expanded to the six-well size plate, then cryopreserved in duplicate aliquots. Further analysis was done on genomic DNA prepared from these expanded clones by 39 and 59 Southern blotting and PCR to confirm homologous recombination. Confirmed clones were injected into an ICR morulae host embryo (Nagy et al., 1993). Five clones produced a high percentage of agouti coat color chimeras which were used for breeding and expansion of the colony for this study.

Olfactory bulb immunohistochemistry
Three-to six-week-old mice were habituated for 1 h and then received either control or ACT stimulation by spotting 10 ml of freshly made 0.01% ACT onto blotting paper placed in an odor cassette for 1 h. Control conditions consisted of placing clean blotting paper into an odor cassette. The mice were sacrificed and olfactory bulbs were collected into cold fixative solution consisting of 1Â PBS, 4% paraformaldehyde for 2 h with gentle shaking. Subsequently, olfactory bulbs were washed with cold 1Â PBS three times and placed into cryoprotective solution consisting of 1Â PBS, 30% sucrose, overnight with gentle shaking (Nutator S0500-VWR). The next day, sunken olfactory bulbs were collected and placed into optimal cutting temperature (OCT) compound and flash frozen with liquid nitrogen. Olfactory bulbs were then sectioned at 20-mm sections using a cryostat. Detection of endogenous TdTomato fluorescence using 554/581 nm (excitation/emission) light was used to confirm glomerulus presence whilst sectioning. For fluorescent labeling of c-Fos-positive periglomerular cells, olfactory bulb sections were immersed in 1Â PBS for 5 min at room temperature. Slides were then blocked with 4% donkey serum (Jackson ImmunoResearch 017-000-121) and 1% Triton X-100 for 1 h at room temperature. Subsequently, slides were incubated with the primary anti c-Fos antibody (1:400, 9F6 Cell Signaling 2250) overnight at 4°C. The next day, slides were washed in 1Â PBS three times at room temperature. Slides were then incubated with secondary donkey anti-rabbit Alexa Fluor 488 (1:200, Jackson Laboratories) diluted in 1Â PBS, 5% skim milk for 1 h. Slides were then washed in 1Â PBS two times, stained with bisbenzimide (1:1,000,000) for 5 min, washed with 1Â PBS, distilled water, and mounted with mowiol. Glomeruli with the greatest density of tdTomato expressing OSN terminals were imaged and periglomeruli cells in the immediate vicinity of the Olfr923 glomerulus were counted. C-Fos signal intensities were quantified as follows. We obtained 200Â magnification Z-stacked images with 2-mm intervals using the Zeiss Axiocam MRm and upright inverted fluorescent microscope. The filter sets used were as follows: Zeiss filter set #38 for Alexa Fluor 488, #43 for tdTomado, and #49 for bisbenzimide.

GloSensor assay
The liquid-stimulation GloSensor assay was conducted as previously described (Kida et al., 2018) with slight modifications. Hana3A cells were plated onto 96-well plates and placed into a 5% CO 2 incubator overnight; 18-24 h after plating, the cells were transfected with 80 ng/ well of Olfr923-pCI plasmid, 5 ng/well of RTP1S and 10 ng/well of GloSensor plasmid (Promega) and placed into a 5% CO 2 incubator overnight; 18-24 h later, the medium was replaced with 50 ml of HBSS (Invitrogen) containing 10 mM HEPES and 1 mM glucose (wash step). This was followed by adding 30 ml of the HBSS containing GloSensor cAMP Reagent (Promega). Plates were kept in a dark place at room temperature for 2 h to equilibrate the cells. ACT and heptanal were diluted individually in DMSO (Sigma-Aldrich) to 100 mM working stocks, and further diluted in the GloSensor buffer (Promega) for dose-response experiments; 30 ml/well of diluted liquid odorants were added and the test plate was immediately inserted in the plate reader. The luminescence in each well was measured at 90-s intervals for six cycles. All luminescence values were normalized by dividing by the value obtained from the wells transfected with the empty pCI vector in the same cycle.

Homology modeling
The protocol follows two previously published methods with an optimal alignment which fulfills constraints provided by .100 mutants from the literature (Charlier et al., 2013;Yu et al., 2015). Four experimental non-olfactory GPCR structures (1U19, 3ODU, 2YDV, and 2LNL) were selected as templates to build Olfr923 by homology modeling using Modeller (Eswar et al., 2006). The percentage of identity for the seven transmembrane domains of Olfr923 with respect to the four templates were as follows: 22.02% with Rhodospin (1U19), 13.68% with CxCR4 (3ODU), 19.87% with A2a (2YDV), and 25.00% with CxCR1 (2LNL). The N-terminal structure was omitted to avoid perturbing the modeling protocol. Five models were obtained and we kept the one that was consistent with several additional structural constraints (no large folded structure in extra-cellular loops should be observed, all TMs and H8 folded as a-helices, a short a-helix structure between TM3 and TM4).

Docking of ligand
ACT and heptanal structures were prepared with the antechamber module of AMBER with AM1-BCC charges. They were docked into the receptor cavity with Autodock Vina with an exhaustiveness parameter of 50. An exhaustiveness of 15 produced the same resulting poses (Trott and Olson, 2010). For each ligand, we considered simulations with the lowest binding free energy for both molecules.

Molecular dynamics simulation
All systems were then embedded in a model membrane made up of POPC lipids solvated by TIP3P water molecules using Maestro (Schrödinger, 2013). The total system was made up of ;77,000 atoms. For all systems in this article, molecular dynamics simulations were performed with sander and pmemd.cuda modules of AMBER with the ff03 force-field for the protein the gaff.lipid force-field for the membrane, and the gaff force-field for the ligands (Case et al., 2014). Bonds involving hydrogen atom were constrained using SHAKE algorithm. Long-range electrostatics interactions were handled with the particle mesh Ewald (PME) method. The cutoff for non-bonded interactions was set to 8 Å. Temperature was maintained constant with a Langevin thermostat with a collision frequency of 2 ps À1 . In addition, a weak coupling anisotropic algorithm with a relaxation time of 1 ps À1 was applied to keep a constant pressure. Snapshots were saved every 2 fs. The detailed workflow of the simulations was detailed in Extended Data Figure 7-1C.

Free energy calculation
The free energy of binding between the odorant and the OR was evaluated using the MM-GBSA method, according the following equation: The ,. corresponds to an average of a given value extracted from the MD trajectory.
Each term of the Equation 1 can be written as: where X corresponds to the complex, the odorant or the OR. The three terms were the gas-phase contribution to the binding energy, the solvation free energy on binding and an entropic contribution, respectively. The gas-phase contribution was the sum of three terms: The first term (DE intra ) was the difference in internal energy (bond, angle, and dihedral energy) of the OR and the odorant between the bound and unbound forms. The two other terms, (DE elec ) and (DE vdw ), were the non-bonded electrostatic and van der Waals energies between the bound and unbound form.
The solvation free energy includes two terms: The first terms DG pol was estimated using the Generalised Born model of Onufriev, Bashford, and Case. The value of the external dielectric constant was set to 80, while the internal one was set to 1. The non-polar term was considered as the solvent accessible surface area (SASA). The terms DG°M M 1 DG solv. represent the DG MM-GBSA component of the free energy of binding. The MM-GBSA analysis was performed on 40 snapshots sampled every 10 ns and covering the whole molecular dynamics simulation.

Results
To understand how mammalian ORs encode odorant concentrations, we conducted high-throughput in vivo screens via phosphorylated ribosomal protein S6 immunoprecipitation followed by RNA-Seq (pS6-IP-Seq). This characterized comprehensive OR responses and changes in OR activation patterns against ACT and TMT ( Fig. 2A) from the starting concentration of 0.01% (v/v) to undiluted 100% (v/v).

Modified pS6-IP-Seq enriched more ORs
We updated the pS6-IP-Seq method as originally described by Jiang et al. (2015), in an effort to profile ORs expressed in active OSNs that respond to lower concentrations of odorants. Figure 1A shows the updated pS6-IP-Seq pipeline. Similar to the Jiang et al. (2015) publication, we habituated and stimulated a mouse for 1 h with 10 ml of odorant at a given starting odorant concentration. We replaced the polyclonal antibody against pS6 with a monoclonal alternative to improve immunoprecipitation consistency. We updated annotations of OR genes to also include UTRs containing Bowtie aligned transcripts (Langmead et al., 2009;Li et al., 2009). As an example, Figure 1B shows coding region (CDS) annotation and updated annotation (CDS 1 UTR) for an ACT OR, Olfr923. Figure 1C also shows Olfr923's transcripts read depth aligned with STAR based on the new OR annotation (for details, see Materials and Methods; Dobin et al., 2013;Bray et al., 2016).
The updated OR annotations now include 1385 UTR annotated ORs, largely based on the work of Ibarra-Soria et al. (2014). To complete this set of OR definitions, we made three additional modifications: reannotation of 58 coding exonal ORs genes to include noncoding UTRs, correction of a misannotated OR (Olfr151) and addition of four previously unannotated ORs from the refGene database (O'Leary et al., 2016). More ORs were significantly enriched (FDR corrected p , 0.05) at 1% and 100% ACT and TMT with the UTRs included gene annotations, compared with ORs reported by Jiang et al. (2015), suggesting a more sensitive detection of OR enrichment with the updated pS6-IP-Seq method (Fig. 2B,D ;.

OR repertoires with different odorant concentrations
An increase in the number of ORs suggests that the updated pS6-IP-Seq method may allow us to identify activated ORs at lower odorant concentrations than what has been previously reported (Jiang et al., 2015). To test this, we stimulated mice with 0.01%, 0.1%, 1%, and 100% (v/v) of ACT and performed the updated pS6-IP-Seq at each condition (n = 3; for details, see Materials and Methods). At the lowest starting concentration, 0.01%, ACT enriched six ORs (FDR corrected p , 0.05; Fig. 2B). The number of ORs recruited increased from six to 115 (Fig. 2B) when the starting odorant concentration was increased from 0.01% to 100%. This shows that more ORs were activated at higher concentrations of odorants (Extended Data Fig. 2-11). Trace amine associated receptors (Taars), a small family of chemosensory receptors activated by volatile amines (Liberles and Buck, 2006) were not significantly enriched by the tested odorants.
To investigate the phylogenetic relationship of ORs responding to an odorant, we plotted enriched ORs on protein sequence-based phylogenetic trees (Fig. 2C). As the odorant concentrations increase, more diverse sets of ORs were enriched (Fig. 2C), consistent with findings previously reported (Jiang et al., 2015).
We performed the same experiment and analyses for another odorant, TMT. The number of ORs recruited increased from three to 68 for TMT (Fig. 2D) when the starting odorant concentration was increased from 0.01% to 100% TMT (Extended Data Fig. 3-3). Similar to ACT, more diverse sets of ORs are enriched as the odorant concentrations increase (Fig. 2E). These findings are consistent with earlier studies based on glomerular activation in the olfactory bulb (Rubin and Katz, 1999;Wachowiak and Cohen, 2001;Fried et al., 2002;Wilson et al., 2017) and complement these studies by identifying ORs that respond to a wide range of concentrations of odorant.
The distribution of OSNs expressing a given OR is confined within continuous and partially overlapping spatial zones along the dorsomedial (dorsal)-ventrolateral (ventral) axis on the olfactory epithelium (Ressler et al., 1993;Vassar et al., 1993;Miyamichi et al., 2005;Tan and Xie, 2018). To test whether active ORs for each odorant show bias in their zonal distributions, we conducted an RNA-Seq using dorsal and ventral olfactory mucosa (for details, see Materials and Methods) and found 412 and 558 ORs to be significantly enriched in dorsal and ventral zones, respectively (  (Table 1) is consistent with the zone index data reported by Tan and Xie (2018; p b , 0.001, Mann-Whitney test). Comparing between the distribution of responsive dorsal and ventral ORs of ACT and TMT at different concentrations, we found that ACT enriched ORs as a group are more dorsally distributed than TMT enriched ORs from low to high odorant concentrations (Fig. 2F

Most sensitive ORs may not be most robust responders at higher concentrations
At 0.01% starting concentration of ACT, Olfr923 was the most significantly enriched OR with the largest log 2 fold change (Log 2 FC) amongst significantly enriched ORs (FDR corrected p a = 1.56 Â 10 À13 , Log 2 FC = 2.96; Fig.  2B), suggesting it is among the most sensitive ORs that recognize ACT. Olfr923 remained significantly enriched from 0.1% to 100% ACT (FDR corrected p value p a , 0.001; Fig. 2B). However, the rank order of Olfr923 in terms of FDR corrected p values and Log 2 FCs went down at 1% and higher concentrations of ACT (Fig. 2B). Olfr983 is the second most significantly enriched OR with the second largest Log 2 FC amongst significantly enriched ORs at 0.01% ACT (FDR corrected p a = 2.55 Â 10 À9 , Log 2 FC = 1.88; Fig. 2B). Similarly, the rank order of Olfr983's activation also went down at higher starting concentrations of ACT (Fig. 2B). In contrast, Olfr907 and Olfr898, ORs with the largest and second largest Log 2 FCs among significantly enriched ORs at 100% ACT, were not significantly enriched at 1% and lower starting concentrations of ACT (Fig. 2B). Olfr30, the most significantly enriched OR at 100% ACT, was also not significantly enriched at 1% and lower starting concentrations of ACT. In addition, Olfr376, the most significantly enriched OR at 1% ACT, was also not significantly enriched at 0.1% and lower concentrations of ACT (Extended Data Fig. 3-3). These results suggest that sensitive ORs responding at lower odorant concentrations may not be the most robust responders at higher odorant concentrations, and vice versa.
We observed similar trends with TMT. The rank order of three significantly enriched ORs at 0.01% TMT (Olfr1395, Olfr165, and Olfr1297) also went down in terms of FDR corrected p values and Log 2 FCs at 1% and 100% starting concentrations of TMT (Fig. 2D). In contrast, Olfr837 and Olfr506, ORs with the largest and the second largest Log 2 FCs among ORs significantly enriched ORs at 100% TMT, were not significantly enriched at 0.1% and 0.01% TMT (Fig. 2D). Lastly, Olfr531 and Olfr376, the most and the second most significantly enriched OR at 100% TMT, were not significantly enriched at 0.01% ( Fig. 2D; Extended Data Fig. 3-4).
To validate pS6 induction for OSNs in a set of significantly enriched ORs, we conducted in situ hybridization and pS6 immunostaining for ORs responsive at 0.01% starting odorant concentrations. Consistent with our pS6-IP-Seq data, both Olfr923 and Olfr983, but not Olfr376expressing OSNs (serving as a control) showed significantly higher pS6 signals at 0.01% ACT compared with no odor controls (Fig. 3A,C e ; Extended Data Fig. 3-1) ( Table 1). Violin plots for Olfr1427 and Olfr901 (pS6-IP-Seq enriched at 0.01% ACT) show significantly higher pS6 signals in staining at 1% and higher starting ACT concentrations (Extended Data Fig. 3-2 f ). In situ hybridization and pS6 immunostaining also confirmed the activation of Olfr1395, Olfr1297, and Olfr165-expressing OSNs towards TMT (Fig. 3B,D g ; Extended Data Fig. 3-1) (Table  1). Consistent with our pS6-IP-Seq data, Olfr376-expressing OSNs did not show significantly higher pS6 signals at 0.01% starting concentration of TMT (Fig. 3D g ; Extended Data Fig. 3-1) (Table 1). As a negative control, we also stained for the pS6 intensity of Olfr1395 when mice were exposed to ACT instead of TMT. Consistent with our pS6-IP-Seq's differential expression data for ACT, Olfr1395 did not show significantly higher pS6 signals at any tested ACT concentrations (Extended Data Fig. 3-2 f ) (Table 1). Together, our data verified a set of sensitive ORs for ACT and TMT.

Correlation between pS6-IP-Seq and in vitro or responses
The in vitro responses of a panel of 500 ORs against ACT based on cAMP-mediated luciferase reporter gene assays in heterologous cells were reported to be correlated with the enrichment of ORs from pS6-IP-Seq (Jiang et al., 2015). Here we compared our new in vivo ACT pS6-IP-Seq dataset to the in vitro ACT dataset reported by Jiang et al. (2015). We first investigated whether significantly enriched ORs grouped using in vivo pS6-IP-Seq (orange circles) showed higher in vitro activation at three different concentrations of ACT (3 mM, 30 mM, and 300 mM), compared with the "not in vivo enriched" OR group (Fig. 4A, gray circles). Additionally, we created a graph for each ACT concentration, with FDR corrected p values on the x-axis, and in vitro responses on the y-axis to provide a sense how FDR corrected p values correlates with in vitro responses (Fig. 4B). Overall, 89 out of the 138 in vivo enriched ORs were tested in vitro and 46 (52%) were activated (Extended Data Fig. 4-1). Consistent with the results of Jiang et al. (2015), we found that the in vivo enriched OR group has a significantly higher in vitro activation from 0.1% to 100% ACT, but there are instances where in vivo enriched ORs do not show in vitro responses, and vice versa. P P P P P P P P P P P P Figure 4. In vivo OR enrichment correlates with the in vitro activation data. A, In vitro Hana3A cell lines cAMP-mediated reporter data for 516 out of the 1088 pS6-IP-Seq mapped intact ORs are classified into in vivo pS6-IP-Seq ACT enriched (orange, FDR corrected p , 0.05 and log 2 fold change . 0) and not enriched groups (gray, FDR corrected p . 0.05 or log 2 fold change , 0). In vitro normalized fold increase in luciferase signals were plotted for each OR groups across four ACT concentrations; 100% is based on the fold of increase of Olfr1126 activated by 300 mM ACT; 0% is determined by the fold of increase of empty rho-pCI vector stimulated with 3 mM ACT. Similarly, we observed that the ORs grouped based on in vitro luciferase assay activation (orange circles) had a significantly higher pS6-IP-Seq fold enrichment at all the ACT concentrations tested, compared with the "not in vitro activated" OR group (Fig. 4B h , gray circles) (Table 1); 46 out of the 105 in vitro activated ORs were also enriched in vivo based on pS6-IP-Seq (Extended Data Fig. 4-1).

Evaluating pS6-IP-Seq versus DREAM techniques
The DREAM technique is another method for the in vivo identification of odor activated ORs. DREAM relies on the transcriptional downregulation of individual ORs after ;5 h of odor exposure (von der Weid et al., 2015). Using DREAM, von der Weid et al. (2015) identified a set of ORs, 22 of which were confirmed by qPCR, to be responsive at 5% ACT.
To understand the relationship between DREAM and our updated pS6-IP-Seq, we focused on the top 45 ORs whose mRNA expressions were downregulated by 5% ACT (von der Weid et al., 2015; Extended Data Fig. 5-2). We remapped the sequence data using our updated OR definitions and performed differential expression analysis using the same bioinformatics pipeline as our pS6-IP-Seq datasets (Extended Data Fig. 5-3; for details, see Materials and Methods). Consistent with the von der Weid et al. (2015) study, all the top 45 ORs were downregulated. Comparing the DREAM data with our pS6-IP-Seq data, many of the top 45 ORs showed positive Log 2 FC values in pS6-IP-Seq ( Fig. 5; Extended Data Fig. 5-1). Overall, 19 of the top 45 ORs showed significant enrichment (i.e., FDR corrected p , 0.05) at least in one of the tested concentrations in our pS6-IP-Seq data. Log 2 FC values were negatively correlated between DREAM at 5% ACT and pS6-IP-Seq from 0.01% to 1% ACT (p values , 0.05 and R 2 range from 0.11 0.25; Fig. 5A i -C i ) (Table 1), while Log 2 FC values were weakly correlated between DREAM at 5% ACT and pS6-IP-Seq at 100% ACT (p = 0.2356 and R 2 = 0.03; Fig.  5D i ) (Table 1). Altogether, our analysis suggests that the degrees of pS6 induction and downregulation of mRNA expression are moderately correlated.

Validating Olfr923 with a gene knock-in approach
As an independent approach to demonstrate in vivo activation of Olfr923 expressing OSNs by ACT, we generated IRES-Cre gene knock-in mice at the Olfr923 locus and crossed them with Rosa26-lox-stop-lox-tdTomato reporter mice ( Fig. 6A; Extended Data Fig. 6-1) to label the Olfr923 glomerulus (Madisen et al., 2010).
To ask whether the Olfr923 glomeruli in the olfactory bulb were activated, we conducted c-Fos immunostaining in periglomerular cells induced by odor stimulation (Guthrie et al., 1993). Mice exposed to 0.01% starting concentration of ACT displayed a significant increase in c-Fos induction in periglomerular cells surrounding the Olfr923 glomeruli compared with no odor controls both in terms of cell counts and percentages of c-Fos-positive cells (p j , 0.05 paired t test, one-tailed; Fig. 6B,C; Extended Data Fig. 6-2), suggesting ACT activates Olfr923 OSNs and also the Olfr923 glomeruli.
The free energy of binding shows the high affinity of Olfr923 for ACT Our in vivo and in situ data consistently showed ACT to be an activator for Olfr923. Lastly, we asked if Olfr923 has evidence of structural properties that are favorable for ACT binding. Heptanal is a non-agonist of Olfr923 as previously reported (Ibarra-Soria et al., 2017). We first aimed to confirm Olfr923 activation in vitro using the GloSensor assay system to monitor the real time response of Olfr923 against ACT and heptanal. ACT or heptanal was dissolved into medium and odor-mediated cAMP induction was monitored in real time. Increases of over time were observed when Olfr923 was exposed to ACT, an indication of Olfr923 activation by ACT, while no significant luminescence increase was observed when Olfr923 was exposed to heptanal (Extended Data Fig. 7-1A). To quantify Olfr923 activations, we analyzed the area under the curve (AUC) of the normalized luminescence and generated dose response curves, confirming differential activation of Olfr923 by ACT and heptanal (Extended Data Fig.  7-1B).
A three-dimensional model of Olfr923 was built by homology modeling (for details, see Materials and Methods). ACT was docked into the canonical binding cavity of Olfr923, and heptanal serves as a negative control. Heptanal is a non-agonist of Olfr923 as previously reported (Ibarra-Soria et al., 2017). This complex was embedded in a membrane model, solvated in water and submitted to multiple molecular dynamics simulations. After several steps of minimization, equilibration and 800 ns of production with no constraints, ACT remained closely associated with the canonical binding cavity (Fig. 7A). The odorant binding cavity is more constricted around ACT, as shown for other non-olfactory GPCR structure bound to agonist (Wingler et al., 2019). However, heptanal exited the binding cavity to enter the lipid bilayer that composes the membrane, suggesting a low ligand-OR binding cavity affinity (Fig. 7B). Next, we computed the free energy of binding (DG binding ) by MMGBSA calculation ( 7C). This method was previously shown to be effective in discriminating OR agonists from non-agonists (Topin et al., 2014). Olfr923 bound to ACT has a favorable negative DG binding (Fig. 7C). In contrast, heptanal has an unfavorable positive DG binding (p k = 0.0371 paired t test, onetailed; Fig. 7C).
In order to identify the residues of Olfr923 involved in the stabilization of ACT and destabilization of heptanal, we decomposed the enthalpy of binding (DG MM-GBSA ) by the contribution of each residue (Fig. 7D). We saw that ACT binding is more strongly stabilized than heptanal with lower enthalpy values for its interaction with residues of the Olfr923 binding cavity. However, residues H165 4.55 , V205 5.39 , and M252 6.42 (Ballesteros-Weinstein numbering is shown in superscript; Ballesteros and Weinstein, 1995) stabilized heptanal binding more than for ACT in our molecular dynamics simulations (Fig. 7D). To identify the positions of the residues involved in the binding of both molecules, we projected them on our three-dimensional model (Fig. 7A,B). ACT engaged in interactions exclusively with residues belonging to the canonical odorant binding cavity of Olfr923 (Man et al., 2004;Katada et al., 2005;de March et al., 2015a), positioning ACT in an effective binding position likely to trigger the receptor activation (Fig. 7A). Heptanal was more mobile into the binding cavity. Interestingly, the residues H165 4.55 , V205 5.39 , and M252 6.42 which participated only in heptanal stabilization are always located outside of the canonical binding cavity (Fig. 7B,D), either in the extracellular part of the receptor or between TM4 and TM5 under the canonical binding cavity.
This explains the tendency of heptanal to exit the canonical binding cavity during our molecular dynamic simulations. Furthermore, it was shown recently that molecular dynamic simulations of ORs can sample active or inactive states of the receptor when bound to an agonist or a non-agonist, respectively (de March et al., 2015b. The most accurate parameter to monitor receptor activation in molecular dynamics simulations is the distance between the intracellular part of TM3 and TM6, particularly between the two protagonists of the ionic lock in GPCR ( Fig. 7E; Altenbach et al., 2008;de March et al., 2015bde March et al., , 2018. The inactive conformation of an OR is maintained by the ionic interaction, called ionic lock, between D 3.49 belonging to the motif D 3.49 RY common to the GPCR superfamily and R 6.30 belonging to the motif R 6.30 xKAFSTCASH specific to the OR family. During activation, conformational changes in the receptor are associated with a break of this ionic interaction and a switch of the intracellular part of TM6 toward the membrane, opening an intracellular cavity for the G-protein coupling (de March et al., 2015b. By observing our multiple molecular dynamics simulations, we noticed that Olfr923 structures bound to ACT show mostly opened ionic locks (.12 Å), suggesting active conformation (Fig. 7E,F). However, Olfr923 structures bound to heptanal exhibit a majority of receptor conformations possessing a close ionic lock (,12 Å), suggesting that the Olfr923 is mainly in its inactive conformation (Fig. 7E,F). This difference in the distribution of receptor conformations between ACT and heptanal shows that our model captured the ligand binding features required to trigger Olfr923 activation. Altogether, our data suggests that the Olfr923 binding cavity has favorable properties for ACT effective binding, and further supports our findings that ACT binds to and activates Olfr923.
Our results from in vivo, in situ, and in silico approaches, together with the reported activation of Olr923 by ACT in vitro (Jiang et al., 2015), suggest that ACT binds and activates Olfr923.

Discussion
In this study, we investigate how odorant concentration information is represented at the receptor level. This comprehensive representation is reported as activated OR across a 10,000-fold starting concentration range. We identified some of the most sensitive ORs for ACT and TMT and validated the activation of Olfr923 by ACT at the glomerular level in the olfactory bulb.

Updated method to identify ORs activated by odors in vivo
To more accurately describe the concentration-dependent combinatorial OR code, we report an updated pS6-IP-Seq method to profile mRNA expression in odor activated OSNs. We used a monoclonal pS6 antibody, streamlined the library preparation methods (Knight et al., 2012;Jiang et al., 2015), and adopted STAR for efficient sequence reads mapping (Dobin et al., 2013). We can now identify more ORs at given odor concentrations and at lower and more physiologically relevant odor concentrations. The updated pS6-IP-Seq can be adopted to identify responsive ORs across a wide range of monomolecular odorants, odorant mixtures, and complex natural odors to further our understanding of peripheral odor coding (Isogai et al., 2018).
Our OR lists overlap with a majority of the ORs identified by the Jiang et al. (2015) publication. Out of the total 86 ACT responsive ORs reported by Jiang et al. (2015), 51 of these ORs were also identified by our current study. The vast majority of the non-overlapping ORs exhibit a positive Log 2 FC, indicating that these ORs were enriched but did not reach our statistical criteria (FDR corrected p , 0.05). It is also likely that the difference in antibodies against pS6 (i.e., polyclonal antiserum in Jiang et al. vs monoclonal antibodies in the current study) contribute some differences (Jiang et al., 2015). Changes in the bioinformatics pipelines also contribute to differences.
Our data support the notion that higher concentration of odorants activates more ORs (Fried et al., 2002;Mainland et al., 2014). In the vast majority of cases, ORs significantly enriched with odor stimulation at a low concentration are also enriched at higher concentrations. However, there are some instances where ORs show significant enrichment only at lower concentrations. The statistical criteria based on the FDR corrected p , 0.05 cutoff does not represent the critical boundary between responsive ORs and non-responsive ORs. It serves as a selection criterion to identify populations of responsive ORs activated backed by statistical confidence. It is also Figure 7. In silico investigation of ACT binding to olfr923. The three-dimensional structure of Olfr923 was built by homology modeling and bound to ACT (red, A) and heptanal (blue, B). Helices are represented in white ribbons and the odorant molecules are shown in Van der Walls volumes. The residues interacting with the odorant molecule during molecular dynamics simulations are highlighted by dotted clouds in the corresponding color (ACT-red, heptanal-blue). Residues specific to heptanal binding, H165 4.55 , V205 5.39 and M252 6.42 (Ballesteros-Weinstein numbering shown in superscript; Ballesteros and Weinstein, 1995), are represented in blue licorice. C, Computed enthalpy DG MM-GBSA, entropy TDS vib , and free energy of binding DG binding for heptanal and ACT bound to Olfr923; pp , 0.05 by paired t test (ACT compared with heptanal). N = 3 simulations per odorant. D, Decomposition per residue of the enthalpy DG MM-GBSA, for ACT (red) and heptanal (blue). E, Tridimensional structure of olfr923 in active and inactive conformations. The intracellular part of TM6 and the two residues involved in the ionic lock are highlighted for the ACT bound system (shade likely that the proportional nature of RNA-Seq contributes to apparent lower fold changes at the higher odorant concentration. With future advancements in next-generation sequencing data collection and analysis, we expect refinement in the list of enriched ORs.

OR populations encoding odor concentrations
Humans perceive the same odor with higher intensity and varying degree of pleasantness as odor concentration increases (Moskowitz et al., 1976;Anderson et al., 2003). In some cases, odor quality and/or pleasantness shift when different concentrations of the same odorant are presented (Gross-Isseroff and Lancet, 1988), although there is no such report for ACT and TMT. Higher concentrations of odorants that result in higher OSN spiking rates (Rospars et al., 2000;Bhandawat et al., 2005;Grosmaitre et al., 2006) are likely to result in higher S6 phosphorylation as shown in the retinal cells (Milner and Do, 2017). This cascade of events, could, in turn, contribute to the neural representation of elevated odor intensity. Higher concentration of odors activate additional, lower affinity ORs at the olfactory epithelium, resulting in more glomerular activations in the olfactory bulb (Fried et al., 2002;Khan et al., 2010). These additionally recruited, low-affinity ORs may drive differences in odor quality and perceived intensity by modulating neuronal firing in higher olfactory areas (Stettler and Axel, 2009). In this study, we identified a set of ORs responding to both low and high concentrations odorant and another set of ORs recruited only at high odorant concentrations. In the future, resources provided in this study will be useful in defining the roles of high-affinity and low-affinity ORs by targeting specific ORs to manipulate in reference to their in vivo activation profiles.

Specific ORs encoding odor identity
Another dimension of smell encoding, besides keeping track of odorant concentrations, is for neurons to also represent concentration differences of a given odor (Parabucki et al., 2019). Despite facing a dynamic environment where they encounter hundreds to thousands of volatile odors at various concentration gradients, humans and animals are capable of maintaining stable odor quality perception across concentrations (Krone et al., 2001;Wachowiak and Cohen, 2001;Uchida and Mainen, 2007;Mainland et al., 2014;Sirotin et al., 2015;Wilson et al., 2017).
What is not yet clear is how odor identities are encoded and maintained despite massive differences in the activated OR repertoires across concentrations. Stable odor ratio information (Uchida and Mainen, 2007), input-output transformation at the level of the olfactory bulb (Storace and Cohen, 2017) and early activated ORs during the first sniff (Wilson et al., 2017) are implicated in concentration invariance. Sensitivity of ORs, as well as abundance and positions of OSNs expressing individual ORs within the olfactory epithelium (zonal expression), can play a role in determining the timing of neuronal activation at the level of glomeruli in the olfactory bulb (Schneider et al., 1966;D'Hulst et al., 2016). Our datasets reporting ORs responding to low concentrations of odorants, in conjunction with a number of OSNs expressing different ORs which correlates with OR mRNA abundance (Ibarra-Soria et al., 2014) together with zonal expression of ORs (Ressler et al., 1993;Vassar et al., 1993;Miyamichi et al., 2005;Tan and Xie, 2018) can provide valuable information in targeting a specific set of ORs to investigate their role in odor coding.

Comparisons of pS6-IP-Seq, DREAM and heterologous expression in measuring or responses
The DREAM technique identifies odor responsive ORs based on a decrease in OR mRNA abundance after odor stimulation (von der Weid et al., 2015). Our analysis suggests that OR responses based on in vitro assays and the two in vivo-based assays (pS6-IP-Seq and DREAM) correlated each other yet they did not show strong linear correlations. For example, Olfr923 is one of the most sensitive and robust ORs responding to ACT based on pS6-IP-Seq, ranking top in terms of Log 2 FC amongst significantly enriched ORs from 0.01% to 1% ACT. Olfr923 ranked 36 th out of the top 45 responsive ORs based on 5% ACT DREAM RNA-Seq fold change in the von der Weid publication, Olfr923 now ranks fifth when our bioinformatics pipeline is used to reanalyze the 5% ACT DREAM RNA-Seq data (von der Weid et al., 2015;Fig. 5;Extended Data Fig. 5-3). When tested in vitro via the luciferase assay, Olfr923 ranked 32nd out of the 105 responsive ORs based on fold luciferase induction with 300 mM ACT (Jiang et al., 2015). Difficulties in functionally expressing certain ORs in heterologous cells partly explains the differences. Future mechanistic understanding of how odorant activation drives S6 phosphorylation or OR mRNA decrease should help explain these differences across the in vivo and in vitro OR response studies. In addition, the abundance of ORs can also be widely varied (Ibarra-Soria et al., 2014). Some ORs with low read counts did not reach the set significance (FDR corrected p . 0.05) despite large fold changes. A deeper sequencing read depth will help determine whether these ORs are bona fide ORs responding to the tested odorants.
In summary, our study comprehensively identifies ORs responding to varying concentrations of odorants and contributes to the understanding of how odorant concentration information is encoded at the receptor level in the peripheral olfactory system. Identifying ORs responding to both low to high odorant concentrations continued of red) and heptanal bound system (blue). Insert, Intracellular view of the receptor. The TM3-TM6 distance monitored in F is shown by a dotted line. F, Distribution of the number of frames in the molecular dynamics simulations sampling TM3-TM6 distances for ACT (red) and heptanal (blue). See Extended Data Figure 7-1 for GloSensor assay in vitro validation and the molecular dynamics workflow. will be valuable in future studies evaluating the contribution of high-affinity and low-affinity ORs towards odor perception.