Spatial and temporal variation at major histocompatibility complex class IIB genes in the endangered Blakiston’s fish owl

Introduction Quantifying intraspecific genetic variation in functionally important genes, such as those of the major histocompatibility complex (MHC), is important in the establishment of conservation plans for endangered species. The MHC genes play a crucial role in the vertebrate immune system and generally show high levels of diversity, which is likely due to pathogen-driven balancing selection. The endangered Blakiston’s fish owl (Bubo blakistoni) has suffered marked population declines on Hokkaido Island, Japan, during the past several decades due to human-induced habitat loss and fragmentation. We investigated the spatial and temporal patterns of genetic diversity in MHC class IIβ genes in Blakiston’s fish owl, using massively parallel pyrosequencing. Results We found that the Blakiston’s fish owl genome contains at least eight MHC class IIβ loci, indicating recent gene duplications. An analysis of sequence polymorphism provided evidence that balancing selection acted in the past. The level of MHC variation, however, was low in the current fish owl populations in Hokkaido: only 19 alleles were identified from 174 individuals. We detected considerable spatial differences in MHC diversity among the geographically isolated populations. We also detected a decline of MHC diversity in some local populations during the past decades. Conclusions Our study demonstrated that the current spatial patterns of MHC variation in Blakiston’s fish owl populations have been shaped by loss of variation due to the decline and fragmentation of populations, and that the short-term effects of genetic drift have counteracted the long-term effects of balancing selection. Electronic supplementary material The online version of this article (doi:10.1186/s40851-015-0013-4) contains supplementary material, which is available to authorized users.


