# The sensory shark: high-quality morphological, genomic and transcriptomic data for the small-spotted catshark *Scyliorhinus canicula* reveal the molecular bases of sensory organ evolution in jawed vertebrates

H. Mayeur, J. Leyhr, J. Mulley, Nicolas Leurs, L. Michel, K. Sharma, R.

Lagadec, J.-M. Aury, O.G. Osborne, P. Mulhair, et al.

► To cite this version:

H. Mayeur, J. Leyhr, J. Mulley, Nicolas Leurs, L. Michel, et al. The sensory shark: high-quality morphological, genomic and transcriptomic data for the small-spotted catshark *Scyliorhinus canicula* reveal the molecular bases of sensory organ evolution in jawed vertebrates. Molecular Biology and Evolution, 2024, 41 (12), pp.msae246. 10.1101/2024.05.23.595469. http://hal-04785835

HAL Id: hal-04785835

https://hal.science/hal-04785835v1

Submitted on 20 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Sensory Shark: High-quality Morphological, Genomic and Transcriptomic Data for the Small-spotted Catshark

## *Scyliorhinus Canicula* Reveal the Molecular Bases of Sensory Organ Evolution in Jawed Vertebrates

Hélène Mayeur ID,<fnref n="1" /> Jake Leyhr ID,<fnref n="2" /> John Mulley ID,<fnref n="3" /> Nicolas Leurs ID,<fnref n="4" /> Léo Michel ID,<fnref n="1" />,†

Kanika Sharma,<fnref n="5" /> Ronan Lagadec ID,<fnref n="1" /> Jean-Marc Aury ID,<fnref n="6" /> Owen G. Osborne ID,<fnref n="3" /> Peter Mulhair ID,<fnref n="7" />

Julie Poulain ID,<fnref n="6" /> Sophie Mangenot ID,<fnref n="6" /> Daniel Mead ID,<fnref n="8" /> Michelle Smith,<fnref n="8" /> Craig Corton,<fnref n="8" />

Karen Oliver,<fnref n="8" /> Jason Skelton,<fnref n="8" /> Emma Betteridge,<fnref n="8" /> Jale Dolucan ID,<fnref n="9" />,<fnref n="10" /> Olga Dudchenko ID,<fnref n="9" />,<fnref n="10" />

Arina D. Omer ID,<fnref n="9" />,<fnref n="10" /> David Weisz ID,<fnref n="9" />,<fnref n="10" /> Erez L. Aiden ID,<fnref n="9" />,<fnref n="10" />,<fnref n="11" /> Shane A. McCarthy,<fnref n="8" /> Ying Sims,<fnref n="8" />

James Torrance ID,<fnref n="8" /> Alan Tracey,<fnref n="8" /> Kerstin Howe ID,<fnref n="8" /> Tobias Baril ID,<fnref n="12" />,† Alexander Hayward ID,<fnref n="12" />

Camille Martinand-Mari ID,<fnref n="4" /> Sophie Sanchez ID,<fnref n="2" />,<fnref n="13" /> Tatjana Haitina ID,<fnref n="2" /> Kyle Martin ID,<fnref n="14" />

Sigrun I. Korsching ID,<fnref n="5" />,* Sylvie Mazan ID,<fnref n="1" />,* Mélanie Debiais-Thibaud ID,<fnref n="4" />,*

## Abstract

Cartilaginous fishes (chondrichthyans: chimeras and elasmobranchs -sharks, skates, and rays) hold a key phylogenetic position to explore the origin and diversifications of jawed vertebrates. Here, we report and integrate reference genomic, transcriptomic, and morphological data in the small-spotted catshark *Scyliorhinus canicula* to shed light on the evolution of sensory organs. We first characterize general aspects of the catshark genome, confirming the high conservation of genome organization across cartilaginous fishes, and investigate population genomic signatures. Taking advantage of a dense sampling of transcriptomic data, we also identify gene signatures for all major organs, including chondrichthyan specializations, and evaluate expression diversifications between paralogs within major gene families involved in sensory functions. Finally, we combine these data with 3D synchrotron imaging and in situ gene expression analyses to explore chondrichthyan-specific traits and more general evolutionary trends of sensory systems. This

Received: July 07, 2024. Revised: November 21, 2024. Accepted: November 21, 2024

© The Author(s) 2024. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution.

This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site -for further information please contact journals permissions@oup.com.

Open Access

<footnote n="1">¹Sorbonne Université, CNRS, Biologie Intégrative des Organismes Marins, BIOM, Banyuls-sur-Mer, France</footnote>

<footnote n="2">²Department of Organismal Biology, Uppsala University, Uppsala, Sweden</footnote>

<footnote n="3">³School of Environmental and Natural Sciences, Bangor University, Bangor, Gwynedd LL57 2UW, UK</footnote>

<footnote n="4">⁴Institut des Sciences de l'Evolution de Montpellier, ISEM, University of Montpellier, CNRS, IRD, EPHE, Montpellier, France</footnote>

<footnote n="5">⁵Institute of Genetics, Faculty of Mathematics and Natural Sciences of the University at Cologne, Cologne 50674, Germany</footnote>

<footnote n="6">⁶Génomique Métabolique, Génoscope, Institut François Jacob, CEA, CNRS, Univ Evry, Université Paris-Saclay, Evry 91057, France</footnote>

<footnote n="7">⁷Department of Biology, University of Oxford, Oxford OX1 3SZ, UK</footnote>

<footnote n="8">⁸Sequencing Department, Wellcome Sanger Institute, Cambridge CB10 1SA, UK</footnote>

<footnote n="9">⁹The Center for Genome Architecture and Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, TX 77030, USA</footnote>

<footnote n="10">¹⁰The Center for Theoretical Biological Physics, Rice University, Houston, TX 77030, USA</footnote>

<footnote n="11">¹¹Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA</footnote>

<footnote n="12">¹²Centre for Ecology and Conservation, University of Exeter, Cornwall TR10 9FE, UK</footnote>

<footnote n="13">¹³European Synchrotron Radiation Facility, Grenoble, France</footnote>

<footnote n="14">¹⁴Department of Earth Sciences, Natural History Museum, London SW7 5BD, UK</footnote>

<footnote>†Present addresses: The Department of Fundamental Neuroscience, The University of Lausanne, Lausanne 1005, Switzerland; Laboratory of Evolutionary Genetics, Institute of Biology, University of Neuchâtel, Neuchâtel 2000, Switzerland.</footnote>

<footnote>*Corresponding authors: E-mails: sigrun.korsching@uni-koeln.de; mazan@obs-banyuls.fr; melanie.debiais-thibaud@montpellier.fr.</footnote>

<footnote>Associate editor: Dr. Jeffrey Townsend</footnote>
approach brings to light, among others, novel markers of the ampullae of Lorenzini electrosensory cells, a duplication hotspot for crystallin genes conserved in jawed vertebrates, and a new metazoan clade of the transient-receptor potential (TRP) family. These resources and results, obtained in an experimentally tractable chondrichthyan model, open new avenues to integrate multiomics analyses for the study of elasmobranchs and jawed vertebrates.

**Key words:** vertebrate evolution, cartilaginous fish genome, olfactory receptors, TRP ion channels, transcriptomics, electroreceptors.

Introduction

Cartilaginous fishes, or chondrichthyans, comprise elasmobranchs (i.e. sharks, skates, and rays) and holocephalans (chimeras) (Janvier 1996), and are a group of major interest for vertebrate evolutionary biologists. As sister group to bony fishes, which contain all major vertebrate model organisms (e.g. mouse, chicken, Xenopus, zebrafish), they occupy an important phylogenetic position for elucidating early evolutionary events at the jawed vertebrate crown node (Kuraku 2021; Tan et al. 2021; Yamaguchi et al. 2021; Zhang et al. 2022; Marlétaz et al. 2023). This position makes them important models for determining genomic bases of phenotypic evolution, especially the origin of ancestral features of jawed vertebrates and subsequent taxon-specific diversification. Such approaches require the availability of high-quality genomic data from across cartilaginous and bony fish diversity, coupled with a comprehensive knowledge of the morphological evolution of traits of interest. In this respect, more than 1,000 species of chondrichthyans remain massively under-represented. However, past few years have seen a growing effort to produce genomic data for this group resulting in a well-assembled genome for a holocephalan, the elephant shark *Callorhinchus milii*, and several chromosome-level assembled genomes for elasmobranch fishes including batooids (smalltooth sawfish *Pristis pectinate*, thorny skate *Amblyraja radiata*, little skate *Leucoraja erinacea*) and seilachians (one Carcharhiniform, the small-spotted catshark *Scyliorhinus canicula* (this study); one lamniform, the great white shark *Carcharodon carcharias*; two orectolobiforms, the epaulette shark *Hemiscyllium ocellatum*, and the white spotted bamboo shark *Chiloscyllium plagiosum*; one

squaliform, the spiny dogfish *Squalus acanthias*), as listed in Table 1. Additional genomes have been sequenced with a lower quality of assembly, and their annotations are also available (Weber et al. 2020; Stanhope et al. 2023; Zhou et al. 2023) notably through the Squalomix project (Hara et al. 2018; Nishimura et al. 2022). Even fewer cartilaginous fishes have been the focus of integrative studies associating genomic or transcriptomic data with different organizational levels from organism to organ, to cell, or successive developmental stages (see Gillis et al. 2022). Reasons for this lack of integrative data include the difficulty of accessing biological material (e.g. protection status for certain species, such as the smalltooth sawfish), biological features (e.g. size and reproduction, such as for the great white shark and whale shark), and more generally the lack of previous physiological/anatomical studies.

This work aims at filling this gap, focusing on sensory systems and using the small-spotted catshark *Scyliorhinus canicula*, an emerging model organism, as the reference species. Understanding how sensory systems diversify during evolution is crucial to apprehend the multiple strategies employed by animals to process environmental information and it typically requires a multi-level approach, from gene to cell and organ. The rise of paired sensory organs, including our eyes, nose, and ears, was a key innovation of early vertebrates, emerging alongside a complex, centralized brain and the prototypical "new head" of vertebrates (Gans and Northcutt 1983). The lateral line system is composed of two components: (i) the mechanosensory nemotasts located just below the epidermis or in a canal running into the dermis, and (ii) the electrosensory ampullary organs, also formed of canals, but filled with a specialized

<table><caption>Table 1 Summary of chromosome-level assembly features available in cartilaginous fishes with focus on the small-spotted catshark *Scyliorhinus canicula* (*), and in the partially assembled genome of the elephant shark *Callorhinchus milii*</caption><thead><tr><th>Order</th><th>Assembly (NCBI)</th><th>Species name</th><th>N</th><th>Genome size</th><th>N > 40 Mb (%)</th><th>Maximum chrom. size</th></tr></thead><tbody><tr><td rowspan="3">Orectolobiforms</td><td>sRhiTyp1.1</td><td>Rhincodon typus</td><td>[1]</td><td>50</td><td>2.9 Gb</td><td>28 (56%)</td><td>178 Mb</td></tr><tr><td>ASM401019v1</td><td>Chiloscyllium plagiosum</td><td>[2]</td><td>51</td><td>3.6 Gb</td><td>35 (68%)</td><td>156 Mb</td></tr><tr><td>sHemOce1.pat.X.cur</td><td>Hemiscyllium ocellatum</td><td>[3]</td><td>52 + XY</td><td>4.0 Gb</td><td>39 (73%)</td><td>167 Mb</td></tr><tr><td>Squaliforms</td><td>ASM3039002v1</td><td>Squalus acanthias</td><td>[4]</td><td>29</td><td>3.7 Gb</td><td>25 (86%)</td><td>266 Mb</td></tr><tr><td>Lamniforms</td><td>sCarCar2.pri</td><td>Carcharodon carcharias</td><td>[5]</td><td>40 + XY</td><td>4.3 Gb</td><td>25 (60%)</td><td>274 Mb</td></tr><tr><td>Carcharhiniforms*</td><td>sScyCan1.1*</td><td>Scyliorhinus canicula</td><td>[6]</td><td>31</td><td>4.2 Gb</td><td>21 (68%)</td><td>313 Mb</td></tr><tr><td rowspan="2">Rajiforms</td><td>Leri_hhj_1</td><td>Leucoraja erinacea</td><td>[7]</td><td>50</td><td>2.2 Gb</td><td>19 (38%)</td><td>170 Mb</td></tr><tr><td>sAmbRad1.1.pri</td><td>Amblyraja radiata</td><td>...</td><td>49</td><td>2.6 Gb</td><td>24 (48%)</td><td>193 Mb</td></tr><tr><td>Rhinopristiforms</td><td>sPriPec2.1.pri</td><td>Pristis pectinata</td><td>...</td><td>46 + X</td><td>2.3 Gb</td><td>20 (42%)</td><td>150 Mb</td></tr><tr><td>Holocephalan</td><td>IMCB_Cmil_1.0</td><td>Callorhinchus milii</td><td>[8]</td><td>na</td><td>0.9 Mb</td><td>na</td><td>na</td></tr></tbody></table>

N: number of chromosomes in the haploid genome. Maximum chromosome size as described in GenBank. Reference publications [1] (Read et al. 2017), [2] (Zhang et al. 2020), [3] (Sendell-Price et al. 2023), [4] (Wagner et al. 2023), [5] (Marra et al. 2019), [6] this work, [7] (Yoo et al. 2022; Marlétaz et al. 2023), [8] (Venkatesh et al. 2007; Venkatesh et al. 2014).
secreted hydrogel which has high conductivity and, as a result, weak electric fields can be detected by sensory cells in the deepest part of the tube (the ampullae). Both components of the lateral line were lost in amniotes during the transition to life out of the water, while amphibians retain them in larval or aquatic adult stages. Cartilaginous fishes display the complete set of paired ectodermal sensory organs (eyes, olfactory organs, ears, electrosensory organs, and lateral line; Collin 2012). However, cartilaginous and bony fish groups diverged over about 420 My ago and cartilaginous fishes diversified into a myriad of forms with a great range of habitats and life history traits (Ebert et al. 2021). This independent evolution is expected to have generated significant variation in their sensory systems, as observed in anatomical aspects of the electrosensory system (reviewed by Collin 2012 and see Newton et al. 2019; Bottaro 2022) and olfactory system (Schluessel et al. 2008; Yopak et al. 2015), as well as in sensory cell types such as photoreceptors (Hart et al. 2020; Yamaguchi et al. 2021).

We have focused on the small-spotted catshark *Scyliorhinus canicula* to provide a reference, laying the groundwork for the analysis of this diversity both within cartilaginous fish and at the origin of jawed vertebrates (Coolen et al. 2008). With easy access across developmental stages, a small body size, non-endangered status, high abundance in the North-East Atlantic, and easy maintenance in aquaria, this species has become a target elasmobranch model for research in European laboratories for more than 20 yr (Mazan et al. 2000; Tanaka et al. 2002; Coolen et al. 2007; O'Neill et al. 2007; Sakamoto et al. 2009; Oulion et al. 2010; Debiais-Thibaud et al. 2011; Rodríguez-Moldes et al. 2011; Mulley et al. 2014; Tulenko et al. 2016; Redmond et al. 2018; Mayeur et al. 2021). However, until recently, a major limitation was the absence of a chromosome-level genome assembly.

Here, we describe the assembled and thoroughly annotated genome of the small-spotted catshark, a wide series of transcriptomic data from adult tissues and embryonic stages, and cutting-edge 3D morphological data, now allowing integrative studies in this species, from the level of individual sequences to the whole organism. We exemplify how these resources can be combined to better understand the molecular bases of evolution, focusing on paired sensory organs. From the RNA-seq data, we identify the most highly expressed genes in each sensory organ. We then consider candidate gene families known to play major roles in sensory organs: opsins and crystallins in the retina and lens respectively; the four main vertebrate olfactory receptor families in the nose; and the transient-receptor potential (TRP) family, essential for diverse sensory functions. For each family, we describe the catshark gene repertoire, examine its evolution in relation to other vertebrates, and analyze tissue-specific expression patterns. Sites of expression are interpreted in the light of new 3D anatomical data and histology, facilitating a deeper understanding of the development and organization of sensory organs in the small-spotted catshark.

These integrative data highlight ancestral jawed vertebrate features of the electrosensory system and evolutionary trends in the olfactory and visual systems of cartilaginous fishes.

## Results

### Genomic and Transcriptomic Resources

#### Genomic Resources

The chromosome-level genomic assembly was obtained from one male small-spotted catshark collected in the Mediterranean Sea (western Gulf of Lion, off Port-Vendres, France). The sequencing strategy relied on a combination of PacBio, 10× Genomics Chromium, Bionano and Hi-C, which resulted in a chromosome-level genome reconstruction including 31 chromosomes (covering 97.7% of the total genome) and 614 unassembled contigs, for a total of 4.22 Gb (Fig. 1a and b). These data are available in NCBI as sScyCan1.1 (RefSeq GCF_902713615.1), and summary statistics for the assembly are presented in supplementary table S1, Supplementary Material online. BlobToolKit (Challis et al. 2020) GC-coverage and cumulative sequence plots are provided in supplementary fig. S1, Supplementary Material online. Twenty-one of these chromosomes are longer than 40 Mb (ranging between 71 and 313 Mb) and should therefore be considered macrochromosomes (Marlétaz et al. 2023), with the remaining ten microchromosomes (which includes the X chromosome) ranging in size between 13 Mb and 30 Mb (supplementary table S2, Supplementary Material online).

The genome of the small-spotted catshark has therefore over 60% macrochromosomes, as observed in most other shark genomes and in contrast to batoids species (Table 1). Macrosyteny is generally well conserved across elasmobranch fishes (Fig. 1c). However, the 31 chromosomes identified here represent a smaller karyotype than that reported for several other chondrichthyan species (Fig. 1d). Large synteny blocks carried by one chromosome in *S. canicula* are often distributed on 2 to 3 chromosomes in *C. carcharias* and *P. pectinata*, suggesting a number of recent chromosome fusions in the catshark or parallel fusions in the 2 other species (e.g. catshark chromosomes 1 or 4, Fig. 1c).

Analysis of GC content (supplementary figs. S1 and S2, Supplementary Material online) suggests a relatively stable GC percentage across chondrichthyan chromosomes, with an AT-bias in short (~1Mb) unplaced scaffolds in *S. canicula*. Repeat sequences are highly abundant in the Mediterranean small-spotted catshark genome, with transposable element sequences accounting for more than two-thirds of total genomic content (68.14%, supplementary table S3, Supplementary Material online, Fig. 1e). LINEs are the most abundant transposable element sequence type by far (34.2%, supplementary table S3, Supplementary Material online) and account for the vast majority of recently and currently active elements, as indicated by low sequence distances of annotated elements to their respective
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Genome assembly of *Scyliorhinus canicula*. a) Hi-C contact map of the assembly, visualized using HiGlass. The absence of off-diagonal signals illustrates the correct assignment of sequence to chromosome-level scaffolds. The remaining unassigned contigs are shown as bars at the bottom and right of the map. b) The BlobToolKit Snailplot shows N50 metrics and BUSCO gene completeness. The main plot is divided into 1,000 size-ordered bins around the circumference with each bin representing 0.1% of the 4.2 Gb assembly. The distribution of scaffold lengths is shown in dark gray with the plot radius scaled to the longest scaffold present in the assembly (314 Mb). Orange and pale-orange arcs show the N50 and N90 scaffold lengths (199 Mb and 98 Mb), respectively. The pale gray spiral shows the cumulative scaffold count on a log scale with white scale lines showing successive orders of magnitude. The blue and pale-blue area around the outside of the plot shows the distribution of GC, AT and N percentages in the same bins as the inner plot. A summary of complete, fragmented, duplicated, and missing BUSCO genes in the vertebrata_odb10 set is shown in the top right. c) Macrosynteny shared between three elasmobranch species. Synteny blocks based on alignment using all coding genes in three species, with 15,540 placed in shared synteny blocks between *Scyliorhinus canicula* and another shark, the great white shark *Carcharodon carcharias*, and 14,832 between *Scyliorhinus canicula* and a batoid, the small-tooth sawfish *Pristis pectinata*. Colors represent conserved synteny blocks and are ordered based on the focal species *Scyliorhinus canicula*. d) Genome evolution across chondrichthyans: species phylogeny (left) consisting of 16 species including 3 outgroup species (see Material and methods). Bar plots show genome size distribution per species, and chromosome count per species (both estimated from genomes assemblies). e) Pie-chart: proportion of the main TE classes in the assembled genome; graph: repeat landscape plot showing the proportion of repeats at different genetic distances (%) to their respective RepeatModeler consensus sequence, genetic distance is calculated under a Kimura 2 parameter model with correction for CpG site hypermutability: lower genetic distances suggest shorter time of divergence.</figcaption></figure>
consensus sequences (Fig. 1e). In addition, many long terminal repeat (LTR) transposable elements (11.7%, supplementary table S3, Supplementary Material online) and DNA transposable elements (6%, supplementary table S3, Supplementary Material online) sequences are also present in the small-spotted catshark genome (supplementary table S3, Supplementary Material online, Fig. 1e). However, both LTR TEs and DNA TEs appear to be decreasing in activity over more recent timescales (Fig. 1e).

A previously sequenced genome was obtained from an Atlantic individual (sex unknown) through Illumina sequencing of genomic DNA libraries (average fragment lengths ranging between 180 bp and 8 Kb) and resulted in a 4.31 Gb draft assembly, consisting of 3,502,619 scaffolds (>200 bp; maximum size 206 Kb; N50 = 6,665 bp). Both assemblies could be aligned over more than 90% of their length (Fig. 2b), with a 98.3% nucleotide identity. Unaligned sequences (about 9% and less than 1% of the Atlantic and Mediterranean genomes, respectively) essentially consist of repetitive sequences, and no coding sequences were identified in these data.

