MHC in newts illuminates the evolutionary dynamics of complex regions in giant genomes

Comparative genomics of the MHC region(s) in newts

We analyzed four publicly available newt genomes (Lissotriton helveticus, L. vulgaris, Pleurodeles waltl, Triturus cristatus) (Fig. 1A), all assembled at the chromosome level. In each genome, there is a continuous MHC region assigned to a specific chromosome. Because the initial mapping of MHC transcripts to genomes indicated that automatic annotation of MHC genes was, in many cases, not accurate, MHC-I and MHC-II genes were annotated manually (Supplemental Data). The following description incorporates key observations informed by expression data and phylogenetic analyses, which are presented in detail later.

Figure 1.

Newt phylogeny and a schematic representation of MHC genes and proteins. (A) Time-calibrated phylogeny of the studied species (Stewart and Wiens 2025); colored dots next to the species name indicate the availability of a reference genome. (B) Schematic representation of the structure of MHC molecules and their genes. MHC-I molecules are composed of a heavy chain with three extracellular domains (α1, α2, and α3), paired with a nonvariable β2-microglobulin (β2m) molecule, which is encoded by a gene outside the MHC region. MHC-II molecules are heterodimers, each consisting of an α and a β chain. Each chain contributes two extracellular domains: α1 and α2 for the α chain, and β1 and β2 for the β chain. The peptide-binding groove in MHC-I is formed by the α1 and α2 domains, whereas in MHC-II it is formed by the α1 and β1 domains. All MHC molecules also include a transmembrane domain and a cytosolic tail. There is a direct correspondence between structural domains and exons, as illustrated on the right-hand side. The first exon (gray box) encodes a leader peptide; each extracellular domain is encoded by a separate exon (indicated as “ex” with the exon number inside the box). The remaining regions—including the transmembrane domain and cytosolic tail—are encoded by a variable number of exons (black boxes), typically ranging from one to four, depending on the molecule type.

Open in new tabDownload PowerPointLink to figure

There are both considerable similarities and striking differences in the organization of the MHC region between the newt genomes (Figs. 2, 3). A core MHC region, hereafter defined as spanning from the first MHC-II to the last MHC-I in the vicinity of APGs or the last APG, is on an orthologous chromosome in all four species, although, because chromosome nomenclature is not standardized across assemblies, it is designated Chr 5 in L. helveticus and L. vulgaris, Chr 1 in T. cristatus, and Chr 6 in P. waltl (Fig. 2). In Lissotriton and Pleurodeles the MHC region is located at the terminal part of the chromosome reaching ∼ 0.5 (L. vulgaris)–9 (P. waltl) Mb from the end (Fig. 2). In T. cristatus, the MHC region is at the distance of >120 Mb from the end of the 2.5-Gb chromosome. Notably, the orientation of the region differs between assemblies, with the extended class II region located terminally/distally in L. vulgaris and T. cristatus and internally/proximally in L. helveticus and P. waltl. The availability of linkage maps of L. vulgaris and T. cristatus (France et al. 2025) allowed us to establish that the MHC region is actually located at the opposite ends of linkage group (LG) 1 in these two species, and a comparison with the P. waltl genome suggests that the region was translocated into a new location in the T. cristatus lineage. The length of the core MHC region varies from 5.5 Mb in P. waltl to 16.7 Mb in L. vulgaris.

Figure 2.

The MHC region in amphibian genome assemblies and on linkage maps. The location and length of the MHC region comprising APG, MHC-II, extended MHC-II gene, and all annotated MHC-I genes located within 100 Mb. MHC genes in the four newt species were annotated manually in this study; manual annotation of MHC in A. mexicanum is from Migalska et al. (2025), and MHC genes in X. laevis are from the publicly available annotation of the v. 10.1 genome assembly. Note that, although in all newt genomes, the MHC region is on the same chromosome, naming and coordinates follow the original genome assemblies; numbers at the ends of chromosomes indicate coordinates (in Mb) in the original assemblies. Comparative linkage maps of L. vulgaris and T. cristatus (France et al. 2025) indicate that the MHC region is located on the opposite ends of LG1 in these species—cM 101 and 15, respectively.

Open in new tabDownload PowerPointLink to figure

Figure 3.