Introduction
Blakiston's fish owl (Bubo blakistoni Seebohm; hereafter 'fish owl'), the world's largest owl, is currently categorized as 'Endangered' on the IUCN Red List of Threatened Species [1]. This owl is non-migratory and endemic to northeastern Asia, comprising two subspecies: B. b. blakistoni in central and eastern Hokkaido and the southern Kuril islands; B. b. doerriesi in continental Asia, including the Russian Far East, northeastern China, and probably North Korea (Figure 1) [2]. Unlike many other owl species, the fish owl is specialized to aquatic prey, mainly freshwater fishes, and thus its habitat is limited to areas close to lakes, rivers, springs, and shoals that do not freeze in winter. In addition, it requires large, old trees with large cavities for nesting [2][3][4]. Over the last several decades, riparian old-growth forests, the preferred habitat for the fish owl, have been decimated due to human land use, causing a rapid decline in fish owl populations.
The current population size of the fish owl is estimated to be 120-150 individuals in Hokkaido, 70-85 in the southern Kuril Islands, and a few thousand in continental Asia [2,3,5]. The fish owl was widespread in Hokkaido until the middle of the 20th century [3,6], but its population has decreased markedly over the last several decades due to human-induced habitat loss and fragmentation, and probably fell to fewer than 100 individuals in the 1970-80s [7]. A project began in the early 1980s aimed at conservation of the fish owl, through the installation of nest boxes and artificial feeding. Thanks to these efforts, the Hokkaido population is now gradually recovering, but the risk of extinction remains high due to the loss of adaptive genetic variation and inbreeding depression in highly fragmented populations.
Recent population genetic studies using selectively neutral markers such as mitochondrial DNA (mtDNA) sequences and microsatellite genotypes revealed low overall genetic diversity in Hokkaido fish owl populations [8,9]. These studies also detected genetic differentiation among fragmented populations, indicating limited movement of owls between areas. In contrast, mtDNA sequence analyses of historical specimens, including taxidermied specimens and archaeological bones, indicated that gene flow had occurred over a broad area of Hokkaido until the middle of the 20th century [9]. Both the mtDNA and microsatellite data showed a sharp reduction in genetic diversity after around 1980 [9]. The low levels of genetic variation and genetic differentiation among subpopulations currently inhabiting Hokkaido are probably due to genetic drift and inbreeding resulting from habitat loss and fragmentation.
The major histocompatibility complex (MHC) is a multigene family that plays a crucial role in the vertebrate immune system. It contains genes coding for cellsurface glycoproteins, the MHC class I and II molecules, each of which is specifically involved in presenting antigen peptides derived from intra-or extracellular pathogens to T cells, initiating the adaptive immune response [10]. These genes are among the most polymorphic in the vertebrate genome, exhibiting high allelic diversity [11]. MHC polymorphism is generated by frequent gene duplications and deletions, intra-and inter-locus recombination or gene conversion, and the accumulation of de novo mutations [12,13]. In addition, the enormous variation in MHC genes is probably maintained by pathogen-driven balancing selection, mediated through either heterozygote advantage or negative frequencydependent selection, and by sexual selection via MHCmediated disassortative mating (reviewed in [14][15][16][17]).
Although balancing selection shapes MHC polymorphism to a great extent over the long term, MHC variation is often substantially reduced in populations that have undergone extreme bottlenecks (e.g., [18][19][20][21]). Habitat fragmentation may further accelerate the reduction in MHC diversity by strengthening the effects of genetic drift (e.g., [22][23][24]). Although the effects of reduced MHC diversity on the long-term viability of populations that undergo bottlenecks has remained unclear (reviewed in [25]), the example of contagious cancer in Tasmanian devils (Sarcophilus harrisii) strongly supports the possibility that populations with low MHC diversity are more vulnerable to outbreaks of novel infectious diseases (reviewed in [26]). Quantifying MHC diversity can thus provide an indirect measure of the immunological fitness (i.e. the potential of resistance to novel infectious diseases) of a population, and should thus be incorporated into studies of endangered species aimed at establishing conservation plans [27,28].
The first comprehensive information on the genomic structure of the MHC in birds came from studies of the domestic chicken, Gallus gallus [29]. The chicken MHC  is more compact than the mammalian MHC; it has only two classical MHC class I and II genes, which is referred as 'the minimal essential MHC' [29]. This compact MHC structure has been reported in various bird lineages (e.g., [30][31][32][33][34][35]), but is not universal in birds. The Japanese quail (Coturnix japonica) appears to have a more complex MHC structure [36] than the chicken, even though the two belong to the same family. A highly complex MHC structure has been reported in passerine species, which commonly have many copies of both class I and II MHC genes and a large number of pseudogenes (e.g., [37][38][39][40]). The considerable variation among bird taxa suggests that the avian MHC structure is evolutionarily labile, probably due to recent gene duplications and pseudogene formation [41].
In the present study, we investigated the genetic diversity of the MHC class IIβ loci in Hokkaido fish owl populations, based on samples collected from 1963 to 2012. We used massively parallel pyrosequencing [42,43] for exhaustive genotyping of the fish owl MHC class IIβ loci. Our aims were (1) to describe the polymorphism at second exon of MHC class IIβ genes in the fish owl, and (2) to elucidate spatial and temporal patterns in the MHC diversity in order to assess how the recent population decline associated with habitat loss and fragmentation have affected variation in these functional genes.  Figure 1). Most samples were from firstyear juveniles, though some were from adults. All blood samples were non-invasively collected from owls by a veterinarian in the activity for conservation of the Blakiston's fish owl by the Ministry of the Environment, Japan. Some drops of the blood samples were preserved in 99.9% ethanol or dried on filter paper, and stored at −20°C. Fibroblasts were cultured from skin tissue samples according to the method described in Nishida et al. [44]. Total genomic DNA was extracted from blood and fibroblasts with the DNeasy Blood & Tissue Kit (Qiagen).

RNA extraction and cDNA synthesis
For expression analysis, total RNA was extracted with the RNeasy Mini Kit (Qiagen) from skin fibroblasts from 16 Hokkaido individuals. First-strand cDNA was synthesized from 1 μg of total RNA using the oligo-(dT) 20 primer with the Superscript III First-strand Synthesis System (Invitrogen). The expression of MHC class II molecules in skin fibroblasts was confirmed through polymerase chain reaction (PCR) using the primers Sf-ex1F (5′-CAC TGG TGG TGC TGG GAG CC-3′) and Tyal-ex3-3′R (5′-AGG CTG ACG TGC TCC ACC TG-3′), which bind in conserved regions of exons 1 and 3 of owl MHC class II loci [32,45].

