Table Of ContentGBE
Contrasted Evolution of the Vomeronasal Receptor
Repertoires in Mammals and Squamate Reptiles
Urszula Brykczynska1,y, Athanasia C. Tzika1,y, Ivan Rodriguez2, and Michel C. Milinkovitch1,*
1LaboratoryofArtiEcial&NaturalEvolution(LANE),DepartmentofGenetics&Evolution,UniversityofGeneva,SciencesIII,Geneva,
Switzerland
2LaboratoryofNeurogenetics&NationalResearchCenterFrontiersinGenetics,DepartmentofGenetics&Evolution,UniversityofGeneva,
SciencesIII,Geneva,Switzerland
yTheseauthorscontributedequallytothiswork.
*Correspondingauthor:E-mail:[email protected].
Accepted:January 18, 2013
Abstract
Thevomeronasalorgan(VNO)isanolfactorystructurethatdetectspheromonesandenvironmentalcues.Itconsistsofsensory
neuronsthatexpressevolutionaryunrelatedgroupsoftransmembranechemoreceptors.ThepredominantV1RandV2Rreceptor
repertoiresarebelievedtodetectairborneandwater-solublemolecules,respectively.Ithasbeensuggestedthattheshiftinhabitatof
earlytetrapodsfromwatertolandisreflectedbyanincreaseintheratioofV1R/V2Rgenes.Snakes,whichhaveaverylargeVNO
associated with a sophisticated tongue delivery system, are missing from this analysis. Here, we use RNA-seq and RNA in situ
hybridization to study the diversity, evolution, and expression pattern of the corn snake vomeronasal receptor repertoires. Our
analyses indicate that snakes and lizards retain an extremely limited number of V1R genes but exhibit a large number of V2R
genes,includingmultiplelineagesofreptile-specificandsnake-specificexpansions.Wefinallyshowthatthepeculiarbigenicpattern
ofV2Rvomeronasalreceptorgenetranscriptionobservedinmammalsisconservedinsquamatereptiles,hintingatanimportantbut
unknownfunctionalroleplayedbythisexpressionstrategy.Ourresultsdonotsupportthehypothesisthattheshifttoavomeronasal
receptorrepertoiredominatedbyV1Rsinmammalsreflectstheevolutionarytransitionofearlytetrapodsfromwatertoland.This
studyshedslightontheevolutionarydynamicsofthevomeronasalreceptorfamiliesinvertebratesandrevealshowmammalsand
squamatesdifferentiallyadaptedthesameancestralvomeronasalrepertoiretosucceedinaterrestrialenvironment.
Key words: vomeronasal organ (VNO), monogenic expression, evolution of sensorial abilities, squamates, snakes, phylogeny.
Herrada and Dulac 1997; Matsunami and Buck 1997; Ryba
Introduction
andTirindelli1997).Previousanalysesoffullysequencedver-
Thevomeronasalorgan(VNO),orJacobson’sorgan,contains tebrate genomes (Grus and Zhang 2007) indicated de novo
an olfactory sensory neuroepithelium enclosed in a cartilagi- appearance and expansion of V1R and V2R subfamilies in
nous or bony capsule in contact with the base of the nasal some lineages, as well as the presence of a significant
cavity. It plays a major role in interindividual interactions number of pseudogenes in others (Grus et al. 2005; Young
(throughthedetectionofpheromonesandkairomones)and etal.2005;ShiandZhang2007;YoungandTrask2007).The
environmental recognition (Houck 2009; Su et al. 2009). remarkable divergence between V1R and V2R gene reper-
Chemosensoryreceptorproteinsarepresentonthemicrovilli toires is thought to reflect the ecological niche and social
of vomeronasal sensory neuron dendrites and interact with habits of a particular species. Indeed, a large expansion of
moleculesintheVNOcavity.Thepredominantchemosensory the V1R family accompanied by a sharp decline in the
receptorsinthevertebrateVNOcomeinthreeclasses:formyl numberofV2Rgeneswasobservedinmammalsincompar-
peptidereceptors(FPRs),V1Rs,andV2Rs.Theseareevolutio- ison to Xenopus and fish species. Together with the notion
narily unrelated families of seven-transmembrane (7TM) that V1Rs and V2Rs might detect volatile and water-soluble
G-protein-coupled receptors (GPCRs) (Dulac and Axel 1995; molecules,respectively,theseresultsledtothehypothesisthat
(cid:2)TheAuthor(s)2013.PublishedbyOxfordUniversityPressonbehalfoftheSocietyforMolecularBiologyandEvolution.
ThisisanOpenAccessarticledistributedunderthetermsoftheCreativeCommonsAttributionNon-CommercialLicense(http://creativecommons.org/licenses/by-nc/3.0/),which
permitsunrestrictednon-commercialuse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited.
GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013 389
GBE
Brykczynskaetal.
this shift, reflected by the ratio of V1R and V2R functional Allanimalsweresexuallymature,unlessindicatedotherwise,
genes,correspondstotheevolutionarytransitionofearlytet- and2–5yearsold.
rapodsfromwatertoland(ShiandZhang2007).
Surprisingly,thevomeronasalchemosensoryreceptorgene
Histology
repertoireofsnakes,avertebratecladethatpossessesavery
largeVNO (Dawley1998),hasneverbeenstudied,probably The VNO of one male adult corn snake was dissected,
becauseoftheunder-representationofnonmammalianspe- fixed in 4% paraformaldehyde (PFA) overnight at 4(cid:2)C, and
ciesamongfullysequencedvertebrategenomes(Milinkovitch decalcified in 0.5M ethylenediaminetetraacetic acid (EDTA)
andTzika2007;Milinkovitchetal.2010).TheVNOofsnakes for 3 days at 4(cid:2)C. Seven-micrometer paraffin sections
isassociatedwithasophisticatedtonguedeliverysystem:The were stained with hematoxylin and eosin using standard
tongue collects chemicals in the environment by means of procedures.
tongueflickingandtransfersthemtothevomeronasalopen-
ingsonthepalate.Severalneurologicalandbehavioralstudies
demonstratedthatsnakeslargelydependontheirvomerona- LabelingofVNOEpitheliumbyRetrogradeNeuronal
salsystemforpreytrailing,capture,andconsumption,aswell Tracing
asforcourtshipandshelterselection(MillerandGutzke1999; The head of a juvenile male corn snake was fixed in 4%
ZuriandHalpern 2003; Huang etal. 2006). Componentsof PFA for 4h at 4(cid:2)C. The skull was opened and
the vomeronasal receptor signal transduction pathway, 1,10-dioctadecyl-3,3,30,30-tetramethylindocarbocyanine per-
including several G proteins and the secondary messenger chlorate (Dil) crystals (Invitrogen) were placed on top of the
triphosphoinositol (IP3), have been identified in the garter AOB and the main olfactory bulb (MOB). The head was
snakeVNO(Luoetal.1994).Inthesamespecies,stimulations embedded in 4% low-melting agarose and placed in 4%
by prey-derived chemoattractants generate both transient PFA at 37(cid:2)C for 6 weeks before it was cut sagittally. The
activation of neurons inthe vomeronasalsensory epithelium VNOwasdissected,embeddedinOCT(optimalcuttingtem-
(Cinellietal.2002)andanincreaseinfiringratesofindividual peratureembeddingmedium),andfrozenondryice.Sixteen-
neurons in the accessory olfactory bulb (AOB), that is, the micrometercoronalsectionswerepreparedfromcryoblocks,
projection site of vomeronasal sensory neurons in the brain counterstainedwithHoechst(Invitrogen),andmountedwith
(Jiangetal.1990). Vectashield(VectorLabs).
Recenttechnologicaladvancesallowfordenovosequen-
cingandassemblyoftranscriptomes(RNA-seq),thatis,even
in the absence of reference genomes (Gibbons et al. 2009; cDNAPreparationandDeepSequencing
Schwartz et al. 2010; Surget-Groba and Montoya-Burgos
Immediately after sacrifice of the animals, the VNO was dis-
2010; Tzika et al. 2011). Here, using RNA-seq and RNA in
sected and stored in RNA-protective solution (PrepProtect,
situ hybridization, we report on the diversity and evolution
MiltenyiBiotech).FollowingtissuedisruptionwithaPolytron
of the corn snake (Elaphe guttata, recently renamed
device (Kinematica), mRNAwas extractedusing themMACS
Pantherophis guttatus) vomeronasal receptor repertoire. We
kit (Miltenyi Biotech). After verification of mRNA integrity
foundalargeV2Rrepertoirecomposedofsmall,butmultiple,
usingaBionanalyser(Agilent),mRNAwaseitherdirectlysub-
snake-specificandreptile-specificexpansionsandonlyafew
mitted to sequencing, or cDNA was synthesized and DSN
V1R receptor genes. We also show punctate expression of
normalized (i.e., enriched for less abundant transcripts using
V2Rs in the corn snake vomeronasal neuroepithelium, an
Kamchatka crab duplex-specific nuclease) using published
expressionpatterncompatiblewiththatobservedinmammals
procedures(Zhulidovetal.2004).NormalizedcDNAsamples
(Martinietal.2001;Silvottietal.2007;IshiiandMombaerts
were polymerase chain reaction (PCR)-amplified (30 cycles),
2011).Ourresultsrefutethehypothesisthatthedistribution
andthe50-polyTcDNAsynthesisadaptorwasremovedusing
of vomeronasal receptors in vertebrates is dictated by the
restrictionenzymes.RNAsamplesandnormalizedcDNAsam-
evolutionary transition from water to land (Shi and Zhang
ples, pooled or individually (supplementary fig. S1 and table
2007) and shed light on the evolutionary dynamics of V1R
S1,SupplementaryMaterialonline),weresubmittedtoRoche/
andV2Rfamiliesinvertebrates.
454(Macrogen)andIlluminasingle-endsequencing(libraries
werepreparedaccordingtostandardprotocols).Thesamples
Materials and Methods sequencedusingRoche/454werethreemales(M1,M2,and
M3), pooledandnotnormalized; threefemales(F1,F2,and
Animals
F3), pooled and normalized; and one juvenile female (F4),
Cornsnakes(Pantherophisguttatus)wereobtainedfromour normalized. The samples sequenced using Illumina were
in-house breeding colony and maintained according to the three females (F1, F2, and F3), pooled and normalized, and
Geneva Canton regulations (authorization 1008/3421/1R). onemale(M4),notnormalized.
390 GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013
GBE
SnakeVomeronasalReceptorRepertoires
ContigAssemblyandDatabaseSearches intoaasequencesandalignedusingtheMAFFTmultiplealign-
mentprogram(Katohetal.2002).Aftermanualeditingofthe
The two types of reads (Roche, 200–500bp, and Illumina,
alignment using JalView (Waterhouse et al. 2009), contigs
114bp) were assembled in two steps (see supplementary
spanning mostly the 30-UTR or with stop codons in the
fig.S1andtableS1,SupplementaryMaterialonline).
coding part of the transcript were removed. The remaining
setof196potentiallyfunctionalV2Rtranscriptsarelistedand
Step1:PreselectionAssembly
their sequences provided, in supplementary file S1, Supple-
Adaptersremovalandqualitytrimmingaswellasassemblyof mentaryMaterialonline.Sequenceswithstopcodonsinthe
all Roche/454 reads (from M1–3, F1–3, and F4) were per- coding sequence are listed and their sequences provided, in
formedwithSeqManNGenv.2(DNASTAR).Readsobtained supplementary file S2, Supplementary Material online.
fromIlluminaweresubjectedtoadapterremovalandquality Visualization of contig statistics and properties (fig. 2) was
trimmingusingtheFASTX-Toolkit(http://hannonlab.cshl.edu/ performed using the R software (www.r-project.org, last
fastx_toolkit,lastaccessedFebruary3, 2013) andassembled accessedFebruary3,2013).
using Velvet with the Oases extension (Zerbino and Birney
EstimationofTranscriptNumberandVariability
2008).Optimalk-merlengths(matchlengthusedbyVelvet)
maydifferfortranscriptswithdifferentabundances(Gibbons ToestimatethenumberofdistinctV2Rtranscriptsandtheir
et al. 2009; Surget-Groba and Montoya-Burgos 2010). To sequence variability, we used the 350 aa of the alignment,
maximizethechancesofidentifyingVRtranscripts,weapplied which correspond to the 9-cysteine domain and the 7TM
an additive k-method as described in Surget-Groba and domain. We divided the alignment into four segments
Montoya-Burgos(2010).K-mervaluesof51,61,71,and81 (fig.2D)and,withineachsegment,wecountedthenumber
wereused,followedbytheremovalofredundantcontigswith ofuniquesequencesamongthosewithlessthan50%missing
CD-HIT-EST(LiandGodzik2006).TheShortReadRpackage data. We aligned the corresponding nucleotide sequences
(www.r-project.org,lastaccessedFebruary3,2013)wasused and counted the number of unique nucleotide sequences.
foranalyzingIlluminadata(Morganetal.2009). Assembled Tocalculatesequencevariability,wedividedtheaaalignment
Roche and Illumina contigs and unassembled Roche reads into eight segments and computed pairwise identity among
(singletons)withlengthabove200bpwereannotatedusing sequenceswithnomissingdata.
BLASTX,implementedinLANErunner(Tzikaetal.2011),with
PCRValidationofAssembledContigs
the following criteria: minimal match length 30 amino acids
(aa),minimalsequenceidentity50%,maximal Evalue0.05, PCR primer pairs were designed for each of the 25 selected
and release 60 of the Ensembl protein database for the fol- contigs. cDNA samples from the VNO of one male and one
lowing species: Anolis carolinensis (anole lizard), Danio rerio femalecornsnakeswerepreparedasdescribedearlier,pooled
(zebrafish),Monodelphisdomestica(opossum),Musmusculus andusedastemplatesforPCR.Productsoftheexpectedsizes
(mouse),Ornithorhynchusanatinus(platypus),Rattusnorvegi- wereamplifiedfor24ofthe25PCRassaysandsubmittedto
cus (rat), and Xenopus tropicalis (western clawed frog). The capillarySangersequencing:23ofthe24sequencescouldbe
anole lizard sequence that we identified as V1R (see below) aligned to the original corresponding contig. Percentages of
wasincludedinthesearch.AssembledIlluminacontigsaswell nucleotide differences were calculated over the validated
asassembled(contigs)andunassembled(singletons>200bp) sequence length (supplementary fig. S2, Supplementary
Roche/454 sequences are available at http://www.reptilian- Material online). Sequences of the 23 primer pairs used for
transcriptomes.org (last accessed February 3, 2013) validation (amplification and sequencing) of the predicted
(supplementary files S6a and b, Supplementary Material sequences are listed in supplementary table S3, Supplemen-
taryMaterialonline.
online).
PotentialPythonV2RGenes
Step2:FinalAssembly
We performed Basic Local Alignment Search Tool (BLAST)
Reads contributing to contigs with a VR hit and singletons analysesusingournewlyidentifiedcornsnakeV2Rtranscripts
with a VR hit (in total 8,048 Roche and 54,191 Illumina asinput,aswellasV2Rsfromothervertebrates,todetectV2R
reads) were pooled and assembled using SeqMan NGen v.2 genesintheBurmesepython(Pythonmolurusbivittatus)draft
(DNASTAR). Assembled contigs were annotated using LANE genome(Castoeetal.2011).OurTBLASTXanalysesrevealed
runner(Tzikaetal.2011)asdescribedearlier.Althoughstrin- thepresenceof216partialsequencesoflikelyV2Rgenesthat
gencyofBLASTXsearchcriteriawaslow,alltheV2Rcontigs werealsocomparedwiththenonredundantNationalCenter
wereidentifiedwithaveryhighconfidence(E<10(cid:3)9).Read for Biotechnology Information (NCBI) database to confirm
alignmentsofpotentialV2Rcontigsweremanuallycorrected their annotation (supplementary file S3, Supplementary
forlonghomopolymererrorsintroducedbyRoche/454pyro- Material online). For each of these shotgun-sequenced frag-
sequencing.PotentialVRtranscriptfragmentsweretranslated mentsofthePythongenome,thelongestopenreadingframe
GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013 391
GBE
Brykczynskaetal.
(ORF)wasidentifiedusingtheORFfinderofNCBI(supplemen- queries against the garter snake database, two potential
tary file S4, Supplementary Material online). Among these, Python V1R sequences were identified as two independent
114 ORFs (longer than 140 aa) were selected and aligned shotgun-sequenced fragments of the species genome
withknownV2Rs.Thisalignmentindicatedthatthemajority (Python V1r1—GenBank: AEQU010364851.1 and Python
of the Python V2R sequences correspond to the variable N- V1rb—GenBank: AEQU010376814.1; from the whole-
terminusasonly34alignedwiththe7TMdomain. genome shotgun sequencing project AEQU000000000.1).
Searching our corn snake VNO transcriptome with TBLASTX
PhylogeneticAnalyses did not reveal any V1R homologs. On the other hand, two
slightlydifferentcopiesofeachPythonV1ra1andV1r2were
aa sequences of the 7TM region (265 aa) of 66 V2R corn
successfully amplified from Pantherophis guttatus genomic
snake contigs spanning at least 60% of the 7TM region
DNA using primers based on the Python sequences (supple-
werealignedwith67representativeV2Rsequencesofzebra-
mentaryfileS7,SupplementaryMaterialonline).Thedistances
fish, frog, anole lizard, mouse, and opossum using MAFFT
(Katoh et al. 2002), followed by manual adjustment using ofthetwocopiestoeachotherandtoseveralvertebrateV1Rs
JalView (Waterhouse et al. 2009). An alignment with the were computed using MetaPIGA-v2.1 (Helaers and Milinko-
additionof15PythonORFswasalsoperformed.GTR20was vitch2010),bothforthenucleotide(withJukes–Cantorcor-
selected as the best aa-substitution model using the Akaike rection)andtheaasequences(withPoissoncorrection). The
InformationCriterionimplementedinMetaPIGA-v2.1(http:// lizard and snake V1Rs thus identified were translated and
www.metapiga.org(lastaccessedFebruary3,2013)[Helaers aligned with their mouse, fish, and frog homologs. The
andMilinkovitch2010]).Maximumlikelihood(ML)treeswere poorlyalignedregionsofthe438-aaMAFFTalignmentwere
computed using the Metapopulation Genetic Algorithm trimmedwiththe“strict”criterionoftrimal(Capella-Gutierrez
(MetaGA [Lemmon and Milinkovitch 2002]) implemented in et al. 2009) implemented into MetaPIGA-v2.1 (Helaers and
MetaPIGA-v2.1 (Helaers and Milinkovitch 2010). We used Milinkovitch 2010), and a final alignment of 230 aa per
probability consensus pruning among four populations of sequence was phylogenetically analyzed as described earlier
four individualseachandestimatedposteriorprobability dis- fortheV2Ralignment.
tribution of possible trees by performing replicated metaGA
searches and stopping when the mean relative error values
mRNAInSituHybridization
among10consecutiveconsensustreesremainedbelow2%.
ML trees were also inferred using RaxML (v7.2.6) with 500 Templatesforprobeswereamplified(primerslistedinsupple-
rapidbootstrapreplicates(Stamatakis2006).Finally,MC3ana- mentary table S3, Supplementary Material online) from a
lyses under a Bayesian framework were performed using mixed sample of VNO cDNA from male and female corn
MrBayes(v3.2)(Huelsenbecketal.2001),andposteriorprob- snakes. Digoxigenin (DIG) and fluorescein (FLUO) probes
abilitieswereestimatedafter2milliongenerations(burnin¼ weresynthesizedaccordingtotheDIGandFLUORNAlabeling
7,000 trees sampled every 100 generations). The RaxML kitsprotocols(Roche).Sixteen-micrometercryosectionswere
topologyisshowninfigure3(brancheswithbootstrapsup- prepared from freshly frozen VNO of male and female corn
portbelow50%arecollapsed),andRaxML,MetaPIGA,and snakes. Sections were fixed in 4% PFA for 20min at room
MrByesbranchsupportsareindicatedabovethebranchesof temperature.Single-anddouble-staininginsituhybridizations
majorclades. were performed using published protocols (Riviere et al.
2009). Hybridizations were carried out at 62(cid:2)C for 14h.
IdentificationofLizardandSnakeV1RSequences Anti-DIG antibody coupled to alkaline phosphatase (AP)
RepresentativeV1Rnucleotidesequencesofzebrafish,opos- and anti-FLUO antibody coupled to horseradish peroxidase
sum, mouse, platypus, rat, and frog were used as input of (POD) were used (Roche). FastRed (Sigma) and biotinyl-tyra-
TBLASTXqueriesagainsttheanolegenome(Ensemblrelease mide(PerkinElmer)wereusedassubstratesforAPandPOD,
64).Thisledtotheidentification(by“best-reciprocalhit”)ofa respectively. Biotinyl-tyramide was detected by streptavidin-
single Anolis sequence (ENSACAG00000025768) as a V1R conjugatedAlexaFluor488(Invitrogen).Sectionswerecoun-
homolog. The presence of the corresponding sequence in terstained with Hoechst (Invitrogen) and mounted with
the genomic DNA of the anole lizard was confirmed by Vectashield (Vector Labs). Low-magnification images were
PCR,followedbySangersequencing.KnownV1Rnucleotide taken with a standard fluorescent microscope (Zeiss). High-
sequencesofvariousvertebratesalongwiththenewlyidenti- magnificationimagesweretakenwithaconfocalmicroscope
fied Anolis V1R gene were also compared (TBLASTX using (Leica).Thenumberofcellsexpressingreceptorsandthetotal
LANE runner) with the fully sequenced genome of the number of cells in the sensory epithelium were manually
Burmese python (Python molurus bivittatus), as well as the counted and used to calculate the percentage of expressing
multiorgan transcriptome of the garter snake (Thamnophis neuronsinsupplementaryfigureS5,SupplementaryMaterial
sirtalis). Although no significant hit was retrieved from the online.
392 GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013
GBE
SnakeVomeronasalReceptorRepertoires
DataAvailability partofthesequence(the7TMandcysteinedomains)(supple-
mentaryfileS1,SupplementaryMaterialonline).Wealsoiden-
All developed tools and supplementary files, Supplementary
tified 39 contigs representing putative pseudogenes as they
Materialonline,areavailableathttp://www.reptilian-transcrip
included stop codons leading to a translation termination
tomes.org(lastaccessedFebruary3,2013).
beforetheendoftheseventhtransmembranedomain(sup-
plementary file S2, Supplementary Material online). We did
Results and Discussion
not observe any large insertions or deletions in any of these
TheCornSnakeVNO putative V2Rs or pseudogenes. To validate our results, we
selected 25 contigs for PCR amplification and sequencing.
The corn snake, suggested as one of the most convenient
Among these, 92% (23/25) could be amplified, and their
reptilian model species (Milinkovitch and Tzika 2007), was
sequenceswereidenticalornearlysowiththecorresponding
selectedfortheinvestigationofthevomeronasalchemorecep-
target sequence (average percentage of nucleotide differ-
tor repertoire of snakes (fig. 1A and B). First, to identify the
ences¼0.012;supplementaryfig.S2,SupplementaryMaterial
sensorypartofthevomeronasalepitheliuminthecornsnake
online).Thesefewdifferencesbetweentargetandamplified
VNO,weretrogradelylabeledsensoryneuronsbyapplyinga
lipophilicneuronaltracer(Godementetal.1987)ontheAOB sequences were attributed to transcriptome sequencing
ofafixedbrain.Figure1C–Eindicatesthatthemajorityofcells errors,amplificationofcloselyrelatedgenefamilymembers,
present in the vomeronasal epithelium are neurons whose or polymorphisms between individuals used for sequencing
axons project to the AOB and whose dendrites protrude andvalidation.Takenasawhole,ourvalidationdemonstrates
intothelumenoftheVNOcavity.ThecornsnakeVNOneu- that the majority of assembled contigs represent genuine
roepithelium is organized in columns (fig. 1B and E), like in transcripts.
othersnakespecies(WangandHalpern1980;Taniguchietal. BoththemeanandmedianlengthsofsnakeV2Rcontigs
2000). The unlabeled cells in the basal zone of the were approximately 500bp (fig. 2A). Hence, most of them
epitheliumlikelyrepresentneuronalstem cells.Attheedges represent a fragment of a full transcript (the approximate
of the epithelium, we observe groups of labeled cells that mean length of a V2R coding sequence in vertebrates is
possiblyrepresentapoolofimmatureneurons,astheirden- 2,100bp). V2R sequences corresponding to the N-terminus
driticprojectionsintothelumenarenotlabeled(arrowheads, areunder-representedinourdatasetforatleasttworeasons:
fig.1D). 1) The low evolutionary conservation of the V2R N-terminal
region makes homology assignment difficult using BLASTX
TheSnakeVNOChemoreceptorTranscriptome searches against the genomes of other vertebrates and 2)
cDNA libraries tend to be biased against 50 sequences. To
We isolated vomeronasal mRNA from four male and four
estimatetheminimumnumberofdistincttranscriptspresent
femaleindividuals,ofwhichallbutonefemaleweresexually
inourassembly,wefirsttranslatedthe196putativeV2Rsnake
mature.SomeoftheproducedcDNAsampleswerenormal-
transcripts into aa sequences. Next, we generated multiple
izedtoenrichforlessabundanttranscripts.Wesubmittedtwo
alignmentsofnucleotide(nt)andaasequenceswithMAFFT
normalized and two non-normalized samples to Roche/454
(Katohetal.2002).Wedividedthealignmentsintofourseg-
andIlluminasequencing(seesupplementarytableS1,Supple-
ments(of228–279ntor76–93aa).Ineachofthesesegments,
mentaryMaterialonline,forsequencingandassemblystatis-
wecountedthenumberofuniquesequences(fig.2D);note
tics). In total, we obtained 343,062 reads of 400bp mean
length (Roche/454) and 54,394,908 reads of 114bp theN-terminushighvariabilityofV2Rsequences.Thethirdand
(Illumina). The high number of Illumina reads allowed for fourthsegmentsofthealignment, locatedinthemorecon-
greater depth coverage, whereas the longer Roche/454 served 7TM domain, indicated that our 196 snake V2R
readsreducedthechancesofchimericassemblies.Following sequencesexpressedinthecornsnakeVNOcorrespondtoa
a two-stepapproach(supplementary fig.S1, Supplementary minimum of 116 distinct transcripts, encoding 109 different
Material online), we first generated contigs separately from proteins.
thetwotypesofreadsandidentified(usingBLASTXagainstall The suggestion of a large and variable V2R repertoire in
protein-coding sequences of seven vertebrate species) the snakesisalsosupportedbyourTBLASTXanalyses(supplemen-
snake contigs likely to encode for vomeronasal receptors. tary file S3, Supplementary Material online) against 1) the
Second,thereadscontributingtothesecontigswerepooled Burmese python (Python molurus bivittatus) draft genome
into a single set of 62,239 sequences that were assembled (Castoe et al. 2011), in which we identified 216 partial
into467mixed-readcontigs.Outofthese,196contigswere sequences of potential V2R genes and 2) a garter snake
identified as fragments of likely functional snake V2R tran- (Thamnophis elegans) multiorgan transcriptome (from brain,
scripts based on high-confidence BLASTX similarity hits gonad,heart,kidney,liver,spleen,andbloodtissuesbutnot
(E<10(cid:3)9) against V2R transcripts from other species (fig. olfactorytissue[Schwartzetal.2010])inwhichweidentified
2C) and the absence of stop codons in the protein-coding 38potentialV2Rgenes.
GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013 393
GBE
Brykczynskaetal.
FIG.1.—TheVNOofthecornsnake.(A)Cornsnake.(B)CoronalsectionofacornsnakeVNOstainedwithhematoxylinandeosin(HE).Fullorgan(left),
scalebar200mm;close-upsofthecolumnarsensoryepithelium(topright)andthezoneincontactwiththeoutsideworld(bottomright),scalebars20mm;
SE,sensoryepithelium;NE,nonsensoryepithelium;L,lumen.(C)Schematicrepresentationofaheadhemi-sectionwithapictureofthefluorescencedetected
intheVNOandthemainolfactoryepithelium(MOE)afterapplicationofaretrogradetracingdyeontotheAOBandmainolfactorybulb(MOB).(D,E)Coronal
sectionoftheVNOwithfluorescence(red)detectedafterapplicationofthetracingdyeontotheAOB.Arrowheadsandarrowsin(D)indicategroupsof
potentiallyimmatureneuronsattheedgesandbase,respectively,oftheneuroepithelium;arrowsanddottedframein(E)indicatedendriticprojectionstothe
lumen.DNAisstainedwithHoechst(blue).Scalebarin(D),200mmandin(E),50mm.
Reptile-SpecificExpansionofV2Rs estimation performed with MrBayes (Huelsenbeck et al.
2001).Thegene-familytreeincludesthethreewell-described
ToinfertheevolutionaryhistoryofsnakeV2Rs,weperformed
mammalianV2RsubfamiliesA,B,andD,aswellasthemore
phylogeneticanalysesofsnakeV2Raasequenceswithhomo-
distantly related subfamily C (Yang et al. 2005; Grus et al.
logsfromvariousvertebratespecies.Amultiplealignmentwas
2007; Young and Trask 2007). Our analyses indicate that
built among a selection of 67 sequences from Anolis caroli-
nensis,Musmusculus,Monodelphisdomestica,Xenopustro- family C is also present in nonmammalian vertebrates but is
picalis, and Danio rerio, as well as 66 of our corn snake representedbyasinglememberbothintheanolelizardand
sequences (only sequences that spanned at least 60% of thecornsnake.TheV2Rgenetreealsoincludesalargereptile-
the 7TM domain were used; supplementary fig.S3, Supple- specific clade whose topology indicates that 1) a significant
mentaryMaterialonline).Notethatthemajorityofthesecorn portionofthesnakeV2Rrepertoirearosebeforethesplitof
snakesequenceshaveanaveragedepthcoverageabove10 thislineagefromothersquamates(i.e.,severalmonophyletic
(fig. 2B). The consensus topology shown in figure 3 is sup- groupsincludebothcornsnakeandanolesequences)and2)
portedbythreedifferentheuristicsforlargephylogenyinfer- snakes experienced further evolution of their V2R repertoire
ence: the metapopulation genetic algorithm (Lemmon and butintheformofmultiplesmallexpansions(ofwhichfiveare
Milinkovitch 2002) implemented in metaPIGA-v2 (Helaers indicatedinfig.3).GiventheshortsizeoftheV2Ralignment
and Milinkovitch 2010), the rapid bootstrap analysis carried (265aa/sequence),increasingthenumberofsequencesinthe
out with RaxML (Stamatakis 2006), and MC3 Bayesian alignment will tend to reduce the robustness of branches
394 GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013
GBE
SnakeVomeronasalReceptorRepertoires
A B C
D E
FIG.2.—TheV2Rrepertoireofthecornsnake.(A)Lengthdistributionof196cornsnakeV2Rcontigs.(B)Distributionoftheaveragereaddepthcoverage
inrelationshiptothelengthof196cornsnakeV2Rcontigs(thoseusedforphylogeneticanalysesaremarkedinred).(C)DistributionofBLASTXEvaluesin
relationshiptothelengthof196snakeV2Rcontigs.(D)Identitystatistics(basedonsegmentsofamultiplealignmentof196snakeV2Rcontigs)foraaandnt
sequencevariabilityanalysesandestimationofthenumberofdistinctV2Rtranscripts;seetextfordetails.(E)AbsoluteandrelativesizesofVRgene
repertoires,basedonthepresent(redframe)andpublishedstudies(Grusetal.2007;ShiandZhang2007;YoungandTrask2007;Date-Itoetal.2008;
Hashiguchietal.2008;Jietal.2009;Youngetal.2010;Alfoldietal.2011).
duringphylogenyinference.Weanywaybuiltasecondalign- and EG154). We also analyzed the expression pattern of
mentwiththeadditionof15PythonputativeV2Rfragments one identified pseudogene (V2R-EG204). At least three
spanningatleast60%ofthe7TMdomain.Phylogeneticana- snakes were analyzed for each receptor gene, and similar
lysesofthislargerdataset(supplementaryfig.S4,Supplemen- patternsofexpressionwereobservedamongindividuals(sup-
tary Material online) yield a topology highly similar to that plementaryfig.S5,SupplementaryMaterialonline).Alltested
shown in figure 3, with Python sequences clustering with transcriptsweredetectedinbothmalesandfemales.Foreach
corn snake sequences (including one Python sequence in ofthefiveV2Rgenes,wedetectedexpressioninthevomer-
family C; supplementary fig. S4, Supplementary Material onasalsensoryneuroepithelium(fig.4)andnotintheMOE.
online). Similar to what has been reported for mice (Herrada and
Itismostlikelythattheclearseparationbetweenmamma- Dulac 1997; Matsunami and Buck 1997; Ryba and Tirindelli
lian, frog, and squamate non-C clades reflects phylogenetic 1997), we observed punctate signals (e.g., V2R-EG107 and
signalratherthanlong-branchattractionor homogenization EG154) corresponding to single neuron expression, with a
through gene conversion (Mallon et al. 2004; Ezawa et al. cytoplasmic localizationofV2RmRNA. Asobservedinother
2006). The latter hypothesis could be investigated through vertebratespecies,theprobabilityofexpressionofagivenV2R
analysis of synteny conservation and patterns of sequence fromnon-Csubfamiliesisgenedependent.Thisprobabilityis
divergenceamongcloselyrelatedparalogs(Sawyer1999). high for V2R-EG108 and lower for V2R-EG107/154/204
(figs.4and5).Wenotethatwecannotexcludethatdistinct,
ExpressionPatternsofV2RSubfamiliesAreConserved
butverycloselyrelated,receptortranscriptswerecodetected
amongVertebrates
bysomeprobes.
TheexpressionofsnakeV2Rgeneswasinvestigatedwithin One striking feature of V2Rs in the mouse VNO is that
situhybridizationsoncornsnakeVNOsections(fig.4)using subfamily-C members are broadly expressed and are coex-
probescorrespondingtoputativelyfunctionalV2Rtranscripts pressed in the same cells with non-C (ABD) V2Rs (Martini
(indicatedwitharrowsinfig.3:V2R-EG001,EG107,EG108, et al. 2001; Silvotti et al. 2007; Ishii and Mombaerts 2011).
GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013 395
GBE
Brykczynskaetal.
FIG.3.—EvolutionaryhistoryofsnakeandothervertebrateV2Rs.MLtree,basedonamultiplealignmentof7TM-domainaasequencesfrom66corn
snakeV2Rsand67representativeV2Rsequencesfromfivevertebratespecies.Branchsupportvalues(underRaxML/MetaPIGA/MrBayes)areindicatedfor
majorclades.Brancheswith<50%supportvaluesarecollapsed.Pinkboxeshighlightfiveofthesnake-specificexpansions.Arrowsindicatesequencesused
forinsituhybridization(figs.4and5).Thetreeisrootedwithtastereceptors(TAS1Rs)asanoutgroup.Scalebar:Meannumberofaasubstitutionspersite.
The functional relevance of this peculiar pan- and coexpres- coexpressed with members of non-C subfamilies, such as
sion pattern, if any, is unknown but could reflect that these V2R-EG154(fig.5B).
receptors work as heterodimers with one monomer widely
expressed and the other specifically expressed. Interestingly, V1RversusV2RReceptorRepertoires
wefindthatthecornsnakesubfamily-Cgene(V2R-EG001)is No V1R sequences were found in our snake vomeronasal
expressedinthemajorityofthevomeronasalsensoryneurons transcriptome,evenunderrelaxedcriteriaofBLASTXsimilarity
(fig. 4 and supplementary fig. S5, Supplementary Material searches. To further investigate the potential lack of V1R
online). To test for monogenic expression of non-C and C genes in corn snakes and other Sauropsida reptiles, we per-
V2Rs, we performed double-staining in situ hybridizations. formed TBLASTX searches against the fully sequenced gen-
Figure 5A indicates that V2R-EG154 and V2R-EG107, that omesoftheanolelizard(Ensemblrelease64)andtheBurmese
is, two close members of the non-C family of snake V2Rs python (Castoe et al. 2011), as well as the multiorgan tran-
(fig. 3), show mutually exclusive transcription. In contrast, scriptomeofthegartersnake(Schwartzetal.2010).Although
the snake gene belonging to subfamily C (V2R-EG001) is no lizard V1R gene is predicted in the Ensembl database
396 GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013
GBE
SnakeVomeronasalReceptorRepertoires
A
B
FIG.5.—MonogenicversusnonmonogenicexpressionofV2Rs.Insitu
hybridizations with antisense RNA probes for (A) two ABD-subfamily
members V2R-EG154 (green) and V2R-EG107 (red), (B) one C (V2R-
EG001;red),andoneABD(V2R-EG154;green)subfamilymembers;coex-
pression of V2Rs is only observed between family C and family A/B/D
members. Left and middle panels show single channels in green and
red,andrightpanelshowsthemergeplusHoechststainingofDNAin
blue.Scalebar,10mm.
V1r1–4)correspondtotwoallelesorcloseparalogsforeach
ofthetwoV1Rs.
Wealsoperformedphylogeneticanalysesofa230-aadata
setincludingrepresentativesofthe12V1Rmousesubfamilies
(Rodriguez et al. 2002), the 19 known frog V1Rs (Date-Ito
et al. 2008), V1Rs from teleost fish (Pfister et al. 2007), the
detectedanolelizardV1Ralongwithadditionalfishandfrog
homologs(asprovidedbyEnsemblrelease66),andthePython
andcornsnakeV1Rsidentifiedabove.TheresultingMLphy-
logenetictreeindicatesthatthesnakeV1Rs(alongwiththeir
FIG.4.—ExpressionofV2Rtranscriptsinthesnakevomeronasalsen- orthologs from lower vertebrates) belong to two distant
sory epithelium. In situ hybridizations (green) of coronal sections of a
lineagesofparalogs(fig.6),oneofwhichincludestheanole
femalecornsnakeVNO,withantisenseRNAprobesforonepseudogene
V1r.Hence,thereptilianV1Rsarenotreptilianacquisitions.In
(V2R-EG204)andfourlikelyfunctionalV2Rtranscripts.DNAisstainedwith
situ hybridizations of corn snake VNO sections revealed
Hoechst(blue).Scalebar,20mm.
vomeronasal transcripts of Pantherophis V1r3 and 4 that
were absent or present depending on the individual exam-
ined, whereas no Pantherophis V1r1 or 2 transcription was
(whereas 39 genes are annotated as V2Rs), our analyses
detected.Thereasonforthisvariabilityisunknownbutmay
recognized as a V1R a single anole previously unidentified becorrelatedwithseasonalexpressionbecausepositivestain-
protein-codingsequence.Weconfirmedthepresenceofthis ing was observed during the mating period but not during
sequence in the anole lizard genome by PCR amplification hibernation (supplementary fig. S6, Supplementary Material
followedbysequencing.Ourinsilicoanalysesalsoidentified online). Additional experiments are required to statistically
twopotentialV1Rs(calledherePythonV1r1andV1r2)inthe determinewhetherthisistrue.
Python draft genome, and their presence was confirmed by Furthermore, we searched our snake VNO transcriptome
PCR amplification from Python genomic DNA. On the other data sets for other GPCR classes. We did not find any
hand, no significant match was obtained by searching the sequenceswithsignificantsimilaritytotraceamine-associated
multiorgan garter snake transcriptome. We then designed receptors (TAARs) (Liberles and Buck 2006) or FPRs (Riviere
PCR primers based on the Python V1R sequences that were etal.2009)bothofwhichareexpressedinthemurineolfac-
used to probe the corn snake genomic DNA. This approach tory system. We, however, found several olfactory receptor
led to the identification of four V1R genes. Based on the (ORs) transcripts (data not shown) that may have originated
distances among corn snake and Python sequences (supple- fromaminorcontaminationofoursamplesbysnakeMOE.To
mentary table S2, Supplementary Material online), it is likely testthispossibility,weusedaPantherophisORprobe(Locus_
that these four corn snake sequences (called Pantherophis 18943—supplementary file S5, Supplementary Material
GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013 397
GBE
Brykczynskaetal.
FIG.6.—EvolutionaryhistoryofsnakeandothervertebrateV1Rs.MLtree,basedona230-aamultiplealignmentofreptilian,mouse,frog,andfishV1Rs.
ThetopologyoftheRaxMLanalysesisshown,andbranchsupportvalues(underRaxML/MetaPIGA/MrBayes)areindicatedformajorclades.Brancheswith
<50%supportvaluesarecollapsed.Thetreeisunrootedandhasbeenarbitrarilyorientedfordisplaypurposes(thereisnooutgroup).Thedashedarrow
indicatestheAnolissingleV1R.Scalebar:Meannumberofaasubstitutionspersite.
online),correspondingtoatranscriptthatwasidentifiedinour Ourdeep-sequencinganalysessuggestthatthecornsnake
sequence data set, and found punctate transcription in the genome contains more than 116 V2R genes. Although we
MOEandnotranscriptionintheVNO(datanotshown),con- cannotexcludethatsomereconstructedtranscriptscombine
sistentwiththeideaofMOEcontamination. polymorphisms among individuals, this number is likely to
Takentogether,ourresultsstronglysuggestthatsquamate reflect a minimal estimation of the snake V2R repertoire
reptiles (snakes and lizards) have a very limited repertoire of because 1) we identified 196 high-quality partial V2R
functionalV1Rs.ThismakestheratioofsquamateV1Rversus sequences that potentially represent distinct transcripts (thus
V2R sequences very similar to the ones observed in amphi- making the snake V2R gene repertoire the second largest
bians andteleosts, incontrast to mammals whose genomes known after Xenopus), 2) we identified 216 potential V2R
encodemoreV1RsthanV2Rs(fig.2E). genes in the Python draft genome, 3) our transcriptome
approach, contrary to genomic surveys, might have missed
poorly transcribed and temporally regulated VR genes, and
Conclusions
4)someVRgenesmaybeexpressedinnonolfactorystructures
Wereportontheidentificationofthevomeronasalreceptor andbeabsentfromolfactoryneurons.Ourphylogeneticana-
repertoireinsnakes(i.e.,thelineagewiththemostdeveloped lyses indicate the presence of reptile-specific and also of
VNO among all vertebrates), as well as the analysis of its snake-specific clades of V2R genes (fig. 3) that may reflect
expression pattern and evolutionary history. Our analyses specializationofspecificsnakereceptorstoparticularligands.
demonstrate that 1) snakes exhibit a large V2R repertoire This specialization is, however, unlikely to be a rule because
but a very limitednumber of V1Rgenesand 2) the peculiar electrophysiological recordings in garter snakes showed that
V2RexpressionpatternobservedintheVNOofmice(mono- individual neurons of the vomeronasal sensory epithelium
genicexpressionforV2RABD-subfamilymembersandbroad respondtomultipleclassesofstimuli,includingdifferentpep-
expressionofC-subfamilymembers)isconservedinsnakes. tides and aa (Inouchi et al. 1993). The snake vomeronasal
398 GenomeBiol.Evol.5(2):389–401. doi:10.1093/gbe/evt013 AdvanceAccesspublicationJanuary24,2013
