The dopamine receptor D₅ gene shows signs of independent erosion in toothed and baleen whales

Synteny analysis

To build the synteny maps for the DRD₅ gene locus in Cetacea and Hippopotamus amphibius (common hippopotamus) several annotated Cetacea genome assemblies were inspected and scrutinized using the NCBI browser, namely Orcinus orca (killer whale; GCF_000331955.2), Lagenorhynchus obliquidens (Pacific white-sided dolphin; GCF_003676395.1), Tursiops truncatus (common bottlenose dolphin; GCF_001922835.1), Delphinapterus leucas (beluga whale; GCF_002288925.1), Neophocaena asiaeorientalis asiaeorientalis (Yangtze finless porpoise; GCF_003031525.1), Lipotes vexillifer (Yangtze River dolphin; GCF_000442215.1), Physeter macrocephalus (sperm whale; GCF_002837175.1) and Balaenoptera acutorostrata scammoni (minke whale; GCF_000493695.1). Bos taurus (cattle; GCF_002263795.1), a fully terrestrial relative of extant cetaceans, was used as reference. Next, (1) in genome assemblies with annotated DRD₅, the following procedure was used: five protein-coding genes, upstream and downstream of the DRD₅ gene, and from the same strand, were collected; (2) if DRD₅ gene annotations were not present, the genomic locus was retrieved using Bos taurus (cattle) DRD₅ flanking genes as reference. The selected anchoring genes to search upstream and downstream DRD₅ neighboring genes were CPEB2 and OTOP1, respectively. Regarding Hippopotamus amphibius (common hippopotamus), the synteny map was built via BLAST searches against the assembled, fragmented and unannotated genome of the same species, available at NCBI (GCA_002995585.1). Bos taurus (cattle) DRD₅ flanking genes were used as reference; using the discontiguous megablast task from blastn, the best BLAST hit (highest alignment identity and query coverage) was retrieved and the coordinates of the alignment in the target genome carefully inspected. The Hippopotamus amphibius (common hippopotamus) synteny map was then built by sorting the genes according to the subject alignment coordinates within genes aligning at the same genomic scaffold.

Sequence retrieval and gene annotation

The NCBI “low-quality protein” (LQ) tag marks RefSeq sequences that have been automatically modified relative to the corresponding genome sequence to correct for possible ORF protein-altering indels or mismatches, assuming that these arise from sequencing errors or genome assembly artefacts. Yet, in several cases these frameshift and nonsense codons correspond to authentic biological modifications leading to shifts or truncations of the predicted ORF, culminating into an abnormal protein amino acid constitution. For this reason, to clarify the functional status of DRD₅ in species presenting LQ tag for this gene, the corresponding genomic regions were directly collected from NCBI.

Concerning species with no annotated DRD₅ genes two scenarios were identified: (1) annotated genomes excluding DRD₅ gene annotation and (2) unannotated genomes. Regarding the first case (i.e., the cetaceans Lipotes vexillifer, the Yangtze River dolphin and Balaenoptera acutorostrata scammoni, the minke whale, as well as other non-cetacean mammals), the DRD₅ genomic locus was retrieved using, as reference, annotated DRD₅ flanking genes from closely related species (for cetaceans, the terrestrial extant sister clade Bos taurus (cattle), regarding non-cetacean mammals see Table S1). If the genomic locus exhibited severe genomic fragmentation (presence of Ns), thus, hindering the retrieval using neighboring genes, blastn searches were conducted against the Whole Shotgun Contigs of the corresponding species via discontiguous megablast task, using as query the DRD₅ coding sequence (CDS) of the same reference species (Table S1). The genomic sequence corresponding to the BLAST hit with the highest alignment identity and query coverage was selected.