Primer design
To develop specific primers for pyrosequencing, complete and partial sequences of the second exon of the MHC class IIβ loci were amplified from several individuals from different populations in Hokkaido and from different years by using the primer sets Sf-ex1F and Tyal-ex3-3′R, and Stri2FC (5′-CMC ACA CAG GGG TTT TCC-3′) and Stri2RC (5′-AAC GYG YGG CCA CGC GCT CA −3′) [30]. Six cDNA and 10 genomic DNA samples were used as PCR templates. PCR amplifications were performed in 25-μl reaction volumes, each containing 0.5 units of high-fidelity PrimeSTAR GXL DNA Polymerase (TaKaRa Bio), 1 × PrimeSTAR GXL buffer, 200 μM each dNTP, 1 μM each primer, and approximately 200 ng of cDNA or 50 ng of genomic DNA. Cycling conditions were 98°C/1 min, 28 × (98°C/ 10 s, 68°C/30 s), 68°C/3 min for Sf-ex1F and Tyal-ex3-3′R; and 98°C/1 min, 28 × (98°C/10 s, 58°C/10 s, 68°C/30 s), 68°C/3 min for Stri2FC and Stri2RC. Because the genomic sequences of the MHC class IIβ region are extraordinarily GC-rich in avian species, we used high annealing temperature to increase the efficiency of amplifications as recommended in Burri et al. [46]. PCR products were cloned using the Zero Blunt TOPO PCR Cloning Kit (Invitrogen), and 24-48 clones per sample were sequenced with the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) and an ABI 3730 DNA Analyzer (Applied Biosystems). From consensus sequences for the sequences obtained, specific primers were designed to amplify part of the second exon of the fish owl MHC class IIβ loci.

Preparation of an amplicon library and pyrosequencing
Pyrosequencing was performed on 242 amplicons obtained from 200 genomic DNA samples and 16 cDNA samples. To generate an amplicon library, fusion primers were designed that contained the GS FLX Titanium primer sequence (A in forward, B in reverse primers) at the 5′ end, followed by a 10-bp multiplex identifier (MID) sequence and the sequence of a gene-specific primer. The MID sequences were chosen from the Extended MID Set (Roche) and were used to distinguish amplicons obtained from different PCR reactions. For genomic DNA samples, newly designed primers BublIIb2F (5′-GAG TGT CAG YAC CTY RAY RG-3′) and BublIIb2R (5′-CTT TCY TCT SCS TGA YGW AGG-3′) were used as forward and reverse gene-specific primers to amplify 203-bp fragments of the second exon of MHC class IIβ loci ([see Additional file 1: Figure S1]). Ten genomic DNA samples were amplified and sequenced twice to estimate the genotyping error. For cDNA samples, we used 2 reverse primers (BublIIb2-3R1, 5′-TTC CAC CTC GGG CGG GAC TTT C-3′; BublIIb2-3R2, 5′-CCT CAC CTT GGG CTG AAC TTT C-3′) that spanned the exon-exon junction between the second and third exons. PCR was conducted twice for each cDNA sample to amplify 213-bp fragments, using the fusion primer pairs containing the sequences of BublIIb2F/BublIIb2-3R1 and BublIIb2F/ BublIIb2-3R2.

