We used the previously described chicken and the Australian ghost shark (C. milii) TLR15 protein sequences (18) as queries to perform exhaustive blast searches on taxa representative of different vertebrate lineages, including mammals, reptiles/birds, amphibians, lungfishes, the coelacanth, ray-finned fish, cartilaginous fish (elasmobranchs and holocephalans), lampreys and hagfish genomes ( Supplementary Table 1 ). Additionally, for holocephalans, searches were also performed in unpublished genomic databases (Castro et al. in prep) covering the three taxonomic orders from Chimaeriformes, namely Callorhinchidae (Callorhinchus millii), Chimaeridae (Chimaera opalescens; Hydrolagus affinis, H. colliei, and H. mirabilis) and Rhinochimaeridae (Harriotta raleighana). All the protein sequences retrieved ( Supplementary Figure 1 , Supplementary Data 1 ) were aligned using Multiple sequence Comparison by Log-Expectation (MUSCLE) as implemented in Geneious Prime (http://www.geneious.com), and the position of indels were adjusted manually. Sequence alignment of putatively functional proteins can be found in Supplementary Figure 1 , while the alignment including putatively non-functional proteins, i.e., pseudogenes, can be found in Supplementary Figures 2A, B .

The genomic region surrounding the TLR15 gene was surveyed for all neighboring genes (up to three genes upstream and downstream of TLR15) in several taxa representative of the different vertebrate lineages, using the gene annotations available on NCBI for each taxon. This analysis allowed comparisons of the gene composition and order across vertebrates and insights into the putative conserved TLR15 synteny reported previously (18). When any of the conserved syntenic genes were absent in some species, we conducted additional blastn and blastx searches (with default parameters) on publicly available and unpublished genomes ( Supplementary Table 1 ). For those species in which no TLR15 ortholog was found, the genomic region between the two flanking genes common across vertebrates (ERLEC1 and GPR75) was subjected to blastx searches against NCBI non-redundant protein sequence database (nr). We performed this approach in different species of cartilaginous fishes, ray-finned fish, the coelacanth, lungfish, amphibians, turtles, crocodilians, birds, squamates, tuatara and mammals ( Supplementary Table 1 ; Figure 1 ).

To clarify the identification of TLR15 as an independent family, the full-length proteins ( Figure 2 ) of vertebrate TLR members of all TLR families were used and aligned (data not shown) using MUSCLE (as described above). A phylogenetic reconstruction of vertebrate TLRs relationships was performed in MEGAX using the ML method, with JTT+G+I as the best-fit amino acid substitution model (determined by MEGAX using ML as statistical method), and mid-point rooting. All positions with less than 95% site coverage were eliminated (i.e., fewer than 5% alignment gaps, missing data, and ambiguous bases were allowed at any position) using the partial deletion option. The final dataset used 352 positions from a total of 1800 positions from 259 amino acid sequences.

Genomic DNA sequences corresponding to the putative TLR15 orthologs were used for prediction of intron-exon boundaries and domain structure, namely signal peptide, ectodomain (ECD), transmembrane region (TM) and Toll/interleukin-1 receptor (TIR) domain, using the SignalP6.0 server (https://services.healthtech.dtu.dk/service.php?SignalP) (21), SMART (http://smart.embl-heidelberg.de/) (22) and TMHMM (https://services.healthtech.dtu.dk/service.php?TMHMM-2.0) (23). Delimitation of LRR motifs was performed using different tools: LRRpredictor (https://lrrpredictor.biochim.ro/) (24), LRRsearch (http://lrrsearch.com/index.php?page=tool) (25) and Conserved Domain Database (www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) (26). At the time of manuscript preparation, the commonly used tool LRRfinder was not available due to technical issues. Thus, we adopted a conservative approach for LRR motif delimitation: only the motifs detected by at least 2 of the above tools were considered. Additionally, we used the I-TASSER webserver (https://zhanggroup.org/I-TASSER/) (27) to predict the three-dimensional (3D) structure of TLR15 in lungfish and holocephalans and the modeled structure was displayed by PyMol (Schrödinger, LLC).

The evolutionary dynamics of amino acid substitution among TLR15 proteins was estimated using the ConSurf algorithm (28). For this we used the TLR15 amino acid alignment ( Supplementary Figure 1 ), and the chicken protein as query sequence. We allowed the algorithm to infer the phylogenetic tree using a maximum likelihood (ML) approach and the best evolutionary substitution model. The conservation scale retrieved is defined from the most variable positions to the most conserved positions ( Supplementary Figure 3 ).

To validate the pseudogene allele retrieved for H. affinis, and to confirm that the large insertion present in holocephalans were indeed part of the transcript (and not an intron), we performed PCR amplification from cDNA of H. affinis. Thus, we designed the primers Forward: 5’ GGAATTCTAGCAACTGAGGAGAAAGAGG 3’ and Reverse: 5’ GAAAGGTCCAGAATTTCAAGAGAG 3’, and the PCR was made with the Multiplex PCR Kit (Qiagen, Hilden, Germany) according to the manufacturer ‘ s protocol. Sequencing was performed on an ABI PRISM 3100 Genetic Analyzer (PE Applied Biosystems) and PCR products were sequenced in both directions ( Supplementary Figure 2B ). Due to sample limitations, we were not able to corroborate by PCR the pseudogenization in H. raleighana. Additionally, to test if the larger insertions in holocephalans were part of the TLR15 gene or correspond to an intron, the gene annotation software AUGUSTUS (v. 3.3.2) (29) and GeneMark (v. 3.61) (30) were used to check for the presence of introns in the holocephalan TLR15 sequences. Results including start and stop codons, exons and introns were produced in gff3 format.

The ratio (ω) of non-synonymous substitutions per non-synonymous sites (dN) over synonymous substitutions per synonymous sites (dS), dN/dS, was used to infer the selection pressures acting on TLR15 (Supplementary Tables 2 and 3). For this we used two ML frameworks, the CODEML program of Phylogenetic Analysis by Maximum Likelihood (PAML) 4.9 package (31, 32), and the HyPhy package implemented in the Datamonkey webserver (33, 34). In CODEML, a neighbour-joining tree of the TLR15 gene constructed in MEGAX (35) (with options: p-distances as the substitution model and complete deletion for gaps/missing data) was used as guide tree to compare the opposing site models M7 vs M8 using Likelihood Ratio Tests (LRT). While M7 (i.e. null model) assumes that ω ratios are distributed among sites according to a beta distribution allowing codons to evolve neutrally or under negative selection, M8 is an extension of the M7 model with an extra class of sites with an independent ω ratio freely estimated from the data allowing positive selection. Both, M7 and M8 models were compared by taking twice the difference in log likelihood between the two models, and the obtained value was assessed with a χ² distribution (df = 2) to test the null model (p<0.05). Amino acids detected as under positive selection were identified using the Bayes Empirical Bayes (BEB) approach, with posterior probability > 95%. BEB is the preferred approach because it accounts for sampling errors in the ML (31, 32, 36–38).

In the datamonkey Web Server, all the methods can take recombination into account. Thus, prior to the selection analysis we used the GARD module (39) to screen our sequences for recombination breakpoints. Since recombination breakpoints were detected for the TLR15 gene, we used the partitioned dataset obtained in GARD as input for the selection models. The nucleotide sequences of the TLR15 were analyzed under four available models: Single Likelihood Ancestor Counting (SLAC), Fixed-Effect Likelihood (FEL), Mixed Effects Model of Evolution (MEME) and Fast Unconstrained Bayesian AppRoximation (FUBAR). The SLAC model is based on the reconstruction of ancestral sequences and counts the number of d_S and d_N changes at each codon position of the phylogeny (40). FEL estimates ratios of d_N to d_S changes for each site in an alignment (40). MEME detects sites evolving under positive selection under a proportion of branches (41). FUBAR detects site-specific selection assuming that the selection pressure for each site is constant along the entire phylogeny (42). For SLAC, FEL and MEME the p-value was set to ≤0.05, while for FUBAR we used a posterior probability ≥ 0.95. For a more conservative approach, and as used previously (43, 44), only sites detected to be under positive selection in more than one ML method were considered.

Gene expression quantification was performed for TLR15 transcripts using bioinformatic mapping of the paired reads with RSEM (RNA-Seq by Expectation Maximization) (version 1.3.1.) by calling rsem-prepare-reference with specific parameter –bowtie2 and the rsem-calculate-expression with default parameters for paired-end reads (45). Read mapping was performed separately for C. milli and H. colliei for which RNAseq data was available from NCBI and from unpublished data (LFC Castro, in prep), using species-specific TLR15 transcripts as references ( Figure 3 ). TLR2 and TLR3 gene expression were also estimated for each species (Supplementary Figures 4, 5, respectively), following the procedure described above, to allow insights into the relative expression levels of TLR15 compared to other constitutively expressed TLRs, and to ensure that the observed read counts were not due to unbalanced representation of genes in the dataset. The final read counts are in transcripts per million (TPM) and fragments per kilobase of transcript per million of fragments mapped (FPKM).

We found single-copy TLR15 gene orthologs in two jawed vertebrate lineages: holocephalans and lungfishes. The retrieved protein sequences share the conserved synteny ( Figure 1 ) and cluster with other TLR15 members in a phylogenetic analysis ( Figure 2 ). Among holocephalans, TLR15 gene orthologs were found in all species studied, covering the three taxonomic families included in the Order Chimaeriformes ( Supplementary Figures 1, 2 ), namely Callorhinchidae (Callorhinchus millii), Chimaeridae (Chimaera opalescens; Hydrolagus affinis, H. colliei, and H. mirabilis) and Rhinochimaeridae (Harriotta raleighana) (accession numbers BK061675 and BK061828, and Supplementary data 1). Interestingly, the careful examination of the identified TLR15 gene sequences indicates that this gene is rendered non-functional in H. raleighana due to the occurrence of early sequence stop codons, while in H. affinis there is one functional allele and one rendered non-functional due to a frameshift mutation ( Supplementary Figures 2A, B ). The results in H. affinis were further confirmed by PCR amplification ( Supplementary Figure 2B ). Likewise, TLR15 orthologs were found in the two lungfish species (46, 47), which are representative of two extant families, Neoceratodontidae (Neoceratodus forsteri) (accession number BK061674) and Lepidosirenidae (Protopterus annectens). In NCBI, the Australian ghost shark C. milii TLR15 gene described here is incorrectly annotated as TLR6 (XM_042333155.1), while in the West African lungfish P. annectens the sequence presented here (Supplementary data 1) is an update (larger sequence encoding a signal peptide) to the incorrectly annotated TLR2 type-2 sequence (XM_044055095.1). The syntenic genomic block where TLR15 resides is conserved in vertebrates, although the flanking gene GPR75 is located elsewhere in all holocephalans analyzed ( Figure 1 ). Additionally, a TLR15-like remnant was also found in tuatara in the conserved syntenic block. This remnant contains 1462 nucleotides spanning the final nine LRR described for chicken, the C terminal, TM domain, and TIR domain; however, it has several frameshift mutations and early stop codons ( Supplementary Data 1 ).

TLRs are classified into seven families according to their ectodomain architecture and phylogenetic criteria (3, 8). In addition to the conserved synteny, our phylogenetic analyses support the identity of TLR15 in holocephalans and lungfish since the retrieved protein sequences cluster within the well-supported TLR15 clade both when using the full-length protein ( Figure 2 ) and the ECD only (data not shown). Our phylogenetic analyses also show, with good support, that TLR15 forms a distinct subfamily within vertebrate TLRs and that it is most closely related to the TLR1 subfamily.

Putative TLR15 orthologs could not be found in mammals, amphibians, the coelacanth, ray-finned fish and elasmobranchs, despite the conserved synteny of the genes surrounding TLR15 ( Figure 1 ). Likewise, our searches for putative TLR15 genes in jawless vertebrates did not retrieve any results. The searches were performed in several species of lamprey and in one hagfish ( Supplementary Table 1 ), but only partial CDS of the syntenic genes could be located ( Figure 1 ; Supplementary Table 4 ).

Motif prediction showed that the TLR15 in holocephalans and lungfish included the typical N-terminal signal peptide, ectodomain (ECD), transmembrane and TIR domains ( Supplementary Figure 6 ).

In birds and reptiles, the TLR15 ECD is generally composed of a N-terminal LRR, a C-terminal LRR and 19 additional LRR motifs (48). Probably due to the different approaches used, in this study we were only able to detect 18 LRRs for the chicken TLR15 ( Supplementary Figures 1 , Figure 4A ). Holocephalans and lungfishes present longer TLR15 proteins when compared to birds or reptiles ( Supplementary Figures 1, 6 ). Indeed, the full TLR15 protein in these two lineages has upwards of 1007 aa while the chicken TLR15 has only 868 aa ( Supplementary Figure 1 ). The major size differences are in the ectodomain region: taking the birds/reptiles TLR15 protein as reference, holocephalans exhibit 231 additional aa between the signal peptide and the first LRR motif described, and both holocephalans and lungfish have interval insertions of ~130 amino acids between the LRR3 and LRR4 ( Supplementary Figures 1, 4 ). Structurally, the exclusive 231 aa region of holocephalan TLR15 modifies the protein conformation ( Figures 4B–D ) by disrupting its horseshoe structure. By using H. affinis cDNA, we were able to amplify part of these 231 aa insertions. Additionally, when we tested this larger insertion with different annotation tools (AUGUSTUS and GeneMarker) they did not infer any intron. In turn, the extra ~130 aa shared by holocephalans and lungfish does not seem to affect the typical horseshoe-shaped solenoid structure ( Figures 4B–F ). Interestingly, when compared to the chicken ortholog, holocephalan and lungfish TLR15 genes had five and four extra LRR motifs (marked with an asterisk in Figures 4B–F ), respectively, with three of the extra LRRs being located in the inserted regions ( Supplementary Figure 1 and Figure 4 ). In contrast, the intracellular region of the protein (TIR domain) responsible for signal transduction is highly conserved across holocephalans, lungfish, birds and reptiles ( Supplementary Figures 1, 3 ).

In this work, we used RNAseq data to analyze the TLR15 expression levels in different tissues of two holocephalans, C. milii and H. colliei ( Figure 3 ). The patterns differed between species whereby TLR15 is expressed in all tissues of C. milii, with higher expression in gills, spleen, thymus and testis, while it is weakly expressed in H. colliei and only in heart, spleen, kidney and testis ( Figure 3 ). When compared to two constitutively expressed TLR2 and TLR3 genes, TLR15 shows relatively low expression in both holocephalans, but more so in H. colliei (Supplementary Figures 4, 5).

Here, we searched for signatures of selection in the functional TLR15 proteins using the corresponding nucleotide sequences from three birds, three crocodilians, four squamate reptiles, two lungfishes and four holocephalans. Our results detected 7 sites under positive selection and 292 sites under negative selection ( Supplementary Tables 2 , 3 , respectively, and Supplementary Figure 1 ). From the 7 positively selected codons (PSC), six are located in the ectodomain. In turn, from de 292 sites under negative selection, 85 sites were detected in the TIR domain, 21 sites were in the C-terminus region and 75 sites were detected on the ECD, specifically in the LRR motifs. With the exception of LRR3 (no sites under selection), all LRR motifs are under negative pressure, with LRR18 exhibiting the highest number of sites under negative selection (7 sites).