Genes in the MHC region(s). Upper left: Genes in the genomic regions corresponding to those in Figure 2; APG, MHC-II, extended class II, and framework genes are indicated as colored blocks; MHC-I genes are shown individually; the width of the rectangles indicates gene lengths; filled rectangles below MHC-I genes show the maximum relative expression. Coordinates were rearranged so that the start/end of COL11A2, the extended class II gene adjacent in genomes to class II, is given coordinate 0. Pseudogenes are marked with asterisks. Upper right: MHC-I/MHC-like (if any) genes located outside the main MHC region. Lower part: The zoomed in “core” MHC region indicated with dashed rectangles in the upper plot. MHC-I, MHC-II and APG are marked and labeled individually, with the width of rectangles corresponding to gene lengths. Maximum relative expression: the maximum fraction of total MHC-I expression (FPKM) attributable to a sequence across all RNA-seq libraries from postmetamorphic individuals examined within a species (Supplemental Tables S3, S4).

Open in new tabDownload PowerPointLink to figure

In all species, multiple MHC-I genes (from three in L. vulgaris to eight in P. waltl), at least some of them highly expressed, are located between the MHC-II region and APGs (Fig. 3). In L. helveticus, most APGs are duplicated with four TAP2, three PSMB8, and two PSMB9 copies, whereas in three other species, each APG is a single copy gene. In both Lissotriton genomes, some highly expressed MHC-I genes are interspersed among APGs. Newt genomes contain also two divergent lineages of MHC-I sequences, which we name MHC-I-like1 and MHC-I-like2 (Fig. 4A; Supplemental Table S1). In both genomes where we found it (P. waltl and L. vulgaris), MHC-I-like1 gene is close to TAP1, and MHC-I-like1 pseudogene in the T. cristatus genome is also in this location. A single or two MHC-I-like2 genes are located between MHC-II and APGs in Lissotriton and Triturus but on a separate chromosome in Pleurodeles.

Figure 4.

Phylogeny of newt MHC in the context of MHC sequences from other tetrapod taxa. The RAxML-NG maximum likelihood trees were constructed from protein sequences under the JTT+G4 amino-acid substitution model. The trees were rooted with zebrafish sequences (not shown). Support values for clades with the minimum bootstrap support of 70% (100 bootstrap replicates) are shown. The tree contains sequences from the G+T data set (as described in Methods) obtained from Pacific Biosciences (PacBio) Iso-Seq, de novo assembly of RNA-seq and protein sequences predicted for manually annotated MHC genes in genome assemblies of L. helveticus, L. vulgaris, P. waltl, and T. cristatus. Newt sequences are color-coded according to the species and species phylogeny is in Figure 1A. Additional salamanders included in the tree are the giant salamander (Andrias: GenBank AGY55962.1, AGY55988.1, AGY56015.1) and the axolotl (Ambystoma: sequences reported by Migalska et al. 2025). The trees include representatives of other major tetrapod groups, taken from GenBank or Ensembl (accessions starting with ENS): anurans (Xenopus: AAA16064.1, AAA16359.1, NP_001090513.1, BAA02842.1), birds (Gallus: ENSGALG00010003817, NP_001231990.1, ENSGALG00010003022), caecilians (Geotrypetes: XP_033779966.1, XP_033779817.1, XP_033779862.1), lepidosaurs (Sphenodon: ENSSPUG00000003482, ENSSPUG00000017679, ENSSPUG00000012680), and mammals (Homo: ENSG00000206503, ENSG00000204592, ENSG00000277263, ENSG00000196126). (A) MHC-I—a single (MHC-I-like1 and MHC-I-like2) or three (MHC-I) sequences randomly picked per newt species are included; for human and Xenopus, both MHC-Ia and MHC-Ib were included. (B) MHC-IIA, and (C) MHC-IIB—three sequences (if available) randomly picked per newt species were included.

Open in new tabDownload PowerPointLink to figure

In all species, some MHC-I genes are located outside of the core MHC region. Notably, MHC-I genes may be located on the same chromosome but >40 Mb (L. helveticus) or even >100 Mb (T. cristatus) from the core MHC region (Figs. 2, 3). In three species, there are some MHC-I genes at the opposite end of the MHC chromosome (L. vulgaris) or on other chromosomes (P. waltl and T. cristatus) (Fig. 3). In T. cristatus, there are four highly expressed MHC-I genes that occupy ∼3 Mb on Chr 5, not accompanied by any other genes usually found in the MHC region (Fig. 3). In the absence of another Triturus assembly or cytogenetic FISH evidence, we cannot completely rule out the possibility that the apparent translocation is an assembly error. However, genes flanking the additional MHC region from both sides map to a continuous genomic region in other newt assemblies, rendering misassembly unlikely.