Using RNA-seq data acquired from individuals of known origin, we estimated the genomic diversity in a Mediterranean and two Atlantic (North Wales and English Channel) catshark populations (n = 9 individuals, n = 3 and n = 2, respectively) using an SNP-based population differentiation analysis. The mean FST value across all variants is 0.044 (supplementary table S2, Supplementary Material online) and the Mediterranean and Atlantic samples are better described with distinct ancestry (K = 2) with minimal admixture (Fig. 2a, c, and d) supporting differentiation between Mediterranean and Atlantic populations, in accordance with previous studies (Gubili et al. 2014; Kousteni et al. 2015; Ramírez-Amaro et al. 2018; Manuzzi et al. 2019; Melis et al. 2023). Significant outliers in the FST distribution were scattered across all chromosomes (Fig. 2e). Three genomic long-range regions (on chromosomes 6, 13, and 17) display high interpopulation and intra-population nucleotide diversity values (boxed zones in Fig. 2e to h). Gene density also was high in these regions while GC content is often low (Fig. 2i and j): scanning through the annotated loci identified gene clusters with apparently a large number of recent duplicates (e.g. histone coding genes on chromosome 13, zinc-finger binding transcription factor coding genes on chromosomes 13 and 17). Our use of an RNA-seq dataset obtained from gonads and early embryonic stages necessarily restrains the interpretation of the FST analysis. Several statistics of this genomic landscape show significant and strong correlation (supplementary fig. S2, Supplementary Material online) supporting the idea of high FST and nucleotide diversity values to come with regions of high gene density and GC content. As recombination rate is known to correlate with gene density and GC content, this genomic landscape may result from several forces, including balancing or divergent selection but also regional variation in recombination and mutation rates.

Transcriptomic Resources

Transcriptomic resources from a reference set of adult tissues and embryonic stages were generated by Illumina sequencing of cDNAs and used to annotate the Mediterranean chromosome-level genome, resulting in the prediction of 24,473 gene models. Alignment of these gene models against the Atlantic genome highlighted a 98.22% sequence identity (Fig. 2b), indicating substantial sequence variation in transcribed sequences between both origins. In order to improve the representativity of this gene model reference and extend the length of reference transcripts notably in 3'UTRs, we generated a second reference, referred to hereafter as *ncbi-utrs*, using a clustering approach and including NCBI collapsed isoforms plus transcriptomic data obtained by the CEL-Seq2 protocol (Hashimshony et al. 2016), which results in 3'UTR enriched sequence data. This *ncbi-utrs* database contains a total of 30,348 sequences (N50 = 5,181 bp, mean length = 2,980 bp) and harbors a 95.7% BUSCO score with 3,208 complete BUSCO groups out of 3,354 searched (69 fragmented and 77 missing). In accordance with previous studies (Oulion et al. 2010), we identified 34 Hox genes in this dataset (supplementary fig. S3, Supplementary Material online, supplementary table S4, Supplementary Material online). We found no evidence for the HoxC cluster known to be present in the elephant shark *Callorhinchus milii* (Ravi et al. 2009) and the zebra shark *Stegostoma tigrinum* (Yamaguchi et al. 2023b). Over 24,473 NCBI sequences, 11,945 (48.8%) were extended with a mean addition of about 1,845 bp per prolonged sequence. This *ncbi-utrs* reference was used for all subsequent analyses. Its mapping onto the Mediterranean and Atlantic genomes suggests excellent coverage for both assemblies, with only 1.88% and 1.32% of the "ncbi-utrs" sequences respectively missing in each one.

We next took advantage of the transcriptomic reference to gain insight into (1) the dynamics of gene expression during organogenesis and (2) expression specificities of a broad range of adult tissues. Concerning organogenesis, we focused on a temporal window encompassing stages 22, 24, 26, 30, and 31, characterized by major morphogenetic changes in the nervous, sensory, cardiovascular, respiratory, excretory, and musculoskeletal systems (Fig. 3a). This covers a long period of embryonic development: one month between stage 22 and stage 30, and another 20 d for stage 31 itself (Ballard et al. 1993). For each of these stages and each *ncbi-utrs* gene model, TPM (Transcripts Per kilobase Millions) and Z-score values were calculated to provide quantified indicators of relative expression levels and expression specificities of genes across this series of embryonic stages (supplementary dataset S1, Supplementary Material online). Pearson correlation analysis for all genes across stages 22 to 31 shows that stages 30 to 31 are most clearly correlated and set apart from earlier stages, suggesting major developmental transitions between stages 26 and 30 (Fig. 3b). To further identify genes showing similarities in their expression dynamic between stages 22 and 31, we conducted a
<figure><img src="image_3.png" /><figcaption>Fig. 2. Atlantic and Mediterranean populations of the small-spotted catshark. a) Principal Component Analysis based on 12,005 SNPs separating Atlantic and Mediterranean samples. b) Proportion of reciprocal coverage between the Atlantic and the Mediterranean small-spotted catshark genomes, and between genomes and the reference gene model used in this study (percentage of the sequences from the source (first column) that can be mapped on the target (first line). c) ADMIXTURE analysis at K = 2. Bars show estimated ancestry proportion of each individual, and colored lines at the top and bottom of the plot indicate the source population (same legend as in A). d) Cross-validation error (CV error; K = 2 to 6; b) indicates that K = 2 was the best fit to the data. e-j) Genome-wide patterns of differentiation (e; Weir and Cockerham's FST, based on 172,078 SNPs, dotted line indicates the mean value for all SNPs in the genome), absolute divergence (f; dxy, dotted line indicates the absolute value for the whole genome), nucleotide diversity in Atlantic (g; πAll) and Mediterranean (h; πMed) populations, gene density (i; annotated loci per Mb) and GC content (j). Mean values for non-overlapping 1 Mb windows are shown as points in alternating colors for each chromosome. Red points in the FST plot show site-level FST values for variants which were significant outliers in the PCAdapt analysis.</figcaption></figure>

hierarchical clustering approach. This analysis was restricted to a set of 10,651 genes, selected for the coherence of their expression changes between successive stages

(expression autocorrelation score across stages > 0.1: supplementary dataset S2, Supplementary Material online). In line with the Pearson correlation analysis, it led
<figure><img src="image_4.png" /><figcaption>(a) Developmental stages</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Pearson correlations</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) GO-term analysis</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) Hierarchical clustering</figcaption></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Fig. 3.** Transcriptomic comparison between stage 22, 24, 26, 30 and 31 embryos. a) Broad morphological characteristics of the stages analyzed (after Ballard et al. 1993). Lateral (stages 22 and 24) or ventral views (stages 26, 30 to 31) of a catshark embryo (stage 22) or its embryonic cephalic region (stages 24 to 31) are schematized for each stage. b) Pairwise Pearson correlations between stages across all genes with sum of TPM > 50. A black box shows the correlation observed between stages 30 and 31. c) Gene groupings obtained by hierarchical clustering of Z-score transformed expression of the 10,651 genes with regionalized expression (autocorrelation > 0.1). The two major clusters are respectively numbered 1 and 2, the latter is further subdivided into four groups, referred to as 2.1, 2.2, 2.3, and 2.4. The corresponding nodes are shown by arrowheads in the dendrogram on the left, and the clusters are delineated by vertical bars on the right. d) Selection of enriched GO-terms in clusters 1, 2.1, 2.2, 2.3, and 2.4, with P-values indicated. Detailed results for the GO-term analysis are shown in **supplementary dataset S3, Supplementary Material online**. C1 to C6 refer to pharyngeal clefts, from anterior to posterior. AP: anterior-posterior; dev: development; diff: differentiation; DV: dorsal-ventral; PD: proximal-distal; pns: peripheral nervous system.

to the identification of two major gene clusters, referred to as clusters 1 and 2 (Fig. 3c; **supplementary dataset S2, Supplementary Material online**). These two clusters respectively contain genes with higher expressions at stages 22 to 24 (cluster 1) and stages 30 to 31 (cluster 2) (Fig. 3c). Cluster 2 is further subdivided into four groups: clusters 2.1 and 2.2 contain genes whose expression increases at stage 26, with peaks at either stage 30 (cluster 2.1) or stage 31 (cluster 2.2), while both other clusters contain genes with a later increase of expression at stage 30 together with highest expression either at stage 30 (cluster 2.3) or at stage 31 (cluster 2.4). Gene ontology (GO)-term analysis reveals very different functional annotations between each of these 5 clusters (Fig. 3d; **supplementary dataset S3, Supplementary Material online**). Genes in cluster 1 (expression repressed at stages 30 to 31) are frequently associated with cell division (DNA replication, cell cycle,

chromosome segregation, and their regulation), pattern specification, and some aspects of neuronal differentiation (neural fate specification and spinal cord motor neurons). As an example, several Hox genes and related long non-coding RNAs belonged to cluster 1 (**supplementary fig. S3, Supplementary Material online**), with peaks of expression either early (for “anterior” Hox genes) or late (more “posterior” Hox genes) in the stages 22 to 26 period.

GO-terms related to different aspects of organogenesis and myelination prevail in clusters 2.1 and 2.2, respectively (genes showing an increase of expression by stage 26). Genes in cluster 2.3 (expression established at stage 30) show a large number of highly significant enrichments in GO-terms related to the generation and differentiation of neurons, the formation of axonal and dendritic projections, the assembly and function of synapses, and the homeostasis of neurotransmitters, suggesting that the
establishment of the basic architecture of the catshark neural circuitry is a major ongoing process at stage 30. GO-terms related to behavioral, cognitive, or sensory processes as well as to the formation of the circulatory system are also significantly enriched at this stage. A significant association with GO-terms related to synaptic and neurotransmitter signaling, as well as sensory perception is maintained for genes in cluster 2.4 (expression peaking at stage 31). Additional GO-terms enriched in this cluster include heart function, the differentiation of smooth and striated muscles, the differentiation of different lymphoid cell types, and hormonal and immune responses, suggesting that these processes become prevailing between stages 30 and 31 (Fig. 3d and supplementary dataset S3, Supplementary Material online). Taken together, these data highlight relationships between specific gene expression dynamics and the initiation of biological or developmental processes in the developmental window analyzed. The functional annotation thus provides a basis to search for underlying developmental regulators.

Concerning the analysis of gene specificities in adult tissues, we focused on organs of the nervous, digestive, excretory, cardiovascular, and reproductive systems, as well as on

a selection of sensory organs and skeletal structures (Fig. 4a). Z-Score values were calculated for each gene in this restricted set of adult tissues (supplementary dataset S4, Supplementary Material online). Pearson correlations for each pair of tissues and across selected genes (Z-score > 1 and TPM > 5 for at least one tissue) highlight positively correlated combinations of tissues, such as neural crest and paraxial mesoderm-derived skeletal components (Meckel's cartilage, vertebrae and chondrocranium), central nervous system organs (brain, eye, spinal cord), striated muscles (hypaxial muscle, heart), mineralized tissues (dental lamina, skin denticles), tissues involved in the maintenance of osmotic homeostasis (rectal gland and kidney), and hematopoietic organs or cells (spleen, blood) (supplementary fig. S4, Supplementary Material online). Some unexpected correlations were retrieved between ovary, pancreas, stomach, and liver. Because these organs all displayed a high proportion of genes with lower levels of expression, we interpret these correlations as an artifact due to shared low expression levels for a large set of genes. Another correlation was observed between gills, dental lamina, and skin denticles that may come from all three tissues being of epidermal origin, in addition to common processes of mineralized

<figure><img src="image_11.png" /><figcaption>**Fig. 4.** Transcriptomic analysis of small-spotted catshark adult tissues. a) Scheme showing the adult tissues included in the transcriptomic analysis. (b), (c), (d), and e) graphs showing Z-score and TPM values for each gene in the eye (b), heart (c), blood (d), and pancreas (e). Each dot represents one gene and the graph is restricted to genes with a Z-score > 1 and TPM > 5 in the tissue considered. GO-terms relevant to known functions of these tissues are indicated below the graphs. Arrows point to a selection of genes which exhibit high Z-scores in the catshark, consistent with selective expressions either in the eye (b), heart (c), blood (d), or pancreas (e) as documented in bony fishes.</figcaption></figure>
structure development shared by skin denticles and teeth (Debiais-Thibaud et al. 2011; Debiais-Thibaud et al. 2019). For all tissues analyzed, we retrieved genes possessing both high Z-scores and high TPM values, albeit with very different distributions along these two parameters depending on the tissue (Fig. 4b to e; supplementary fig. S4, Supplementary Material online). Among those with the highest Z-scores, we identified many known markers of homologous organs or cell populations in bony fishes (Fig. 4b to e; see also supplementary fig. S4, Supplementary Material online). For organs still poorly characterized, such as the rectal gland (a salt-secreting osmoregulatory organ), placoid skin denticles (dermal scales), or the ampullae of Lorenzini, we establish lists of candidate gene signatures, which provide a starting point to examine their molecular functions (supplementary fig. S4, Supplementary Material online). As expected, GO-term enrichment analysis for each one of these tissues led to the identification of terms relevant for known functions of these tissues (Fig. 4; supplementary dataset S5, Supplementary Material online). In summary, despite some possible paralog misidentifications, unavoidable in automatic annotations, these transcriptomic data yield coherent molecular blueprints of selected organs and provide a chondrichthyan reference to explore ancient organ molecular signatures across jawed vertebrate as well as chondrichthyan specificities.

## Characterization of Sensory Organs and Cells

*The Anatomy of Sensory Organs in Three Dimensions*

To provide a comprehensive morphological reference of sensory organs and their 3D organization, we analyzed the head of a late embryo (pre-hatching stage, 8 cm total length) by combining diffusible iodine-based contrast enhancement with propagation phase-contrast synchrotron radiation micro-computed tomography (DICE-PPC-SRµCT [Leyhr et al. 2023]). The reconstructed scan images were manually segmented to produce 3D renderings of the soft tissues (Fig. 5a to f). The volume of the whole head was measured to approximately 283 mm³. The eyes, nasal cavities, neuromast canals of the lateral line, and ampullae of Lorenzini (including their canals) altogether comprise 14% of the head total volume in this late embryonic stage (Fig. 5a to c), measuring 22.8 mm³, 8.6 mm³, 0.9 mm³, and 6.0 mm³, respectively. The sensory surface area of the ampullae of Lorenzini (approximated by excluding the canal surface, Fig. 5f) was 58.2 mm², retinas were 46.6 mm² (Fig. 5c) and olfactory epithelia (surface of the olfactory rosettes, Fig. 5e) were 146.3 mm². The ampullae of Lorenzini therefore occupy a volume comparable to the olfactory organs, and have a sensory surface area comparable to that of the visual system.

The 3D visualization of the 219 left ampullae of Lorenzini and their canals (Fig. 5a to c) recovers the superficial ophthalmic (SO), buccal (BUC), and mandibular (MAN) clusters (following Norris 1929; Rivera-Vicente et al. 2011). We can refine this clustering to 17 subclusters defined by distinct canal orientations and/or surface pore

locations within the SO (10) and BUC (7) clusters (Fig. 5b and c). Through these subclusters, sensory ampullae within both the SO and BUC clusters are connected to surface pores on the anterior, dorsal, and lateral/ventrolateral surfaces of the head, with the greatest concentration of pores on the anteroventral surface of the snout (Fig. 5a to c, supplementary fig. S5, Supplementary Material online). The 3D identification of the canal, the ampullae, and their sensory alveoli, with correspondence to histological observations is similar to data obtained in adult individuals despite the difference in size (i.e. 150 to 200 µm alveolar diameter at pre-hatching and juvenile stage (Fig. 6) versus 750 to 800 µm alveolar diameter in adults [Crooks and Waring 2013]). The number of alveoli was correlated with the size of the ampullae in our sample, suggesting the number of alveoli observed in a pre-hatching might not be representative of the adult condition.

The number of observed olfactory lamellae (~33 per rosette, Fig. 5e) is consistent with the number in adult individuals (Ferrando et al. 2017) which is low relative to most other chondrichthyans (Ferrando et al. 2017; Dymek et al. 2021; Clark et al. 2022). Similar to the adult olfactory organ, each olfactory rosette is partially separated by the central raphe into two compartments: one located dorsoanterior, the other ventroposterior. A large number of nerve bundles from individual lamellae run either lateral or medial to the olfactory rosette, cross the chondrocranium border through fibrous fenestrae, and merge into two short olfactory nerves (a medial and a lateral component), which reach the olfactory bulb medially and laterally in the anterior braincase (arrowhead, Fig. 5f, h, and k, supplementary figs. S6 and S7, Supplementary Material online, supplementary datasets S6 and S7, Supplementary Material online). The catshark olfactory bulb is bipartite, with two separate zones for olfactory nerve layer and glomerular layer, but a common granule cell layer (supplementary fig. S7, Supplementary Material online, supplementary datasets S6 and S7, Supplementary Material online). The virtual sections confirm that medial nerve obtains input from medial lamellae, whereas the lateral nerve collects the output from the more laterally situated lamellae of the olfactory epithelium, each from both halves of the olfactory organ (dorsoanterior and ventroposterior, supplementary datasets S6 and S7, Supplementary Material online).

Three major branches of the anterior lateral line nerve can be delineated: the superficial ophthalmic, maxillary (or buccal), and mandibular (following Boord and Campbell 1977). By visualizing the anterior lateral line nerves and the sensory ampullae without their canals (Fig. 5f), it is apparent that the SO, BUC, and MAN clusters are innervated by the superficial ophthalmic, maxillary, and mandibular nerves respectively, with no visible overlap or exceptions. One projection of the superficial ophthalmic nerve extends ventrally (black arrow, Fig. 5f), and one maxillary nerve extends anteriorly (white arrow, Fig. 5f) but these interact exclusively with the anterior lateral line neuromasts, not the nearby BUC or SO electro sensory ampullae, respectively (similar to observations in *Squalus acanthias* [Norris and Hughes 1920; Norris
<figure><img src="image_12.png" /><figcaption>**Fig. 5.** 3D organization of head sensory organs in a pre-hatchling *Scyliorhinus canicula*. a) Lateral view of whole-head surface rendering, with dermal denticles in white and the neuromast canal (NM) pores in black. Surface pores for the ampullae of Lorenzini (AoL) canals are displayed in shades of red. b) Rendering of subclusters of AoL as determined by canal orientations, shown in a range of different colors (neuromast canals in white). c) Lateral 3D rendering of sensory organs. Shades of red = AoL with canals colored for three clusters as defined by the location of the ampullae: orange = superficial ophthalmic (SO); red = buccal (BUC); old rose = mandibular (MAN); white = neuromast canals (NM); green = olfactory rosette (OR); blue = retina; dark blue = lens. d) Medial rendering of the anterior snout (anterior to the left), isolating a single neural projection of the superficial ophthalmic nerve (light green) to both the neuromast canals (white) and several AoL with canals and sensory ampullae (orange = internal cast, and its superficial transparent yellow lining = epithelium). e) Anteroventral rendering of the isolated left OR. f) Lateral rendering displaying the interaction between sensory organs and the nervous system. Lilac = brain, yellow = internal fibrous surface of chondrocranium, light green = superficial ophthalmic nerve (son), light blue = maxillary nerve (mn), orange = mandibular nerve (mbn), white = neuromast canals, red = ampullae of the AoL system. White arrowhead indicates olfactory nerve foramen leading to the olfactory bulb (ob), black asterisk indicates optic nerve foramen, white arrow indicates anterior projection of the maxillary nerve, black arrow indicates ventral projection of the superficial ophthalmic nerve. g, h) Virtual parasagittal and horizontal histological sections with left sensory organs colored as in (c), white arrowhead locates axon projections of the olfactory neurons emerging from the sensory epithelium (compare to the K panel). The dotted line in h) indicates the position of the sagittal section displayed in (g). i) Parasagittal histological section in plane corresponding to (g). tel = telencephalon, di = diencephalon, rh = rhombencephalon, in = inner ear, Mc = Meckel's cartilage. j) Histological detail in the neuromast canals and AoL. k) Histological detail in the olfactory epithelium showing secondary folding (sf) of the sensory epithelium on each lamella and axon projections of the olfactory neurons (oln, white arrowhead) emerging from the sensory epithelium.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 6. Histological (HES) staining of the sensory ampullae of Lorenzini in a juvenile catshark in their longitudinal (a) and transverse b) sections. Gene expression patterns for *otof* (c, d), *slc1a3* (e, f) and *vglut3* (g, h) in similar orientations. (i, j) 3D renderings of a sensory ampulla (red = hydrogel-filled cavity, transparent yellow = epithelium) in similar orientations.</figcaption></figure>

1929]). Other anterior lateral line nerve bundles innervate both the neuromast canal and ampullae of Lorenzini sensory bulbs (Fig. 5d). Virtual histological sections (Fig. 5g and h) provide comparable morphological detail to traditional stained histological sections (Fig. 5i to k). Overall, the description of sensory organs in a pre-hatching stage resembles the adult situation despite the overall smaller size of organs, suggesting that sensory organs mostly display their final functional organization prior to hatching.

## The Sensory Epithelium of the Ampullae of Lorenzini

Several candidate genes for the ampullae of Lorenzini have been identified by differential expression in RNA-seq data in another shark, *Scyliorhinus retifer* (Bellono et al. 2017) and in a bony fish, the paddlefish *Polyodon spathula* (Modrell et al. 2017). Similar to the study by Bellono et al. (2017), we found that the *parvalbumin alpha* (*pvalb*) gene is one of the most highly expressed transcripts with high specificity to the ampullae of Lorenzini (**supplementary dataset S8, Supplementary Material online**). In addition, we uncovered the expression of several ion channels previously associated with electroreception in elasmobranch and/or paddlefish electrosensory cells (**supplementary dataset S9, Supplementary Material online**): the L-type voltage-gated calcium channel *Cav1.3* (*cacna1d*, XM_038812023.1), the calcium-gated potassium channel *BK* (*kcnma1*, XM_038783287.1), and 2 voltage-gated potassium channel subunits (*kcnab3*, XM_038786707.1; *kcnai*, XM_038811203.1 [**supplementary dataset S10, Supplementary Material online**]). Two other ion channels involved in signal transduction in these sensory cells are also found to be expressed in the ampullae of Lorenzini tissue, albeit without strong specificity (Z-score<1: *cacnb2* XM_038798048.1 and *atp1b2* XM_038786982.1). Only *kcnab3* displayed both specific and strong expressions in the ampullae of Lorenzini (Z-score>5, **supplementary dataset S8, Supplementary Material online**). We found strong expression

(TPM >13) in the ampullae of Lorenzini (**supplementary dataset S9, Supplementary Material online**) for several transcription factors, including *mafa* (XM_038808833.1), *atoh1* (XM_038793060.1), *pou4f3* (XM_038796138.1), and *six2* (XM_038813278.1) that were previously associated to lateral line organs (Modrell et al. 2017; Minařík et al. 2024). Finally, 2 genes known to code for presynaptic proteins of electrosensory cells in bony fish (Modrell et al. 2017), *Otoferlin* (*otof*) and *Vglut3* (*slc17a8*), had both very high specificity and levels of expression (**supplementary dataset S8, Supplementary Material online**). Another highly expressed and specific transcript coding for a glutamate transporter, *slc1a3*, was also identified (**supplementary dataset S7, Supplementary Material online**). We confirmed that expression of these last 3 genes was specifically located in cells of the sensory epithelium of the ampullae of Lorenzini, excluding the central cup (see Fig. 6).

