Signature of adaptive evolution in olfactory receptor genes in Cory’s Shearwater supports molecular basis for smell in procellariiform seabirds

Olfactory receptors (ORs), encoded by the largest vertebrate multigene family, enable the detection of thousands of unique odorants in the environment and consequently play a critical role in species survival. Here, we advance our knowledge of OR gene evolution in procellariiform seabirds, an avian group which relies on the sense of olfaction for critical ecological functions. We built a cosmid library of Cory’s Shearwater (Calonectris borealis) genomic DNA, a model species for the study of olfaction-based navigation, and sequence OR gene-positive cosmid clones with a combination of sequencing technologies. We identified 220 OR open reading frames, 20 of which are full length, intact OR genes, and found a large ratio of partial and pseudogenes to intact OR genes (2:1), suggestive of a dynamic mode of evolution. Phylogenetic analyses revealed that while a few genes cluster with those of other sauropsid species in a γ (gamma) clade that predates the divergence of different avian lineages, most genes belong to an avian-specific γ-c clade, within which sequences cluster by species, suggesting frequent duplication and/or gene conversion events. We identified evidence of positive selection on full length γ-c clade genes. These patterns are consistent with a key role of adaptation in the functional diversification of olfactory receptor genes in a bird lineage that relies extensively on olfaction.

The sense of olfaction is used by vertebrates to process information about their surroundings, playing a critical role in survival and reproduction as it facilitates the recognition of suitable food, predators and prey, mates, and territories 1,2 . Odour perception begins with the binding of odorant molecules to G protein-coupled receptors (GPCRs) primarily expressed in the nasal olfactory epithelium 3 . These receptors, encoded by a complex multi-gene family containing large numbers of olfactory receptor (OR) genes 4 , allow discrimination of a vast array of unique odorants.
Phylogenetic analyses of olfactory receptor genes in model vertebrates enabled the identification of nine different monophyletic classes of receptors (α, β, δ, ε, γ, ξ, η, θ and κ) divided into Type I and Type II 5,6 ( Fig. 1). Type I is the most diverse, and includes Class I (α, β, δ, ε, ξ and θ groups) and Class II (γ). Particularly the α and γ groups have experienced extremely large copy number expansions and because they occur in tetrapods they are thought to mediate the recognition of airborne odorants 5,7 (Fig. 1). The remaining Type I and Type II (which only includes group η) are mainly present in fish and amphibian genomes and are thought to detect water-soluble odorants 5 . The avian species surveyed in greater genomic depth (the chicken, Gallus gallus and the zebra finch Taeniopygia guttata) mostly have OR genes belonging to group γ, but also a few α and θ 6,8,9 . Comparative studies have shown that OR gene copy number can vary by several orders of magnitude among vertebrate groups, and even among closely related lineages [10][11][12] . Although the proportion of OR pseudogenes is characteristically high (20-60% 11 ), the relative numbers of intact, pseudo-(with an in-frame stop codon) www.nature.com/scientificreports www.nature.com/scientificreports/ The most comprehensive hybrid assembly of the OR-containing cosmids had a cumulative length of 7.4 Mb. It consisted of 399 unique scaffolds with a cumulative length of 2.6 Mb and 22 K degenerate contigs totaling 4.8 Mb (Supplementary Info). The final set of EVm models supported by ORFs consisted of 1,501 genes encoding proteins with a mean length of 135 amino-acids (range of 75-1215 amino acids). Of these, 317 were found in the assembly scaffolds and the remaining in the degenerate contigs.
We identified 220 Cory's shearwater full and partial length OR genes among the 1,501 gene models, ~ 60% of the OR gene number found in the genome assembly of the Northern fulmar, the only procellariidae with a draft genome available 6,40 , and 35% of the OR gene number found in the genomes of G. gallus and T. gutatta. Many additional small ORFs were discarded due to size (<75 amino acids), despite significant homology with other avian OR genes, so this is likely an underestimate of the total number of olfactory receptor genes in Cory's shearwaters.
In total, 45 OR genes were in the main scaffolds, where ORFs ranged between 297-936 bp (with mean length ± standard deviation of 678.3 ± 249.3 bp). The remaining 175 OR genes were found in the degenerate contigs, and the ORFs ranged in length between 225-945 bp (mean ± SD 394.9 ± 118.6 bp).
Structure and proportion of potentially functional shearwater oR genes. Of all OR genes identified, only 56 (~ 25%) were not at the edge of a contig. Of these, 20 (36%) were complete genes, 15 (27%) were pseudogenes and 21 (38%) were partial sequences (missing a start and/or stop codon and no in-frame stop codons). In all except four cases there was only one OR locus per scaffold, even though almost 80 scaffolds were above 10 Kb (Fig. S1b). In the scaffolds with more than one OR locus, ORFs were between 44 and 1,517 bp apart, and at least one was always a pseudogene. Of the 220 OR genes, 85% were annotated as putatively belonging to the γ-c clade (family 14 6 ) ( Table 1), 18 of which are complete. The remaining ORs belong to families already identified in birds and reptiles except for one gene, with highest sequence similarity to mammal genes (24-like family gene) ( Table 1).
Full-length γ-c shearwater OR genes include a single exon with a length of 312 amino acids. The logo view of the alignment suggests that the distribution of genetic variation is not uniform along the length of these genes (Fig. 2). The characteristic, conserved olfactory receptor-specific motifs are present, such as the MAYDRYVAIC in the transmembrane domain (TM) 3 -intracellular domain (IC) 2 boundary, the FSTC(LP)H at the end of IC3, and the conserved cysteine residues in the extra-cellular (EC) 1 and EC2 domains and prolines at TM7 (Fig. 2), which are important for the conformation stability of the receptor 41 . The most variable positions are found in the TM 4 and TM5 and to a lesser extent in TM3 and TM6 (Fig. 2).
phylogenetics of shearwater oR gene sequences. Phylogenetic analyses suggest that there are at least two distinct OR gene clades in the shearwater (Fig. 3a), in agreement with the functional annotation. A few OR genes cluster with those of the other sauropsid species in a γ clade that pre-dates the divergence of the different bird lineages. However, most chicken, zebra finch, fulmar and shearwater OR genes belong to the avian specific, well-supported γ-c clade, within which sequences cluster by avian family (Fig. 3a). The fulmar γ-c genes cluster with the shearwater genes in a paraphyletic Procellariidae clade, some clustering in a fulmar specific clade, suggesting that they are younger than the divergence of the two procellariids. When all 220 sequences are included in the analysis, some shearwater and fulmar γ-c sequences appear to be evolving at significantly faster rates than  Table 1. Functional annotation and location of the 220 Olfactory receptor genes found in the Cory's Shearwater genome. Genes that could not be reliable assigned to an OR sub-family with the functional annotation pipeline belong to an "unknown" family. *Genes were classified as Intact if they had start and stop codons, no premature stops or frameshift mutations; Partial if they had an incomplete coding region and were not at the edge of a contig; Pseudogene if they had the full coding region but had internal stop or frameshifts, Truncated if they were located at the edge of a contig.
www.nature.com/scientificreports www.nature.com/scientificreports/ those from the other species, as revealed by the very long branch lengths (Fig. 3a). This pattern may be due to a large number of pseudogenes evolving free of selective constraints, as it does not persist in analysis including only intact OR genes (Fig. 3b). In the analysis with the intact genes, the monophyly of the avian γ-c clade remains well supported as is the clade including all shearwater genes (bootstrap value of 92%).

