Analysis of the vomeronasal organ transcriptome reveals variable gene expression depending on age and function in rabbits

Genomics 113 (2021) 2240–2252 Available online 18 May 2021 0888-7543/© 2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Original Article Analysis of the vomeronasal organ transcriptome reveals variable gene expression depending on age and function in rabbits P.R. Villamayor a , b , D. Robledo c , C. Fern´ andez a , J. Gull´ on d , L. Quintela e , P. S´ anchez-Quinteiro b , * , P. Martínez a a Department of Zoology Genetics and Physical Anthropology, Faculty of Veterinary, University of Santiago de Compostela, Lugo, Spain b Department of Anatomy, Animal Production and Clinical Veterinary Sciences, Faculty of Veterinary, University of Santiago de Compostela, Lugo, Spain c The Roslin Institute and Royal (Dick) School of Veterinary Studies, University of Edinburgh, Midlothian, UK d Conejos Gallegos, COGAL SL, Rodeiro, Pontevedra, Spain e Department of Animal Pathology, Faculty of Veterinary, University of Santiago de Compostela, Lugo, Spain ARTICLE INFO Keywords: Vomeronasal organ Gene expression Age Sex VNO-receptors Rabbits RNA-Seq ABSTRACT The vomeronasal organ (VNO) is a chemosensory organ specialized in pheromone detection that shows a broad morphofunctional and genomic diversity among mammals. However, its expression patterns have only been well-characterized in mice. Here, we provide the first comprehensive RNA sequencing study of the rabbit VNO across gender and sexual maturation stages. We characterized the VNO transcriptome, updating the number and expression of the two main vomeronasal receptor families, including 128 V1Rs and 67 V2Rs. Further, we defined the expression of formyl-peptide receptor and transient receptor potential channel families, both known to have specific roles in the VNO. Several sex hormone-related pathways were consistently enriched in the VNO, highlighting the relevance of this organ in reproduction. Moreover, whereas juvenile and adult VNOs showed significant transcriptome differences, male and female did not. Overall, these results contribute to understand the genomic basis of behavioural responses mediated by the VNO in a non-rodent model. 1. Introduction Many animals rely on chemical communication to regulate social and reproductive interactions between conspecifics [1]. This is usually mediated by pheromones, chemical cues mainly detected by two multigenic families, vomeronasal type-1 and type-2 receptors (V1Rs and V2Rs), which are expressed in the neuroepithelium of the vomeronasal organ (VNO) [2,3], a specialized structure located in the nasal cavity and containing vomeronasal sensory neurons (VSNs) [4]. The threshold for detecting some of these chemicals is remarkably low, near 10 −11 M, placing VSNs among the most sensitive chemodetectors in mammals [5]. The VNO displays unique anatomical, physiological, biochemical and genetic features depending on the species [6]. Semiochemical perception and signal amplification and transduction from the VNO to the central circuits of the brain has been mostly studied in rodents [7,8]. However, the molecular and cellular rationale underlying vomeronasal function is not fully understood, likely due to the wide range of molecules detectable by the VNO –from proteins or peptids to steroids, major Abbreviations: 2MB2, 2-methyl-but-2-enal; Aas, Amino acids; Ano, Anoctamin; AOB, Accesory olfactory bulb; Arg, Arginase; DE, Differential expression; DEGs, Differentially expressed genes; Dlk1, Delta like non-canonical Notch ligand 1; Dll4, Delta-like 4 signaling gene; Esr, Estrogen receptor; FC, Fold change; FDR, False discovery rate; FPRs, Formyl peptide receptors; Fshr, Follicle-stimulating hormone receptor; GnRH, Gonadotropin releasing hormone; GO, Gene ontology; MHC, Major Histocompatibility Complex; KEGG, Kyoto Encyclopedia of Genes and Genomes; LCN, Lipocalins; Lhcgr, Luteinizing hormone/choriogonadotropin receptor; MOE, Main olfactory epithelium; Mup4, Major urinary protein 4; Mups, Major urinary proteins; M&M, Material and Methods section; NCBI, National Center for Biotechnology Information; Oxtr, Oxytocin receptor; PCA, Principal component analyses; PE, Paired-end; Pgr, Progesterone receptor; Pgrmc, Progesterone receptors membrane component; RIN, RNA integrity number; RNA-Seq, RNA sequencing; RT-PCR, Reverse transcription polymerase chain reaction; SBT, 2-s-butyl-4,5- dihydrothiazole; SRA, Short Read Archive; TPM, Transcripts Per Million; Trp, Transient receptor potencial channel; Trpc2, Transient receptor potential cannel 2; V1Rs, Vomeronasal receptors type 1; V2Rs, Vomeronasal receptors type 2; VNO, Vomeronasal organ; VRs, Vomeronasal receptors; VSNs, Vomeronasal sensory neurons. * Corresponding author. E-mail address: [email protected] (P. S´ anchez-Quinteiro). Contents lists available at ScienceDirect Genomics journal homepage: www.elsevier.com/locate/ygeno https://doi.org/10.1016/j.ygeno.2021.05.007 Received 8 December 2020; Received in revised form 23 April 2021; Accepted 14 May 2021 Genomics 113 (2021) 2240–2252 2241 histocompatibility complex (MHC) molecules and some small metabolites [9,10,11,12]–, and also to the diverse family receptor repertoire found in VSNs (V1R, V2R, formyl-peptide and MHC receptors) [13,14,15,16]. Additionally, the number of V1Rs and V2Rs genes greatly varies across mammalian genomes, from non-functional genes in macaques to several hundred in rodents [17,18,19,20]. The genomic characteristics of the VNO have been studied through comparative genomics, comparing genome assemblies of different species [21], and through gene expression of specific vomeronasal receptor (VR) genes or subsets of genes by reverse transcription polymerase chain reaction (RT-PCR), in situ hybridization and microarrays [22,23,24,25]. However, only a few studies have addressed the analysis of the whole VNO transcriptome in a handful of species (mouse, zebrafish and bats; [26,27,28]). In particular, the mice vomeronasal organ has been studied under different experimental conditions, showing variation in vomeronasal receptors expression between strains [29], with pregnancy [30] and under sex-separation [31], suggesting that VNO features may vary not only among species but also under different conditions within the same species. Rabbits are a species pertaining to the order Lagomorpha, phylogenetically close to Rodentia. However, since pheromones are speciesspecific signals, each species should be considered independently according to its reproductive pattern and behavioural priorities [32]. Indeed, the rabbit is the only mammalian species in which a mammary pheromone −2-methyl-but-2-enal (2 MB2)- has been detected [33], becoming one of the best studied models of pheromonal communication in mammals [34]. Additionally, rabbits are a farm species and more recently they have also become a common pet, being the third preferred pet worldwide after dogs and cats [35]; thus understanding pheromone perception by the olfactory subsystems may contribute towards the implementation of pheromone-based therapies for improving both animal production and welfare. In fact, the Rabbit Appeasing Pheromone –2 MB2 analog– has been commercialized as a method to reduce stress, increase reproductive efficiency and improve animal welfare [36], and some behavioural data regarding the rabbit mammary pheromone have been collected [37,38,39]. These studies, along with the comprehensive study of the vomeronasal system at anatomical level [40,41], highlight the relevance of chemocommunication in this species. The rabbit VNO was reported to be highly developed with many specific morphological features [40], but surprisingly it has not been considered in most vomeronasal phylogenetic studies. Indeed, Grus et al. [21] showed that the complete repertoire of rabbit vomeronasal receptors has still to be determined and to date 160 V1Rs have been reported in the rabbit genome [18] and 37 V2Rs have been depicted in a phylogenetic tree of mammals by Francia et al. [13]. Additionally, despite the well-characterized VR subfamilies in mice [42], no data regarding VR subfamilies and clades in rabbits have been found. Molecular-genetic approaches have greatly contributed to increase our understanding of vertebrate pheromonal communication [32], but to our knowledge the rabbit VNO has not yet been studied. Here, we carried out an RNA sequencing (RNA-Seq) approach to: i) characterize the rabbit VNO transcriptome; ii) define the vomeronasal receptors expressed in the rabbit VNO iii) identify vomeronasal-specific genes by comparing the vomeronasal transcripts repertoire to other rabbit tissues; and iv) evaluate differences in gene expression both between males vs females and juveniles vs adults. We obtained a broad outlook of the rabbit vomeronasal gene repertoire, identifying VNO-specific genes, and observing a considerably fluctuation in their overall expression profile, especially depending on age but not on gender, despite the enrichment of several sex hormone-related pathways supporting the pivotal role of the VNO in reproductive behaviour. Altogether, these results suggest a great flexibility of the vomeronasal expression in rabbits and may contribute to understand the diverse physiological mechanisms underlying the vomeronasal organ function. 2. Material & methods 2.1. Experimental design and sampling We employed 24 animals, separated in 12 juveniles and 12 adults (40 and 180 days, respectively) including the same number of males and females within each group. All animals pertained to a commercial hybrid -Hyplus strains PS19 and PS40 for female and male, respectively- and were maintained on a farm (COGAL SL, Rodeiro, Spain) under the same temperature conditions (18–24 ◦C), dark-light cycles of 12:12 h and ad libitum drinking and feeding. All individuals were humanely sacrificed by an abattoir of the same company, in accordance with the current legislation, and their heads were separated from the carcasses in the slaughtering line. The whole VNOs were immediately dissected out after opening the lateral walls of the nasal cavity and removing the nasal turbinates. The hidden location of the VNOs in rabbits under the mucosal covering the ventral part of the nasal septum, and their anatomical structure have been previously described in Villamayor et al. [40]. Previous experience was critical to obtain the organs in a short period of time (<5 min after death) essential for obtaining the appropriate RNA quality (see below). The two VNOs of each animal were collected and immediately stored in Trizol and kept in ice (~4 ◦C). Due to the double bone and cartilage envelope of the rabbit VNO, the tissue was homogenized using a mixer to guarantee the whole tissue sample was soaked by Trizol. After 20 min, all samples were stored at −80 ◦C for further RNA extraction. 2.2. Transcriptomic analysis RNA extraction was performed using Trizol followed by DNase treatment in accordance with the manufacturer’s instructions. RNA quantity and integrity were evaluated in a NanoDrop® ND-1000 spectrophotometer (NanoDrop® Technologies Inc.) and in a 2100 Bioanalyzer (Agilent Technologies), respectively. The 24 RNA samples displayed RNA integrity number (RIN) values >7.6, so appropriate for library construction and sequencing. Samples were barcoded and prepared for sequencing by Novogene (Cambridge, UK) on an Illumina Nova-Seq 6000 150 bp PE run. The quality of the sequencing output was assessed using FastQC v.0.11.7 (https://www.bioinformatics.babraham.ac.uk/projects/fast qc/). Quality filtering and removal of residual adaptor sequences were performed on read pairs using Trimmomatic v.0.39 [43]. Specifically, residual Illumina-specific adaptors were clipped from the reads, the read trimmed if a sliding window average Phred score over five bases was <20, and only reads where both pairs had a length longer than 50 bp post-filtering were retained. Filtered reads were aligned against the rabbit genome (OryCun2.0; [44]) and assigned to genes based on the latest annotation of the rabbit genome [45] using STAR v.2.7.0e [46] two-pass mode and the following parameters: the maximum number of mismatches for each read pair was set to 10% of trimmed read length, and minimum and maximum intron lengths were set in 20 bases and 1 Mb, respectively. Gene count data were used to calculate gene expression and estimate differential expression (DE) using the Bioconductor package DESeq2 v.1.28.1 [47] in R v.3.6.2 (R [48]). Size factors were calculated for each sample using the ‘median of ratios’ method and count data were normalized to account for differences in library depth. Normalized reads were used as a measure of expression and were calculated by taking the average of the normalized count values, dividing by size factors, taken over all samples. This corresponds to the ‘basemean’ obtained with DESeq2. Next, gene-wise dispersion estimates were fitted to the mean intensity using a parametric model and reduced towards the expected dispersion values. Finally, differential gene expression was evaluated using a negative binomial model that was fitted for each gene, and the significance of the coefficients was assessed using the Wald test. The Benjamini-Hochberg false discovery rate (FDR) correction for multiple P.R. Villamayor et al. Genomics 113 (2021) 2240–2252 2242 tests was applied, and transcripts with FDR <0.05 were considered differentially expressed genes (DEGs). Hierarchical clustering and principal component analyses (PCA) were performed to assess the clustering of the samples and identify potential outliers over the general gene expression background. The R packages “pheatmap”, “PCAtools” and “EnhancedVolcano” were used to plot heatmaps, principal component analysis (PCA) and volcano plots, respectively. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the David 6.8 Bioinformatic Database [49,50] for all conditions. Additionally, to deal with STAR outputs multi-mapping reads, we further employed kallisto [51], which pseudoaligns reads to a reference, producing a list of transcripts that are compatible with each read while avoiding alignment of individual bases. We compared the results obtained by both software, STAR and kallisto, especially regarding vomeronasal receptor genes, as they are highly prone to multi-mapping due to sequence homology. Note that STAR was used by default in our analysis and therefore it is only specifically mentioned when such data is compared to the kallisto output. Several protein coding genes in the rabbit genome were not annotated in Ensembl, and were blasted against the protein database Swiss-Prot ([52]; E-value <10 −5 ). The new identified VR genes were further validated against Pfam 34.0 model using HMMER [53]. Some specific genes belonging to formyl peptide (FPRs) and transient receptor potential channel (Trp) families, known to play an important role in the mice VNO, were not found in the OryCun 2.0; thus, mice sequences were carefully explored by blasting them against the latest version of the rabbit genome OryCun 3.0 (https://www.ncbi.nlm.nih. gov/assembly/GCA_013371645.1, 2020). It should be noted that we did not use this genome version as the reference in our study because it has not been annotated and many of the transcripts have indels (insertions / deletions) that shift the reading frame. After personal communication (Leif Anderson, Upssala University), authors suggested to use OryCun2.0 meanwhile OryCun4.0 is not available, despite Ory- Cun 3.0 is a useful tool to explore gene content in the rabbit genome. Finally, we also blasted our data against the mice transcriptome (Mus musculus, assembly GRCm39) (Genome Reference Consortium, https: //www.ncbi.nlm.nih.gov/assembly/GCF_000001635.27/, 2020), employing kallisto, in order to detect potential unidentified genes in the rabbit genome. 2.3. Vomeronasal receptors. Comparison among different rabbit transcriptome tissues and phylogenetic evaluation Vomeronasal receptor gene expression was compared to that of other seven rabbit tissues (hindbrain, forebrain, ovary, testis, liver, heart and kidney) available in the Rabbit Expression Atlas (https://www.ebi.ac. uk/gxa/experiments/E-MTAB-6782/Supplementary%20Information; [54]). Genes of the VNO were considered functionally relevant for the VNO when at least more than 90% of the VNO samples analyzed (measured as Transcripts Per Million, TPM) doubled the expression in more than 98% of the atlas tissue samples. The evolutionary history of the rabbit VR repertoire was evaluated against the mouse VR repertoire. OrthoFinder [55,56] was used to infer orthogroups and gene trees for all VR genes found in the rabbit and the mouse genome. The R package ggtree [57] was used to draw the phylogenetic tree representing the evolutionary history of the rabbit and mouse VR1 and VR2 repertoires. 2.4. Data availability The RNA-seq raw reads analyzed in the current study are available on the National Center for Biotechnology Information (NCBI) Short Read Archive (SRA) database under accession number PRJNA720622. 3. Results 3.1. RNA sequencing output A total of ~843 million paired-end (PE) reads were generated, for an average of ~35 million reads for each of the 24 samples analyzed. After filtering, ~798 million reads were retained (94.7%; >33 million reads per sample). On average, 76.6% of the PE filtered reads were aligned to the available rabbit genome (OryCun2.0; [44]) at unique (70.0%) or multiple (6.6%) genomic positions. These reads were assigned to the 29,587 genes identified in the rabbit genome. A total of 7690 out of the 22,675 protein coding genes in the rabbit genome were not annotated in Ensembl, and were blasted against the protein database Swiss-Prot ([52]; E-value <10 −5 ), being annotated 2670 additional genes. 3.2. Vomeronasal transcriptome Given the lack of previous molecular and genomic knowledge about the rabbit VNO, we first characterized its global transcriptome. A total of 19,482 genes (65.8%) with at least one normalized read were identified, and 14,584 (49.3%) presented ≥20 normalized reads (Table 1); this is the threshold that we established to consider a gene as expressed based on the distribution of the number of normalized reads per gene (Supplementary file 1). GO term and KEGG pathway enrichment analyses of the expressed genes (≥20 normalized reads) (Supplementary file 2) revealed the overrepresentation of several interesting functional categories such as PI3K- Akt signaling, an important pathway for survival and proliferation of VSNs [58], and MAPK, involved in the activation of VSNs [59]. Importantly, the reproductive function of the VNO was represented with terms such as gonadotropin releasing hormone (GnRH) signaling and estrogen signaling pathway (p <0.05).To refine the analyses of the rabbit VNO transcriptome constitution, we considered different gene expression ranges: i) low expression (20–100 normalized reads), ii) moderate expression (100–1000 normalized reads) and iii) high expression (≥ 1000 normalized reads) (Table 1). Additionally, since STAR does not account for multi-mapping reads, which might lead to biased gene expression estimates, we compared the results obtained by STAR with those from kallisto (Supplementary files 3 and 4), which accounts for multi-mapping reads. Both software provided similar results, especially for the higher expression categories (Table 1). Functional enrichment in the three different subgroups (Supplementary file 2) showed the enrichment of the KEGG pathway terms neuroactive ligand-receptor interaction and calcium signaling pathway in the 20–100 list, which fits with the low expression needed by VRs to trigger signal transduction. In VSNs, such signal is mediated by G-pro- tein coupled mechanisms, and accordingly, terms such as heterotrimeric G-protein complex or G-protein coupled receptor signaling pathway were also found in the GO 20–100 list (p <0.05). Similarly, other terms such as presynaptic membrane, neuronal cell body membrane and glutamatergic synapse, may underlie the ‘basic needs’ of the VSNs. Moderate and high expression genes were instead mainly enriched in reproduction pathways such as circadian rhythm (involved in animals seasonal reproduction), GnRH, estrogen, oxytocin, and prolactin Table 1 Gene expression in the rabbit vomeronasal organ. Range of expression STAR software Kallisto software Genes Genes Whole transcriptome 29,587 21,268 ≥1 19,482 17,300 ≥20 14,584 13,976 Low expression (20−100) 3019 2725 Moderate expression (100−1000) 8573 8559 High expression (≥1000) 2601 2692 P.R. Villamayor et al. Genomics 113 (2021) 2240–2252 2243 signaling pathways, together with p53 signaling pathway, whose loss leads to a significant decrease of fertility [60] (p <0.05) (Supplementary file 2). This is consistent with the wide expression of reproductionrelated genes found in the rabbit VNO, such as nine estrogen-related genes and eight progesterone-related genes (Supplementary file 5). The case of GnRH is especially interesting since several genes controlling GnRH production were moderately expressed in the rabbit VNO such as Gli3, related to GnRH-1 migration to the central nervous system [61]. Dmxl2 and Lgr4, whose deficiency results in decreased fertility or delayed puberty due to fewer GnRH neurons, were also moderately expressed in our samples [62,63]. Recently, Taroc et al. [64] identified the transcription factor Isl1 sparsely expressed in postnatal VNO neurons, but they suggested a dispensable role for Isl1 in GnRH neurons. Accordingly, we found this gene lowly expressed in rabbits. However, the two GnRH receptors (GNRH1 and GNRHR) showed very low expression (1–10 normalized reads), below the established expression threshold (see Material and Methods section, M&M). Additionally, three prolactin-related genes were moderately or highly expressed, and 15 genes related to androgens were also moderately expressed on average. Conversely, oxytocin receptor (Oxtr), luteinizing hormone/choriogonadotropin receptor (Lhcgr) and folliclestimulating hormone receptor (Fshr) showed lower expression or were not detected. We also found 34 genes related to spermatogenesis, but their role in the VNO remains unknown. Finally, we detected 11 genes involved in circadian rhythm (Supplementary file 5). 3.3. Vomeronasal receptors We were especially interested in understanding the particular expression of the two main vomeronasal receptor multigene families (V1Rs and V2Rs). To date, 141 VRs (125 V1Rs and 16 V2Rs) have been found in the rabbit genome (OryCun 2.0). Here, by assessing the nonannotated protein coding genes of our rabbit VNO transcriptome (see M&M), we were able to identify 58 new VRs (4 V1Rs and 54 V2Rs). These new genes were further validated against Pfam 34.0 model using HMMER [53], and 54 out of 58 genes (3 V1R and 51 V2R) showed protein domains consistent with vomeronasal receptors (Supplementary file 6). In sum, we identified a total of 195 VR genes in the rabbit genome, which is towards the upper range among mammalian VR gene repertoires (from >250 in mice to ~8 in dogs), including 128 V1Rs and 67 V2Rs genes (Supplementary file 4 and 7). The latter figure is rather surprising because we identified up to 51 new V2Rs, a family nearly degenerated in most mammals. Additionally, following the same criteria as [65], we classified the total VR gene repertoire into 175 intact VRs (encoding >300 amino acids (aas)) and 20 partial VRs (>100 aas and < 300 aas). Most vomeronasal receptor genes showed relatively low expression, especially those pertaining to the family V1Rs (Table 2). Nonetheless, a large dynamic range of expression was observed in the VRs identified (Fig. 1). The highest expressed gene was VN2R1 (ENSOCUG00000002314) with 308.20 normalized reads, while 14 V2R genes showed moderate expression levels (100–1000 normalized reads). However, most V2Rs were lowly expressed (20–100 normalized reads), and the majority of V1Rs genes showed less than 20 normalized reads on average (Supplementary file 7). The mean expression for V1Rs was 11.56 ±1.14 whereas for V2Rs was 56.84 ±8.04. This is not surprising because most V1Rs are thought to be expressed only in one or two specific subsets of neurons, being tuned up by a cognate ligand, while at least some V2Rs can be co-expressed in the same neurons [66]. Accordingly, since each neuron expresses different VRs, it is reasonable that the expression of each individual VR in the whole VNO is relatively low. To explore the evolution of the rabbit VR repertoire, we performed a phylogenetic analysis (Fig. 2) comparing the Ensembl VR repertoire of the rabbit (Oryctolagus cuniculus, OryCun 2.0), including our newly identified VR, with that of the mouse (Mus musculus, GRCm39) –pseudogenes discarded–. We found species-specific expansions of ancestral genes in the rabbit VR repertoire which are generally distributed in a few rabbit-specific clades. Additionally, those VR clades seem to cluster at specific genomic regions, suggesting tandem or local duplications. For example, we found 30 V1R rabbit genes in chromosome 9. Regrettably, most rabbit VR genes were located in scaffolds, and although many of them are also grouped, a new rabbit genome assembly will be necessary to confidently study the genomic architecture and evolutionary history of this gene family. A third class of VNO receptors includes 5 out of the 7 members of the formyl peptide receptor family -Fpr-rs1 (also known as Fpr3), Fpr-rs3, Fpr-rs4, Fpr-rs6 and Fpr-rs7, previously reported in mice [67,16]. We did not find these genes in the OryCun 2.0 or in OryCun 3.0 rabbit genome assemblies (see M&M), instead all mice VNO Fprs sequences matched to a single genomic location in the rabbit genome (OryCun 3.0), which corresponds to the rabbit FPR2 gene, thus suggesting that mice VNO Fprs might represent an expansion of the ancestral FPR2 ortholog. Furthermore, we aligned our rabbit VNO RNA-Seq against the last version of the mice transcriptome (GRCm39) and none of the reads matched against the mice VNO Fpr genes. On the other hand, we found that not only FPR2 but also FPR1 are expressed at low levels in the rabbit VNO. While expression of these receptors in the mouse vomeronasal neurons has not been described [16], they are expressed in other cell types and tissues more related to immunity [68], and this is consistent with the presence of lymphoid tissue in the VNO [69]. 3.4. Vomeronasal-specific gene expression To assess the expression pattern of the genes found in the rabbit vomeronasal organ in a broader scenario, we compared their expression to those reported in hindbrain, forebrain, ovary, testis, liver, heart and kidney using publicly available RNA-Seq data [54] (Supplementary file 8). We found 429 genes predominantly expressed only in the VNO. Among them, 80 vomeronasal receptors (45 V1Rs and 35 V2Rs) were identified as VNO-specific, but other 71 VRs (47 V1Rs and 24 V2Rs; also considering those which did not match against the Pfam model) were not present in the rabbit expression atlas, and therefore their VNO- specificity could not be assessed. Additionally, some VRs VNO- specific, especially V1Rs, showed high expression in the adult testis (Supplementary Figs. 1 and 2); this is not surprising considering the pervasive expression described in this tissue [70], and does not imply that these genes have a functional role in the testis. We analyzed the expression of the remaining 48 ‘non VNO-specific’ VRs (37 V1Rs and 11 V2Rs), and interestingly their overall expression was higher in the VNO than in any other analyzed tissue (Supplementary file 8; Supplementary Figs. 3 and 4). However, in some VNO samples their expression was 0, and in general the low expression of these genes explains why they were not detected as VNO-specific. Additionally, some VRs were slightly expressed in other tissues, especially in testis, ovary and brain. No expression of VR genes was found in heart, liver, and kidney (Supplementary Figs. 1-4). The transient receptor potential channel 2 gene (Trpc2) was also exclusively found in the VNO. Additionally, six other genes of the MHC, Table 2 VR gene expression in the rabbit vomeronasal organ. Range of expression STAR software Kallisto software VR V1R V2R VR V1R V2R Whole transcriptome 195 128 67 195 128 67 ≥1 176 112 64 176 114 62 ≥20 66 22 44 54 20 33 Low expression (20–100) 52 22 30 46 20 26 Moderate expression (100–1000) 14 0 14 8 0 8 High expression (≥1000) 0 0 0 0 0 0 P.R. Villamayor et al. Genomics 113 (2021) 2240–2252 2244 either from class I or class II, as well as several immunoglobulins and other immunity-related genes such as lipocalins (LCN) were predominantly expressed in the VNO. In more detail, the major urinary protein 4 (Mup4) and two lipocalin genes (ENSOCUG00000026763, ENSOCUG00000028189; both of them coded as LCN1) were among the most expressed genes in the VNO (>164,191 normalized reads) and considered VNO-specific, thus suggesting the importance of these proteins in the VNO function. Looking at the molecular function, the term pheromone receptor activity was enriched in vomeronasal-specific genes, thus supporting pheromone-detection as the VNO principal role. Additionally, some biological process and cellular component terms were related to the immune system (i.e. immune response, MHC class I protein complex, MHC class II protein complex), as expected considering the wide range of immune-related genes found in the list (Supplementary file 8). This is consistent with previous studies where genes from immune response and chemosensory receptor classes were found as the most represented genes in the VNO [29]. 3.5. Differential expression of the vomeronasal organ between females and males We studied the transcriptomic differences in the VNO between animals of different gender, aiming to discover the functional basis of the reproductive role of the VNO in rabbits. The vomeronasal organ transcriptome showed very small differences between females and males as shown in the PCA (Fig. 3A) and volcano plot (Fig. 3B). Only 12 significant DEGs were detected, 6 more expressed in females and 6 more in males (Supplementary file 9), none of them pertaining to the VR families or representing a relevant role in reproductive-hormonal circuitries or sexual dimorphism. 3.6. Differential expression of the vomeronasal organ between juveniles and adults Similarly, we studied the difference between the vomeronasal transcriptomes of juvenile and adult rabbits. The vomeronasal organ transcriptome showed sharp differences between juveniles and adults, and the two groups clustered separately in the PCA (Fig. 4A). A total of 3061 DEGs were detected between juveniles and adults (Fig. 4B), 1376 upregulated (higher expression in adults) and 1685 down-regulated (higher expression in juveniles) (Supplementary file 10). Some genes of the VR gene repertoire were more expressed in juveniles, including six and one genes of the V2R and V1R families, respectively. The remaining VR genes did not show differential expression, but many of them displayed log2 FC either >1 or <1, thus suggesting that their lack of significance might be due to their low expression levels or sampling error, and in fact it was not possible to calculate the Benjamini-Hochberg pvalue correction for 70 VR genes (see Supplementary file 10, DE_VR and DE_VR_all). In fact, the seven VR DEGs were among the most expressed VR genes, all with more than 20 normalized reads. This almost doubles the mean expression of V1Rs and likely explains why only one V1R was detected as differentially expressed. Additionally, a few reproduction-related genes (i.e. progesterone receptor Pgr, and estrogen receptor Strn) and others related to vomeronasal function such as Sema3F, Kirrel2, Anoctamin 2 (Ano2) and Fig. 1. Heatmap of vomeronasal receptor genes. Heatmap showing the patterns of expression of vomeronasal receptors (V1Rs and V2Rs) in the VNO of juvenile and adult rabbits. Gender has not been included since no significant DE genes were detected between males and females. P.R. Villamayor et al. Genomics 113 (2021) 2240–2252 2245 Robo2 were either up- or down-regulated. Remarkably, Trpc2, essential in VNO signal transduction, was more expressed in adults, whereas Trpc1 was more expressed in juveniles. Other groups such as arginase (Arg1, Arg2) or Notch-related genes (Delta-like 4 signaling gene (Dll4), Delta like non-canonical Notch ligand 1 (Dlk1)) also showed differences between ages (Supplementary files 10 and 5). Overall, a set of 32 genes related to vomeronasal activity and function and showing highly significant differences between juveniles and adults enabled the classification of all individuals by age in our experiment with full confidence (Fig. 5). GO term and KEGG pathway enrichment analyses were performed for both lists of up-regulated (more expressed in adults) and downregulated (more expressed in juveniles) genes separately (Supplementary file 11). Some reproductive-related terms such as estrogen, prolactin and GnRH signaling pathways were found in any of the two lists (p value <0.05) (Supplementary file 9). In juveniles, terms such as the previously mentioned MAPK and PI3K-Akt signaling, as well as p53 signaling, relevant to ensure faithful development and reproduction [71], were significantly enriched (p value <0.05). Instead, adults showed the enriched term myelin sheath, which together with other nominally significant terms such as receptor activity, neuronal cell body and myelination in peripheral nervous system (p value <0.05), fit the expected functional role of the adult vomeronasal organ. All in all, the VNO transcriptomes of juveniles and adults show notable differences, and the results suggest that the VNO is involved in modulating reproductive functions from early life-stages at least from day 40 after birth, just post-weaning. 4. Discussion This study provides a comprehensive gene expression analysis of the rabbit vomeronasal organ using RNA-Seq, characterizing for the first time the rabbit VNO transcriptome. Remarkably, the two main VR families (V1Rs and V2Rs) showed very low expression, especially V1Rs. Even though our study addresses the whole VNO and not specific cell types such as VSNs, the low expression of VRs as well as their wellknown expression in VSNs would be consistent with the hypothesized expression of each V1R in a single sensory neuron in mice [72,73]. We also identified genes significantly more expressed in the VNO than in any other tissue in the rabbit gene expression atlas [54]. While the VNO transcriptomes of juveniles and adults have been compared before in mice [74], the study was focused on lncRNA and not on VR expression. Our study is the first to compare VR expression in the VNO between adults and juveniles in mammals, and has revealed considerable fluctuation in the vomeronasal repertoire depending on age. Finally, whereas male and female VNO transcriptomes did not show significant differences, we found several sex hormone-related pathways enriched in the overall rabbit VNO, confirming previous results in mice. 4.1. The unique gene repertoire of the rabbit vomeronasal organ The vomeronasal organ acts as an interface between the immune and nervous systems and carries out many functions such as pheromone Fig. 2. Phylogenetic tree of the rabbit and mouse V1R and V2R. Fig. 3. Comparison between female and male VNO transcriptomes. A) Principal Components Analysis (PCA) showing the clustering of the rabbit VNO RNA-Seq samples, with the samples colored according to their sex. B) Volcano plot showing the differential expression test between female and male VNO RNA-Seq samples. Each point in the plot represents a gene, with its log2 fold change (FC) in the x-axis and its log10 p-value in the y-axis. Genes are classified in four categories depending on their FC and FDR corrected p-value: i) grey =p-value >0.01 and log2 FC between −0.5 and 0.5; ii) green =p-value >0.01 and log2 FC <−0.5 or >0.5; iii) blue =p-value <0.01 and log2 FC between −0.5 and 0.5; and iv) red =p-value <0.01 and log2 FC <−0.5 or >0.5). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) P.R. Villamayor et al. Genomics 113 (2021) 2240–2252 2246 perception and signal transduction, receptor activity and neuron maintenance, allowing the regulation of a wide range of behaviours, including individual recognition, reproduction, aggression [75,76,77,13] and sick conspecific avoidance behaviour [78]. This is mediated by highly qualified sensory receptors expressed in the neuroepithelium of the VNO, whose high specificity results in very low expression. Due to the various functions of the VNO, it is reasonable that its transcriptome shows wide intrinsic variability and fluctuates over time depending on the function at each life-stage and social interaction. a) Vomeronasal-type receptors The rapid evolution of the two main vomeronasal receptor superfamilies -V1Rs and V2Rs- in vertebrate lineages places them among the gene families showing the broadest variation in gene number [79,80]. For instance, snakes and lizards exhibit a large number of V2R genes, but have retained an extremely limited number of V1Rs [81]. Instead, mammals tend to have a highly developed VNO with an expanded V1R repertoire, as it occurs in platypus (283 genes), mouse (239), rat (109) and rabbit (here updated to 128). A total of 121 and 37 intact V2Rs (here updated to 67) were found in mouse and rabbit, respectively, whereas no functional V2R genes have been reported in species such as Old World monkeys or humans [13]. The loss of these receptors in this lineage is confirmed by the high number of pseudogenes and a vestigial and presumably non-functional VNO [82], even though the latter is still a matter of controversy [83,84,85]. Additionally, VRs are very dynamic genes with rapid rates of gene duplication, pseudogenizations and lineage-specific expansions [86], mostly studied in rodents [87,88]. We found that rabbits have undergone several expansions from ancestral genes that are not shared with mice. Consequently, the vomeronasal gene repertoire at both genomic and transcriptomic levels must be considered in each species independently. The VR repertoire of a species may also vary in different conditions. In rabbits, seven VRs were more expressed in juveniles than in adults. Although we expected greater differences, the low expression of these genes makes it difficult to reach significance and accordingly, more important biological differences in VR gene expression between juveniles and adults cannot be ruled out. Moreover, VRs tended to be more expressed in juveniles in our study, even those that did not reach statistical significance. Considering that the VNO has been long known to elicit innate or unconditioned responses to specific chemical signals [89], these results might reflect the continuous exposure to new stimuli during the first stages of life, leading to a full activation of the organ in 40 days old animals. To our knowledge, this is the first time the expression of VRs between juveniles and adults has been studied in mammals, and more studies focused on the VR repertoire are needed to clarify the role of these receptors at different life-stages, not only in rabbits but also in other species. Fig. 4. Comparison between juvenile and adult VNO transcriptomes. A) Principal Components Analysis (PCA) showing the clustering of the rabbit VNO RNA-Seq samples, with the samples colored according to their age. B) Volcano plot showing the differential expression between juvenile and adult VNO RNA-Seq samples. Each point in the plot represents a gene, with its log2 FC in the x-axis and its log10 p-value in the y-axis. Genes are classified in 4 categories depending on their FC and FDR corrected p-value: i) grey =p-value >0.01 and log2 FC between −0.5 and 0.5; ii) green =p-value >0.01 and log2 FC <−0.5 or >0.5; iii) blue =p-value <0.01 and log2 FC between −0.5 and 0.5; and iv) red =p-value <0.01 and log2 FC <−0.5 or >0.5). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 5. Heatmap of genes of interest showing differential expression between juvenile and adult VNO. Heatmap showing the expression pattern of genes differentially expressed between juvenile and adult VNO RNA-Seq samples and with putative VNO-related function. P.R. Villamayor et al. Genomics 113 (2021) 2240–2252 2247 On the other hand, a comparison of the rabbit VR gene expression among different tissues has revealed 80 VRs VNO-specific and 48 VRs ‘non VNO-specific’. The latter were not considered VNO-specific following the criteria employed in the analysis rather than biological reasons, and in both cases VR genes showed overall higher expression in the VNO than in other tissues. This is consistent with previous studies in mice, where the vast majority of VR genes were only expressed in the VNO [25]. However, we found some VNO-specific VRs, especially V1Rs, also expressed in adult testis, which is consistent with previous studies in mice [90] and swine [65]. Some VRs, mainly those belonging to the ‘non VNO-specific’ group, were also found lowly expressed in both gonads and brain. Expression in brain and bulb was previously detected in mice [25]. No other tissues in any species have shown VR expression apart from the V1Rs expressed in the main olfactory epithelium (MOE) of goats, mice, lemur (also showed VN2R2 expression) and humans [91,82,92]. Interestingly, Ibarra-Soria et al. [26] found one V2R, Vmn2r29, expressed in MOE, being its expression higher than the median olfactory receptor expression, suggesting previous unrecognized mechanisms of pheromone detection in the MOE. In rabbits, the MOE transcriptome has not yet been assayed, and further studies are needed to obtain a map of the expression of its chemosensory receptors, verifying whether VRs play a role in the main olfactory network. b) Formyl peptide receptors The formyl peptide receptor family was first described in immune cells [93]. In mammals, FPR1 and FPR2 are believed to play an important role in innate immune response and are expressed by immune cells such as granulocytes and monocytes [94,95]. These two genes are expressed in the rabbit VNO transcriptome. Despite our study does not work at cellular resolution, considering the VNO tissue heterogeneity (blood vessels, immune cells, glands, etc.), one could speculate that FPR1 and FPR2 might be expressed in immune cells and not in vomeronasal neurons. Mammalian FPRs have a complex evolutionary history [96]; while a single copy of FPR1 has been described, FPR2 has undergone dynamic episodes of duplication events widely varying among species, which points towards the frequent neofunctionalization of this subfamily. In fact, in Rodentia a third family of vomeronasal receptors, Fprs, exclusively expressed in vomeronasal tissue extracts has been reported [67,16], representing an expansion of FPR2 that shows variation between species (for instance, Fpr-rs4, 6 and 7 paralogs are specific to Muridae species). While VNO Fprs specific function is not clear, it is widely thought that they play a role in VNO pathogen sensing [97]. Indeed, Bufe et al. [98] have recently determined that bacterial MgrB peptide activates Fpr3 in a subset of VSNs and drives avoidance behaviour, representing the only functional study of any Fpr in the VNO to date. There was a wide consensus on the rodent-specific nature of the VNO Fprs expansion, but until recently most studies had been restricted to a reduced number of rodents and primates [99,100]. Lately several new FPR2 paralogs have been descriped in Chiroptera and Perissodactyla [96], which suggests that FPR2 expansions are not restricted to Rodentia. The relation of these non-rodent FPR2 paralogs with the VNO remains unknown and further studies are needed to determine if they also have a functional role in the VNO. In any case, in rabbits FPR2 has apparently not undergone any duplication, and the only two FPR genes found in the rabbit VNO transcriptome were FPR1 and FPR2. c) Transient receptor potential channels Pheromone signals enter the vomeronasal organ and activate essential components of the signal transduction machinery. Even though our understanding of the vomeronasal signaling cascade remains incomplete, the need of a specific transient receptor potential channel (Trp) to guarantee signal transduction seems to be clear [89]. Trpc2 has been detected in VSNs, including microvilli, which represent the primary signal transduction compartment [101,102]. While we found expression of this gene only in the VNO when compared to the rabbit atlas expression, it is also expressed in testis, sperm, the dorsal root ganglion and the brain of mice [103], and in rat thyroid cells [104]. Moreover, it has also been detected in olfactory sensory neurons within the MOE of mice and lemur [91,105]. Since the MOE was not part of the rabbit expression atlas, further studies regarding Trpc2 expression would aid to elucidate whether this gene is also involved in olfactory signaling through the main olfactory pathway in rabbits. We found several Trp genes with differential expression between adults and juveniles. Trpc2 is more expressed in the VNO of adult rabbits, while three other Trp transcripts -Trpm7, Trpm8 and Trpc1- not previously reported as expressed in the mammalian VNO, are more expressed in juveniles. We speculate that, despite signal-transduction mediated by Trpc2 may be enhanced in the VNO of mature animals (6 months old), as suggests the up-regulation of this gene, the higher expression of Trpc1 in juveniles suggests that it might assume Trpc2 function in immature animals. Alternatively, a given member of the Trp family or even a combinatorial Trp expression may trigger or be involved in a specific effect depending on the perceived chemical, as supported by the increased expression of Trpm7 and Trpm8 in juveniles. In short, different members of the Trp family might be implicated in pheromone-signaling in both juveniles and adults, and their combinatorial function in the VNO could be significantly different from their function when expressed individually in other cell types. Additionally, despite not detected in the OryCun 2.0, we found Trpm5 and Trpm4 in the OryCun 3.0, and alignment against the mice transcriptome of the RNA-Seq reads that did not map to the rabbit assembly suggests that both genes might be moderately expressed in the rabbit VNO. Trpm5 was found in the mice non-sensory VNO but not in VSNs [106,107], and despite its function in the VNO is still unclear, its presence in the supporting cell layer of mice MOE has revealed its protective role for maintaining the olfactory function [108]. Trpm4 has also been described in the VNO of other mammals, and has been involved in vomeronasal circuitry in mice, showing dimorphic expression in male and female VSNs depending on female estrous cycle [106]. Nonetheless, an improved rabbit genome assembly is required to confirm whether these genes exist in rabbit, as well as to evaluate where exactly they are expressed at cellular resolution level and whether Trpm4 also shows sexually dimorphic expression in rabbits. d) Major urinary proteins. Mup4 in the rabbit nasal mucosa The lipocalin family is known to be involved in immune response and pheromone transport [109], and it contains major urinary proteins (Mups), which are directly involved in vomeronasal stimuli and communication of information in urine-derived scent marks [110,111]. Mups genomic complexity is consistent with their species-specific function in mammals [112]. To date, expression of these proteins has only been found in rodents, especially in urine, liver and exocrine glands [113,114,115,116]. However, their expression was also reported in the VNO [117]; for instance, Mup4 is known to be expressed in mice vomeronasal mucosa and it is highly specific for the male mouse pheromone 2-s-butyl-4,5-dihydrothiazole (SBT), which promotes aggression among males while inducing synchronized estrus in females [118,119]. Here we provided the first evidence that several lipocalins, including Mup4, belong to the most expressed genes in the rabbit VNO transcriptome and are VNO-specific. This is consistent with a previous study from Ibarra- Soria et al. [26] where Mup4 was found among the most expressed genes in the mice VNO transcriptome. Additionally, Broad and Keverne [22] found Mup 1–3 and 5 in the glandular region of the mice VNO, thus suggesting a role of nasal Mups in binding pheromones and likely presenting them to their receptors. We did not find such proteins in the rabbit VNO transcriptome, and despite P.R. Villamayor et al. Genomics 113 (2021) 2240–2252 2248 more in-depth studies are needed to discard its presence, it is possible that such Mups are rodent-specific. Finally, vomeronasal Mups may provide an additional means of selectively modulating the activation of vomeronasal receptors, and more species-specific studies are needed to disclose the vomeronasal Mups repertoire of each species, as well as the mechanisms underlying pheromone-receptor activity. 4.2. VNO-mediated reproductive behaviour does not depend on VNO transcriptome sex differences Males and females exhibit striking sexually dimorphic behavioural responses to pheromones [120,121,110], but how those signals are processed by the sensory systems and central circuits needs to be elucidated. It could be hypothesized that behavioural differences might have its origin on sex-specific expression of sensory receptors; however, we hardly found differences between sexes in the rabbit VNO transcriptional repertoire (only 12 DE genes, none of them VRs neither involved in essential VNO function), which supports previous findings where VR expression was found largely identical between males and females in mice [29,26]. To our knowledge, a single mice sensory receptor, Olfr692, has been reported to show a sexual dimorphic expression in a subpopulation of VSNs in response to pup odors [122]. We did not find that gene in the rabbit genome, but the mice sequence gave a significant homology against the rabbit OR52W1, which showed very low expression in our VNO transcriptome. While post-translational changes as well as variation in signal transduction mechanisms could explain sex behavioural dimorphism, it is more likely that pheromones are equally detected by both sexes but interpreted in different ways due to sexually dimorphic central circuits [123]. Indeed, previous studies have found a significant sexual dimorphism in the density of the accessory olfactory bulb (AOB) principal cells of rabbits [41] and rats [124], as well as in the medial amygdala [125, 126]. The important role of hormonal state in behavioural modulation should also be considered. Accordingly, several key hormone receptors are known to be expressed in pheromone-responsive nodes within the brain, contributing to sexually dimorphic social behaviour [123]. Internal hormones can also trigger behavioural changes through the modulation of vomeronasal sensory neurons before information is delivered to the brain. For instance, progesterone, a sex-specific hormone which is mainly produced during female diestrus, acts directly in a subset of female VSNs expressing the progesterone receptor membrane component 1 (Pgrmc1), and inhibits their ability to detect male pheromones (Mups). Therefore, female mice are “blind” to male Mups in diestrus, but they display attraction towards the same ligands in estrus, when progesterone is present at basal levels [127]. The VNO is known to play an important role in reproductive function (i.e sexual attraction or maternal behaviour) [128]. Our rabbit VNO transcriptomic analysis showed several enriched GO terms and KEGG pathways as well as hormone receptor genes closely related to the hormone/reproductive physiology of the VNO itself. Specifically, progesterone receptors membrane component 1 and 2 (Pgrmc1 and Pgrmc2) showed high expression in the rabbit VNO. Future studies on the activation/silencing of VSNs under progesterone exposure across rabbits ovulation cycle might help understand the mechanisms underlying state-specific behaviour in this species. Additionally, we found progesterone receptor (Pgr) lowly expressed in our samples, which is consistent with previous studies of the VNO in pregnant female mice [30]. This gene was also more expressed in adult rabbits, suggesting a broader function of this receptor in mature animals. Moreover, some estrogen receptors are found in male and female mice VSNs [129] and their expression depends on the animal condition. For instance, the estrogen receptor 1 (Esr1) showed a prominent expression in the VNO of pregnant female mice [30]. In rabbits, estrogen receptors showed poor expression, but further studies considering rabbit females at different stages –estrus, diestrus, and pregnancy– may help understand their importance in estrogen-detection in this species. Recently, oxytocin has also been involved in regulation of VNO activity, reducing pup-direct aggression of male mice [130]. In our case, we found functional enrichment of the oxytocin signaling pathway (pvalue <0.05), but the oxytocin receptor (OXT) was hardly expressed in the rabbit VNO. More specifically, Gli3, a gene that controls GnRH production, was moderately expressed in the rabbit VNO. This is consistent with its expression in the apical portions of the developing mice VNO in maturing rabbits and could be related to exerting control of pubertal onset and fertility [61]. Additionally, we found moderately expressed genes such as Dmxl2 and Lgr4, which might be related not only to fertility onset but also to its maintenance throughout reproductive life [62,63]. Dmxl2 has been previously involved in mice olfactory mucosa and Dmxl2-knockouts have shown deficiencies in olfactory signal transmission [63]. However, this is the first time in which its expression in the VNO is shown. Interestingly, Dmxl2 exerts its principal function in the testes at the onset of puberty and was previously found highly expressed in spermatogonia and spermatocytes, therefore playing a dual role in olfactory information and first wave spermatogenesis [131]. Moreover, Lgr4 is strongly expressed in the VNO during development in mice, and its deficiency impairs GnRH neuron development and delays puberty [62]. GnRH neurons migrate along axons of cells that reside within the VNO to the forebrain, and therefore the expression of Lgr4 in the VNO may be related to the GnRH pathway. Accordingly, we found this gene moderately expressed in the rabbit VNO. Other genes of this family such as Lgr5 and Lgr6 were found moderately expressed in the rabbit VNO, and Lgr6 was more expressed in juveniles than in adults. Although Lgr5 has been previously found in MOE and taste stem cells [132], to our knowledge this is the first time in which the expression of both genes is reported in the VNO. Altogether, we have detected genes connected to several sex hormone-related pathways in the rabbit VNO, thus confirming that this phenomenon is not restricted to mice. Its biological relevance should be further explored in a behavioural context to define the potential role of sex hormones in modulating VNO-mediated behaviours in rabbits. 4.3. Insights into vomeronasal connectivity Sensory systems are responsible of relaying information from the environment to the central nervous system. Specifically, the VNO is implicated in transforming chemical cues into electrical signals, in which the ion channels anoctamin 1 and 2 (Ano1 and Ano2) play a fundamental role [133]. These genes contributing to vomeronasal signal amplification [134] have been found expressed in the rabbit VNO. Additionally, Ano2 and Ano10 were more expressed in adults than in juveniles, an observation not previously reported in any species. The functional maturation and connectivity of basal VSNs in mice is driven by Smad4 through the canonical BMP/TGFβ signaling pathway [135]. This gene showed the highest expression (>773 normalized reads) among the 8 belonging to Smad family in the rabbit VNO transcriptome. We also found Smad5, a gene not previously described as expressed in the VNO, more expressed in adults. Additionally, 21 BMP and 22 TGFβ genes were expressed in the rabbit VNO, seven and three of them differentially expressed between juveniles and adults, respectively. The high expression of BMP4 is quite relevant since it plays an important role in defining neurogenic fate in the developing vomeronasal system [136]. Additionally, Tfap2 family, specifically Tfap2e, is expressed in mice basal VSNs and plays an important role in the maintenance of their identity [137]. We found this gene lowly expressed in the rabbit VNO. Other genes of this family such as Tfap2a and Tfap2c were moderately expressed in our samples, the latter being overexpressed in adults. Despite the plasticity of olfactory sensory neurons has been more broadly studied, the VSNs of the VNO also regenerate throughout life [138]. In rabbits, the constant regeneration and differentiation of VSNs is reflected by stem cell and Notch-related enriched terms. In fact, we found 11 Notch genes differentially expressed, with Delta-like 4 signaling gene (Dll4) more expressed in adults, and delta like non- P.R. Villamayor et al.