MHC genotyping
SFF Tools 2.8 (Roche) was used to assign the processed reads to respective individuals based on the MID sequences in the forward and reverse fusion primers, and seq_crumbs 0.1.8 (available from http://bioinf.comav.upv.es/seq_crumbs/) was used to trim and filter sequence reads based on quality and sequence length.
Allele detection and genotyping were performed with the recently developed pipeline ngs_genotyping 0.9.0 ( [47]; available from https://github.com/enormandeau/ ngs_genotyping). This pipeline uses an iterative procedure in three successive steps. The first step generates putative allele sequences for each individual. The second step combines and strengthens these into the global alleles for all individuals. The third step then genotypes each individual. Prior to these steps, sequence reads were iteratively cleaned and aligned by means of the MUSCLE algorithm [48]. Throughout the text, unique sequence variants are referred to as 'alleles' for convenience, although this is not strictly correct, as these sequences represent multiple loci.
The following parameters were used with a hierarchical clustering analysis in the first step of the pipeline to filter out sequencing errors and to detect putative alleles: minimum internal branch length, 0.06; minimal proportion threshold to define the cluster, 0.02. Likewise, the following parameters were used to detect global consensus alleles in the second step: minimum internal branch length, 0.06; minimum number to define the cluster, two sequences. In the third step, the number of sequence reads of each global allele was counted for each individual by using the BLASTn algorithm. Each individual was then genotyped according to the minimal proportion threshold of each allele. Threshold values ranging from 0.01 to 0.05 were used, depending on natural breaks in the number of reads per individual for each allele. The above settings of the threshold values in the pipeline meet the two-PCR criterion that is standard in MHC studies [49].

Data analysis
Phylogenetic relationships among the fish owl MHC class IIβ alleles were reconstructed by Bayesian inference and maximum likelihood (ML) analyses, using the programs MrBayes 3.2.2 [50] and GARLI 2.01 [51]. The best-fit model of nucleotide substitution was selected based on the Bayesian information criterion implemented in jModelTest 2.1.1 [52]. MrBayes analyses were performed with two parallel runs of 20 million generations each and using one cold and three heated Markov chains, with sampling every 1000 steps. Convergence and stationarity of the chains were confirmed by the average standard deviation of split frequencies (<0.01). The first 25% of trees were discarded as burn-in, and a 50% majority-rule consensus tree was constructed from the remaining trees. An ML bootstrap search with 1000 pseudo-replicates were perfomed with GARLI. For these analyses, sequences of the MHC class IIβ alleles of northern goshawk (Accipiter gentilis; EF370917) and Eurasian black vulture (Aegypius manachus; EF370890) were also included as outgroups. An additional analysis was conducted to reconstruct the relationships between MHC class IIβ alleles identified in the fish owl and those in other owl species (GenBank accession nos. EF370927-370928, EF370930-370946, and EF641223-EF641262 [30,32]).
A Z-test for selection was performed to detect the signature of historical selection on the fish owl MHC class IIβ loci by using MEGA 5.2.2 [53] with the modified Nei-Gojobori method and Jukes-Canter correction. The average rates of synonymous (d S ) and non-synonymous substitution (d N ) were calculated for all sites, as well as for positions encoding amino acids in the putative peptidebinding region (PBR), and the remaining positions (putative non-PBR). The location of the PBR was inferred from the molecular structure of human leukocyte antigen class II, HLA-DR1 [54]. In addition, Bayesian inference was performed with omegaMap 0.5 [55] to detect amino acid sites under positive selection. This program uses a coalescentbased model for detecting natural selection in the presence of recombination, and uses the reversible-jump Markov chain Monte Carlo to perform Bayesian inference of both the d N /d S ratio (ω) and the recombination rate (ρ). The omega_model was set to independent, with the rho_model set to variable and the rho_block set to 3. To reduce computation time, two independent random subsamples of 200 alleles each were generated from all populations. Three independent runs of 2 million steps each were performed for each subsample. Convergence and stationarity of the runs were assessed by calculating the effective sample size, which exceeded 200 for all estimated parameters. The first 50,000 steps were discarded as burnin. Codons were considered to be under positive selection if the posterior probability of ω > 1 exceeded 0.95 in both independent subsample run sets.
The detection of potential recombinant sequences in our data set was carried out using a set of seven nonparametric detection methods implemented in RDP4 beta 4.27 software [56]: RDP, GENECONV, MaxChi, Chimaera, BootScan, SiScan, and 3Seq. The analysis was performed with default settings for the various detection methods, and the Bonferroni-corrected P value cutoff was set at 0.05. Recombination events were accepted when detected with at least three of the seven detection methods. Additionally, the web-based service GARD (genetic algorithm for recombination detection) [57] was used to detect recombination breakpoints.
Population-level allelic richness was calculated through 1000 bootstrap replicates with a constant sample size (n = 4). Nucleotide diversity (π) within individuals and populations were calculated using MEGA. To assess the MHC diversity within and between populations, R 3.0.1 [58] with the package ecodist 1.2.7 [59] was used to calculate Jaccard distances from a binary presence/absence matrix for each allele in pairwise comparisons of individuals. Jaccard dintance is the measure of dissimilarity between sample sets, and is defied as one minus ratio of the size of the intersection and union of the sample sets [60]. Non-metric multidimensional scaling (NMDS) was performed to visualize the relationships among individual MHC genotypes based on Jaccard distances.
To assess the genetic differentiation among populations, global and population pairwise G ST values [61] were estimated from allele frequencies, using the mmod package version 1.2.1 [62] in R. The significance of the G ST values were evaluated with the permutation tests (1000 replicates). Finally, an analysis of molecular variance (AMOVA) [63] was conducted to evaluate the spatial and temporal patterns of genetic variations, using the ade4 package version 1.6-2 [64] in R. An AMOVA was performed with two different genetic distance matrixes calculated in pairwise comparison of individuals: 1) Jaccard distance, and 2) average number of nucleotide substitutions per site between individuals, calculated with the Kimura 2-parameter model. The total genetic variance was partitioned into three hierarchical levels: among populations, between periods within populations, and within periods. The significance of the variance components were assessed with the permutation tests (1000 replicates).

