Coadaptation of the chemosensory system with voluntary exercise behavior in mice.

Ethologically relevant chemical senses and behavioral habits are likely to coadapt in response to selection. As olfaction is involved in intrinsically motivated behaviors in mice, we hypothesized that selective breeding for a voluntary behavior would enable us to identify novel roles of the chemosensory system. Voluntary wheel running (VWR) is an intrinsically motivated and naturally rewarding behavior, and even wild mice run on a wheel placed in nature. We have established 4 independent, artificially evolved mouse lines by selectively breeding individuals showing high VWR activity (High Runners; HRs), together with 4 non-selected Control lines, over 88 generations. We found that several sensory receptors in specific receptor clusters were differentially expressed between the vomeronasal organ (VNO) of HRs and Controls. Moreover, one of those clusters contains multiple single-nucleotide polymorphism loci for which the allele frequencies were significantly divergent between the HR and Control lines, i.e., loci that were affected by the selective breeding protocol. These results indicate that the VNO has become genetically differentiated between HR and Control lines during the selective breeding process. Although the role of the vomeronasal chemosensory receptors in VWR activity remains to be determined, the current results suggest that these vomeronasal chemosensory receptors are important quantitative trait loci for voluntary exercise in mice. We propose that olfaction may play an important role in motivation for voluntary exercise in mammals.

Introduction

Chemical senses are involved in many aspects of behavior. Olfaction is especially important for controlling such intrinsically motivated behaviors as food-seeking, social interactions, and reproductive- and fear-driven behaviors [1]. An ethologically relevant cue is sensed by chemosensory receptors expressed in the sensory organs, which activate a specific neural circuitry for behavioral motivation and induces an appropriate behavioral output in a specific context [2-5]. Comparative functional studies involving insect model species propose a model wherein molecular evolution of chemosensory receptors is sufficient to induce changes in neural circuit activity and behavioral patterns [6]. Thus, ethologically relevant cues, receptors, neural circuitries, and behavioral habits are likely to evolve together (coadapt) in response to natural and sexual selection.

One olfactory organ, the vomeronasal organ (VNO), occurs in some amphibians, squamates, and some mammals, including rodents. The VNO is known to detect intraspecific signals known as pheromones that trigger behavioral and physiological changes in receivers [2]. Pheromones are non-volatile peptides and small molecular weight compounds that are excreted in such fluids as urine and tears. These molecules are taken up from the environment to the VNO by direct contact and activate the vomeronasal sensory neurons (VSNs) [7, 8]. Generally, each VSN expresses a member of the vomeronasal receptor families: type 1 vomeronasal receptors (Vmn1rs), type 2 vomeronasal receptors (Vmn2rs), and formyl peptide receptors (Fprs), with some exceptions [9-13]. The signals detected by these receptors in the VSNs are axonally sent to glomerular structures and synaptically transmitted to the postsynaptic neurons, also known as mitral-tufted cells, in the accessory olfactory bulb (AOB) [14, 15]. The signals are then processed in the amygdala and hypothalamus, which induce the animal’s instinctive behavioral responses and endocrinological changes [2, 16, 17].

Rapid evolution of the receptor genes is a pronounced feature of the vomeronasal system [18-31]. Different species of animals have divergent family members of vomeronasal receptor genes [23, 26–28, 32, 33]. Even within the Mus musculus (house mouse) species complex, variation in the coding sequence is frequently observed [18]. Moreover, the abundance of receptor genes expressed in the VNO varies even among different inbred mouse strains [34]. Distributions of single nucleotide polymorphisms (SNPs) observed in lab-derived strains are non-random, and correlated with vomeronasal receptor phylogeny as well as genomic clusters [18]. These observations led us to hypothesize that selective breeding for a behavior that is modulated by chemosensory signals would induce an alteration in genomic clusters of vomeronasal receptors that are potentially involved in the behavior.

