We primarily relied on publicly available genomic resources to retrieve MHC class IIA and IIB sequences across different clades of non-passerine birds (GenBank numbers and sequences provided in Appendix 1 , 2 , respectively). Because the ancient MHC-IIB divergence was originally resolved in owls (24, 25), we generated a consensus MHC-IIB query from the orthologous DAB1 and DAB2 exon 3 sequences (282 bp) retrieved from a true owl species, the tawny owl Strix aluco (GenBank nos. KJ162540 and KJ162542, respectively). Because no owl MHC-IIA sequences were available in the Nucleotide Database of the National Center for Biotechnology Information (NCBI, Bethesda, MD, USA), we generated a consensus MHC-IIA query from the available MHC-IIA exon 3 sequences retrieved from a procellariform species, the Leach’s storm petrel Hydrobates leucorhous (GenBank nos. MN061391-MN061408). Both queries were BLASTed against 370 genomes of non-passerine birds available in the Genome NCBI Database (as accessed in November 2022), including a considerable number of long-read assemblies, which greatly improves the resolution of genomic features (32). Due to long divergence times between different non-passerine clades, we used the blastn algorithm, which was designed to search for short queries of moderate similarity in a cross-species framework. We only retained contigs containing complete or nearly complete (>90% of query coverage) exon 3 sequences. We also searched the Nucleotide NCBI Database for avian MHC-IIA and MHC-IIB exon 3 sequences obtained via targeted sequencing and these searches primarily focused on avian orders with low availability of genomic data. We extracted exon 3 sequences from all the contigs and aligned them by order in Geneious v.10.0.5 software (Biomatters Ltd., Auckland, New Zealand). Alignments from all palaeognath orders were merged into a single alignment due to low sample size within each order, while alignments from neognath orders were kept separate. All alignments were trimmed to full codons and we also removed the first (5’-end) and the last (3’-end) codons, because they were missing in some alignments, retaining a uniform 276 nt length of exon 3 across all orders. We removed highly similar (>98% of pairwise nucleotide similarity) sequences within each species (ca. 40% of all originally aligned sequences) and we only retained alignments with at least three MHC-IIA and MHC-IIB sequences originating from at least two species per order for the downstream analyses. Based on this selection criteria, we retained a single combined alignment for Palaeognathae (n = 4 orders) and 17 alignments for different orders of non-passerine Neognathae birds. Our final dataset comprised 184 MHC-IIA and 268 MHC-IIB exon 3 sequences originating from 235 species (156 and 208 species for MHC-IIA and IIB, respectively) and representing 71 avian families. The wider phylogenetic coverage of MHC-IIB than MHC-IIA sequences was primarily attributed to their better availability in the Nucleotide NCBI Database (MHC-IIA are largely understudied in birds), rather than to better retrieval rate from genomic resources.

To resolve clusters of orthologous MHC-IIA sequences within our alignments we used Bayesian tree-building algorithms, as implemented in MrBayes v.3.2.6 (33). Since orthology of MHC-IIB sequences was originally resolved in owls (24), we started phylogenetic analysis of MHC-IIA exon 3 sequences from this avian order. Prior to phylogenetic analysis, selection of an appropriate model of nucleotide substitution was performed in Mega v.6.0 software (34) using maximum likelihood and the lowest Bayesian Information Criteria (BIC). The Chinese giant salamander Andrias davidianus class IIA gene was used as outgroup (Genbank no. KF611869). Phylogenetic analysis was conducted using two independent chains, each consisting of one million iterations, burin-in length of 250 000 steps/iterations, and sampling every 500^th tree, yielding an expected sample size of 1500 retained trees per chain. The analyses were run using the Molecular Clock with Uniform Branch Lengths models and the default uninformative flat priors. We also tested models with Unconstrained and Exponential Branch Lengths across selected avian orders, but they provided qualitatively similar tree topologies (results not shown). The average effective sample size for model parameters was 1453 ± 17 [SE]. Convergence of independent chains was confirmed by the average standard deviation of split frequencies (ASDSF) approaching zero (35) and the potential scale reduction factors (PSRF) approaching one (36). Extended majority-rule consensus topologies were determined using Geneiuos v.10.0.5.

