- Research article
- Open Access
Convergent reduction of V1R genes in subterranean rodents
BMC Evolutionary Biologyvolume 19, Article number: 176 (2019)
Vomeronasal type 1 receptor genes (V1Rs) are expected to detect intraspecific pheromones. It is believed that rodents rely heavily on pheromonal communication mediated by V1Rs, but pheromonal signals are thought to be confined in subterranean rodents that live in underground burrows. Thus, subterranean rodents may show a contrasting mode of V1R evolution compared with their superterranean relatives.
We examined the V1R evolution in subterranean rodents by analyzing currently available genomes of 24 rodents, including 19 superterranean and 5 subterranean species from three independent lineages. We identified a lower number of putatively functional V1R genes in each subterranean rodent (a range of 22–40) compared with superterranean species (a range of 63–221). After correcting phylogenetic inertia, the positive correlation remains significant between the small V1R repertoire size and the subterranean lifestyle. To test whether V1Rs have been relaxed from functional constraints in subterranean rodents, we sequenced 22 intact V1Rs in 29 individuals of one subterranean rodent (Spalax galili) from two soil populations, which have been proposed to undergo incipient speciation. We found 12 of the 22 V1Rs to show significant genetic differentiations between the two natural populations, indicative of diversifying selection.
Our study demonstrates convergent reduction of V1Rs in subterranean rodents from three independent lineages. Meanwhile, it is noteworthy that most V1Rs in the two Spalax populations are under diversifying selection rather than relaxed selection, suggesting that functional constraints on these genes may have retained in some subterranean species.
Olfaction has been widely believed to trigger essential animal behaviors such as mate choice, food location and predator avoidance . In most tetrapod vertebrates, two olfactory systems are present: the main olfactory system (MOS) and the vomeronasal system (VNS) . Through distinct signal transduction pathways, the two olfactory systems act in concert to detect odorants . The intraspecific pheromones are species-specific odor molecules that trigger sexual and social behaviors in conspecific individuals , which are more likely to evoke the VNS chemoreceptors that show an evolutionary pattern with a species-specific manner; By contrast, the MOS chemoreceptors show an evolutionary pattern with a relatively conserved manner [3, 5]. There are at least three families of receptor genes that are expressed in the VNS, including vomeronasal type 1 receptor genes (V1Rs), vomeronasal type 2 receptor genes (V2Rs), and formyl peptide receptor genes (FPRs) [6,7,8]. Of the three gene families, V1Rs are of particular interest for evolutionary analysis within a phylogenetic framework, because the V1R gene family in mammals shows dramatic among-species variation in size, from 0 in two bat species, dolphin and rhesus macaque to 283 in platypus [9,10,11,12]. Earlier studies have reported that five surface-dwelling rodents (superterranean rodents) have an expanded V1R repertoire and suggested that V1Rs play an important role in pheromonal communication in these animals [12, 13]. Indeed, mutant mice lacking 16 intact V1R genes showed a reduced level of maternal aggression and male sexual drive . To our knowledge, however, all studies published so far were focused on V1R gene repertoires in superterranean rodents, little is known about evolution of V1R repertoires in subterranean rodents.
Subterranean rodents live in underground burrows and spend most of their life below the soil surface . Due to the limited space in burrows, chemical communication was thought to be confined in subterranean rodents . However, there was evidence supporting that subterranean rodents may use pheromonal communication . To study the evolution of V1Rs in subterranean rodents, we examined V1R evolution by analyzing currently available genomes of 24 rodents. We next sequenced 22 intact V1R genes in 29 individuals of a subterranean rodent species (Spalax galili) from two soil populations, with the aim of testing functional importance of V1Rs in natural populations of subterranean rodents.
We identified V1R genes in 24 currently available rodent genome sequences, including five subterranean and 19 superterranean species (Fig. 1). A total of 4763 V1R genes were identified in the present study, including 2244 putatively functional V1R genes and 2519 pseudogenes (Additional file 1: Table S1). We divided putatively functional genes into two categories: intact genes and partial genes. The former contain an intact open reading frame (ORF) and a complete coding region, whereas the latter retain an intact ORF and a partial coding region due to incomplete sequencing [17,18,19,20,21]. These 2244 putatively functional V1R genes include 1983 intact genes (mean 83, median 88) and 261 partial genes (mean 11, median 9) (Additional file 1: Table S1). The total number of putatively functional V1Rs in each species varied from 22 in the naked mole rat (Heterocephalus glaber) to 221 in the mouse (Mus musculus) (Additional file 1: Table S1). Notably, subterranean rodents (mean 24, median 22) tend to have a lower number of intact V1R genes than their superterranean relatives (mean 98, median 96) (Fig. 1). The same tendency was also observed when both intact and partial V1Rs were considered (Additional file 1: Table S1). The reduction of V1R genes has occurred repeatedly in five subterranean rodents from three different lineages (Fig. 1), suggestive of convergent reduction. In addition, we constructed a gene tree (Additional file 2: Figure S1) using all intact V1Rs from the 24 species of rodents. Our phylogenetic analyses revealed that different subterranean lineages with similar sizes of small gene repertoires tend to retain different gene groups (Additional file 2: Figure S1).
Phylogenetically independent contrast (PIC) analysis
To determine whether the reduction of V1R genes has resulted from the subterranean lifestyle shared by the five species, we performed phylogenetically independent contrast (PIC) analysis based on the V1R genes identified from 24 rodent species (Additional file 1: Figure S2, Table S2). Due to the tendency of fewer V1R genes in subterranean rodents, we coded the lifestyle of a rodent as 0 (subterranean species) or 1 (superterranean species). We converted the lifestyle codes and the V1R gene numbers into phylogenetically independent contrasts (PICs)  and subsequently conducted a regression analysis between lifestyle code contrasts and gene number contrasts. This analysis revealed a significant positive correlation between the PICs of the intact V1R gene numbers and those of lifestyle codes (Spearman’s ρ = 0.424, P = 0.044; Fig. 2a). The same trend was identified when the PICs of the total numbers of putatively functional V1Rs (intact and partial genes) were correlated with the PICs of lifestyle codes (Spearman’s ρ = 0.432, P = 0.040; Fig. 2b). We also performed PIC analysis based on the V1R genes identified from 18 species that were sequenced on the Illumina platform and obtained similar results (Additional file 1: Figure S3). These findings strongly suggest that the lifestyle impacts the repertoire size of functional V1Rs, with a smaller size in subterranean rodents relative to their superterranean relatives (Fig. 1, Additional file 1: Table S1).
Population genetic differentiation of V1Rs in Spalax galili
To characterize the evolution of those survival V1Rs in subterranean rodents, we performed population genetic analyses of V1Rs between two natural populations of one subterranean rodent species (the blind mole rat, S. galili) (Fig. 3a). Our samples were collected in a recent study , including 29 individuals of S. galili from two soil populations (16 from the basalt and 13 from the chalk) (Fig. 3a). From the publicly available genome of S. galili (GenBank assembly GCF_000622305.1) , we identified 23 putatively functional V1R genes with complete coding regions and intact ORFs. Twenty-two of the 23 genes were successfully amplified and sequenced in the 29 individuals of S. galili (Fig. 3b, Additional file 1: Table S3). To test whether these V1R genes have undergone relaxed selection, we compared V1R genes with neutrally evolved sequences, which were taken from a recent study . These sequences that were assumed to be under neutral evolution were 18 randomly selected noncoding regions, which were sequenced in the same 29 mol rats . We identified single nucleotide polymorphisms (SNPs) in each locus from each animal but did not observe any SNPs that are fixed within either population, suggesting that the two populations diverged very recently. In line with this finding, genomic evidence has estimated a divergence time of 0.2–0.4 million years ago between the two populations .
The intra-population genetic variations of 22 V1R genes and 18 noncoding regions are presented in Additional file 1 Table S4 and Table S5, respectively. We found the average number of nucleotide differences per site (π) at the 18 noncoding regions to be same in both populations (0.0014) (Additional file 1: Table S5). By contrast, V1Rs showed an overall higher π in the chalk population (0.0018) than in the basalt population (0.0014) (Additional file 1: Table S4), although the difference is not statistically significant (P = 0.24, two-tailed paired t-test). After comparing with the noncoding regions, we found V1Rs to have a same π (0.0014) in the basalt population but an overall higher π (0.0018) in the chalk population (P = 0.45, two tailed Mann-Whitney U test) (Additional file 1: Tables S4, S5). We also calculated the Watterson’s polymorphism per site (θ) in V1Rs and noncoding regions, and obtained similar results to the estimates of π (Additional file 1: Tables S4, S5). Tajima’s neutrality test  was used to assess whether a V1R gene or a noncoding region fit the neutral theory model. The rejection of the null hypothesis of neutrality suggests a non-neutral process of sequence evolution, including purifying selection, positive selection, and demographic changes . We found zero V1R and four noncoding regions that showed significantly positive or negative Tajima’s D (Additional file 1: Table S4, S5), suggesting all V1Rs fit the neutral theory model and the majority of noncoding regions also fit. In the basalt population, one noncoding region showed significantly positive D, whereas two noncoding regions showed significantly negative D (Additional file 1: Table S5). In the chalk population, only one noncoding region showed a significantly negative D (Additional file 1: Table S5). However, we did not identify a significant difference between the basalt (3/18 = 16.7%) and chalk (1/18 = 5.6%) populations in the fraction of noncoding regions with significantly positive or negative D (P > 0.6, Fisher’s exact test). In either population, the fraction of loci with significantly positive or negative D in V1Rs (Additional file 1: Table S4) is not significantly different from that in noncoding regions (Additional file 1: Table S5) (P = 0.23 for the basalt and P = 1 for the chalk population, Fisher’s exact test). We also performed Fu and Li’s neutrality test  and found similar results to Tajima’s test. Neither population shows a significant difference between noncoding regions and V1Rs in the fraction of loci with a significantly different Fu and Li’s D (Additional file 1: Tables S4, S5) (P = 0.16 for the basalt and P = 0.65 for the chalk population, Fisher’s exact test).
The genetic differentiations between the chalk and basalt populations were measured as the fixation index (termed FST). Pairwise FST statistics between the two populations are provided in Fig. 3b. We found that 14 V1Rs are significantly differentiated between populations (P < 0.05 after false discovery rate (FDR) adjustment) (Fig. 3b). Only one noncoding region was found to be significantly differentiated between populations (P < 0.05 after FDR adjustment) . The percentage of significantly differentiated loci is statistically greater for V1Rs (14/22 = 63.6%) than for noncoding regions (1/18 = 5.6%) (P < 0.001, Fisher’s exact test), suggesting that positive selection may have shaped the V1R evolution. We also compared the frequency of each nonsynonymous SNP in the chalk population with that in the basalt population, with the aim to detect significantly differentiated nonsynonymous SNPs that may be of functional importance. Briefly, we used Fisher’s exact test to examine the difference of the frequency of each nonsynonymous SNP between the two populations. Nonsynonymous SNPs with adjusted P-values less than 0.05 were considered significantly differentiated nonsynonymous SNPs (Additional file 1: Table S6). We found 12 of 22 V1Rs to have at least 1 significantly differentiated nonsynonymous SNP, with the number ranging from 1 to 6 (Additional file 1: Figure S4). With one exception, the 12 V1Rs in this analysis could also be identified by the FST analysis (Fig. 3b). The exception is V1R6, which showed overall genetic differentiation between populations (Fig. 3b) but did not have nonsynonymous SNPs (Additional file 1: Figure S4). We further examined the distributions of these significantly differentiated nonsynonymous SNPs on each V1R gene by predicting V1R transmembrane protein topologies (Additional file 1: Table S6). Our prediction revealed that most of these SNPs occurred in putative ligand binding sites [28,29,30]. These findings unambiguously demonstrate nonsynonymous substitutions contribute to the genetic and functional divergence of V1Rs in the two natural populations of the subterranean rodent species.
Our data show that the number of functional V1Rs was markedly reduced in 5 subterranean rodents from 3 different lineages. A previous study has identified a positive correlation between the number of intact V1Rs and the complexity of vomeronasal organ (VNO) morphology across mammals . Our results are consistent with the known information about VNO morphology of studied rodents. For example, Heterocephalus glaber, one of the five subterranean rodents in this study, has a small and growth-deficient VNO , while other studied superterranean rodents, like Mus musculus, Microtus ochrogaster, Rattus norvegicus and Chinchilla laniger, have a well-developed VNO [32,33,34,35]. Combined our genetic data with previous morphological evidence, we concluded that pheromonal olfaction mediated by V1Rs is commonly reduced in subterranean rodents. In addition, the small V1R repertoire in subterranean rodents can not be compensated by a large V2R repertoire, because the majority of the V2R repertoires in Spalax galili and H. glaber are characterized by pseudogenes . Thus, we speculate that the function of VNS in subterranean rodents is less important than their superterranean relatives. By contrast, the MOS, which mainly expresses olfactory receptors (ORs), was expected to play a more important role in African mole-rats. Specifically, the naked mole rat possesses 1401 functional ORs, which are around three times more than its superterranean relative, the guinea pig . Thus, ORs seem to have shown an evolutionary pattern distinct from V1Rs in rodents.
The observation of reduction or complete loss of V1R repertoire in mammals may have resulted from a number of ecological factors, such as a sensory tradeoff between vision and pheromonal olfaction in hominids and Old World monkeys , and an aquatic lifestyle in dolphins . Why did the reduction of V1R repertoire occur in subterranean rodents? Previous studies suggested that chemical signals are relatively less important than vibrational signals in these animals . In other words, living in the underground burrow confined the dispersal of chemical signals, whereas vibrational signals can spread beyond the confines of underground burrow . We thus infer that selective pressures acting on V1R repertoire may have been relaxed because of ineffective chemical communication in subterranean rodents. In addition, we also identified the V1R gene repertoires in four species of soricomorphans, including one subterranean species (Condylura cristata, star-nosed mole) and three superterranean relatives (Erinaceus europaeus, Solenodon paradoxus, Sorex araneus). We found that the star-nosed mole possesses the fewer intact V1R genes (n = 15) compared with its superterranean relatives (n = 56 in E. europaeus and n = 54 in S.araneus), which is consistent with the overall pattern in rodents (Fig. 1), although the Hispaniolan solenodon (S. paradoxus) carries the fewest number of intact V1Rs (n = 6) (Additional file 1: Table S1).
The substantial reduction of functional V1Rs does not represent nonfunctional pheromonal olfaction in subterranean rodents, because considerable evidence has suggested that pheromone cues are used in subterranean rodents such as Spalax, Heterocephalus, ctenomyids, and bathyergids . Indeed, subterranean rodents have been assumed to rely on chemical cues, because life underground deprives them of visual cues, which otherwise represent an important type of sensory stimuli in mammals . In addition, pheromones are relatively easy to synthesize and are also long-lasting cues , which may provide an inexpensive mode of communication for subterranean rodents coping with multiple stressors such as darkness, energetics, hypoxia, and hypercapnia . Our population genetic data include intra-population polymorphisms and inter-population variations. The intra-population polymorphisms in V1Rs are similar to those in noncoding regions, suggesting that V1R evolution is largely neutral within either population of the blind mole rats. In accordance with our study, earlier studies in humans , mice [26, 42], and other mammals  also revealed that genetic drift plays an important role in V1R evolution. By contrast, our inter-population analysis (FST) showed significantly greater genetic differentiation of V1Rs than neutrally evolving noncoding regions, with the fraction of significantly differentiated loci much higher for V1Rs (14/22) than for noncoding regions (1/18). Furthermore, a number of significantly differentiated nonsynonymous SNPs were observed in the significantly differentiated V1Rs between populations, implying that diversifying selection may also play an important role in V1R functional changes. Thus, the microevolution of V1Rs in S. galili is at least in part shaped by natural selection, supporting the hypothesis that V1Rs may be involved in reproductive isolation of S. galili. In addition, several studies have suggested that pheromonal olfaction mediated by the VNS is involved in Spalax reproductive isolation: S. galili exhibits mate choice by choosing mates with similar genetically determined odors ; A behavioral test demonstrated that reproductive isolation associated with olfaction plays a major role in Spalax speciation across Israel ; Olfactory discrimination was observed to serve as a reproductive isolating mechanism in Spalax speciation ; Pheromones identified from the urine of male mole rats were found to be involved in sexual attractance . Together, our study shows that functional V1Rs in subterranean rodents are significantly reduced in number and may remain important in some species of subterranean rodents.
Our study suggests that the subterranean or superterranean lifestyle shaped the evolution of V1Rs in rodents. Subterranean rodents lived in underground burrows have a smaller V1R repertoire than their superterranean relatives. However, we noticed that the small V1R repertoire of Spalax galili are under diversifying selection between two populations, indicative of functional significance of these genes. It would be interesting to test whether V1Rs remain significant genetic differentiation in other populations or species of subterranean rodents.
Identification of V1R repertoires
A total of 24 rodent genome sequences were retrieved from the National Center for Biotechnology Information database (http://www.ncbi.nlm.nih.gov/). Of them, five species are subterranean, while the remaining 19 species are superterranean (Additional file 1: Table S1). The species names and genome accession numbers are given in Additional file 1: Table S1. Because vertebrate V1Rs are single-exon genes encoding vomeronasal type 1 receptors with seven transmembrane domains , we generally followed two previous studies [11, 20] to identify the V1R gene repertoires in the 24 rodents. Briefly, TblastN searches  were conducted on the genome sequences using phylogenetically diverse V1Rs that were obtained from previous studies, including full-length V1R genes from the mouse, rat, cow, dog, two fishes, and one ancV1R of Hippopotamus amphibius [11, 12, 48,49,50] as queries, with a cutoff e-value of 1e-5. The blast hits shorter than 300 bp were discarded, and the remaining blast hits were extended in both directions to acquire the start or stop codons. All full-length candidate genes were confirmed for the presence of seven-transmembrane domains using the TMHMM  or HMMTOP method . A putative gene was considered to be real when the BLASTN searches against the non-redundant database of GenBank yielded the best hit being the known V1R gene. All candidate genes were classified into three categories: intact genes, partial genes and pseudogenes. Intact genes contain at least 270 codons, a putative start and stop codon, and an intact ORF; Partial genes contain more than 300 nucleotides, a putative start or a stop codon, and a truncated ORF due to incomplete sequencing or poor genomic assembly. Pseudogenes contain more than 300 nucleotides and an interrupted ORF because of nonsense or frameshift mutations. We conducted BLASTP searches  against non-redundant database of GenBank, and chose the best hit as the basis of nomenclature of rodent V1Rs. By contrast, we named the 22 V1R genes from S. galili numerically with the order in which they were identified from the genome, because 16 of them were annotated as V1R4 in GenBank and 4 of them as V1R1 (Additional file 1: Table S3). For simplicity and convenience, we used V1R1-V1R22 to name the 22 V1R genes in Spalax galili. All intact V1R genes newly identified in this study were provided in Additional file 3. The information about gene location and orientation on the scaffolds for each V1R was provided in the Additional file 4. The phylogenetic tree (Additional file 2: Figure S1) of all intact V1R sequences was reconstructed using the Bayesian method  with 6 million generations, and hTAS2R1 (GenBank accession number: NM_019599.2) was used as an outgroup.
Phylogenetically independent contrast analysis
To assess the potential impact of the lifestyle categories (subterranean or superterranean lifestyle) on the evolution of V1Rs in rodents, we conducted a regression analysis of V1R gene number against the lifestyle category following a recent study . Briefly, we coded each rodent as 0 (subterranean species) and 1 (superterranean species) by assuming that subterranean rodents may contain fewer V1Rs than their superterranean relatives. We conducted a phylogenetically independent contrast (PIC) analysis using the package APE (Analyses of Phylogenetics and Evolution) . The PIC method was used to quantify the correlation between V1R gene number and the lifestyle after eliminating the confounding effects of phylogenetic inertia. The input tree (Additional file 1: Figure S2) is the established species tree , and the branch lengths of the input tree were estimated by divergence times (Additional file 1: Table S2). The regression analysis was performed after converting the lifestyle codes and V1R gene numbers into PICs. The correlation between the lifestyle and V1R gene number was assessed by the nonparametric Spearman’s rank correlation coefficient (ρ).
DNA sequencing and population genetic analysis
By searching the published genome sequence of the blind mole rat Spalax galili (GenBank assembly: GCF_000622305.1) , we identified 23 full-length V1Rs with an intact open reading frame, suggesting that they are putatively functional. A suite of primers (Additional file 1: Table S3) were designed to amplify all V1R genes except for the ancV1R based on the corresponding flanking sequences. A total of 29 individuals of S. galili, 16 from the basalt soil and 13 from the chalk soil, were captured alive near Rehaniya, Upper Galilee Mountains in northern Israel (33°02.5′N, 35°29.2′E) in a previous study . Animals were euthanized by injecting Ketaset CIII at 5 mg/kg of body weight. Muscle tissues were sampled and then stored in 95% (vol/vol) ethanol to be kept at − 80 °C until DNA extraction. Genomic DNAs were isolated from the muscle tissue using the DNeasy Blood and Tissue Kit (Qiagen). PCRs (polymerase chain reactions) were processed with high-fidelity KOD-Plus-Neo DNA polymerase (Toyobo), and details of PCR amplication were described previously [17,18,19, 21, 56]. PCR products were purified by the QIAquick PCR Purification Kit (Qiagen) and sequenced directly from both strands. A total of 1276 sequences that were newly obtained were deposited in GenBank with accession numbers MF118921-MF120196, and we included both alleles for each gene in the data deposition. All sequences were edited and aligned by MEGA6 . Nucleotide diversity π, a measure of genetic variation, is defined as the average number of nucleotide differences per site . Another important parameter of genetic diversity is Watterson’s θ, which was estimated by the number of segregating sites . Tajima’s D and Fu and Li’s D* [25, 27] are classical statistical tests to assess whether DNA sequences evolve neutrally. The significance of Tajima’s D test or Fu and Li’s D* test was conducted by 10,000 replicates of coalescent simulation. The population genetic differentiation of each gene between the basalt and chalk populations was measured by the fixation index (termed FST) . The nearest-neighbor statistic(Snn) test  was conducted to evaluate the significance level of the genetic differentiation between the two populations. The estimates of significance level were corrected for multiple comparisons using the FDR adjustment. All these parameters and statistical tests of population genetic differentiation were computed and performed using DnaSP v5.10 .
Availability of data and materials
All DNA sequences of Spalax V1Rs generated in this study can be found under the GenBank accession numbers: MF118921-MF120196. The accession numbers for all other genetic data used in this study has been included in the manuscript and its additional files.
False discovery rate
Formyl peptide receptor gene
Main olfactory system
Open reading frame
Polymerase chain reaction
Phylogenetically independent contrast
- Snn test:
Nearest-neighbor statistic test
Single nucleotide polymorphism
Vomeronasal type 1 receptor
Vomeronasal type 2 receptor
Doty RL. Mammalian olfaction, reproductive processes, and behavior. New York: Academic Press; 1976.
Dulac C, Torello AT. Molecular detection of pheromone signals in mammals: from genes to behaviour. Nat Rev Neurosci. 2003;4:551–62.
Grus WE, Zhang J. Distinct evolutionary patterns between chemoreceptors of 2 vertebrate olfactory systems and the differential tuning hypothesis. Mol Biol Evol. 2008;25:1593–601.
Wyatt TD. Pheromones and animal behaviour: communication by smell and taste. New York: Cambridge University Press; 2003.
Zhang XH, Zhang XM, Firestein S. Comparative genomics of odorant and pheromone receptor genes in rodents. Genomics. 2007;89:441–50.
Matsunami H, Buck LB. A multigene family encoding a diverse array of putative pheromone receptors in mammals. Cell. 1997;90:775–84.
Riviere S, Challet L, Fluegge D, Spehr M, Rodriguez I. Formyl peptide receptor-like proteins are a novel family of vomeronasal chemosensors. Nature. 2009;459:574–7.
Mombaerts P. Genes and ligands for odorant, vomeronasal and taste receptors. Nat Rev Neurosci. 2004;5:263–78.
Young JM, Kambere M, Trask BJ, Lane RP. Divergent V1R repertoires in five species: amplification in rodents, decimation in primates, and a surprisingly small repertoire in dogs. Genome Res. 2005;15:231–40.
Grus WE, Shi P, Zhang J. Largest vertebrate vomeronasal type 1 receptor gene repertoire in the semiaquatic platypus. Mol Biol Evol. 2007;24:2153–7.
Grus WE, Shi P, Zhang YP, Zhang J. Dramatic variation of the vomeronasal pheromone receptor gene repertoire among five orders of placental and marsupial mammals. P Natl Acad Sci U S A. 2005;102:5767–72.
Young JM, Massa HF, Hsu L, Trask BJ. Extreme variability among mammalian V1R gene families. Genome Res. 2010;20:10–8.
Wang GD, Shi P, Zhu ZH, Zhang YP. More functional V1R genes occur in nest-living and nocturnal terricolous mammals. Genome Biol Evol. 2010;2:277–83.
Del Punta K, Leinders-Zufall T, Rodriguez I, Jukam D, Wysocki CJ, Ogawa S, Zufall F, Mombaerts P. Deficient pheromone responses in mice lacking a cluster of vomeronasal receptor genes. Nature. 2002;419:70–4.
Lacey EA, Patton JL, Cameron GN, editors. Life underground: the biology of subterranean rodents. Chicago: The University of Chicago Press; 2000.
Heth G, Todrank J. Assessing chemosensory perception in subterranean mole rats: different responses to smelling versus touching odorous stimuli. Anim Behav. 1995;49:1009–15.
Feng P, Zheng JS, Rossiter SJ, Wang D, Zhao H. Massive losses of taste receptor genes in toothed and baleen whales. Genome Biol Evol. 2014;6:1254–65.
Hong W, Zhao H. Vampire bats exhibit evolutionary reduction of bitter taste receptor genes common to other bats. Proc Roy Soc B. 2014;281:20141079.
Lu Q, Wang K, Lei F, Yu D, Zhao H. Penguins reduced olfactory receptor genes common to other waterbirds. Sci Rep. 2016;6:31671.
Wang K, Zhao H. Birds generally carry a small repertoire of bitter taste receptor genes. Genome Biol Evol. 2015;7:2705–15.
Zhao H, Li J, Zhang J. Molecular evidence for the loss of three basic tastes in penguins. Curr Biol. 2015;25:R141–2.
Felsenstein J. Phylogenies and the comparative method. Am Nat. 1985;125:1–15.
Li K, Hong W, Jiao H, Wang GD, Rodriguez KA, Buffenstein R, Zhao Y, Nevo E, Zhao H. Sympatric speciation revealed by genome-wide divergence in the blind mole rat Spalax. P Natl Acad Sci U S A. 2015;112:11905–10.
Fang X, Nevo E, Han L, Levanon EY, Zhao J, Avivi A, Larkin D, Jiang X, Feranchuk S, Zhu Y. Genome-wide adaptive complexes to underground stresses in blind mole rats Spalax. Nat Commun. 2014;5:3966.
Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.
Park SH, Podlaha O, Grus WE, Zhang J. The microevolution of V1r vomeronasal receptor genes in mice. Genome Biol Evol. 2011;3:401–12.
Fu YX, Li WH. Statistical tests of neutrality of mutations. Genetics. 1993;133:693–709.
Yoder AD, Chan LM, dos Reis M, Larsen PA, Campbell CR, Rasoloarison R, Barrett M, Roos C, Kappeler P, Bielawski J, et al. Molecular evolutionary characterization of a V1R subfamily unique to Strepsirrhine primates. Genome Biol Evol. 2014;6:213–27.
Schwartz TW. Locating ligand-binding sites in 7TM receptors by protein engineering. Curr Opin Biotechnol. 1994;5:434–44.
Baldwin JM. Structure and function of receptors coupled to G proteins. Curr Opin Cell Biol. 1994;6:180–90.
Smith TD, Bhatnagar KP, Dennis JC, Morrison EE, Park TJ. Growth-deficient vomeronasal organs in the naked mole-rat (Heterocephalus glaber). Brain Res. 2007;1132:78–83.
Wilson KCP, Raisman G. Age-related changes in the neurosensory epithelium of the mouse vomeronasal organ: extended period of postnatal-growth in size and evidence for rapid cell turnover in the adult. Brain Res. 1980;185:103–13.
Weiler E, McCulloch MA, Farbman AI. Proliferation in the vomeronasal organ of the rat during postnatal development. Eur J Neurosci. 1999;11:700–11.
Maico LM, Burrows AM, Mooney MP, Siegel MI, Bhatnagar KP, Smith TD. Size of the vomeronasal organ in wild Microtus with different mating strategies. Acta Biol Hung. 2003;54:263–73.
Oikawa T, Shimamura K, Saito TR, Taniguchi K. Fine structure of the vomeronasal organ in the Chinchilla (Chinchilla Laniger). Exp Anim. 1994;43:487–97.
Francia S, Silvotti L, Ghirardi F, Catzeflis F, Percudani R, Tirindelli R. Evolution of spatially coexpressed families of type-2 vomeronasal receptors in rodents. Genome Biol Evol. 2014;7:272–85.
Hughes GM, Boston ESM, Finarelli JA, Murphy WJ, Higgins DG, Teeling EC. The birth and death of olfactory receptor gene families in mammalian niche adaptation. Mol Biol Evol. 2018;35:1390–406.
Zhang J, Webb DM. Evolutionary deterioration of the vomeronasal pheromone transduction pathway in catarrhine primates. P Natl Acad Sci U S A. 2003;100:8337–41.
Nevo E. Mosaic evolution of subterranean mammals: regression, progression, and global convergence. New York: Oxford University Press; 1999.
Shorey HH. Pheromones. In: Sebeok TA, editor. How animals communicate. Bloomington: Indiana University Press; 1977. p. 137–63.
Nozawa M, Kawahara Y, Nei M. Genomic drift and copy number variation of sensory receptor genes in humans. P Natl Acad Sci U S A. 2007;104:20421–6.
Kurzweil VC, Getman M, Green ED, Lane RP, Sequencing NC. Dynamic evolution of V1R putative pheromone receptors between Mus musculus and Mus spretus. BMC Genomics. 2009;10:74.
Tzur S, Todrank J, Juergens A, Nevo E, Heth G. Odour-genes covariance within a natural population of subterranean Spalax galili blind mole rats. Biol J Linn Soc. 2009;96:483–90.
Heth G, Nevo E. Origin and evolution of ethological isolation in subterranean mole rats. Evolution. 1981;35:259–74.
Nevo E, Bodmer M, Heth G. Olfactory discrimination as an isolating mechanism in speciating mole rats. Experientia. 1976;32:1511–2.
Menzies RA, Heth G, Ikan R, Weinstein V, Nevo E. Sexual pheromones in lipids and other fractions from urine of the male mole rat, Spalax ehrenbergi. Physiol Behav. 1992;52:741–7.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.
Suzuki H, Nishida H, Kondo H, Yoda R, Iwata T, Nakayama K, Enomoto T, Wu J, Moriya-Ito K, Miyazaki M, et al. A single pheromone receptor gene conserved across 400 My of vertebrate evolution. Mol Biol Evol. 2018;35:2928–39.
Shi P, Zhang J. Comparative genomic analysis identifies an evolutionary shift of vomeronasal receptor gene repertoires in the vertebrate transition from water to land. Genome Res. 2007;17:166–74.
Saraiva LR, Korsching SI. A novel olfactory receptor gene family in teleost fish. Genome Res. 2007;17:1448–57.
Krogh A, Larsson B, von Heijne G, Sonnhammer ELL. Predicting transmembrane protein topology with a hidden markov model: application to complete genomes. J Mol Biol. 2001;305:567–80.
Tusnady GE, Simon I. The HMMTOP transmembrane topology prediction server. Bioinformatics. 2001;17:849–50.
Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19:1572–4.
Paradis E, Claude J, Strimmer K. APE: analyses of phylogenetics and evolution in R language. Bioinformatics. 2004;20:289–90.
Blanga-Kanfi S, Miranda H, Penn O, Pupko T, DeBry RW, Huchon D. Rodent phylogeny revised: analysis of six nuclear genes from all major rodent clades. BMC Evol Biol. 2009;9:71.
Wang K, Hong W, Jiao H, Zhao H. Transcriptome sequencing and phylogenetic analysis of four species of luminescent beetles. Sci Rep. 2017;7:1814.
Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.
Nei M, Li WH. Mathematical model for studying genetic variation in terms of restriction endonucleases. P Natl Acad Sci U S A. 1979;76:5269–73.
Watterson GA. On the number of segregating sites in genetical models without recombination. Theor Popul Biol. 1975;7:256–76.
Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984;38:1358–70.
Hudson RR. A new statistic for detecting genetic differentiation. Genetics. 2000;155:2011–4.
Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.
We thank Yi Wang and Feng-Juan Mu for technical assistance in the laboratory.
This work was supported by the National Natural Science Foundation of China (91331115, 31722051). The funders were not involved in the design of the study; the collection, analysis and interpretation of data; the writing of the manuscript and any decision concerning the publication of the paper.
Ethics approval and consent to participate
The field study and other animal ethics had been explicitly approved by the Ethics Committees of University of Haifa and Wuhan University. All procedures involving the mole rats were included in the main protocol approval, and no additional approval is required. All samples were lawfully acquired and their use was conformed to national and local laws and regulations, including the 1994 Animal Protection Law (Israel’s major animal welfare law; http://www.chai.org.il/en/compassion/legislation_protection.htm), Animal Experimentation Law (1994; https://www.aaalac.org/intlRefs/Israel/Israel-Law.pdf), and Wildlife Protection Law (1955; http://www.sviva.gov.il/English/Legislation/Documents/Wildlife%20Protection%20Laws%20and%20Regulations/WildlifeProtectionLaw1955.pdf).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Numbers of V1R genes identified from the genome assemblies of 28 mammals. Twenty-four species of rodents and four species of soricomorphans (shaded) were included. The subterranean species were shown in bold. Table S2. Node ages in the species tree depicted in the Additional file 1: Figure S2. Table S3. Primers used for amplification. Table S4. Intra-population genetic variations of 22 intact V1R genes in the two soil populations of Spalax galili. Table S5. Intra-population genetic variations of the 18 noncoding regions in the two soil populations of Spalax galili. Table S6. All significantly differentiated nonsynonymous SNPs between populations in each of the 22 V1R genes. Ten V1Rs were left blank and shown in gray due to the lack of significantly differentiated nonsynonymous SNPs. Figure S2. The species tree of 24 rodents in this study. The letters (A-W) indicate the nodes of the tree. Branch lengths were estimated from the divergence times among species (Table S2). Figure S3. PIC analysis based on the V1Rs identified from 18 species that were sequenced on the Illumina platform, with the aim to avoid biases resulting from different sequencing platforms. Figure S4. Number of significantly differentiated nonsynonymous SNPs between populations in each of the 12 V1Rs. (DOCX 545 kb)
Figure S1. The Bayesian tree based on the alignment of 1983 intact V1Rs. (PDF 3369 kb)
Nucleotide sequences of intact V1R genes identified from 28 mammalian genomes in the present study. These genes were named following the best hits of BLASTP searches. (DOCX 683 kb)
Genomic locations for all annotated V1R genes studied. (TXT 579 kb)