Voluntary wheel running (VWR) is an intrinsically motivated behavior, and even wild mice run on a wheel placed in nature [35]. Notably, VWR is one of the most widely studied behaviors in laboratory rodents [36-38]. Individual differences in VWR are highly repeatable on a day-to-day basis, the trait is heritable within outbred populations of rodents, and genes and genomic regions associated with VWR are being identified [39]. Moreover, some of the underlying causes of variation in VWR have been elucidated, in terms of both motivation and ability for voluntary exercise [37, 40, 41]. Importantly, a previous study demonstrated that the presence of conspecific urine increased VWR activity level in adult wild-derived mice [42], suggesting that external chemosensory cues also have a modulatory role in VWR activity.

We have established 4 independent, artificially evolved mouse lines by selectively breeding individuals showing high VWR activity (High Runners; HRs), along with 4 independent, non-selected Control lines over 88 generations [43, 44]. Briefly, all 4 HR lines run ∼2.5–3.0-fold more revolutions per day as compared with the 4 Control lines [45, 46]. Studies of mice allowed access to clean wheels or those previously occupied by a different mouse revealed that HRs show higher sensitivity to previously-used wheels and display greater alteration in daily wheel running activities than the Controls [47]. This result suggests that selective breeding for high running activity was accompanied by altered sensitivity to other individuals, suggesting a potential coadaptation of the chemosensory system with voluntary wheel running.

In this study, we examined whether selective breeding for VWR has differentiated the vomeronasal receptor genes between HR and Control lines. We found that a repertoire of receptor genes was differentially expressed between the VNO of HR and Control lines, which resulted from reduction or increase of specific vomeronasal receptor-expressing cells in the VNO of HR lines. We also found that this gene expression change was partially due to the genetic alteration upon selective breeding for VWR, suggesting a relationship between high running activity and the function of the VNO in HR lines. Taken together, our results indicate vomeronasal receptors as quantitative trait loci (QTL) for voluntary exercise behavior in mice.

Results

Differential expression of vomeronasal receptors in the VNO of HR and Control lines

To examine the impact of selective breeding for VWR activity on receptor gene expression in the VNO, we conducted transcriptome analysis of the VNO from both males and females of HR and Control lines. For each sex and replicate line of the 4 HR and 4 Control lines, total RNA samples were prepared, each consisting of the combined VNOs from 3 individual mice. After RNA sequencing, we identified 132 differentially expressed (DE) genes in the HR line group compared to the Control line group (Fig 1). We identified 19 vomeronasal receptor genes in the DE gene set, which belong to either the Fpr, Vmn1r or Vmn2r family (Fig 1, shown in red). Of the 19 DE receptor genes, the Reads Per Kilobase Million (RPKM) of Fpr3, Fpr-rs4, Vmn2r8, Vmn2r9, Vmn2r11, Vmn2r96, Vmn2r98, Vmn2r102, Vmn2r107, Vmn2r110, and Vmn2r114 were significantly up-regulated, while Vmn1r188, Vmn1r236, Vmn2r13, Vmn2r15, Vmn2r16, Vmn2r23, Vmn2r45, and Vmn2r99 were significantly down-regulated in the VSNs of HR lines compared to Control lines (Fig 2A). The RPKM of olfactory marker protein (OMP), which is abundantly and exclusively expressed in all mature VSNs in the VNO [48], was not different between HR and Control lines (data not shown), indicating that receptor gene expression changes were not due to variation in VSN number. The log₂ fold change of normalized expression of the DE genes varied from -3.3 to 2.1 (Fig 2B). Fpr-rs4 and Vmn2r16 showed the largest upregulation and downregulation, respectively. To examine sexually-biased expressions, we separated RPKMs of each DE receptor gene into sex and linetype groups (S1 Fig) and performed a two-way ANOVA test to examine sex and linetype differences. As shown in S2 Fig, p values for sexually-biased expression and interaction were above 0.05 in all DE receptor genes, demonstrating that there are no sex-dimorphic expression changes between HR and Control lines. These results suggest that expression of the chemosensory receptor genes is differentially regulated in the VSNs between HR and Control lines without a sex difference.

Fig 1

RNAseq analysis of the vomeronasal organs of High Runner and Control mice.

Heatmap of differentially expressed (DE) genes between HR and Control lines. DE vomeronasal receptors are shown in red. Fpr3 (analyzed in Fig 5) is highlighted with a red underline.

Fig 2

RPKM comparison of DE vomeronasal receptors between High Runner and Control mice.