## The Olfactory Epithelium: Supporting Cells Versus Olfactory Receptors

Olfactory sensory neurons (OSNs) detect odorant molecules and pheromones by means of G protein-coupled receptors (GPCRs). Extensive research has uncovered the genetic variability underlying the adaptive evolution of these olfactory receptor families in OSNs, from their origin and evolution in vertebrates to more recent events of gene family expansion or contraction (Saraiva and Korschning 2007; Grus and Zhang 2009; Bear et al. 2016; Dieris et al. 2021; Policarpo et al. 2022). Examination of the 50 most abundantly and specifically expressed genes from the "brain & olfactory epithelium" library (**supplementary dataset S8, Supplementary Material online**) did not identify genes encoding GPCRs, but highlights 2 genes known to possess an olfactory role in bony vertebrates, *trpc2* and *s100z*, and 2 duplicates of *moxD2* that we named *moxd2.1* and *moxd2.2*. The vertebrate *moxd2* gene encodes
<figure><img src="image_14.png" /><figcaption>Fig. 7. Histological (HES) staining of the olfactory epithelium in a juvenile catshark. a) with higher magnification (a2, a3) in transverse section. Gene expression patterns for *moxd2.1* (b), *s100z* (c) and *trpc2* (d) in parallel sections. ecm: extracellular matrix of the basal conjunctive tissue; lu: lumen of the olfactory rosette; oln: olfactory nerve; dotted line: basal lamina between the sensory epithelium and the underlying conjunctive tissue; asterisk in panel a3: putative ionocyte (large clear cell); open arrowhead: putative sensory cell bodies; black arrow head: putative supporting cells.</figcaption></figure>

a monooxygenase named "Monooxygenase Dopamine-Beta-Hydroxylase Like 2", of unknown but possible olfactory function (Hahn et al. 2007). Our genome screening and phylogenetic reconstruction recovered one *moxd2* gene in the great white shark genome and none in other available chondrichthyan genomes, but the two duplicates in the small-spotted catshark have a time of divergence older than that of the last common ancestor of these two shark species (supplementary fig. S8, Supplementary Material online). Several duplicates were also found in *Petromyzon marinus*, *Protopterus annectens*, and *Erpetoichthys calabaricus*, but all duplication events appeared specific to each one of these three branches (supplementary fig. S8, Supplementary Material online). Both *moxd2* copies were found located in a single chromosomal region in *Protopterus annectens* and *Scyliorhinus canicula* (Chr8.part0 and Chr16, respectively), separated by other genes. The two genes found between the catshark *moxd2* duplicates (namely *emg1* [XM_038773333] and *peptidase inhibitor 1G-like* [XM_038773822]) have orthologs on the Chr8.part0 of *Protopterus annectens* but the former is located more than 300Mb upstream and the latter 1.5Gb downstream of the *moxd2* genes. Little similarity between these loci in these two species and the phylogenetic relations of these sequences (supplementary fig. S8, Supplementary Material online) support a scenario of

parallel tandem duplication events in these two lineages. Both catshark *moxd2* paralogs showed identical expression patterns in the most superficial cell layer of the olfactory epithelium, known to consist of supporting cells (Fig. 7a and b and supplementary fig. S9, Supplementary Material online for *moxd2.2* expression; Ferrando et al. 2010).

Phylogenetic analysis identified a single *s100z* gene in each tested vertebrate species including the catshark (supplementary fig. 10, Supplementary Material online), and genomic data support conserved synteny between human and catshark for the genes neighboring *s100z*. The *s100z* gene encodes a protein that belongs to the S100 family of calcium-binding proteins and is known to be expressed in a large population of microvillous OSNs in zebrafish (Oka et al. 2012). Pseudogenization of *s100z* in mammals was found in species with a reduced or absent vomeronasal organ, which houses the microvillous OSNs (Hecker et al. 2019). Here we show that *s100z* expression is located in an intermediate layer of the catshark olfactory epithelium (Fig. 7c), known to host the sensory cell bodies (Ferrando et al. 2017; Syed et al. 2023). A similar expression pattern was observed for the *trpc2* gene (Fig. 7d), which is known to be expressed in microvillous OSNs in the catshark (Syed et al. 2023). The *s100z* cell population seems less contiguous compared to *trpc2*-expressing cells and may delineate a neuronal subpopulation.
## Rise and Fall of Paralogs in Sensory Gene Families in Chondrichthyans

We examined the evolutionary dynamics, gene-wise and expression-wise, of several candidate gene families known to play major roles in the respective sensory organs: the four main olfactory receptor families for the olfactory epithelium, opsins and crystallins for the eye, and the transient-receptor potential (TRP) ion channel family essential for diverse sensory functions. For each family, we established the catshark gene repertoire, showed its

evolution in relation to other vertebrates using systematic phylogenetic analyses for unambiguous resolution of orthology relationships, and analyzed tissue-specific expression patterns for selected family members.

### Small Repertoire Sizes of Olfactory Receptor Gene Families

The olfactory sense of vertebrates employs four major olfactory receptor families for the detection of odors. They are the odorant receptors (OR), the trace amine-associated receptors (TAAR), the olfactory receptors class A-related (ORA) also

<figure><img src="image_15.png" /><figcaption>**Fig. 8.** Gene expression patterns for *or3* (a-c), *taar2* (d-f) and *v2rl1* (g-m) in transverse sections of a juvenile catshark (except in h, i: stage 31 embryo) in the olfactory epithelium (OE; a, d, g) or non-olfactory sites: ampullae of lorenzini (AoL; b, c), c: canal of an AoL, se: sensory epithelium of an AoL; differentiated e) and undifferentiated retina (h), ac: amacrine cells; gc: ganglion cells, hc & bc: horizontal and bipolar cells; INL: inner nuclear layer, IPL: inner plexiform layer, ONL: outer nuclear layer, OPL: outer plexiform layer, pr: photoreceptor cell layer; posterior brain (f, i), arrowhead: posterior hypothalamus, in: inner ear, v: ventricle of the mesencephalon; taste bud (j); epiderm and neuromast (NM; k), sb: scale bud; gills (l) and muscle tissue (m). Black arrowheads point to scattered epithelial cells positive for gene expression; open arrow heads point to endothelial (i) or muscle (m) cells.</figcaption></figure>
called vomeronasal receptor type 1 (V1R), and the vomero-
nasal receptors type 2 (V2R, in bony fish, also called olfactory
receptor class C, OlfC). These receptor families were first iden-
tified in mammalian species before being found later in bony
fish, and all four have also been detected in cartilaginous and
jawless fishes (Nei et al. 2008; Sharma et al. 2019; Dieris et al.
2021; Kowatschew and Korsching 2021; Syed et al. 2023). The
ORs are the largest olfactory receptor family of bony verte-
brates (e.g. 1,948 genes in the elephant *Loxodonta africana*
(Niimura et al. 2014)), but the catshark has only 2 canonical
OR and 6 OR-like genes (Syed et al. 2023). The ORA family
shows large gene expansions in tetrapods (here named
V1Rs, Young et al. 2010), the TAAR family is massively ex-
panded in teleost fish (Korsching 2020), similar to the
TAAR-like (TARL) family in lampreys (Dieris et al. 2021),
but these families consist of only 3, 3, and 2 genes respectively
in the catshark (Syed et al. 2023). The restricted diversity of
orthologs in these gene families in chondrichthyans was
shown to be linked to the loss of function of these receptors
in the adult olfactory epithelium (Syed et al. 2023). Here we
examined the potential alternate sites of expression of these
receptors.

No specificity of expression toward the “brain & olfac-
tory epithelium” tissue was detected for either of the
two canonical OR genes, or1 and or2, nor the six OR-like
genes (or3 to or9, excluding or7 that is not a predicted
gene in the current version of the genome, but see Syed
et al. 2023), supporting the loss of olfactory function in
the catshark. For or3 we have confirmed the absence of ex-
pression in the olfactory epithelium by in situ hybridiza-
tion (Fig. 8a). However, six of the eight OR genes showed
highly specific expression in other tissues in the transcrip-
tomic data (supplementary fig. 1, Supplementary Material
online1). Among these, or3 and or5 had their highest levels
of expression in the ampullae of Lorenzini. In line with this,
or3 expression was detected in the cells lining the canals
but not in the sensory zones of the ampullae (Fig. 8b
and c). Transcripts of or2 are broadly identified across tis-
sues but over-represented in the testis (supplementary fig.
11, Supplementary Material online). Both or1 and or4 are
expressed with high specificity in the spleen, and or9 in
the stomach (supplementary fig. 11, Supplementary
Material online). The specific albeit low expression of or9
in the stomach may be related to an expression restriction
to a small subpopulation of cells (supplementary fig. 11,
Supplementary Material online), conceivably involved in
sensing metabolic status (Kitamura et al. 2014). The two
remaining genes, orG and or8, show overall weak and rather
unspecific expression and may exert their function in other
tissues or developmental stages than those examined here
(supplementary fig. 11, Supplementary Material online).

Similarly, none of the genes of the ORA, TAAR, or
TARL show a high expression specificity for the “brain &
olfactory epithelium” in our transcriptomic data
(supplementary fig. 11, Supplementary Material online),
but some exhibit strikingly specific expression in other tis-
sues, e.g. *taar2* in the eye and *tarl4* in the spiral intestine. In
situ hybridization with a *taar2* probe showed dense

labeling in several brain regions and layered expression in
the eye at both borders of the inner and the outer nuclear
layer (layers that include photoreceptors, horizontal and
bipolar cells, amacrine cells, and ganglion cells, Fig. 8e
and f). We also detected sparse labeled cells in the juvenile
olfactory epithelium, suggesting a conserved function in
olfaction, restricted to the juvenile stage of development
(Fig. 8d) as *taar2* expression was not detected in adult ol-
factory tissues (Syed et al. 2023).

V2R is the main olfactory receptor family in the catshark (Sharma et al. 2019; Syed et al. 2023), and here we examined the 32 canonical v2r genes annotated in the gene model as well as three v2r-like genes. Twenty-one v2r genes show expression specificity in the “brain & olfactory epithelium” sample (supplementary fig. 12, Supplementary Material online), consistent with their known expression in adult OSNs (Syed et al. 2023). Six genes show enriched expression in “brain & olfactory epithelium” but their maximal Z-score is in another tissue (v2r2, v2r5, v2r9, v2r11, v2r29, and v2r33) and four genes are not enriched in the “brain & olfactory epithelium” sample (v2r4, v2r10, v2r12, and v2r28). These data suggest that the large majority of v2r genes indeed have expression restricted to the olfactory sensory neurons, while others are additionally or exclusively expressed in other organs suggesting additional functions. For instance, a group of three closely related genes, v2r9, v2r10, and v2r11 (Syed et al. 2023), also have enriched expression in the rectal gland (supplementary fig. 12, Supplementary Material online). This suggests that the v2r9, v2r10, v2r11 clade may have undergone neofunctionalization in the rectal gland, possibly in ionic homeostasis (Burger and Hess 1960). Along the same line, three phylogenetically distant genes (Syed et al. 2023), v2r2, v2r5, and v2r12 are highly enriched in the testis, suggesting three independent neofunctionalization events. Two closely related genes, v2r28 and v2r29, have enriched expression in some embryonic stages (22 to 26) and “brain & olfactory epithelium” tissue (supplementary fig. 12, Supplementary Material online), although both with low levels of expression. Stages 22 to 26 correspond to early organogenesis and indeed the olfactory placode, visible already at stage 20, matures into a pit during these stages (Ballard et al. 1993), suggesting these genes may have acquired a function in the early development of the sensory olfactory placode. Several pairs of paralogous genes show very different expression specificities, suggesting that expression regulation is highly dynamic in this family (e.g. v2r9/v2r10, and v2r33/v2r35a). As a consequence, despite the v2r gene family showing the typical expansion seen for olfactory receptor gene families, their function is not restricted to the olfactory system and may involve such different sites as testis and rectal gland. Finally, in the v2r-like gene family, only v2r14 has a higher level of expression in the “brain & olfactory epithelium” tissue, but also displays enriched expression in gills (supplementary fig. 11, Supplementary Material online). The other two genes, v2rl1 and v2rl3, have their higher expression level first in gills, with v2rl1 also more highly
expressed in the esophagus, skin denticles, and dental lamina (supplementary fig. 11, Supplementary Material online). Because the dental lamina, gills, and anterior esophagus are places where taste buds are found, we wished to better localize gene expression in relation to taste buds. *In situ* hybridization for *v2rl1* and *v2rl4* showed virtually identical sites of expression, but no specific expression in taste buds, rather in cells scattered throughout the epithelium (Fig. 8j and supplementary fig. 13, Supplementary Material online). These genes also display expression in several brain regions (Fig. 8i and supplementary fig. 13, Supplementary Material online). Expression is also detected in the eye, intervertebral tissue, some chondrocytes, some muscle, and gills, illustrating a wide array of expression sites during embryonic development and in the juvenile stage (Fig. 8 and supplementary fig. 13, Supplementary Material online). Unexpectedly, *v2rl1* and *v2rl4* expression was also detected in sparse cells in the olfactory epithelium that may be OSNs (Fig. 8g and supplementary fig. 13, Supplementary Material online). Taken together, our study points to several independent heterotopic shifts of expression in the *v2r* and *v2r-like* families, with noteworthy diverse expression patterns of the *v2r-like* gene family, which suggests a broad range of functions in this clade of receptors.

## Contracted Visual and Non-visual Opsin Gene Families

Visual opsins are light-activated receptors that are conserved in all animal groups. Together with their sister groups of non-visual opsins and non-canonical opsins, they belong to a wider group of GPCRs that are involved in signal transduction in animal cells (Koyanagi et al. 2021). Vertebrate opsins are classically categorized into several opsin subfamilies (Terakita 2005) and many lineage-specific duplications and

losses were identified over the course of vertebrate diversification (Rennison et al. 2012).

In the small-spotted catshark genome, we identified two genes for visual opsins (supplementary fig. 14, Supplementary Material online): the one encoding Rhodopsin (*rho*) with extremely high and exclusive expression in the eye (supplementary fig. 15, Supplementary Material online), and the second visual opsin gene *rh2*, which has very low levels of expression in all sampled tissues (supplementary fig. 15, Supplementary Material online). *In situ* hybridization of *rh2* showed expression in the poorly differentiated retina of a stage 31 embryo (Fig. 9c) but not in a juvenile (Fig. 9b), supporting a restricted transient expression in early retina progenitors. Rhodopsin therefore appears as the single visual opsin involved in light-sensing in the adult small-spotted catshark. The other three known jawed vertebrate visual opsin genes (*sws1*, *sws2*, and *lws/mws*, expressed in vertebrate cone cells) are absent in the catshark: *sws* genes are lost in all chondrichthyans while *lws/mws* loss is more recent as it is still found in some shark species (Yamaguchi et al. 2021). In addition to the isolated presence of the Rhodopsin visual pigment in the adult catshark retina, only three G-protein coding genes were recovered in our top-50 Eye list (supplementary dataset S8, Supplementary Material online): *gnat* (G-protein, alpha transducing activity polypeptide 1), *gngt1* (G-protein, gamma transducing activity polypeptide 1) and *gnb1* (G-protein G(I)/G(S)/G(T) subunit beta-1) that are considered specific to signal transduction in rod cells only (Peng et al. 1992), congruent with the previous finding that the small-spotted catshark adult eye only has rod cells (Bozzanao et al. 2001) with the single rhodopsin pigment.

Jawed vertebrates also have three canonical non-visual opsins: Pinopsin, VA-opsin, and Parietopsin. The small-

<figure><img src="image_16.png" /><figcaption>**Fig. 9.** a) Histological (HES) staining of a transverse section in the head of a juvenile catshark: general view (a1) with inset (a2) on the retina histology. Gene expression patterns for *rh2* (b, c) and *pinopsin* (d, e; arrowhead points to expression in the pineal stalk): in a juvenile (b, d); in the developing retina of a stage 31 embryo (c); in an adult catshark (e). Scales in microns.</figcaption></figure>
spotted catshark has all but the parietopsin gene (supplementary fig. 14, Supplementary Material online), similar to other chondrichthyans (Yamaguchi et al. 2021). In the catshark, the gene coding for VA-opsin is only found expressed at very low levels but specifically in the eye, coherent with the putative light-sensing function of VA-opsin in retinal interneurons (Jenkins et al. 2003). The gene coding for Pinopsin is expressed quite ubiquitously at low levels with its zone of highest expression in the ampullae of Lorenzini (TPM = 8.5; Z-score = 2.88; supplementary fig. 15, Supplementary Material online). However, in situ hybridizations failed to detect the expression of *pinopsin* either in the ampullae of Lorenzini or in the retina of a juvenile individual (Fig. 9d2). Instead, a highly specific expression was detected in the pineal organ of both a juvenile and an adult specimen (Fig. 9d1 and e) in accordance with photoreceptive abilities of the pineal gland in vertebrates (Okano and Fukada 2001).

Non-canonical opsins include several subfamilies (Terakita 2005). In the Melanopsin subfamily, jawed vertebrates have up to two genes (*Opn4* and *Opn4x*). Only *opn4* is found in the small-spotted catshark (supplementary fig. 16, Supplementary Material online), similar to other elasmobranchs (Yamaguchi et al. 2021). The catshark *opn4* shows strong and specific expression in the eye (supplementary fig. 15, Supplementary Material online), supporting a conserved function in vision (Díaz et al. 2017). Within the Encephalopsin/TMT (teleost multiple tissue)-opsin subfamily (Kato et al. 2016), the catshark genome includes the encephalopsin (*opn3*) gene and a TMT1 gene (*tmt1*, supplementary fig. 17, Supplementary Material online), similar to other elasmobranchs (Yamaguchi et al. 2021). Encephalopsin expression is highest in the testis, then in embryonic stages and the “brain & olfactory epithelium” samples (supplementary fig. 15, Supplementary Material online). Encephalopsin is named after its expression that was first described in brain tissues (Blackshaw and Snyder 1999), and extensive brain expression is known at embryonic stages (Davies et al. 2021). A function of opsins in sperm thermotaxis was demonstrated in the mouse, where Encephalopsin was the most highly expressed of all opsins in spermatozoa (Pérez-Cerezales et al. 2015). Both these expression patterns appear conserved in the catshark, suggesting sensing capabilities by opsins in catshark sperm cells. On the other hand, the catshark *tmt1* (multiple-tissue opsin) is primarily expressed in the eye but at low levels, and its second highest site of expression is the spinal cord (supplementary fig. 15, Supplementary Material online). In the Neuropsin subfamily (Opn5/6), up to 7 paralogs were identified in jawed vertebrates (Yamashita 2020) but only one was found intact in the catshark genome (supplementary fig. 18, Supplementary Material online), which is less than most other elasmobranchs (Yamaguchi et al. 2021). Expression of this gene, *opn5L1b*, is found in skin denticles, testis, eye, and spinal cord, although with low levels of expression (supplementary fig. 15, Supplementary Material online). From synteny comparisons with other

<table><caption>**Table 2** Identified opsin gene complement (bold) in the small-spotted catshark with confirmed gene losses ('lost') or pseudogenization ('pseudogene')</caption><thead><tr><th></th><th></th><th># Of<br>genes</th><th>Sites of expression</th></tr></thead><tbody><tr><td>Canonical opsins</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="2">Visual opsins</td><td>rho</td><td>1</td><td>Eye</td></tr><tr><td>rh2</td><td>1</td><td>Transient expression in the<br>embryonic retina*</td></tr><tr><td rowspan="3">Pinopsin</td><td>sws1</td><td>0</td><td>lost</td></tr><tr><td>sws2</td><td>0</td><td>lost</td></tr><tr><td>lws/mws</td><td>0</td><td>lost</td></tr><tr><td></td><td>pinopsin</td><td>1</td><td>Juvenile and adult pineal<br>gland *</td></tr><tr><td>VA-opsin</td><td>va-opsin</td><td>1</td><td>Eye</td></tr><tr><td>Parietopsin</td><td>...</td><td>0</td><td>lost</td></tr><tr><td rowspan="2">Melanopsin</td><td>opn4</td><td>1</td><td>Eye</td></tr><tr><td>opn4x</td><td>0</td><td>lost</td></tr><tr><td rowspan="2">Encephalopsin</td><td>opn3</td><td>1</td><td>Testis, embryonic stages,<br>brain and olfactory<br>epithelium</td></tr><tr><td>tmt1</td><td>1</td><td>Eye, spinal cord</td></tr><tr><td rowspan="5">Neuropsin</td><td>tmt2</td><td>0</td><td>lost</td></tr><tr><td>tmt3</td><td>0</td><td>lost</td></tr><tr><td>opn5m</td><td>P</td><td>pseudogene</td></tr><tr><td>opn5m2</td><td>0</td><td>lost</td></tr><tr><td>opn5L1a</td><td>0</td><td>lost</td></tr><tr><td rowspan="4">Peropsin</td><td>opn5L1b</td><td>1</td><td>Skin denticles, testis, eye, and<br>spinal cord</td></tr><tr><td>opn5L2</td><td>0</td><td>lost</td></tr><tr><td>opn5L2c</td><td>0</td><td>lost</td></tr><tr><td>rrh</td><td>1</td><td>Eye</td></tr><tr><td rowspan="2">Photoisomerases</td><td>rgr1</td><td>1</td><td>Eye and testis</td></tr><tr><td>rgr2</td><td>P</td><td>pseudogene</td></tr></tbody></table>

Sites of expression are defined as sites of over-represented expression in the RNA-seq data (Z-score > 1) or confirmed sites of expression by in situ hybridization (*).

elasmobranch *opn5m* genes, a pseudogenized sequence of *opn5m* was identified suggesting a recent loss of this opsin in the catshark lineage (supplementary fig. 18, Supplementary Material online).

Two other opsin subfamilies are involved in the regulation of photoreception in vertebrates: the Peropsin (Rrh) subfamily involved in regulation of all-trans-retinol uptake from photoreceptors (Cook et al. 2017), and the Photoisomerase (RGR, or Opn6) subfamily involved in the modulation of retinaldehyde levels in retinal cells (Díaz et al. 2017). We found one gene for each subfamily: the catshark *rrh* (supplementary fig. 19, Supplementary Material online) is expressed with specificity in the eye while the catshark *rgr1* expression was highest in the eye and testis (supplementary fig. 15, Supplementary Material online). A remnant sequence of a pseudogenized *rgr2* is identifiable in the catshark locus where *rgr2* is found in other elasmobranchs (supplementary fig. 20, Supplementary Material online).

