Allelic diversity and patterns of selection at the major histocompatibility complex class I and II loci in a threatened shorebird, the Snowy Plover (Charadrius nivosus)

Background Understanding the structure and variability of adaptive loci such as the major histocompatibility complex (MHC) genes is a primary research goal for evolutionary and conservation genetics. Typically, classical MHC genes show high polymorphism and are under strong balancing selection, as their products trigger the adaptive immune response in vertebrates. Here, we assess the allelic diversity and patterns of selection for MHC class I and class II loci in a threatened shorebird with highly flexible mating and parental care behaviour, the Snowy Plover (Charadrius nivosus) across its broad geographic range. Results We determined the allelic and nucleotide diversity for MHC class I and class II genes using samples of 250 individuals from eight breeding population of Snowy Plovers. We found 40 alleles at MHC class I and six alleles at MHC class II, with individuals carrying two to seven different alleles (mean 3.70) at MHC class I and up to two alleles (mean 1.45) at MHC class II. Diversity was higher in the peptide-binding region, which suggests balancing selection. The MHC class I locus showed stronger signatures of both positive and negative selection than the MHC class II locus. Most alleles were present in more than one population. If present, private alleles generally occurred at very low frequencies in each population, except for the private alleles of MHC class I in one island population (Puerto Rico, lineage tenuirostris). Conclusion Snowy Plovers exhibited an intermediate level of diversity at the MHC, similar to that reported in other Charadriiformes. The differences found in the patterns of selection between the class I and II loci are consistent with the hypothesis that different mechanisms shape the sequence evolution of MHC class I and class II genes. The rarity of private alleles across populations is consistent with high natal and breeding dispersal and the low genetic structure previously observed at neutral genetic markers in this species.


Background
The genes of the major histocompatibility complex (MHC) are crucial for the immune response in vertebrates [1,2]. Their encoded proteins are involved in presenting antigen derived from pathogens to immune cells, which then trigger a cascade of immune responses [3,4]. Because of their functional importance and the direct link between MHC diversity, fitness and individual behaviour [3,5], the MHC has been the subject of ecological and evolutionary studies ranging from assessing individual survival and mate choice to the processes of speciation and practical conservation management [6][7][8][9][10]. Adaptive genes, which are directly associated with individual fitness, are important for population viability and hence conservation [2,11]. The loss of adaptive genetic diversity has been associated with an increase in the risk of extinction, especially in species with low population sizes [11]. The maintenance of MHC diversity is crucial for pathogen resistance, which represents one of the principal selective forces impacting fitness and long-term survival of endangered species [2,12].
MHC genes display the highest degree of polymorphism within vertebrate genomes [5,13]. Pathogenmediated selection results in positive selection and the substitution of amino acids in the codons of the peptide-binding region (PBR), as well as balancing selection including heterozygote advantage, frequencydependent selection and fluctuating selection [2,4,14]. Recently, a number of studies showed that the evolutionary dynamics of the MHC genes is driven by high rates of recombination, duplication and conversion [15][16][17]. Through these processes populations can respond to a great number of antigens [10,17]. The MHC genes are divided into two principal classes: class I, which is responsible for immune defence against intracellular pathogens such as viruses, and class II, which is responsible for dealing with extracellular pathogens such as bacteria and nematodes [2].
The number of MHC genes varies between and within species [10]. In mammals, MHC genes are organized into a dense genomic region and are characterized by complex organization and many pseudogenes, leading to extraordinary genetic diversity. For example, in humans approximately 9000 class I alleles and 3000 class II alleles have been described [18]. In birds, the structure and organization of the MHC region varies not only between, but also within the same family [9]. Some groups, such as chicken Gallus gallus and some birds of prey, have an extraordinarily compact MHC region (coined as the "minimal essential" MHC, [19,20]). However, other galliform species have duplications, leading to many MHC alleles [21,22]. In contrast, in Passeriformes, the MHC shows a complex architecture, and is frequently composed of multiple expressed loci and pseudogenes [1,23,24]. Other groups of birds, such as the Charadriiformes, appear to have a diversity and complexity intermediate between chicken and Passeriformes [25,26]. The differences in the number and organization of the MHC genes in vertebrates might be best explained by different evolutionary dynamics in the birth and death of genes [27]. Here, new genes are generated by duplication, with some daughter copies conserving their function while others are inactivated or eliminated from the population [10,27].
Within the order Charadriiformes, study of the MHC have so far been restricted to three families: Alcidae, Laridae and Scolopacidae. These studies revealed considerable differences in the diversity and organization of the Charadriiform MHC [16,25,26,28,29]. To investigate which mechanisms generate and promote MHC evolution and diversity [9], studies in more phylogenetically distinct families are required.
We investigated the diversity and organization of the MHC in the Snowy Plover Charadrius nivosus, a member of the Charadriidae [30]. Until recently, Snowy Plovers were lumped with Kentish Plovers Charadrius alexandrinus and considered to be a cosmopolitan species, but the analysis of genetic and phenotypic traits has shown that the two species are separate and that Snowy Plovers are characterized by low genetic diversity at neutral genetic markers [31][32][33]. Three Snowy Plover lineages are commonly recognized and distinguishable with genetic markers: Western Snowy Plovers (C. n. nivosus) in North America, Cuban Snowy Plovers (C. n. tenuirostris) in the Caribbean, and the Peruvian Snowy Plover (C. n. occidentalis) in South America [31]. Coastal populations in particular have declined throughout the range and now require active conservation management [34][35][36][37]. Here, we develop markers for MHC class I and class II genes to examine adaptive genetic diversity in Snowy Plovers. Our study has four aims: First, we characterize functionally essential regions of the MHC class I and class II loci to provide novel genetic markers for studying adaptive diversity in this species. Second, we investigate the genetic diversity (number of segregating sites and nucleotide diversity), evolutionary distance and type of selection acting on MHC class I and class II alleles. Third, we compare the diversity of MHC classes I and II with the respective diversity observed in other Charadriiformes. Finally, we describe the pattern of diversity across the geographic range and test for the presence of private alleles within Snowy Plover populations.