(A) Scatter plots showing the RPKM of DE vomeronasal receptor genes in each line of HR or C mice. Black bars represent medians. (B) Bar graph denoting log2 fold change of the relative expression of DE vomeronasal receptor genes between the HR and Control lines.

Fig 5

RNA scope in situ hybridization analysis of a DE vomeronasal receptor gene in the VNO.

(A) Images showing RNAscope-derived fluorescent signals for Fpr3 (left) and Gnao1 (middle) transcripts. In merged images (right), Fpr3 and Gnao1 are shown in red and yellow, respectively, together with DAPI staining (blue). Upper and lower panels show representative images from the VNO of a Control (line 5) line and a HR (line 3) line, respectively. (B) A scatter plot showing the percentage of Fpr3 signals in vomeronasal sensory neurons (VSNs) per VNO slice for each line of mice. n = 12 slices in 2–3 mice per line. Black bars represent medians. (C) The mean percentage of Fpr3 signals in VSNs in Control and HR lines. Each dot indicates the mean of one line. Black bars represent medians. * indicates p < 0.05 in Mann-Whitney test. (D) Representative images showing Fpr3 signals (red) observed in the VNO of Control and HR lines of mice. DAPI signals are shown in blue. (E) The mean of Fpr3 signal intensity (arbitrary unit, AU) per VSN in the Control and HR lines. Each dot indicates the mean of one line. Black bars represent medians.

Interestingly, 15 out of 19 DE receptor genes are co-localized in specific vomeronasal receptor clusters. Vmn2r8, Vmn2r9, Vmn2r11, Vmn2r13, Vmn2r15, and Vmn2r16 are localize in a ~1 Mb vomeronasal receptor cluster on chromosome 5, and Fpr3, Fpr-rs4, Vmn2r96, Vmn2r98, Vmn2r99, Vmn2r102, Vmn2r107, Vmn2r110, and Vmn1r236 are localized in a ~3 Mb vomeronasal receptor cluster on chromosome 17 (Fig 3). Vmn2r23, Vmn2r45, Vmn1r188, and Vmn2r114 exist on chromosome 6, 7, 10, and 17, respectively (Fig 3), and they are surrounded by vomeronasal receptors which are not differentially expressed between HR and Control lines. These results suggest the correlation of differential expression of vomeronasal receptors and genomic chromosomal locations, as well as other unknown causes.

Fig 3

Chromosomal locations of DE vomeronasal receptors.

A plot showing the log2 fold change of 336 vomeronasal receptor genes in their relative positions along the chromosomes. Each circle represents a member of Vmn1r, Vmn2r, or Fpr family genes. DE vomeronasal receptors are shown in filled circles.

Accumulation of all-or-none SNPs in a vomeronasal receptor cluster

We then hypothesized that differences between HR and Control lines in vomeronasal receptor gene expression may be associated with differences in genomic sequences, especially allele frequencies between HR and Control lines caused by the selective breeding. Previous genome-wide SNP analysis detected 152 out of 25,318 variable SNP loci for which allele frequencies are significantly different between HR and Control lines after correction for multiple comparisons [49]. As explained in the previous paper, the differentiation in allele frequencies for these 152 loci cannot be attributed to random genetic drift. Of the 152 SNP loci, we particularly focused on 61 loci that are fixed for the same allele in all 4 replicate HR lines but not fixed in any of the 4 replicate Control lines, or vice versa (which we term “all-or-none SNPs”, S1 Table). The 61 SNP loci are not randomly distributed throughout the genome (S3A Fig). The majority of them (59 of 61) exist as a member of groups of 3 or more which are located in close proximity on the genomic chromosomes (S3B Fig). As a result, only 11 all-or-none SNP clusters are observed in the genome (S3A Fig)