Genotyping of fish owl MHC class IΙβ
Amplification primers with the complete MID sequences were identified in 193,393 reads. Final genotypes were based on 139,832 reads; the mean coverage (number of reads per amplicon) was 577.8 reads (SD = 307.7, range = 34-1851). For conservative genotyping, individuals with fewer than 200 reads were excluded from the analysis. No correlation between the number of reads and the number of alleles observed per individual was detected in the remaining samples ([see Additional file 1: Figure S2]), indicating that the coverage was sufficient for reliable genotyping. Of the 10 replicated pairs that were run twice to estimate genotyping error, nine pairs had > 200 reads for both replicates, and the same alleles were found in these pairs. In all, 174 individuals were genotyped (Table 1; Additional file 2).

Phylogenetic relationships among MHC alleles and allelic expression patterns
From the partial sequences (203 bp) of MHC class IIβ exon 2 among 174 fish owls, 19 unique variants (alleles) were identified. We named these alleles Bubl-DAB*01-19 (only the abbreviation Bubl is used hereafter), following the nomenclature proposed by Klein et al. [65]. Sequences of all alleles were deposited in the DNA Data Bank of Japan (DDBJ) under accession nos. LC007937-LC007955 (sequence alignment in Additional file 3).
The alleles we detected showed high degrees of nucleotide (94/203 sites variable) and amino acid (43/67 sites variable) polymorphism. The number of alleles per individual ranged from eight to 16, indicating there are at least eight MHC class IIβ loci in this species. Nine of 19 alleles were also found in 16 cDNA samples (Figure 2 and [see Additional file 1: Table S1]). Up to eight alleles were detected in cDNA samples from single individuals, indicating that at least four loci are expressed in skin fibroblasts.
Additional file 1: Figure S3 shows the phylogenetic relationships among MHC class IIβ alleles from the fish owl and other owls. This pattern is one of trans-species polymorphism [66]: alleles in one species are most closely related to alleles in other species, rather than grouping by species. No alleles were shared between the fish owl and other owl species based on the nucleotide sequences or amino acid sequeneces.

Patterns of selection and recombination
Because the allele Bubl12 contains an internal stop codon ([see Additional file 1: Figure S4]), this allele was excluded from the Z-test and the omegaMap analysis. A highly significant excess of non-synonymous over synonymous substitutions was observed for codons in the putative peptide-binding regions (PBR), whereas synonymous substitutions were more frequent than non-synonymous substitutions in non-PBR codons ( Table 2). The similar patterns were obtained in the analyses perfomed only with expressed alleles and that only with unexpressed alleles, although more synonymous and non-synonymous mutations   were detected in the comparisons between the unexpressed alleles than between the expressed alleles ( Table 2). The omegaMap analysis indicated that 13 of 67 amino acid positions are under balancing selection ([see Additional file 1: Figure S5]). These positions lined up with 10 of 19 putative PBR residues defined by alignment to human HLA-DR1 ([see Additional file 1: Figure S4]). In addition, 2 of the putative non-PBR residues had a signature of balancing selection. No evidence of recombination were found within exon 2 sequences of the fish owl MHC class IIB alleles in either RDP or GARD analyses.