For species without annotated genomes (i.e., Balaenoptera bonaerensis (Antarctic minke whale), Eschrichtius robustus (gray whale), Balaena mysticetus (bowhead whale), Sousa chinensis (Indo-pacific humpbacked dolphin), as well as Hippopotamus amphibius (hippopotamus)) genomic sequences were retrieved through blastn searches in the corresponding genome assembly using the Bos taurus (cattle) DRD₅ CDS as query. For each species, the best genomic scaffold corresponded to the BLAST hit with the highest query coverage and identity value. Due to the presence of a fragmented genomic region in the DRD₅ gene annotation of T. truncatus (common bottlenose dolphin), the same BLAST search procedure was carried out for this species, using as target the Whole Genome Shotgun contig dataset of the same species.

Collected genomic sequences were further imported into Geneious Prime 2019 (www.geneious.com) and the DRD₅ gene CDSs manually annotated for each species. Briefly, using the built-in map to reference tool with the highest sensibility parameter selected, the reference single-exon DRD₅ gene, 3′ and 5′ UTR flanked, was mapped against the corresponding genomic sequence of the in-study species. Aligned regions were further carefully screened for ORF abolishing mutations including frameshift mutations and in-frame premature stop codons. For Cetacea and Hippopotamus amphibius (common hippopotamus) DRD₅ gene annotation, Bos taurus (cattle) DRD₅ was selected as reference. Regarding non-cetacean mammals DRD₅ annotation, different references were chosen according to the phylogenetic relationships between reference and test species (Table S1). Finally, the identified mutations were next validated using raw sequencing reads retrieved from at least two independent genomic NCBI SRA projects (when available). Briefly, blastn searches were conducted against the selected SRA projects (Supplemental Materials 1–3) using as query the nucleotide sequence containing the mutation. BLAST hits were imported into Geneious Prime 2019 and mapped against the manually annotated sequences, using the built-in map to reference tool, to confirm the presence of the identified mutation.

Phylogenetic and d_N/d_S analyses

The predicted cetacean and Hippopotamus amphibius (common hippopotamus) DRD₅ CDSs, as well as the DRD₅ CDSs of Homo sapiens (human) and Bos taurus (cattle) available at NCBI were translation aligned in Geneious Prime 2019 using the Blosum62 substitution matrix. The alignment (1,440 bp) was manually inspected and predicted DRD₅ CDSs involved in alignment suffered prior inspection, with frameshift insertions and deletions being omitted and stop codons recorded as missing. The alignment was posteriorly exported for phylogenetic analysis. A Maximum likelihood phylogenetic tree, performed in PhyML3.0 server (Guindon et al., 2010), was produced with the best sequence evolutionary model being determined using the built-in smart model selection (HKY85 +G) (Lefort, Longueville & Gascuel, 2017). Node support was inferred based on 1,000 bootstrap replicates. The resulting phylogenetic tree newick file is available as File S1 and was subsequently imported for visualisation in FigTree (Supplemental Material 4) (Rambaut & Drummond, 2012).

PAML 4.6 (Yang, 2007) was used to estimate the ratio (ω) of the nonsynonymous substitution rate (d_N) to the synonymous substitution rate (d_S) using the produced phylogenetic tree (Supplemental Material 4). Codeml with the branch model was implemented to estimate d_N/d_S ratios (ω) for seven different branch categories: Mysticeti, Odontoceti, Functional, a category corresponding to the common cetacean branch, other concerning the stem Mysticeti branch, an additional one regarding the stem Odontoceti branch (excluding Physeter macrocephalus, the sperm whale) and finally, a category corresponding to the Physeter macrocephalus (sperm whale) ancestor branch (Table 1). Mysticeti and Odontoceti branch categories comprised all predicted Mysticeti and Odontoceti DRD₅ sequences, respectively. The Functional branch category comprised Hippopotamus amphibius (common hippopotamus), Bos taurus (cattle) and Homo sapiens (human) DRD₅ CDSs. For each of the seven different branch categories, likelihood ratio tests (following χ² distribution) were conducted to compare the estimated ratio (ω) of the nonsynonymous substitution rate (d_N) to the synonymous substitution rate (d_S) determined by PAML (the alternative hypothesis), with the expected ratio (ω) of the nonsynonymous substitution rate (d_N) to the synonymous substitution rate (d_S) according to a neutral model of evolution (ω = 1) (the null hypothesis). All PAML analyses were run with the CodonFreq set to three and codon sites with ambiguous data (including gaps and missing data) were included in the analyses.