Interestingly, 8 of the 61 all-or-none SNP loci are located in a ~3 Mb interval on chromosome 17 that contains clusters of Vmn1rs (14), Vmn2rs (21), and Fprs (7) (Fig 4). Strikingly, 9 out of the 19 DE vomeronasal receptors are located in this all-or-none SNP cluster. Five of the all-or none SNPs are localized near the differentially expressed vomeronasal receptors (Table 1): SNP ID rs29503987 at 30 kb downstream and 20 kb upstream of Fpr3 and Fpr-rs4, respectively, rs33447983 at 8.4 kb downstream of Vmn2r99, rs6224641 at 30 kb downstream of Vmn2r99, rs33649277 at 44 kb upstream of Vmn2r102, rs29522462 at 524 bp downstream of Vmn2r110, and rs33120398 at the intron 1 of Vmn2r110. Other SNPs, such as rs33463529, are also closely located near vomeronasal receptors, though changes in expression of the nearby receptors were not observed. These results strongly suggest that changes in vomeronasal receptor gene expression between HR and Control lines are at least partially caused by changes in allele frequencies at multiple loci in response to selective breeding for VWR activity.

Fig 4

A genomic cluster containing both SNPs and DE vomeronasal receptors.

A genomic cluster containing DE vomeronasal receptor genes in the mouse chromosomes 17. Vomeronasal receptors in red and blue indicate up- and down-regulations, respectively. Non-DE genes are shown in black. Purple arrowheads indicate locations of SNPs that are significantly differentiated between HR and Control groups [49], as shown in Table 1.

Table 1

Locations of SNPs on the chromosome 17 that are significantly differentiated between HR and Control groups.

	SNP ID	SNP location (mm10)	HR / Cont	Gene (within 50kb)	strand	up/dw	Distance (bp)
1	rs29503987	18001459	C / A	Fpr-3	+	dw	29935
				Fpr-rs4	+	up	20274
				Vmn2r124	+	up	48025
2	rs33375308	18210739	T / C	Vmn2r92	+	dw	25561
3	rs33447983	19403011	G / T	Vmn2r99	+	dw	8421
4	rs6224641	19424358	A / G	Vmn2r99	+	dw	29768
5	rs33649277	19616228	G / A	Vmn2r101	+	dw	3911
				Vmn2r102	+	up	44171
6	rs29522462	20573305	A / C	Vmn2r109	-	up	8549
				Vmn2r110	-	dw	524
7	rs33120398	20587484	C / T	Vmn2r109	-	up	22728
				Vmn2r110	-	intron 1
				Fpr-rs3	-	dw	36362
8	rs33463529	20779567	G / A	Vmn1r227	+	dw	43441
				Vmn1r228	-	up	2066
				Vmn1r229	+	up	34928

The rest of the DE vomeronasal receptors, which are located in another vomeronasal receptor cluster on chromosome 5, as well as other solitary ones, are not surrounded by SNP loci that are significantly differentiated between HR and Control lines [49]. Therefore, expression changes of the receptor genes may be mediated by SNPs that remain polymorphic in both lines, or by different mechanisms.

Differential number of Fpr3-expressing VSNs in the VNO of HR and Control lines

To determine the significance at the cellular level in the VNO of the DE chemosensory receptor genes, we chose one representative gene to determine whether there are differences in the number of receptor-expressing VSNs, or alternatively, differences in transcript abundance in each receptor-expressing VSN. We performed in situ hybridization to detect Fpr3 in the VNO of 2–3 individual mice from each of the 4 HR and 4 Control lines, together with a probe for the Gαo (Gnao1). Expression of Fpr3 is ~3 times higher in HR lines compared to Control lines in RNAseq analysis (Fig 2B). Although Fpr3-expressing VSNs were observed in the VNO of both HR and Control lines, the number of Fpr3-expressing VSNs in each VNO slice varied among lines (Fig 5A and 5B). In 3 Control lines (line 1, 4, and 5), Fpr3 signal was barely observed in each VNO tissue slice, while Control line 2 had a significantly higher number of Fpr3-expressing VSNs (Fig 5B). This result was consistent with the RNAseq data, in which the amount of Fpr3 transcripts in line 2 was higher than other Control lines (Fig 2A). We consistently observed multiple Fpr3-expressing VSNs in most of the VNO tissue slices from the 4 HR lines and Control line #2 (Fig 5B and 5C). Generally, there were significantly more Fpr3-expressing VSNs in the HR versus the Control lines (Fig 5C, Mann-Whitney test (p < 0.05)). The fluorescent intensity derived from Fpr3 gene transcripts in each VSN was not distinguishable between the VNO tissues from HR and Control lines (Fig 5D and 5E). Taken together, these results demonstrate that the different expression levels of the chemosensory receptors result from changes in the number of the receptor-expressing VSNs. The results observed here are consistent with a previous finding in the main olfactory system, which has a similar monoallelic pattern of receptor gene expression in each sensory neuron [50]: a positive correlation between tissue RNA levels of olfactory receptor genes and numbers of OSNs expressing the receptors. Thus, the number of VSNs expressing specific sets of chemosensory receptors are differentially regulated after selective breeding for VWR.

