Our full updated maximum likelihood phylogeny of the OTR/VTR gene family, including all 88 analyzed species and 519 sequences, has been deposited in Figshare (see Supplementary Material section). A smaller phylogeny including 55 species (Table 1) is shown in Figures 1 and 2. The smaller phylogeny, as well as the alignments underlying both phylogenies, have also been deposited in Figshare. Our conclusions are supported by both phylogenies. We have used the gene nomenclature VTR1A (AVPR1A), VTR1B (AVPR1B), OTR (OXTR), VTR2A (AVPR2), VTR2B and VTR2C.

The updated phylogenies are consistent with our previously published analyses (12, 15), showing 6 ancestral subtypes VTR2A, VTR2B and VTR2C (Figure 1) as well as VTR1A, VTR1B and OTR (Figure 2) diverging early in jawed vertebrate evolution. As in previous analyses, the VTR2C genes form a paraphyletic group, with the VTR2C sequences from teleost fishes (Teleostei) as well as some shark (Selachimorpha) VTR2C sequences clustering at the base of the VTR2A, VTR2B and the main VTR2C branches (Figure 1). Nonetheless, conserved synteny as well as exon/intron structures and amino acid sequence features support the identity of this seemingly basal branch as VTR2C. All other OTR/VTR subtypes are well-supported and, in general, follow the known and accepted phylogeny of jawed vertebrate groups. As in previous analyses, the coelacanth (Latimeria chalumnae) OTR/VTR sequences in some cases diverge from the expected phylogenetic position, see the VTR1B and OTR branches in Figure 2. Other discrepancies in our current phylogenies include the following:

Aside from these minor discrepancies, both phylogenies are generally well-supported and well-resolved, and allowed us to draw more detailed conclusions about the repertoires of OTR/VTR subtype genes in different vertebrate lineages than before.

To extend our previous analyses of jawed vertebrate chromosomal regions for the OTR/VTR family, we have added information from newer high-coverage genome assemblies, particularly genomes that are distantly related to those that were available when we published our previous report describing the four-fold chromosomal similarity (12). At that time, we focused our analyses on chromosomal regions in the human, chicken, spotted gar and zebrafish genomes, and produced a phylogeny that also included OTR/VTR sequence information from other jawed vertebrates, namely mouse (Mus musculus), gray short-tailed opossum (Monodelphis domestica), Carolina anole lizard (Anolis carolinensis), Western clawed frog, coelacanth, several teleosts, and the holocephalan cartilaginous fish elephant shark.

Here, we have added OTR/VTR sequence information from a larger number of vertebrate species, as well as chromosomal regions from the polypteriform fish reedfish, representing an even earlier-diverging branch of the ray-finned fishes, the cartilaginous fish thorny skate as well as a more recent high-quality assembly of the Western clawed frog genome. We focused on species with key phylogenetic and evolutionary positions as well as high-coverage genome assemblies mapped to chromosomes or pseudochromosomes. In particular, we benefited from the recent advancements in the sequencing of cartilaginous fish genomes (32–35) as well as the advancement of the Vertebrate Genomes Project (36). This study includes two additional mammalian species, including platypus (Ornithorhynchus anatinus), nine additional avian species representing a wide segment of the extant avian diversity, nine additional reptilian species, four additional amphibian species, the gray bichir, reedfish, American paddlefish (18) and sterlet (17), which occupy key positions at the base of ray-finned fish evolution, ten additional teleost species, seven additional cartilaginous fish species representing all three major lineages, as well as the cyclostomes inshore hagfish (Eptatretus burgeri) (37), Pacific lamprey (Entosphenus tridentatus) (38) and sea lamprey (Petromyzon marinus). The full list of species, including genome assembly and database information has been deposited in Figshare (see Supplementary Material section), and a subset is shown in Table 1.

Updating the data with newly available assemblies also allowed us to assign chromosome locations to previously unassigned genes, as for the Western clawed frog genes and several of the three-spined stickleback (Gasterosteus aculeatus) genes, or to correct previously incomplete gene predictions. We also found that many of the NCBI gene models we identified for this study matched our manually curated predictions from the previous study (12), lending some validity to the manual curation process we had employed. All chromosome locations and gene/sequence prediction notes are available in Figshare (see Supplementary Material).