Overall, these data show a strong conservation of all aspects of photoreception in the small-spotted catshark despite a restricted gene repertoire (10 opsins out of the 26 paralogs identified in jawed vertebrates (see Yamaguchi et al. 2021, Table 2).
<figure><img src="image_17.png" /><figcaption>**Fig. 10.** Simplified representation of the phylogenetic relationships between the small-spotted catshark gamma-crystallins, with the single gamma-S crystallin gene (*crygs*) as the outgroup (a) mapped on the gene organization along chromosome 2 (b); gene names appear as the last 5 digits of their genomic ID XM_0387xxxxx. TPM values in three selected samples (embryo stage 30, embryo stage 31, and adult eye) for each gene copy c) with color code indicated the associated Z-score (vermillion: Z-score > 4; orange: 4 > Z-score > 1). Histological staining of a transverse section through the eye of a juvenile catshark (d1 general view, d2 close-up on lens cells) and in situ hybridization with a gamma-crystallin (*cryg1*) RNA probe in the corresponding sample e) and in a stage 31 embryo (f). The lens epithelial cells are labeled with an arrowhead, and the elongated lens fiber cells are marked with a star: both cell types are positive for *cryg* RNA in situ hybridization.</figcaption></figure>

## Amplified Crystallin Gene Family

Vertebrate crystallins evolved with the vertebrate camera eye. Alpha-crystallins evolved from heat-shock proteins

while beta- and gamma-crystallins evolved from calcium-binding proteins that gave rise to an ancestral chordate beta-/gamma-crystallin gene (Kappé et al. 2010; Slingsby
et al. 2013; Cvekl et al. 2017). These proteins are highly expressed in lens fiber cells to ensure the transparency of, and a refraction gradient in, the lens of the vertebrate eye. Their genes were shown to have undergone rapid evolution in sequence and gene number in several vertebrate lineages in relation to ecological diversification (Slingsby et al. 2013; Malik et al. 2021) and are also linked to eye pathologies (Panda et al. 2016).

Out of the top 50 genes in our eye transcriptome (supplementary dataset S8, Supplementary Material online), sixteen were annotated as crystallin coding genes. Of these, two were alpha-crystallins (cryaa and cryab) and five were beta-crystallins (cryba1, cryba2, crybb1, crybb2, and crybb3), making this gene repertoire similar to other vertebrates (Slingsby et al. 2013) with probable high functional conservation (Ghahghaei et al. 2009). The remaining nine crystallin genes were categorized by the automatic annotation as gamma-crystallins. We searched the catshark genome for additional similar genes and identified a total of forty-one gamma-crystallin genes, showing a much higher diversity than described in mammals (Slingsby et al. 2013). To further understand the evolution of this clustered gene family, we generated a phylogenetic reconstruction of their relationships that identified well-supported vertebrate clades for one gamma-S and one gamma-N crystallin gene, both used as outgroups (supplementary fig. 21, Supplementary Material online, and Weadick and Chang 2009). The catshark gamma-S gene (crygs) was located 2 Mb away from a dense cluster of all 40 other genes located along a 1Mb portion of chromosome 2 that also includes an eya-4-like gene (Fig. 10b). An additional non-coding RNA with sequence similarity to the surrounding gamma-crystallin genes was identified on the same cluster, suggesting recent pseudogenization. Our phylogeny gives evidence for multiple series of duplication at different time scales. In particular, it shows a large chondrichthyan-specific clade and a large bony vertebrate-specific clade suggesting much of gene duplication has happened after the divergence of these two lineages. (supplementary fig. 21, Supplementary Material online). For both taxa, there are several copies in a given species that group together, showing recently generated duplicates in all lineages (the four human genes coding for gamma-crystallins A, B, C, and D group together, as a sister clade to thirty copies in the amphibian Xenopus tropicalis (supplementary fig. 21, Supplementary Material online). For the chondrichthyan clade of gamma-crystallins, we generated a nucleotide-sequence phylogenetic reconstruction for better resolution in the more recent nodes (supplementary fig. 22, Supplementary Material online): we uncovered a clade of 26 duplicates in the elephant shark Callorhinchus mili; a clade of 20 duplicates in the skate Amblyraja radiata; a clade of 16 duplicates in the great white shark; and a clade of 30 duplicates in the catshark Scyliorhinus canicula (lighter colored boxes in Fig. 10b). These 30 copies constitute a monophyletic group and therefore originated from duplication events that

occurred after the divergence from the great white shark lineage, less than 180 My ago (following TimeTree of Life, Kumar et al. 2022). Other duplicates, also located on the same 40 gene cluster, were closer to great white shark sequences, so they duplicated earlier in their common lineage (deep blue boxes in Fig. 10b). Species-related clades of gamma-crystallins in bony fishes (supplementary fig. 21, Supplementary Material online) suggest this “duplication hotspot” is a shared characteristic of jawed vertebrates. We compared the expression patterns of all copies in our transcriptomic data: two copies generated by the most ancient duplication events show a nonspecific pattern of expression, but all others had expression specific to the adult eye and/or embryonic stages 30 and 31 (Fig. 10c). Notably, most copies generated by more ancient duplications had an adult eye-biased expression, while most of the more recent copies had expression biased toward the embryonic stages, sometimes in addition to strong expression also in the adult eye (Fig. 10c). We confirmed the location of expression of these recent copies in the embryonic and juvenile lens cells by in situ hybridization (Fig. 10d to f). Tests for positive selection in the clade of these recently duplicated gamma-crystallin genes in the catshark were non-significant. Taken together, these data suggest that ancestral copies were expressed mostly in the adult eye, but that duplications, followed by partitioning of the ancestral expression pattern led to specialized sites of expression in the embryonic eye without positive selection for new protein functions.

## A New Vertebrate Transient-Receptor Potential Channels Clade

The superfamily of Transient-Receptor Potential (TRP) channels are classified in vertebrates as eight families belonging to two groups: group 1 with Trpa (Ankyrin family), Trpc (Canonical family), Trpm (Melastatin family), Trpn (nomp-C family), Trpv (Vanilloid family) and a more recently described Trpvl (Vanilloid-Like family); group 2 with Trpp (Polycystic family, that also includes Pkd1 and related paralogs, and Pkd2 and related paralogs) and Trpml (Mucolipin family) (Venkatachalam and Montell 2007; Peng et al. 2015; Himmel and Cox 2020). TRPs are structurally defined as ion channel subunits that form pores when organized in tetramers. All are permeable to cations but some are selective to calcium ions (Venkatachalam and Montell 2007). TRPs are considered polyvalent actors of sensation because of the wide diversity of physical and chemical stimuli that may activate their opening: in mammals, they were shown to function in vision, taste, and chemesthetic sensation, olfaction, hearing, touch, nociception, and thermo- and osmosensation (reviewed in Venkatachalam and Montell 2007). They function in the plasma membrane of excitable cells (e.g. depolarizing sensory cells) but also calcium signaling in non-excitable cells by modulating calcium movement through the plasma membrane or through intracellular membranes enclosing calcium reservoirs (reviewed in Gees et al. 2010).
<figure><img src="image_18.png" /><figcaption>**Fig. 11.** a) Phylogenetic relationships between the Trpc, Trpv, and Trpa groups and the identification of a new metazoan Trpw clade: Trpc is taken as an outgroup and robust sister-relationship is recovered between the metazoan clades of Trpa and Trpw, excluding the chordate Trpv clade. Branch support is indicated by a black dot if the posterior probability value is >0.88. Species names: Ar: *Amblyraja radiata*; Cc: *Carcharodon carcharias*; Cm: *Callorhinchus milii*; Dm: *Drosophila melanogaster*; Dr: *Danio rerio*; Lo: *Lepisosteus oculatus*; Pa: *Protopterus annectens*; Pm: *Petromyzon marinus*; Sc: *Scyliorhinus canicula*; Xt: *Xenopus tropicalis*. b) scheme of the ankyrin and trans-membrane (TM) protein domains.</figcaption></figure>

c) TPM values for a selection of tissues (red underline for Z-score > 4, light orange for 4 > Z-score > 1). d) expression of *trpw1* in the telencephalon of a juvenile (arrowheads point to individual cells). e) expression of *trpw1* in the sensory layer (open arrowheads) of the olfactory epithelium of a juvenile. (f, h, l, o) Histological (HES) staining of transverse sections of a juvenile catshark through: the olfactory epithelium (f: general view with h as a close-up on the non-sensory epithelium), the gills l), and anterior vertebra (o); putative ionocytes in the olfactory epithelium and gill epithelium are indicated with open arrowhead. Gene expression patterns for *trpv8.2* in the olfactory epithelium (G and close up in (i), in an ampulla of Lorenzini (j), in the gills k), and the vertebra n) of a stage 31 embryo, and in the developing gills m) and vertebra p) of a juvenile catshark. Isolated cells with gene expression are indicated with a black arrowhead; c: ampulla canal (see Fig. 4); cc: central cup (see Fig. 4); nc: notochord, nt: neural tube; se: ampulla sensory epithelium (see Fig. 4).
Our extensive survey in the small-spotted catshark genome identified most of the expected TRP genes, but we also uncovered a previously unknown vertebrate clade of TRP. This TRP class was until now only identified in protostomes and non-bilaterian species and named Trpv1 because in previous analyses it grouped as a sister clade to the Trpv group (Peng et al. 2015; Himmel and Cox 2020). However, we recovered a sister relationship with Trpa with strong support (Fig. 11a, supplementary fig. 23, Supplementary Material online): for this reason, we suggest changing the name of the class to Trpw. In the Trpw clade, one metazoan group is recovered (including sequences from protostomes and cnidarians) where an ortholog herein named Trpw1 is also identified in Xenopus tropicalis, Protopterus annectens, Danio rerio, and Lepisosteus oculatus, in addition to their chondrichthyan orthologs. Three additional catshark copies group with one Callorhinchus milii and three Petromyzon marinus orthologs, although low support in the clade and lack of information from other jawed vertebrate clades do not allow further identification of the duplication history within this clade. The catshark paralogs in this clade were herein named trpw2.1, trpw2.2, and trpw2.3 (Fig. 11a, supplementary fig. 23, Supplementary Material online). Their predicted proteins include the expected 6 transmembrane motifs that are involved in the ion pore structure, and all have 4 repeats of an ankyrin domain which makes them structurally more similar to the Trpc family (Fig. 11b, Huffer et al. 2020). In the RNA-seq data, the newly identified trpw1 gene exhibited its highest expression in the “brain & olfactory epithelium” (Fig. 11c and supplementary fig. 24, Supplementary Material online). In situ hybridization on sections of a juvenile showed scarce punctuated signal in the telencephalon (Fig. 11d) and a highly specific expression in the intermediate cell layer of the olfactory epithelium, known to mostly consist of the sensory neuron cell bodies (Fig. 11e).

More generally in the group 1 TRPs, we uncovered all Trpc paralogs (Trpc1 to 7, supplementary fig. 25, Supplementary Material online) as well as the Trpn and Trpa paralogs (supplementary fig. 26, Supplementary Material online and supplementary fig. 27, Supplementary Material online, respectively). We recovered most of the diversity of Trpm genes (jawed vertebrate clades Trpm1 to Trpm7) although Trpm8 sequences appear to be lost in all chondrichthyan genomes (supplementary fig. 28, Supplementary Material online). The Trpv clade has an overall greater diversity of paralogs in chondrichthyans than in other jawed vertebrates (Morini et al. 2022). However, the Trpv1 and Trpv2 sister paralogs are restricted to sarcopterygian species, and only one gene was found orthologous to these duplicates in chondrichthyans and in actinopterygians, herein named trpv1/2 (supplementary fig. 28, Supplementary Material online, following Morini et al. 2022). The Trpv5/6 sequences have diversified in parallel through lineage-specific duplications in jawed vertebrates (Flores-Aldama et al. 2020; Morini et al. 2022), with a total of 3 chondrichthyan Trpv5/6 copies (herein named trpv5/G.1, trpv5/G.2, and trpv5/G.3; supplementary

fig. 29, Supplementary Material online). In our transcriptomic data (supplementary fig. 24, Supplementary Material online), the three catshark Trpv5/6 paralogs have their highest expression (Z-score > 1) in gills, but also “brain & olfactory epithelium” and seminal vesicle (trpv5/G.2) or kidney (trpv5/G.1), supporting the hypothesized function in ion transport for secretion (Flores-Aldama et al. 2020). We also identified the 2 additional jawed vertebrate paralogs that are the sister groups to the Trpv5/6 clade (Trpv7 and Trpv8 following Morini et al. 2022), one of them being duplicated in chondrichthyans (paralogs named trpv8.1 and trpv8.2 in the catshark; supplementary fig. 29, Supplementary Material online). All 3 genes have their highest expression in the gills, but both trpv8.1 and trpv8.2 also have high expression in the ampullae of Lorenzini and “brain & olfactory epithelium” tissues, suggestive of sensory functions (supplementary fig. 24, Supplementary Material online). As these paralogs were only recently identified (Morini et al. 2022), we chose to characterize the expression pattern of one of them (trpv8.2): its expression was pervasive in the gill epithelium (Fig. 11i and k), congruent with expression by cells involved in ion transport. It is also expressed in the olfactory epithelium with expression in few cells of the sensory epithelium (thickened zone of the lamellae) but more cells in the non-sensory part of the epithelium (thinner epithelium between two lamellae, Fig. 11f and g), similar to the known distribution of ionocytes (Dymek et al. 2021). Expression in the ampullae of Lorenzini also suggested expression in a scarce population of cells that seem to be located in the zone intermediate between the central cup and the sensory epithelium (Fig. 11h). Further cellular characterization will be necessary to identify putative ion transport in this zone. We also detected expression of this trpv8.2 gene in the neural system (developing and juvenile spinal cord) and in a fibrous zone of the vertebral body (Fig. 11l and n).

In the TRP group 2, we recovered all three paralogs of the Trpml (Mucolipin) genes, and identified a recent duplication in the catshark lineage that generated two duplicates of trpml3 (named trpml3.1 and trmp13.2; see supplementary fig. 30, Supplementary Material online). The expression of trpml3.1, but not of trpml3.2, appeared higher in the ampullae of Lorenzini and the “brain & olfactory epithelium” tissues (supplementary fig. 24, Supplementary Material online). We also identified catshark orthologs for each three known jawed vertebrate Trpp paralogs (pkd2, pkd2L1, and pkd2L2, supplementary fig. 31, Supplementary Material online with gene nomenclature following England et al. (2017)), with only pkd2 showing over-expression in the ampullae of Lorenzini (supplementary fig. 24, Supplementary Material online). Despite poor current knowledge of the evolution of the Pkd1-related group (England et al. 2017), we uncovered a catshark ortholog for each of the jawed vertebrate paralogs pkd1, pkd1L1, pkdrej, and pkd1L2 genes (supplementary figs. S32 and S33, Supplementary Material online) but none showed biased expression toward sensory organs (supplementary fig. 24, Supplementary Material online). Two additional duplicates were found in several chondrichthyans, which we named pkd1L4.1 and pkd1L4.2
because they were found both outside of the known Pkd1L2 and Pkd1L3 clades (supplementary fig. 33, Supplementary Material online). *pkd1L4.2* had very low levels of expression, while *pkd1L4.1* showed higher expression in a range of tissues, including the eye, ampullae of Lorenzini, and spinal cord (supplementary fig. 24, Supplementary Material online).

Overall, the proper identification of paralogs in TRP gene families and their associated pattern of expression in RNA-seq data highlight several potential actors of sensory and ion-trafficking functions that could not have been identified from a simple analogy with what is known in classical model organisms.

## Discussion

### High-Quality Datasets of Diverse Nature Fuel Integrative Analyses of Morphological and Molecular Evolution

The catshark has attracted attention as an experimental chondrichthyan model for more than a century (Coolen et al. 2008), but the lack of a high-quality genome has been a major limitation to its development as a model. Here we report the availability of a chromosome-level genome assembly and transcriptomic data across a broad array of adult and embryonic tissues. We integrate these data with morphological descriptions and in situ analyses to gain new insights into the evolution of sensory systems in jawed vertebrates. This approach opens the field for a plethora of cellular and organ-level analyses in this important model organism.

### Genome Evolution of “Big Genome” Organisms

Although long neglected because of their large genomes, elasmobranchs have recently seen a number of chromosome-level genomes released. These resources now allow comparative genomic analyses within this group and with other jawed vertebrates, helping to identify the genetic bases of morphological or physiological diversifications in the early stages of vertebrate evolution. The catshark genome is one of the largest sequenced in chondrichthyans but has comparatively fewer chromosomes. This was suggested to reflect chromosomal organization traits of ancestral jawed vertebrates (Marlétaz et al. 2023), but may also result from secondary fusions during the recent evolution of Carcharhiniforms, a hypothesis that could be further tested by synteny comparisons across a denser sampling of selachian genomes. The transposable element (TE) content in the catshark represents two-thirds of the total genome, with a strong representation of LINE retrotransposons (34%). These results are very similar to those reported for other shark genomes, such as the great white shark (Marra et al. 2019), hammerhead shark, and mako shark (Stanhope et al. 2023), and also for the little skate genome (Marlétaz et al. 2023) where TEs account for more than half of total genomic content, and where LINEs account for approximately one-third of total genomic content. Expansion of LINEs accounts for larger genomes in other lineages of vertebrates (Tan et al. 2021), and this process

