RNA-mediated symmetry breaking enables singular olfactory receptor choice

Mice

Mice were treated in compliance with the rules and regulations of the Institutional Animal Care and Use Committee of Columbia University under protocol number AABG6553. Mice were euthanized using CO₂ followed by cervical dislocation. Both male and female mice were used for experiments. All experiments were performed on dissected olfactory epithelium tissue or on dissociated cells prepared from whole olfactory epithelium tissue. This study used several mouse lines (Mus musculus) on mixed C57BL/6J and 129 backgrounds. For Dip-C, H3K27ac HiChIP and liquid Hi-C, cells expressing the OR P2 were obtained by crossing tetO-P2-IRES-GFP mice to Gng8(gg8)-tTA mice⁴² and sorting GFP⁺ cells from dissociated MOE. For Dip-C, Gng8tTA>tetO-P2 and Mor28-IRES-GFP²⁵ mice were crossed to CAST/EiJ mice (Jax strain 000928) to generate F₁ hybrids where known single-nucleotide polymorphisms could be used for haplotype imputation. For the Hi-C data shown in Supplementary Fig. 8a–f, horizontal basal cell and INP analyses were performed on previously published Hi-C data⁴, iOSNs were isolated by performing Hi-C on heterozygous Atf5-IRES-RFP⁴³ OMP-IRES-GFP mice, sorting RFP⁺GFP⁻ cells, and GFP⁺ cells from OMP-IRES-GFP mice²⁵ were used to isolate mOSNs. For ATAC-seq, RNA-seq and in situ Hi-C results shown in Fig. 4g–n, Mor28-IRES-cre²⁵, Rosa26(LSL-tdTomato/+)⁴⁴, OMP-ires-tTA and tetO-P2 alleles were crossed to create mice heterozygous for all alleles. For immunofluorescence, Hi-C and RNA-seq, tetO-P2(nc) mice were generated by performing CRISPR/non-homologous end-joining on heterozygous tetO-P2 embryos with the following guide targeting the 5′ region of the P2 CDS (5′-GGGAAACTGGACAACTGTCA-3′). Verification of frameshift was done by performing TIDE analysis on PCR amplicons of the unmutated and mutated tetO-P2 sequence from gDNA of F₁ pups of founder mice and stock tetO-P2 mouse lines. For immunofluorescence and RNA-seq, tetOM71(nc) mice were generated by first assembling a tetOM71(nc)-IRES-GFP construct made by performing an NEB HiFi assembly using an M71(nc)-IRES-GFP gene block made with Integrated DNA Technologies (IDT,  and a pTRE Tight tetO-Fv2E-Perk plasmid (gift from H. Shayya). The M71 CDS was rendered non-coding by changing the 11th amino acid to a stop codon and mutating all in-frame methionine codons to another missense codon that would result in few modifications to RNA secondary structure, thereby preventing any in-frame translation. NheI restriction digest released a fragment containing the tetOM71(nc) construct, which was used for pronuclear injection in B6CBAF1 zygotes. Tail biopsy and PCR were used to identify founder mice containing the transgene; these were crossed to Omp-irestTA⁴⁵ animals to screen for both germline transmission and tTA-dependent transgene expression in mOSNs. tetOM71-LacZ mice⁴⁶ and tetO-GFP mice were also crossed to OMP-tTA and/or Gng8tTA drivers for immunofluorescence and RNA-seq experiments. For all experiments, mice were between 5 and 12 weeks of age.

Fluorescence-activated cell sorting

Cells were prepared for FAC sorting as previously described⁴ by dissociating olfactory epithelium tissue with papain for 40 min at 37 °C according to the Worthington Papain Dissociation System. Cells were washed twice with cold PBS before being passed through a 40 μm strainer. Live (DAPI-negative) fluorescent cells were collected for RNA-seq and liquid Hi-C. Alternatively, for Hi-C and HiChIP, cells were fixed for 10 min in 1% formaldehyde in PBS at room temperature, quenched with glycine and washed with cold PBS before sorting of fluorescent cells. For Dip-C, cells were fixed in 2% formaldehyde in PBS at room temperature for 10 min, inactivated with 1% bovine serum albumin (BSA) and washed with cold 1% BSA in PBS before sorting of fluorescent cells. All cells were sorted on a Beckman Coulter Low Flow Astrios EQ.

Olfactory epithelium immunofluorescence

Immunofluorescence assays were performed as previously described⁴³. In brief, dissected MOEs were fixed in 4% (w/v) paraformaldehyde in PBS for 1 h at 4 °C and then washed three times for 10 min each time in PBS. Olfactory epithelia were decalcified overnight at 4 °C in 0.5 M EDTA (pH 8) and washed again in PBS. MOEs were cryoprotected overnight at 4 °C in 30% (w/v) sucrose in PBS, embedded in OCT, frozen over an ethanol/dry ice slurry and stored at −80 °C until sectioning. To ensure full coverage of the MOE, tissue was serially sectioned in the coronal plane, moving from the flat posterior surface to the anterior surface. Six slides were prepared with four sections per slide, of 15 mm sections collected on slides starting at the moment when turbinate 3 separated from the dorsalmost aspect of the epithelium⁴⁷. Slides were frozen at −80 °C until the day of staining experiments, when they were thawed, washed for 5 min in PBS and postfixed for 10 min at room temperature in 4% (v/v) formaldehyde (Thermo Fisher) in PBS. Tissue was then washed three times (5 min each time, in PBS + 0.1% Triton X-100 (Sigma)) and blocked for 1 h at room temperature in 4% (v/v) donkey serum (Sigma) + 1% Triton X-100 in PBS. Primary antibodies against GFP (chicken anti-GFP ab13970, 1:2,000), P2 (Olfr17 antibody were raised in guinea pig, 1:2,000), M71 (1:3,000)¹¹ and/or LacZ (abcam ab4761, 1:16,000) were diluted in block solution and used for incubation overnight at 4 °C. The following day, sections were washed, incubated with secondary antibodies (Jackson Immunoresearch, 1:500 in block solution) for 1 h at room temperature, washed again and mounted using VECTASHIELD Vibrance (Vector Labs) mounting medium. Images were rendered with ImageJ 2.0.0.

In situ Hi-C, liquid Hi-C and H3K27ac HiChIP

In situ Hi-C and liquid Hi-C

In situ Hi-C was performed exactly as previously described⁴. The liquid Hi-C protocol²⁶ was integrated into our Hi-C protocol to perform liquid Hi-C in OSNs. In brief, MOE was dissociated from gg8-tTA>tetO-P2 mice, and 400,000 GFP⁺ cells were sorted as described above per condition per replicate, with three biological replicates per time point. After sorting, cells were pelleted at 600g, for 10 min at 4 °C, and resuspended in 300 μl chilled lysis buffer (50 mM Tris pH 7.5, 0.1% Igepal, 150 mM NaCl, protease inhibitor in water). Samples were then pelleted for 7 min at 700g and 4 °C and then resuspended in 105 μl DpnII-MasterMix (DpnII Buffer, 250 U DpnII) and placed on a preheated thermomixer at 37 °C with shaking at 900 rpm for 5 min, 30 min or 60 min. Samples were immediately placed on ice for 10 min after predigestion. For 0 min liquid Hi-C, after lysis, cells were immediately processed for fixation. For fixation, samples were diluted into 1% formaldehyde in PBS, rotated on a rotisserie for 10 min at room temperature and quenched with 1/10 volume of 1.25 M glycine. Samples were pelleted at 2,500g, for 5 min at 4 °C, washed with PBS and then resuspended in nuclear permeabilization solution (as described in the in situ Hi-C protocol). All subsequent steps and the library preparation were performed as previously described⁴. Samples were sequenced paired-end 50 bp or 100 bp on Illumina NextSeq 550, Illumina NovaSeq2000 or Illumina NextSeq2000. Three biological replicates were created for all liquid Hi-C experiments; once libraries had been confirmed to be similar, they were merged. Heatmaps were generated from merged cooler files, and Welch’s two-sample t-tests on CSS scores were performed on unmerged replicates.

H3K27ac HiChIP

The HiChIP protocol was given by the Chang laboratory and integrated into our Hi-C protocol for H3K27ac HiChIP on OSNs²⁸. MOE from 5–7 gg8-tTA>tetO-P2 mice were dissociated to obtain 4 million GFP⁺ cells per replicate, for a total of two replicates. Cells were processed according to the in situ Hi-C protocol with the following exceptions: nuclei were digested for only 2 h instead of overnight, and complete nuclei digestion was verified by running reverse cross-linked digested nuclei on a DNA agarose gel. After ligation, nuclei were pelleted at 2,500g, for 5 min at 4 °C, and stored overnight at −20 °C. The next day, nuclei were resuspended in 130 μl of HiChIP nuclear lysis buffer (50 mM Tris pH 7.5, 10 mM EDTA, 1% sodium dodecyl sulfate, protease inhibitor in water) and sheared on a Covaris S220 with the following parameters: duty cycle, 2%; PIP, 140; cycles/burst, 200; time, 4 min. After shearing, samples were precleared, immunoprecipitation was performed with 1 µg H3K27ac antibody per 4 million cell input (Abcam GR323193701) and libraries were prepared exactly as previously described²⁸. Samples were sequenced paired-end 50 bp on an Illumina NextSeq2000.

In situ Hi-C, liquid Hi-C and HiChIP alignment and data preprocessing

Alignment and data preprocessing were performed exactly as previously described²². In brief, reads were aligned to the mm10 genome using the distiller pipeline ( requirements: java8, nextflow and Docker); uniquely mapped reads (mapq > 30) were retained, and duplicate reads were discarded. Contacts were then binned into matrices using cooler⁴⁸. Data pooled from two to three biological replicates were analysed, after the results of analyses of individual replicates had been confirmed to be similar.

RNA-seq

RNA extraction and library preparation

All RNA-seq experiments were performed under RNA clean conditions. For RNA-seq, live cells were sorted into RNase-free PBS, pelleted at 600g, for 5 min at 4 °C, then resuspended in 500 μl TRIzol, flash-frozen in liquid nitrogen and stored overnight at −80 °C. RNA extraction was performed the next day. TRIzol suspensions were thawed on ice, 1/5 V of 1-bromo-3-chloropropane was added, and tubes were shaken vigorously to combine phases. Phases were allowed to separate for 2 min at room temperature, then tubes were centrifuged at 10,500 rpm, for 15 min at 4 °C, in an Eppendorf centrifuge C5424R. We collected the upper aqueous phase and transferred to a new tube. Then, 1/2 V of isopropanol and 1 μl of linear polyacrylamide (Sigma Aldrich 56575) were added, the tube was inverted to mix the contents, and RNA was allowed to precipitate for 10 min at room temperature. Tubes were centrifuged for 10 min at 10,500 rpm and 4 °C. The supernatant was removed, and 1 V of 75% ethanol was added to the pellet, which was dislodged by flicking the tube. Tubes were centrifuged for another 5 min, at 10,500 rpm and 4 °C. Ethanol was removed, and tubes were allowed to air dry for 5 min until the pellet turned clear. Next, we added 26 μl of RNase-free water, 3 μl of Ambion DNase I 10× buffer and 1 μl of DNase I (AM2222) to remove all DNA and incubated tubes at 37 °C for 30 min. RNA was purified by a 1.5× AMPure bead clean-up, measured on a nanodrop and used as the input for library preparation with a SMARTER Stranded Total RNA-Seq Kit – Pico Input Mammalian v2 (TaKaRa Bio USA). OMP-tTA>tetO-GFP, gg8-tTA>tetO-GFP and two gg8-tTA>tetO-P2 libraries were prepared with the TruSeq kit. However, mOSN samples were compared with both OMP-tTA>tetO-GFP (TruSeq prep) and OMP-IRES-GFP (TaKaRa Bio USA), which label the same neurons, and produced the same results (Extended Data Fig. 10c,g–i). Libraries were sequenced on either a NextSeq550 or a NextSeq2000 and were sequenced to a targeted coverage of approximately 25 million reads. All RNA-seq experiments were performed with two to three biological replicates.

RNA-seq data processing and analysis

Data processing and analysis was performed as previously described¹². In brief, adaptor sequences were removed from raw sequencing data with CutAdapt. RNA-seq reads were aligned to the mouse genome (mm10) using STAR⁴⁹. SAMtools was used to select uniquely aligning reads by removing reads with alignment quality alignments below 30 (-q 30). RNA-seq data were analysed in R with the DESeq2 package⁵⁰. For MA plots, DESeq2 normalized gene counts were compared between control and knockout mice, and significantly changed genes were identified with an adjusted P value cutoff of 0.05. DESeq2 normalized counts were used to examine expression levels of genes (Extended Data Fig. 2a–c). Principal component analysis on all genes except Olfr genes was performed on RNA-seq datasets, to separate cells according to their developmental cell stage (Extended Data Fig. 10b).

ATAC-seq

ATAC-seq library preparation

ATAC-seq libraries, data processing and bigwig generation were performed exactly as previously described¹². In brief, cells were pelleted (500g, 5 min, 4 °C) and then resuspended in lysis buffer (10 mM Tris-HCl, pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% IGEPAL CA-630). Nuclei were immediately pelleted (1,000g, 10 min, 4 °C). Pelleted nuclei were resuspended in transposition reaction mix prepared from Illumina Nextera reagents (for 50 μl: 22.5 μl water, 25 μl 2× TD buffer, 2.5 μl Tn5 transposase). The volume of the Tn5 transposition reaction was scaled to the number of cells collected: 1 μl mix per 1,000 cells. If fewer than 10,000 cells were collected by FACS, 10-μl-scale reactions were performed. Transposed DNA was column purified using a Qiagen MinElute PCR cleanup kit (Qiagen). The transposed DNA was then amplified using barcoded primers and NEBNext High Fidelity 2× PCR Master Mix (NEB). Amplified libraries were purified using Ampure XP beads (Beckman Coulter) at a ratio of 1.6 μl of beads per 1 μl of library and eluted in 30 μl of elution buffer (10 mM Tris-HCl pH 8, 0.1 mM EDTA). Libraries were sequenced on either a NextSeq550 or a NextSeq2000 and were sequenced to a targeted coverage of approximately 25 million reads.

ATAC-seq data processing

Adaptor sequences were removed from raw sequencing data with CutAdapt, and reads were aligned to the mouse genome (mm10) using Bowtie2. Default settings were used, except that a maximum insert size of 1,000 (-X 1,000) was allowed for ATAC-seq. PCR duplicate reads were identified with Picard and removed with SAMtools. SAMtools was used to select uniquely aligning reads by removing reads with alignment quality alignments below 30 (-q 30). For ATAC-seq, regions of open chromatin were identified by running HOMER peak calling in ‘region’ mode, with a fragment size of 150 bp and a peak size of 300 bp. For ATAC-seq signal tracks, the results of replicate experiments were merged, and HOMER was used to generate 1 bp resolution signal tracks normalized to a library size of 10,000,000 reads. Reads were shifted 4 bp upstream to more accurately map the Tn5 insertion site. Reads were extended to the full fragment length, as determined by paired-end sequencing. Bigwigs were visualized with the Integrated Genome Browser 9.0.0.

Dip-C generation

Dip-C and data preprocessing

Cas mice were crossed to gg8-tTA>tetO-P2-IRES-GFP or Mor28-IRES-GFP heterozygous F₁ hybrids. Dip-C and data preprocessing were performed exactly as previously described²² and following the quality control metrics as previously described¹³, with the following exceptions. Each Dip-C library was sequenced on a single lane of an Illumina NovaSeq 6000. Reads were trimmed with CutAdapt v.1.17, and Dip-C libraries were aligned with BWA 0.7.17. Haplotype-imputed single-cell contacts were generated using the dip-c package ( requirements: hickit r291 and k8-Linux K8: 0.2.5-r80. We excluded cells that had fewer than around 400,000 contacts, a low contact-to-read ratio, or high variability in three-dimensional structure across computational replicates. Overall, the median number of contacts across nuclei was 715,690 contacts per cell for 74 cells for Mor28-IRES-GFP Dip-C and 694,462 contacts per cell for 84 cells for gg8-tTA>tetO-P2-IRES-GFP Dip-C, for a total of 161 cells. Three-dimensional reconstruction of Dip-C models was performed in PyMOL 2.5.3 as previously described²¹.

DNA FISH

Oligopaint probes specific for 20 kb encompassing the 30 most interacting GIs (based on bulk Hi-C results) and for the P2 locus were generated using oligominer scripts ( Sections of the MOE were fixed, denatured and hybridized as previously described^51,52. Imaging was performed using the Vutara VXL at the Zuckerman Institute Imaging Platform.

Multiome generation

Purification of nuclei

Nuclei must be purified under RNA clean conditions. A cell suspension of mouse MOE was obtained from an adult mouse following the dissociation conditions previously described¹². Cell pellets were immediately resuspended in 300 μl of cold RNAse-free lysis buffer (10 mM Tris-HCl, pH 7.4, 10 mM NaCl, 3 mM MgCl₂, 0.1% IGEPAL CA-630), and nuclei were pelleted in an Eppendorf 5810R centrifuge at 1,000g for 10 min at 4 °C. Nuclei were resuspended in 500 μl 10× homogenization buffer (100 mM Trizma base, 800 mM KCl, 100 mM EDTA, 10 mM spermidine trihydrochloride, 10 mM spermidine tetrahydrochloride in double-distilled H₂O), and the pH was adjusted to 9–9.4 with NaOH. Instructions for preparation of homogenization buffer can be found in Zhang et al.⁵³. RNAse inhibitor (NEB MO314L) was added, followed by 500 μl 82% OptiPrep solution (4.1 ml OptiPrep solution (Sigma Aldrich D1556-250ML), 25 μl 1 M CaCl₂, 15 μl 1 M magnesium acetate, 50 μl 1 M Tris pH 8, 810 μl water), and the mixture was placed on ice. Then, 1 ml homogenate was carefully added on to 1 ml of 48% OptiPrep solution (2.4 ml OptiPrep solution, 800 μl 1 M sucrose, 25 μl 1 M CaCl₂, 15 μl 1 M magnesium acetate, 50 μl 1 M Tris pH 8, 1,710 μl water) and spun down in a precooled swinging bucket centrifuge (Eppendorf 5810R) at 32,00g for 20 min at 4 °C, with acceleration 5/9 and deceleration 0/9 (no break)⁵⁴. The supernatant was aspirated and disposed of without dislodging the pellet. The pellet was air-dried and resuspended in 500 µl PBS diluted with 0.04% BSA with RNAse inhibitor. Cell concentration was measured for accurate loading into the 10× pipeline. Two independent multiomes were generated from a 12 week old (Fig. 1, wild-type background) and a 5-week-old mouse (Extended Data Fig. 1; gg8-tTA>tetO-P2(nc) background) and analysed separately. Both multiomes produced the same findings.

10x Genomics scATAC and scRNA library generation

Joint scRNA-seq and scATAC-seq libraries were prepared in collaboration with the Columbia Genome Center using the 10x Genomics Single Cell Multiome ATAC + Gene Expression kit according to the manufacturer’s instructions. Both 10X Single-Cell Expression (GEX) and ATAC libraries were sequenced to around 350 million reads on an Illumina NovaSeq 6000 150PE.

Generation of aligned multiome data

Raw sequencing data were demultiplexed with cellranger-arc mkfastq and aligned with cellranger-arc count. An mm10 fasta file and a custom GTF with extended OR annotations⁵⁵ were used to generate a reference package for alignment with cellranger-arc mkref. Our multiome contained an estimated 8,856 cells (12,936 cells for independent replicates; Extended Data Fig. 1) from the MOE, with a median of 2,671 high-quality ATAC fragments per cell (median 9,078 high-quality ATAC fragments per cell for independent replicates; Extended Data Fig. 1) and a median of 1,316 GEX genes per cell (1,006 GEX genes per cell for independent replicates; Extended Data Fig. 1). All multiome data were analysed in R v.4.1.3 using packages Signac v.1.6.0 and Seurat v.4.1.0.

Molecular dynamics simulations of GI hubs in OSNs

To investigate the symmetry-breaking mechanism of GI hubs occurring in OSNs, classical molecular dynamics simulations were used⁵⁶. Each hub was made of three distinct polymers, modelled as standard self-avoiding-walk strings composed of N = 30 beads. Each polymer was equipped with three binding sites, located in the central region. Polymer ends in a specific hub were anchored to the vertices of a hexagon (Fig. 5c) to ensure hub specificity and spatial separation between the polymers in the hub. Other geometries (for instance, triangular) gave similar results. Binding sites could attractively interact with binders with an affinity E_P and binder total concentration c. In addition, binders could interact among themselves with affinity E_B. For the sake of simplicity, polymer bead and binders had the same diameter σ and mass m, which were both set to 1 (dimensionless units)⁵⁶. All particles interacted with a repulsive Lennard–Jones (LJ) potential to take into account their excluded volume, with diameter σ and energy scale ε = 1k_BT, where T is the temperature and k_B is the Boltzmann constant. Between two consecutive beads of a polymer, a finite extensible nonlinear elastic⁵⁶ potential was used, with length constant R₀ = 1.6σ and elastic constant K = 30k_BT/σ², as previously described⁵⁷.

The interactions among binders, as well as the interactions between binders and binding sites, were modelled as a truncated, shifted LJ potential⁵⁷: \({V}_{{\rm{LJ}}}(r)=4\varepsilon \,\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^{6}-{\left(\frac{\sigma }{{R}_{{\rm{int}}}}\right)}^{12}+{\left(\frac{\sigma }{{R}_{{\rm{int}}}}\right)}^{6}\right]\) for R_int < 1.3σ and 0 otherwise, where r is the distance between particle centres, and ε, sampled in the range 8–12 k_BT, regulates the interaction intensity. The affinities E_B shown in Fig. 5c,d correspond to the minimum of V_LJ. For the sake of simplicity, the interaction between binder and binding sites was kept constant (E_P = 3.5k_BT). To map the length scale σ in physical units, we equalized the average interhub distance of nearest neighbouring hubs with the median interhub distance of ∼2 μm; this was estimated by measuring the average inter-GI distance in Dip-C nuclei, which was 33.4 p.r., obtaining σ = 60 nm. Binder concentrations were computed as previously described⁵⁷, using c = N_B/VN_A, where N_B is the number of binders, V is the volume (in litres) of the simulation box and N_A is the Avogadro number.

The system was in contact with a thermal bath at temperature T; therefore, positions evolved according to the Langevin equation⁵⁸, with the following standard parameters: friction coefficient ζ = 0.5, temperature T = 1 and timestep dt = 0.012 (ref. ⁵⁷). Integration was performed with a velocity Verlet algorithm using the LAMMPS software⁵⁹. The simulation was performed in a cubic box (linear size D = 64σ) with boundary periodic conditions to avoid finite size effects. For each parameter setting, we performed ten independent simulations. The system was initialized with polymers in random self-avoiding-walk states and binders randomly located in the simulation box and then equilibrated up to 10⁸ time = steps. Configurations were logarithmically sampled up to the equilibrium sampling frequency, that is, every 10⁵ timesteps.

Phase diagram and symmetry-breaking dynamics

The phase diagram was obtained by considering several different combinations of system control parameters, that is, binder self-interaction affinity E_B and binder concentration c. Symmetry-breaking events were called if, at equilibrium, a large and stable aggregate of binders in a GI hub was detected. To this end, we performed standard hierarchical clustering applied directly to the coordinates of binders, using their Euclidean distance as a metric⁶⁰. Clustering was performed using the linkage function from the Python package scipy.cluster. Then, a distance threshold R_thr = 1.3σ (as large as the attractive LJ distance cutoff) was set, and a cluster was defined as the set of binders whose cophenetic distance was lower than R_thr.

To study the dynamics of symmetry-breaking events associated with the formation of a stable cluster in a single GI hub, we considered system configurations from the starting state to the equilibrium state. For each sampled timestep, we applied the clustering procedure described above and then selected the largest clusters, that is, those containing the highest fractions of binders. We then used averaging over independent runs to generate the curves shown in Fig. 5d.

Statistics

All statistical analyses used Welch’s two-sample t-test. All averages are reported as mean ± s.e.m. In plots with error bars, points are centred on the mean, and error bars indicate the s.e.m.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.