We constructed a full maximum likelihood phylogeny based on all 519 identified OTR/VTR sequences from 88 species (see Supplementary Material), as well as a smaller phylogeny shown in Figures 1 and 2 based on the 55 representative species in Table 1. Both phylogenies support the emergence of the VTR1A, VTR1B, OTR, VTR2A, VTR2B and VTR2C genes early in vertebrate evolution, prior to the radiation of jawed vertebrates. We could identify all OTR/VTR subtype genes in cartilaginous fishes, except for VTR2A. We do not interpret this fact as support for the view that VTR2A arose in a bony fish (Osteichthyes) ancestor, millions of years after the emergence of VTR2A and VTR2B. Assuming an even rate of sequence evolution, our phylogenies suggest that the VTR2A, VTR2B and VTR2C genes arose concomitantly, and thus that the loss of VTR2A in cartilaginous fishes is a subsequent and independent event. This view is also supported by our chromosomal analyses of conserved synteny, discussed below.

Our suggested evolutionary scenario for the emergence and evolution of the OTR/VTR genes is shown in Figure 8, based on our updated phylogeny, species distribution and conserved synteny analyses. An ancestral OTR/VTR gene underwent a local gene duplication, giving rise to a VTR1/OTR ancestral gene and a VTR2 ancestral gene located on the same chromosome. This first local duplication took place a in chordate ancestor after the divergence of tunicates from the lineage leading to vertebrates, at the earliest, or in a vertebrate ancestor prior to the first of two WGD events, at the latest. Subsequently, these two ancestral genes duplicated in the 1R and 2R WGD events, as shown in Figure 9, after which one VTR1/OTR duplicate and one VTR2 duplicate were lost, leaving the ancestral jawed vertebrate subtype genes VTR1A, VTR1B and OTR on one side and VTR2A, VTR2B and VTR2C on the other. Due to the original local duplication, VTR1A and VTR2C as well as OTR and VTR2B remain located on the same chromosomes in many jawed vertebrate genomes.

Subsequently, there were several differential gene losses in several jawed vertebrate lineages. VTR2A was possibly lost from cartilaginous fishes, followed by the putative loss of VTR1B and VTR2B in sharks, rays and skates (Elasmobranchii) producing the smallest repertoire of OTR/VTR genes in any of the studied lineages. This is possibly related to the urea storage-based osmoregulatory system in cartilaginous fishes, which is unique among jawed vertebrates and whose endocrine control remains largely unknown (41). VTR2A also seems to have been lost from the lineage leading to Acipenseriformes, including the American paddlefish and sterlet, as well as from a common ancestor of birds and crocodilians. VTR1B was likely lost from the lineage leading to neopterygian fishes, including the holosteans spotted gar and bowfin as well as the teleosts. VTR2B was lost from the lineage leading to tetrapods after the divergence of the coelacanth lineage, while VTR2C was lost at least three times Sarcopterygii - from the coelacanth, mammals, and the lineage leading to caenophidian snakes.

The OTR/VTR genes duplicated further in the early teleost 3R WGD event, except for VTR1B which had been lost from neopterygian fishes (Figure 9). However, none of the analyzed teleost species preserve duplicate VTR2C genes. Subsequently VTR2C seems to have been lost independently several times from different teleost lineages. Likewise, several of the analyzed teleost species have lost either one or the other of the two VTR2B duplicates. We could also detect additional duplicates of VTR1Aa, OTRa, OTRb, VTR2Aa, VTR2Ba and VTR2Bb in the brown trout (Salmo trutta) genome (see Supplementary Material) that likely arose in the salmonid-specific 4R WGD event (42). The OTR/VTR genes also duplicated in the independent 3R WGD events in the lineages leading to the American paddlefish and the sterlet, except for VTR2A which likely had been lost in a common ancestor. The American paddlefish 3R has been dated at 50 MYA (18) and the sterlet 3R has been estimated at both 180 MYA (17) and 21.3 MYA (19).