Discussion

In this study, we utilized a unique animal model: 4 replicate mouse lines that have been experimentally evolved by selectively breeding individuals showing high VWR activity (HR lines), along with their 4 independent, non-selected Control lines maintained over 88 generations [43]. The HR and Control lines provide a strong model for determining the contribution of genetics to voluntary-exercise related traits [44, 51]. In addition to the exercise ability-related genetic adaptations found after selective breeding [37, 44], several changes at the level of the central nervous system have also been identified, which contribute to elevation of VWR for HR mice [37, 40, 52]. Through SNP mapping analysis (S1 Fig), we found that 3 of the 61 all-or-none SNP loci that were fixed in all 4 replicate HR lines (but none of the 4 replicate Control lines) were located in a genomic cluster exclusively containing T-box genes on Chromosome 5. These genes are associated with GO terms of “bundle of His development,” “embryonic forelimb morphogenesis,” “cardiac septum morphogenesis,” “ventricular septum development,” and “cardiac muscle cell differentiation”. Indeed, compared with their 4 non-selected Control line counterparts, mice from the 4 replicate HR lines have been shown to have increased ventricular mass [45, 53–55], as well as altered cardiac functions [55-57]. Thus, the genome-wide SNP analysis of HR and Control lines of mice [49] could robustly identify QTL associated with voluntary exercise behavior.

Vomeronasal receptors are among the most rapidly evolving genes in vertebrates [18-31]. Different taxonomic groups have divergent family members of vomeronasal receptor genes [21, 23, 26–28, 32, 33], and the abundance of receptor genes expressed in the VNO is different even among inbred mouse strains [34]. Moreover, many of the mouse pheromones identified as ligands for vomeronasal receptors show strain specificity. For example, expression of the male pheromone ESP1 is only observed in a few inbred strains, although males of wild-derived strains all secrete abundant ESP1 peptide into their tears [58]. Likewise, expression of juvenile pheromone ESP22 is missing in some inbred strains [59]. Major urinary proteins (MUPs) are potential ligands for vomeronasal receptors, and all male mice of a given inbred strain secrete identical MUP members, whereas wild-derived mice each exhibit a unique profile of emitted MUPs [60]. Thus, pheromones and receptors in the vomeronasal system may have evolved in response to various environmental changes, including domestication, which resulted in alteration of coding sequences and expression patterns.

Considering the extensive evolution of receptor genes, selective breeding for a chemosensory-mediated behavior is an attractive alternative approach to reveal the functions of vomeronasal receptors. VWR activity of a mouse strain that recently derived from the wild has been shown to be increased by urinary chemosignals from other individuals [42]. Therefore, if the function of the VNO is involved in the modulation of VWR activity, then we would expect that selective breeding for high VWR activity should impact vomeronasal receptors. Indeed, we found that expression levels of several vomeronasal receptor genes as well as a few SNPs near the DE receptor genes were different between HR and Control lines. Although the role of each DE receptor in VWR activity needs to be determined in future studies, the current results suggest that vomeronasal chemosensory receptors could be important QTLs for voluntary exercise in mice.

