In this study, characterization of great tit Mhc class I exon 3 was adopted and large-scale, high-resolution genotyping was carried out by the use of 454 pyrosequencing method. A total of 862 alleles with varying sizes and high nucleotide diversity were detected in 857 great tits, demonstrating that this species has a polymorphic, multilocus Mhc region; tests of historical selection implied this polymorphism was mainly maintained by balancing selection. Full analysis of sequences through a stepwise variant validation procedure allowed reliable typing of 871 samples and duplicates confirmed the repeatability of the genotyping method. Lastly, the effective Mhc repertoire of each great tit was described by clustering alleles with similar antigen binding affinities into 17 supertypes.
Our characterization effort proved to be effective as it enabled us to define the allelic groups of interest and the sequences to prevent (pseudogenes) during 454 pyrosequencing. High diversity among the potentially functional alleles led to the design of a degenerate primer pair, however attention was paid to minimize amplification of alleles with stop codons. Consequently, functional sequences constituted the major fraction of the reads (87%) that were retained following variant validation procedure, whereas the number of reads belonging to putatively non-functional genes were much lower. This highlights the importance of carrying out a thorough background work prior multilocus typing, not only to attain an initial understanding of gene and allelic diversity, but also to design the ideal primers that target the expressed alleles only. For instance, Zagalska-Neubauer et al.  used 454 pyrosequencing to type the Mhc class II region of the collared flycatcher and the experiment successfully generated a mean of 541 reads per individual. However the coverage was insufficient for reliable genotyping since the primers extensively amplified the putative Mhc pseudogenes, greatly decreasing the coverage of expressed alleles.
We retrieved 76% of the alleles identified during the characterization stage. These were randomly distributed in the phylogeny and represented each allelic group, suggesting that we successfully amplified the allelic groups of interest. Few of the alleles that were not retrieved belonged to Group 1 or the pseudogene group as expected. The rest of the alleles identified at the initial characterization stage were most likely PCR-generated errors that were not differentiated at the characterization stage as they were either 1 basepair different from more common alleles or they could be explained as a chimera of more common alleles. Lenz and Becker  showed that PCR generated artefacts could constitute up to 25 percent of the alleles when no approach is taken to reduce artefact formation. Although we adjusted the PCR conditions during amplification for 454 genotyping, we used the standard conditions during characterization. Hence it would not be surprising to have generated a relatively high frequency of artefacts during the cloning –sequencing processes.
We used the 16 regions of the Pico Titer Plate and pooled 96 great tit amplicons in each region to maximize the number of individuals sequenced in the experiment. However from the total of 1536, only 871 samples (56%) were successfully genotyped following the variant validation procedure. The efficiency of the experiment could have been improved by optimising the sample number for a full plate run . In amplicon sequencing, the depth of coverage is tightly linked to the number of amplicons being pooled. Therefore by pooling a smaller number of amplicons, we could have genotyped more individuals. However, it is rather difficult to determine the optimal number of samples to use in a 454 run, especially in cases where the gene copy number is initially unknown. Alternatively, our method could have been improved by calculating the concentration of each amplicon using a Nanodrop, in order to assure the pooling of equimolar quantities . This way the variation in read number per sample could have been minimized and possibly more individuals would have passed the minimum coverage threshold. Due to the error-prone sequencing technology of NGS, we applied a five-step variant validation procedure to differentiate real alleles from PCR/sequencing artefacts. It can be argued that some rare alleles might have been missed in few individuals (especially in the ones that had read numbers slightly exceeding 200) while these steps were being applied. However we believe this is unlikely to be a very large effect given that no significant correlation between read number and allele number per individual was found. Moreover the consistency of the genotyping method and the variant validation criteria was confirmed by the high repeatability (0.94) of duplicates, and by the substantial increase in the repeatability measures between the 3rd and 5th step of variant validation procedure.
Still, of the 12 duplicates only six had identical genotypes following variant validation. Therefore it is plausible to suggest that the quality of the genotyping method could have been improved by modifying the variant validation procedure. As mentioned in the methods section, we tried to maintain a balance between the quality and size of the dataset whilst setting up the genotyping criteria; hence a better experimental design would have maximized both measures. The inconsistency among samples in six of the duplicates was mostly due to removal of real alleles, because the variants in question had low frequency on per individual basis and were highly similar to more common variants. Although the removal of rare alleles poses a problem for accurate genotyping, the application of supertype classification alleviated the issue to some extent, since highly similar alleles were grouped into the same supertype and the contribution of single alleles became irrelevant following supertyping. Hence, of the 12 duplicates eight had identical genotypes after supertyping, no correlation was found between the read number and supertype number per individual and the agreement between duplicates was calculated as 0.96. To date few Mhc studies carried out an empirical assessment of genotyping error via running duplicates in non-model vertebrates , therefore it is hard to assess how good the agreement is, but we believe that in the context of ecological studies, this level of repeatability is sufficient.
Four great tit cDNA libraries were constructed from mRNA sequences extracted from blood, and the libraries were used to identify transcribed allelic clusters. However a few cDNA sequences were retrieved from each allelic cluster, so it was not possible to isolate any group as more functionally important. Surprisingly, pseudogene alleles that bore stop codons were also found in the cDNA libraries. This emphasized the fact that classical sequencing methods are limited and even problematic if used for expression analysis, as they only inform on the presence or absence of transcripts. Such data are insufficient, because even non-functional genes could be transcribed in low levels unless a mutation at the promotor region halts the transcription completely. Therefore an approach based on expression level quantification (for instance using real-time PCR) should be adopted in order to differentiate highly expressed, hence functionally important alleles . Assessment of gene expression through cDNA coverage has started to find wider application in NGS methods .
Alternatively, the presence of non-functional cDNA sequences can also be explained by genomic DNA contamination during RNA extraction. However we believe this is unlikely, because the method used for RNA purification (RNeasy Protect Animal Blood) included a DNase digestion step. Although mRNA data did not inform us on allelic expression, the phylogenetic tree of the class I alleles and the historical selection tests significantly improved our understanding of functional alleles. The Group 1 sequences formed a monophyletic cluster with pseudogene alleles and codon based Z-test of selection suggested none of the putative antigen binding sites were positively selected. Still we detected a weak and non-significant signal for the positive selection of Group 1 ABS. Such weak signals are expected, as it has been shown to take around 19–74 million generations for a positive dN-dS signal to disappear even in the absence of selection . Using a likelihood ratio modelling approach three sites were found to be under positive selection in Group 1, yet these were not within highly variable sites. In contrast all the putative antigen-binding amino acids were under positive selection in Group 2, in addition to three highly variable sites. Moreover Z –test of selection suggested both ABS and PSS were under great selective pressure. These results implied that Group 2 sequences were actively involved in antigen recognition and have been under selection from a variety of parasites; whereas Group 1 sequences were non-functional, hence the observed variation possibly accumulated as a result of relaxed selection.
In total 755 alleles were considered functional and individual great tits possessed 9 to 32 putative expressed alleles. A discrepancy in the number of alleles per individual is common in Mhc studies across taxa, due to variation in gene copy number within species and allele-sharing between loci [24, 31, 61, 62]. We cannot exclude the possibility that some alleles might have been missed at the PCR stage, although we believe this to be relatively unlikely considering that we gathered extensive sequence information during characterization, and amplifications were performed using a degenerate primer pair. Similarly sized alleles did not form monophyletic groups implying that 3 bp insertions and deletions occurred several times in the history of the gene, independent of each other. The important contribution of indels to Mhc genomic diversity has already been shown in chicken . Moreover a comparative study between human and chimpanzee Mhc proposed indels as the major driving force for genomic divergence . Alternatively, this pattern might also be an effect of recombination and gene conversion, processes shown to be frequent in passerine Mhc systems [20, 23].
The alleles of Group 2 were clustered into 17 functional supertypes based on the physicochemical properties of their PSS. We adopted a bioinformatic approach that was described by Doytchinova and Flower  to determine peptide specificity; however an experimental approach would have been ideal. For instance, a recent study established computational antigen-binding prediction algorithms based on empirical datasets, and showed an evolutionary advantage for allele pairs that are more divergent in recognizing a broader range of potential antigens . Moreover the biological relevance of grouping Mhc alleles with similar antigen binding affinities into supertypes is supported by a growing number of human and non-human primate studies [33, 35, 66, 67]. In this study we utilized a newly proposed method for identifying genetic clusters (reviewed in ). This method is an advance on previous methods, as it does not require arbitrary clustering decisions , but uses K-means clustering algorithm and model selection approach to compute associated summary statistics.
We found evidence for the presence of at least 16 functional loci in the great tit. This constitutes the highest number of expressed Mhc class I alleles/loci identified in a passerine species. However we believe it is likely that this complexity is not specific to great tits and that similar patterns can be found in other passerines. Studies on great reed warbler and scarlet rosefinch have already revealed high polymorphism and the existence of more than 5–6 functional Mhc class I loci in these species, although these studies had lower sample sizes (248 and 120 respectively) and used conformation-based mutation detection methods that rely on physical separation of alleles, without providing allele sequence information [20, 43]. The utilization of indirect typing methods and motif-specific primers can be a rewarding approach in species where functional alleles are well described and confined to an allelic subset. However in species where the Mhc structure is complex and the number of co-amplifying alleles is high, like in the great tit, it is impossible to differentiate variants based on their migratory patterns and it is inevitable to underestimate the allelic diversity as a result of peak overlapping . Moreover, because indirect typing methods require cloning-sequencing effort to reveal the nucleotide information, hundreds of clones must be sequenced to obtain a good estimate of allelic composition in such complex systems . Therefore it is feasible to suggest that the employment of 454-pyrosequencing is ideal for passerine species that harbour many Mhc loci and high allelic diversity, as it has the potential to reveal the extent of Mhc complexity. In a recent study 454 technology was used for Mhc genotyping Atlantic cod and it was shown that the species possess more than 100 Mhc class I loci, a greatly expanded number compared to other teleost fish . This extraordinary expansion of class I genes was explained by the absence of class II loci; this might represent a compensatory mechanism adapted by the Atlantic cod immune system. This work illustrates the great flexibility in Mhc genomic organization among closely related species.