appears to have shaped elasmobranch genomes for more than 270My (last common ancestor of elasmobranch fishes in Time Tree of Life (Kumar et al. 2022)). Further studies on the evolution of cell size in fossil elasmobranch and stem chondrichthyans may help identify the timing of genome expansion, following a strategy previously used in the study of salamander genomes (Laurin et al. 2016). Such large genomes may have ecological consequences as they were suggested to impact cell size, cell division time, and therefore reproductive and developmental features of the species (Cavalier-Smith 1978; Fung and Bergmann 2023). Sharks and rays have evolved diverse modes of reproduction (from egg-laying to several forms of live-bearing strategies (Katona et al. 2023) that may be shaped to be compatible with these constraints. The availability of two sets of genomic sequences and RNA-seq data coming from two populations of the small-spotted catshark allowed us to describe a population genomic landscape for the first time in this species. The identified outlier regions may contain the genetic bases of diverging morphotypes between Atlantic and Mediterranean populations (as identified in growth patterns or tooth shape (Berio et al. 2022). Further studies using genomic re-sequencing data and a larger sample size could delineate these with more confidence.

### Integrating Scales of Observation: Genome, Molecules, Cells, Organs

Genome-based technologies have opened novel perspectives to explore gene regulatory networks, cellular phenotype, organ development, and physiology model organisms that are neither genetically nor surgically amenable. Our integration of 3 very different methods (synchrotron imaging, histological observations, and densely sampled RNA-seq data) exemplifies the potential of integrative approaches combining morphological observations at different resolutions and scales with gene-centered analyses for a comprehensive view of a complete system. The combination of high-resolution 3D reconstruction of organs at nearly cellular resolution using synchrotron imaging and 2D in situ hybridization analyses on histological sections thus gives access to gene expression features with a cellular resolution in the complete organ context. This allowed for instance to pinpoint the cell heterogeneity in the olfactory epithelium, with only a small subset of cells expressing *trpv8.2*. Our bulk RNA-seq data themselves provide an unparalleled depth of molecular characterization of organs, which opens perspectives on at least two levels of analysis. Firstly, our large diversity of twenty-five different adult tissues leads to unbiased, broad-scale identification of organ-specific gene expression signatures. This expression information is an important complement to systematic comparative genomic approaches, aimed at identifying gene gains and losses relevant to morphological and physiological diversifications. The organ signatures thus provide putative gene identifiers of highly specialized cell types to be further characterized using advanced organ-focused approaches such as single-cell RNA-seq or spatial transcriptomics. Such data are crucial to explore the
evolution of cell types within organs. Similarly, the dense transcriptomic sampling of a developmental window when major organogenesis processes take place leads to a fine-grained evaluation of the timing of underlying cell differentiation progression and concomitant gene expression dynamics. Second, the dense RNA-seq data combined with the availability of a chromosome-level genome assembly allows comprehensive pictures of the structural and functional evolution of candidate multi-gene families of interest. On the one hand, the chromosome-level genome assembly allows unambiguous identifications of gene identities through additional phylogenetic analyses (as automated annotation occasionally leads to debatable ortholog identification) and helps identify gene losses or/and gene pseudogenization (e.g. opsins) using synteny arguments. On the other hand, transcriptomic characterizations of a large sampling of tissues and embryonic stages provide expression information for each individual gene family member leading to hypotheses of neo- or subfunctionalization.

## Multi-pronged Biological Datasets in the Catshark to Understand the Evolution of Sensory Systems in Chondrichthyans and in Early Jawed Vertebrates

### Visual System: Opposing Trends in Two Major Gene Families

Using these new resources, we identified and characterized a large locus in chondrichthyan genomes made of a cluster of tens of gamma-crystallin genes. In all aquatic species, many of these genes grouped together as species-specific clades in our phylogenetic reconstruction, making them recent duplicates (in an amphibian, several actinopterygians [also see Lemopoulos and Montoya-Burgos 2022], and chondrichthyans). Considering these gene clusters as the product of many recent duplications that happened independently in all sampled lineages, we argue that our data point to an ancient duplication hotspot in the genomes of jawed vertebrates. Shared by chondrichthyans and bony fishes, this locus of continuous duplication of gamma-crystallins may have been active already in early jawed vertebrates, most probably also associated with a high degeneration rate as an ancestral condition leading to the loss of any 1-to-1 orthologous gene between lineages. Consistent with this hypothesis, we identified in the catshark cluster a putative pseudogene, producing a non-coding RNA still showing extensive expression in the adult eye with a nucleotide sequence similar to other gamma-crystallin genes. We did not detect positive selection within the most recent clade of gamma-crystallins in the catshark, suggesting that maintenance of large complements of similar genes is driven by selection on the quantity of the produced proteins, rather than by positive selection on sequence variants. However, such a gene cluster with high rates of duplication (and probably pseudogenization) generates a large number of genes coding for lens proteins that may then be subject to subfunctionalization and adapt quickly to ecological constraints of vision in deep or shallow waters (Weadick and Chang 2009). Such

gene clusters may participate in population divergence, which should be further tested with adequate sequence data (e.g. eye RNA-seq data to generate SNP comparison). A reduced necessity for large amounts of these proteins may have led to the loss of most gamma-crystallin genes in terrestrial vertebrates (Plotnikova et al. 2007; Cvekl and Eliscovich 2021).

An opposing trend was observed in the opsin gene family. With this high-quality genome of the small-spotted catshark, we confirmed the secondary loss of many opsin genes in the last common ancestor of elasmobranchs (sws1, sws2, opn4x, tmt2, tmt3, and several clades of the neuropsin family), which was previously suggested on the basis of incomplete genomic data (Yamaguchi et al. 2021). In addition, we identified a previously undescribed catshark rh2 gene with restricted expression in the undifferentiated embryonic retina. The general opsin gene loss identified for the elasmobranch lineage (where rh2, one of the sws genes, opn4x, tmt, VA-opsin, and pinopsin were lost) might be compared to the losses associated with the “nocturnal bottleneck” hypothesized in mammals during their early diversification (reviewed by Gerkema et al. 2013). These co-occurring losses of opsin genes (not observed in the elephant shark [Yamaguchi et al. 2021]) should now be more densely mapped onto the chondrichthyan phylogeny to test if they could have resulted from a single event, e.g. raising the hypothesis of a deepwater habit for the common elasmobranch ancestor, or rather from gradual, progressive losses. The reduction of visual opsin genes in elasmobranchs seems to have been balanced not only by adaptive mutations in rhodopsin but also by structural adaptation of the rod cells (Bozzanao et al. 2001; Yamaguchi et al. 2023a). The catshark transcriptomic dataset further showed a strong specificity toward adult eye expression in some but not all of both visual and non-visual opsin genes, the latter being known to be expressed in the eyes of other jawed vertebrates (e.g. Upton et al. 2021; Salgado et al. 2022). Expression of opsins in other tissues raises questions on additional functions of these proteins, including at embryonic stages, as previously suggested (Liebert et al. 2022). A curious observation is the co-expression of rgr1 and encephalopsin (Opn3) in the testis. A putative thermal

<table><caption>**Table 3** Summary of the TRP genes with higher expression in sensory organs</caption><thead><tr><th></th><th></th><th>Ampullae of Lorenzini</th><th>Eye</th><th>'brain & olfactory epithelium'</th></tr></thead><tbody><tr><td rowspan="4">Group1</td><td>Trpc</td><td>...</td><td>trpc4, 5</td><td>trpc2*, trpc3, trpc5</td></tr><tr><td>Trpm</td><td>trpm G</td><td>trpm1, 3</td><td>trpm2*, trpm3</td></tr><tr><td>Trpv</td><td>trpv8.1,<br><u>trpv8.2*</u></td><td>trpv9</td><td>trpv5/G.2, trpv3,<br>trpv8.1, trpv8.2*</td></tr><tr><td>Trpw</td><td>...</td><td>...</td><td>trpw1*</td></tr><tr><td rowspan="2">Group2</td><td>Trpml</td><td>trpml3.1</td><td>...</td><td>trpml3.1</td></tr><tr><td>Trpp</td><td>pkd2</td><td>pkd1L4.1</td><td>...</td></tr></tbody></table>

Bold denotes Z-score > 4; others: Z-score between 1 and 4. Asterisk: genes validated by in situ hybridization; underlined: gene expression incongruent with sensory cell expression.
sensation of sperm cells through expression of opsin proteins, as proposed in mammals, should also be investigated in elasmobranch fishes, in particular regarding the processes of sperm cell guidance to, and storage in, the spermatheca of females (Griffiths et al. 2012).

Finally, we identified distinct combinations of Trp gene expression in several sensory organs sampled in the small-spotted catshark, including in the eye (Table 3). These results suggest considerable conservation in visual organs between the lineages of chondrichthyans and sarcopterygians, since a similar combination of Trps was reviewed for the mammalian retina, both in the photoreceptors and associated cells involved in transducing the signal (Pickard and Sollars 2011). Further comparisons with actinopterygians will be necessary to infer if this situation is conserved from the jawed vertebrate's last common ancestor with putative secondary loss in actinopterygians, or rather derived from convergent evolution.

## Olfactory System: Neofunctionalization in Three Minor and One Major Olfactory Receptor Families

We show a morphologically and molecularly well-developed olfactory epithelium already at late embryonic stages and described it in 3D as very similar between the pre-hatching and adult stages despite the size difference (Schluessel et al. 2010). We found no anatomical evidence for a separate accessory olfactory system such as those observed for both lamprey and lungfish, in this pre-hatching stage (Ren et al. 2009; Nakamuta et al. 2012). If this is representative of the adult condition, then the catshark is similar to actinopterygian fish, rather than to lungfish. This might suggest that a single olfactory epithelium is the ancestral state for vertebrates and that both accessory olfactory systems described in the lamprey and sarcopterygians (lungfish and tetrapods) have evolved independently. The anatomical arrangement of this olfactory epithelium and its projection through two olfactory nerves to the olfactory bulb is consistent with previous reports at earlier catshark stages (Quintana-Urzainqui et al. 2014) appears similar to a recent report for adult bamboo shark (Camilieri-Asch et al. 2020) and seems to be widespread in the chondrichthyan lineage (Yopak et al. 2015). Since all molecular markers investigated show similar expression in all lamellae, the division in two olfactory nerves does not appear to reflect a functional subdivision of the olfactory epithelium.

Our analysis of highly enriched expression yielded several genes with specific expression patterns in the olfactory epithelium. Trpc2 is a component of the signal transduction cascade in olfactory microvillous neurons of bony vertebrates (Liman et al. 1999; Stowers et al. 2002; Liberles 2014) and has often been used as a marker gene for these neurons (e.g. see Sato et al. 2005; Hecker et al. 2019; Syed et al. 2023). The high expression levels of *trpc2* are consistent with its presence in (nearly) all olfactory receptor neurons in the catshark (our results and Syed et al. 2023); the other main population of olfactory sensory neurons, ciliated neurons, are absent in sharks (Ferrando and Gallus 2013). The expression of *trpc2* in the olfactory epithelium suggests the

jawed vertebrate-wide conservation of Trpc2 function in chemo-sensation. The catshark gene *trpw1*, a vertebrate Trpw gene clade first described in this work, shows a very similar expression in cells of the neuronal layer, albeit less frequent than *trpc2*-expressing cells, consistent with *trpw1* expression labeling a subpopulation of microvillous neurons. In addition, the calcium-binding protein S100z has been suggested as another marker of microvillous receptor neurons in bony vertebrates (Kraemer et al. 2008; Oka et al. 2012; Hecker et al. 2019), and our results suggest it to be shared by all jawed vertebrates. Moreover, *s100z*, together with *trpw1*, expression appears to label fewer cells than *trpc2* suggesting that different subpopulations of microvillous neurons may occur in the catshark.

Moxd2 is a monoxygenase related to dopamine beta-hydroxylase that originated in vertebrates as a duplication of an older *moxd* gene. The function of Moxd proteins is unknown in any species, but its relation to dopamine-beta-hydroxylase (DBH) and the presence of *Moxd2* expression in mouse medial olfactory epithelium (array data in Su et al. 2004) have led other authors to suggest a function in neurotransmitter conversion or modification of incoming odorant molecules (Hahn et al. 2007). We note that high expression levels of *moxd2* (similar to *s100z* levels) are present also in a recent transcriptome of zebrafish olfactory epithelium (gene name *moxd1*, Calfún et al. 2016). Our transcriptome data show that enrichment of *moxd2* in the olfactory epithelium appears to be a general jawed vertebrate feature. Moreover, we identify for the first time the cell type expressing *moxd2* genes in the catshark. Both *moxd2* paralogs are abundantly expressed in the supporting (sustentacular) cells of the catshark olfactory epithelium (apical layer), but not in the neuronal nor basal cell layer. It is currently unknown whether the *moxd2* genes of bony vertebrate species share the same expression in supporting cells. However, supporting cells were shown to metabolize odorants by a two-step process involving monoxygenases in phase I (P450 cytochromes, Chen et al. 1992), and several transferases in phase II (olfactory uridine 5'-diphospho- glucuronosyl transferase, Lazard et al. 1991; and sulfotransferase, Miyawaki et al. 1996). Moxd2 may then act in the supporting cells as a phase I enzyme to begin metabolization and thereby removal of odorants from the catshark olfactory epithelium. However, *moxd2* genes were apparently lost in parallel in a variety of lineages (several chondrichthyans, our study; humans, Hahn et al. 2007; catarrhines and whales, Kim et al. 2014) suggesting functional redundancy with other monoxygenases. Conservation of two *moxd2* paralogs with similar expression patterns in the small-spotted catshark should be better understood by testing potential differences in the function of each protein.

Sharks have a small olfactory receptor repertoire compared to many bony vertebrates; e.g. the main olfactory receptor family of bony vertebrates, the ORs, is contracted in sharks and consists of only 2 OR genes plus seven OR-like genes (Syed et al. 2023). In bony vertebrates OR genes are expressed by ciliated neurons, so the loss of OR genes and
loss of OR expression in the olfactory epithelium is consist-
ent with the putative secondary loss of ciliated neurons in
sharks (Ferrando and Gallus 2013). We find that the major-
ity of OR genes are most highly expressed in other organs,
suggesting function in the immune system (spleen) or
other sensory organs (ampullae of Lorenzini). Similarly,
three even smaller families of olfactory receptors, TAARs,
TARLs (olfactory in cyclostomes only), and ORAs, are
not enriched in olfactory tissue and all of them have high-
er expression levels in different other tissues. ORAs are
supposed to be ancestral olfactory receptors, as they
are partially olfactory in cyclostomes (Kowatschew and
Korsching 2022) and olfactory in bony fish; neofunctiona-
lization outside of the olfactory system in chondrichthyans
therefore involved the loss of the ancestral function. The
catshark V2R-like genes appear to have a broad spectrum
of functions both in the olfactory epithelium and else-
where, including other sensory organs such as the eye. In
contrast, nearly all V2R genes (considered the main olfac-
tory receptor family of sharks [Syed et al. 2023]) show un-
ique enrichment in the “brain & olfactory epithelium”
tissue, consistent with other expression data (Syed et al.
2023). This supports specific function in olfaction only for
V2R genes, although additional and sometimes exclusive
expression of some V2Rs in non-olfactory tissues suggest
some degree of neofunctionalization (e.g. in ionic homeo-
stasis), even in the only family considered to be dedicated
to olfactory function in the catshark.

To summarize, the catshark (and possibly elasmobranchs
in general) appear to use primarily a single gene family for
the detection of odors (the V2Rs) out of the four canonical
vertebrate olfactory receptor gene families. This represents
a nearly inverse situation compared to the jawless vertebrate
lamprey, whose V2Rs are not expressed in the olfactory epi-
thelium, in contrast to ORs, TAAR-like genes, and V1Rs
(Kowatschew and Korsching 2022), highlighting the need
for a focused examination of the evolution of their function
in chondrichthyans. The expression patterns of several genes
in the olfactory neuron layer are suggestive of subpopulations
in these sensory cells, and axons coming from medial or lat-
eral lamellae project through 2 separate olfactory nerves:
both these observations support a structural complexity
that should be uncovered by further functional studies of
the shark olfactory system.

Electrosensory System: Catshark Provides Insight into
Ancestral State of Electroreception

Embryonic stages of the small-spotted catshark already ex-
hibit adult features in both the 3D spatial organization of
the ampullae of Lorenzini as well as in the histological dif-
ferentiation of its multiple cell types (i.e. sensory, support-
ing, and canal zones). Nervous connections between the
different clusters of ampullae and the anterior lateral line
nerves appear to coalesce with sensory neurons connected
to the sensory neuromasts of the lateral line canals.
Previous analyses of genes involved in electrosensory sys-
tems in vertebrates have focused on electrophysiology
(e.g. in elasmobranch fishes, Bellono et al. 2017) or

differential gene expression between the skin and ampullary
organ tissue (e.g. in elasmobranch fishes, Bellono et al. 2017;
and in bony fish, Modrell et al. 2017). Our analysis of tran-
scriptomic data supports the identification of glutamatergic
and ribbon-synapse components in electrosensory cells of
the catshark, similar to what was proposed for the paddlefish
(Modrell et al. 2017), with expression of glutamate transpor-
ters, L-type voltage-gated calcium channels (Ca,1 channels)
and the transmembrane protein Otoferlin (Modrell et al.
2017; Manchanda et al. 2021). In the paddlefish, Modrell
et al. (2017) reported the expression in ampullary organs
but not neuromasts of two voltage-gated potassium channel
subunit genes, *kcna5* and *kcnb3*, and one transcription fac-
tor, *mafa*. We detected high and specific expression of the
orthologous potassium voltage-gated channel *kcnb3* gene
in the catshark ampullae of Lorenzini, together with high
and specific expression of *mafa*, suggesting both the tran-
scriptional regulators and the terminal differentiation genes
of electrosensory cells show conserved ancestral characteris-
tic for all jawed vertebrates.

Our results also suggest that some TRPs might be ex-
pressed in sensory organs but with an ion-transport func-
tion instead of a sensory function (e.g. *trpv8.2* in the
ampullae of Lorenzini and the olfactory epithelium). The
ampullae of Lorenzini display their own combination of ex-
pressed TRPs (Table 3), making the members of this gene
family excellent candidates for further description of signal
transduction in the electrosensory system of vertebrates.

Conclusion

Despite the paramount importance of chondrichthyans in
vertebrate evolution, their poor suitability to genetic ap-
proaches due to their large size and long generation times
have been a major limitation in their development as bio-
logical models. However, the availability of an increasing
number of chromosome-level genomes in this clade marks
its entry in the post-genomic era, with novel perspectives
in multiple thematic fields, including not only comparative
genomics but also population genetics, physiology, neurosci-
ence, and evolutionary developmental biology. Combinations
of epigenomics, single-cell, and spatial transcriptomics are
currently revolutionizing our understanding of the cellular di-
versity of embryonic and adult organs. These techniques are
fully applicable to chondrichthyans, with large organ size
even becoming advantageous in spatially resolved ap-
proaches. Due to its experimental accessibility and key phylo-
genetic position, the small-spotted catshark is an excellent
model to understand the origin and development of cell
type and organ complexity during vertebrate evolution.

**Materials and Methods**

Sampling and Sequencing
RNA-seq Sampling and Sequencing

RNAs were directly extracted from tissues flash-frozen in
liquid nitrogen at the time of their dissection, following
the Sigma TRI-reagent® extraction protocol: one whole embryo was used for each standard embryonic stage 12, 22, 24, 26, 30, and 31 (Ballard et al. 1993), and a list of twenty-five different adult tissues were sampled from three individuals (one female, two males), all of these being of Mediterranean origin, fished off the coast of Banyuls-sur-mer, Fr. RNA libraries were sequenced on Illumina NovaSeq 6000 (150bp, paired-end). Raw data were deposited in NCBI (Bioproject: PRJEB35946; BioSample accessions: 21694013 to 21694062; see supplementary table S5, Supplementary Material online).

## Mediterranean Genome Sequencing and Annotation

Tissues used for DNA extraction were sampled from the male individuals in the sample list (supplementary table S5, Supplementary Material online). The assembly sScyCan1.1 is based on 63× PacBio data (CLR sequencing), 43× 10× Genomics Chromium data, and BioNano data generated at the Wellcome Sanger Institute, and 17× Hi-C data generated at the Baylor College of Medicine. The assembly process included the following sequence of steps: initial PacBio assembly generation with Falcon-unzip (Chin et al. 2016); retained haplotig identification with purge_dups; 10× based scaffolding with scaff10x; BioNano hybrid-scaffolding; Hi-C-guided scaffolding by the Aiden lab as part of a collaboration with the DNA Zoo Consortium (www.dnazoo.org) using 3D-DNA and Juicebox Assembly Tools (Durand et al. 2016; Dudchenko et al. 2017, 2018); Arrow polishing; and two rounds of FreeBayes polishing (Garrison and Marth 2012). Finally, the assembly was analyzed and manually improved using gEVAL (Chow et al. 2016; Howe et al. 2021), with BlobToolKit (Challis et al. 2020) used to detect possible sequence contamination. Chromosome-scale scaffolds confirmed by the Hi-C data were named in order of size. The interactive contact maps for the Hi-C-guided genome assembly step can be found at https://www.dnazoo.org/assemblies/scyliorhinus_canicula.

## Atlantic Genome Sequencing and Assembly

For Illumina libraries sequenced at Genoscope, 250 ng of genomic DNA extracted from an adult individual fished in Roscoff (France; sex unknown) were sonicated using the S2 Covaris instrument (Covaris, Inc., USA). Sheared DNA was used for Illumina library preparation by a semi-automatized protocol. Briefly, end repair, A tailing and Illumina compatible adaptors (BiooScientific) ligation were performed using the SPRIWorks Library Preparation System and SPRI TE instrument (Beckmann Coulter), according to the manufacturer protocol. A 300 to 600 bp size selection was applied in order to recover the most of fragments. DNA fragments were amplified by 12 cycles of the polymerase chain reaction (PCR) using Platinum Pfx Taq Polymerase Kit (Life Technologies) and Illumina adapter-specific primers. Librarie was purified with 0.8× AMPure XP beads (Beckmann Coulter), and size-selected on 1.5% agarose gel around 600 bp. After library profile analysis by Agilent 2100 Bioanalyzer (Agilent Technologies,

USA) and qPCR quantification (MxPro, Agilent Technologies, USA), the library was sequenced using 101 base-length read chemistry in a paired-end flow cell on the Illumina sequencer (Illumina, USA).

The mate-pair libraries were prepared using the Nextera Mate Pair Sample Preparation Kit (Illumina, San Diego, CA). Briefly, genomic DNA (4 µg) was simultaneously enzymatically fragmented and tagged with a biotinylated adaptor. Tagmented fragments were size-selected (3 to 5; 5 to 8 and 8 to 11 Kb) through regular gel electrophoresis and circularized overnight with a ligase. Linear, non-circularized fragments were digested, and circularized DNA was fragmented to 300- to 1000-bp size range using Covaris E210. Biotinylated DNA was immobilized on streptavidin beads, end-repaired, then 3'-adenylated, and Illumina adapters were added. DNA fragments were PCR-amplified using Illumina adapter-specific primers (15 cycles) and then purified. Finally, libraries were quantified by qPCR, and library profiles were evaluated using an Agilent 2100 bioanalyzer (Agilent Technologies, USA). Each library was sequenced using 100 base-length read chemistry on a paired-end flow cell on the Illumina HiSeq2000 (Illumina, USA). Raw data were deposited in European Nucleotide Archive with project reference PRJEB80353, and are available online at https://www.ebi.ac.uk/ena/browser/view/PRJEB80353.

The reads were submitted to an initial cleaning step to trim sequences and remove low-quality or very short reads and then corrected using Musket (Liu et al. 2013) with a 17 kmer size. 180 base pair (bp) paired-ends reads were joined using fastq-join (https://code.google.com/p/ea-utils/wiki/FastqJoin) with default parameters. Assembly and scaffolding of paired-end sequences was conducted using CLC (http://www.clcbio.com/products/clc-assembly-cell/) with a 61 bp kmer size and only solving bubbles shorter than 500 bp. Mate pair sequences were introduced in additional scaffolding and gap closing steps, respectively using SSPACE (Boetzer et al. 2011) and Gapcloser (Luo et al. 2012).

## Genome Structure and Synteny
### Species Tree Inference

In order to compare evolution of genome structure across chondrichthyans, we obtained thirteen genome assemblies from this group labeled as chromosome-level assemblies from NCBI. Three representative outgroup species with chromosome-level assemblies were also included in the analysis (supplementary table S6, Supplementary Material online). Species tree reconstruction of all 16 species was carried out using a data set of 1,005 genes annotated with BUSCO v5.1.2 (Manni et al. 2021) using the Vertebrata gene set. BUSCO genes present in each of the species in our dataset were extracted, aligned, and trimmed using the busco2phylo pipeline (https://github.com/istevens17/busco2phylo-nf), which employed MAFFT v7.4 (Katoh 2005) and trimAl v1.4 (Capella-Gutiérrez et al. 2009) for gene alignment and trimming. All trimmed alignments were concatenated into a supermatrix using PhyKIT (Steenwyk et al. 2021), and species
tree inference was carried out with IQ-TREE v2.0 (Minh et al. 2020), applying the LG model with a gamma distribution with 4 categories. Species trees were plotted, alongside the inferred genome size and chromosome number per species (inferred from the genome assemblies) using Toytree (Eaton 2020).

## GC Variation Across Genomes

GC content was calculated at three levels; (i) genome-wide, measured using the complete genome assembly, (ii) coding regions, measured using the nucleotide versions of the trimmed BUSCO gene alignments, and (iii) third codon position, also inferred using the BUSCO gene alignments. These three values were plotted against the species tree in supplementary fig. 34, Supplementary Material online using ggtree (Yu et al. 2017). For each species, GC content of each chromosome above 1Mb in size was plotted against chromosome size, to inspect for patterns of correlation between GC content and chromosome size in Chondrichthyes (supplementary fig. 35, Supplementary Material online).

## Genome Macrosynteny

Macrosynteny of the focal species *Scyliorhinus canicula* was compared to two other species in Chondrichthyes (*Carcharodon carcharias* and *Pristis pectinata*) using the MCscan pipeline in the jcvi software: github.com/tanghaibao/jcvi (Tang et al. 2008). Briefly, this pipeline uses BLAST to carry out a pairwise synteny search between species using the coding sequences and gene locations. Conserved synteny blocks containing 30 genes or more were retained and visualized, with *Scyliorhinus canicula* set as the focal reference species.

## Repeat Annotation

The sequences of transposable elements (TEs) in the small-spotted catshark genome were annotated with the Earl Grey TE annotation pipeline (version 1.2, github.com/TobyBaril/EarlGrey) as described previously (Baril and Hayward 2022). In brief, the pipeline was used firstly to identify known TEs from the vertebrate subset of Dfam (release 3.4) and RepBase (release 20,181,026) (Jurka et al. 2005; Hubley et al. 2016). Next, de novo TEs were identified and the boundaries of consensus sequences were extended using an automated “BLAST, Extract, Extend” process (Platt et al. 2016). Any redundant sequences were removed from the consensus library before the genome assembly was annotated using a final combined library, which consisted of the known and de novo TE libraries. TE annotations were subsequently processed to remove any overlaps and also to defragment broken annotations before the final TE quantification was performed.

## Hox Cluster Analysis

Hox protein sequences from a range of chondrichthyans (supplementary table S7, Supplementary Material online) and humans were downloaded from NCBI, and used as queries in tBLASTn searches of the *Scyliorhinus canicula* gene model dataset (Mayeur et al. 2021) and the genome assembly.

We also used chondrichthyan and human protein sequences of the HoxC flanking genes *tspan31*, *ddx23*, *rnd1*, *adcyG*, *smug1*, and *copz1* (Yamaguchi et al. 2023b) to identify the syntenic region in the *S. canicula* assembly.

## Mediterranean and Atlantic Population Differentiation

To estimate the population differentiation between Mediterranean and Atlantic catshark populations, we used RNA-seq data from three NCBI bioprojects: PRJEB36280 (Mediterranean), PRJNA255185 (North Wales, UK), and PRJNA504730 (English channel, Brittany, France). We used the Nextflow (v.23.10.0; Di Tommaso et al. 2017) f-core rnavar (v.1.1.0; Ewels et al. 2020) pipeline for read mapping and SNP calling, which is based on the Genome Analysis Toolkit (GATK) RNA-seq variant calling best-practice workflow (github.com/gatk-workflows/gatk4-rnaseq-germline-snpsnindels). We used the “--star_tropass” option to run STAR two-pass read mapping (Dobin et al. 2013) and ran read sets of different lengths separately so that optimal STAR indices could be produced for each read length (additionally, the following non-default options were used: --remove_duplicates, --skip_baserecalibration, -skip_variantannotation, and --generate_gvcf). Variant calls from all individuals were then combined and jointly genotyped using the COMBINEGVCFs and GENOTYPEGVCF programs from the GATK software suite (v.4.2.6.1; McKenna et al. 2010). To produce unbiased estimates of population genetic statistics, we called genotypes at both variant and invariant sites. All sites were filtered to a maximum missing data percentage of 15% (i.e. two individuals missing) and a minimum read depth of 5 using VCFtools. Variant and invariant sites were then separated and variant sites were additionally filtered to include only SNPs with minor allele frequency >0.1 and minimum genotype quality of 30 with VCFtools, before the sites were recombined with the concat function from BCFtools (Li 2011). Fixation index (Weir and Cockerham's *F*<sub>ST</sub>), absolute divergence (*d*<sub>XY</sub>), and nucleotide diversity for each population (π) were calculated in non-overlapping 1Mb windows with pixy (Korunes and Samuk 2021). Additionally, per-site Weir and Cockerham's *F*<sub>ST</sub> was calculated using VCFtools. Gene density and GC content were calculated in non-overlapping 1Mb windows using BEDTools (Quinlan and Hall 2010). Outlier SNPs potentially under divergent selection between Atlantic and Mediterranean populations were identified using pcadapt (Luu et al. 2017) with *K* = 2, LD clumping (size = 500, thresh = 0.1) and the “componentwise” method. Genome-wide statistics were visualized with qqman (Turner 2018). Principal Component Analysis (PCA) was performed using PLINK 2.0 (v.2.00a3LM; Chang et al. 2015), after thinning SNPs to one per 100 Kb to reduce the impact of genetic linkage. Population structure was estimated using ADMIXTURE (v.1.3.0; Alexander et al. 2009) using the same thinned set of SNPs as the PCA analysis and 5-fold cross-validation. The code used in this section is available at github.com/ogosborne/Scyliorhinus_canicula_popgen.
## Reference Gene Model

The reference gene model “ncbi-utrs” previously published in Mayeur et al. (2021) was used for RNA-seq mapping. Its construction aimed at maximizing sequence predictions in 3' UTRs, and it was conducted as follows using custom R and PERL scripts (supplementary text, Supplementary Material online). First, an “isoform-collapsed” version of gene models was constructed by selecting the most supported isoform of each gene, as defined in NCBI Annotation Release 100 of the *Scyliorhinus canicula* genome (sScyCan1.1 version). All NCBI predicted transcripts were then mapped back by BLASTn (identity = 100% and e-value < 1E-06) to this ‘isoform-collapsed’ reference in order to identify for each gene the isoform with the longest 3’UTR. When detected, the additional 3’UTR sequence was appended to the sequence of the “isoform-collapsed” reference. A second round of 3’UTR expansion was conducted from the resulting gene model set by including divergent isoforms identified with lower BLASTn supports (identity < 100% and e-value < 1E-06). A final iteration of 3’UTR expansion was carried out using a transcriptome assembled from available public NCBI SRA reads as well as locally obtained RNA-seq data as a query (listed in supplementary table S8, Supplementary Material online). The transcriptome assembly used in this final step was conducted using a de novo RNA-Seq assembly pipeline (Cabau et al. 2017) and isoform-collapsed using Corset and SuperTranscripts (Davidson and Oshlack 2014; Davidson et al. 2017). In this final step, the 3’UTR appending step was carried out using only BLASTn hits with 100% identity with the subject, e-value < 1E-06.

## Transcriptome Mapping and Gene Expression Ranking

The paired-end RNA-seq reads from adult and embryonic *S. canicula* tissue samples were mapped onto the reference gene model ncbi-utrs using BWA-MEM (Li 2013). After sorting and converting the resulting alignments to BAM files using SAMtools, counts for each library were extracted. Count data was then normalized as Transcripts Per Kilobase Millions (TPMs) and used as such for quantification of expression (supplementary dataset S9, Supplementary Material online). If a tissue had several replicates (max 3 replicates), TPM values were averaged, leading to a final number of 31 tissues compared. For each gene, Z-scores were then calculated to evaluate expression over-represented in a tissue, as $Z = (x - \mu) / \sigma$ where for any sample: *x* is the observed TPM value in a given tissue, *μ* is the mean TPM value of all tissues, and *σ* is the standard deviation of all tissues. We generated 3 Z-score tables: one encompassing all tissues and stages (supplementary dataset S10, Supplementary Material online), another calculated only with TPM values in embryonic stages (supplementary dataset S1, Supplementary Material online), and one calculated only with TPM values in adult tissues (supplementary dataset S4, Supplementary Material online). For the analysis of gene expression

dynamics during the developmental window spanning stages 22 to 31, we used Moran’s index as a measure of temporal autocorrelation. Moran’s indexes and their statistical support were calculated for each gene with the R package spdep and its function moran.test following a data formatting step with ade4. To identify genes with high level and high specificity of expression in a given organ, we calculated a score integrating expression level and tissue specificity as follows: $score = \ln(TPM + 1) * Z$. For each tissue, we could then rank the genes of interest and focused on the 50 genes with the highest score.

## Phylogenetic Reconstructions

Phylogenetic analyses of gene families were conducted by retrieving protein sequences for all ohnologs (orthologs and paralogs) from a set of species including actinopterygians (*Erpetoichthys calabaricus*, *Lepisosteus oculatus*, *Danio rerio*, *Oryzias latipes*), amniotes (*Homo sapiens* or *Mus musculus*, *Gallus gallus* or *Anolis carolinensis*), an amphibian (*Xenopus tropicalis*), a lungfish (*Protopterus anneticens*) and chondrichthyans (*Scyliorhinus canicula*, *Carcharodon carcharias*, *Amblyraja radiata*, *Callorhinchus milii*), as well as cyclostomes (*Petromyzon marinus*, *Eptatretus burgeri*) when available. Homologous sequences from the amphioxus *Branchiostoma floridae* were also included as an outgroup reference. This search was conducted using Genomicus, controlled by synteny data when available (Muffato et al. 2010; Louis et al. 2013, 2015) and TBLASTN (e-value < 10⁻¹⁰) searches against Genbank or relevant predicted cDNA databases, controlled by reciprocal BLAST. The protein sequences were aligned using SeaView v4.6.4 (Gouy et al. 2010), MUSCLE v3.7 (Robert C Edgar 2004a, 2004b), and MAFFT v7 (Katoh et al. 2019). Gene tree inferences were performed using maximum likelihood with IQ-Tree (Minh et al. 2020) and ModelFinder to find the model of best fit (Kalyaanamoorthy et al. 2017). Branch supports were given as percentages of 1000 SH-aLRT branch tests and 1000 UltrafastBootstrap tests (Minh et al. 2013).

## 3D Reconstruction of Synchrotron Data

Iodine staining and synchrotron scanning of the head of a *S. canicula* hatchling (8cm total length) were performed as previously described (Leyhr et al. 2023). Briefly, the head was fixed in 4% paraformaldehyde, rinsed in 1× phosphate-buffered saline, and gradually dehydrated in 25%, 50%, 75%, and 100% ethanol. Staining with iodine was performed overnight in a 1% iodine solution in 100% ethanol according to the previously published protocol for I₂E staining (Metscher 2009). Iodine staining was combined with synchrotron phase contrast (DICE-PPC-SRµCT) to highlight the organs of interest. After staining, the shark head was transferred to 96% ethanol before being scanned at 3 µm isotropic voxel size using the Hierarchical Phase-Contrast Tomography (HiP-CT) protocol developed by Walsh et al. (2021) to optimize the contrast of the images and reduce the noise. The experiment was performed at the beamline BM05 of
the European Synchrotron Radiation Facility—Extremely Brilliant Source (ESRF-EBS) in France (proposal LS-3021). A PCO Edge CLHS 4.2 camera was used with a magnifying optic based on a Canon 68 mm F/D 2.8 Super Macro lens coupled to a 250-µm-thick LuAG:Ce scintillator to produce images with a 3-µm-voxel size. The beam was filtered with 2.3 mm of aluminum and 8 5-mm bars of silicon dioxide, thereby resulting in a detected average energy of 78 keV. To maximize the phase contrast, the specimen was placed at a distance of 1.45 m from the detector. Tomographic data were reconstructed using filtered back-projection algorithm with the software PyHST2 (Mirone et al. 2014) coupled to a modified version (Sanchez et al. 2012) of single-distance phase retrieval (Paganin et al. 2002). Reconstructed jpeg2000 image stacks were imported into VGStudio MAX version 3.5.1 (Volume Graphics, Germany) for manual segmentation and rendering. Organ volumes and surface areas were calculated using the Porosity/Inclusion Analysis module. Only the left side of the head was manually segmented, and bilateral organ volumes are reported as double the left values, assuming bilateral symmetry. The dataset used in this study for 3D virtual microanatomy/histology is available at the ESRF repository (https://doi.org/10.15151/esrf-dc-1947663452) and additional 3D visuals will be accessible in the Paleontology database (https://paleo.esrf.eu/).

## Histology and *in situ* Hybridization

The head of embryos at stage 31 (late stage of organogenesis) and stage 34 (preceding hatching) and of a recently hatched juvenile (6 wk post-hatching) were sampled from individuals previously euthanized, fixed in paraformaldehyde 4% in phosphate-buffered saline 1× (PBS), and preserved in 100% ethanol (after graded dehydration in ethanol/PBS) for other studies. These samples were embedded in paraffin and cut at 10 µm thickness. In situ hybridizations with RNA probes were performed as described previously (Lagadec et al. 2015) followed by nuclear counterstaining with Nuclear Fast Red (Sigma) and mounting in Eukitt. DNA matrices were synthesized by Twist Bioscience with probe sequences as given in supplementary table S9, Supplementary Material online. For crystallin genes, we designed two probes (against XM_038788722.1 and XM_038788824.1), both supposed to hybridize all of the recently duplicated crystallin gene copies as the nucleotide sequences differ only by about 5%. More distant copies (XM_038788959.1, XM_038790385.1, XM_038788706.1) differ from the probe sequence of about 15%. These two probes gave equal expression patterns. The probes for vomeronasal receptor type 2 genes *v2rl1* and *v2rl4* were not considered to cross-react (about 65% distance in the alignment between the probe and the potentially cross-reacting gene sequence). Similar dissimilarity was observed for the two copies of Monooxygenase DBH Like 2 genes *moxd2.1* and *moxd2.2* (60% distance one way, 75% the other way). Histological staining with Hematoxylin-Eosin-Saffron was made by an automatLeica autostainer X with standard protocol, at the

RHEM platform in Montpellier on adjacent sections of the same individual (10 µm thick). All slides were scanned with a Nanozoomer S210 (Hamamatsu).

## Supplementary Material

Supplementary material is available at *Molecular Biology and Evolution* online.

## Acknowledgments

The work was funded by the Centre National de la Recherche Scientifique, by Agence Nationale de la Recherche contract n°ANR-16-CE13-0013-02 to S.M. It was supported by a PhD fellowship to L.M. (Ministère de la Recherche, ED515). T.B. was supported by a studentship from the Biotechnology and Biological Sciences Research Council-funded South West Biosciences Doctoral Training Partnership (BB/M009122/1). A.H. was supported by a Biotechnology and Biological Sciences Research Council (BBSRC) David Phillips Fellowship (BB/N020146/1). E.L.A. was supported by grants from the Welch Foundation (Q-1866), an NIH Encyclopedia of DNA Elements Mapping Center Award (UM1HG009375), a US-Israel Binational Science Foundation Award (2019276), the Behavioral Plasticity Research Institute (NSF DBI-2021795), NSF Physics Frontiers Center Award (NSF PHY-2210291), and an NIH CEGS (RM1HG011016-01A1). Genome assembly was performed in association with the DNA Zoo Consortium (www.dnazoo.org), which acknowledges support from Illumina, IBM, and Pawsey Supercomputing Center. S.S. was awarded a grant from the Swedish research Council–Vetenskaprådet (2019-04595) to fund the workstation and software used for 3D virtual microanatomy and histology investigations.

We thank Ram Reshef, University of Haifa, for giving access to transcriptomic data prior to public release and the Sanger Institute for providing access to the catshark genome prior to public release. We thank Dr. Paul Tafforeau and Dr. Kathleen N. Dollman at the European Synchrotron Radiation Facility for their assistance in the data acquisition and tomographic reconstruction. Beamline at BM05 was granted as part of proposal LS3021 (awarded to S.S., T.H., J.L., and M.D.-T.). We thank Dr. Vincent Dupret and Vishruth Venkataraman for their assistance in delineating the catshark cranial nerves. We acknowledge ABIMS and Genotoul for access to bioinformatic facilities, and UMR7232 Service de Bio-Informatique BSBII for HM's support. We thank Michaël Fuentes and Pascal Romans and the Centre de Ressources Biologiques Marines of the Banyuls Oceanological Observatory (OOB) for providing catshark specimens, and Marion Coolen for the initial sampling of an Atlantic catshark. Catshark infrastructures of OOB and collaborations were funded by EMBRC-France and EMBRC-Europe. To generate the histology data, we acknowledge the "Réseau d'Histologie Expérimentale de Montpellier"-RHEM facility supported by SIRIC Montpellier Cancer Grant INCa_Inserm_DGOS_12553, the
European regional development foundation and the Occitanie region (FEDER-FSE 2014-2020 Languedoc Roussillon), REACT-EU (Recovery Assistance for Cohesion and the Territories of Europe), IBISA and Ligue contre le cancer for histology technics.

## Data Availability

The data underlying this article are available in the article, and in its online supplementary material, and in InDores, at https://doi.org/10.48579/PRO/B8OJ7Y.

## References

Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009:19(9): 1655-1664. https://doi.org/10.1101/gr.094052.109.

Ballard WW, Mellinger J, Lechenault H. A series of stages for development of *Scyliorhinus canicula* the lesser spotted dogfish (Chondrichthyes: Scyliorhinidae). J Exp Zool. 1993:267(3): 318-336. https://doi.org/10.1002/jez.1402670309.

Baril T, Hayward A. Migrators within migrators: exploring transposable element dynamics in the monarch butterfly, danaus plexippus. Mob DNA. 2022:13(1):5. https://doi.org/10.1186/s13100-022-00263-5.

Bear DM, Lassance J-M, Hoekstra HE, Datta SR. The evolving neural and genetic architecture of vertebrate olfaction. Curr Biol. 2016:26(20): R1039-R1049. https://doi.org/10.1016/j.cub.2016.09.011.

Bellono NW, Leitch DB, Julius D. Molecular basis of ancestral vertebrate electroreception. Nature. 2017:543(7645):391-396. https://doi.org/10.1038/nature21401.

Berio F, Bayle Y, Baum D, Goudemand N, Debiais-Thibaud M. Hide and seek shark teeth in random forests: machine learning applied to *Scyliorhinus canicula* populations. PeerJ. 2022:10:e13575. https://doi.org/10.7717/peerj.13575.

Blackshaw S, Snyder SH. Encephalopsin: a novel mammalian extra-retinal opsin discretely localized in the brain. J Neurosci. 1999:19(10):3681-3690. https://doi.org/10.1523/JNEUROSCI.19-10-03681.1999.

Boetzer M, Henkel CV, Jansen HJ, Butler D, Pirovano W. Scaffolding pre-assembled contigs using SSPACE. Bioinformatics. 2011:27(4): 578-579. https://doi.org/10.1093/bioinformatics/btq683.

Boord RL, Campbell CBG. Structural and functional organization of the lateral line system of sharks. Am Zool. 1977:17(2):431-441. https://academic.oup.com/icb/article/17/2/431/163641.

Bottaro M. Sixth sense in the deep-sea: the electrosensory system in ghost shark Chimaera monstrosa. Sci Rep. 2022:12:9848. https://doi.org/10.1038/s41598-022-14076-2.

Bozzanao A, Murgia R, Vallerga S, Hirano J, Archer S. The photo-receptor system in the retinae of two dogfishes, *Scyliorhinus canicula* and *Galeus melastomus*: possible relationship with depth distribution and predatory lifestyle. J Fish Biol. 2001:59(5): 1258-1278. https://doi.org/10.1111/j.1095-8649.2001.tb00190.x.

Burger JW, Hess WN. Function of the rectal gland in the spiny dogfish. Science. 1960:131(3401):670-671. https://doi.org/10.1126/science.131.3401.670.

Cabau C, Escudié F, Djari A, Guiguen Y, Bobe J, Klopp C. Compacting and correcting Trinity and oases RNA-Seq de novo assemblies. PeerJ. 2017:5:e2988. https://doi.org/10.7717/peerj.2988.

Calfún C, Domínguez C, Pérez-Acle T, Whitlock KE. Changes in olfactory receptor expression are correlated with odor exposure during early development in the zebrafish (*Danio rerio*). Chem Senses. 2016:41(4):301-312. https://doi.org/10.1093/chemse/bjw002.

Camilieri-Asch V, Shaw JA, Yopak KE, Chapuis L, Partridge JC, Collin SP. Volumetric analysis and morphological assessment of the ascending olfactory pathway in an elasmobranch and a teleost

using diceCT. Brain Struct Funct. 2020:225(8):2347-2375. https://doi.org/10.1007/s00429-020-02127-1.

Capella-Gutiérrez S, Silla-Martinez JM, Gabaldón T. Trimal: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics. 2009:25(15):1972-1973. https://doi.org/10.1093/bioinformatics/btp348.

Cavalier-Smith T. Nuclear volume control by nucleoskeletal DNA, selection for cell volume and cell growth rate, and the solution of the DNA C-value paradox. J Cell Sci. 1978:34(1):247-278. https://doi.org/10.1242/jcs.34.1.247.

Challis R, Richards E, Rajan J, Cochrane G, Blaxter M. BlobToolKit-interactive quality assessment of genome assemblies. G3. 2020:10(4):1361-1374. https://doi.org/10.1534/g3.119.400908.

Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience. 2015:4(1):7. https://doi.org/10.1186/s13742-015-0047-8.

Chen Y, Getchell ML, Ding X, Getchell T V. Immunolocalization of two cytochrome P450 isozymes in rat nasal chemosensory tissue. Neuroreport. 1992:3(9):749-752. https://doi.org/10.1097/0001756-199209000-00007.

Chin C-S, Peluso P, Sedlazeck FJ, Nattestad M, Concepcion GT, Clum A, Dunn C, O'Malley R, Figueroa-Balderas R, Morales-Cruz A, et al. Phased diploid genome assembly with single-molecule real-time sequencing. Nat Methods. 2016:13(12):1050-1054. https://doi.org/10.1038/nmeth.4035.

Chow W, Brugger K, Caccamo M, Sealy I, Torrance J, Howe K. gEVAL-a web-based browser for evaluating genome assemblies. Bioinformatics. 2016:32(16):2508-2510. https://doi.org/10.1093/bioinformatics/btw159.

Clark A, Porter M, Meredith T. Morphometric analysis of the elasmobranch olfactory rosette. J Morphol. 2022:283(11):1464-1477. https://doi.org/10.1002/jmor.21514.

Collin SP. The neuroecology of cartilaginous fishes: sensory strategies for survival. Brain Behav Evol. 2012:80:80-96. https://doi.org/10.1159/000339870.

Cook JD, Ng SY, Lloyd M, Eddington S, Sun H, Nathans J, Bok D, Radu RA, Travis GH. Peropsin modulates transit of vitamin A from retina to retinal pigment epithelium. J Biol Chem. 2017:292(52): 21407-21416. https://doi.org/10.1074/jbc.M117.812701.

Coolen M, Menuet A, Chassoux D, Compagnucci C, Henry S, Lèvèque L, Da Silva C, Gavory F, Samain S, Wincker P, et al. The dogfish *Scyliorhinus canicula*: a reference in jawed vertebrates. CSH Protoc. 2008:2008:pdb.emo111. https://doi.org/10.1101/pdb.emo111.

Coolen M, Sauka-Spengler T, Nicolle D, Le-Mentec C, Lallemand Y, Da Silva C, Plouhinec J-L, Robert B, Wincker P, Shi D-L, et al. Evolution of axis specification mechanisms in jawed vertebrates: insights from a chondrichthyan. PLoS One. 2007:2(4):e374. https://doi.org/10.1371/journal.pone.0000374.

Crooks N, Waring CP. A study into the sexual dimorphisms of the Ampullae of Lorenzini in the lesser-spotted catshark, *Scyliorhinus canicula* (Linnaeus, 1758). Environ Biol Fishes. 2013:96(5): 585-590. https://doi.org/10.1007/s10641-012-0048-8.

Cvekl A, Eliscovich C. Crystallin gene expression: insights from studies of transcriptional bursting. Exp Eye Res. 2021:207:108564. https://doi.org/10.1016/j.exer.2021.108564.

Cvekl A, Zhao Y, McGreal R, Xie Q, Gu X, Zheng D. Evolutionary origins of Pax6 control of crystallin genes. Genome Biol Evol. 2017:9(8):2075-2092. https://doi.org/10.1093/gbe/exv153.

Davidson NM, Hawkins ADK, Oshlack A. SuperTranscripts: a data driven reference for analysis and visualisation of transcriptomes. Genome Biol. 2017:18(1):148. https://doi.org/10.1186/s13059-017-1284-1.

Davidson NM, Oshlack A. Corset: enabling differential gene expression analysis for de novoassembled transcriptomes. Genome Biol. 2014:15(7):410. https://doi.org/10.1186/s13059-014-0410-6.

Davies WIL, Sghari S, Upton BA, Nord C, Hahn M, Ahlgren U, Lang RA, Gunhaga L. Distinct opsin 3 (Opn3) expression in the developing nervous system during mammalian embryogenesis.
eNeuro. 2021;8(5):ENEURO.0141-21.2021. https://doi.org/10.1523/ENEURO.0141-21.2021.

Debiais-Thibaud M, Oulion S, Bourrat F, Laurenti P, Casane D, Borday-Birraux V. The homology of odontodes in gnathostomes: insights from Dlx gene expression in the dogfish, *Scyliorhinus canicula*. BMC Evol Biol. 2011;11(1):307. https://doi.org/10.1186/1471-2148-11-307.

Debiais-Thibaud M, Simon P, Ventéo S, Muñoz D, Marcellini S, Mazan S, Haitina T. Skeletal mineralization in association with type X collagen expression is an ancestral feature for jawed vertebrates. Mol Biol Evol. 2019;36(10):2265-2276. https://doi.org/10.1093/molbev/msz145.

Díaz NM, Morera LP, Tempesti T, Guido ME. The visual cycle in the inner retina of chicken and the involvement of retinal G-protein-coupled receptor (RGR). Mol Neurobiol. 2017;54(4):2507-2517. https://doi.org/10.1007/s12035-016-9830-5.

Dieris M, Kowatschew D, Korsching SI. Olfactory function in the trace amine-associated receptor family (TAARs) evolved twice independently. Sci Rep. 2021;11(1):7807. https://doi.org/10.1038/s41598-021-87236-5.

Di Tommaso P, Chatzou M, Floden EW, Barja PP, Palumbo E, Notredame C. Nextflow enables reproducible computational workflows. Nat Biotechnol. 2017;35(4):316-319. https://doi.org/10.1038/nbt.3820.

Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15-21. https://doi.org/10.1093/bioinformatics/bts635.

Dudchenko O, Batra SS, Omer AD, Nyquist SK, Hoeger M, Durand NC, Shamim MS, Machol I, Lander ES, Aiden AP, et al. De novo assembly of the *Aedes aegypti* genome using Hi-C yields chromosome-length scaffolds. Science. 2017;356(6333):92-95. https://doi.org/10.1126/science.aal3327.

Dudchenko O, Shamim MS, Batra SS, Durand NC, Musial NT, Mostofa R, Pham M, Hilaire BGS, Yao W, Stamenova E, et al. The Juicebox Assembly Tools module facilitates de novo assembly of mammalian genomes with chromosome-length scaffolds for under $1000. bioRxiv 254797. https://doi.org/10.1101/254797, 28 January 2018, preprint: not peer reviewed.

Durand NC, Robinson JT, Shamim MS, Machol I, Mesirov JP, Lander ES, Aiden EL. Juicebox provides a visualization system for Hi-C contact maps with unlimited zoom. Cell Syst. 2016;3(1):99-101. https://doi.org/10.1016/j.cels.2015.07.012.

Dymek J, Muñoz P, Mayo-Hernández E, Kuciel M, Żuwala K. Comparative analysis of the olfactory organs in selected species of marine sharks and freshwater batoids. Zool Anz. 2021;294:50-61. https://doi.org/10.1016/j.jcz.2021.07.013.

Eaton DAR. Toytree: a minimalist tree visualization and manipulation library for Python. Methods Ecol Evol. 2020;11(1):187-191. https://doi.org/10.1111/2041-210X.13313.

Ebert DA, Dando M, Fowler S. Sharks of the world. Princeton University Press; 2021.

Edgar RC. MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics. 2004a:5(1):113. https://doi.org/10.1186/1471-2105-5-113.

Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004b:32(5):1792-1797. https://doi.org/10.1093/nar/gkh340.

England SJ, Campbell PC, Banerjee S, Swanson AJ, Lewis KE. Identification and expression analysis of the complete family of zebrafish pkd genes. Front Cell Dev Biol. 2017;5:5. https://doi.org/10.3389/fcell.2017.00005.

Ewels PA, Peltzer A, Fillinger S, Patel H, Alneberg J, Wilm A, Garcia MU, Di Tommaso P, Nahnsen S. The nf-core framework for community-curated bioinformatics pipelines. Nat Biotechnol. 2020;38(3):276-278. https://doi.org/10.1038/s41587-020-0439-x.

Ferrando S, Gallus L. Is the olfactory system of cartilaginous fishes a vomeronasal system? Front Neuroanat. 2013;7:37. https://doi.org/10.3389/fnana.2013.00037.

Ferrando S, Gallus L, Gambardella C, Ghigliotti L, Ravera S, Vallarino M, Vacchi M, Tagliafierro G. Cell proliferation and apoptosis in the olfactory epithelium of the shark Scyliorhinus canicula. J Chem Neuroanat. 2010;40(4):293-300. https://doi.org/10.1016/j.chemneu.2010.08.004.

Ferrando S, Gallus L, Ghigliotti L, Amaroli A, Abbas G, Vacchi M. Clarification of the terminology of the olfactory lamellae in Chondrichthyes. Anat Rec. 2017;300(11):2039-2045. https://doi.org/10.1002/ar.23632.

Flores-Aldama L, Vandewege MW, Zavala K, Colenso CK, Gonzalez W, Brauchi SE, Opazo JC. Evolutionary analyses reveal independent origins of gene repertoires and structural motifs associated to fast inactivation in calcium-selective TRPV channels. Sci Rep. 2020;10(1):8684. https://doi.org/10.1038/s41598-020-65679-6.

Fung HF, Bergmann DC. Function follows form: how cell size is harmed for developmental decisions. Eur J Cell Biol. 2023;102:151312. https://doi.org/10.1016/j.ejcb.2023.151312.

Gans C, Northcutt RG. Neural crest and the origin of vertebrates: a new head. Science. 1983;220(4594):268-273. https://doi.org/10.1126/science.220.4594.268.

Garrison E, Marth G. 2012. Haplotype-based variant detection from short-read sequencing, arXiv, arXiv:1207.3907, preprint: not peer reviewed.

Gees M, Colsoul B, Nilius B. The role of transient receptor potential cation channels in Ca2+ signaling. Cold Spring Harb Perspect Biol. 2010;2(10):a003962. https://doi.org/10.1101/cshperspect.a003962.

Gerkema MP, Davies WIL, Foster RG, Menaker M, Hut RA. The nocturnal bottleneck and the evolution of activity patterns in mammals. Proc Biol Sci. 2013;280(1765):20130508. https://doi.org/10.1098/rspb.2013.0508.

Ghahghaei A, Rekas A, Carver JA, Augusteyn RC. Structure/function studies of dogfish alpha-crystallin, comparison with bovine alpha-crystallin. Mol Vis. 2009;15:2411-2420. https://doi.org/molvis/v15/a257.

Gillis JA, Bennett S, Criswell KE, Rees J, Sleight VA, Hirschberger C, Calzarette D, Kerr S, Dasen J. Big insight from the little skate: *Leucoraja erinacea* as a developmental model system. In: Goldstein B, Srivastava M, editors. Current topics in developmental biology. London: Academic Press; 2022. Chapter 21, Vol. 147, p. 595-630.

Gouy M, Guindon S, Gascuel O. SeaView version 4: a multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Mol Biol Evol. 2010;27(2):221-224. https://doi.org/10.1093/molbev/msp259.

Griffiths AM, Jacoby DMP, Casane D, McHugh M, Croft DP, Genner MJ, Sims DW. First analysis of multiple paternity in an oviparous shark, the small-spotted catshark (*Scyliorhinus canicula* L.). J Hered. 2012;103(2):166-173. https://doi.org/10.1093/jhered/esr112.

Grus WE, Zhang J. Origin of the genetic components of the vomeronasal system in the common ancestor of all extant vertebrates. Mol Biol Evol. 2009;26(2):407-419. https://doi.org/10.1093/molbev/msn262.

Gubili C, Sims DW, Veríssimo A, Domenici P, Ellis J, Grigoriou P, Johnson AF, McHugh M, Neat F, Satta A, et al. A tale of two seas: contrasting patterns of population structure in the small-spotted catshark across Europe. R Soc Open Sci. 2014;1(3):140175. https://doi.org/10.1098/rsos.140175.

Hahn J, Jeong S, Lee B. Inactivation of MOXD2 and S100A15A by exon deletion during human evolution. Mol Biol Evol. 2007;24(10):2203-2212. https://doi.org/10.1093/molbev/msm146.

Hara Y, Yamaguchi K, Onimaru K, Kadota M, Koyanagi M, Keeley SD, Tatsumi K, Tanaka K, Motone F, Kageyama Y, et al. Shark genomes provide insights into elasmobranch evolution and the origin of vertebrates. Nat Ecol Evol. 2018;2(11):1761-1771. https://doi.org/10.1038/s41559-018-0673-5.

Hart NS, Lamb TD, Patel HR, Chuah A, Natoli RC, Hudson NJ, Cutmore SC, Davies WIL, Collin SP, Hunt DM. Visual opsin diversity in sharks and rays. Mol Biol Evol. 2020;37:811-827. https://doi.org/10.1093/molbev/msz269.
Hashimshony T, Senderovich N, Avital G, Klochendler A, de Leeuw Y, Anavy L, Gennert D, Li S, Livak KJ, Rozenblatt-Rosen O, et al. CEL-Seq2: sensitive highly-multiplexed single-cell RNA-Seq. Genome Biol. 2016;17(1):77. https://doi.org/10.1186/s13059-016-0938-8.

Hecker N, Lächelle U, Stuckas H, Giere P, Hiller M. Convergent vomeronasal system reduction in mammals coincides with convergent losses of calcium signalling and odorant-degrading genes. Mol Ecol. 2019;28(16):3656-3668. https://doi.org/10.1111/mec.15180.

Himmel NJ, Cox DN. Transient receptor potential channels: current perspectives on evolution, structure, function and nomenclature. Proc Biol Sci. 2020;287(1933):20201309. https://doi.org/10.1098/rspb.2020.1309.

Howe K, Chow W, Collins J, Pelan S, Pointon D-L, Sims Y, Torrance J, Tracey A, Wood J. Significantly improving the quality of genome assemblies through curation. Gigascience. 2021;10(1):giaa153. https://doi.org/10.1093/gigascience/giaa153.

Hubley R, Finn RD, Clements J, Eddy SR, Jones TA, Bao W, Smit AFA, Wheeler TJ. The Dfam database of repetitive DNA families. Nucleic Acids Res. 2016;44(D1):D81-D89. https://doi.org/10.1093/nar/gkv1272.

Huffer KE, Aleksandrova AA, Jara-Oseguera A, Forrest LR, Swartz KJ. Global alignment and assessment of TRP channel transmembrane domain structures to explore functional mechanisms. Elife. 2020;9:e58660. https://doi.org/10.7554/eLife.58660.

Janvier P. Early vertebrates. New York: Oxford University Press Inc; 1996.

Jenkins A, Muñoz M, Tarttelin EE, Bellingham J, Foster RG, Hankins MW. VA opsin, melanopsin, and an inherent light response within retinal interneurons. Curr Biol. 2003;13(15):1269-1278. https://doi.org/10.1016/S0960-9822(03)00509-8.

Jurka J, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Wallichiewicz J. Repbase update, a database of eukaryotic repetitive elements. Cytogenet Genome Res. 2005;110(1-4):462-467. https://doi.org/10.1159/000084979.

Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods. 2017;14(6):587-589. https://doi.org/10.1038/nmeth.4285.

Kappé G, Purkiss AG, van Genesen ST, Slingsby C, Lubsen NH. Explosive expansion of βγ-crystallin genes in the ancestral vertebrate. J Mol Evol. 2010;71(3):219-230. https://doi.org/10.1007/s00239-010-9379-2.

Kato M, Sugiyama T, Sakai K, Yamashita T, Fujita H, Sato K, Tomonari S, Shichida Y, Ohuchi H. Two opsin 3-related proteins in the chicken retina and brain: a TMT-type opsin 3 is a blue-light sensor in retinal horizontal cells, hypothalamus, and cerebellum. PLoS One. 2016;11(11):e0163925. https://doi.org/10.1371/journal.pone.0163925.

Katoh K. MAFFT version 5: improvement in accuracy of multiple sequence alignment. Nucleic Acids Res. 2005;33(2):511-518. https://doi.org/10.1093/nar/gki198.

Katoh K, Rozewicki J, Yamada KD. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 2019;20(4):1160-1166. https://doi.org/10.1093/bib/bbx108.

Katona G, Szabó F, Végvári Z, Székely T, Liker A, Freckleton RP, Vági B, Székely T. Evolution of reproductive modes in sharks and rays. J Evol Biol. 2023;36(11):1630-1640. https://doi.org/10.1111/jeb.14231.

Kim DS, Wang Y, Oh HJ, Lee K, Hahn Y. Frequent loss and alteration of the MOXD2 gene in catarrhines and whales: a possible connection with the evolution of olfaction. PLoS One. 2014;9(8):e104085. https://doi.org/10.1371/journal.pone.0104085.

Kitamura A, Tsurugizawa T, Uematsu A, Uneyama H. The sense of taste in the upper gastrointestinal tract. Curr Pharm Des. 2014;20(16):2713-2724. https://doi.org/10.2174/13816128113199990569.

Korsching SL. Olfaction. In: Currie S, Evans DH, editors. The physiology of fishes. 5th ed. Boca Raton (FL): CRC Press; 2020. p. 191-202.

Korunes KL, Samuk K. PIXY: unbiased estimation of nucleotide diversity and divergence in the presence of missing data. Mol Ecol Resour. 2021;21(4):1359-1368. https://doi.org/10.1111/1755-0998.13326.

Kousteni V, Kasapidis P, Kotoulas G, Megalofonou P. Strong population genetic structure and contrasting demographic histories for the small-spotted catshark (*Sclyiorhinus canicula*) in the Mediterranean sea. Heredity (Edinb). 2015;114(3):333-343. https://doi.org/10.1038/hdy.2014.107.

Kowatschew D, Korsching SL. An ancient adenosine receptor gains olfactory function in bony vertebrates. Genome Biol Evol. 2021;13(9):evab211. https://doi.org/10.1093/gbe/evab211.

Kowatschew D, Korsching SL. Lamprey possess both V1R and V2R olfactory receptors, but only V1Rs are expressed in olfactory sensory neurons. Chem Senses. 2022;47:bjac007. https://doi.org/10.1093/chemse/bjac007.

Koyanagi M, Saito T, Wada S, Nagata T, Kawano-Yamashita E, Terakita A. Optogenetic potentials of diverse animal opsin: parapinopsin, peropsin, LWS bistable opsin. Adv Exp Med Biol. 2021;1293:141-151. https://doi.org/10.1007/978-981-15-8763-4_8.

Kraemer AM, Saraiva LR, Korsching SL. Structural and functional diversification in the teleost S100 family of calcium-binding proteins. BMC Evol Biol. 2008;8(1):48. https://doi.org/10.1186/1471-2148-8-48.

Kumar S, Suleski M, Craig JM, Kasprowicz AE, Sanderford M, Li M, Stecher G, Hedges SB. TimeTree 5: an expanded resource for species divergence times. Mol Biol Evol. 2022;39(8):msac174. https://doi.org/10.1093/molbev/msac174.

Kuraku S. Shark and ray genomics for disentangling their morphological diversity and vertebrate evolution. Dev Biol. 2021;477:262-272. https://doi.org/10.1016/j.ydbio.2021.06.001.

Lagadec R, Laguerre L, Menuet A, Amara A, Rocancourt C, Péricard P, Godard BG, Celina Rodicio M, Rodriguez-Moldes I, Mayeur H, et al. The ancestral role of nodal signalling in breaking L/R symmetry in the vertebrate forebrain. Nat Commun. 2015;6(1):6686. https://doi.org/10.1038/ncomms7686.

Laurin M, Canoville A, Struble M, Organ C, de Buffrénil V. Early genome size increase in urodeles. C R Palevol. 2016;15(1-2):74-82. https://doi.org/10.1016/j.crpv.2014.12.006.

Lazard D, Zupko K, Poria Y, Net P, Lazarovits J, Horn S, Khen M, Lancet D. Odorant signal termination by olfactory UDP glucuronosyl transferase. Nature. 1991;349(6312):790-793. https://doi.org/10.1038/349790a0.

Lemopoulos A, Montoya-Burgos JI. Whole genome assembly of the armored loricarid catfish *Ancistrus triradiatus* highlights herbivory signatures. Mol Genet Genomics. 2022;297(6):1627-1642. https://doi.org/10.1007/s00438-022-01947-6.

Leyhr J, Sanchez S, Dollman KN, Tafforeau P, Haitina T. Enhanced contrast synchrotron X-ray microtomography for describing skeleton-associated soft tissue defects in zebrafish mutants. Front Endocrinol (Lausanne). 2023;14:1108916. https://doi.org/10.3389/fendo.2023.1108916.

Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011;27(21):2987-2993. https://doi.org/10.1093/bioinformatics/btr509.

Li H. 2013. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM, arXiv, arXiv:1303.3997, preprint: not peer reviewed.

Liberles SD. Mammalian pheromones. Annu Rev Physiol. 2014;76(1):151-175. https://doi.org/10.1146/annurev-physiol-021113-170334.

Liebert A, Pang V, Bicknell B, McLachlan C, Mitrofanis J, Kiat H. A perspective on the potential of opsin as an integral mechanism of photobiomodulation: it's not just the eyes. Photobiomodul Photomed Laser Surg. 2022;40(2):123-135. https://doi.org/10.1089/photob.2021.0106.

Liman ER, Corey DP, Dulac C. TRP2: a candidate transduction channel for mammalian pheromone sensory signaling. Proc Natl Acad
Sci U S A. 1999:96(10):5791-5796. https://doi.org/10.1073/pnas. 96.10.5791.

Liu Y, Schröder J, Schmidt B. Musket: a multistage k-mer spectrum-based error corrector for illumina sequence data. Bioinformatics. 2013;29(3):308-315. https://doi.org/10.1093/bioinformatics/bts690.

Louis A, Muffato M, Roest Crollius H. Genomicus: five genome brows-ers for comparative genomics in eukaryota. Nucleic Acids Res. 2013;41 part 2:Database issue):D700-D705. https://doi.org/10.1093/ nar/gks1156.

Louis A, Nguyen NTT, Muffato M, Roest Crollius H. Genomicus up-date 2015: KaryoView and MatrixView provide a genome-wide perspective to multispecies comparative genomics. Nucleic Acids Res. 2015:43 part 2:Database issue):D682-D689. https://doi.org/ 10.1093/nar/gku1112.

Luo R, Liu B, Xie Y, Li Z, Huang W, Yuan J, He G, Chen Y, Pan Q, YunjieL, et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience. 2012:1(1):18. https://doi.org/10.1186/2047-217X-1-18.

Luu K, Bazin E, Blum MGB. Pcadapt: an R package to perform genome scans for selection based on principal component analysis. Mol Ecol Resour. 2017;17(1):67-77. https://doi.org/10.1111/ 1755-0998.12592.

Malik A, Almaharfi HA, Khan JM, Hisamuddin M, Alamery SF, Haq SH, Ahmed MZ. Protection of ζ-crystallin by α-crystallin under thermal stress. Int J Biol Macromol. 2021:167:289-298. https://doi.org/10.1016/j.ijbiomac.2020.11.183.

Manchanda A, Bonventre JA, Bugel SM, Chatterjee P, Tanguay R, Johnson CP. Truncation of the otoferlin transmembrane domain alters the development of hair cells and reduces membrane docking. Mol Biol Cell. 2021;32(14):1293-1305. https://doi.org/ 10.1091/mbc.E20-10-0657.

Manni M, Berkeley MR, Seppey M, Simão FA, Zdobnov EM. BUSCO update: novel and streamlined workflows along with broader and deeper phylogenetic coverage for scoring of eukaryotic, pro-karyotic, and viral genomes. Mol Biol Evol. 2021;38(10): 4647-4654. https://doi.org/10.1093/molbev/msab199.

Manuzzi A, Zane L, Muñoz-Merida A, Griffiths AM, Verissimo A. Population genomics and phylogeography of a benthic coastal shark (*Scyliorhinus canicula*) using 2b-RAD single nucleotide polymorphisms. *Biol J Linn Soc.* 2019;126(2):289-303. https:// doi.org/10.1093/biolinnean/bly185.

Marlétaz F, de la Calle-Mustienes E, Acemel RD, Paliou C, Naranjo S, Martínez-García PM, Cases I, Sleight VA, Hirschberger C, Marcet-Houben M, et al. The little skate genome and the evolutionary emergence of wing-like fins. Nature. 2023;616(7957): 495-503. https://doi.org/10.1038/s41586-023-05868-1.

Marra NJ, Stanhope MJ, Jue NK, Wang M, Sun Q, Pavinski Bitar P, Richards VP, Komissarov A, Rayko M, Kliver S, et al. White shark genome reveals ancient elasmobranch adaptations associated with wound healing and the maintenance of genome stability. Proc Natl Acad Sci U S A. 2019:116(10):4446-4455. https://doi. org/10.1073/pnas.1819778116.

Mayeur H, Lanoizelet M, Quillien A, Menuet A, Michel L, Martin KJ, Dejean S, Blader P, Mazan S, Lagadec R. When bigger is better: 3D RNA profiling of the developing head in the catshark *Scyliorhinus canicula*. Front Cell Dev Biol. 2021;9:744982. https://doi.org/10. 3389/fcell.2021.744982.

Mazan S, Jaillard D, Baratte B, Janvier P. Otx1 gene-controlled morphogenesis of the horizontal semicircular canal and the origin of the gnathostome characteristics. Evol Dev. 2000;2(4): 186-193. https://doi.org/10.1046/j.1525-142x.2000.00062.x.

McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9): 1297-1303. https://doi.org/10.1101/gr.107524.110.

Melis R, Vacca L, Cariani A, Carugati L, Cau A, Charilaou C, Di Crescenzo S, Ferrari A, Follesa MC, Hemida F, et al. Commercial

sharks under scrutiny: baseline genetic distinctiveness supports structured populations of small-spotted catsharks in the Mediterranean sea. Front Mar Sci. 2023;10:1050055. https://doi. org/10.3389/fmars.2023.1050055.

Metscher BD. MicroCT for comparative morphology: simple staining methods allow high-contrast 3D imaging of diverse non-mineralized animal tissues. BMC Physiol. 2009:9(1):11. https:// doi.org/10.1186/1472-6793-9-11.

Minařík M, Modrell MS, Gillis JA, Campbell AS, Fuller I, Lyne R, Micklem G, Gela D, Pšenička M, Baker CVH. Identification of multiple transcription factor genes potentially involved in the development of electrosensory versus mechanosensory lateral line organs. Front Cell Dev Biol. 2024;12:1327924. https://doi. org/10.3389/fcell.2024.1327924.

Minh BQ, Nguyen MAT, von Haeseler A. Ultrafast approximation for phylogenetic bootstrap. Mol Biol Evol. 2013;30(5):1188-1195. https://doi.org/10.1093/molbev/mst024.

Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, von Haeseler A, Lanfear R. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol. 2020;37(5):1530-1534. https://doi.org/10.1093/ molbev/msaa015.

Mirone A, Brun E, Gouillart E, Tafforeau P, Kieffer J. The PyHST2 hybrid distributed code for high speed tomographic reconstruction with iterative reconstruction and a priori knowledge capabilities. Nucl Instrum Methods Phys Res B. 2014;324:41-48. https://doi. org/10.1016/j.nimb.2013.09.030.

Miyawaki A, Homma H, Tamura H, Matsui M, Mikoshiba K. Zonal distribution of sulfotransferase for phenol in olfactory sustentacular cells. EMBO J. 1996;15(9):2050-2055. https://doi.org/10. 1002/j.1460-2075.1996.tb00558.x.

Modrell MS, Lyne M, Carr AR, Zakon HH, Buckley D, Campbell AS, Davis MC, Micklem G, Baker CV. Insights into electrosensory organ development, physiology and evolution from a lateral line-enriched transcriptome. Elife. 2017:6:e24197. https://doi.org/10. 7554/eLife.24197.

Morini M, Bergqvist CA, Asturiano JF, Larhammar D, Dufour S. Dynamic evolution of transient receptor potential vanilloid (TRPV) ion channel family with numerous gene duplications and losses. Front Endocrinol (Lausanne). 2022;13:1013868. https://doi.org/10.3389/fendo.2022.1013868.

Muffato M, Louis A, Poisnel C-E, Crollius HR. Genomicus: a database and a browser to study gene synteny in modern and ancestral genomes. Bioinformatics. 2010;26(8):1119-1121. https://doi. org/10.1093/bioinformatics/btq079.

Mulley JF, Hargreaves AD, Hegarty MJ, Heller R, Swain MT. Transcriptomic analysis of the lesser spotted catshark (*Scyliorhinus canicula*) pancreas, liver and brain reveals molecular level conservation of vertebrate pancreas function. BMC Genomics. 2014;15(1): 1074. https://doi.org/10.1186/1471-2164-15-1074.

Nakamura S, Nakamura N, Taniguchi K, Taniguchi K. Histological and ultrastructural characteristics of the primordial vomeronasal organ in lungfish. Anat Rec. 2012;295(3):481-491. https://doi. org/10.1002/ar.22415.

Nei M, Niimura Y, Nozawa M. The evolution of animal chemosensory receptor gene repertoires: roles of chance and necessity. Nat Rev Genet. 2008;9(12):951-963. https://doi.org/10.1038/ nrg2480.

Newton KC, Gill AB, Kajiura SM. Electroreception in marine fishes: chondrichthyans. J Fish Biol. 2019:95:135-154. https://doi.org/ 10.1111/jfb.14068.

Niimura Y, Matsui A, Touhara K. Extreme expansion of the olfactory receptor gene repertoire in African elephants and evolutionary dynamics of orthologous gene groups in 13 placental mammals. Genome Res. 2014;24:1485-1496. https://doi.org/10.1101/gr. 169532.113.

Nishimura O, Rozewicki J, Yamaguchi K, Tatsumi K, Ohishi Y, Ohta T, Yagura M, Niwa T, Tanegashima C, Teramura A, et al. Squalomix: shark and ray genome analysis consortium and its data sharing
platform. F1000Res. 2022:11:1077. https://doi.org/10.12688/f1000research.123591.1.

Norris HW. The distribution and innervation of the ampullae of lor-enzini of the dogfish, squalus acanthias. Some comparisons with conditions in other plagiostomes and corrections of prevalent errors. J Comp Neurol. 1929:47(3):449-465. https://doi.org/10.1002/cne.900470306.

Norris HW, Hughes SP. The cranial, occipital, and anterior spinal nerves of the dogfish, *Squalus acanthias*. J Comp Neurol. 1920:31(5):293-402. https://doi.org/10.1002/cne.900310502.

Oka Y, Saraiva LR, Korsching SI. Crypt neurons express a single V1R-related ora gene. Chem Senses. 2012:37(3):219-227. https://doi.org/10.1093/chemse/bjr095.

Okano T, Fukada Y. Photoreception and circadian clock system of the chicken pineal gland. Microsc Res Tech. 2001:53(1):72-80. https://doi.org/10.1002/jemt.1070.

O'Neill P, McCole RB, Baker CVH. A molecular analysis of neurogenic placode and cranial sensory ganglion development in the shark, *Sclyiorhinus canicula*. Dev Biol. 2007:304(1):156-181. https://doi.org/10.1016/j.ydbio.2006.12.029.

Oulion S, Debiais-Thibaud M, d'Aubenton-Carafa Y, Thermes C, Da Silva C, Bernard-Samain S, Gavory F, Wincker P, Mazan S, Casane D. Evolution of hox gene clusters in gnathostomes: insights from a survey of a shark (*Sclyiorhinus canicula*) transcriptome. Mol Biol Evol. 2010:27(12):2829-2838. https://doi.org/10.1093/molbev/msq172.

Paganin D, Mayo SC, Gureyev TE, Miller PR, Wilkins SW. Simultaneous phase and amplitude extraction from a single de-focused image of a homogeneous object. J Microsc. 2002:206(1): 33-40. https://doi.org/10.1046/j.1365-2818.2002.01010.x.

Panda AK, Nandi SK, Chakraborty A, Nagaraj RH, Biswas A. Differential role of arginine mutations on the structure and functions of α-crystallin. Biochim Biophys Acta. 2016:1860(1 Pt B): 199-210. https://doi.org/10.1016/j.bbagen.2015.06.004.

Peng G, Shi X, Kadowaki T. Evolution of TRP channels inferred by their classification in diverse animal species. Mol Phylogenet Evol. 2015:84:145-157. https://doi.org/10.1016/j.ympev.2014.06.016.

Peng YW, Robishaw JD, Levine MA, Yau KW. Retinal rods and cones have distinct G protein beta and gamma subunits. Proc Natl Acad Sci U S A. 1992:89(22):10882-10886. https://doi.org/10.1073/pnas.89.22.10882.

Pérez-Cerezales S, Boryshpolets S, Afanzar O, Brandis A, Nevo R, Kiss V, Eisenbach M. Involvement of opsins in mammalian sperm thermotaxis. Sci Rep. 2015:5(1):16146. https://doi.org/10.1038/srep16146.

Pickard GE, Sollars PJ. Intrinsically photosensitive retinal ganglion cells. In: Nilius B, Amara SG, Gudermann T, Jahn R, Lill R, Offermans S, Petersen OH, editors. Reviews of physiology, biochemistry and pharmacology. Berlin, Heidelberg: Springer Berlin Heidelberg; 2011. Vol. 162, p. 59-90.

Platt RN, Blanco-Berdugo L, Ray DA. Accurate transposable element annotation is vital when analyzing new genome assemblies. Genome Biol Evol. 2016:8(2):403-410. https://doi.org/10.1093/gbe/evw009.

Plotnikova OV, Kondrashov FA, Vlasov PK, Grigorenko AP, Ginter EK, Rogaev EI. Conversion and compensatory evolution of the γ-crystallin genes and identification of a cataractogenic mutation that reverses the sequence of the human CRYGD gene to an ancestral state. Am J Hum Genet. 2007:81(1):32-43. https://doi.org/10.1086/518616.

Policarpo M, Bemis KE, Laurenti P, Legendre L, Sandoz J-C, Rétaux S, Casane D. Coevolution of the olfactory organ and its receptor repertoire in ray-finned fishes. BMC Biol. 2022:20(1):195. https://doi.org/10.1186/s12915-022-01397-x.

Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010:26(6):841-842. https://doi.org/10.1093/bioinformatics/btq033.

Quintana-Urzainqui I, Rodríguez-Moldes I, Candal E. Developmental, tract-tracing and immunohistochemical study of the peripheral

olfactory system in a basal vertebrate: insights on Pax6 neurons migrating along the olfactory nerve. Brain Struct Funct. 2014;219(1):85-104. https://doi.org/10.1007/s00429-012-0486-2.

Ramírez-Amaro S, Picornell A, Arenas M, Castro JA, Massuti E, Ramon MM, Terrasa B. Contrasting evolutionary patterns in populations of demersal sharks throughout the western Mediterranean. Mar Biol. 2018:165(1):3. https://doi.org/10.1007/s00227-017-3254-2.

Ravi V, Lam K, Tay B-H, Tay A, Brenner S, Venkatesh B. Elephant shark (*Callochlorhinus milli*) provides insights into the evolution of hox gene clusters in gnathostomes. Proc Natl Acad Sci U S A. 2009:106(38):16327-16332. https://doi.org/10.1073/pnas.0907914106.

Read TD, Petit RA, Joseph SJ, MdT A, Weil MR, Ahmad M, Bhimani R, Vuong JS, Haase CP, Webb DH, et al. Draft sequencing and assembly of the genome of the world's largest fish, the whale shark: Rhincodon typus smith 1828. BMC Genomics. 2017:18(1):532. https://doi.org/10.1186/s12864-017-3926-9.

Redmond AK, Macqueen DJ, Dooley H. Phylotranscriptomics suggests the jawed vertebrate ancestor could generate diverse helper and regulatory T cell subsets. BMC Evol Biol. 2018:18(1):169. https://doi.org/10.1186/s12862-018-1290-2.

Ren X, Chang S, Laframboise A, Green W, Dubuc R, Zielinski B. Projections from the accessory olfactory organ into the medial region of the olfactory bulb in the sea lamprey (*Petromyzon marinus*): a novel vertebrate sensory structure. J Comp Neurol. 2009:516(2):105-116. https://doi.org/10.1002/cne.22100.

Rennison DJ, Owens GL, Taylor JS. Opsin gene duplication and divergence in ray-finned fish. Mol Phylogenet Evol. 2012:62(3): 986-1008. https://doi.org/10.1016/j.ympev.2011.11.030.

Rivera-Vicente AC, Sewell J, Tricas TC. Electrosensitive spatial vectors in elasmobranch fishes: implications for source localization. PLoS One. 2011:6(1):e16008. https://doi.org/10.1371/journal.pone.0016008.

Rodríguez-Moldes I, Carrera I, Pose-Méndez S, Quintana-Urzainqui I, Candal E, Anadón R, Mazan S, Ferreiro-Galve S. Regionalization of the shark hindbrain: a survey of an ancestral organization. Front Neuroanat. 2011:5:16. https://doi.org/10.3389/fnana.2011.00016.

Sakamoto K, Onimaru K, Munakata K, Suda N, Tamura M, Ochi H, Tanaka M. Heterochronic shift in hox-mediated activation of social hedgehog leads to morphological changes during fin development. PLoS One. 2009:4(4):e5121. https://doi.org/10.1371/journal.pone.0005121.

Salgado D, Mariluz BR, Araujo M, Lorena J, Perez LN, Ribeiro RL, Sousa JF, Schneider PN. Light-induced shifts in opsin gene expression in the four-eyed fish *Anableps anableps*. Front Neurosci. 2022:16:995469. https://doi.org/10.3389/fnins.2022.995469.

Sanchez S, Ahlberg PE, Trinajstic KM, Mirone A, Tafforeau P. Three dimensional synchrotron virtual paleohistology: a new insight into the world of fossil bone microstructures. Microsc Microanal. 2012:18(5): 1095-1105. https://doi.org/10.1017/S1431927612001079.

Saraiva LR, Korsching SI. A novel olfactory receptor gene family in teleost fish. Genome Res. 2007:17(10):1448-1457. https://doi.org/10.1101/gr.6553207.

Sato Y, Miyasaka N, Yoshihara Y. Mutually exclusive glomerular innervation by two distinct types of olfactory sensory neurons revealed in transgenic zebrafish. J Neurosc. 2005;25(20):4889-4897. https://doi.org/10.1523/JNEUROSCI.0679-05.2005.

Schluessel V, Bennett MB, Bleckmann H, Blomberg S, Collin SP. Morphometric and ultrastructural comparison of the olfactory system in elasmobranchs: the significance of structure-function relationships based on phylogeny and ecology. J Morphol. 2008:269:1365-1386. https://doi.org/10.1002/jmor.10661.

Schluessel V, Bennett MB, Bleckmann H, Collin SP. The role of olfaction throughout juvenile development: functional adaptations in elasmobranchs. J Morphol. 2010:271(4):451-461. https://doi.org/10.1002/jmor.10809.

Sendell-Price AT, Tulenko FJ, Pettersson M, Kang D, Montandon M, Winkler S, Kulb K, Naylor GP, Phillyll A, Fedrigo O, et al. Low
mutation rate in epaulette sharks is consistent with a slow rate of evolution in sharks. Nat Commun. 2023:14(1):6628. https://doi.org/10.1038/s41467-023-42238-x.

Sharma K, Syed AS, Ferrando S, Mazan S, Korschling S. The chemosensory receptor repertoire of a true shark is dominated by a single olfactory receptor family. Genome Biol Evol. 2019:11(2):398-405. https://doi.org/10.1093/gbe/evz002.

Slingsby C, Wistow GJ, Clark AR. Evolution of crystallins for a role in the vertebrate eye lens. Protein Sci. 2013:22(4):367-380. https://doi.org/10.1002/pro.2229.

Stanhope MJ, Ceres KM, Sun Q, Wang M, Zehr JD, Marra NJ, Wilder AP, Zou C, Bernard AM, Pavinski-Bitar P, et al. Genomes of endangered great hammerhead and shortfin mako sharks reveal historic population declines and high levels of inbreeding in great hammerhead. iScience. 2023:26(1):105815. https://doi.org/10.1016/j.isci.2022.105815.

Steenwyk JL, Buida TJ, Labella AL, Li Y, Shen X-X, Rokas A. PhyKIT: a broadly applicable UNIX shell toolkit for processing and analyzing phylogenomic data. Bioinformatics. 2021:37(16):2325-2331. https://doi.org/10.1093/bioinformatics/btab096.

Stowers L, Holy TE, Meister M, Dulac C, Koentges G. Loss of sex discrimination and male-male aggression in mice deficient for TRP2. Science. 2002:295(5559):1493-1500. https://doi.org/10.1126/science.1069259.

Su AI, Wiltshire T, Batalov S, Lapp H, Ching KA, Block D, Zhang J, Soden R, Hayakawa M, Kreiman G, et al. A gene atlas of the mouse and human protein-encoding transcriptomes. Proc Natl Acad Sci U S A. 2004:101(16):6062-6067. https://doi.org/10.1073/pnas.0400782101.

Syed AS, Sharma K, Policarpo M, Ferrando S, Casane D, Korschling S. Ancient and nonuniform loss of olfactory receptor expression renders the shark nose a *de facto* vomeronasal organ. Mol Biol Evol. 2023:40(4):msad076. https://doi.org/10.1093/molbev/msad076.

Tan M, Redmond AK, Dooley H, Nozu R, Sato K, Kuraku S, Koren S, Phillippy AM, Dove AD, Read T. The whale shark genome reveals patterns of vertebrate gene family evolution. Elife. 2021:10: e65394. https://doi.org/10.7554/eLife.65394.

Tanaka M, Münsterberg A, Anderson WG, Prescott AR, Hazon N, Tickle C. Fin development in a cartilaginous fish and the origin of vertebrate limbs. Nature. 2002:416(6880):527-531. https://doi.org/10.1038/416527a.

Tang H, Bowers JE, Wang X, Ming R, Alam M, Paterson AH. Synteny and collinearity in plant genomes. Science. 2008:320(5875): 486-488. https://doi.org/10.1126/science.1153917.

Terakita A. The opsins. Genome Biol. 2005:6(3):213. https://doi.org/10.1186/gb-2005-6-3-213.

Tulenko FJ, Augustus GJ, Massey JL, Sims SE, Mazan S, Davis MC. Hoxd expression in the fin-fold compartment of basal gnathostomes and implications for paired appendage evolution. Sci Rep. 2016;6(1):22720. https://doi.org/10.1038/srep22720.

Turner SD. Qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. J Open Source Softw. 2018:3(25):731. https://doi.org/10.21105/joss.00731.

Upton BA, Díaz NM, Gordon SA, Van Gelder RN, Buhr ED, Lang RA. Evolutionary constraint on visual and nonvisual mammalian opsins. J Biol Rhythms. 2021:36(2):109-126. https://doi.org/10.1177/0748730421999870.

Venkatachalam K, Montell C. TRP channels. Annu Rev Biochem. 2007:76(1):387-417. https://doi.org/10.1146/annurev.biochem.75.103004.142819.

Venkatesh B, Kirkness EF, Loh Y-H, Halpern AL, Lee AP, Johnson J, Dandona N, Viswanathan D, Tay A, Venter JC, et al. Survey sequencing and comparative analysis of the elephant shark (*Callorhinchus milii*) genome. PLoS Biol. 2007;5(4):e101. https://doi.org/10.1371/journal.pbio.0050101.

Venkatesh B, Lee AP, Ravi V, Maurya AK, Lian MM, Swann JB, Ohta Y, Flajnik MF, Sutoh Y, Kasahara M, et al. Elephant shark genome

provides unique insights into gnathostome evolution. Nature. 2014;505(7482):174-179. https://doi.org/10.1038/nature12826.

Wagner CI, Kopp MEL, Thorburn J, Jones CS, Hoarau G, Noble LR. Characteristics of the spiny dogfish (*Squalus acanthias*) nuclear genome. G3. 2023:13(9):jkad146. https://doi.org/10.1093/g3journal/jkad146.

Walsh CL, Tafforeau P, Wagner WL, Jafree DJ, Bellier A, Werlein C, Kühnel MP, Boller E, Walker-Samuel S, Robertus JL, et al. Imaging intact human organs with local resolution of cellular structures using hierarchical phase-contrast tomography. Nat Methods. 2021:18(12):1532-1541. https://doi.org/10.1038/s41592-021-01317-x.

Weadick CJ, Chang BSW. Molecular evolution of the lens crystallin superfamily: evidence for a retained ancestral function in N crystallins. Mol Biol Evol. 2009:26(5):1127-1142. https://doi.org/10.1093/molbev/msp028.

Weber JA, Park SG, Luria V, Jeon S, Kim H-M, Jeon Y, Bhak Y, Jun JH, Kim SW, Hong WH, et al. The whale shark genome reveals how genomic and physiological properties scale with body size. Proc Natl Acad Sci U S A. 2020:117(34):20662-20671. https://doi.org/10.1073/pnas.1922576117.

Yamaguchi K, Koyanagi M, Kuraku S. Visual and nonvisual opsin genes of sharks and other nonosteichthyan vertebrates: genomic exploration of underwater photoreception. J Evol Biol. 2021;34(6):968-976. https://doi.org/10.1111/jeb.13730.

Yamaguchi K, Koyanagi M, Sato K, Terakita A, Kuraku S. Whale shark rhodopsin adapted to deep-sea lifestyle by a substitution associated with human disease. Proc Natl Acad Sci U S A. 2023a:120(13):e2220728120. https://doi.org/10.1073/pnas.2220728120.

Yamaguchi K, Uno Y, Kadota M, Nishimura O, Nozu R, Murakumo K, Matsumoto R, Sato K, Kuraku S. Elasmobranch genome sequencing reveals evolutionary trends of vertebrate karyotype organization. Genome Res. 2023b:33(9):1527-1540. https://doi.org/10.1101/gr.276840.122.

Yamashita T. Unexpected molecular diversity of vertebrate nonvisual opsin Opn5. Biophys Rev. 2020:12(2):333-338. https://doi.org/10.1007/s12551-020-00654-z.

Yoo D, Park J, Lee C, Song I, Lee YH, Yun T, Lee H, Heguy A, Han JY, Dasen JS, et al. Little skate genome provides insights into genetic programs essential for limb-based locomotion. Elife. 2022:11: e78345. https://doi.org/10.7554/eLife.78345.

Yopak KE, Lisney TJ, Collin SP. Not all sharks are "swimming noses": variation in olfactory bulb size in cartilaginous fishes. Brain Struct Funct. 2015;220(2):1127-1143. https://doi.org/10.1007/s00429-014-0705-0.

Young JM, Massa HF, Hsu L, Trask BJ. Extreme variability among mammalian V1R gene families. Genome Res. 2010;20(1):10-18. https://doi.org/10.1101/gr.098913.109.

Yu G, Smith DK, Zhu H, Guan Y, Lam TT. GGTREE: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol. 2017;8(1):28-36. https://doi.org/10.1111/2041-210X.12628.

Zhang Y, Gao H, Li H, Guo J, Ouyang B, Wang M, Xu Q, Wang J, Lv M, Guo X, et al. The white-spotted bamboo shark genome reveals chromosome rearrangements and fast-evolving immune genes of cartilaginous fish. iScience. 2020;23(11):101754. https://doi.org/10.1016/j.isci.2020.101754.

Zhang Z, Sakuma A, Kuraku S, Nikaido M. Remarkable diversity of vomeronasal type 2 receptor (OlfC) genes of basal rayfinned fish and its evolutionary trajectory in jawed vertebrates. Sci Rep. 2022;12(1):6455. https://doi.org/10.1038/s41598-022-10428-0.

Zhou J, Liu A, He F, Zhang Y, Shen L, Yu J, Zhang X. Draft genome of white-blotched river stingray provides novel clues for niche adaptation and skeleton formation. Genomics Proteomics Bioinformatics. 2023:21(3):501-514. https://doi.org/10.1016/j.gpb.2022.11.005.