ORF-disrupting mutations of DRD₅ in Cetacea

For cetacean species with annotated genomes, we started by examining the DRD₅ gene locus, including neighboring genes, to verify and elucidate the orthology of the annotated and non-annotated genes and outline the genomic regions to be inspected (Fig. S1). All analyzed loci were found to be conserved, including in both Lipotes vexillifer (Yangtze River dolphin) and Balaenoptera acutorostrata scammoni (minke whale), which lacked previous DRD₅ gene annotations (Fig. S1). Subsequent manual annotation of all collected cetacean genomic sequences revealed DRD₅ gene erosion across all analyzed species, except in Balaenoptera acutorostrata scammoni (minke whale), for which the DRD₅ coding status could not be accessed due to fragmentation of the 5′ end of the respective genomic region (presence of sequencing gaps (Ns)). In detail, a conserved 2-nucleotide deletion was detected for the full set of Odontoceti examined species, except for Physeter macrocephalus (sperm whale) that presented a premature stop codon near the middle of the gene and a single nucleotide insertion close to the end of the gene (Fig. 1). The 2-nucleotide deletion alters the reading frame, leading to a drastic change in downstream amino acid composition. Additionally, non-conserved mutations were found in Lipotes vexillifer (Yangtze River dolphin), which presented two premature stop codons, N. a. asiaeorientalis (Yangtze finless porpoise) that presented a premature stop codon close to the middle of the ORF, and D. leucas (beluga whale) with a single nucleotide deletion near the 5′ end of the gene (Fig. 1). D. leucas (beluga whale) presumed DRD₅ sequence presented another noticeable feature. A massive and abrupt alignment identity decrease observed when aligning Bos taurus (cattle) DRD₅ gene against the genomic target region of D. leucas (beluga whale). The alignment identity drop was noted approximately in the middle of the complete alignment length for this species, suggesting that the DRD₅ gene sequence is interrupted, and further supporting pseudogenization in this species. Regarding Odontoceti, at least one ORF-abolishing mutation was validated using available genomic SRAs experiments for all studied species, excluding S. chinensis (Indo-pacific humpbacked dolphin) and Lipotes vexillifer (Yangtze River dolphin) for which no genomic sequencing runs were available at the NCBI SRA database (Supplemental Material 1).

Regarding the Mysticeti suborder, a conserved single nucleotide deletion was detected in all species except Balaenoptera acutorostrata scammoni (minke whale) (Fig. 1). A non-conserved 2-nucleotide insertion was found also found in Balaenoptera bonaerensis (Antartic minke whale) and an insertion of one nucleotide was detected at Eschrichtius robustus (gray whale) near the 5′ end of the DRD₅ sequence (Fig. 1). Again, at least one ORF-abolishing mutation was validated using genomic SRA experiments for all analyzed species (Supplemental Material 1). Importantly, no conserved mutations were detected between most Odontoceti, Physeter macrocephalus (sperm whale, Odontoceti) and Mysticeti lineages, suggesting that three pseudogenization events occurred independently after their evolutionary divergence (Fig. 1). To increase the robustness of our analysis, we further scrutinized the genome of the extant sister clade of the Cetacea, the Hippopotamidae and were able to predict a fully functional CDS for DRD₅ in Hippopotamus amphibius (common hippopotamus), supporting the loss of DRD₅ after Cetacea diversification.

d_N/d_S analyses support independent inactivation events of DRD₅ in Cetacea lineages