positive selection in full length γ-c shearwater oR genes.
To determine the selection regimen governing the evolution of the intact γ-c shearwater genes, we considered two data partitions as GARD found evidence of one breakpoint at alignment position 246, which falls in the putative TM6 domain (Fig. 2). Based on the two-partition dataset of all full-length OR genes, the different methods found that the overall d N /d S was never significantly different from 1, varying between ω = 0.552 (SLAC) and ω = 1.044 (REL). However, positive selection was identified on individual site codons of the full-length shearwater γ-c genes, consistent with a key role of selection in the functional diversification of these genes in the shearwater (Fig. 4a, Table 2). The number of sites with evidence of positive selection varied with the method used, with amino-acid positions 109 and 274 (respectively in TM3 and TM7 domains) identified by all individual methods as well as by the integrative analysis ( Table 2). In addition, positive selection was identified in amino-acid positions 111, 114, 209, 210 and 219 by at least two methods. These five sites are also located in TM domains (Fig. 4a). Critically, amino acids 109 and 274 have side chains oriented towards the center of the receptor (Fig. 4b), suggesting that, in this species, they may be part of the odorant binding pocket.

Discussion
This work provides one of the first detailed molecular insights into the evolution of chemoreception in an avian group for which olfaction has been shown to play a vital role, by specifically studying sequence diversity and evolution of olfactory receptors in the Cory's shearwater genome.
Several studies suggest that mammalian and avian OR repertoires are influenced by life histories 6,10,12,18-20 . Procellariiform seabirds have anatomical characteristics that suggest a strong reliance on a functional olfactory system, such as one of the largest olfactory bulb size to brain size (OB) ratio among extant avian species and specialized tube-like nostrils 42 . It has been suggested that there is a positive correlation between OB ratio, OR gene repertoire and olfactory function 10,15,21,25,28,29 . In this context, one might expect shearwaters to have more OR genes than species with smaller OB ratio and less obvious reliance on olfaction for survival, such as the chicken and zebra finch, in which 674 and 688 isolated OR genes have been identified, respectively 9 . However, the procellariiform Snow petrel was found to have among the lowest estimated number of OR genes (n = 212) among nine surveyed avian species, despite its much larger OB ratio 25 . Likewise, the Northern fulmar was found to only have an average number of OR genes in a genomic survey that included 48 phylogenetic diverse avian lineages 6 . The evidence from Cory's shearwater is inconclusive in this regard. While we identified with certainty 220 OR genes (similar to the numbers found in other petrels but about a third of the total number estimated in the chicken), our study was limited by the maximum number of cosmid clones considered for sequencing, and possibly by the specificity of the probes used to screen the library, which were restricted to known OR motifs. Furthermore, putative OR sequences mapping to the edge of assembly contigs and for which the translated sequence was shorter than 75 amino acids were discarded. Finally, the 220 shearwater OR genes were identified in only a very small (albeit biased toward OR-containing) fraction of the genome (7.4 Mb) compared with the other species for which the whole genome was analyzed. On the other hand, we found 17 intact γ-c genes in the shearwater, a larger number than the 2 previously found in the fulmar, despite possibly having missed characterization of the full repertoire. www.nature.com/scientificreports www.nature.com/scientificreports/ This raises the possibility that a smaller olfactory receptor repertoire is compensated by a larger number of intact genes, or possibly a broader receptor binding pocket specificity through the use of other sensory receptors 27,43 . This issue should be revisited once a high quality assembled Cory's shearwater genome becomes available. www.nature.com/scientificreports www.nature.com/scientificreports/ A remarkable finding of this study is that the vast majority (85%) of the OR genes identified in the shearwater belong to the γ-c clade. If the proportion of γ-c genes among all identified shearwater OR genes is representative of their real frequency in the genome, then the olfactory receptor family in this species is mostly composed of young paralogous loci, reflecting a relatively recent expansion of highly similar genes putatively linked to the recognition of volatile odorants. Although the γ-c clade seems to be as abundant in the shearwater as in the chicken (63%) and the zebra finch (80%), it was found to be not as abundant in other avian species, including the Northern fulmar 6 , although low genome assembly quality may have compromised the ability to characterize large expansions of recently duplicated genes in some species 44 . In two falcon species the γ-c clade was found missing, possibly reflecting reliance on other sensory cues 26 . Nevertheless, the presence of this striking γ-c diversity in the shearwater genome, and its higher proportion of intact genes relative to what has been found in most species surveyed so far 16 , is an example of how ecological adaptation is likely contributing to an enhanced olfactory function in Procellariiformes.
Lineage-specific duplications and losses of olfactory receptor genes is characteristic of the evolution of the OR gene family in vertebrates and have resulted in a large range of OR gene repertoire size across lineages 5,10-12,20 . The co-occurrence of γ-c intact, partial and pseudogenes in the shearwater OR repertoire, sometimes in the same contig, agrees with the extreme dynamic nature of the evolutionary models proposed to govern the genes encoding OR and other vertebrate multi-gene families 14,44,45 . In the context of the olfactory receptor evolution, the presence of many paralogs facilitates the evolution of novel gene functions that allow the recognition of new odors, under a scenario of neo-functionalization, but may also result in large number of pseudogenes due to functional redundancy and consequent relaxation of selective constraints [45][46][47] .
A hallmark of the olfactory receptor structure is a highly variable pattern of amino acid sequence conservation determining the identity and specificity of odorants that are recognized by each receptor. Although it is possible that some ligand interactions will be mediated by conserved amino acids 48,49 , peaks of polymorphism have been consistently detected on the same transmembrane domains in different species, namely in the TM3, TM5, TM6 and TM7 16 , similar with what was found in this study. If these domains do indeed form part of the odorant binding pocket in seabirds, as it has been shown in other vertebrate lineages 22 , it is not surprising to find evidence of diversifying selection at these sites, since allelic diversity would confer a selective advantage by allowing individuals to recognize a wider spectrum of odors.
Lineage-specific patterns of OR gene diversity have been found to be tightly linked to the habitat and foraging ecology of the different groups 6,18,21 . Procellariiform seabirds are pelagic species that forage over vast areas of a seemingly featureless sea expanse. However, these species rely on extremely fine tuned sensory cues to recognize heterogeneity in the marine environment, ensuring foraging efficiency over areas of patchy productivity. Exactly what odors are used to build an odoriferous seascape remains to be studied, but the large expansion of the γ-c clade, with its high levels of polymorphism and underlying patterns of adaptive evolution may contribute to the necessary discriminatory power to be used in a genetic mechanism linking the olfactory system with navigation, foraging and homing and perhaps even with kin recognition and other critical social interactions in this and other avian systems.
Much remains to be learnt about the structure and functioning of sensory receptors in birds 6,11,12,16 , and the isolation and characterization of olfactory receptor loci in the Cory's shearwater is a significant step towards understanding the role of this fascinating multi-gene family in birds and, more generally, in the genomic underpinnings of ecological adaptation.