In all four genomes, MHC-II genes are in the same order (Fig. 3): MHCIIB-IIA-BRD2-DMA-DMB and span from <1 Mb (P. waltl) to >3 Mb (L. vulgaris). In the L. vulgaris genome, there are two MHC-IIA and three MHC-IIB genes. The extended class II region is immediately adjacent to MHC-II in all species and spans from ∼ 5 Mb in Pleurodeles to >10 Mb in T. cristatus (Fig. 3).

Finally, we examined the immediate genomic context of the core MHC region, to see whether gene categories neighboring the MHC class I and II regions in mammals can be found in this group. Several framework genes (see Methods) were found in the vicinity of MHC in some species (Fig. 3). These are non-MHC genes interspersed in the eutherian MHC class I region (Amadou 1999), which are useful reference points for comparative analysis in mammals but not necessarily beyond (Belov et al. 2006). Additionally, some MHC class III genes were found within or close to the core MHC genes in T. cristatus, although we note that precise localization of Class III is beyond the scope of the present article.

Diversity of MHC sequences in newts

The MHC gene family frequently undergoes rapid expansions and contractions, leading to variation in gene content and blurred orthology among closely related species and even between haplotypes within species. As a result, a single reference genome per species is unlikely to capture the full extent of MHC diversity. Comprehensive sampling—both across individuals within species and across the phylogeny—is therefore key to understanding the evolution of MHC genes. To assess the diversity of MHC sequences in newts, we complemented the sequences annotated in genomes (hereafter referred to as data set G) with full-length coding sequences obtained from transcriptomes of, altogether, seven species (Supplemental Tables S2–S4; Fig. 1A). The coding sequences were obtained mostly through Pacific Biosciences (PacBio) Iso-Seq, supplemented by a selection of well-supported transcripts assembled from Illumina RNA-seq data. Translated sequences were clustered within each species to reduce redundancy (details are in Methods), and the longest protein from each cluster was selected as a cluster representative. This set of combined genome- and transcriptome-derived sequences will hereafter be referred to as data set G+T (Supplemental Data).

Phylogenies of MHC-I and MHC-II from the data set G+T are shown in Figure 4 and Supplemental Figures S1 and S2. The most striking feature of the MHC-I tree is the monophyly of all newt sequences and the presence among them of two well-supported clusters separated from the rest of the tree by long branches (Fig. 4A; Supplemental Fig. S1). These clusters, MHC-I-like1 and MHC-I-like2, are characterized by large deletions in the α1 and α2 domains, respectively (Fig. 5). MHC-I-like proteins, in particular MHC-I-like2, commonly lack conserved amino acids at the anchor residues,those important for anchoring the termini of antigenic peptides, that are conserved in the classical MHC-I of most taxa (Fig. 5; Supplemental Fig. S1; Supplemental Table S5). MHC-I-like1 sequences were detected in five species (I. apestris, L. boscai, L. vulgaris, P. waltl, and T. marmoratus), in each, most likely a single gene (Supplemental Table S6). We did not detect MHC-I-like1 sequences in either transcriptome data or in the genome assembly of L. helveticus, whereas in the T. cristatus genome, there is an apparent MHC-I-like1 pseudogene, which suggests the loss of functional MHC-I-like1 in these two species. MHC-I-like2 sequences, which, as long branches in the phylogeny indicate, evolve quickly, were detected in all species, and they do not cluster according to the species phylogeny (Fig. 4A; Supplemental Fig. S1; Supplemental Table S6). The number of variants per individual indicates multiple genes in some species. The remaining sequences on the MHC-I tree did not generally form well-supported clusters, although there was a tendency to group by species or genus (Fig. 4A; Supplemental Fig. S1). Given that MHC domains have different functions (Fig. 1B), they may be subject to different evolutionary pressures, including positive/purifying selection and gene conversion. Therefore, we additionally reconstructed separate phylogenies for MHC-I α1-3 domains (Supplemental Figs. S3–S5). These trees did not depart substantially from the full-length protein phylogenies, although branch lengths, also those leading to MHC-I-like1 and MHC-I-like2, were much shorter for α3 than for the other two domains. Additionally, in α3, the MHC-I-like2 clade was not strongly supported (Supplemental Fig. S5).

Figure 5.

Protein alignments of α1–α3 domains, conserved and positively selected sites. Representative sequences of MHC-I, MHC-I-like1 (L1), and MHC-I-like2 (L2) proteins predicted from newt genomes were aligned to the Xenopus laevis allele with resolved crystal structure (Xela-UAAg, PDB: 6A2B) and representative, highly expressed MHC-I sequences of Ambystoma mexicanum (Amme-04 and Amme-05) (Migalska et al. 2025). Positively selected sites identified in newts by FUBAR are shaded in gray; conserved anchoring residues are marked with orange stars. Species abbreviations: (Lihe) L. helveticus, (Livu) L. vulgaris, (Plwa) P. waltl, (Trcr) T. cristatus.