One of the important remaining questions is how the vomeronasal system modulates VWR behavior in HR lines. One study measuring patterns of brain activity using c-Fos immunoreactivity revealed multiple areas in the brain that appear to be associated with motivation for VWR in HR lines [61]. These areas include brain nuclei known to be motivation-related, such as the prefrontal cortex, medial frontal cortex, and nucleus accumbens (NAc) [61]. In addition, it was recently shown in mice that VNO-mediated signals regulate the mesolimbic dopaminergic system, especially by upregulating the ventral tegmental area (VTA)-NAc circuit, and that they enhance reproductive motivation in mice [15]. Thus, it is possible that the VNO-mediated chemosensory signals also upregulate VWR activity by stimulating the VTA-NAc circuit. Moreover, one of the hypothalamic targets of the vomeronasal system, the medial preoptic area (MPOA), has been shown to regulate wheel-running activity in a hormone-dependent manner [62-65]. It is therefore also conceivable that the VNO-mediated chemosensory signals upregulate VWR by directly activating MPOA neurons. Combined with these previous observations, we propose that chemosensory signals detected by the VNO activate specific areas of the central nervous system that contribute to VWR activity. Future studies are expected to reveal the role of the VNO in modulating physical exercise and other voluntary behaviors in rodents.

Materials and methods

Animals

The experimental procedures were approved by the UCR Institutional Animal Care and Use Committee and were in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. The VNOs studied were from 12-week old male and female mice of 4 lines selected for high voluntary wheel running and 4 Control lines. The studied mice were derived from generation 88 of a replicated selective breeding experiment for increased voluntary wheel running behavior that began with a base population of the outbred Hsd:ICR strain [43]. Wheel revolutions were recorded in 1-minute intervals continuously for 6 days, and mice were selected within-family for the number of revolutions run on days 5 and 6. In each selected HR line, the highest-running male and female within 10 individual families were selected per generation and each mouse was mated to a mouse from another family, within its line. This within-family selection regimen minimized inbreeding such that the effective population size was approximately 35 in each line [43]. In the Control lines, one female and one male within each family were chosen at random, though full sibling mating was again prevented. The mice in the present study were neither full nor half-siblings.

RNA sequencing

The VNO tissues were harvested from 3 male and 3 female mice from each of the 4 HR and 4 Control lines, immediately transferred to RNA later (Sigma-Aldrich), then stored at −80˚C until use for RNA-seq. VNO tissues from the same sex and line of mice were pooled and homogenized in Trizol Reagent (Life Technologies, Carlsbad, CA) and processed according to the manufacturer's protocol. Trizol-purified RNA samples were quantified using Qubit1 2.0 (Life Technologies). The integrity of isolated RNA was measured by the 28S/18S rRNA analysis using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara CA) with RNA Nano chip (Agilent Technologies, Palo Alto, CA). Samples had RNA integrity number values of at least 8.30. Using the Ultra II Directional RNA Library Prep kit (NEB), each RNA sample was depleted of ribosomal RNA and used to prepare an RNA-seq library tagged with a unique barcode at the UCR IIGB Genomics Core. Libraries were evaluated and quantified using Agilent 2100 Bioanalyzer with High Sensitivity DNA chip, then sequenced with the Illumina NextSeq 500 system (Illumina, San Diego, CA, USA) and 75nt-long single-end reads were generated at the UCR IIGB Genomics Core. A total of 8 libraries (4 HR lines and 4 Control lines) were multiplexed and sequenced in a single lane which yielded ~11,000 M reads, averaging ~1,400 M reads per sample.

The RNA-seq data files are available in the National Center for Biotechnology Information Gene Expression Omnibus (GEO) database (accession identifier GSE146644).

Differential gene expression analysis

The analysis compared the transcriptome profiles from both males and females of the HR and Control lines of mice. Quality control of the sequence reads included a minimum average Phred score of 30 across all positions using FastQC. Sequencing reads were aligned to the mouse reference genome (GRCm38/mm10), using STAR aligner ver. 2.6.1d [66] with an increased stringency unforgiving any of mismatches per each read (‘-outFilterMismatchNmax 0’). Any reads that map to multiple locations in the genome are not counted (‘-outFilterMultimapNmax 1’) since they cannot be assigned to any gene unambiguously. In order to determine the differentially expressed (DE) genes, generated BAM files were accessed with Cuffdiff [67], a program included in Cufflinks. Cuffdiff reports reads per kp per million mapped reads (RPKM), log₂ fold change, together with p-value, and adjusted p-values (q-values). After Benjamini-Hochberg false discovery correction, genes with q-values less than 0.05 were considered as DE genes. To examine sexual differences, RPKM of DE vomeronasal receptor genes in male and female HR and Control mice were subjected to two-way ANOVA tests using Prism (GraphPad). Notably, p-values for linetype-biased expression of 4 DE receptor genes (Vmn2r13, Vmn2r23, Vmn2r45, and Vmn2r107) were > 0.05 due to differences of statistical tests: the negative binomial regression in DE gene detection and the general linear model in two-way ANOVA.

