Spatial and temporal variation at major histocompatibility complex class IIB genes in the endangered Blakiston’s fish owl
Zoological Letters volume 1, Article number: 13 (2015)
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.
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.
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.
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 . 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) . 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-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 . 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 . Both the mtDNA and microsatellite data showed a sharp reduction in genetic diversity after around 1980 . 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 cell-surface 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 . These genes are among the most polymorphic in the vertebrate genome, exhibiting high allelic diversity . 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 frequency-dependent selection, and by sexual selection via MHC-mediated disassortative mating (reviewed in [14-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-21]). Habitat fragmentation may further accelerate the reduction in MHC diversity by strengthening the effects of genetic drift (e.g., [22-24]). Although the effects of reduced MHC diversity on the long-term viability of populations that undergo bottlenecks has remained unclear (reviewed in ), 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 ). 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 . 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’ . This compact MHC structure has been reported in various bird lineages (e.g., [30-35]), but is not universal in birds. The Japanese quail (Coturnix japonica) appears to have a more complex MHC structure  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-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 .
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.
Materials and methods
Sampling and DNA extraction
From 1963 to 2012, blood or skin tissue samples were collected from 200 individuals of the fish owl from five geographically isolated populations across Hokkaido: Shiretoko (SR), Konsen (KS), Akan (AK), Daisetsu (DS), and Hidaka (HD) populations (Figure 1). Most samples were from first-year 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. . 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].
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′) . 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. . 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.
PCR conditions were 98°C/1 min, 28 × (98°C/10 s, 58°C/10 s, 68°C/30 s), and 68°C/3 min. PCR amplicons were purified by using the MinElute PCR Purification Kit (Qiagen). After quantification by agarose-gel electrophoresis, the purified amplicons were pooled in approximately equimolar quantities. The amplicon library was commercially sequenced on a 1/4 Titanium Pico Tire Plate with the GS FLX Titanium Sequencing Kit XLR70 (Roche) at Hokkaido System Science Co. (Sapporo, Japan).
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 (; 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 . 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 .
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  and GARLI 2.01 . The best-fit model of nucleotide substitution was selected based on the Bayesian information criterion implemented in jModelTest 2.1.1 . 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  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 peptide-binding 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 . In addition, Bayesian inference was performed with omegaMap 0.5  to detect amino acid sites under positive selection. This program uses a coalescent-based 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 burn-in. 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 non-parametric detection methods implemented in RDP4 beta 4.27 software : 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)  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  with the package ecodist 1.2.7  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 . 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  were estimated from allele frequencies, using the mmod package version 1.2.1  in R. The significance of the G ST values were evaluated with the permutation tests (1000 replicates). Finally, an analysis of molecular variance (AMOVA)  was conducted to evaluate the spatial and temporal patterns of genetic variations, using the ade4 package version 1.6-2  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. . 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 : 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, P = 0.949; median of pairwise G ST = −0.0012 [see Additional file 1: Table S2]) or in the period of 1993–2002 (global G ST = −0.0010, P = 0.927; median of pairwise G ST = −0.0002), whereas slight but significant genetic differences were detected among poplations in the period of 2003–2012 (global G ST = 0.0012, P = 0.039; median of pairwise G ST = 0.0022). 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 .
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 . 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.  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 , passerine , and wallaby MHC genes . 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 ).
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 . 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  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 . A microsatellite analysis also detected relatively high genetic diversity in this population . 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.
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 . 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.
Availability of supporting data
Sequences of the MHC class IIβ alleles are available in GenBank/EMBL/DDBJ accession nos. LC007937–LC007955.
Sequence alignment and genotyping data are included in the additional files.
BirdLife International 2013. Bubo blakistoni. The IUCN red list of threatened species. version 2014.2. [http://www.iucnredlist.org]
Slaght JC, Surmach SG. Biology and conservation of Blakiston’s fish owls (Ketupa balakistoni) in Russia: a review of the primary literature and assessment of the secondary literature. J Raptor Res. 2008;42:29–37.
Takenaka T. Distribution, habitat environments, and reasons for reduction of the endangered Blakiston’s fish owl in Hokkaido, Japan. PhD thesis. Sapporo: Hokkaido University, Department of Environmental Earth Science; 1998.
Slaght JC, Surmach SG, Gutiérrez RJ. Riparian old-growth forests provide critical nesting and foraging habitat for Blakiston’s fish owl Bubo blakistoni in Russia. Oryx. 2013;47:553–60.
Takenaka T. Shimafukuro. In: Shiretoko Museum, editor. Birds in Shiretoko. Sapporo: The Hokkaido Shinbun press; 1999. p. 80–125.
Hayashi Y. Past and present distribution of Blakiston’s fish owl (Ketupa blakistoni) in Hokkaido, Japan –based upon museum specimens–. J Yamashina Inst Ornithol. 1999;31:45–61.
Brazil MA, Yamamoto S. The status and distribution of owls in in Japan. In: Meyburg B-U, Chancellor RD, editors. Raptors in the Modern World. Berlin, Germany: World Working Group on Birds of Prey; 1989. p. 389–402.
Omote K, Nishida C, Takenaka T, Masuda R. Temporal changes of genetic population structure and diversity in the endangered Blakiston’s fish owl (Bubo blakistoni) on Hokkaido Island, Japan, revealed by microsatellite analysis. Zoolog Sci. 2012;29:299–304.
Omote K, Nishida C, Takenaka T, Saito K, Shimura R, Fujimoto S, Sato T, and Masuda R. Recent fragmentation of the endangered Blakiston’s fish owl (Bubo blakistoni) population on Hokkaido Island, northern Japan, revealed by mitochondrial DNA and microsatellite analyses. Zoolog Lett. in press.
Klein J. Natural History of the Major Histocompatibility Complex. New York: John Wiley & Sons; 1986.
Garrigan D, Hedrick PW. Perspective: detecting adaptive molecular polymorphism: lessons from the MHC. Evolution. 2003;57:1707–22.
Ohta T. Role of diversifying selection and gene conversion in evolution of major histocompatibility complex loci. Proc Natl Acad Sci U S A. 1991;88:6716–20.
Nei M, Rooney AP. Concerted and birth-and-death evolution of multigene families. Annu Rev Genet. 2005;39:121–52.
Bernatchez L, Landry C. MHC studies in nonmodel vertebrates: what have we learned about natural selection in 15 years? J Evol Biol. 2003;16:363–77.
Sommer S. The importance of immune gene variability (MHC) in evolutionary ecology and conservation. Front Zool. 2005;2:16.
Piertney SB, Oliver MK. The evolutionary ecology of the major histocompatibility complex. Heredity (Edinb). 2006;96:7–21.
Spurgin LG, Richardson DS. How pathogens drive genetic diversity: MHC, mechanisms and misunderstandings. Proc R Soc B Biol Sci. 2010;277:979–88.
Bollmer JL, Vargas FH, Parker PG. Low MHC variation in the endangered Galápagos penguin (Spheniscus mendiculus). Immunogenetics. 2007;59:593–602.
Zhu L, Ruan X-D, Ge Y-F, Wan Q-H, Fang S-G. Low major histocompatibility complex class II DQA diversity in the Giant Panda (Ailuropoda melanoleuca). BMC Genet. 2007;8:29.
Castro-Prieto A, Wachter B, Sommer S. Cheetah paradigm revisited: MHC diversity in the world’s largest free-ranging population. Mol Biol Evol. 2011;28:1455–68.
Babik W, Kawalko A, Wójcik JM, Radwan J. Low major histocompatibility complex class I (MHC I) variation in the European bison (Bison bonasus). J Hered. 2012;103:349–59.
Biedrzycka A, Radwan J. Population fragmentation and major histocompatibility complex variation in the spotted suslik, Spermophilus suslicus. Mol Ecol. 2008;17:4801–11.
Říčanová Š, Bryja J, Cosson J-F, Gedeon C, Choleva L, Ambros M, et al. Depleted genetic variation of the European ground squirrel in Central Europe in both microsatellites and the major histocompatibility complex gene: implications for conservation. Conserv Genet. 2011;12:1115–29.
Strand TM, Segelbacher G, Quintela M, Xiao L, Axelsson T, Höglund J. Can balancing selection on MHC loci counteract genetic drift in small fragmented populations of black grouse? Ecol Evol. 2012;2:341–53.
Radwan J, Biedrzycka A, Babik W. Does reduced MHC diversity decrease viability of vertebrate populations? Biol Conserv. 2010;143:537–44.
Belov K. Contagious cancer: lessons from the devil and the dog. BioEssays. 2012;34:285–92.
Hughes AL. MHC polymorphism and the design of captive breeding programs. Conserv Biol. 1991;5:249–51.
Ujvari B, Belov K. Major histocompatibility complex (MHC) markers in conservation biology. Int J Mol Sci. 2011;12:5168–86.
Kaufman J, Milne S, Göbel TWF, Walker BA, Jacob JP, Auffray C, et al. The chicken B locus is a minimal essential major histocompatibility complex. Nature. 1999;401:923–5.
Alcaide M, Edwards SV, Negro JJ. Characterization, polymorphism, and evolution of MHC class II B genes in birds of prey. J Mol Evol. 2007;65:541–54.
Ekblom R, Sæther SA, Jacobsson P, Fiske P, Sahlman T, Grahn M, et al. Spatial pattern of MHC class II variation in the great snipe (Gallinago media). Mol Ecol. 2007;16:1439–51.
Burri R, Hirzel HN, Salamin N, Roulin A, Fumagalli L. Evolutionary patterns of MHC class II B in owls and their implications for the understanding of avian MHC evolution. Mol Biol Evol. 2008;25:1180–91.
Hughes CR, Miles S, Walbroehl JM. Support for the minimal essential MHC hypothesis: a parrot with a single, highly polymorphic MHC class II B gene. Immunogenetics. 2008;60:219–31.
Kikkawa EF, Tsuda TT, Sumiyama D, Naruse TK, Fukuda M, Kurita M, et al. Trans-species polymorphism of the Mhc class II DRB-like gene in banded penguins (genus Spheniscus). Immunogenetics. 2009;61:341–52.
Wang Z, Zhou X, Lin Q, Fang W, Chen X. Characterization, polymorphism and selection of major histocompatibility complex (MHC) DAB genes in vulnerable Chinese egret (Egretta eulophotes). PLOS ONE. 2013;8:e74185.
Shiina T, Hosomichi K, Hanzawa K. Comparative genomics of the poultry major histocompatibility complex. Anim Sci J. 2006;77:151–62.
Anmarkrud JA, Johnsen A, Bachmann L, Lifjeld JT. Ancestral polymorphism in exon 2 of bluethroat (Luscinia svecica) MHC class II B genes. J Evol Biol. 2010;23:1206–17.
Balakrishnan CN, Ekblom R, Völker M, Westerdahl H, Godinez R, Kotkiewicz H, et al. Gene duplication and fragmentation in the zebra finch major histocompatibility complex. BMC Biol. 2010;8:29.
Zagalska-Neubauer M, Babik W, Stuglik M, Gustafsson L, Cichoń M, Radwan J. 454 sequencing reveals extreme complexity of the class II Major Histocompatibility Complex in the collared flycatcher. BMC Evol Biol. 2010;10:395.
Sepil I, Moghadam HK, Huchard E, Sheldon BC. Characterization and 454 pyrosequencing of Major Histocompatibility Complex class I genes in the great tit reveal complexity in a passerine system. BMC Evol Biol. 2012;12:68.
Hess CM, Edwards SV. The evolution of the major histocompatibility complex in birds. Bioscience. 2002;52:423–31.
Babik W, Taberlet P, Ejsmond MJ, Radwan J. New generation sequencers as a tool for genotyping of highly polymorphic multilocus MHC system. Mol Ecol Resour. 2009;9:713–9.
Galan M, Guivier E, Caraux G, Charbonnel N, Cosson J-F. A 454 multiplex sequencing method for rapid and reliable genotyping of highly polymorphic genes in large-scale studies. BMC Genomics. 2010;11:296.
Nishida C, Ishijima J, Ishishita S, Yamada K, Griffin DK, Yamazaki T, et al. Karyotype reorganization with conserved genomic compartmentalization in dot-shaped microchromosomes in the Japanese mountain hawk-eagle (Nisaetus nipalensis orientalis, Accipitridae). Cytogenet Genome Res. 2013;141:284–94.
Burri R, Niculita-Hirzel H, Roulin A, Fumagalli L. Isolation and characterization of major histocompatibility complex (MHC) class II B genes in the Barn owl (Aves: Tyto alba). Immunogenetics. 2008;60:543–50.
Burri R, Promerová M, Goebel J, Fumagalli L. PCR-based isolation of multigene families: lessons from the avian MHC class IIB. Mol Ecol Resour. 2014;14:778–88.
Pavey SA, Sevellec M, Adam W, Normandeau E, Lamaze FC, Gagnaire P-A, et al. Nonparallelism in MHCIIβ diversity accompanies nonparallelism in pathogen infection of lake whitefish (Coregonus clupeaformis) species pairs as revealed by next-generation sequencing. Mol Ecol. 2013;22:3833–49.
Edgar RC. MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics. 2004;5:113.
Babik W. Methods for MHC genotyping in non-model vertebrates. Mol Ecol Resour. 2010;10:16228–33.
Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61:539–42.
Zwickl DJ. Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasets under the maximum likelihood criterion. In: Ph.D. thesis. Austin: The University of Texas; 2006.
Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9:772.
Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28:2731–9.
Brown JH, Jardetzky TS, Gorga JC, Stern LJ, Urban RG, Strominger JL, et al. Three-dimensional structure of the human class II histocompatibility antigen HLA-DR1. Nature. 1993;364:33–9.
Wilson DJ, McVean G. Estimating diversifying selection and functional constraint in the presence of recombination. Genetics. 2006;172:1411–25.
Martin DP, Lemey P, Lott M, Moulton V, Posada D, Lefeuvre P. RDP3: a flexible and fast computer program for analyzing recombination. Bioinformatics. 2010;26:2462–3.
Pond SLK, Posada D, Gravenor MB, Woelk CH, Frost SDW. GARD: a genetic algorithm for recombination detection. Bioinformatics. 2006;22:3096–8.
R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2012 [http://www.R-project.org/]
Goslee SC, Urban DL. The ecodist package for dissimilarity-based analysis of ecological data. J Stat Softw. 2007;22:1–19.
Jaccard P. Nouvelles recherches sur la distribution florale. Bull Soc Vaud Sci Nat. 1908;44:223–70.
Nei M. Analysis of gene diversity in subdivided populations. Proc Natl Acad Sci U S A. 1973;70:3321–3.
Winter DJ. MMOD: an R library for the calculation of population differentiation statistics. Mol Ecol Resour. 2012;12:1158–60.
Excoffier L, Smouse PE, Quattro JM. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics. 1992;131:479–91.
Dray S, Dufour AB, Chessel D. The ade4 package-II: Two-table and K-table methods. R News. 2007;7:47–52.
Klein J, Bontrop RE, Dawkins RL, Erlich HA, Gyllensten UB, Heise ER, et al. Nomenclature for the major histocompatibility complexes of different species: a proposal. Immunogenetics. 1990;31:417–26.
Klein J. Origin of major histocompatibility complex polymorphism: the trans-species hypothesis. Hum Immunol. 1987;19:155–62.
Whitlock MC. G’ST and D do not replace F ST. Mol Ecol. 2011;20:1083–91.
Burri R, Salamin N, Studer RA, Roulin A, Fumagalli L. Adaptive divergence of ancient gene duplicates in the avian MHC class II β. Mol Biol Evol. 2010;27:2360–74.
Dawkins R, Leelayuwat C, Gaudieri S, Tay G, Hui J, Cattley S, et al. Genomics of the major histocompatibility complex: haplotypes, duplication, retroviruses and disease. Immunol Rev. 1999;167:275–304.
Siddle HV, Deakin JE, Coggill P, Whilming LG, Harrow J, Kaufman J, et al. The tammar wallaby major histocompatibility complex shows evidence of past genomic instability. BMC Genomics. 2011;12:421.
Hughes AL, Nei M. Pattern of nucleotide substitution at major histocompatibility complex class I loci reveals overdominant selection. Nature. 1988;335:167–70.
Hughes AL, Yeager M. Natural selection at major histocompatibility complex loci of vertebrates. Annu Rev Genet. 1998;32:415–35.
Alcaide M, Edwards SV, Negro JJ, Serrano D, Tella JL. Extensive polymorphism and geographical variation at a positively selected MHC class II B gene of the lesser kestrel (Falco naumanni). Mol Ecol. 2008;17:2652–65.
Alcaide M, Muñoz J, Martínez-de la Puente J, Soriguer R, Figuerola J. Extraordinary MHC class II B diversity in a non-passerine, wild bird: the Eurasian Coot Fulica atra (Aves: Rallidae). Ecol Evol. 2014;4:688–98.
We thank M. H. Dick and M. T. Kimura for invaluable comments and editing the draft manuscript, R. Burri and one anonymous reviewer for insightful comments and suggestions. We also thank Kushiro Zoo and the Institute for Raptor Biomedicine Japan for supplying samples. Samples were originally obtained through conservation activities for Blakiston’s fish owl by the Ministry of the Environment, Japan. The study was supported in part by a Grant-in-Aid for Scientific Research (No. 23310164) from the Japan Society for the Promotion of Science (JSPS) to RM, a JSPS Fellowship (No. 13 J01622) to KO, and a grant (No. D-1201) from the Environment Research and Technology Development Fund of the Ministry of the Environment, Japan.
The authors declare that they have no competing interests.
TIK designed the study, performed MHC genotyping, analyzed the data, and wrote the manuscript. KO assisted in laboratory work. CN performed the cell culture. TT made field collections during conservation work. KS and SF provided samples. RM gave the guidance through the study. All authors read and approved the final manuscript.
About this article
Cite this article
Kohyama, T.I., Omote, K., Nishida, C. et al. Spatial and temporal variation at major histocompatibility complex class IIB genes in the endangered Blakiston’s fish owl. Zoological Lett 1, 13 (2015). https://doi.org/10.1186/s40851-015-0013-4