Since phylogenetic analysis of MHC-IIA sequences in owls provided support for the presence of two lineages of orthologous sequences (henceforth referred to as DAA1 and DAA2; see results for details), we have used the same protocols as described above to resolve MHC-IIA orthology in all other avian orders. Blakiston’s fish owl Bubo blakistoni sequences (contig GenBank nos. BJCB01040766 and BJCB01033660, respectively; genome assembly GenBank no. GCA_004320225) were added to each alignment as reference DAA1 and DAA2 sequences. Since all amino acid differences between consensus owl DAA1 and DAA2 sequences concentrated in the upstream (5’-end) 90 nt region of exon 3, we expected this part of the sequence to yield the highest power to resolve orthology of MHC-IIA genes in other bird orders. In fact, the same pattern was described for MHC-IIB exon 3, where amino acid differences between DAB1 and DAB2 sequences also concentrated in the upstream region (24, 26). Taking all this into account, all phylogenetic analysis was based on this short upstream fragment (90 nt) of MHC-IIA exon 3. An assignment of MHC-IIA sequences to either DAA1 or DAA2 gene lineage was based on their monophyletic clustering with reference DAA1 and DAA2 sequences. Since nucleotide sequences provided poor cluster support (<0.50) in Psittaciformes, we used amino acid sequences to resolve orthology in this clade.

The evolutionary history of both orthologous gene lineages was reconstructed along the bird tree topology developed by Prum et al. (37) and visualized using the Interactive Tree Of Life (iTOL) v.5.0 webserver (38). We also tested whether phylogenetic resolution of DAA1 and DAA2 was retained across the remaining part of MHC class IIA exon 3. For this purpose, we conducted phylogenetic analysis of the central and downstream region (186 nt) of exon 3 in three selected avian orders with extensive coverage (at least ten available sequences per order), all showing evidence for the presence of DAA1 and DAA2 in the analysis of the upstream region of exon 3 (Charadriiformes, Procellariiformes, and Coraciiformes).

Although the history of the ancient duplication at the MHC-IIB has already been tracked throughout the evolution of birds (26), we aimed to assign all the retrieved MHC-IIB sequences to either DAB1 or DAB2. Similarly to MHC-IIA, we focused on the upstream (5’-end) 90 nt region of exon 3 in the phylogenetic inference. We used Bayesian phylogenetic analysis algorithms, as described above. Tawny owl sequences (GenBank nos. KJ162540 and KJ162542) were added to each alignment as DAB1/DAB2 reference and the Chinese giant salamander was used as outgroup (Genbank no. KF723002).

Our phylogenetic analyses were only used to identify the occurrence (presence or absence) of gene lineages at both MHC-IIA and MHC-IIB (DAA1/DAB1 and DAA2/DAB2) and we did not aim to quantify gene copy number within each lineage. The occurrence of gene lineages was inferred on the level of bird orders (Neognathae) or across orders (Palaeognathae). We did not aim to characterize gene lineage occurrence at the species level, as technical artefacts (sequences missing in low-quality genome assemblies) could be easily confused with a true absence (loss) of either gene lineage.

From the contigs containing either MHC-IIA or MHC-IIB exon 3 sequences, for which orthology was successfully resolved via phylogenetic inference, we also retrieved sequences of exon 2 encoding the peptide-binding α1 and β1 domains of the MHC-II molecule. In total, we retrieved 157 MHC-IIA and 215 MHC-IIB full length (255 and 270 nt, respectively) sequences of exon 2 across the entire dataset. All retrieved sequences encoded putatively functional proteins (lacked stop codons or frame-shift mutations).

All retrieved MHC-IIA and IIB exon 2 and 3 sequences were used to quantify and compare sequence diversity and selection between membrane-proximal (α2 and β2) and peptide-binding (α1 and β1) domains of DAA1/DAA2 and DAB1/DAB2. For this purpose, we prepared combined alignments across all avian orders for each exon (separately for each gene and gene lineage) and we first used DnaSP v6.10.03 (39) to quantify sequence diversity (as measured with nucleotide diversity and the total number of mutations). Second, we removed putative recombinant sequences from the alignments, since under high recombination rate no unique tree topology can describe evolutionary history of entire sequences (40), which may lead to biases in selection estimates. We used RDP v.4.80 software (41) to assess recombination signal using seven different algorithms (all available except for the most computationally demanding LARD approach). All analyses were run using default settings and recombination events were identified if supported by at least three different approaches. An inference of recombination signal at each exon within each gene lineage was conducted separately for each avian order (except of palaeognath birds which were analysed jointly). In total, 11.3% sequences showed recombination signal (1.4-18.5% sequences per exon per gene lineage) and they were all removed from the final alignments prior to the analyses of selection patterns. However, because inferring recombination between species (within orders) may produce false positives (e.g. homoplasy could yield a similar signal as recombination), we also re-ran all selection inferences using the entire dataset (no putative recombinant sequences removed). The estimates of selection (dN and dS, see below for details) for each exon from each gene lineage were highly consistent between the two datasets (average Pearson correlation coefficient r = 0.94 ± 0.02 [SE]) and, thus, we concluded that any possible removal of false positive recombinant sequences from the alignments was unlikely to introduce any major bias in the results. Our final alignments (after removal of recombinant sequences) mostly consisted of single sequences per species and only ca. 5% sequences were identified as multiple replicates within species (max. 3 sequences per species). To test if the presence of multiple sequences per species biased our selection inferences, we randomly retained a single sequence per species in each alignment and re-ran the analyses. The dN and dS estimates were highly consistent between the two approaches (average r = 0.989 ± 0.003 [SE]), thus providing support for no biases.