In a previously published study (12), we identified 34 neighboring gene families that showed patterns of conserved synteny in the chromosomal regions harboring OTR/VTR genes. For most of these neighboring gene families, 23 out of 34, phylogenetic analyses using two different methods (neighbor-joining and maximum likelihood) supported gene duplications in the time window of the 1R and 2R WGD events before the jawed vertebrate radiation. For the remaining 11 neighboring gene families identified in our previous study (12), a subset supported duplications in the time window of 1R/2R in one of the two phylogenetic methods, and yet another subset is compatible with duplications in 1R/2R but lacked an appropriate outgroup to relatively date the time window of duplications. Thus, we concluded that the chromosomal regions harboring these genes likely arose through 1R and 2R as well. This has been the gold standard for detecting the patterns of four-fold paralogy that resulted from 1R and 2R for nearly two decades, and we have been able to elucidate the evolution of many jawed vertebrate gene families through these ancient WGDs in this way. Notably, not all gene families that expanded prior to the jawed vertebrate radiation show signs of having duplicated through 1R and 2R. Out of the four jawed vertebrate genes that encode the growth hormone and prolactin family of hormones - GH, SL, PRL1 and PRL2 - only PRL1 and PRL2 arose through either of the 1R or 2R WGDs (43).

In the present study, we re-analyzed seven representative neighboring gene families (Table 1) out of the original 34 in the human, chicken, spotted gar and zebrafish genomes and added novel chromosomal data from the reedfish, thorny skate and Western clawed frog genomes. These seven gene families were chosen based on three criteria: 1) their clear and well-supported phylogenies in our previous study; 2) the presence of member genes on all four paralogous chromosome regions that were identified; and 3) their representation across the length of these paralogous regions. With these seven neighboring gene families, we could detect the same pattern of four-fold conserved synteny around the OTR/VTR genes in the Western clawed frog, reedfish and thorny skate genomes (Figures 5, 6). This provides us with several key relative time points for the four-fold paralogy: at the divergence between cartilaginous fishes and bony fishes, relatively soon after the WGDs, at the base of ray−finned fishes, close to the divergence between ray-finned and lobe-finned fishes, and at the base of the tetrapod radiation. We could not identify the VTR2A-bearing chromosome region in the thorny skate genome (Figure 5). A possible explanation for this is that the thorny skate genome assembly provided by the Vertebrate Genomes Project was found to have a very high repeat content (54 %), which resulted in a large number of assembly gaps, even after manual modifications to the assembly pipeline and curation of the assembly to increase its quality (36). It is possible that the VTR2A-bearing chromosome region, and the VTR2A gene, are missing from the thorny skate genome assembly for this reason. The analyses of the white shark (Carcharodon carcharias), whale shark (Rhincodon typus), brownbanded bambooshark (Chiloscyllium punctatum) and cloudy catshark (Scyliorhinus torazame) genomes show a similarly high repeat content of approximately 50 % of the genome assemblies (33, 34). The repeat content of the elephant shark genome is also relatively high, at approximately 42 % [see Supplementary Note 4 in (33)]. We could not identify the VTR2A gene from any of the available cartilaginous fish genomes, and it is possible that this reflects the difficulty of integrating its surrounding chromosomal region into high-quality genome assemblies due to the presence of repeat sequences.

There have been several attempts at reconstructing the ancestral genome prior to the early vertebrate WGDs and to map the twice-duplicated ancestral chromosomes onto the genomes of extant vertebrate species. In the reconstruction by Nakatani et al. (9), the paralogous chromosome regions we have detected closely match the D0-D3 paralogous regions (marked in Figure 9) that arose through the quadruplication of ancestral linkage group D [see Figure 4 in (9)]. In the reconstruction by Putnam et al. (10), our paralogous chromosome regions match the four paralogous chromosomes that arose from the quadruplication of the ancestral chordate linkage group 13 [see Figure 3 in (10)]. The more recent reconstructions by Sacerdot et al. (44) and Simakov et al. (45) are also compatible with our findings. In the reconstruction by Sacerdot et al. (44), the paralogous regions we have detected closely match the reconstructed pre-1R chromosome 4 [see Figures 5 and 6 in (44)], and in the reconstruction by Simakov et al. (45) our paralogous regions in the human, chicken, Western clawed frog and spotted gar closely match the reconstructed chordate linkage group E [see Figure 2 in (45)]. The same pattern of four-fold paralogy has thus been detected by four independent studies, using different reconstruction methods, in addition to ours. Furthermore, the paralogous chromosome regions in teleost genomes, although heavily rearranged, match chromosome regions that have been shown to have originated in the teleost 3R WGD event in several independent studies (31, 46, 47).