MHC class I: exon 3
We discarded 32 samples because they did not pass the quality filters or because they did not have the minimum number of reads per amplicon. Among the 218 remaining samples from eight populations, we found a total of 40 alleles. Five alleles ('long alleles'; Chni-UA*12 to 16) showed a 3-bp insertion at position 178-180 (amino acid position 60 in the alignment) in the α2 domain, retaining the correct reading frame. Meanwhile, the remaining 35 alleles ('short alleles'; Chni-UA*01 to UA*11 and Chni-UA*17 to UA*40) did not show this insertion. The alignment of exon 3 displayed five of the eight highly conserved amino acid (peptide main-chain) sites in birds (amino acids: TKWYY, Fig. 1). Alleles Chni-UA*01 to UA*10 had an amino acid substitution at the second of these conserved sites (N for K), whereas in alleles Chni-UA*17 to UA*21 the fourth site was substituted (C for Y). Four of the 18 highly conserved intra-and interdomain contacts described in vertebrates [38] were present, and none of these showed a polymorphism.

Allelic diversity
We found two to seven alleles per individual (x ± SD: 3.70 ± 0.92), which suggested that we obtained alleles from up to four loci, assuming heterozygosity. We detected up to three non-classical alleles in 216 of 218 samples. The number of alleles did not differ between individuals across populations (Table S1). However, we found differences in the number of alleles per population, with populations Nayarit and Puerto Rico showing fewer alleles than the other populations (Table S2). Long alleles were less common among the 218 individuals genotyped, with only 23 individuals displaying one or two long alleles, whereas the short alleles were present in all individuals. Chni-UA*21 was the most common allele and detected in 173 individuals (83.5%), followed by Chni-UA*20 in 102 individuals (49.2%) and Chni-UA*30 in 69 individuals (33.3%). All individuals genotyped had at least two alleles. Most individuals (48.3%) had four alleles, 24.6% had three alleles, 13% had five alleles, 11.6% had two alleles, 2.4% had six alleles and 0.9% had seven alleles.

Diversity and inference of selection
The average nucleotide diversity (π) for the complete sequence was similar among the three lineages, ranging from 0.05 to 0.07. Populations at Nayarit and Puerto Rico showed the highest levels of nucleotide diversity, the nucleotides distance and average nucleotide diversity at PBR sites, whereas populations at Utah, San Quintin and Ceuta had the lowest levels of nucleotide diversity (both complete sequence and PBR sites) and nucleotide distance ( Table 1). Within all populations, PBR sites showed higher diversity than non-PBR sites, suggesting balancing selection at these sites (Table 1). Bayesian analysis of selection (FUBAR: Fast Unconstrained Bayesian AppRoximation) identified six sites (sites 7, 23, 25, 58, 66 and 69) under positive selection and 11 sites (sites 1, 6,16,27,28,33,41,46,64,80 and 85) that displayed diversity patterns consistent with negative, purifying selection (Fig.  1a). Sites with purifying selection were exclusively located in the non-PBR region. Also, differences in nonsynonymous substitution rate and synonymous substitution rate (dN/dS) suggested stronger positive selection at PBR sites in comparison with non-PBR sites ( Table 1). Among genetic lineages, C. n. tenuirostris (Puerto Rican population) showed the lowest level of positive selection. Using GARD (Genetic Algorithm for Recombination Detection) we found no evidence for recombination among the 40 alleles. Table 1 Diversity at MHC class I exon 3 in the Snowy Plover (Charadrius nivosus). Segregating sites of amino acids (S aa ), average nucleotide diversity (π), evolutionary distance of nucleotide sequences (d nt ) and amino acid sequences (d aa ). The measures of diversity, and the synonymous and non-synonymous substitution rates, were calculated for the complete sequence (All), and separately for the PBR and non-PBR sites within each population. The number of genotyped individuals per population (N) and the number of samples that passed the quality filters (n) are shown. To have comparable measures of diversity we randomly selected 15 allele sequences for each population, with the exception of Nayarit and Puerto Rico, where we found only 13 and 9 alleles, respectively