Open in new tabDownload PowerPointLink to figure

The number of MHC-II sequences in the G+T data set was much smaller than that of MHC-I, and their phylogenies showed simpler patterns (Fig. 4B,C; Supplemental Fig. S2). In the case of nonclassical MHC-II, DMA phylogeny was poorly resolved, whereas the DMB tree generally reflected species phylogeny, with one or two sequences per species (Supplemental Fig. S2). Shorter DMA branch lengths indicate slower evolution compared to the DMB. Phylogenies of classical MHCIIA and MHCIIB showed considerable similarity to each other—there was a somewhat distinct Pleurodeles lineage and three divergent clusters, each of them containing sequences from all the remaining three genera. This pattern caused MHCIIA and MHCIIB phylogenies to depart drastically from species phylogeny. We did not find any evidence for the recently described MHC W type genes, that have both MHC-I and MHC-II features (Okamura et al. 2021), in newt genomes or transcriptomes.

The number of MHC genes and pseudogenes

All studied newts have multiple MHC-I genes (Supplemental Table S6). It is, however, not straightforward to estimate their number. Excluding MHC-I-like, there are from six (L. vulgaris) to 11 (T. cristatus) genes in the genome assemblies. The number of sequences retrieved from transcriptomes of individuals with Iso-Seq data available (these estimates include both Iso-Seq and Illumina RNA-seq G+T data set sequences in these individuals) varied from 10 in P. waltl to 15 in I. alpestris, and in all cases was larger than the number of genes detected in the assembly (Supplemental Table S6). This is not surprising as divergent alleles found in heterozygous individuals could have been assigned into separate clusters. Also, the number of MHC-I pseudogenes and gene fragments differs between species, ranging from eight in P. waltl to 23 in L. helveticus, with an overwhelming majority of pseudogenes and gene fragments within or in the vicinity of the MHC region (Fig. 3; Supplemental Table S6).

In each genome, there are single DMA and DMB genes and, except in L. vulgaris, single MHC-IIA and MHC-IIB genes (Fig. 3). We did not identify DM pseudogenes in any genome, and generally, there were none or few MHC-II pseudogenes, the only exception being 21 MHC-IIA pseudogenes in the T. cristatus genome. In contrast to MHC-I, MHC-II pseudogenes were located mostly outside the MHC region (Supplemental Table S7).

Ontogenetic and tissue expression profiles, polymorphism, and identification of putative nonclassical MHC-I genes

Ontogenetic profiles of MHC expression can shed light on the timing of adaptive immunity onset in newts. Moreover, when combined with expression pattern in multiple tissues and polymorphism data, they can help to distinguish classical MHC-Ia from nonclassical MHC-Ib genes. Overall, changes in expression of MHC-I and MHC-II were synchronized throughout development (Fig. 6). In P. waltl, there is a steady increase of MHC expression through the larval stages, whereas in the remaining six species, MHC expression in the larval period and immediately after metamorphosis is at least an order of magnitude lower than in adults (Fig. 6), with a sharp increase visible only in mature individuals. The ontogenetic expression profile of MHC-I-like genes is very different; they are expressed at much lower level, and their expression remains low throughout the larval period and metamorphosis and in adults (Fig. 6).

Figure 6.

Overall MHC expression through ontogeny. Expression in fragments per kilobase of transcript length per million reads mapped (FPKM) of all MHC-I, MHC-II (MHC-IIA+MHCIIB), and MHC-I-like1 and MHC-I-like2 sequences through the ontogeny. Points are expression values (summed over all transcripts of a given MHC category) for RNA-seq libraries (Supplemental Table S4), and bars are averages across libraries for each species and stage. Species abbreviations: (Ia) I. alpestris, (Lb) L. boscai, (Lh) L. helveticus, (Lv) L. vulgaris, (Pw) P. waltl, (Tc) T. cristatus, (Tm) T. marmoratus.

Open in new tabDownload PowerPointLink to figure