Analysis of all-or-none SNPs

SNP data in supplemental table 7 (Data_S7) in Xu and Garland (2017) [49] were used for this analysis. SNPs that separate all 4 HR and 4 Control lines (which we term all-or-none SNPs) were selected (S1 Table) and mapped onto mouse genome (NCBI37/mm9) using UCSC Genome Browser (https://genome.ucsc.edu). We noticed that most of the all-or-none SNPs occurred in groups. Thus, we mapped those SNP clusters onto genomes (S3A Fig). Each cluster was defined—and + 0.1 Mb from the first and last SNP, respectively, observed in a specific location of the genome. Information of coding genes in each SNP cluster were extracted (S3B Fig). For some clusters, Gene-to-GO mappings was performed with PANTHER (http://pantherdb.org).

RNAscope in situ hybridization

Female mice (11 Controls and 10 HRs) were utilized for this analysis. The animals were intracardially perfused with 4% Paraformaldehyde in Phosphate Buffered Saline (PBS). VNOs were dissected from perfused animals and fixed overnight. The VNO samples were decalcified in EDTA pH 8.0 for 48 hours, then cryoprotected in 15% sucrose in PBS followed by 30% sucrose in PBS. All samples were ultimately embedded in optimal cutting temperature (OCT) medium (Electron Microscopy Sciences) above liquid nitrogen and sectioned at 20 μm using Leica CM3050S Cryostat. The cross sections analyzed were from the VNO regions with clearly discernible two crescent shapes. They were collected approximately 120 μm apart from each other spanning approximately 360 μm of the VNO containing regions from each mouse. Slides were stored at -80°C until use for in situ hybridization staining.

RNA detection in VNOs were performed with ACD RNAscope® control and target GNAO1 (ACD # and FPR3 (ACD #503451) using RNAscope® Multiplex Fluorescent Reagent Kit v2 (ACD# 323100) Assay. Probe binding was detected with Akoya Biosciences’ Opal 690 (FP1497001KT) and 570 (FP1488001KT) Dyes at 1:750 dilution in RNAscope TSA Buffer. Nuclear staining was visualized with DAPI (EMS #17989–20). Images were acquired at 20X or 40X magnification on Zeiss Axio Imager.M2, and FPR3-positivity was quantified with a proprietary script using QuPath software [68]. Fluorescent intensity was measured by Fiji software. 4–8 slices in each animal were examined. Mann-Whitney test was used to examine statistical significance in Fig 5C and 5E.

RPKM comparison of DE vomeronasal receptors between males and females in High Runner and Control mice.

Scatter plots showing the RPKM of DE vomeronasal receptor genes in males and females from each line of HR or Control mice. Black bars represent medians.

(TIF)

Click here for additional data file.

Two-way ANOVA tests of DE vomeronasal receptors.

A table showing p-values for interactions, sex differences, and linetype differences in the RPKM of DE vomeronasal receptor genes in two-way ANOVA analyses. n.s., *, **, ***, and **** represent not significant, p < 0.05, p < 0.01, p < 0.001, and p < 0.0001, respectively.

(TIF)

Click here for additional data file.

Analysis of all-or-none SNP loci in HR and Control lines of mice.

(A) A schematic diagram showing the relative positions of loci containing 1 or more all-or-none SNPs. Blue triangles indicate non-chemosensory clusters, and a red triangle indicates clusters containing only chemosensory (vomeronasal) receptors. (B) A table showing chromosomal location and length of each all-or-non SNP cluster, and the number of SNPs and genes within the clusters. The row highlighted in red is the cluster containing only the vomeronasal receptor genes.

(TIF)

Click here for additional data file.

Genomic locations, p-value by the mixed model approach [49], and allele frequencies of the 61 all-or-on SNP loci.

(XLSX)

Click here for additional data file.