An important corollary of the whole-genome duplication scenario is that the gene duplicates initially most likely had not only identical coding sequences, but also identical regulatory regions. Thus, the OTR, VTR1A and VTR1B genes started as “identical triplets”, and likewise for the three VTR2 genes. Subsequently, the OTR, VTR1A and VTR1B genes diverged functionally from each other resulting in three receptors with unique roles for each one, either by sub-functionalization (dividing the mother gene’s functions between the three daughter genes) or neo-functionalization (acquisition of novel functions). This would explain why they have been retained in the genomes of many vertebrate lineages (Figure 8). By the same token, the conservation of the three VTR2 genes in diverse vertebrate lineages (Figure 8) suggests that they too have evolved unique functional roles. However, little has been done to investigate the functions of V2B and V2C receptors. Yamaguchi et al. (48) characterized some functional aspects of elephant shark and teleost V2B and V2C receptors and discovered that they elicit an intracellular Ca²⁺ cascade on activation, like V1A, V1B and OTR and unlike V2A, which engage the adenylyl cyclase/cAMP signal transduction pathway. It is also known that the chicken V2 receptor, which is a V2C, primarily signals via Ca2+ as a second messenger (49). Interestingly, the VTR2B gene appears to have been lost in a common ancestor of tetrapods, suggesting that its role became dispensable in this lineage, possibly in relation to the transition from water to land. It will be interesting to identify its unique role among other vertebrates. Furthermore, there seems to have been an asymmetric loss of V2C from a mammalian ancestor and V2A from an archosaurian ancestor of birds and crocodilians (Figure 8). This might be related to the fact that the kidneys in birds and crocodilians do not contribute significantly to water reabsorption. They rely instead on the intestine and well-developed salt glands for water/electrolyte balance. In domestic chickens the V2C receptor is primarily expressed in the shell gland and the brain (49). It will likewise be interesting to identify the unique role of VTR2C among other vertebrates, especially those that have both VTR2A and VTR2C genes, like turtles (Testudines).

Similarly, it will be interesting to find out how VTR1B could be lost in a neopterygian ancestor of bowfin, gars and teleosts but retained in earlier-diverging ray-finned fish lineages, like bichirs and reedfish as well as sturgeon.

Our analyses of the expanded IL3 regions of teleost V2A receptors also have functional implications. We found that these IL3 regions contain several relatively well-conserved Tyr, Ser and Thr residues, as well as a high degree of intrinsic disorder. These sequence features indicate a degree of functionality as well as similarities to sequence features seen in IL3 regions from other GPCRs. It has been shown that GPCRs can be phosphorylated on intracellular Tyr, Ser and Thr residues by various receptor kinases, mediating the modulation of intracellular signaling cascades as well as GPCR attenuation/desensitization, endocytosis and intracellular trafficking (50). The low functional activity of the relatively short IL3 of the human V2A receptor has been demonstrated through mutagenesis studies (51, 52), however for many GPCRs IL3 contains important sites for downstream interactions with G proteins, receptor kinases and arrestins, for example. It has also been shown that the IL3 regions of GPCRs in particular harbor intrinsically disordered stretches that facilitate such interactions (53–55), as well as GPCR oligomerization (56).

Recent studies of the evolution of OTR/VTR genes that included cyclostome sequences have questioned the double WGD evolutionary scenario and have proposed other schemes. One study concluded that the analyses supported one WGD event that was followed by “chromosomal or segmental duplications and translocations” in the jawed vertebrate lineage (13). This conclusion was repeated more recently in a book chapter from the same group (14). However, their analyses did not include neighboring gene families. Besides, it is well known that lamprey sequences tend to cluster in ambiguous positions in vertebrate phylogenies, partly due to a unique pattern of sequence evolution at the protein level as a result of GC-rich codons, but also possibly due to asymmetric gene retention after the WGD events between cyclostomes and jawed vertebrates (57, 58). It is also unclear when jawless vertebrates diverged from the branch leading to jawed vertebrates with respect to the two rounds of WGD (59). This makes early vertebrate branching events more difficult to interpret.