In both G and G+T data sets, only some sequences are highly expressed at any larval stage or in adults. Whereas in Triturus and perhaps also P. waltl and L. vulgaris, not more than two MHC-I genes are highly expressed in adults, in I. alpestris, L. boscai, and L. helveticus, several genes appear to have similarly high expression (Fig. 7; Supplemental Fig. S6). Moreover, many sequences highly expressed in at least one tissue/ontogenetic stage had conserved amino acids in most (eight or nine out of nine) anchor residues—one of the hallmarks of MHC-Ia molecules (Supplemental Fig. S1; Supplemental Table S5). We did not detect pronounced, several-fold differences between adult tissues in relative expression of particular MHC-I sequences, with the exception of P. waltl Plwa-03 and Plwa-04 which have much higher relative expression in skin (tailtip) and T. marmoratus Tm-13 with much higher expression in lungs than in other tissues (Supplemental Fig. S7).

Figure 7.

Expression through ontogeny of MHC-I genes annotated in the genomes (data set G). Relative expression, that is, the fraction of all RNA-seq reads in a library mapping to all MHC-I genes and pseudogenes in the assembly that mapped to a particular gene sequence, is shown. Box plots show medians, interquartile, and total ranges. Gray stripes indicate genes that show higher relative expression in larvae than in adults—those may be nonclassical (MHC-Ib) genes. Species abbreviations: (Lh) L. helveticus, (Lv) L. vulgaris, (Pw) P. waltl, (Tc) T. cristatus.

Open in new tabDownload PowerPointLink to figure

A striking pattern of ontogenetic changes in the relative expression of some MHC-I genes was detected in all seven species. During the larval period, when the overall MHC expression is low, the relative expression of some MHC-I genes was much higher than in adults (Fig. 7; Supplemental Fig. S6). Such a pattern suggests a possible functional divergence from the role played by classical MHC genes and hints at a role in larval immunity or development. The sequences with high relative expression in larvae of Triturus and Pleurodeles are easily identified visually in Figure 7 and Supplemental Figure S6, whereas the assignment of sequences to this category in other species is less straightforward. Still, such sequences do not form separate clades and occur in different parts of the phylogeny (Supplemental Fig. S8). Therefore, regardless of difficulties in defining the exact set of such sequences, it is clear that they have emerged repeatedly and are probably undergoing fast turnover on evolutionary time scales.

For the four species with assembled genomes, we estimated polymorphism of MHC-I genes by mapping variants obtained in previous amplicon-based studies (Palomar et al. 2021; Gaczorek et al. 2023) to the genome and calculating the number of variants mapping to a particular MHC-I gene (Supplemental Tables S8 and S9; Supplemental Fig. S9). Combining expression patterns with these polymorphism estimates (Fig. 8) suggests that sequences highly expressed in adults, which are also highly polymorphic, represent classical (MHC-Ia) genes. Sequences with high relative expression in larvae, lower expression in adults and/or substantial differences in expression among adult tissues, which also show low polymorphism (Fig. 8), likely represent nonclassical (MHC-Ib) genes. The status of the remaining genes, which have low relative expression in all examined stages and tissues, as well as generally low polymorphism, is unclear.

Figure 8.

Polymorphism of MHC-I genes annotated in the genomes. Bar heights show polymorphism of each gene expressed relative to the most highly polymorphic gene in the species; MHC-I-like genes are not shown because reliable polymorphism estimates were not available. Polymorphism was estimated based on mapping amplicon variants detected in population samples to genes, as described in the text. Maximum relative expression: the maximum fraction of total MHC-I expression (FPKM) attributable to a sequence across all RNA-seq libraries from postmetamorphic individuals within a species. Species abbreviations: (Lh) L. helveticus, (Lv) L. vulgaris, (Pw) P. waltl, (Tc) T. cristatus.

Open in new tabDownload PowerPointLink to figure

Overall, in all species, we see three broad categories of MHC-I genes: (i) putative classical MHC-Ia, expression of which increases in ontogeny, is very high in adults, and they exhibit high polymorphism; (ii) putative nonclassical (MHC-Ib) with low polymorphism, higher relative expression in larval stage and/or tissue-specific expression in adults; and (iii) either nonclassical or nonfunctional, with low expression throughout ontogeny and tissues and, more often than not, low polymorphism. The distribution of these three categories on phylogeny indicates fast turnover and rapid evolution/diversification (Supplemental Fig. S8).

Evolution of the MHC-I gene family in newts

The results presented above indicate a dynamic evolution of MHC-I family in newts, with multiple gene duplications and losses. To obtain a more quantitative picture of these dynamics, we performed gene tree- species tree reconciliation to infer the number of duplications and losses and identify clusters of orthologs. Reconciliation of the MHC-I gene tree with the species tree was performed separately for data sets G and G+T. Each data set produced a single event history—for G, 26 duplications and 10 losses, and for G+T, 91 duplications and 66 losses. All species have experienced lineage-specific duplications and losses, although the exact numerical results in Supplemental Table S10 should be treated with caution because some analyzed sequences may be derived from the same gene. However, we were not able to reliably reconstruct groups of orthologous sequences because, within the same event history, a large number of nodes in the reconciled tree could have been swapped without changing the overall number of events.