Comparison and phylogenetic relationships with other Charadriiformes
Average nucleotide diversity and number of segregation sites were lower than those reported from Red Knot, Icelandic Black-tailed Godwit and Red-billed Gull ( Table 2; Figure S1). Nevertheless, the dN/dS proportion at PBR sites was higher for the Snowy Plover, indicating stronger positive selection at these sites in comparison with the other three charadriiform species. As in Red Knots and Red-billed Gulls, we found alleles with putatively non-classical functions among the Snowy Plover alleles. These non-classical alleles (Chni-UA*17 to UA*21) formed a well-defined cluster in the phylogenetic network (Fig. 2a) and showed a low level of polymorphism (Figs. 1a and 2a). Non-classical alleles were present in 216 of 218 samples. These alleles clustered together with the Red-billed Gull non-classical allele Lasc-UDA*03 in the neighbour-joining tree (Fig. 3a).

MHC class II: exon 2
We discarded 36 samples, as they did not pass the quality filters. In total, we found six alleles across 214 individuals of eight populations. We found no more than two alleles per individual ( x ± SD: 1.45 ± 0.50), which suggests that we genotyped only one locus. The most common allele was Chni-DAB*01, which we detected in 171 individuals (79.9%), whereas the least common allele was Chni-DAB*03, present in four individuals (1.9%). Most individuals (55.6%) were homozygous.

Diversity and inference of selection
The nucleotide diversity (π) for subspecies nivosus and occidentalis were similar, with the PBR sites showing higher diversity in comparison to the non-PBR sites, suggesting that balancing selection could be acting at PBR sites ( Table 3). The Bayesian site-by-site test in FUBAR suggested positive selection for one site (site 82), and purifying selection for two sites (sites 17 and 73, Fig. 1b). Sites with purifying selection were exclusively located in the non-PBR region. The analysis of rates of changes dN/dS indicated positive selection for the PBR sites in comparison to non-PBR sites for the nivosus and occidentalis subspecies ( Table 3). All tenuirostris samples were homozygous. GARD identified one recombination breakpoint, which was located at position 222. When we re-ran a codon-specific model for the non-recombinant fragment of the sequences, our results remained unchanged.

Diversity and phylogenetic relationships within the Charadriiformes
The average nucleotide diversity of Snowy Plovers (π = 0.06) was well within the range of diversity observed in other Charadriiformes (range 0.02 to 0.15, Table 4; Figure S1). This was also true for the nucleotide diversity at PBR sites ( Table 4). The phylogenetic network and the neighbor joining tree showed two distant allele groups for MHC Class II (Figs. 2b and 3b). The neighbour-joining tree showed that MHC class II alleles of Snowy Plovers (Chni-DAB) were located on a different branch than most other charadriiform MHC II alleles identified so far (Fig. 3b).

Geographic pattern of MHC diversity
We found that 25% of alleles were private alleles for MHC class I and 16.7% private alleles for MHC class II. For MHC class I, Snowy Plover populations from Utah, Ceuta and Puerto Rico showed private alleles (Fig. 4a). For Utah, four (Chni-UA*03, 04, 31 and 37) of 30 alleles were private, for Ceuta, four (Chni-UA*25, 28, 32 and 33) of 34 alleles were private and for Puerto Rico, two (Chni-UA*05 and 17) of nine alleles were private. Other populations lacked private alleles (notably the Peruvian population of the subspecies occidentalis) but showed similar numbers of alleles (San Quintin: 23 alleles, Nayarit: 13 alleles, Texcoco: 19 alleles, Florida: 17 alleles and Table 2 Comparison of genetic diversity at MHC class I exon 3 in four different charadriiform species. Segregating sites of amino acids (S aa ), nucleotides (S nt ), and average nucleotide diversity (π), from 21 sequences randomly chosen by each of these species. Diversity indices and the synonymous (dS) and non-synonymous (dN) substitution rates were calculated for the complete sequences (All), the PBR and the non-PBR sites  Peru: 18 alleles). With the exception of the island population of Puerto Rico (lineage tenuirostris), private alleles generally occurred in low frequencies (Fig. 4). At MHC class II, there were no private alleles within lineages tenuirostris or occidentalis. Only the Ceuta population had a private allele (Chni-DAB*03), and we found all six alleles in this population. The other populations had three to five alleles, except for Puerto Rico, where we detected only one allele (Fig. 4b).