The most recent report on the evolution of OTR/VTR genes favored one WGD event followed by large segmental duplications in the jawed vertebrate lineage giving rise to VTR1A and VTR1B as well as VTR2B and VTR2C (4). This study included both lamprey and hagfish sequences. The analyses did include neighboring gene families, but comparisons were only performed between lamprey and various jawed vertebrates, not among the jawed vertebrates. Thus, the conclusion was most likely affected by the distant relationship to the cyclostomes, which makes conserved synteny difficult to interpret. We were able to identify the same sea lamprey and inshore hagfish sequences as the previously published study by Theofanopoulou et al. (4), and included them in our full phylogeny. However, based on the known set of problems when trying to infer the orthology between cyclostome and jawed vertebrate sequences, we don’t consider these results conclusive, and much caution should be taken before trying to reconcile the ambiguous positions of the cyclostome sequences within phylogenies with any scenario of OTR/VTR gene evolution. Furthermore, recent reports from ancestral genome reconstructions have suggested that cyclostomes share only the first WGD event (1R) with jawed vertebrates (27, 45) whereupon the cyclostome lineage underwent a genome triplication (hexaploidization) (27), further accentuating the complications related to comparisons between cyclostomes and jawed vertebrates. The study by Nakatani et al. from 2021 (27), nevertheless firmly concluded that the lineage leading to jawed vertebrates has undergone two WGD events. Aside from asymmetric gene duplicate retention after 1R, it is possible that further OTR/VTR gene duplicates could have arisen through different mechanisms in the two lineages: the second WGD in jawed vertebrates and hexaploidization in cyclostomes. We are currently carrying out further careful study of the chromosome regions in the most recent sea lamprey and inshore hagfish genome assemblies with the aim to resolve the evolutionary history of the cyclostome OTR/VTR genes.

Another reason why Theofanopoulou et al. (4) proposed that jawed vertebrates may have undergone only one WGD followed by several segmental duplications was probably that they could identify five, rather than four, chromosome segments showing conserved synteny [see Figure 2 in (4)]. However, their two blocks of genes on chromosome 3 are actually one and the same. Only three gene families are present in both blocks, and all three were analyzed incorrectly: 1) The ARF and ARL8 genes do not form a quartet that arose in vertebrates, in fact ARF and ARL8 existed much earlier in evolution; 2) The ATP2B genes form a quartet on the four human chromosomes, 1, 3, 12 and X (Figure 5), however the authors have omitted the ATP2B1 gene on chromosome 12, misplaced the ATP2B3 gene on chromosome 3 when it is in fact located on chromosome 12, and included ATP6AP1, which is not part of this gene family; and 3) the CNTN3, -4, and -6 genes arose through a local gene triplication after 2R and the large distance to the CNTN3 gene is due to a recent rearrangement in the primate lineage. Thus, all analyzed gene families in the study by Theofanopoulou et al. (4) nicely support chromosomal quadruplication, hence two WGD events.

Yet another reason that Theofanopoulou et al. (4) proposed segmental duplications in jawed vertebrates is the lack of one of the four paralogous chromosome regions (bearing VTR2A) in the elephant shark genome, indicating that it might have arisen in a bony fish (Osteichthyes) ancestor. Our extended analysis of seven additional species confirms that this chromosome region indeed seems to be missing from cartilaginous fishes. However, as we discuss above, we cannot reject the possibility that these regions are missing from cartilaginous fish genome assemblies due to their high repeat content. We could not identify the VTR2A-bearing chromosome region in avian genomes either. Large chromosomal regions can indeed be lost, as has been shown for HOX clusters in teleosts after the 3R event.

Our sequence-based phylogenies as well as conserved synteny analysis strongly support an origin of all four paralogous chromosome regions before the jawed vertebrate radiation. The distinct four-fold paralogy is preserved in diverse jawed vertebrate lineages and consists of gene families that diverged early in vertebrate evolution. This provides a clear case for gene duplication through the 1R and 2R WGD events.