Spatial and temporal patterns of MHC diversity
The variablity in the number of alleles per individual differed among populations (Figure 3), although there is little difference in the median number of alleles per indiviual among populations ( Table 1). The SR population shows highest allelic richness in the all periods (Table 1), and all of the 19 alleles were observed in this population during the period of 2003-2012. A relatively high alleric richness was detected in HD in the period of 2003-2012 (Table 1). Although DS shows moderate allelic richness (Table 1), alleles Bubl16 and Bubl17, both of which were found before 2002, were not found thereafter. Conversely, Bubl19 was found only in individuals sampled after 2003 (Figure 2). The estimates of alleleic richness were decreased in KS and AK in the last several decades, and only 11 or 12 alleles were observed in these populations in the period of 2002-2012 (Table 1). Some alleles whose frequencies were initially low were probably lost in these populations ( Figure 2). Meanwile, nucleotide diversities were relateively higher in the KS and AK populations, indicating that remained alleles in these populations were highly diverged from each other. Interestingly, these alleles were also maintained in high frequencies in the other populations ( Figure 2). Figure 4A shows the distribution of Jaccard distances among individuals within populations, calculated from samples collected from 1993 to 2012. Mean Jaccard values were 0.244 (SR), 0.023 (KS), 0.123 (AK), 0.188 (DS), 0.245 (HD), and 0.220 (All). The peaks of the distributions of Jaccard distance were near zero for KS and AK ( Figure 4A), indicating that most individuals in these populations had a similar set of alleles. An NMDS plot shows that genotypic diversity was lower in these populations than in the other populations ( Figure 4B). No significant genetic differences were detected among populations in the period of 1963-1992 (global G ST = −0.0014,   Hierachical analysis by AMOVA revealed that substantial amount of genetic variation was attiributed to differences among populations (Table 3). An AMOVA also detected weak but partially significant genetic differences among periods whitin populations (Table 3). Interestingly, AMOVA revealed a high proportion of variance explained between populations, while G ST values were very small. This may be because G ST only considers allele frequencies, and maximal G ST is reduced for highly variable markers such as typically microsatellites, but also MHC [67].

MHC class IIβ genes in Blakiston's fish owl
Our results show that the Blakiston's fish owl genome contains at least eight MHC class IIβ loci, among which at least four are expressed in skin fibroblasts. Although the number of individuals used in the expression analysis was insufficient to confirm the expression status of all alleles, some alleles isolated from genomic DNA samples were not detected in the cDNA samples from the same individuals ( Figure 2, [see Additional file 1: Table S1]), suggesting that these alleles are pseudogenes or have different functions. Indeed, one of these alleles is clearly non-functional, as it contains an internal stop codon. The other alleles not detected in the cDNA samples are possibly expressed at very low levels or in tissues other than skin fibroblasts. Further study is required to clarify the actual number of gene copies and expression status of the fish owl MHC class IIβ genes.  The evolution of the genomic structure of avian MHC is of particular interest, as it varies considerably among taxa [41]. To our knowledge, the number of MHC class IIβ genes in the fish owl suggested by our results is the greatest in any non-passerine bird. Besides the fish owl, Burri et al. [45] reported that barn owl (Tyto alba) has two duplicated MHC class IIβ genes, by means of a Southern blot analysis. Although the number of MHC class IIβ genes has not yet been investigated in detail in other owls, the cDNA based study revealed that the strigiform owls have two or more gene copies [32,68]. Our results suggest significant variation in the number of gene copies in owl MHC class IIβ. It is unknown why the fish owl has a large number of MHC class IIβ genes. Although the reasons for the taxonomic trends in avian MHC structure are yet unclear, gene duplications may have occurred at a higher frequency in some lineages than others due to the activity of endogenous retroviral elements, as has been suggested for primate [69], passerine [38], and wallaby MHC genes [70]. Future comparative genomic studies of the fish owl MHC structure using more powerful approaches, such as whole-genome sequencing or screening BAC libraries, will provide new insights into the evolutionary history of the avian MHC. From an eco-immunological point of view, the complexity of MHC genes may be related to the complexity of pathogens in the environment. As the fish owl specializes on aquatic prey, it may be exposed to different suites of pathogens than most other owls, which feed on terrestrial prey. However, little is known about pathogens infecting the fish owl, and further study is necessary.
The Z-test and the omegaMap analysis revealed that the ratio of non-synonymous to synonymous substitutions was significantly higher than expected in PBRs, based on the neutral model, but not in non-PBRs (Table 2; [see Additional file 1: Figure S5]), a pattern consistent with balancing selection [71,72]. The pattern of trans-species polymorphism further supports the conclusion that balancing selection acted in the past ([see Additional file 1: Figure S3]). These results are concordant with most other studies reporting selective signatures for MHC genes (reviewed in [14]).

Effect of population reduction and fragmentation on fish owl MHC diversity
Although balancing selection might maintain MHC polymorphisms over the long term, strong genetic drift in populations that undergo a bottleneck and fragmentation can counteract the effects of balancing selection [25]. Our results indicate that a recent habitat loss and fragmentation has substantially lowered MHC diversity in the fish owl. Using massively parallel pyrosequencing, we identified 19 alleles among 174 individuals from Hokkaido populations. Assuming that there are at least eight loci, the level of genetic variation in MHC class IIβ loci is low in the fish owl populations in Hokkaido, as other species often show higher levels of MHC class IIβ allelic diversity. For example, Alcaide and colleagues [73] detected 103 alleles at a single locus among 121 individuals of the lesser kestrel (Falco naumanni). High levels of allelic diversity at MHC genes in wild bird populations were also reported in the recent studies using next-generation sequencing approach [39,40,74]. The low variation detected in the fish owl populations was not a consequence of limited sampling, as our population sample covered most families currently living in Hokkaido.
Our results also revealed that the MHC diversity differed among geographically isolated populations in Hokkaido. The highest allelic diversity was observed in the SR population (Table 1), in which all 19 alleles identified from all Hokkaido populations combined were observed during the last decade (Table 1 and Figure 2). NMDS analysis also detected high inter-individulal MHC diversity within SR populaion ( Figure 4). These results correlate with the large size of the SR population; approximately half the fish owls in Hokkaido inhabit the Shiretoko Peninsula [3]. A microsatellite analysis also detected relatively high genetic diversity in this population [8]. A moderate level of MHC diversity was observed in the DS population, which is the second largest Hokkaido. MHC diversity was also relatively high in the HD population, despite its small size. Interestingly, this is not congruent with microsatellite data [8,9], where the HD population showed the lowest genetic diversity. This discrepancy suggests that analyses of neutral genetic markers alone inadequately quantify genetic variation. In contrast, inter-individulal MHC diversity was very low in the AK population and nearly lacking in the KS population ( Figure 4). In the past two decades the number of alleles per individual in these populations was as high as 11 or 12 (Figure 3), which corresponded to the total number of alleles observed in these populations ( Table 1). Assuming that there are at least eight loci, this result indicates that more than half the MHC class IIβ loci are homozygous in individuals in the AK and KS populations. Both populations are relatively small in size, and during the 1980s there may have been as few as one to several breeding pairs in each population. As we detected more alleles in individuals in these populations sampled before 1992 (Table 1 and Figure 2), the currently low MHC variation appears to be due to genetic drift associated with recent habitat loss and fragmentation.
Finaly, we did not investigate samples from a continental population in the present study. As the fish owl habitat is better preserved and population sizes are larger on the continent than in Hokkaido, higher MHC diversity may have been maintained in these population.
Further analysis on a continental population may lead to a better understanding of the extent of MHC diversity in this species.

Conclusions
Our study demonstrated substantial spatial variation in MHC diversity in the current fish owl populations on Hokkaido. Analysis of mtDNA sequences from historical specimens suggested that the genetic diversity was higher before the 1980s, and that before the middle of 20th century, there was gene flow over broad areas of Hokkaido [9]. The current spatial pattern in MHC diversity may thus have been shaped by loss of variation due to genetic drift in the fragmented populations. Low MHC variation among individuals, as in the KS and AK populations, should make these populations more vulnerable to epidemics of novel pathogens. The effects of low MHC diversity on the viability of fish owl populations remains unclear, however, as information on the abundance of pathogens and prevalence of diseases in wild fish-owl populations is scarce. Based on the results of recent population genetic studies, genetic rescue of the fish owl populations by translocation of outbred individuals has been proposed, and the results of our study will be useful in establishing such translocation plans. Although translocations could increase the adaptive genetic variation, a detailed screening of pathogens should be done before any attempts at mixing different populations are considered.