Materials and Methods
Sampling and oR gene isolation. A blood sample from a female Cory's shearwater was collected in Selvagem Grande Island, and kept in 95% ethanol until laboratory analyses. Genomic DNA (~1500 μg) was extracted with a phenol-chloroform standard protocol and used to build a cosmid library (built with SuperCos1 Cosmid Vector Kit, Stratagene, La Jolla, CA, USA, outsourced to America Pharma Source LLC). The library had inserts that varied in size between 30-42 Kb and a titer equivalent to 4X of estimated genome size. Since there is no genome size estimates for Procellariiformes, we assumed 1,36 Gb based on the average of all avian genome sizes available 50 , a value similar to the genome assembly size for Northern fulmar of 1,135Gb 40 . Cosmid clones were screened using Southern Blotting with three probes, two targeting primarily sequences from the γ-c OR clade and one sequences from the non γ-c OR clade 25 . However, it is noteworthy that all probes targeted highly similar, evolutionarily conserved coding motifs that are present across all OR gene families, as implied by the cross-hybridization results below. Probes were generated as described elsewhere 25 . A ~ 500 bp fragment was amplified in the shearwater, which was cloned and sequenced to confirm olfactory receptor identity (Supplementary Info). We picked 96 positive clones, approximately a third were positive for all three probes, and the remaining for either one or two probes.
Sequencing, assembly and annotation of oR gene-positive cosmids. The sequencing strategy aimed at generating long assemblies while minimizing costs and avoiding chimeras. Sanger sequencing generated ~ 900 bp long reads from cosmid ends which anchored each cosmid sequence. Depth of coverage was obtained from an equi-molar pool of non-tagged cosmids, using Illumina HiSeq 2000. 34 million Illumina 100 bp-long reads were used to build multiple assemblies at different depth of coverage using Velvet 51   www.nature.com/scientificreports www.nature.com/scientificreports/ pooled into ten subsets that minimized within pool clone similarity. For each pool, DNA was sheared to 3Kb and a MID 454 GS FLX Titanium Paired End library was built and sequenced. A total of 541,133 reads were obtained for the 10 pools, with an average read length of 311 bp, resulting in an additional 20X coverage of each of the 96 cosmids. Hybrid assemblies were built with the Celera Assembler v.7 52 using all 454 data and varying numbers of Illumina reads. Further analyses proceeded on the most comprehensive hybrid assembly.
All scaffold and degenerate contigs (the latter contain repetitive sequences that could not be unambiguously assembled with the main scaffolds) were processed through the Institute for Genome Sciences' eukaryotic structural annotation pipeline, which was configured to use a combination of tools, including ab initio gene finders, protein aligners and model combiners. Gene finders including GeneMark-ES (GM-ES) 53 , Augustus v2.7 54 and Maker 55 were run using chicken evidence data, and GeneID v1.4 56 and GlimmerHMM v3.0.1 57 used human and zebra finch informant data, respectively. Exonerate's protein2genome option 58 was run with a score cutoff of 500 and a cutoff of 75% to identify homologs to known chicken (n = 37,258 protein models) and zebra finch (n = 17,773) proteins. Finally, Evidence Modeler (EVm) r2012-06-25 59 was used to combine evidence from spliced protein alignments and ab initio gene predictions into weighed consensus gene models.
Proteins were searched with Hidden Markov Models (HMMs) generated with HMMbuild using default settings and calibrated with HMMcalibrate. One HMM was built from full length OR genes (mostly from G. gallus and T. guttata), and the second was based on an alignment of partial OR gene sequences from eight different avian species 25 . The HMMs were validated by searching UniProt (Universal Protein Resource) using HMMsearch (EMBL_EBI) with a cut off value of 1e −10 . A set of open reading frames (ORFs) was generated with ENSEMBL's getorf which was searched using the new HMMs, with default settings. ORFs and EVm models matching HMMs were compared and uniqued, and the assembly location of each of the resulting genes determined. Proteins generated by EVm were annotated with a custom functional annotation pipeline, implemented within a reusable Ergatis pipeline, where input proteins are searched using several tools in parallel including HMMer3 60 NCBI-BLASTP (cutoff 1e-6), SignalP 61 , ScanProsite 62 and TMHMM 63 . The HMM database used for HMMer3 is custom-made with TIGRFAMs, PFAM, among other custom records. Proteins which match the HMMs above trusted cut off levels inherit these annotations, or else BLAST results were used. phylogeny of shearwater oR genes. Phylogenetic trees were estimated for two different amino acid (AA) alignments, one including all OR shearwater sequences of length > 75AA (n = 220) (Supplementary Dataset 1) and the other only the full-length ORs (n = 20) (Supplementary Dataset 2). The alignments included complete OR sequences from chicken, zebra finch and green anole (Anolis carolinensis) from 9 and all non-redundant Northern fulmar (Fulmarus glacialis) proteins from Protein database of GenBank. We use a non-OR rhodopsin-like family GPCR gene (Adenosine receptor A2b from the chicken, zebra finch and emperor penguin (Aptenodytes forsteri), respective GenBank accession numbers NP_990418, XP_002198489, XP_009276449) as outgroup 21 . Sequences were aligned with MAFFT 7.12 64 , using E-INS-i parameters, and the alignments were manually curated in Mesquite v.2.75 65 . Approximate maximum-likelihood trees were reconstructed with FastTree v2.1.10 66 under the JTT evolutionary model 67 , implemented in the CIPRES Science Gateway 68 . Nodal support was evaluated with Shimodaira-Hasegawa (SH) tests with 1,000 bootstrap replicates 66 . estimates of positive selection on clade γ-c. The ratio of non-synonymous to synonymous substitutions (ω = d N /d S ) is used as an indicator of the strength of selection on protein-coding genes. The presence of recombination or gene conversion, frequent in multi-gene families, precludes the use of methods that estimate selection based on a single phylogenetic tree 69,70 . We estimated d N /d S per codon with four codon-based maximum likelihood methods available in HYPHY 71 , as implemented in the Datamonkey webserver 72 . We accounted for the possibility of recombination using the GARD module (also implemented in Datamonkey 72,73 to infer breakpoints in the alignment and generate tree topologies for each of the non-recombinant data partitions. We then used SLAC, FEL, REL 70 , MEME 74 methods, which differ in assumptions for the distribution of rates across sites or lineages, and an integrative approach, which considers all codon sites detected by each method, to infer selection within each partition. A significance level of P < 0.1 was used for SLAC, FEL and MEME, and a Bayes factor >50 for REL 70 . WEBLOGO (http://weblogo.berkeley.edu/logo.cgi, 75 ) was used to generate a sequence logo with the amino acid alignment of all γ-c complete shearwater olfactory receptor sequences, to facilitate visualization of sequence variation for the putative transmembrane (TM), intracellular (IC) and extracellular (EC) domains. The height of each letter corresponds to the relative frequency of each amino acid at a given position, reflecting sequence conservation.
Chimera 76 (http://www.cgl.ucsf.edu/chimera/index.html) was used to display the location of the polymorphic sites on the secondary structure of a shearwater complete γ-c OR. ethical statement. The volume of the blood sample collected from the female shearwater was 150 μl (approximately 0.2% vol/weight). Sampling guidelines were approved by the Instituto da Conservação da Natureza e da Biodiversidade (ICNB) and by the Parque Natural da Madeira (currently Instituto das Florestas e Conservação da Natureza), Portugal under permit: 2/2012S. Sampling further followed the requirements of the Directive 2010/63/EU of the European Parliament and of the council for the protection of animals used for scientific purposes.

Data availability
Cory's Shearwater cosmid library shotgun Illumina HiSeq and 454 data as well as the assembled scaffold and degenerate contigs and OR gene sequences will be deposited in NCBI under BioSample accession number SAMN13721836. For additional information contact corresponding author (MCS).