Discovery of All Three Types in Cartilaginous Fishes Enables Phylogenetic Resolution of the Origins and Evolution of Interferons

Interferons orchestrate host antiviral responses in jawed vertebrates. They are categorized into three classes; IFN1 and IFN3 are the primary antiviral cytokine lineages, while IFN2 responds to a broader variety of pathogens. The evolutionary relationships within and between these three classes have proven difficult to resolve. Here, we reassess interferon evolution, considering key phylogenetic pitfalls including taxon sampling, alignment quality, model adequacy, and outgroup choice. We reveal that cartilaginous fishes, and hence the jawed vertebrate ancestor, possess(ed) orthologs of all three interferon classes. We show that IFN3 groups sister to IFN1, resolve the origins of the human IFN3 lineages, and find that intronless IFN3s emerged at least three times. IFN2 genes are highly conserved, except for IFN-γ-rel, which we confirm resulted from a teleost-specific duplication. Our analyses show that IFN1 phylogeny is highly sensitive to phylogenetic error. By accounting for this, we describe a new backbone IFN1 phylogeny that implies several IFN1 genes existed in the jawed vertebrate ancestor. One of these is represented by the intronless IFN1s of tetrapods, including mammalian-like repertoires of reptile IFN1s and a subset of amphibian IFN1s, in addition to newly-identified intron-containing shark IFN1 genes. IFN-f, previously only found in teleosts, likely represents another ancestral jawed vertebrate IFN1 family member, suggesting the current classification of fish IFN1s into two groups based on the number of cysteines may need revision. The providence of the remaining fish IFN1s and the coelacanth IFN1s proved difficult to resolve, but they may also be ancestral jawed vertebrate IFN1 lineages. Finally, a large group of amphibian-specific IFN1s falls sister to all other IFN1s and was likely also present in the jawed vertebrate ancestor. Our results verify that intronless IFN1s have evolved multiple times in amphibians and indicate that no one-to-one orthology exists between mammal and reptile IFN1s. Our data also imply that diversification of the multiple IFN1s present in the jawed vertebrate ancestor has occurred through a rapid birth-death process, consistent with functional maintenance over a 450-million-year host-pathogen arms race. In summary, this study reveals a new model of interferon evolution important to our understanding of jawed vertebrate antiviral immunity.