Gene conversion may lead to the homogenization of sequences of different members of gene families, resulting in concerted evolution that obliterates the signal of orthology. It is also an important mechanism generating new MHC variants by transferring small sequence segments between alleles and genes. To see whether concerted evolution via interlocus conversion could have contributed to the poor resolution of MHC-I phylogeny and to the tendency for MHC sequences to group by species, we tested for gene conversion, estimated the length of conversion tracts, and checked whether some domains were involved in the process more often than others. For MHC-I-like1 and -like2 which were analyzed together and with all species pooled, the signal of gene conversion was detected only for domain α3, and the tracts were shorter than in MHC-I. In MHC-I, conversion was detected in all species, with a tendency to be more common in α3, which may contribute to the low sequence divergence seen in this domain (Supplemental Fig. S5). We detected both conversion tracts contained within a single domain (i.e., within a single exon), and tracts spanning more than one exon. Because more of the latter were detected when mismatches within tracts were allowed (Supplemental Table S11) and because tracts are expected to be relatively short so that they should not span across long newt introns, we suspect some signal may represent false positives. Nevertheless, the process of gene conversion appears widespread in newt MHC-I and may lead to concerted evolution of MHC genes within species.

Selective pressures across MHC-I phylogeny

We were interested whether MHC-I and two MHC-I-like lineages experienced different selective pressures. In particular, we tested, using RELAX and FUBAR: (i) whether a change in the strength of selection occurred along the long branches separating MHC-I-like1 and MHC-I-like2 from other sequences in MHC-I phylogeny, and (ii) whether the number and location of codons under positive/purifying selection differ among the three categories of genes.

There was evidence of selection relaxation for MHC-I-like1: a highly significant (P-values < 0.001) RELAX test result for the branch leading to MHC-I-like1 cluster (k: 0.356, LRT: 401.7), and for the cluster itself (k: 0.476, LRT: 42.7), as compared to the MHC-I sequences. Accordingly, a direct, codon-based selection test (FUBAR) found only a few codons under positive or negative selection among MHC-I-like1 sequences (three and seven, respectively—posterior probability >0.9). The second divergent clade, MHC-I-like2, was associated with a slight relaxation of selection intensity in the branch leading to the clade (P-value: 0.021, k: 0.857, LRT: 5.3) but not in the cluster itself (P-value: 0.112). A direct selection test found 12 sites under positive and 50 sites under negative selection, which was considerably fewer compared to MHC-I sequences (30 and 122, respectively). Only two of the sites under positive selection overlapped between the two groups (Fig. 5). The selection tests were performed only on the three extracellular domains of MHC-I (α1, α2, α3), because the divergence of the remaining parts (i.e., the connecting, transmembrane, and cytosolic portion of the molecules) rendered the alignment unreliable. Roughly half of the residues in the MHC-I were under negative selection, compared to ∼20% in the case of the MHC-I-like2 sequences. Most of the positively selected sites in MHC-I were located in the α1 and α2 domains (14 and 10, respectively) and largely coincided with or were immediately adjacent to residues forming peptide biding pockets as identified in Xenopus (Fig. 5). MHC-I-like2 sequences had many fewer positions under positive selection, scattered through the extracellular domains (four, six and three in α1, α2, and α3, respectively) (Fig. 5).

The overall concordance between the positively selected sites identified in the sequences forming the major MHC-I clade (Fig. 4A; Supplemental Fig. S1) and the peptide binding pockets identified in Xenopus (Fig. 5) supports the classification of these sequences as the putative MHC-Ia (and perhaps some MHC-Ib). In contrast, markedly relaxed selection in MHC-I-like1 suggests a trajectory leading to pseudogenization, whereas, in the case of MHC-I-like2, change in selective pressure strength and location of codons under positive selection hints to alterations in function, perhaps specialization.

Structural modeling of MHC-I-like sequences