Discussion
In this study, we developed new MHC markers that amplified with high success exon 3 of MHC class I and exon 2 of MHC class II in the Snowy Plover, a shorebird species of high conservation concern. This allowed us to examine adaptive genetic diversity at this important locus in the family Charadriidae and provides novel and useful markers for future studies in other Charadrius species. Surveying both MHC classes in more than 200 individuals from eight populations, we report differences in allelic diversity across both MHC classes, with nearly seven times as many alleles at MHC class I than at MHC class II. The genotypic variation in individuals suggests that the markers amplified four highly similar loci for MHC class I, as we registered up to seven alleles at MHC class I, whereas a maximum of two alleles per individual is consistent with amplification of a single locus at MHC class II. Our network analysis suggests the presence of classical and non-classical genes among the amplified MHC class I loci. Non-classical loci have evolved from classic MHC genes to perform specific tasks in the immune recognition [39]. Five alleles (Chni-UA*17 to 21), formed a welldefined cluster in the phylogenetic network, which is characteristic for non-classical alleles. Non-classical alleles have also previously been described in other shorebirds (Red-billed Gull: [29]; Red Knot: [25]). Consistent with expectations derived from other nonclassical alleles, these five alleles showed a low level of polymorphism [40,41]. Four alleles were only differentiated by synonymous substitutions, whereas the fifth allele, found only in the Puerto Rico population, deviated in a single amino acid substitution. We observed that two of three non-classical alleles Lasc-UCA*01 and Caca-UA*35 clustered with the classical alleles of their respective species Red-billed Gull and Red Knot. This suggests that they may have evolved through recent duplications that have occurred after the divergence of the species in this tree. By contrast, Lasc-UDA*03 and the Snowy Plover non-classical alleles may have been created through a more ancient duplication event.
We observed a second well-defined cluster in our phylogenetic network at MHC class I; this cluster corresponds to five alleles (Chni-UA *12 to 16) that shared a 3-bp insertion. The observed insertion kept the reading frame intact, presumably preserving their function. Contrary to the non-classical alleles, these five alleles with a 3-bp insertion show a higher nucleotide diversity, including several non-synonymous changes between these five alleles, suggesting these alleles have classical functions. In birds, insertions are less frequent than deletions in MHC genes, and it has been suggested that insertions have a reduced adaptive advantage [42]. However, some insertions may contribute to adaptive MHC variation due to possible changes in the PBR [42,43]. Unlike the MHC class I, the six alleles at MHC class II formed an undifferentiated cluster, with both synonymous and non-synonymous substitutions present. Together, with the observation of no more than two alleles per individual, we concluded that our markers amplified a single classical MHC II locus. When we compared the nucleotide diversity at the PBR sites, we found that the MHC class I showed moderately higher values of diversity (0.16 ± 0.03) than at MHC class II (0.12 ± 0.03). This result may be a consequence of having amplified four loci at MHC I and only one locus at MHC II.
For both MHC classes, we observed that the number of sites under negative (purifying) selection was higher than sites under positive (diversifying) selection when we considered the full sequence. As PBR sites are interacting directly with the antigens derived from pathogens, it is expected that these sites are subject to stronger positive diversifying selection than non-PBR sites. Consistent with this, we found that sites with signatures of purifying selection were predominantly non-PBR sites, whereas five of the six sites under diversifying selection in the MHC class I and in one site at MHC class II corresponded to PBR sites. Similar patterns were observed in another shorebird, the Icelandic Black-tailed Godwit [26], and shown in a recent comparative study of selection at the avian MHC [42].
Our finding that more sites showed signatures of positive selection at MHC class I than at MHC class II (See figure on previous page.) Fig. 3 Neighbour-joining tree of amino acid sequences for MHC class I and II from different Charadriiformes. For MHC class I, four species are represented (Snowy Plover (Chni), Black-tailed Godwit (Lili), Red Knot (Caca) and Red-billed Gull (Lasc)) and (b) MHC class II represented by eight species (Snowy Plover (Chni), Black-tailed Godwit (Lili), Ruff (Phpu), Great Snipe (Game), Black-legged Kittiwake (Ritr), Razorbill (Alto), Common Murre (Uraa) and Atlantic Puffin (Frar)). Long alleles are indicated by a white square, non-classical alleles with a black square. The homologous MHC sequences of the chicken Gallus gallus (Gaga-HQ141386 and Gaga-BLB2) served as an outgroup differs from a recent global analysis of selection at the avian MHC [42]. Here, the authors suggested that the pressure of extracellular pathogens is higher in nonpasserines, resulting in a stronger signature of selection for the MHC class II in non-passerines than in passerines. The order Charadriiformes does not seem to fit to this proposed rule. Consistent with Minias et al. [42], gulls and most of the alcids show a strong signature of selection at MHC class II and weaker selection at MHC class I, but the paraphyletic group of shorebirds (plovers and sandpipers) instead shows a pattern more similar to passerines [16,25,26,28,29,44]. Several morphological or ecological variables may explain this discrepancy. First, body size may be related to parasite abundance, as larger hosts may provide a greater variety of niches and, in turn, support a higher number of parasites than smaller birds [45]. Snowy Plovers are relatively small birds with a mean body mass of 38-50 g [32,46]. Second, the selective pressures imposed by parasites may be habitat dependent. Although aquatic birds (mainly non-passerines) show a more diverse parasite community than their terrestrial counterparts [45], there is a difference between freshwater and saltwater habitats. Snowy Plovers inhabit the shores of alkaline water bodies, such as salt lakes, salt evaporation ponds and sandy beaches [46]. These saline habitats are typically considered to have a lower abundance of extracellular parasites [47][48][49], which would be consistent with the observed low diversity at MHC class II. In general, shorebirds show a low prevalence in intracellular pathogens, although viral infections (West Nile virus; [50], Newcastle disease virus [50], avian influenza; [51]), avian haemosporidians (Plasmodium and Haemoproteus spp.; [52]) and bacterial infections (Mycobacterium; [53]) have all been reported. Furthermore, Snowy Plovers inhabit low latitudes, where the diversity of intracellular Table 3 Diversity at the MHC class II (segregating sites of amino acids (S aa ) and average diversity of nucleotides (π)) and evolutionary distance (nucleotide sequences (d nt ) and amino acids (d aa )) of the alleles of the exon 2 in the Snowy Plover (Charadrius nivosus). The measures of diversity and the synonymous and non-synonymous substitution rates were calculated for the complete sequences (All), the PBR and the non-PBR sites for each population. The number of genotyped individuals per population (N) and number of samples that passed the quality filter (n) are shown  (17) All pathogens, as well as of their vectors, is expected to be high [48,52]. In addition to the copy number variation for the fragment amplified at MHC I, the high abundance of intracellular pathogens in the tropics may contribute to the high allelic diversity at MHC class I. Further research is needed to determine whether Charadriiformes themselves show unusual group-specific variation in intra-and extracellular pathogens, or whether other life-history, ecological and evolutionary differences explain the observed differences in signatures of selection. Across continental populations (subspecies nivosus and occidentalis), we found a similar number of alleles and nucleotide diversity. This result is consistent with the high gene flow observed among Snowy Plover populations [31,33,54]. At MHC I only, we found that the Nayarit population had a lower number of alleles (13 alleles), although this may reflect the sample size, as we only genotyped ten birds from this population. Among continental populations, private alleles were generally rare, and when present, occurred at low frequencies, suggesting that these represent rare alleles rather than geographically distinct alleles. By contrast, the island population of Puerto Rico (lineage tenuirostris) showed two privates alleles at MHC class I at moderate frequencies (Fig. 4): these alleles may represent true geographic variants. Interestingly, the tenuirostris population was monomorphic at MHC class II. Together with the lower positive selection at PBR sites at MHC class I, this might suggest that pathogen pressure is weaker in this population. The low genetic diversity may be a common feature of island populations that are thought to be exposed to fewer pathogens than continental pathogens [55]. For example, among shorebirds, Icelandic Black-tailed Godwits did not show positive selection at MHC class I PBR sites [26]. Other biogeographic features could contribute to the observed differences in genetic diversity. In contrast to most other populations, tenuirostris inhabits the Atlantic Ocean and the diversity of pathogens is thought to be lower in the Atlantic Ocean than the Pacific Ocean [56]. However, the other Atlantic population located in Florida did not show lower diversity than Pacific populations. Demography may also have played a role in shaping MHC diversity. A recent study showed that all Snowy Plover lineages went through a bottleneck within the last 1000 years, with particularly strong effects observed in C. n. tenuirostris [33]. A similar pattern of loss of adaptive diversity has been observed in other bird populations subject to recent bottlenecks (see [57][58][59]).
Comparing allelic diversity across species, we found that diversity at MHC class I in Snowy Plovers (40 alleles across four loci amplified) was similar to other Charadriiformes species (Red-billed Gull [29]; (38 alleles), Red Knot [25]; (36 alleles) and Black-tailed Godwit [26]; (47 Table 4 Comparison of genetic diversity at MHC class II (segregating sites of amino acids (S aa ) and nucleotides (S nt ) and average nucleotide diversity (π)) in eight species of Charadriiformes, from six alleles randomly selected for each species (except for Black-tailed Godwit and Ruff where only four alleles were available). The measures of diversity and the synonymous and non-synonymous substitution rates were calculated for the complete sequences (All), the PBR and the non-PBR sites alleles)). By contrast, MHC class II (six alleles) showed a lower allelic diversity than other Charadriiformes (Great Snipe -50 alleles, [28]; Marbled Murrelet Brachyramphus marmoratus -27 alleles, [60]). Despite the copy number variation, the nucleotide diversity at MHC class I in Snowy Plovers was lower than in other Charadriiformes (Table 2). However, the observed values were within the range of other non-passerines (such as the Humboldt Penguin Spheniscus humboldti -0.03 to 0.04 and Magellanic Penguin Spheniscus magellanicus -0.05 to 0.06, [61]; or grouse specieson average 0.05, [58]). For MHC class II, we found intermediate nucleotide diversity in comparison to other Charadriiformes (Marbled Murrelet -0.08; [60]) and other nonpasserines (Magellanic Penguin -0.02, [61]; Eurasian Coot Fulica atra -0.11, [62]; and Black Grouse Tetrao tetrix -0.11, [63]).