Interferons orchestrate host antiviral responses in jawed vertebrates. They are categorized into three classes; IFN1 and IFN3 are the primary antiviral cytokine lineages, while IFN2 responds to a broader variety of pathogens. The evolutionary relationships within and between these three classes have proven difficult to resolve. Here, we reassess interferon evolution, considering key phylogenetic pitfalls including taxon sampling, alignment quality, model adequacy, and outgroup choice. We reveal that cartilaginous fishes, and hence the jawed vertebrate ancestor, possess(ed) orthologs of all three interferon classes. We show that IFN3 groups sister to IFN1, resolve the origins of the human IFN3 lineages, and find that intronless IFN3s emerged at least three times. IFN2 genes are highly conserved, except for IFN-γ-rel, which we confirm resulted from a teleost-specific duplication. Our analyses show that IFN1 phylogeny is highly sensitive to phylogenetic error. By accounting for this, we describe a new backbone IFN1 phylogeny that implies several IFN1 genes existed in the jawed vertebrate ancestor. One of these is represented by the intronless IFN1s of tetrapods, including mammalian-like repertoires of reptile IFN1s and a subset of amphibian IFN1s, in addition to newly-identified intron-containing shark IFN1 genes. IFN-f, previously only found in teleosts, likely represents another ancestral jawed vertebrate IFN1 family member, suggesting the current classification of fish IFN1s into two groups based on the number of cysteines may need revision. The providence of the remaining fish IFN1s and the coelacanth IFN1s proved difficult to resolve, but they may also be ancestral jawed vertebrate IFN1 lineages. Finally, a large group of amphibian-specific IFN1s falls sister to all other IFN1s and was likely also present in the jawed vertebrate ancestor. Our results verify that intronless IFN1s have evolved multiple times in amphibians and indicate that no one-to-one orthology exists between mammal and reptile IFN1s. Our data also imply that diversification of the multiple IFN1s present in the jawed vertebrate ancestor has occurred through a INTRODUCTION Antiviral immunity in jawed vertebrates is directed by interferons released by host cells in response to viral pathogens (1,2). Interferons are members of the class II α-helical cytokines along with interleukin (IL)-10,−19,−20,−22,−24, and −26 (hereafter called the IL-10 family), and are categorized into three classes, denoted as type I [e.g., IFN-α, β, κ, etc. in human [amongst others] and chicken], II (IFN-γ), and III (IFN-λs; known as IL-28A, IL-28B and IL29 in humans) interferons (hereafter called IFN1, IFN2, and IFN3), based on their receptors, genomic location, and sequence/structural homology (2,3). Roles beyond antiviral immunity have recently come to light for interferons, and IFN2 has been shown to contribute mainly toward defense against bacterial (especially mycobacteria), parasitic and fungal pathogens, leaving IFN1 and IFN3 as the main antiviral cytokines (1,2,(4)(5)(6)(7).
The evolutionary relationships between the three interferon classes, as well as intra-class evolutionary histories, have received considerable attention, but have proven difficult to resolve. The origins of the IFN3 lineage are particularly contentious. While some early studies suggested that the IFN3 and IFN1 lineages diverged in tetrapods, with teleost fishes possessing IFN1/3like molecules ( Figure 1A) (8)(9)(10), other studies suggested that teleosts possessed either IFN3 (11,12) or IFN1 orthologs (13). Later work, incorporating protein structures, showed that the teleost molecules were indeed IFN1s (14), and suggested that IFN3s likely emerged early in vertebrate evolution following whole genome duplication events (15)(16)(17). This scenario was supported by the discovery of IFN3 receptor homologs (along with those of IFN1 and IFN2) in cartilaginous fishes (18,19). However, other structure-based studies concluded that IFN3 is either part of the IL-10 family (specifically IL-22 or IL-19) ( Figure 1B) (3,(20)(21)(22), or sister to IFN1s ( Figure 1C) (15). A recent study has also revived the idea that IFN3 may have emerged from within IFN1 in the tetrapod ancestor ( Figure 1A) (23). Crucially, as no study has yet identified either the root of the class II α-helical cytokine family or orthologs of IFN3 genes outside tetrapods, none of these three hypotheses has been firmly rejected.
Evolutionary histories within each interferon class also remain unclear. For example, IFN2 is typically considered the most conserved interferon, however a tandem duplicate (IFN-γ-rel) has been found in some teleost fish-species (24), and phylogenetic analyses have failed to clarify whether this an ancient jawed vertebrate gene lost in other lineages (24)(25)(26), or teleost-specific (27). Multiple IFN3 genes often exist in individual species, but IFN3s are thought to be tetrapod-specific (28). However, very few studies have specifically focused on IFN3 evolution across vertebrates (28). The evolution of IFN1 genes, while better studied, also appears to be the most complicated. IFN1 genes are often present as lineage-specific clusters; for example, with the exception of IFN-κ, the IFN1 molecules of humans are evidently not directly orthologous with those of chickens (29)(30)(31). Clusters of lineage-specific IFN1s have also been observed in teleost fishes (32), classified as belonging to fish-specific group 1 or group 2 based on cysteine patterns (having two and four conserved cysteines, respectively) in the mature peptide sequence (9), and in amphibians (23). In fact, some phylogenetic analyses place all IFN1 sequences from mammals, teleosts, and amphibians into lineage-specific clades (9,(33)(34)(35), supporting a scenario where IFN1s evolve through concerted evolution (36,37). This would imply that high-turnover, lineage-specific gene gain and loss events, and/or gene conversion are major driving forces of IFN1 evolution (37,38). This is consistent with functional data, where individual genes appear to be specialized for defense against specific viruses. However, some studies have found phylogenetic relationships between IFN1s that are more difficult to interpret (23).
Poor resolution of interferon phylogenies hinders our ability to infer the history of evolutionary events including retro(trans)position, intron gains and losses, and changes in disulphide bridge structure. Amniote IFN1s are intronless and are classically thought to have arisen as a result of a retro(trans)position early in amniote evolution, as fish and amphibian interferons were found to contain four introns (2,9,19,30,33,39). Recent studies have revealed that both introncontaining and intronless IFN1 genes also exist in amphibians, leading to two competing hypotheses to explain the origins of intronless amphibian IFN1 genes; (i) they arose from the same event as the amniote intronless IFN1s (23,35), or (ii) they arose during independent retro(trans)position events (40). Intronless IFN3 genes also exist in mammals and amphibians (23), but whether they resulted from a single event or not remains to be tested. Similarly, two-and four-cysteine containing IFN1s exist in mammals and teleosts but it is thought that two-cysteine containing IFN1s emerged independently in each lineage; with intronless mammal and intron-containing teleost two-cysteine containing IFNs having lost a different cysteine pair (and hence disulphide bridge) from an ancestral four-cysteine containing IFN1 (19). Better resolution of IFN1 and IFN3 evolution could help determine both the frequency and timing of emergence of such features.
The primary amino acid sequence of interferons are short and rapidly evolving, both characteristics expected to promote phylogenetic error (41). Short alignments may have insufficient phylogenetic information to infer relationships between sequences and are more prone to stochastic errors. On the other hand, rapidly evolving sequences can be difficult to align and may induce systematic errors, resulting in long branch attraction (LBA) (42). Homoplasy (i.e., convergence due to hidden substitutions) is the best studied cause of LBA, and has previously been acknowledged as a concern when inferring immune gene phylogenies (41). Fortunately, it can often be counteracted by breaking long branches with additional taxa (43)(44)(45), applying site-heterogeneous models of evolution (42,46), removing fast-evolving sites (47), and/or identifying the best outgroup (48)(49)(50)(51), in addition to using outgroupfree methods (41,(52)(53)(54)(55) to root the tree. Compositional heterogeneity, resulting from differing codon usage preferences among sequences under comparison can also lead to LBA through non-phylogenetic similarity between lineages (56,57), and can be remedied by applying time-heterogeneous models of evolution (58,59), or removing compositionally biased sites or sequences (56,57,60,61). Other sources of systematic error have been identified (e.g., heterotachy, heteropecilly, nonindependence of sites), but are either less well studied or thought to have a less important effect on tree topology (62). Attempting to account for multiple sources of error, or applying several errorattenuating strategies at once is thought to improve phylogenetic accuracy (49,63,64), and this has proven successful for immune genes in the past (41,51).
Here, taking account of important phylogenetic considerations overlooked in past studies, we infer the origins and evolutionary history of interferons using a dataset that incorporates unprecedented sampling of both species and interferon diversity. Our findings offer a substantially overhauled model of interferon evolution and provide insights into the varied issues that hinder such studies, which have broader implications for immune gene phylogenetic analysis.

Homolog Identification and Characterization
TBLASTN (65) searches were carried out against a denselysampled set of genomes spanning chordate phylogeny (Table S1). An e-value cut-off of 10 was used in all searches, and sequences with either >75% identity compared to the query sequence [a set of known phylogenetically diverse IFNs were used for all searches, while known sequences (Table S2) from closely related species were also applied on an ad hoc basis], and/or with a top BLASTP (66) hit against an interferon in the NCBI non-redundant protein database, were retained for further analysis. To increase taxon sampling of cartilaginous fishes beyond elephant shark [until recently the only cartilaginous fish species with a sequenced genome (18,67)], transcriptome datasets for the small spotted-catshark were also analyzed from this lineage (41,68,69). Gene predictions were performed where a protein sequence was not already available, with the FGENESH+ webserver, using parameters for the closest related species available, and using either the blast hit or query as the homologous sequence (70). Structural homology prediction was performed through the Phyre2 protein structure prediction webserver using the "intensive" search option (71). Assessment of evolutionary conservation of sites necessary for IFN-λ3-receptor binding was achieved through visual comparison of multiple sequence alignment.

Multiple Sequence Alignment
All multiple sequence alignments were generated using PRANK, which has been shown to improve inference of insertions and deletions compared to other alignment approaches (72). This should help avoid alignment of non-homologous sites, reducing the potential for phylogenetic error. Manual curation was also performed (e.g., positions with no homologous amino acid in other classes were removed when examining inter-class relationships). Due to the rapidly evolving nature of IFN1s a set of high-quality, known sequences (Table S2) was used to build a base alignment before adding additional sequences from transcriptome and draft genome datasets, which may be truncated and more error prone, to the IFN1 dataset. Prior to analyzing this dataset, the PRANK alignment process was bootstrapped using GUIDANCE, which identifies sites that are not consistently aligned (73,74). Site alignments present in <93% of the GUIDANCE replicates were then removed to avoid use of unreliably aligned sites in phylogenetics (73,74). The "-add" function of MAFFT was then used with the L-INS-i approach to add new sequences to this high-quality core alignment (75). Alignment positions present in only a single species were then removed. See Data S1 for all multiple sequence alignment files.

Phylogenetic Analyses
All maximum likelihood phylogenetic analyses and model selection were performed in IQ-tree v1.6.7 (76). The Bayesian information criterion was used for model selection using IQtree's ModelFinder (77), and 1,000 ultrafast bootstrap replicates were generated to provide branch support values (78). IQ-tree was also used to detect compositionally biased sequences, using the built-in χ 2 test (71).
Outgroup-free rooted phylogenetic analyses were performed using a relaxed clock model, that permits root inference while accommodating rate variation among different tree branches (52). We have previously applied this approach to root other fast-evolving immune gene families (41,51,54,79,80), and it appears to work consistently for such datasets, except in the face of extreme rate asymmetry (41). This analysis was performed in BEAST v1.8.3 (81) applying an uncorrelated lognormal relaxed molecular clock model (52), a Yule speciation prior (82,83), and the best-fit amino acid substitution model (as inferred with IQ-tree). Two Markov chain Monte Carlo (MCMC) chains were run until effective samples sizes (>200) and convergence were sufficient, as assessed in Tracer v1.6 (84). Maximum clade credibility trees were generated in RootAnnotator (85).
Bayesian phylogenetic analyses incorporating outgroups were performed in PhyloBayes v4.1b (86), which also permits testing of site-heterogeneous models. Two MCMC chains were run for each analysis, until convergence was reached and effective sample sizes were sufficient for all statistics. This was assessed using the bpcomp (maxdiff < 0.3) and tracecomp (effsize > 50, and rel_diff < 0.3) programs within the PhyloBayes package (86).

Site-Heterogeneous Models, Cross-Validation, and Posterior Predictive Analyses
Site-heterogeneous models typically allow for better detection of homoplasy by accommodating site-specific evolutionary constraints in phylogenomic datasets (42). Such models have been applied with the objective of generating more reliable immune gene phylogenies (87,88), and have recently been shown to be capable of better explaining the site-specific evolutionary processes of aligned immune gene datasets (41). However, such models do not always provide a better fit for short alignments, and their relative fit cannot always be compared to standard models with the commonly applied information criteria. As such we used 10-fold cross-validation, as implemented in PhyloBayes (42,86), to compare the relative fit of a range of site-heterogeneous mixture models to the bestfitting standard model for the IFN1 dataset [JTT+Ŵ (89,90)]. The models tested included the infinite mixture model CAT (46), empirical derivations of CAT (C10/20/30/40/50/60) with limited numbers of site-categories, intended for gene family phylogenies (91), as well as an alternative site-heterogeneous model, WLSR5 (92), and a three-matrix substitution model, UL3, that loosely accommodates evolutionary process differences between structural features (91). Cross-validation relies on randomly partitioning the alignment into equal sized subsamples (10 here, as the analysis is 10-fold), before one of these subsamples is used for validation to test the model, while the rest are combined as a training set. This process is then repeated, using each of the other subsamples as the validation dataset, and then the average results are used for comparison against other models. We ran each individual chain (i.e., one chain for each of the 10 training sets for each model) for 1,000 points, using the first 100 as burn-in (i.e., 10 chains for each model tested).
In addition to assessing relative model fit, posterior predictive simulations (PPS) were also performed in PhyloBayes to determine if the model applied could adequately describe the real data for the tested statistic (42,86). This approach consists of generating simulated data under the model in question, for comparison against the observed (i.e., real) data. Here, PPS was used to investigate the ability of models to account for homoplasy and compositional heterogeneity across lineages in the IFN1 dataset (42,56,57,60). The compositional heterogeneity test was used to generate a second IFN1 dataset by identifying and removing sequences that deviate significantly from the assumption of homogeneity, measured at Z-scores < −2 and >2 (the default in PhyloBayes, which is slightly more inclusive than P < 0.05 cut-off). All PPS analyses were specifically performed under JTT+Ŵ, as this model should be the most susceptible to error compared to the tested mixtures, so we viewed this as more conservative. Finally, we also tested a time-heterogeneous (59), and a site-and time-heterogeneous (58) model for the IFN1 dataset, but these analyses failed to reach convergence despite running continuously for more than three months each.

Testing Exacerbation of Potential Errors in Interferon Phylogenetics Analysis
Multiple approaches were tested to induce phylogenetic error (93) in the rapidly-evolving IFN1 family to better explain the discrepancy in the results of past studies, as well as the difficulty in inferring IFN1 evolutionary history. This included applying a more distantly related outgroup in place of the closest related outgroup to root the tree (50,51,94), inferring the phylogeny under less well-fitting substitution models (50,51,54,93,95), including sequences that introduce significant compositional bias in the analysis (63,64), as well as sequence removal to lengthen target branches and increase the potential for LBA (93,96).

A Cartilaginous Fish IFN-λ
Reciprocal BLAST searches of a multi-tissue small-spotted catshark transcriptome (41) revealed a putative cartilaginous fish IFN3 sequence. Characterization of this sequence by multiple sequence alignment against human interferon sequences and IL-10 (as a representative of the broader IL-10 family) support this assignment (Figure 2A and Figure S1). Additionally, the signature disulfide bridge-forming Cys pair were present at the C-terminus of this sequence (Figure 2A and Figure S1) (21,22). Interestingly, the most important receptor binding sites of human IFN-λ are poorly conserved in this sequence (22); including Phe158, which is vital for human IFN-λ3-receptor interaction (22). However, this may not preclude antiviral functionality of catshark IFN3, considering that this residue is FIGURE 2 | (A) Analysis of key residues in the catshark IFN-λ sequence compared to a set of other interferons and IL-10. Sequences are represented as cartoon bars, which are relatively scaled according to amino acid sequence length. Arrows denote the end of the signal peptide region, while disulphide bridges are shown as connected regions underneath each cartoon bar, with "C" in the C-terminal region being an unpaired Cys from the characteristic C-terminal disulphide bridge of IFN3s (22). Above the bar the most important residues for IFN-λ3 receptor binding are shown (22). Residues filled in black are conserved, whereas residues filled white are not well conserved, and gray-filled residues involve conserved replacements (e.g., K→ R). The bar over the VXXQ motif of catshark IFN-λ indicates that this is not aligned perfectly to human IFN3s, while the star indicates that mutation of this residue abolishes binding in human IFN-λ3 (22). See Figure S1 for full alignment. (B) Relaxed clock (uncorrelated lognormal) rooted class II α-helical cytokine family phylogeny under JTT + I + τ and a Yule speciation prior. The tree is rooted at the best supported root position. Root posterior probabilities (RPP) are shown for branches with a non-negligible probability (i.e., posterior probability <0.05) of being the root. Posterior probabilities are also shown for key nodes, and clades representing individual family members, or the entire IL-10 family have been collapsed to emphasize deep relationships within the family.
also not conserved in human IFN-λ4 (Figure 2A and Figure S1), which appears to be capable of binding the IFN3 receptor (97). Further, our preliminary analyses suggest that catshark IFN-λ is involved in antiviral defense (unpublished data). Submission to the Phyre2 protein structure prediction server, an approach which has previously been employed to aid orthology assignment of fast-evolving immune genes in cartilaginous fishes (98), also indicated a best-structural match of the putative catshark sequence to mammalian IFN-λs ( Figure S2). Finally, phylogenetic analysis (see next section), verified this assignment; indeed, catshark IFN-λ forms a clade with tetrapod IFN3s, to the exclusion of other class II α-helical cytokines, with maximal support (posterior probability [PP] = 1.00; Figure 2B, Figure S3). The existence of a cartilaginous fish IFN3 allows us to unequivocally reject the hypothesis that IFN3s emerged by duplication from IFN1s in the tetrapod ancestor ( Figure 1A) (8,10,23).

Deep Relationships Within the Class II α-Helical Cytokine Family
To understand the evolutionary relationships between the interferon classes and other class II α-helical cytokines, and identify the closest outgroups to best infer within-class interferon relationships, we performed a phylogenetic analysis of the full class II α-helical cytokine family. A relaxed clock model was used to root this tree, as the deeper phylogenetic origins of the family, and thus potential outgroups, are not known. Our analysis supports a sister group relationship between IFN1 and IFN3 (PP = 0.91; Figure 2B and Figure S3), while on the other side of the tree root (root posterior probability [RPP] = 0.82; Figure 2B and Figure S3), a monophyletic IL-10 family is sister to IFN2 (PP = 0.9; Figure 2B and Figure S3). These findings reject the hypothesis that IFN3 is part of the IL-10 family (14,21,22) (Figure 1B), while the root placement suggests that the deepest divergence in the class II α-helical cytokines separates the main antiviral interferons from the rest of the family, consistent with the model of class II α-helical cytokine evolution proposed by Siupka et al. (15). A second root position placing IFN1 as sister to the rest of the family could not be rejected however, although this was only very weakly supported (PP = 0.05; yielding a 16:1 weighting in favor of the best root) ( Figure 2B). These results concur with the conclusion that the IFN1/3s previously identified outside tetrapods are in fact true IFN1s, consistent with their structural and functional features (14). Further, by supporting a sister group relationship between IFN1 and IFN3 (Figure 1C), our findings indicate that IFN1 and IFN3 can be used as reciprocal outgroups in phylogenetic analyses, enabling outgroup-rooted IFN1 and IFN3 phylogenies without the inclusion of more distant, and potentially-biasing, outgroups like IFN2 and/or the IL-10 family (which by the same rationale can be applied as outgroups for each other).

IFN2 Evolution Indicates That IFN-γ-Rel Is Teleost-Specific
IFN2 is the most structurally conserved of the interferon classes and is thought to have the simplest evolutionary history. Despite being present in single copy across most of vertebrate phylogeny, an additional gene, IFN-γ-rel (24), is present in tandem to IFNγ in some teleosts. Phylogenies in some previous studies suggest IFN-γ-rel could be an ancient lineage, lost from other vertebrates (24)(25)(26). Not all phylogenies support this however (27), and it has also been suggested that IFN-γ-rel arose through duplication of IFN-γ during teleost evolution (19). Here, we tested this using a PRANK alignment of IFN-γ sequences spanning jawed vertebrate phylogeny, as well as the best-fitting substitution model, and the most closely related outgroup, the IL-10 family. This phylogenetic analysis maximally supported IFN-γ-rel as sister to teleost IFN-γ (Ultrafast Bootstrap [UFBOOT] = 100%) (Figure 3). Together with its absence outside of teleosts, this indicates that IFN-γ-rel resulted from a teleost-specific tandem gene duplication.

Divergence of Human IFN-λs and Convergent Intron Loss in IFN3 Evolution
In light of the newly discovered catshark IFN-λ, and identification of the IFN1 family as the closest outgroup, we reassessed the evolutionary history of the IFN3s (Figure 4 and Figure S4). Along with a selection of known tetrapod IFNs, BLAST analyses of genomes spanning vertebrate phylogeny revealed putative new amphibian and reptile IFN3 sequences, but consistent with previous studies failed to identify coelacanth or teleost IFN3s. Phylogenetic analysis of this dataset revealed that catshark IFN-λ falls sister to all tetrapod IFN3s (UFBOOT = 91%) (Figure 4 and Figure S4), while within tetrapods, amphibians and amniotes form separate sister clades (UFBOOT = 88%) (Figure 4 and Figure S4). Our analyses verified the presence of intronless IFN3s in amphibians (23). Strikingly however, we found that the intronless IFN3s of amphibians and mammals emerged independently, and that intronless IFN3s have evolved at least twice during amphibian evolution (UFBOOT = 100%) (Figure 4 and Figure S4), and hence at least three times throughout vertebrate evolution. Within amniotes, our results are largely consistent with those of Chen et al. (28), as we find that reptiles have at least two IFN3 lineages. In our analyses these lineages form clades with mammalian IFN-λ4 (though reptiles are not monophyletic) (UFBOOT = 97%) and mammalian IL-28/29 (UFBOOT = 82%), suggesting that they are orthologous, and that the IL-28/29 and IFN-λ4 lineages split in the amniote ancestor (Figure 4 and Figure S4). The reptile IL-28/29-like gene appears to have been duplicated in the ancestor of archelosaurians (turtles, birds, and crocodiles) (UFBOOT = 97%), while the human IL-28 and IL-29 lineages appear to have been duplicated in placental mammals (UFBOOT = 98%), with IL-28A and IL-28B later splitting during primate evolution (Figure 4 and Figure S4).

Accounting for Phylogenetic Errors to Generate a Reliable IFN1 Tree
The evolutionary history of IFN1s has been studied intensively and many very different tree topologies generated. However, previous studies have not intentionally accounted for any of the major known sources of phylogenetic error. To help counter this we first applied two data-centric approaches designed to combat phylogenetic errors (43-45, 48-50, 94, 96). First, we applied only the closely related IFN3 as an outgroup, and increased taxon sampling by identifying new IFN1s from a dense sample of genomes across vertebrate phylogeny. This revealed hundreds of new IFN1 sequences (Table S2), which we subsampled prior to phylogenetic analyses, keeping only sequences above 100 amino acids in length (except for cartilaginous fishes and Japanese eel, where sequences of 50 or more amino acids were retained, given the paucity of data available for these species and their important evolutionary placement within jawed vertebrates and teleosts), and removing highly similar sequences within species from densely sampled lineages to reduce computational burden without negatively affecting deeper nodes in the tree.
Next, we tested the utility of site-heterogeneous phylogenetic mixture models to resolve the IFN1 phylogeny; such models have been shown to offer an improved fit to many datasets (42,49,95), including immune genes (41), as well as being more resistant to LBA artifacts (42) by accounting for evolutionary process variation among sites. As such we ran PhyloBayes analyses under the best-fitting standard model, JTT+Ŵ, as well as under a variety of site-heterogeneous mixture models (46,91,92,99). Unfortunately, none of these runs reached convergence. However, this is not an uncommon occurrence in PhyloBayes analyses of difficult datasets, and can be remedied by identifying and removing error causing sequences or branches (49).
Interrogation of the MCMC chains for the JTT+Ŵ analyses revealed that effective sample sizes were sufficient, but that individual runs had become "stuck" at different log-likelihoods (the lesser of which must represent a local optimum), and at different tree topologies (i.e., PhyloBayes bpcomp "maxdiff = 1"). Trying to resolve this objectively, while also reducing systematic error, we used PPS to detect and remove sequences that deviated from the assumption of compositional heterogeneity (60). We then re-ran the phylogenetic analyses using this reduced "compositionally homogeneous" dataset (CHOM) and found that runs now converged for JTT+Ŵ ( Figure 5A and Figure S5), as well as all the tested mixture models. Before examining the resultant tree topologies however, we sought to gain further understanding of the difference between analyzing the full and CHOM datasets, and to determine if the site-heterogeneous models might provide a better fit than JTT+Ŵ.
Interrogation of the compositional heterogeneity PPS results for the full dataset showed consistency between chains, with almost fully overlapping sets of sequences identified as biased in both chains, implying that sequence removal (i.e., identification of compositional bias) was not affected by lack of convergence ( Figure 5B; Table S3), and therefore should be reliable.
An additional consideration is that removal of sequences may serve to lengthen branches in the phylogenetic tree, reducing the ability to detect hidden substitutions and increasing the potential of LBA artifacts. To determine if this was an issue we performed PPS analyses to test whether JTT+Ŵ could adequately accommodate homoplasy in both the full and CHOM datasets (42). While a greater level of homoplasy was both observed and predicted in the full dataset, this was adequately predicted by JTT+Ŵ (Figure 5C), implying that it should not be a major source of error (including topological inconsistency between chains), in either dataset.
Finally, using Bayesian cross-validation (42) we determined that JTT+Ŵ was in fact better-fitting than any of the tested mixture models (Figure 5D), perhaps due to the short alignment length, and so only this tree was used to make evolutionary inferences.

Birth-Death Evolution of Multiple IFN1 Genes Since the Jawed Vertebrate Ancestor
Phylogenetic analyses of the CHOM dataset under JTT+Ŵ (Figure 5A and Figure S5) suggest a new paradigm for IFN1 evolution. The resultant tree indicates that duplication and loss events have occurred frequently since the origins of IFN1s (Figure 5A and Figure S5). This fits a rapid birth-death model of evolution (100,101), as proposed for salmonid IFN1s (32), rather than the concerted evolution model (i.e., IFN1 expansions and contractions are confined to specific lineages) implied by many previous phylogenies. Multiple features of the tree topology support this scenario, including the presence of an amphibian lineage (red star in Figure 5A) that falls sister to all other IFN1s (PP = 0.65), intimating this lineage existed in the jawed vertebrate ancestor and has since been lost in other jawed vertebrates (Figure 5A; Figure S5). In addition, we identified a new IFN1 locus in elephant shark (Table S2), the intron-containing genes of which form the sister group to a clade (PP = 0.82) containing all amniote intronless interferons (PP = 0.74), suggesting that these genes are orthologous and that other jawed vertebrate classes have lost orthologs of this gene (Figure 5A and Figure S5). This means that the intron-containing IFN1s of teleosts and amphibians which gave rise to the idea that the retrotransposition event occurred in the amniote ancestor are in fact paralogous to this lineage, and as such are not informative on this point. Within the intronless amniote interferon clade, reptiles possess lineagespecific expansions comparable to those seen in mammals, and may retain ancient amniote IFN1s lost in mammals ( Figure 5A). Elsewhere in the tree, we find that IFN-f, previously found only in teleosts, is likely also present in amphibians (PP = 0.74), and possibly cartilaginous fishes (PP = 0.53), where a lineagespecific expansion has occurred (green star in Figure 5A). Thus IFN-f appears to be an ancient jawed vertebrate IFN1 lineage, secondarily lost in amniotes. Importantly, this suggests that the ray-finned fish IFN1s may not be monophyletic and that the two groups defined by cysteine structures (9) may not have a phylogenetic basis (Figure 5A and Figure S5). Despite this, non-IFN-f ray-finned fish IFN1s group together (PP = 0.59), though, due to a polytomy, they do not appear to be identifiably orthologous to IFN1s from any other jawed vertebrate lineage (Figure 5A and Figure S5). The relationships of coelacanth IFN1s are similarly unresolved (Figure 5A and Figure S5). This result seems most consistent with both lineages representing lineage-specific expansions of ancient jawed vertebrate IFN1 genes lost in other jawed vertebrates, though the support for any showed that ∼50 IFN sequences introduced significant potential for compositional bias, which were removed to minimize branching artifacts (i.e., forming the CHOM dataset used for part (A). (C) PPS also shows that JTT adequately predicts homoplasy in both the full and CHOM datasets. (D) Model selection via 10-fold Bayesian cross-validation indicates that the site-homogenous JTT model fits the data better than a range of site-heterogeneous mixture models. (E) IFN1 topology is highly sensitive to both dataset bias and methodological error: (i) the full (i.e., compositionally heterogeneous) dataset places the cartilaginous fish group otherwise identified as sister to amniote IFNs in a monophyletic group with the amphibian sequences that form the sister group to all other IFN1s (see also Figures S6, S7), while (ii) and (iii) show that less well fitting models (see also Figure S8) and distant outgroup taxa (full tree in Figure S9) result in evolutionarily irreconcilable root placement.
Frontiers in Immunology | www.frontiersin.org scenario is low. Taken together, these results imply that several distinct IFN1 genes existed in the jawed vertebrate ancestor and have undergone rapid birth-death evolution since, meaning that ancient interferon genes are sometimes retained in only one or very few extant descendant taxa, while at the same time lineage-specific interferon expansions and contractions are common.

Sensitivity of IFN1 Family to Phylogenetic Error
Major differences were observed between our results and those of previous studies. While we believe this is a result of improved methodology, we attempted to formally test this by performing experiments designed to exacerbate error potential in phylogenetic analyses (93). First, given that sequences displaying compositional bias contributed to non-convergence of PhyloBayes analyses, we instead built the full IFN1 phylogeny using alternative software (BEAST). This produced a similar topology to that obtained for the PhyloBayes CHOM analysis, however, the cartilaginous fish lineage that fell sister to the intronless amniote IFNs in the CHOM PhyloBayes analysis was instead placed sister to the amphibian sequences that fell sister to all other IFN1s (PP = 0.5) (Figure 5E and Figure S6). As the CHOM dataset does not deviate from the assumption of compositionally homogeneity, we considered this result to be an error induced by compositional bias. To further explore stability of this cartilaginous fish lineage in the CHOM dataset, we pruned sequences contributing to nearby branches to lengthen this branch, but this did not perturb its placement in the tree ( Figure S7). Second, we examined tree topologies generated under the less well-fitting mixture models (Figure S8). Even for the second best-fitting model, UL3, this resulted in major issues with root placement (Figure 5E and Figure S8), suggesting an extremely non-parsimonious evolutionary scenario. Third, a similar outcome was observed when more distantly related IFN2 was applied as the outgroup instead of IFN3 (Figure 5E and Figure S9). Collectively, these results suggest that IFN1 phylogeny is highly sensitive to methodological and sampling errors.

Intronless IFN1s Emerged in the Tetrapod Ancestor and Multiple Times in Amphibians
Since performing our IFN1 analyses, recent studies have identified new sequences not present in our dataset that may be relevant to IFN1 evolution. However, given the large compute time of performing all our PhyloBayes analyses (i.e., including cross-validation and PPS), it was not practical to rerun these with the addition of the new sequences (102). Instead, we decided upon a reasonable compromise; given that the best-fit model in PhyloBayes was a standard site-homogeneous model, we added the relevant sequences to our alignment and ran this under JTT+Ŵ in IQ-tree, with compositionally biased sequences (as identified by the χ 2 test implemented in IQ-tree) removed.
This dataset (hereafter called EXT) included the recently identified fish IFN-h (103), as well as additional IFNs (both intron-containing and intronless) from amphibians. The EXT phylogenetic tree (Figure 6A and Figure S10) is generally consistent with that of the CHOM analysis, except that non-IFN-f ray-finned fish IFN1s, coelacanth IFN1s, and IFN-f form a weakly supported clade (UFBOOT = 52%) rather than a polytomy (i.e., PP < 0.5 in PhyloBayes analyses) (Figures 5A, 6A; Figures S5, S10). Within this clade, non-IFNf ray-finned fish IFN1s fall sister to IFN-f (UFBOOT = 79%) with the coelacanth IFN1s being sister to both (Figure 6A and Figure S10). Because the IFN-f clade includes cartilaginous fish (UFBOOT = 79%) and amphibian (UFBOOT = 86%) sequences, this is consistent with non-IFN-f ray-finned fish IFN1s being the only surviving lineage of an interferon gene that was present in the jawed vertebrate ancestor (Figure 6A and Figure S10). A similar evolutionary scenario can thus be applied to coelacanth IFN1s, but support for this is weaker (UFBOOT = 52%) (Figure 6A and Figure S10). The newly included IFN-h falls within the clade of non IFN-f teleost IFN1s (UFBOOT = 93%), and as such does not alter the backbone IFN1 phylogeny (Figure 6A and Figure S10). Similarly, despite being placed differently in past analyses (23,35), we find that almost all of the recently identified amphibian IFN1s (23,35) fall into the clade of amphibian sequences (UFBOOT = 100%) that is sister to all other IFN1s (UFBOOT = 88%), in the CHOM analysis (Figures 5A, 6A; Figures S5, S10). Within this clade, the deepest split falls between intronless Xenopus IFN1s, and a clade containing intron-containing Xenopus and Nanorana parkeri sequences, as well as intronless N. parkeri sequences, confirming the recently discovered independent origins of intronless IFN1s in these species (40) (Figure 6B and Figure S11).
Strikingly, a small number of intronless amphibian IFN1s were nested within the mammal and reptile IFN1 clade, falling sister to a clade containing only reptile sequences (UFBOOT = 89%) (Figure 6A and Figure S10). This suggests that orthologs of amniote intronless IFN1s are present in amphibians and arose in the ancestor of tetrapods. Within this intronless tetrapod clade, two additional ancient reptile lineages are also present, one of which forms the sister group to all mammalian IFN1s ( Figure 6A and Figure S10), while the other forms the sister to all other intronless IFN1s (i.e., both former reptile clades, and their mammalian and amphibian intronless counterparts) (UFBOOT ≥ 74%) (Figure 6A and Figure S10). This is consistent with a birth-death model of evolution, where reptiles have retained genes from three ancient intronless lineages that were present in the ancestor of tetrapods, but with amphibians and mammals retaining only one of these each, before the onset of independent lineage-specific diversifications. Intriguingly, an amphibian interferon containing a single intron falls sister to the group of cartilaginous fish IFN1s that were sister to the intronless amniote IFN1s in the CHOM analysis (UFBOOT = 84%) and together they fall sister to the intronless tetrapod IFN1s (UFBOOT = 47%). If accurate, this suggests that these cartilaginous fish genes are paralogous rather than orthologous to mammalian IFN1s, as both clades contain amphibians, further increasing the number of IFN1s likely present in the jawed vertebrate ancestor (Figure 6A and Figure S10).

No One-to-One Orthology Relationships Between Mammal and Reptile IFN1s
It has long been recognized that the IFN-α and IFN-β genes of human and chicken are not orthologous (29). In contrast, the recently discovered chicken IFN-κ is purportedly an ortholog of mammalian IFN-κ (31). Interestingly, our CHOM and EXT IFN1 datasets, which greatly expanded taxon sampling in reptiles, failed to find evidence for orthology between IFN-κ genes of mammals and reptiles, but did not include the lineage containing chicken IFN-α because this was compositionally biased (Figures S5, S10; Table S3). Similarly, a lone amphibian sequence containing a single intron grouped together with the cartilaginous fish sequences that fall sister to the tetrapod intronless interferon clade. As this sequence would, more parsimoniously, be expected to group with the intronless IFN1s we performed more focused phylogenetic analyses to examine this finding. Our analyses included the cartilaginous fish and amphibian sequences that fell sister to this group in the EXT analysis ( Figure 6A and Figure S10), but not more distantly related IFN1s to avoid biases introduced by distant outgroups. We also reinstated sequences, including chicken IFN-α, that were excluded from CHOM and EXT due to compositional bias. Interestingly, in this instance the amphibian sequence sister to cartilaginous fish in CHOM and EXT grouped with the intronless IFN1s of other amphibians (UFBOOT = 68%), away from the cartilaginous fish sequences (UFBOOT = 100%). This, far more parsimonious scenario, verifies the cartilaginous fish sequences as orthologs of the intronless tetrapod IFNs (Figure 7A and Figure S12). No evidence for orthology between any mammalian and reptile IFN1s was observed in this analysis. If rooted with the cartilaginous fish sequences, the results are also consistent with reptile genomes harboring ancient tetrapod intronless interferon lineages lost in mammals ( Figure 7A and Figure S12). Finally, an unrooted analysis (i.e., excluding cartilaginous fish and amphibian sequences) recovered independent mammal and reptile clans, further supporting the lack of orthology between and reptile and mammalian IFNs ( Figure 7B and Figure S12).

Group 1, but Not Group 2, Ray-Finned Fish IFN1s Are Monophyletic
Our IFN1 phylogenies consistently showed that IFN-f is not a member of the ray-finned fish-specific IFN1s (Figures 5A, 6A;  Figures S5, S10). This suggests that IFN1 classification based on conserved cysteine pairs may not have a phylogenetic basis. For example, group 2 IFNs (IFN-b, IFN-c, and IFN-f) do not form a clade despite all having two conserved cysteine pairs in the mature peptide. To better explore this, we performed a focused phylogenetic analysis (Figure 8 and Figure S13) of the remaining ray-finned fish-specific IFN1s that formed a clade in our CHOM and EXT analyses, using IFN-f as an outgroup. This placed the root between the remaining group 2 and group 1 members, in agreement with past hypotheses of fish IFN1 evolution, except for IFN-f (UFBOOT ≥ 77%) (Figure 8 and Figure S13). The group 2 members IFN-b and IFN-c, fell sister to each other (UFBOOT = 95%), while within group 1, IFN-a and IFN-h form a sister group (UFBOOT = 59%), with IFN-d (UFBOOT = 45%) and IFN-e (UFBOOT = 77%) forming successive sister groups (Figure 8 and Figure S13). Thus, our phylogenetic analyses reject the monophyly of group 2 (two pairs of conserved cysteines), due to the independent origins of IFN-f, but not of group 1 (one pair of conserved cysteines) ray-finned fish IFN1s.

DISCUSSION
The origins and evolutionary relationships between, and within, interferon subtypes have proven difficult to resolve. Here, with greatly increased taxon sampling and careful application of alignment and phylogenetic methodology, we overhaul our current understanding of the origins and relationships of the three IFN classes. Our findings also provide a significant step forward compared to previous work in understanding the mode and tempo of intra-class IFN evolution.
A notable study finding was our identification of a cartilaginous fish IFN3 gene, revealing that both IFN3 ligands and receptors existed in the jawed vertebrate ancestor, helping to resolve the deep relationships within the class II α-helical cytokines. We found that the four major lineages of this gene superfamily (i.e., IFN1, IFN2, IFN3, and the IL-10 family) diverged by multiple gene duplications [or genome duplication (15)] in quick succession in the ancestor of jawed vertebrates. We also revealed that the antiviral interferons, IFN1 and IFN3, are likely sister groups, with IFN2 being sister to the IL-10 family, similar to the model proposed by Siupka et al. (15). These results reject both of the other proposed hypotheses of IFN3 origins; (i) that tetrapod IFN3 genes evolved from IFN1s (8,10,23) and (ii) that IFN3 is a member of the IL-10 family (which is based on structural homology) (14,21,22). Structural similarity between IFN3 and the IL-10 family can be explained if these features were ancestral within the class II α-helical cytokines and secondarily lost in the IFN1 and IFN2 lineages. Importantly, unraveling the early evolution of class II α-helical cytokines also allowed us to objectively choose the best outgroups to test ingroup relationships for each of the IFN classes for the first time. This, along with other improvements in phylogenetic approach, made it possible for us to resolve some of the discrepancies noted in previous studies.
Our findings corroborate the conserved nature of IFN2 genes (which are not predominantly antiviral interferons) compared to IFN1 and IFN3 (24,27,104). By incorporating the closest outgroup, including cartilaginous fish IFN-γ (18), and better accounting for insertions and deletions at the alignment stage (72), we found strong support for teleost-specific origins of IFNγ-rel by tandem duplication as proposed previously (19). Thus we can now reject the possibility that this represents an ancestral jawed vertebrate gene that was lost in other groups (24)(25)(26). Applying a similar approach to IFN3 evolution, we were able to delineate the evolution of the major IFN3 gene lineages found in humans for the first time; with the IL-28/29 ancestor diverging from IFN-λ4 in the amniote ancestor, and the IL-28 and IL-29 lineages splitting in the ancestor of placental mammals.
Our results confirm that inferring the evolutionary relationships between IFN1 family members is difficult. IFN1 phylogeny is highly sensitive to several confounding factors, including model inadequacy, distant outgroups, and limited taxon sampling. The short length and rapid evolution of IFN1s may also have driven stochastic errors and resulted in some weakly supported branches in our phylogenetic trees. Importantly, we observed consistency in our analyses that were designed to minimize systematic error (i.e., applying best-fit models and outgroups, and exclusion of compositionally biased sequences), both of which are factors that may be indicative of accuracy, even in the face of weak support (105,106). By accounting for phylogenetic error, and considering consistency across our datasets, we reconstructed a strongly supported scenario of IFN1 evolution where several IFN1 genes existed in the jawed vertebrate ancestor. These genes subsequently underwent extensive lineage-specific gene duplication and loss events. Central to this finding is our unprecedented taxon sampling, which allowed us to identify ancestral jawed vertebrate genes that have become very taxonomically confined due to multiple loss events. Our data imply that while IFN1s often undergo lineage-specific expansions, they can also be lost many times in parallel, generating extreme cases of "elusive" genes (i.e., genes which are difficult to detect because of recurrent loss or biases in generating assembled genomes) (107) and hidden paralogy (i.e., where differential loss results in paralogs presenting as orthologs) (108,109). A key example of this is the discovery of intron-containing cartilaginous fish orthologs of intronless tetrapod IFN1s, which revealed that intron-containing IFN1s of ray-finned fishes are paralogous, rather than orthologous, to the intronless tetrapod IFN1s. This means that the retrotransposition event giving rise to intronless tetrapod IFN1s may have occurred as early as in the ancestor of bony fishes (indicating loss of intronless IFN1s from ray-finned fishes and coelacanth), or as late as in the most recent common ancestor of extant tetrapods (indicating loss of intron-containing IFN1s from ray-finned fishes and coelacanth). Either way, this lineage, which is remarkably expanded in amniotes, has been lost from teleosts and coelacanth. Together these findings imply that IFN1 molecules, like some other immune genes, evolve via a rapid birth-death evolutionary process, and have done so at least since the jawed vertebrate ancestor (100,101,110). This is consistent with a scenario where IFN1 genes have maintained their antiviral function for over 450 million years by evolving rapidly, in terms of both substitutions and gene gain and loss, due to the host pathogen arms race with viruses.
Analyses focused on the evolution of ray-finned fish IFN1s revealed that their group 1 (one conserved cysteine pair), but not group 2 (two conserved cysteine pairs), interferons are monophyletic. Our findings suggest that group 2 should be split into two groups. The first consisting of IFN-b and IFNc (together these form the sister group to group 1), for which we suggest the group 2 name be retained. And the second, consisting only of ray-finned fish IFN-f (although IFN-f appears to be present in at least amphibians and cartilaginous fishes also), which we propose be referred to as group 3. Interestingly, group 1 and group 2 IFN1s use different interferon receptors in zebrafish (11,111), however zebrafish lack IFN-f, and as such it may be that IFN-f (now group 3) may have a different receptor to both group 1 and group 2. If this proved to be the case, analyses of receptor use may also help verify the assignment of amphibian and cartilaginous fish IFN-f. Importantly, although ray-finned fish group 1 and group 2 IFN1s are sister to each other, and seem to be derived from an ancestral jawed vertebrate IFN1 that has been lost in all other species, our results suggest that the ancestor of both groups possessed two conserved cysteine pairs. Based on the presence of the two conserved cysteine pairs across the IFN1 CHOM and EXT trees, our results are also consistent with the ancestral IFN1 possessing two disulphide bridges and four introns (9,19,32).
Similarly focusing on amniote IFN1 evolution we found that several intronless IFN1 genes existed in the tetrapod ancestor, with extensive IFN1 repertoires present in extant reptiles. In fact, as more ancestral tetrapod IFN1s appear to have been retained in reptiles, they evidently have even greater IFN1 diversity than mammals. Our analyses incorporated a greater breadth of mammals and reptiles than previous studies, including aquatic and/or semi-aquatic lineages, and had a more appropriate outgroup, but do not support one-to-one orthology of any mammalian or reptile IFN1s. This confirms non-orthology between human and chicken IFN-α and IFN-β (29,112), while rejecting orthology of chicken and mammal IFN-κ (31).
Emergence of intronless interferons is more common in the IFN1 and IFN3 families than previously thought, consistent with intronless interferons bestowing an evolutionary advantage over those harboring introns (39). Our results suggest that both of the models (19,23,35,40) put forth previously for the origins of amphibian intronless IFN1s are correct, with some emerging multiple times independently within amphibians, and others resulting from the same event that gave rise to amniote IFN1s. Strikingly, we also found that intronless amphibian IFN3s have emerged at least twice and independently from those of mammals on both occasions. Interestingly, amphibians also possess by far the most diverse set of IFN1s, including those which form part of the intronless tetrapod IFN1 group, the IFN-f group, and those in the sister group to all other IFN1s. Given this highly diverse repertoire of antiviral IFN1s and propensity for retrotransposition (or at least gross loss of introns), it is tempting to speculate a link to their morphology (e.g., permeable skin involved in terrestrial cutaneous respiration) or developmental life-history (e.g., aquatic tadpoles undergo metamorphosis to become terrestrial adults), especially as unique interferon responses have been observed between their distinctive stages of life (112,113).
Lastly, our study indicates that a new nomenclature system is required to describe IFN1s to avoid relying on awkward (as applied here) or inaccurate descriptions. We have not attempted to formulate one here, as it is likely to be a substantial undertaking and will require input and agreement from several parties.

DATA AVAILABILITY
All datasets generated and analyzed for the study are included in the manuscript and the Supplementary Files.

AUTHOR CONTRIBUTIONS
AR, JZ, and HD conceived the study. AR performed sequence similarity searches, designed and performed phylogenetic analyses, and drafted the manuscript and figures. JZ and HD performed IFN1 searches for cartilaginous fishes. All authors contributed to and approved study design and the final manuscript.

ACKNOWLEDGMENTS
PhyloBayes analyses were performed using the University of Aberdeen's Maxwell high performance computing cluster. AR was supported by a University of Aberdeen Center for Genome-Enabled Biology and Medicine Ph.D. studentship. DM received support from BBSRC Institutional Strategic Programme funding (grant number: BBS/E/D/20002172). Silhouettes in Figure 7 were obtained from http://phylopic.org; all of which are public domain except for the anole lizard silhouette, which was created by Ghedo and T. Michael Keesey (license: https://creativecommons.org/ licenses/by-sa/3.0/).