MHC-I-like1 and MHC-I-like2 sequences exhibit characteristic deletions within the α1 and α2 domains (Fig. 5), observed in transcripts and confirmed in the available genomes. MHC-I-like1 molecules lack a signal peptide and have a truncated α1 domain (missing 23 N-terminal amino acids). Lack of a signal peptide likely precludes transport to the endoplasmic reticulum, where MHC-I is folded and loaded with antigens. Thus, this modification alone suggests departure from the classical function. In contrast, MHC-I-like2 molecules possess a short N-terminal sequence resembling a signal peptide, an intact α1 domain, and substantial deletions within the α2 domain. To evaluate the structural impact of these alterations on the peptide-binding domains, we performed in silico modeling using AlphaFold 3. Presumably classical MHC-I proteins, included as references, were modeled with high confidence, exhibiting predicted local distance difference test (pLDDT) scores >90 in the α1–α2 domains. In comparison, the predicted α1–α2 domains of MHC-I-like molecules showed lower confidence, with overall pLDDT scores ranging from >70 to <50 in the most uncertain regions. Most of the MHC-I-like1 models preserved the overall MHC-I fold, including the two antiparallel α-helices, but lacked the first β-strand and loop of the α1 domain (Supplemental Fig. S10A), although one variant showed a closed conformation without open pockets (L. vulgaris conformation 2) (Supplemental Fig. S10A). Structural homology searches identified MHC and MHC-like molecules from other vertebrates as the closest matches, but none of the molecules contained deletions in the discussed region, and thus the role of the deletion remains unknown. We speculate that a molecule with such an extensive deletion in the core would not fold altogether, but this speculation remains to be experimentally tested. In MHC-I-like2 molecules, deletions within the α2 domain affect regions of the peptide-biding groove that, in classical MHC-Ia molecules, anchor the C-termini of bound peptides. Two β-strands forming the groove floor were either missing or significantly shortened, and a substantial portion of the α2 helix (including the conserved R146 and W147 residues) was also absent (Supplemental Fig. S10B). Structural homology searches again primarily returned MHC and MHC-like molecules; these molecules however, lacked the discussed deletion. A notable exception was the viral MHC-I homolog m144 (from murine cytomegalovirus, PDB: 1U58), which shares similar deletions in the α2 domain. Unlike m144, which lacks a well-defined groove, the predicted MHC-I-like2 structures displayed seemingly an open conformation with a spacious, predominantly lipophilic pocket—this may, however, reflect the bias of AlphaFold towards classical MHC structures. Functionally, m144 is a molecular decoy mimicking host MHC-I and expressed by the virus to prolong its survival by preventing natural killer cell activation (Natarajan et al. 2006). Although this does not offer clues to the function of MHC-I-like2 in vertebrate hosts, it does support the plausibility of its folding pattern.

Genome assemblies