Conclusions
We developed novel MHC markers to amplify the PBR exon 3 of MHC class I and PBR exon two of MHC class II for the threatened Snowy Plover. These are the first markers for MHC in the family Charadriidae and we anticipate that they will be of high utility for studying MHC in other plover species. Overall, genetic diversity at MHC in Snowy Plovers was low to moderate and likely to be shaped by past demographic processes such as bottlenecks and island colonization. In line with population genetic studies, we find that there is limited genetic differentiation attributable to geographic variation, consistent with the high gene flow observed in this species. Contrasting differences in the allelic diversity between MHC class I and class II indicate stronger positive selection at MHC class I than at MHC class II. These differences may reflect variation in exposure to intracellular and extracellular pathogens [42], but further studies are needed to confirm this.

Population sampling
We collected blood samples of 250 adult and juvenile Snowy Plovers from eight populations across North, Central and South America between 2006 and 2016 (Table 1). We captured adults in funnel traps placed on nests during incubation or used mist nets at other times using the methods described in Székely et al. [64]. After capturing and banding the individuals, we took 20 to 75 μl blood with heparinized capillaries from the brachial vein. We captured chicks at or near to the nest a few hours after hatching and took~20 to 50 μl blood from the tarsal vein. We stored the blood in 1 ml of Queen's lysis buffer [65] at 4°C or pure ethanol at room temperature until further processing.

