Etienne Loire

Sergio Tusso

Pierre Caminade

Dany Severac

Pierre Boursot

Guila Ganem

Carole M Smadja
email: carole.smadja@univ-montp2.fr

§ These

Do changes in gene expression contribute to sexual isolation and reinforcement in the house mouse?

Keywords: assortative mate preference, speciation, hybrid zone, RNA-seq, olfaction, vomeronasal receptors

Expression divergence, rather than sequence divergence, has been shown to be important in speciation, particularly in the early stages of divergence of traits involved in reproductive isolation. In the two European subspecies of house mice, Mus musculus musculus and M. m. domesticus, earlier studies have demonstrated olfactory-based assortative mate preference in populations close to their hybrid zone. It has been suggested that this behaviour evolved following the recent secondary contact between the two taxa (~3,000 years ago) in response to selection against hybridisation. To test for a role of changes in gene expression in the observed behavioural shift, we conducted a RNA-sequencing experiment on mouse vomeronasal organs. Key candidate genes for pheromone-based subspecies recognition, the vomeronasal receptors, are expressed in these organs. Overall patterns of gene expression varied significantly between samples from the two subspecies, with a large number of differentially expressed genes between the two taxa. In contrast, only ~200 genes were found repeatedly differentially expressed between populations within M. m. musculus that did, or did not displayassortative mate preferences (close to, or more distant from the hybrid zone, respectively), with an overrepresentation of genes belonging to vomeronasal receptor family 2. These receptors are known to play a key role in recognition of chemical cues that handle information about genetic identity. Interestingly, four out of five of these differentially expressed receptors belong to the same phylogenetic cluster, suggesting specialisation of a group of closely-related receptors in the recognition of odorant signals that may allow subspecies recognition and assortative mating.

Introduction

Understanding the genetic basis of speciation remains a key challenge in evolutionary biology [START_REF] Seehausen | Genomics and the origin of species[END_REF], in particular determining the genetic architecture of differentiation and reproductive isolation [START_REF] Nosil | Genomic divergence during speciation: causes and consequences[END_REF], the types of genes involved [START_REF] Nosil | The genes underlying the process of speciation[END_REF], and the types of genetic changes responsible for the evolution of reproductive isolation between diverging taxa [START_REF] Hoekstra | The locus of evolution: Evo devo and the genetics of adaptation[END_REF]. Although evidence is rapidly accumulating for the role of protein coding and structural changes in adaptation and reproductive isolation [START_REF] Stapley | Adaptation genomics: the next generation[END_REF], it has also long been argued that gene expression differences may be of prime importance in species differentiation [START_REF] King | Evolution at two levels in humans and chimpanzees[END_REF]Carroll et al. 2005), allowing rapid divergence [START_REF] Wolf | Nucleotide divergence vs. gene expression differentiation: comparative transcriptome sequencing in natural isolates from the carrion crow and its hybrid zone with the hooded crow[END_REF]. Moreover, measures of gene expression provide an important bridge between genotype and phenotype [START_REF] Huestis | From Gene Expression to Phenotype in Insects: Nonmicroarray Approaches for Transcriptome Analysis[END_REF]. In recent years, the study of gene expression has been greatly facilitated by high-throughput sequencing-based methods such as RNA-seq [START_REF] Mortazavi | Mapping and quantifying mammalian transcriptomes by RNA-Seq[END_REF][START_REF] Wang | RNA-Seq: a revolutionary tool for transcriptomics[END_REF], contributing to understand the importance of expression variation in both local adaptation and speciation.

Gene expression changes and the evolution of reproductive isolation