The patterns of selection were quantified using the relative rate of nonsynonymous nucleotide substitutions at non-silent sites (dN), which are a target of selection, to synonymous substitutions at silent sites (dS), which are presumably neutral (42). In general, sustained positive (diversifying) selection promotes changes in the protein sequence and, thus, nonsynonymous substitutions accumulate faster than synonymous substitutions (dN/dS > 1). In contrast, negative (purifying) selection suppresses protein changes, so nonsynonymous substitutions are removed and accumulate at a slower rate than synonymous substitutions (dN/dS < 1). To infer codon-specific signature of pervasive, i.e. constant across the entire tree (affecting all avian lineages) positive and negative selection we used Fast Unconstrained Bayesian AppRoximation (FUBAR) (43), as implemented in the Datamonkey v.2.0 webserver (44). This approach is based on a Markov Chain Monte Carlo (MCMC) routine, which is relatively robust against model misspecifications and leaves the distribution of selection parameters essentially unconstrained (43). Mixed Effect Model of Evolution (MEME) approach (45) was used to infer signature of episodic positive selection, which is apparent at a proportion of tree branches (i.e. affecting only some avian lineages). We also used Contrast Fixed Effects Likelihood (contrast-FEL) method (46) to identify codons showing differences in selection signature (dN/dS ratio) between the two gene lineages (across both domains of either MHC-IIA or IIB gene). All these analyses were run using default settings and input trees inferred directly from the alignments. Finally, we used branch-specific FEL approach (47) to test if positive selection was associated with the early post-duplication divergence of DAA1 and DAA2 exon 3, as previously shown for DAB1 and DAB2 (25). For this purpose, we searched for the signature of positive selection at the branch separating DAA1 from DAA2 lineages. Although this analysis was conducted for the entire exon 3, the input tree was based on the upstream 90 nt region, which provided a good resolution of DAA1 and DAA2 clusters. Codons (sites) under selection were identified based on posterior probabilities >0.95 (FUBAR) and P values <0.05 (MEME, FEL, and contrast-FEL).

The dN/dS ratios (as estimated with FUBAR) were compared between genes (MHC-IIA and IIB), gene lineages (DAA1/DAA2 and DAB1/DAB2), exons within each gene lineage (exon 2 and 3), and regions within each exon. We distinguished two regions within exon 3: i) the upstream region (90 nt), which was used in phylogenetic analysis to resolve orthology and is primarily involved in interactions between the two membrane-proximal (α2 and β2) domains of the MHC-II molecule, and ii) the remaining central and downstream region (186 nt), which is more involved in interactions with the peptide-binding domains, but also forms a conserved hydrophobic pocket interacting with CD4 T cell coreceptor (β2 domain only). Positions of codons involved in different types of interactions (between membrane-proximal domains, with peptide-binding domains, or with CD4 coreceptor) were inferred based on the crystallographic analyses of MHC-II molecules in humans and chickens Gallus gallus (48–50). Within exon 2 we also distinguished two classes of codons: i) the peptide-binding sites (PBS), which are directly involved in interactions with antigens, and ii) the non-peptide binding sites (non-PBS) comprising all the remaining codons. Positions of putative PBS codons were inferred based on the crystallographic analyses of MHC-II molecules in humans (51).

We used general linear models (GLM) to test for differences in selection signature (dN/dS ratio) between genes, gene lineages, exons, and exonic regions. Codon-specific estimates of dN and dS we log-transformed to improve normality (final skewness <1) and log dN/dS ratio was entered as a response variable in the models. First, separate GLMs were run for both exons (2 and 3) and included gene, gene lineage, and exonic region as fixed factors. We also included all possible two-way interactions between fixed factors, but they were removed from the final models if non-significant. Second, to test for the differences in selection between both exons we ran two separate GLMs for MHC-IIA and IIB genes, including exon identity as a fixed factor. All GLMs were run in Statistica v.13.0 (Tibco Software Inc., Palo Alto, CA, USA). We also used G tests to detect differences in the proportion of positively and negatively selected sites between gene lineages of MHC-IIA and IIB. Finally, we used the RELAX approach (52) available at the Datamonkey webserver to test for the relaxation or intensification of selection between gene lineages. Although this method does not explicitly test for positive selection, it uses likelihood ratio tests (LRTs) to identify shifts in the stringency of selection between two subsets of branches (herein corresponding to both gene lineages).