DNA extraction
We isolated genomic DNA using the ammonium acetate precipitation method [66]. We checked the integrity of the DNA using a 0.8% agarose gel stained with SYBRsafe (Invitrogen). We measured DNA concentration using either a fluorometer (FLUOstar OPTIMA) or Nanodrop ND800 (Thermo Fisher Scientific).

Primer design for the MHC loci
We designed primers to capture the most polymorphic PBR sites in exonic regions for both MHC class I and MHC class II genes in Snowy Plovers. For MHC class I, we initially used the primers MHCI-int2F [20] and MHCI-ex3R [67] to isolate exon 3 of non-passerine birds. We undertook polymerase chain reactions (PCRs) in a total volume of 20 μl containing 12 μl Multiplex PCR Master Mix (MM, Qiagen), 4 μl Q-Solution (Qiagen), and 1 μl of each primer (1 μM) and 2 μl DNA (~15 ng/μl). The PCR program started with an initial denaturation step at 95°C for 15 min, followed by 30 cycles at 94°C for 30 s, 56°C for 90 s and 72°C for 90 s, and a final elongation step at 72°C for 10 min. For MHC class II, we first used primers MHC05 [68] and 305 [69], and the primers 306 [69] and RapEx3CR [70] to capture introns 1 and 2, and parts of exons 1 and 3, respectively. We ran PCRs in a total volume of 20 μl, containing 16 μl MM, 1 μl of each primer (1 μM) and 2 μl of DNA (~15 ng/μl). The PCR program consisted of one cycle at 95°C for 3 min, followed by 30 cycles at 94°C for 30 s, 60°C for 90 s and 72°C for 90 s, followed by a final step at 72°C for 10 min. All PCRs were run on a thermocycler PTC-225 DNA Tetrad Engine. We visualized PCR products using an agarose gel at 1.5% stained with ethidium bromide. For MHC class II we obtained multiple bands and subsequently cut out the visible bands of the expected size and extracted the amplified fragment using the QIAquick Gel Extraction Kit (QIAGEN). For MHC class I we only observed a single band, and the product did not require gel excision. We cleaned up MHC class I and II amplicons with ExoSap and sequenced the products using the BigDye terminator v.3.1 chemistry (Applied Biosystems) on an ABI 3730 automated sequencer (Applied Biosystems). For each MHC locus we aligned the sequences of six individuals using CodonCode Aligner 5.0.2 (CodonCode Corporation). We confirmed the identity of the sequences through blast hits in Gen-Bank (NCBI). We then designed new primers Chni-Ex2F 5′-GAACTGCCTCCCTGCACAAA-3′ and ChCR-Ex2R 5′-TTCCCCGGGGAAATGTTCT-3′ to amplify the complete exon 2 for the MHC class II; and ChCR_MHC I_Ex2aF 5′-GGGTCTGTGCCCCACT-3′ for use with primer MHCI-ex3R 5′-CTCACCTTTCCTCTCCAG-3′ [41] to amplify the complete exon 3 of MHC class I.