We described MHC region(s) in chromosome-scale genome assemblies of P. waltl (obtained from the GenBank database [https://www.ncbi.nlm.nih.gov/genbank/] under accession number GCA_031143425.1) (Brown et al. 2025) and three newt species sequenced and assembled by the Darwin Tree of Life Project: L. helveticus (GCA_964261635.1), L. vulgaris (GCA_964263255.1), and T. cristatus (GCA_964204655.1). The analyses were based on the primary haplotypes under the indicated accession numbers.

RNA-seq

We analyzed transcriptomes of the following newt species: I. alpestris, L. boscai, L. helveticus, L. vulgaris, P. waltl, T. cristatus, and T. marmoratus (Fig. 1A). Separate libraries were prepared for each tissue/body part of each individual (Supplemental Table S4). We did not analyze thymus specifically because, given differences in larval size across species and developmental stages, it was often too small to yield sufficient RNA for library preparation with the chosen methodology. To maintain consistency across samples and taxa, we therefore collected larger tissue fragments. Thymus-derived expression data would be highly informative and represent an important avenue for future studies. Details about samples and laboratory procedures are in Supplemental Tables S2 and S4 and Supplemental Methods.

Full-length transcript sequencing

For P. waltl, we used the available PacBio Iso-Seq data obtained from the spleen RNA of the individual used for genome assembly (Brown et al. 2025). For the remaining six species, Iso-Seq libraries from intestine (a tissue where diverse immune cells are abundant) RNA of a single individual per species were prepared and sequenced on the PacBio Revio platform by Macrogen. Note that, for each of these samples, short-read Illumina data were also generated. Details are in Supplemental Methods.

Identification of MHC transcripts

To identify MHC sequences across newt species, we used axolotl and P. waltl MHC proteins as queries to search Iso-Seq transcriptomes. Candidate transcripts were filtered and clustered to reduce redundancy, and representative sequences were aligned and curated. To capture additional diversity that might be absent from Iso-Seq data, we supplemented these analyses with de novo assemblies of short-read RNA-seq data. Although some of the short-read–derived sequences may be artifactual, they provide complementary evidence for additional MHC diversity. Details are in Supplemental Methods.

The genomic organization of the MHC region

The genomic location of the “core” MHC region and any disparate genome fragments that may contain MHC genes was established by BLASTN searches with cluster representatives obtained as described above. To annotate genes both within these regions and in their vicinity, we adopted two routes; protein-coding genes other than MHC-I and MHC-II were automatically annotated, whereas MHC genes were annotated manually. Details are in Supplemental Methods.

To ease comparative analysis, we classified all the annotated functional genes following categories used by Belov et al. (2006): (i) MHC Class I (putative classical or nonclassical); (ii) MHC Class II (including classical MHC-II, DM, and a non-MHC, yet closely linked with Class II region BRD2 gene); (iii) Antigen Processing Genes (which include PSMB8, PSMB9, TAP1, and TAP2); (iv) Extended class II (from COL11A2 to ZBTB22, if present, including TAPBP); (v) MHC Class III (if present); and (vi) Framework region genes (conserved, non-MHC genes known to be interspersed in the eutherian MHC class I region [Amadou 1999], found in the vicinity of the core MHC region in nonplacental mammals [Belov et al. 2006], and more dispersed in nonmammalian genomes). These categories are well suited to describe noneutherian MHC organization while maintaining clear correspondence to the gold standard of the human MHC nomenclature. For clarity of visualization, we show only synteny-informative Framework genes (i.e., present in more than one species) and omit genes that are not mapping to human MHC.

MHC-I polymorphism

One of the key features distinguishing classical MHC genes from nonclassical is their high polymorphism. To assess polymorphism of the MHC-I genes annotated in the genomes, we used exon 2 amplicon sequencing data from previous studies (Palomar et al. 2021; Gaczorek et al. 2023). Polymorphism was measured as the number of alleles mapped to each MHC-I gene, with corrections applied for alleles mapping to multiple genes, and standardized relative to the most polymorphic gene. Details are in Supplemental Methods. This approach is approximate, as alleles may not always be most similar to the reference sequence of their genes of origin due to the high divergence between alleles and interlocus recombination, the genes of origin may be missing from the reference assembly, and alleles may map equally well to multiple genes. Still, even such an approximate measure should provide meaningful information about polymorphism.

Phylogenies of MHC proteins

Clustering of MHC proteins, following the procedure described above, was repeated to include also the proteins predicted in genomes. When available, the longest proteins predicted in the genome were used as cluster representatives. The sequences formed the G+T dataset which was subject to phylogenetic analyses. Details are in Supplemental Methods.

MHC-I anchor residues

MHC-I anchor residues are defined as amino acids at nine residues important for anchoring the termini of antigenic peptides that are conserved in classical MHC-I of most taxa (Kaufman et al. 1994; Sammut et al. 1999; Almeida et al. 2021). The number of residues with conserved amino acids may be helpful in distinguishing classical and nonclassical MHC-I sequences, and we calculated it for our sequences.

MHC-I and MHC-II expression across tissues and developmental stages

MHC-I and MHC-II expression across tissues and developmental stages was quantified by mapping RNA-seq data to genome assemblies (where available) and to representative full-length MHC coding sequences (details are in Supplemental Methods). Expression was calculated as fragments per kilobase of transcript length per million reads mapped (FPKM), and relative expression was defined as the fraction of total MHC-I expression attributable to a given sequence in each library. For cases requiring a single expression value (e.g., for visualization), we used the maximum relative expression across all RNA-seq libraries from postmetamorphic individuals of a species.

Gene tree–species tree reconciliation

Gene tree-species tree reconciliation was performed in Notung (Durand et al. 2006) to infer the number of duplications and gene losses and identify clusters of orthologs in MHC-I. Details are in Supplemental Methods.

Gene conversion

Gene conversion was tested using coding DNA sequences of the G+T data set after removing sequences assembled by Trinity (Grabherr et al. 2011) from short RNA-seq reads, as some of them could have been computationally generated chimeras that may produce a false recombination signal. Details are in Supplemental Methods.

Selective pressures acting on MHC-I and MHC-I-like genes

Selective pressures acting on MHC-I and MHC-I-like genes were investigated using HyPhy (Kosakovsky Pond et al. 2020), focusing on whether selection pressures differed along major phylogenetic branches and between clades, and on the identification of codons under positive or purifying selection. Details are in Supplemental Methods.

Structural modeling

Structural modeling of MHC-I-like molecules and representative classical MHC-I molecules from each species, together with β2m, was performed using AlphaFold3 (Abramson et al. 2024). The resulting models were inspected for structural features and compared to known proteins to assess homology. Details are in Supplemental Methods.