Phylogenetic analysis of MHC class IIA exon 3 sequences in Strigiformes provided evidence for two distinct clusters of sequences yielding very high consensus support ( Figure 1 ). Since the sequences clustered by their apparent orthology rather than by species and since more than one sequence retrieved per genome indicated gene duplication, we considered these clusters to represent two MHC class IIA gene lineages (DAA1 and DAA2).

The clustering of MHC-IIA exon 3 sequences (upstream region) in the basal avian clade Palaeognathae revealed the presence of both DAA1 and DAA2 ( Appendix 3 ). The DAA1 and DAA2 clusters were monophyletic with moderately high to very high posterior probability ( Appendix 2 , 3 ). Assignment of particular sequences to each cluster revealed the presence of both DAA1 and DAA2 in Tinamiformes, whereas the other orders showed the presence of either DAA1 (Casuariiformes) or DAA2 (Rheiformes and Apterygiformes). Phylogenetic analysis of MHC-IIA exon 3 sequences in non-passerine neognath birds revealed the presence of both DAA1 and DAA2 clusters in 10 out of 17 orders ( Figure 2 ). In five orders we recorded only the DAA2 gene cluster and these included both early-branching (Galliformes and Anseriformes) and more derived non-passerine avian clades (Gruiformes, Sphenisciformes, and Accipitriformes) ( Figure 2 ). Whereas the absence of DAA1 sequences in Galliformes and Anseriformes suggests that the loss of DAA1 occurred prior to the split of Galloanserae into these two major lineages, the other clades were scattered across the avian phylogeny, suggesting multiple independent DAA1 losses ( Figure 2 ). At the same time, only the DAA1 cluster was detected in the two most derived non-passerine clades (Falconiformes and Psittaciformes) ( Figure 2 ), indicating a single major transition from DAA1/DAA2 to DAA1 in the evolution of avian MHC-IIA. Bayesian posterior probabilities for monophyletic DAA1 or DAA2 clusters were in most cases high or very high (>0.90; Appendix 3 ).

Whereas the upstream region of MHC class IIA exon 3 provided a good phylogenetic resolution of DAA1 and DAA2 genes across diverse avian orders, this signal was not retained across the remaining part of exon 3 (central and downstream region). Phylogenetic analysis of these sequences in three selected bird orders (Charadriiformes, Procellariiformes, and Coraciiformes) revealed mixed clustering of DAA1 and DAA2 and relatively high consensus support for many mixed DAA1/DAA2 clusters, allowing no reliable resolution of orthology ( Appendix 4 ), and suggesting concerted evolution across the central and downstream regions of exon 3.

Phylogenetic analysis of MHC-IIB exon 3 sequences revealed identical patterns of duplication as found for MHC-IIA ( Appendix 3 ), suggesting coevolution of these genes throughout the evolution of birds ( Figure 2 ). Phylogenetic analysis of MHC-IIB also generally agreed with the results previously presented by Goebel et al. (26), although our analyses lent support for the presence of both DAB1 and DAB2 in some clades, which were previously reported to lack either of gene lineages (Pterocliformes, Otidiformes, Coraciiformes). It was noteworthy that positions of codons most divergent between the two gene lineages (showing amino acid differences between consensus sequences of orthologous gene lineages inferred across all avian orders) clearly overlapped with positions of sites responsible for interactions between the membrane-proximal domains of MHC-IIα and β subunits ( Figure 3 ). Six out of eight divergent codons of the α2 domain overlapped or were adjacent to sites putatively involved in interactions with β2 domain, and five of them showed signal of positive selection along the branch separating DAA1 and DAA2, indicating that the early post-duplication divergence of the α2 domain was likely driven by positive selection. Consistent with this pattern, seven out of eight divergent codons at β2 domain overlapped with sites putatively involved in interactions with α2 domain ( Figure 4 ). There was a minimal overlap of the divergent codons with sites responsible for other types of interactions (with the peptide-binding domain or CD4 coreceptor) ( Figures 3 , 4 ), indicating that functional divergence of MHC-IIA and IIB gene lineages is primarily manifested in the interactions between α2 and β2 domains.