We previously suggested a simplified nomenclature for the family of six OT/VT receptors in vertebrates, V1A, V1B, OTR, V2A, V2B and V2C, based on both phylogenetic analysis and chromosomal location (11, 12). This nomenclature has so far been used in a number of publications (13, 48, 60–64). Recently, Theofanopoulou et al. (4), suggested naming the OTR/VTR genes VTR1A, VTR1B, OTR, VTR2A, VTR2B and VTR2C, foregoing the letter A for arginine-vasopressin currently used in the accepted gene names. We support this proposal and have used this nomenclature throughout the present study. However, Theofanopoulou et al. (4) have interchanged the VTR2A and VTR2C designations based on an inferred chronological order of emergence of the genes, suggesting that VTR2A (which they call VTR2C) and its chromosomal region arose later, in a bony fish ancestor. They base this order on the species distribution of VTR2A genes and their attempt to assign the lamprey V2R-type genes to specific subtypes. However, as we have shown, the VTR2A, VTR2B and VTR2C genes arose concomitantly in the 1R/2R WGD events, not in sequence, and the lamprey V2R-type genes have not yet been confidently assigned to specific subtypes. For these reasons we strongly suggest retaining the original ABC designations as they a) have already been used consistently in the literature for nearly a decade; and b) the genes arose concomitantly in the 1R/2R WGD events, not in a chronological sequence, as Theofanopoulou et al. (4) propose.

The OTR/VTR gene models/sequences identified in our previous study (12) were updated, where applicable, to the corresponding data from the latest available genome assemblies in the National Center for Biotechnology Information (NCBI) Genome database (http://www.ncbi.nlm.nih.gov/genome/) (65). The corresponding chromosomal location data were also updated. OTR/VTR gene models/sequences from additional vertebrate species were identified in the NCBI database using the same BLAST-based method as in our previous study (12), and added to the dataset. All identified OTR/VTR gene models/sequences, with their corresponding location data and database identifiers are available in Figshare (see Supplementary Material). Species, genome assembly and database information are shown in Table 1.

The identified OTR/VTR amino acid sequences were aligned using the MUSCLE algorithm (66) in AliView version 1.6 (67). The alignment was edited manually to adjust poorly aligned sequence stretches with respect to conserved motifs, exon boundaries and incomplete sequences. Missing sequence stretches caused by gaps in the genome assemblies were replaced with X-positions. Partial sequences are marked with asterisks (*) in the final alignment. The curated alignment was used to create a maximum-likelihood phylogeny using IQ-TREE (68) version 1.6.12 (August 15, 2019) supported by IQ-TREE’s non-parametric UltraFast Bootstrap (UFBoot) method (100 iterations) (69). The best-fit amino acid substitution model and substitution parameters were estimated using IQ-TREE’s ModelFinder with the -m TEST option (70), and the proportion of invariant sites was optimized by using the –opt−gamma−inv option. The phylogeny was rooted with the common octopus (Octopus vulgaris) cephalotocin and octopressin receptor sequences (71).

The following neighboring gene families from (12), are included: ATP2B (ATPase plasma membrane Ca²⁺ transporting protein B), CACNA1-L (Voltage-gated Ca²⁺ channel alpha-1 subunit, L-type), CAMK1 (Calcium/calmodulin dependent protein kinase 1 family), L1CAM (L1-related cell adhesion molecules), PLXNA (Plexin A family), SRGAP (SLIT-ROBO Rho GTPase activating protein), and SYP (Synaptophysin/synaptoporin). Briefly, these and other neighboring gene families were identified by comparing the gene composition of the chromosomal regions bearing visual opsin, OTR/VTR, CACNA1-L, GNAT and GNAI genes in the human genome. Those gene families with member genes on at least two chromosomal regions were selected for further analysis [see Methods in (12)].

Gene models and locations for all neighboring genes were updated to the corresponding data from the latest available genome assemblies in the NCBI Genome database. In addition, neighboring family gene model and location data from the Western clawed frog, reedfish and thorny skate were added in order to identify the paralogous OTR/VTR-bearing chromosome regions in these key species.

For the conserved synteny analysis of 3R-generated paralogous regions in the Asian arowana (Sleropages formosus), the conserved neighboring gene pairs were identified through Genomicus v100.01 (72) and the time window of gene pair emergence was inferred from the corresponding Ensembl gene trees (73). Where the orthology relationships to the zebrafish genes was clear in the Ensembl gene trees, gene names with “a” and “b”-designations have been suggested. Where there was only one zebrafish ortholog, this was designated the “a”-duplicate.