There are many examples of differential gene expression causing adaptive phenotypic changes (e.g., Chan et al. 2010a;[START_REF] Mcbride | Evolution of mosquito preference for humans linked to an odorant receptor[END_REF], and large-scale transcriptomics studies have shown heritable expression differences between diverging populations (e.g., [START_REF] Wolf | Nucleotide divergence vs. gene expression differentiation: comparative transcriptome sequencing in natural isolates from the carrion crow and its hybrid zone with the hooded crow[END_REF][START_REF] Uebbing | Divergence in gene expression within and between two closely related flycatcher species[END_REF][START_REF] Davidson | Gene Regulatory Evolution During Speciation in a Songbird[END_REF]. How these components of gene expression can be associated with reproductive isolation is still a largely unexplored question. The major advances on this point come from studies investigating the role of gene regulation in the evolution of postzygotic isolation and, in particular, how misexpression in hybrids can reflect Dobzhansky-Muller incompatibilities (e.g. [START_REF] Haerty | Gene regulation divergence is a major contributor to the evolution of Dobzhansky-Muller incompatibilities between species of Drosophila[END_REF][START_REF] Ortiz-Barrientos | Gene expression divergence and the origin of hybrid dysfunctions[END_REF]Dion-Cote, Renaut, Normandeau, & Bernatchez, 2014;[START_REF] Davidson | Gene Regulatory Evolution During Speciation in a Songbird[END_REF][START_REF] Mack | Gene regulation and speciation in house mice[END_REF]. A recent study also identified significant changes in expression associated with a postmating-prezygotic reproductive barrier between two Solanum plant species, that involve effector molecules possibly responsible for pollen tube arrest [START_REF] Pease | Molecular mechanisms of postmating prezygotic reproductive isolation uncovered by transcriptome analysis[END_REF]. As far as premating isolation is concerned, evidence is mostly found in phytophagous insect species that are reproductively isolated due to host plant preference and specialisation. For example, recent studies have shown expression differences between populations specialised on different plants, involving genes related to host selection (odorant receptors, salivary proteins, detoxification genes) (e.g. cactus-specialised populations of Drosophila mettleri [START_REF] Hoang | Transcriptional variation associated with cactus host plant adaptation in Drosophila mettleri populations[END_REF] Rhagoletis host races: [START_REF] Ragland | Differences in performance and transcriptome-wide gene expression associated with Rhagoletis (Diptera: Tephritidae) larvae feeding in alternate host fruit environments[END_REF]; pea aphid host races: [START_REF] Eyres | Differential gene expression according to race and host plant in the pea aphid[END_REF]). These differences in expression have therefore been proposed to play a key role for adaptive genetic divergence in traits causing reproductive isolation during ecological speciation [START_REF] Pavey | The role of gene expression in ecological speciation[END_REF][START_REF] Etges | No boundaries: genomes, organisms, and ecological interactions responsible for divergence and reproductive isolation[END_REF]. Despite growing evidence for the role of gene regulation in mate choice and sexual selection [START_REF] Wilkinson | The locus of sexual selection: moving sexual selection studies into the post-genomics era[END_REF], support for expression changes contributing to sexual isolation is relatively scarce, and again concerns insect species. Examples include sex-specific gene expression changes underlying shifts in sex-pheromone preferences between populations of the mosquito Anopheles gambiae (Cassone et al. 2008), the silkmoth Bombyx mori [START_REF] Fujii | Sex-linked transcription factor involved in a shift of sexpheromone preference in the silkmoth Bombyx mori[END_REF] and D. melanogaster [START_REF] Bailey | Female Drosophila melanogaster gene expression and mate choice: the X chromosome harbours candidate genes underlying sexual isolation[END_REF], and a change in expression of a fatty acid synthase responsible for cuticular hydrocarbon production inducing sexual isolation between D. serrata and D. birchii (Chung et al. 2014).

Gene expression and the evolution of assortative mate preference in a house mouse hybrid zone

In this study, we investigate the role of expression divergence in the evolution of olfactorybased assortative mate preference between two subspecies of the house mouse (Mus musculus) in order to gain insights into the genetic basis of sexual isolation in this mammalian system. The two European subspecies of the house mouse, M. m. musculus and M. m. domesticus have diverged in allopatry for 0.5 million years in the Indian sub-continent [START_REF] Boursot | Origin and radiation of the house mouse: mitochondrial DNA phylogeny[END_REF], before meeting secondarily in Europe around 5,000 years ago (and around 3,000 years ago in Denmark) (Cucchi et al. 2012). They now form a hybrid zone stretching from Norway to the Black Sea. Assortative mate preference between the two subspecies has been demonstrated in populations at the edges of the hybrid zone in Denmark in both subspecies, and in both males and females, but appears particularly strong in M. m. musculus populations [START_REF] Smadja | Subspecies recognition in the house mouse: a study of two populations from the border of a hybrid zone[END_REF], 2005;Smadja et al. 2004a;[START_REF] Ganem | Variation in mate preference across a house mouse hybrid zone[END_REF]. These results suggest a degree of sexual isolation between the two taxa. In contrast, populations of both subspecies occurring in distant allopatry or slightly further away from the Danish hybrid zone do not display any directional mate preference [START_REF] Smadja | Asymmetrical reproductive character displacement in the house mouse[END_REF][START_REF] Smadja | Seeking signatures of reinforcement at the genetic level: a hitchhiking mapping and candidate gene approach in the house mouse[END_REF]. This pattern of reproductive character displacement has been confirmed in the Czech part of the hybrid zone [START_REF] Bimova | Reinforcement selection acting on the European house mouse hybrid zone[END_REF]. Since selection against hybrids is acting in the mouse hybrid zone (reduced hybrid fertility: [START_REF] Britton-Davidian | Postzygotic isolation between the two European subspecies of the house mouse: estimates from fertility patterns in wild and laboratory-bred hybrids[END_REF][START_REF] Albrechtova | Sperm-related phenotypes implicated in both maintenance and breakdown of a natural species barrier in the house mouse[END_REF][START_REF] Turner | Reduced male fertility is common but highly variable in form and severity in a natural jouse mouse hybrid zone[END_REF][START_REF] Bhattacharyya | Mechanistic basis of infertility of mouse intersubspecific hybrids[END_REF][START_REF] Turner | Genome-Wide Mapping In A House Mouse Hybrid Zone Reveals Hybrid Sterility Loci And Dobzhansky-Muller Interactions[END_REF]; sexual selection against hybrids: [START_REF] Latour | Sexual selection against natural hybrids may contribute to reinforcement in a house mouse hybrid zone[END_REF]; hybrid microbiome dysfunction: [START_REF] Wang | Analysis of intestinal microbiota in hybrid house mice reveals evolutionary divergence in a vertebrate hologenome[END_REF] a reinforcement process has been hypothesised in this system, by which premating barriers could have evolved as a response to selection against hybrids [START_REF] Bimova | Reinforcement selection acting on the European house mouse hybrid zone[END_REF][START_REF] Smadja | Asymmetrical reproductive character displacement in the house mouse[END_REF][START_REF] Smadja | Seeking signatures of reinforcement at the genetic level: a hitchhiking mapping and candidate gene approach in the house mouse[END_REF].

A current challenge is to determine the genetic basis of assortative mate preferences between the two subspecies and of the behavioural shift within each subspecies between populations distant from or close to the hybrid zone. Previous studies have shown that chemical signals present in urine are sufficient to elicit discrimination and assortative mate preferences between the two subspecies [START_REF] Smadja | Subspecies recognition in the house mouse: a study of two populations from the border of a hybrid zone[END_REF], 2005;Smadja et al. 2004a;[START_REF] Ganem | Variation in mate preference across a house mouse hybrid zone[END_REF][START_REF] Bimova | Reinforcement selection acting on the European house mouse hybrid zone[END_REF], suggesting a role of the olfactory system in this behavioural shift. In mammals, two organs are responsible for odorant signal recognition: the main olfactory epithelium (MOE), primarily involved in recognition of airborne chemicals present in the environment (odorants and pheromones), and the accessory olfactory organ or vomeronasal organ (VNO), known to play a key role in both volatile and non-volatile kairomone and pheromone recognition [START_REF] Isogai | Molecular organization of vomeronasal chemoreception[END_REF]Chamero, Leinders-Zufall, & Zufall, 2012;[START_REF] Liberles | Mammalian Pheromones[END_REF][START_REF] Bear | The Evolving Neural and Genetic Architecture of Vertebrate Olfaction[END_REF]. Recognition in these organs involves chemoreceptors (olfactory receptors (ORs) in the MOE and vomeronasal receptors (VRs) in the VNO) that are organised in very large multigene families in the mouse genome (more than 1,200 ORs and 530 VRs, Ibarra-Soria, [START_REF] Ibarra-Soria | The genomic basis of vomeronasal-mediated behaviour[END_REF][START_REF] Ibarra-Soria | The Olfactory Transcriptomes of Mice[END_REF][START_REF] Zhang | The olfactory receptor gene superfamily of the mouse[END_REF]. Although duplications, mutations, and adaptive evolution certainly play an important role in the evolution of chemoreceptors and their associated behaviours [START_REF] Horth | Sensory genes and mate choice: Evidence that duplications, mutations, and adaptive evolution alter variation in mating cue genes and their receptors[END_REF][START_REF] Karn | A Candidate Subspecies Discrimination System Involving a Vomeronasal Receptor Gene with Different Alleles Fixed in M. m. domesticus and M. m. musculus[END_REF][START_REF] Smadja | Seeking signatures of reinforcement at the genetic level: a hitchhiking mapping and candidate gene approach in the house mouse[END_REF], variation in expression at these genes can also strongly influence mouse behaviour [START_REF] Keller | The main and the accessory olfactory systems interact in the control of mate recognition and sexual behavior[END_REF][START_REF] Zhang | High-Throughput Microarray Detection of Vomeronasal Receptor Gene Expression in Rodents[END_REF][START_REF] Bear | The Evolving Neural and Genetic Architecture of Vertebrate Olfaction[END_REF]. A recent study has demonstrated that changes in expression patterns of vomeronasal receptors can induce changes in patterns of mate preference between laboratory strains of the house mouse [START_REF] Broad | The Post-Natal Chemosensory Environment Induces Epigenetic Changes in Vomeronasal Receptor Gene Expression and a Bias in Olfactory Preference[END_REF], confirming the key role of the vomeronasal organ and its receptors in the recognition and transduction of olfactory signals used in mate choice, as well as the possible link between variation in VR expression and variation in mate preference in mice.

Predictions and experimental design

We hypothesized that the observed behavioural divergence between populations close to the hybrid zone and displaying strong assortative mate preferences (hereafter 'Choosy' populations) and populations further away from the hybrid zone and not displaying any directional mate preferences (hereafter 'Non-Choosy' populations) could originate from divergence in expression of some vomeronasal receptors. To test this hypothesis, we conducted deep RNA-sequencing of whole vomeronasal organs, an approach that has already been validated for detection of expression variation at VRs [START_REF] Ibarra-Soria | The Olfactory Transcriptomes of Mice[END_REF].

We compared adult female and male mice originating from different 'Non-Choosy' and 'Choosy' populations of M. m. musculus and M. m. domesticus. Estimates of differential expression allowed us to identify specific sets of vomeronasal receptors and quantify differences in receptor expression potentially related to the shift toward assortative mate preference in the hybrid zone.

MATERIAL AND METHODS

Samples

Mice used in the study Patterns of gene expression were studied in adult mice that were born in the wild, trapped in Jutland, Denmark in October 2010 in several sites (indoor farms and other human dwellings) and then maintained in the laboratory under controlled conditions before being behaviourally tested and euthanised for dissection. As it is known that postnatal odour exposure can influence patterns of mate choice and VNO gene expression of adult mice [START_REF] Broad | The Post-Natal Chemosensory Environment Induces Epigenetic Changes in Vomeronasal Receptor Gene Expression and a Bias in Olfactory Preference[END_REF]Chamero et al. 2012;Cadiou et al. 2014), this experimental strategy allowed us to study behavioural and expression patterns characterising wild animals, while controlling for some factors (odorant environment, age range, diet) just before dissection that could inflate variance in expression among samples. Mice were trapped in three distinct sampling areas, characterised by populations with distinct mate preference behaviours (Figure S1): (1) the border of the hybrid zone on the M. m. musculus side (50 km North to the genetic centre of hybrid zone (defined in [START_REF] Raufaste | Inference of selection and migration in the danish house mouse hybrid zone[END_REF])) where strong assortative mating has been documented ('musculus Choosy' samples) (Smadja et al. 2004b[START_REF] Smadja | Seeking signatures of reinforcement at the genetic level: a hitchhiking mapping and candidate gene approach in the house mouse[END_REF][START_REF] Smadja | Asymmetrical reproductive character displacement in the house mouse[END_REF]Ganem et al. 2008), (2) another area in Denmark but further away from the hybrid zone, where mate preference was not assortative ('musculus Non-Choosy' samples) [START_REF] Smadja | Seeking signatures of reinforcement at the genetic level: a hitchhiking mapping and candidate gene approach in the house mouse[END_REF]; and (3) the border of the hybrid zone on the M. m. domesticus side (50 km South to the genetic centre of the hybrid zone), where weak assortative mating was documented ('domesticus' Choosy samples) [START_REF] Smadja | Asymmetrical reproductive character displacement in the house mouse[END_REF][START_REF] Ganem | Variation in mate preference across a house mouse hybrid zone[END_REF][START_REF] Latour | Sexual selection against natural hybrids may contribute to reinforcement in a house mouse hybrid zone[END_REF]). The M. m. domesticus samples were used as a control for the mapping of M. m. musculus reads against the house mouse reference genome, which is primarily of M. m. domesticus origin, and for interspecific comparisons. We decided to focus on the M. m. musculus side for comparisons between 'Choosy' and 'Non-Choosy' populations, because the strongest behavioural contrasts were observed there (Smadja et al. 2004b;[START_REF] Smadja | Asymmetrical reproductive character displacement in the house mouse[END_REF][START_REF] Ganem | Variation in mate preference across a house mouse hybrid zone[END_REF]. Two populations per geographical area were included as biological replicates, each composed of several trapping sites (six populations in total).

Mice were brought back to the laboratory in Montpellier (France) and maintained in quarantine under controlled and homogenized conditions for two months (same facility room, odorant environment, 12/12 photoperiod, and diet -food ad libitum). Only mice trapped as young adults were selected for this study, so that the ages of individuals participating to the behavioural and expression analyses were likely to have been within the age range at which VR expression in the VNO is optimal (between 90 and 120 days old, [START_REF] Zhang | High-Throughput Microarray Detection of Vomeronasal Receptor Gene Expression in Rodents[END_REF]. During the period in the laboratory, mice were tested for sex odour preferences to confirm assortative preferences in the 'musculus Choosy' group and an absence of directional preference in the 'musculus Non Choosy' group (protocol and results reported in [START_REF] Smadja | Seeking signatures of reinforcement at the genetic level: a hitchhiking mapping and candidate gene approach in the house mouse[END_REF]. We also confirmed weak assortative mate preferences in the 'domesticus Choosy' group [START_REF] Latour | Sexual selection against natural hybrids may contribute to reinforcement in a house mouse hybrid zone[END_REF]. We then selected eight males and eight females per population, resulting in a total of 96 adult male and female mice taking part to the RNA-seq experiment. The average mate preference of RNA-seq samples did not deviate from the average mate preference assessed in larger population samples. Geographical, genetic and behavioural characteristics of the mice used in this study are summarised in Table S1 (the method used to estimate hybridisation index of each mouse is described in detail in [START_REF] Latour | Sexual selection against natural hybrids may contribute to reinforcement in a house mouse hybrid zone[END_REF][START_REF] Smadja | Seeking signatures of reinforcement at the genetic level: a hitchhiking mapping and candidate gene approach in the house mouse[END_REF] and summarised in Table S1).

Vomeronasal organ dissection and RNA extraction

For each of the 96 female and male mice, the VNO was dissected rapidly after death by cervical dislocation. The VNO was immediately immersed in RNAlater and stored at -80C.

Total RNA was extracted using Macherey-Nagel XS kits or the phenol-chloroform based TRIzol procedure. RNA extracts were subsequently quantified and checked for quality with Nanodrop to assess sample concentration and purity and with Agilent Bioanalyzer to assess sample integrity. To constitute samples representative of each of the six populations and sexes, and to obtain enough material for library preparation, equimolar mixes from the eight males or eight females per population were pooled and again checked for quantity and quality using the Agilent Bioanalyzer and Nanodrop, to obtain 3 to 5 μg per pool (200ng/μL).

Trapping and transport of mice were performed under permission from the French Ministry of Agriculture (authorization n°C34-265). Housing took place in the facilities of the University of Montpellier and followed the ethical and animal welfare guidelines of the University of Montpellier and the regional ethical commitee.

Transcriptome sequencing and raw data analysis

Library preparation and sequencing RNA-seq libraries were prepared using the TruSeq RNA sample preparation kit (Illumina) and quantified using Agilent DNA 1000 Labchip kits (Agilent Technologies, Palo Alto, CA). The libraries were then sequenced on an Illumina HiSeq 2000 platform using a paired-end 100bp read-length protocol (3 libraries per lane, 4 lanes in total). Library preparation and sequencing were carried out at the MGX -Montpellier Genomix platform (Institut de Génomique Fonctionnelle-Institut de Génétique Humaine, Montpellier, France).

Filtering and mapping

Raw reads were filtered in order to discard reads with low quality sequence. The remaining reads were trimmed to subtract primer sequences and low quality 3' ends using cutadapt [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. The package RSEM [START_REF] Li | RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome[END_REF], with default mapping parameters, was used to map the reads against the house mouse reference genome and its known transcript annotation database (http://www.ensembl.org/, version GRCm38.p2) and to estimate gene expression abundance. Since the reference genome used for mapping corresponds to an inbred strain of the laboratory mouse (C57BL/6J) (primarily of M. m. domesticus origin), we compared the distribution of expression estimations between wild M. m. domesticus and M. m. musculus samples to detect any potential mapping biases among samples.

Identification of Vomeronasal Receptor and Olfactory Receptor transcripts

Coding sequences belonging to family 1 (V1R) and family 2 (V2R) of vomeronasal receptors [START_REF] Yang | Genomic organization and sequence analysis of the vomeronasal receptor V2R genes in mouse genome[END_REF][START_REF] Zhang | Comparative genomics of odorant and pheromone receptor genes in rodents[END_REF] as well as coding sequences of ORs were extracted with the BioMart tool from Ensembl (http://www.ensembl.org/). Information on chromosome location, genomic coordinates and strand for each gene was according to the NCBI-Build-37 mouse assembly (mm9).

Expression analyses

To compare the amount of gene expression among different classes of genes (e.g. VR genes versus non-VR genes), we normalised coverage values using the Bioconductor package CQN (conditional quantile normalization) (version 1.2.0) [START_REF] Hansen | Removing technical variability in RNA-seq data using conditional quantile normalization[END_REF], which accounts for possible bias due to library size, gene length and GC content. For comparisons of gene expression patterns among samples and differential expression analyses, read counts were normalised using the Bioconductor package edgeR (version 2.13) (McCarthy et al. 2012) implemented for the software R v.3.0.1 (http://www.r-project.org/) and the method of trimmed mean of M values (TMM) [START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF], which normalises for the library size and for the total number of expressed transcripts in each sample. A gene was considered as expressed in one sample if the number of reads was higher than 1 count per million (CPM).

Multidimensional scaling (MDS) plots were generated from the 500 genes maximizing the distances between points, to visualise similarity of global gene expression between samples.

Coverage at the gene level was obtained by summing up transcript coverage values associated with a given gene. To assess whether global patterns of gene expression vary between subspecies, mate preference profile groups ('musculus Choosy' versus 'musculus Non-Choosy') and sexes, we first verified equal dispersions among samples and then performed permutation-based MANOVA analyses [START_REF] Anderson | A new method for non-parametric multivariate analysis of variance[END_REF]) using the PERMDISP and PERMANOVA programs of the Vegan package Adonis [START_REF] Dixon | VEGAN, a package of R functions for community ecology[END_REF] with 1000 permutations.

Two models were used: the effects of the explanatory variables "Subspecies" and "Sex" were tested on the whole dataset (counts ~ Subspecies + Sex + Subspecies*Sex) and the effects of the explanatory variables "Choosiness" (two levels: Choosy and Non-Choosy) and "Sex" were tested on M. m. musculus samples only (counts ~ Choosiness + Sex + Choosiness*Sex). All explanatory variables were defined as fixed factors in our models. Differential expression analyses were conducted using edgeR (version 2.13) (McCarthy et al. 2012). The analyses compared the two subspecies ('domesticus Choosy' versus 'musculus Choosy' samples) or the two groups of M. m. musculus mice showing distinct behavioural profiles ('musculus Choosy' versus 'musculus Non-Choosy'). Since males and females display similar patterns of mate preference within each population (e.g., Smadja et al. 2004b[START_REF] Smadja | Seeking signatures of reinforcement at the genetic level: a hitchhiking mapping and candidate gene approach in the house mouse[END_REF][START_REF] Smadja | Asymmetrical reproductive character displacement in the house mouse[END_REF], we considered four biological replicates per categorial group (two sexes x two populations). Only genes expressed in all samples but one of the same categorical group were included, leading to a different number of genes considered for each comparison. Because simulation analyses have shown the tendency of edgeR to identify a high number of false differentially expressed genes [START_REF] Soneson | A comparison of methods for differential expression analysis of RNA-seq data[END_REF], the significance level was reduced to 0.01. Additionally, a false discovery rate (FDR) method was used to correct the threshold of P-values for multiple testing [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF].

The outcome of this analysis was a list of up-and down-regulated differentially expressed (DE) genes.

Enrichment analyses

Functional annotation of the significantly differentially expressed genes was performed using the functional classification program DAVID (Version 6.8) [START_REF] Huang | Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources[END_REF].

Enrichment analyses were performed for DE genes between the two subspecies and for DE genes between 'Choosy' and 'Non-Choosy' samples within M. m. musculus. Functional categories were considered significantly enriched when FDR < 0.05 and we analysed functional annotation clusters, which group similar functional categories together, to identify enriched biological functions. In the case of VR genes, we ran an independent hypergeometric test, as some of the genes are not annotated as VR genes in the database. In all cases, the reference background genome included only expressed genes for each comparison, instead of using the whole reference mouse genome.

Phylogenetic and clustering analyses

To map DE genes onto the phylogenies of vomeronasal receptor gene family 1 (V1R) and 2 (V2R), we first reconstructed these gene phylogenies using available coding sequences of V1R and V2R (191 and 110 sequences respectively). Sequences were aligned using MACSE (Multiple Alignment of Coding Sequences; [START_REF] Ranwez | MACSE: Multiple Alignment of Coding SEquences Accounting for Frameshifts and Stop Codons (WJ Murphy[END_REF], with manual adjustments. jModelTest v0.1.1 (Posada 2008) identified GTR+G+I (General Time Reversible + Gamma + Invariable sites) as the model of sequence evolution that best described the observed pattern of sequence variation for both V1R and V2R, and this model was used to perform maximum likelihood phylogenetic analyses using RAxML version 8 [START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF]). We then reported positions of DE genes identified previously between the two subspecies and between 'Choosy' and 'Non-Choosy' samples within M. m. musculus on the phylogenetic trees of the V1R and V2R families.

Since the phylogenetic analyses indicated a close relationship among five V2R genes that were found to be differentially expressed between the 'Choosy' and 'Non-Choosy' samples within M. m. musculus, we assessed the probability of phylogenetic clustering using a randomisation procedure. We randomly sampled five receptor genes in the phylogenetic tree of vomeronasal receptor family 2 and computed the total genetic distance among them using the ape R package [START_REF] Paradis | APE: Analyses of Phylogenetics and Evolution in R language[END_REF]). We performed this random sampling 500,000 times to obtain an expected distribution of genetic distance. The observed value of genetic distance among the five DE V2R genes was then compared to this distribution. Since all VR genes belong to only two different gene families, they share a significant amount of genetic identity. Since mapping was based on sequence identity, we were concerned that reads mapping to one truly differentially expressed gene may also map to other VR paralogs with close sequence identity, leading to false positive results in the clustering analysis.

Therefore, for each VR gene found to be differentially expressed in the 'musculus Choosy' versus 'musculus Non-Choosy' comparison, we retrieved the 10 closest paralogs with the program cd-hit-est (http://weizhong-lab.ucsd.edu/cd-hit/). Sequences were aligned using the program MACSE [START_REF] Ranwez | MACSE: Multiple Alignment of Coding SEquences Accounting for Frameshifts and Stop Codons (WJ Murphy[END_REF] and positions that were specific to the differentially expressed genes were recorded. Then, using mapping information, we retrieved the coverage at each position for each of the DE genes in the relevant samples, and checked that positions specific to individual genes were covered at the same level as shared positions.

RESULTS

Variation in gene expression among samples

All samples were sequenced successfully, except one from the 'Musculus Non-Choosy' group ('Musculus Non-Choosy population2 Females', Table S1) due to failure in library construction. For the eleven sequenced samples, the percentage of mapped reads per sample was between 79% and 85%, and the total number of transcripts per sample varied from 33457 to 39770 (Table S2). Table S3 reports the raw and normalised counts for all genes and all samples.

The MDS plots, representing variation in overall gene coverage among samples, indicated a clear separation along dimension 1 between samples of the two subspecies (Figure 1A) and some degree of separation between 'Choosy' and 'Non-Choosy' samples within M. m. musculus (Figure 1B). Multivariate analysis of variance showed that the overall gene expression levels differed significantly between M. m. musculus and M. m. domesticus samples but not between sexes (perMANOVA Adonis test: F subspecies = 6.423; P = 0.004; F sex = 0.804; P = 0.547; F subspecies x sex = 0.322; P = 0.916). The trend for a separation between 'Choosy' and 'Non-Choosy' samples within M. m. musculus was not statistically significant (perMANOVA Adonis test: F choosiness = 2.424; P = 0.089; F sex = 0.782; P = 0.537; F choosiness x sex = 0.506; P = 0.835).

The number of expressed VR genes varied slightly among samples (Table S4a) but remained comparable (range between 197 to 252), with no significant difference between samples from the most divergent groups, i.e. the two subspecies (Wilcoxon rank sum test, P = 0.886).

We also detected expression of Olfactory Receptors (OR) in the vomeronasal organ samples (Tables S3 andS4). The mean levels of expression of VR and OR genes were lower than the mean coverage for the rest of the genes (mean CPM VR = 3.414; mean CPM OR = 0.804; mean CPM other_genes = 27.198), with ORs particularly lowly expressed (Table S4b and Figure S2). In contrast, the mean of median CPM values across all samples was higher for VRs than for the other genes (mean median VR = 0.57; mean median other_genes = 0.29) (Figure S2).

Finally, we report expression of some mouse lipocalins in our VNO samples (Table S3). Seven Odorant Binding Proteins (OBP) were expressed in all samples (including Obp3-p annotated as a pseudogene in the mouse reference genome [START_REF] Stopková | Mouse Lipocalins (MUP, OBP, LCN) Are Coexpressed in Tissues Involved in Chemical Communication[END_REF]). Among the 12 Major Urinary Proteins (MUP) annotated as functional genes and pseudogenes in the reference genome and that belong to Group A (i.e. the ancestral group of MUPs for which mapping information is reliable, [START_REF] Logan | Species Specificity in Major Urinary Proteins by Parallel Evolution[END_REF], we detected expression of three MUPs in all samples (Mup4, Mup5 and to a lesser extent Mup6), and four other MUPs were lowly expressed, mostly in M. m. musculus samples (Mup3, Mup20 (also known as Darcin), Mup-ps19 and Mup-ps20).

Differential expression among samples

Differential expression analyses were performed on the comparison between the two subspecies and the comparison between 'Choosy' and 'Non-Choosy' samples within M. m. musculus. DE genes showed log fold changes from 1 to 10 and average coverage between 1 and 12 logCPM (counts per million) (Figure 2). We found a large number of genes (2624) differentially expressed between M. m. domesticus and M. m. musculus samples (Figure 2A, Table S5a). In contrast, we only found 236 genes differentially expressed between the 'Choosy' and 'Non-Choosy' M. m. musculus samples (Figure 2B, Table S5b). For this comparison, there were more differentially expressed genes up-regulated in the 'Choosy' samples (156 genes) than the 'Non-Choosy' samples (90 genes) (Table S5b), and in general the values of log fold change were also higher in the 'Choosy' samples (Figure 2B).

We identified 90 VR genes and one OR gene (Olfr460) differentially expressed between the two subspecies (Table S5a). We also identified three lipocalins that were differentially expressed between the two taxa (Table S5a): two MUP genes from Group A (Mup6 and Mup20) expressed mostly in M. m. musculus samples, and one OBP gene (Obp6) only expressed in M. m. domesticus samples. More generally, functional analyses revealed an enrichment in 35 functional categories (GO terms) in the set of genes differentially expressed between the two subspecies, with the most enriched clusters of annotation terms corresponding to signal-related functions, vomeronasal receptors of type 2 (V2Rs) and other G-protein-coupled receptors (Table S6a). For the comparison between 'Choosy' and 'Non-Choosy' samples within M. m. musculus, we identified 34 enriched functional categories, with the most enriched clusters of annotation terms corresponding to immune-related functions, signal-related functions and sensory perception of smell (Table S6b). Among the genes found to be differentially expressed between 'Choosy' and 'Non-Choosy' samples, we found five V2R genes (V2r17, V2r69, V2r70, V2r73 and V2r75), which corresponds to a significant enrichment of this category of genes (hypergeometric test: P < 0.01 with 5 out of 140 in a sample of 236 DE genes). Among the five DE V2Rs, four were up-regulated in the 'Choosy' samples and one in the 'Non-Choosy' samples (Figure 3). Differentially expressed genes between 'Choosy' and 'Non-Choosy' samples within M. m. musculus also included two OR genes (Olfr15 and Olfr449) but no members from family 1 of the VR genes (V1R) (Table S5b).

Phylogenetic relationships among differentially expressed vomeronasal receptor genes

When plotting the differentially expressed genes on the vomeronasal receptor phylogenies, we observed that genes differentially expressed between the two subspecies mapped almost everywhere on the V1R and V2R phylogenies (Figure 4A and4B). In contrast, genes differentially expressed between 'Choosy' and 'Non-Choosy' samples only mapped to family 2 of vomeronasal receptors (Figure 4B). One DE VR gene, V2r17, belongs to clade 6 of the V2R family. Interestingly, the other four DE VR genes belong to the same phylogenetic clade, clade 3 (genes V2r69, V2r70, V2r73 andV2r75). The phylogenetic distance between the five DE VR genes was found to be significantly smaller than expected by chance (random sampling test, observed phylogenetic distance = 3.69; p < 0.004 and Figure S3).

The fact that closely related VR genes were found to be differentially expressed between the 'Choosy' and 'Non-Choosy' samples within M. m. musculus could reflect a bias in mapping due to sequence similarity among closely-related paralogs. Indeed, these four VR belong to the same genomic cluster on chromosome 7, and are probably tandem-repeat paralogs with relatively high sequence similarity. However, we found no evidence for such a bias. Although the coverage along the gene was variable, there was no difference between the coverage of conserved and variable regions in any gene, and there were numerous variable sites with high coverage (Figure S4). This was true even when DE genes were compared to their most closely-related paralogs (same phylogenetic cluster). For example, analysis of coverage information along V2r70 (differentially expressed) and its closest relative V2r74 (ENSMUSG00000090774) (not differentially expressed) showed prominent differences in coverage between these two genes for both variable and conserved regions, which indicates that the mapping was able to differentiate reads from different but closely related VR genes.

DISCUSSION

Changes in expression are expected to play an important role in the evolution and divergence of traits involved in reproductive isolation. Here we have used RNA-seq on whole vomeronasal organs collected from wild mice sampled in Denmark to test for the role of changes in expression of vomeronasal receptors in the behavioural shift toward olfactorybased assortative mating observed in the hybrid zone between the two European subspecies of the house mouse. Overall expression levels did not vary significantly between males and females in our dataset, a result in agreement with previous comparisons of male and female mouse VNO transcriptomes [START_REF] Ibarra-Soria | The Olfactory Transcriptomes of Mice[END_REF]. Since assortative mate preferences are displayed by both males and females in the studied mouse populations (Smadja et al. 2004b;[START_REF] Smadja | Asymmetrical reproductive character displacement in the house mouse[END_REF], we did not analyse differential expression between males and females further. Overall expression levels differed between the two subspecies, and we found a very large number of genes differentially expressed between the two taxa (more than 2,000 in total, and 90 VRs), consistent with the degree of divergence between the two taxa [START_REF] Phifer-Rixey | Genome-Wide Patterns of Differentiation Among House Mouse Subspecies[END_REF]. In contrast, only 236 genes were found repeatedly differentially expressed between the M. m. musculus populations showing contrasting patterns of mate preference. Among these DE genes between 'Choosy' and 'Non-Choosy' samples, we found an overrepresentation of genes belonging to family 2 of vomeronasal receptors, known to play a key role in pheromone recognition in mice (Chamero et al. 2012).

Interestingly, four out of five of these differentially expressed VRs belong to the same phylogenetic cluster. In the following paragraphs, we focus on discussing these findings related to expression divergence between populations of M. m. musculus showing contrasting mate preference behaviours.

How confident can we be that differences in expression are associated with differences in behaviour?

Comparative transcriptomics has previously been used to reveal candidate genes for behavioural differences or to reconstruct preference-gene networks (reviewed in [START_REF] Wilkinson | The locus of sexual selection: moving sexual selection studies into the post-genomics era[END_REF]. In our study, how confident can we be that these differences in gene expression are, at least in part, associated with the observed behavioural shift toward assortative mating in the populations nearer to the hybrid zone? First, the use of several biological replicates for each of the two types of natural populations known to differ in their mate preference patterns minimizes the risk of observing population-specific differences that could reflect different regimes of selection or demographic events among populations.

Second, we sampled M. m. musculus populations with contrasting mate preference behaviours but are geographically close in Denmark. The different populations thus share a very recent common ancestor, have similar low rates of introgression from M. m. domesticus and live in similar environments. Our aim was to minimize any other possible natural sources of divergence between these populations and therefore maximize the chance that genetic/expression differences reflect, at least in part, the known behavioural difference between the groups. Since VR expression at the adult stage is known to be influenced by postnatal (from birth to 19 days old) exposure to odours (environment, parental imprinting) at least at the kin recognition level [START_REF] Broad | The Post-Natal Chemosensory Environment Induces Epigenetic Changes in Vomeronasal Receptor Gene Expression and a Bias in Olfactory Preference[END_REF], we decided to study wild-caught mice to ensure assessment of patterns of gene expression as close as possible to natural populations, i.e. potentially influenced by a combination of innate and imprinted mechanisms. However, to homogenise replicates further and avoid any possible variation in gene expression due to differing physiological conditions, we maintained wildcaught mice in the laboratory before dissection to control for the main factors known to influence expression in general and VR expression in particular (diet, temperature, photoperiod, the olfactory environment, age range). Although we could not strictly control for the age of these wild-caught individuals, we selected only mice trapped as young adults expected to fall in the age range corresponding to the peaked plateau of VR expression in adult mice [START_REF] Zhang | High-Throughput Microarray Detection of Vomeronasal Receptor Gene Expression in Rodents[END_REF]. This ensured that age-related variation in VR expression would not strongly influence our results.

Considering that the experiment was performed using wild animals from natural populations, we cannot completely exclude other sources of variation in the observed profiles of gene expression. However, we consider that this series of experimental settings, combined with a targeted approach to candidate genes likely to be involved in the behaviour of interest, offers optimal conditions to detect association between assortative mate preferences and differential gene expression.

Identity and functions of differentially expressed genes

Our main objective was to assess patterns of expression of vomeronasal receptors.

However, we identified several other classes of genes showing repeated patterns of differential expression between 'Choosy' and 'Non-Choosy' samples. Several of these genes have physiological functions related to signal transportation, perception and transduction.

Enrichment analyses showed an overrepresentation of genes involved in regulation of odorant selectivity in the extracellular space and olfactory lumen. This includes genes coding for androgen-binding proteins (ABP), previously suggested to act as subspecific signals in mice [START_REF] Laukaitis | Salivary androgen-binding protein (ABP) mediates sexual isolation in Mus musculus[END_REF][START_REF] Bimova | Reinforcement selection acting on the European house mouse hybrid zone[END_REF], and some organic anion transporters (OAT) expressed in the mouse olfactory mucosa. The latter are potentially involved in odour signalling by binding and transporting organic anions present in urine [START_REF] Kaler | Olfactory mucosa-expressed organic anion transporter, Oat6, manifests high affinity interactions with odorant organic anions[END_REF]).

Other enriched categories include genes involved in regulation of signal transduction (cyclicnucleotide-gated channel (CNG) gene family; [START_REF] Francia | Vomeronasal Receptors and Signal Transduction in the Vomeronasal Organ of Mammals[END_REF]. We also identified two olfactory receptors differentially expressed between the two behavioural groups. Although ORs are primarily expressed in the MOE, our results are congruent with a recent study reporting expression of ORs in mouse VNOs [START_REF] Ibarra-Soria | The Olfactory Transcriptomes of Mice[END_REF]. We did not detect differential expression of MUPs (from Group A) or OBPs between 'Choosy' and 'Non-Choosy' samples, but could confirm previous findings evidencing their expression in VNO tissues [START_REF] Utsumi | Expression of major urinary protein genes in the nasal glands associated with general olfaction[END_REF][START_REF] Ibarra-Soria | The Olfactory Transcriptomes of Mice[END_REF][START_REF] Stopková | Mouse Lipocalins (MUP, OBP, LCN) Are Coexpressed in Tissues Involved in Chemical Communication[END_REF] and expression divergence at some MUP and OBP genes between the two subspecies (although not affecting the same genes as in [START_REF] Stopková | Mouse Lipocalins (MUP, OBP, LCN) Are Coexpressed in Tissues Involved in Chemical Communication[END_REF]. These lipocalins expressed in the VNO may play a role in sequestering volatile pheromones and possibly transporting them to vomeronasal receptors [START_REF] Sharrow | Pheromone binding by polymorphic mouse major urinary proteins[END_REF][START_REF] Stopková | Mouse Lipocalins (MUP, OBP, LCN) Are Coexpressed in Tissues Involved in Chemical Communication[END_REF]. Interestingly, the diversity of genes and functions identified as differentially expressed between 'Choosy' and 'Non-Choosy' samples suggests that the evolution of assortative mating involves a complex network of VNO genes that will need further investigation. Additionnally, genes expressed in other olfactory tissues (MOE) or involved in the integration of the signal in the brain, which we did not study here, may also play an important role in the expression of assortative mate preference.

As far as VR are concerned, we identified five receptor genes (V2r17, V2r69, V2r70, V2r73 and V2r75) that showed repeated differential expression between the two behavioural groups. Since in most cases only one functional receptor is expressed per neuron in the vomeronasal organ (Dalton & Lomvardas 2015;Ishii & Mombaerts 2011), the observed expression divergence at individual VRs could result from differences in transcript level per neuron [START_REF] Young | Odorant receptor expressed sequence tags demonstrate olfactory expression of over 400 genes, extensive alternate splicing and unequal expression levels[END_REF], but also from differences in the numbers of neurons expressing a given receptor, a parameter that shapes the olfactory map in the accessory olfactory bulb and provides to the brain the signal for a behavioural response [START_REF] Mombaerts | Genes and ligands for odorant, vomeronasal and taste receptors[END_REF][START_REF] Bear | The Evolving Neural and Genetic Architecture of Vertebrate Olfaction[END_REF]). Since we used RNA-seq on whole vomeronasal organs, our experiment does not allow us to distinguish these two levels of expression, but is informative about the overall pattern of variation and divergence at individual VRs that likely influence odour perception.

Four of the differentially expressed VR were found to be up-regulated in the 'Choosy' samples and one was up-regulated in the 'Non-Choosy' samples. A majority of VR genes upregulated in the 'Choosy' samples is consistent with the idea that some receptors are recruited in the zone of contact with the other subspecies to allow the expression of assortative mate preferences.

In contrast to differentially expressed VRs between the two subspecies belonging to both families of VRs (family 1 and family 2), expression differentiation within M. m. musculus between 'Choosy' and 'Non-Choosy' samples only affects members of family 2 (V2R). This corresponds to a significant enrichment of this category of genes among all DE genes.

Traditionally, V2Rs were thought to be involved primarily in pheromone recognition (as opposed to V1Rs that were thought to be involved in the detection of the physiological status of other animals). Recent findings have shown that individual V2Rs are specifically tuned to detect the nature of the signal emitter (e.g. conspecifics, heterospecifics, individual recognition), whereas individual V1Rs are activated by cues from multiple species [START_REF] Isogai | Molecular organization of vomeronasal chemoreception[END_REF]. Additionally, V2R-positive Vomeronasal Sensory Neurons detect members of several large peptide or protein families (in particular MHC class I peptides and some Major Urinary Proteins (MUPs) secreted in urine) that are thought to inform the receiver about the genetic identity of a signaler [START_REF] Hurst | Individual recognition in mice mediated by major urinary proteins[END_REF][START_REF] Leinders-Zufall | MHC class I peptides as chemosensory signals in the vomeronasal organ[END_REF][START_REF] Papes | The Vomeronasal Organ Mediates Interspecies Defensive Behaviors through Detection of Protein Pheromone Homologs[END_REF]Chamero et al. 2012). In agreement with this, urines have been shown to diverge in their composition of MUPs between the two subspecies [START_REF] Hurst | Molecular heterogeneity in major urinary proteins of Mus musculus subspecies: potential candidates involved in speciation[END_REF][START_REF] Stopková | Species-specific expression of major urinary proteins in the house mice (Mus musculus musculus and Mus musculus domesticus)[END_REF] and between 'Choosy' versus 'Non-Choosy' samples within M. m. musculus [START_REF] Hurst | Molecular heterogeneity in major urinary proteins of Mus musculus subspecies: potential candidates involved in speciation[END_REF], which establishes MUPs as potential candidates acting as subspecific signals and potentially recognised by V2Rs. Since our study also confirms previous findings that indicate that some MUPs are also expressed in the VNO, a future challenge will be to understand better the respective role of VNO-expressed MUPs and MUPs expressed in the liver and excreted in urine in mouse communication.

Phylogenetic clustering and functional specialisation of differentially expressed vomeronasal receptors

We found that four out of five of the differentially expressed V2Rs belong to the same phylogenetic cluster (V2r69, V2r70, V2r73 and V2r75 in clade 3 of the V2R family). Although these four DE genes are in the same genomic cluster on chromosome 7 and share a relatively high degree of sequence similarity, we have evidence to exclude the possibility of a methodological artefact. We did not find any evidence for a bias when aligning these DE genes with other closely-related paralogs and comparing coverage information in specific versus shared positions along these genes. Moreover, some of the closely-related paralogs belonging to the same genomic cluster as the DE genes did not show significant differential expression. If this pattern of phylogenetic clustering of DE VRs is true, it suggests that the evolution of assortative mate preference in M. m. musculus could involve the specialisation of a group of closely-related receptors in the clade 3 of V2Rs. This hypothesis is consistent with the study by [START_REF] Isogai | Molecular organization of vomeronasal chemoreception[END_REF] which showed that distinct mouse VR receptor subfamilies have evolved towards the specific recognition of certain animal groups (heterospecific signals from predators or sympatric closely-related species/subspecies, male or female conspecific signals) or chemical structures (volatile steroids, non-volatile proteins).

Interestingly, the VNOs analysed in their study were from M. m. domesticus (laboratory mice), and when tested using M. m. musculus scents, several V2Rs were activated including V2r17 (DE gene in our study) and its close relatives and V2r66 (belonging to the same clade 3 as the four VRs differentially expressed in our study). Other receptors identified as differentially expressed in our study could not be individually functionally characterised.

Based on these results, these authors and others (Chamero et al. 2012;[START_REF] Wynn | Genomic variation in the vomeronasal receptor gene repertoires of inbred mice[END_REF] hypothesized a role of these VRs in reproductive isolation. An interesting follow-up from our work would be to apply the methods developed by Isogai et al. to test whether the same set of closely-related V2Rs, identified by these authors as activated by M. m. musculus scents when expressed in M. m. domesticus VNOs, would be activated by M. m. domesticus stimuli when expressed by M. m. musculus VNOs, as our RNA-seq results suggest.

Evolution of expression changes and the role of selection

Regulatory changes in genes involved in reproduction are thought to be prime targets for divergence during speciation. They are expected to evolve under natural and/or sexual selection. Previous studies had shown that expression divergence in testis-associated genes seems to establish during a later phase of the speciation process when reproductive isolation is complete [START_REF] Voolstra | Contrasting evolution of expression differences in the testis between species and subspecies of the house mouse[END_REF]Bryk, Somel, Lorenc, & Teschke, 2013). In our study, we provide evidence for expression divergence in VNO-related genes between populations within the M. m. musculus subspecies, in a context of secondary contact with the subspecies M. m. domesticus. This result may suggest that expression divergence at these genes has evolved rapidly in the 'Choosy' populations from the border of the hybrid zone. Whatever the source of expression divergence (new mutations, standing variation, hybridisation), changes in expression related to the evolution of assortative mate preferences are expected, in the context of reinforcement, to evolve rapidly under diverging selection (natural and sexual selection against hybridisation in the house mouse system, [START_REF] Britton-Davidian | Postzygotic isolation between the two European subspecies of the house mouse: estimates from fertility patterns in wild and laboratory-bred hybrids[END_REF][START_REF] Latour | Sexual selection against natural hybrids may contribute to reinforcement in a house mouse hybrid zone[END_REF]. Although accumulating evidence suggests a predominantly neutral model of gene expression evolution (e.g. in house mice, [START_REF] Staubach | A test of the neutral model of expression change in natural populations of house mouse subspecies[END_REF], positive selection has also been shown to shape changes in expression for some genes and lineages (e.g. [START_REF] Wittkopp | Regulatory changes underlying expression differences within and between Drosophila species[END_REF]Jeukens & Bernatchez, 2012;Chapman et al. 2013).

In a recent study, we identified several genomic regions surrounding VR genes in the house mouse with signatures of positive selection specifically in populations from the border of the hybrid zone that display assortative mate preference (same samples as the 'Choosy' M. m. musculus samples analysed in the present study). Interestingly, one of these genomic regions corresponds to the downstream and putative regulatory region of the cluster containing the candidate VRs identified in the present study, in the vicinity of V2r75 [START_REF] Smadja | Seeking signatures of reinforcement at the genetic level: a hitchhiking mapping and candidate gene approach in the house mouse[END_REF]. Another study on wild-derived strains of the two subspecies has suggested that V2Rs identified by [START_REF] Isogai | Molecular organization of vomeronasal chemoreception[END_REF] as detecting subspecific cues are prone to evolve under positive selection, in contrast to V2Rs detecting other classes of signals [START_REF] Wynn | Genomic variation in the vomeronasal receptor gene repertoires of inbred mice[END_REF]. Analysing genome-wide patterns of nucleotide polymorphism and divergence in the same biological samples as the ones used in the present study should allow us to confirm whether some differentially expressed vomeronasal receptors or their cis-regulatory regions evolve under divergent selection in the populations displaying assortative mate preferences.

To conclude, this study suggests that changes in expression at key receptors for pheromone recognition may play an important role in the evolution of assortative mating and reinforced sexual isolation in the house mouse. Evidence for differential expression at some closelyrelated vomeronasal receptors between populations within M. m. musculus that do or do not display assortative mate preferences (close to or more distant from the hybrid zone, respectively) implicates them in this behavioural shift. This is consistent with previous findings suggesting the specialisation of this group of receptors in the recognition of odorant signals that may allow subspecies recognition and assortative mating.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0

Figure 1 :Figure 2 :

 12 Figure 1: Multi-Dimensional Scaling (MDS) plots indicating variation in gene expression

Figure 4 :

 4 Figure 4: Phylogenetic trees of the two families of vomeronasal receptors (VR) (A: family 1;

Figure S1 :

 S1 Figure S1: Map indicating the geographical origins of the mice used in this study.

Figure S2 :

 S2 Figure S2: Violin plots of log-transformed normalised gene coverage.

Figure S3 :

 S3 Figure S3: Distribution of phylogenetic distances among random sampling of five vomeronasal receptor genes.

Figure S4 :

 S4 Figure S4: Coverage per position (bp) for the five differentially expressed VR genes in the comparison between 'Choosy' and 'Non-Choosy' M. m. musculus samples.

Table S1 :

 S1 Sample description (separated Excel file)

Table S2 :

 S2 General mapping and gene expression statistics (separated Excel file)

Table S3 :

 S3 Raw and normalised (CPM) counts for all genes and all samples (separated Excel

	file, sheet 1: ALL_GENES_raw_counts; sheet 2: ALL_GENES_normalised_counts, sheet 3:
	OR_normalised_counts; sheet 4: VR_normalised_counts; sheet 5: MUP_normalised_counts;
	sheet 6: OBP_normalised_counts)

Table S4 :

 S4 Variation in mean count and in expressed gene number among samples (separated Excel file, sheet 1: Table S4a: number and percentage of expressed genes per sample; sheet 2: Table S4b: average counts per sample)

Table S5 :

 S5 List of differentially expressed genes (separated Excel file, sheet 1: TableS5a:Differentially expressed genes between "Choosy" Mus musculus musculus and "Choosy" M. m. domesticus; sheet 2: TableS5b: Differentially expressed genes between "Choosy" and "Non-Choosy" samples within Mus musculus musculus)

Table S6 :

 S6 Functional enrichment results (separated Excel file, sheet 1: TableS6a: enriched functional categories in differentially expressed genes between 'musculus Choosy' and 'domesticus Choosy' samples; sheet 2: TableS6b: enriched functional categories in differentially expressed genes between 'musculus Choosy' and 'musculus Non-Choosy'

	samples)

ACKNOWLEDGEMENTS

We are very grateful to Ingrid Jakob, former engineer at Centre des Sciences du Goût et de l'Alimentation (UMR 1324), University of Burgundy, Dijon, France, for her advice on dissecting mouse vomeronasal organs. We would also like to thank very much Josette Catalan, Marco Perriat-Sanguinet, Yamin Latour for their help in collecting mice in Denmark and the Danish farmers for their hospitality, as well as Giveskud zoo personnel. We thank Janice Britton-Davidian, Josette Catalan and Marco Perriat-Sanguinet for their help in dissecting mice. We thank the platform "High throughput qPCR" of the University of Montpellier for having given us access to an Agilent instrument to assess RNA sample quality. Finally, we thank Roger Butlin for his comments on a previous version of this manuscript. The research leading to these results has received funding from the European Union's Seventh Framework Programme [FP7/2007-2013] -Marie Curie European Reintegration Grant (ERG), under Grant agreement n°PERG06-GA-2009-251008, as well as from the Agence Nationale pour la Recherche (ANR) under Grant agreement No. 2010BLAN171401-AssortMate. This is publication ISEM 2017-114.

Behavior Genetics, 42, 461-471. Bryk J, Somel M, Lorenc A, Teschke M (2013) Early gene expression divergence between allopatric populations of the house mouse (Mus musculus domesticus). Ecology and evolution, 3, 558-68. Cadiou H, Aoudé I, Tazir B et al. (2014) Postnatal odorant exposure induces peripheral olfactory plasticity at the cellular level. The Journal of neuroscience : the official journal of the Society for Neuroscience, 34, 4857-70. Carroll SB, Grenier JK, Weatherbee SD (2005) From DNA to diversity : molecular genetics and the evolution of animal design. Blackwell Pub. Cassone BJ, Mouline K, Hahn MW et al. (2008) Differential gene expression in incipient species of Anopheles gambiae. Molecular ecology, 17, 2491-504. Chamero P, Leinders-Zufall T, Zufall F (2012) From genes to social communication: molecular sensing by the vomeronasal organ. Trends in Neurosciences, 35, 597-606. Chan YF, Marks ME, Jones FC et al. (2010a) Adaptive Evolution of Pelvic Reduction in Sticklebacks by Recurrent Deletion of a Pitx1 Enhancer. Science, 327, 302-305. Chan YF, Marks ME, Jones FC et al. (2010b) Adaptive Evolution of Pelvic Reduction in Sticklebacks by Recurrent Deletion of a Pitx1 Enhancer. Science, 327, 302-305. Chapman MA, Hiscock SJ, Filatov DA (2013) Genomic Divergence during Speciation Driven by Adaptation to Altitude. Molecular Biology and Evolution. Chung H, Loehlin DW, Dufour HD et al. (2014) A Single Gene Affects Both Ecological Divergence and Mate Choice in Drosophila. Science, 343, 1148-1151. Cucchi T, Auffray J-C, Vigne J-D (2012) History of house mouse synanthropy and dispersal in the near east and europe: a zooarchaeological insight. In: Evolution of the house mouse (eds Macholan M, Baird SJE, P M, Pialek J), pp. 65-93. Cambridge University Press. Dalton RP, Lomvardas S (2015) Chemosensory receptor specificity and regulation. Annual

DATA ACCESSIBILITY

Raw reads have been archived in the NCBI SRA public database under the bioproject identifier PRJNA356041 and through the biosample identifiers SAMN06212561 to SAMN06212563 and SAMN06212565 to SAMN06212572. Expression values (raw and normalised counts) for each gene in each sample are provided in Table S3.

Relative expression per gene

Supplementary information

Discussion:

Robustness of differential expression inferences

RNA-seq on pooled samples is a cost-effective strategy widely used to assess variation in gene expression among experimental conditions or biological groups while producing transcriptome-wide data (reviewed in Todd et al. 2016). This pooling design offers an advantage by increasing sample size into the analysis, and when working with samples that may exhibit a large amount of biological variation, analysing several sample pools, rather than a few single samples per condition, also lessens the impact of single aberrant samples