Pervasive negative (purifying) selection was identified as the key evolutionary force acting at α2 and β2 domains encoded by MHC-II exon 3. Between 37% and 44% sites of the membrane-proximal domains showed signature of negative selection ( Table 1 ; Figures 3 , 4 ). In contrast, no more than three sites were found to be under pervasive positive (diversifying) selection, while an additional 2-5 sites were under episodic positive selection, meaning that on average only 5.7% sites at these domains showed any signature of positive selection ( Table 1 ; Figures 3 , 4 ). The signature of selection was unevenly distributed along exon 3 sequence, as dN/dS ratio was significantly lower in the upstream region used for phylogenetic analysis (stronger negative selection) compared with the remaining central/downstream regions (F_1,364 = 4.17, P = 0.042). Stronger negative selection in the upstream region was apparent across both pairs of MHC class IIA and IIB gene lineages, as indicated by non-significant region-gene (F_1,361 = 0.79, P = 0.37) and region-lineage (F_1,361 = 0.01, P = 0.91) interactions. We found no significant differences in the mean strength of selection (dN/dS ratio) at the membrane-proximal domains between MHC class IIA and IIB genes (F_1,364 = 0.74, P = 0.39) or between gene lineages (F_1,364 = 0.04, P = 0.83; gene-lineage interaction: F_1,361 = 0.07, P = 0.78) ( Figure 5 ). Also, the proportion of positively and negatively selected sites did not differ between gene lineages (DAA1 vs. DAA2: G = 0.08, P = 0.78; DAB1 vs. DAB2: G = 0.02, P = 0.89) and there was no evidence for a relaxation or intensification of selection at DAA1 vs. DAA2 (LR = 0.00, P = 1.00) and DAB1 vs. DAB2 (LR = 0.88, P = 0.35). Despite similar overall strength of selection, we identified codon-specific differences in dN/dS ratio at the membrane-proximal domains of both gene lineages. These differences were more apparent at the DAA genes, where we identified 12 codons showing evidence of different selection signature between DAA1 and DAA2. In contrast, we identified only four codons showing different selection signature between DAB1 and DAB2 membrane proximal domains. Most of these codons (87.5%) showed higher dN/dS ratio at DAA1/DAB1 than DAA2/DAB2.

As expected, the peptide-binding domains (encoded by exon 2) showed significantly stronger signature of positive selection than membrane-proximal domains, both at MHC-IIA (F_1,352 = 8.01, P = 0.005) and MHC-IIB (F_1,362 = 36.79, P < 0.001) genes ( Figure 5 ). Nucleotide diversity was on average 45% higher at the peptide-binding than membrane-proximal domains ( Table 1 ). Significantly stronger positive selection (higher dN/dS) was detected at β1 than α1 domain (F_1,347 = 14.83, P < 0.001) ( Figure 5 ). There were only two sites under pervasive positive selection at the α1 domains of DAA1 and DAA2, although additional 8 (DAA1) and 11 (DAA2) codons showed signature of episodic positive selection ( Table 1 ; Figures 6 , 7 ). In contrast, 17 (DAB1) and 23 (DAB2) codons were found to be under either pervasive or episodic selection at β1 domain ( Figures 6 , 7 ). The number of negatively selected sites was similar between α1 and β1 domains (range: 25-29 at each gene lineage) ( Table 1 ; Figures 6 , 7 ). Most of the positively selected sites (PSS) either overlapped with the putative PBS or were located at positions neighbouring to the PBS (a total of 65% and 78% PSS at α1 and β1 domain, respectively) ( Figures 6 , 7 ). There were no differences in the mean strength of selection (dN/dS ratio) at the peptide-binding domains between gene lineages (F_1,347 = 0.10, P = 0.76; gene-lineage interaction: F_1,346 = 0.11, P = 0.74) ( Figure 5 ) and the proportion of positively and negatively selected sites was also similar between lineage (DAA1 vs. DAA2: G = 0.03, P = 0.87; DAB1 vs. DAB2: G = 0.61, P = 0.43). We found no support for either relaxation or intensification of selection at the peptide-binding domains of DAA1 versus DAA2 (LR = 2.74, P = 0.10) and DAB1 vs. DAB2 (LR = 2.41, P = 0.12). However, similarly to the membrane-proximal domains, we found codon-specific differences in dN/dS ratio at the peptide-binding domains of both gene lineages. Specifically, contrasting selection signature was found at 10 codons for DAA1/DAA2 and 12 codons for DAB1/DAB2 gene lineages. The direction of these differences was highly variable between the codons, showing increased dN/dS at either DAA1 (4 sites) or DAA2 (6 sites) and either DAB1 (4 sites) or DAB2 (8 sites). Interestingly, 77.3% of codons with contrasting selection signature were located at positions directly overlapping or neighbouring with the putative PBS.