To further strengthen the mutational analyses, we next carried out a d_N/d_S approach. The d_N/d_S analyses for species with functional DRD₅ gene (Functional branches category) rejected the null hypothesis (ω = 0.09, p = 0.030, <0.05), confirming that, as expected, these have evolved under purifying selection (ω < 1), favoring the conservation of DRD₅ in the tested lineages (Table 1). In contrast, concerning Mysticeti and Odontoceti branch categories, the performed analyses failed to reject the null hypothesis (ω = 0.386, p = 0.692 for Odontoceti branches category and ω = 0.460, p = 0.291 for Mysticeti branches category, respectively) suggesting that DRD₅ evolved neutrally in both lineages, following gene inactivation events (Table 1). Moreover, evidences of purifying selection were found for the common Cetacean branch category (ω = 0.135, p = 0.001), suggesting the existence of a functional DRD₅ in the common ancestor to all cetaceans (Table 1). Purifying selection was also detected for the Odontoceti stem branch category (excluding Physeter macrocephalus, the sperm whale) (ω = 0.213, p = 0.045). In contrast the Physeter macrocephalus (sperm whale) ancestor branch category, does not present a d_N/d_S ratio (ω) statistically significantly different from 1 (ω = 0.320, p = 0.460) (Table 1). Together with the mutational evidences (Fig. 1), these results suggest two independent inactivation events within the Odontoceti. Regarding the Mysticeti stem branch category, selection analysis suggested neutral evolution (ω = 1.174, p = 0.905) (Table 1). However, mutational evidence supports a conserved inactivation event in the common ancestor of Mysticeti. Taken together, our results support at least two independent DRD₅ inactivation events within Cetacea.

Other mammalian species displaying DRD₅ LQ tags have a coding gene

Initial analysis revealed the presence of at least one ORF-abolishing mutation in Ovis aries (sheep), Phascolarctos cinereus (koala), Bison bison bison (plains bison) and Ochotona princeps (American pika). These are in some cases suggestive of gene inactivation and not sequencing artefacts (Emerling et al., 2017; Lopes-Marques et al., 2017). Manual annotation, including SRA validation, unveiled sequencing reads supporting the absence of disruptive ORF mutations, rebutting each inactivation mutation and suggesting that DRD₅ is, in fact, coding in these species (Supplemental Material 2). Regarding Myotis davidii (vesper bat), the fragmentation of the genomic region (Ns) flanked by upstream and downstream DRD₅ neighboring genes impeded us to infer the DRD₅ coding status in this species. Interestingly, regarding Ochotona princeps (American pika) mutational SRA validation, two scenarios were observed: approximately 50% of aligned reads supported the presence of a premature stop codon in the DRD₅ gene of this species, with the remaining set of aligning reads supporting the absence of a premature stop codon in the same species. This suggest that a polymorphic loss event might have occurred in this species, with non-functional and functional alleles currently segregating in the correspondent population.

Chrysochloris asiatica presents a non-functional DRD₅ gene

Next, we examined other mammalian species with annotated genome yet lacking DRD₅ gene annotations. Results were inconclusive regarding the coding status of Microcebus murinus (gray mouse lemur), Condylura cristata (star-nosed mole), J. jaculus (lesser Egyptian jerboa) and Erinaceus europaeus (western European hedgehog) DRD₅, due to the fragmentation (Ns) of the genomic region flanked by the upstream and downstream DRD₅ neighboring genes, and the unavailability of whole genome shotgun contigs spanning our target gene. For Elephantulus edwardii (Cape elephant shrew) we were able to deduce a fully functional DRD₅ CDS. Curiously, Chrysochloris asiatica (Cape golden mole) presented a single nucleotide insertion in the 5′ end of the gene (validated by genomic SRA, see Supplemental Material 3), suggesting that DRD₅ might be pseudogenized in this species.