Amplification and sequencing
We adopted a two-step PCR protocol to amplify the PBRs of both MHC classes and enable multiplexing. For PCR1 we created new oligonucleotides adding the Illumina overhang sequencing adapters F 5′-TCTACACG TTCAGAGTTCTACAGTCCGACGATC-3′ and R 5′-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT -3′ to the MHC primer sequences (following Campbell et al. [71]). We performed PCRs in a total volume of 10 μl that contained; for MHC class I, 3.5 μl MM, 1.25 μl Q-solution, 1 μl of each primer (10 μM), 1 μl DNA (25 ng/μl) and 2.25 μl water, and for MHC class II, 4 μl MM, 1 μl of each primer (10 μM), 1 μl DNA (25 ng/μl) and 3 μl water. We used the same PCR programs as before but reduced the number of cycles to 28 to minimize the impact of chimera formation [72]. We then checked 1 μl of the products on an agarose gel and cleaned up the remainder of the solution using 8 μl (concentration of 0 We determined the concentration of the PCR2 product using a fluorometer (FLUOstar OPTIMA) and 2 μl of the product. We pooled samples from eight individuals by taking 20 ng per sample and cleaning up the multiplexed PCR products with AMPure XP beads, as described above, with volumes adjusted to a 50 μl solution (concentration of 0.5 X). We used a TapeStation 4200 (Agilent Genomics) to confirm that there were no primer dimers present in the purified samples. We then quantified the PCR products using a qPCR with the KAPA library quantification kit (KAPA Biosystem) using 10 μl reaction volume (8 μl of SYBR Master Mix and 2 μl of template/standard or control), with the program: 95°C for 5 min, followed by 35 cycles of 94°C for 30 s and 60°C for 45 s. We then pooled equimolar amounts per library, preparing six libraries in total, quantified the concentration of the pool with a Qubit (ThermoFisher Scientific, Waltham, USA) and submitted 4 nM per library for sequencing using 250 bp paired-end (500 cycles) Illumina sequencing on the MiSeq (Illumina Inc., San Diego, CA, USA) in six separate runs at the Sheffield Diagnostic Genetics Service.

Processing of data and MHC alleles validation
For the raw MiSeq data processing, we used the Amplicon Sequencing Analysis Tool (AmpliSAT) web server [73]. This tool is divided into different modules that allow the merging, cleaning and assignment of genotypes. First, we used AmpliMERGE with the FLASH algorithm to merge the pair-end reads. Then, we used AmpliCLEAN to filter out low-quality reads (<Q30 score and < 270 bp). After running AmpliCHECK with the default parameters for Illumina sequences, we retained all the remaining reads with lengths of 350 (325) ± 5 bp for MHC class I (II). Finally, we used AmpliSAS to demultiplex, cluster and filter the retained reads using default parameters for Illumina data for clustering, a minimum read depth per amplicon of 2000 and merging minimally different sequences to the dominant sequences when they differed by less than 3 bp and had ≤25% of the read depth in comparison with the dominant sequences [73]. Sequences that differed by 1-3 bp from the dominant sequences, but had more than 25% of the read depth, were classified as 'subdominants' and formed a new cluster representing a putative allele [26,73]. We discarded all sequences that had a frequency of less than 5%, and those identified as chimera sequences. The minimum amplicon depth was set to 150 reads and the maximum amplicon depth set to 5000 reads due to computational limitations [73].

Allele validation
We blasted all putative alleles from Illumina sequencing to the GenBank (NCBI) nucleotide database to examine their similarity to known MHC alleles from other species. The alleles were named Chni-UA*01 to UA*40 (MHC class I) and Chni-DAB*01 to DAB*06 (MHC class II), following the nomenclature suggested by Klein et al. [74].

Diversity analysis and tests for selection
We used MEGA 7.0.18 [75] for initial diversity and selection analysis. First, we aligned the putative alleles using the MUSCLE algorithm [76] implemented in MEGA. We manually checked indel sites and curated alignments in order to preserve triplets within exons. We then estimated the number of segregating amino acid sites (S aa ), nucleotide diversity (π), evolutionary distance for nucleotide sequences (d nt ) and evolutionary distance for amino acid (d aa ), for each of the eight populations included in the study. For the evolutionary distance analyses, we inferred the appropriate substitution models based on the best-fit model (using AIC C ) using JModeltest 2.1.10 [77]. For MHC class I, we employed the Kimura two-parameter model [78] with a gamma distribution, setting the transition rate α to 0.9 for nucleotide sequences, and using the p-distance model with uniform rates to assess amino acid sequences distances. For MHC class II, we implemented the Jukes-Cantor + G model with a gamma distribution, setting the substitution rate α to 0.8 for nucleotide sequences, and we calculated amino acid sequence distances using a p-distance model with uniform rates.
We calculated standard errors of the mean evolutionary distances from 1000 bootstrap replicates. Positive selection, as a response to the selection imposed by pathogens, will lead to an excess of non-synonymous over synonymous substitutions in the PBR (ω = dN/dS > 1). To assess the impact of positive selection on nucleotide diversity, we compared ω at PBR sites and non-PBR sites. We inferred PBR sites based on previously documented transcripts (MHC class I: [79], MHC class II: [80]). We calculated synonymous and non-synonymous substitution rates through the modified Nei-Gojorobi method [81] with Jukes-Cantor correction. Also, we tested site-by-site selection applying Fast Unconstrained Bayesian AppRoximation (FUBAR; http://www.datamonkey.org/fubar, [82]. As recombination is frequent in MHC genes and recombination may lead to overestimation of the number of positively selected sites [62], we tested for evidence of recombination in our MHC alignments using Genetic Algorithm for Recombination Detection (GARD; http:// www.datamonkey.org/gard, [83]).

Phylogenetic diversity and relationships
We compared nucleotide diversity at MHC loci across the Charadriiformes using data available at GenBank. For MHC class I, we obtained sequences from three other charadriiform species: Red Knot Calidris canutus [25], Icelandic Black-tailed Godwit Limosa limosa islandica [26] and the Red-billed Gull Larus scopulinus [29]. As only sequences from 21 alleles were available for Red-billed Gull, we randomly drew sequences of 21 alleles from each species in this comparison to obtain a comparable sample size for the diversity estimate (Table  S3). For MHC class II, we obtained data from seven further species: Common Murre Uria alge [44], Razorbill Alca torda [44], Atlantic Puffin Fratercula artica [44], Black-legged Kittiwake Rissa tridactyla [44], Great Snipe Gallinago media [28], Ruff Philomachus pugnax [16] and Black-tailed Godwit Limosa limosa [16]. As we had only six putative alleles in Snowy Plover, we capped the number of sequences to six per species that we randomly drew for this comparison (Table S3). We evaluated the amino acid (S aa ) and nucleotide (S nt ) segregation sites, nucleotide diversity (π), as well the synonymous (dS) and non-synonymous (dN) substitution rates, using the same parameters described above. We visualized phylogenetic relationships between the MHC class I and class II alleles in Snowy Plovers through the Neighbor-net algorithm implemented in SplitsTree 4.14.6 [84]. This method allows deduction of alternative phylogenetic histories and model incompatibilities in the dataset that may lead to conflicting phylogenetic signals because of duplication, recombination and gene conversion, which are all common in MHC genes [20,25]. Finally, we inferred the phylogenetic relationships among charadriiform MHC alleles using a Neighbour-Joining Tree with Maximum Likelihood implemented in MEGA 7.0.18. We calculated branch support through 1000 bootstrap replications.
Additional file 1: Figure S1. Graphical representation of nucleotide diversity and dN/dS ratios for MHC class I and MHC class II. Silhouettes represent Charadriiform species, i.e. MHC class I (from top to bottom: Red Knot, Black-tailed Godwit, Red-billed Gull and Snowy Plover) and MHC class II (from top to bottom: Black-tailed Godwit, Ruff, Black-legged Kittiwake, Atlantic Puffin, Common Murre, Snowy Plover, Great Snipe and Razorbill). Arrows represent strength of selection in both MHC classes for the Snowy Plover in comparison to the other species. Table S1. Results of generalized linear models testing the differences in the number of alleles per individual between populations for the MHC class I and class II in the Snowy Plover. Table S2. Results of generalized linear models testing the differences in the number of alleles between populations for the MHC class I and class II in the Snowy Plover. Table S3. Sequences ID list randomly drawn for the MHC class I and MHC class II species comparison. Consent for publication Not applicable.