Sampling, and RNA and DNA extraction

Samples for RNA extraction were obtained from two individuals each of B. blakistoni and the Eurasian eagle owl (Bubo bubo), and one individual of the Ural owl (Strix uralensis). Blood samples were collected with the anticoagulant heparin and precipitated at 4°C within 12 hours after sampling. Leukocytes were isolated from the fresh blood and cultured at 39°C for 72 hours under 5% CO₂ in RPMI1640 medium supplemented with 20% fetal bovine serum (FBS), 90 µg/ml of phytohemagglutinin, 3 µg/ml of concanavalin A, and 10 µg/ml of lipopolysaccharide as mitogens. RNA extractions were performed with the RNeasy Plus Mini Kit (Qiagen). Extracted RNA samples were stored at -80°C until use.

Samples for genomic DNA extraction were obtained from 10 individuals of B. blakistoni; two individuals each of B. bubo, the great horned owl (Bubo virginianus), and the snowy owl (Bubo scandiacus); and one individual each of the tawny fish owl (Bubo flavipes) and S. uralensis. DNA extractions were performed with the DNeasy Blood & Tissue Kit (Qiagen).

Sequencing of MHC class I and II genes

To determine long sequences from RNA transcripts from MHCI and MHCII genes, we synthesized complimentary DNA (cDNA) libraries from total RNA from B. blakistoni, B. bubo, and S. uralensis by reverse transcription with the 3′-Full RACE Core Set (TaKaRa). Target genes were amplified by PCR using unique forward primers and a common reverse primer (the 3′-site adaptor primer; 3′-Full RACE Core Set), and PrimeSTAR GXL DNA Polymerase (TaKaRa). The forward primer MHCI-ex2F (5′-CGC TAC AAC CAG ASC RRS G-3′) (Alcaide et al., 2008) was used to amplify MHCI alleles, and the forward primers AvesEx1-F1 (5′-ACT GGT GGC ACT GGT GGY GC-3′) and AvesEx1-F2 (5′-GCA CTG GTG GYG CTG GGA GC-3′) (Burri et al., 2014) were used for semi-nested PCR to amplify MHCIIβ alleles. PCR products were subjected to agarose-gel electrophoresis and visualized by ethidium bromide staining under ultraviolet illumination, and candidate bands were excised and purified with the QIAEX II Gel Extraction Kit (Qiagen). Purified PCR products were cloned into the plasmid vector pBluescript II SK (+), and then transformed into Escherichia coli strain DH5α. Plasmids were purified with the QIAprep Spin Miniprep Kit (Qiagen) and used as templates for Sanger sequencing. Sequencing was performed with the BigDye Terminator v3.1 Cycle Sequencing Kit and an ABI 3730 DNA automated sequencer (Thermo Fisher Scientific) using M13 forward and reverse primers. Nucleotide sequences were aligned by using MEGA 5.0 software (Tamura et al., 2011).

We utilized sequence data from Kohyama et al. (2015), who obtained partial sequences (203 bp) of MHCIIβ exon 2 from 174 B. blakistoni individuals sampled on Hokkaido Island, by massively parallel pyrosequencing. In addition, sequences were obtained from one individual of B. flavipes, and two individuals each of B. bubo, B. virginianus, and B. scandiacus, by the same method as in Kohyama et al. (2015). Briefly, pyrosequencing was performed on PCR products obtained from genomic DNA by using the primer set BublIIb2F (5′-GAG TGT CAG YAC CTY RAY RG-3′) and BublIIb2R (5′-CTT TCY TCT SCS TGA YGW AGG-3′). The amplicon library was commercially sequenced on a 1/4 Titanium Pico-Titer Plate with the GS FLX Titanium Sequencing Kit XLR70 (Roche) at Hokkaido System Science Co. (Sapporo, Japan).

We reconstructed phylogenetic relationships among MHCIIβ alleles by Bayesian inference using BEAST 1.6.2 software (Drummond et al., 2006; Drummond and Rambaut, 2007). Phylogenetic trees of MHCIIβ alleles were reconstructed separately for exon 2 and exon 3 sequences. In owl species, MHCIIβ exon 3 alleles show two distinct lineages, DAB1 and DAB2 (Burri et al., 2014). The cloned sequences, whose allelic lineages were defined, were used for phylogenetic analysis, which included MHCIIβ sequences from 14 other owl species obtained from GenBank (Figs. 1A and 1B).

Quantifying copy numbers of MHC genes

To estimate the copy numbers of MHCI and MHCIIβ genes, quantitative real-time PCR (qPCR) was used on genomic DNA, with a double-strand DNA-binding dye as the reporter. Primer sets were designed for MHCI exon 2 (MHCI-Ex2-F3, 5′-CGC AAA TCA CCA AGA GGA AG-3′; MHCI-Ex2-R1, 5′-GCC CGT AGC TCA CGT ATT TC-3′) and MHCIIβ exon 3 (MHCII-Ex3-F5, 5′-CGA GGT GAA GTG GTT CCA GA-3′; MHCIIEx3-R1, 5′-TGG TTT CCA GCA TCA CCA G-3′), based on sequences conserved in B. blakistoni, B. bubo, and S. uralensis. The glyceraldehyde 3-phosphate dehydrogenase (GAPDH) gene was used as a single-copy reference-gene control for the qPCR, based on copy number validation in avian species (Criscuolo et al., 2009). The primer set GAPDH-F2 (5′-CCA TCA CAG CCA CAC AGA AG-3′) and GAPDH-R2 (5′-TTA GCA GCC CCA GTA GAT GC-3′) was designed for GAPDH exon 3, based on sequences from GenBank: Strix aluco (KF201575), Taeniopygia guttata (NM_001198610), Falco cherrug (XM_014287574), Nipponia nippon (XM_009461700), and Gallus gallus (V00407). The three genes were analyzed by qPCR for B. blakistoni, B. flavipes, B. bubo, B. virginianus, B. scandiacus, and S. uralensis; qPCR reactions were performed with the KAPA SYBER Fast qPCR Kit (Kapa Biosystems). Reactions were performed forth times with primer sets MHCI, MHCIIβ, and GAPDH and 10 µl of KAPA SYBER FAST qPCR Master Mix. Real-time amplifications were performed with the 7300 Real-Time PCR System (Applied Biosystems); cycling conditions were 95˚C for 10 min; 40 cycles of 95˚C for 5 sec and 60˚C for 30 sec; and a following dissociation stage. To generate a reference curve for estimating the amplification efficiency, qPCR was run on serial dilutions (4.8 ng, 2.4 ng, 1.2 ng, 0.6 ng, 0.3 ng, and 0.15 ng of template DNA per reaction).

The number of MHC genes was estimated by comparing the threshold cycle value (Ct) with that of the internal single-gene control (GAPDH) measured from the same DNA sample, using the following formula: gene number = 2^-ΔCt where ΔCt = Ct^{target gene} - Ct^control.

Mapping of MHC loci on chromosomes

To determine the locations of B. blakistoni and B. bubo MHC genes on chromosomes, fl uorescence in situ hybridization (FISH) was performed as described by Matsuda and Chapman (1995). Cloned plasmids containing MHCI (1035–1036 bp) or MHCIIβ (990 bp) sequences were used as probes. In addition, human 18 S and 28 S rRNA gene sequences, pHr21Ab (5.8 kbp) and pHr14E3 (7.3 kbp), respectively, provided by the Japanese Cancer Research Resource Bank, Japan, were used as probes, as these rRNA genes are closely linked with one another and with MHC genes on the same microchromosome in some bird species, such as chicken (Delany et al., 2009) and turkey (Chaves et al., 2007). In a preliminary study, we confi rmed that the human 18 S and 28 S rDNA probes could detect the microchromosomes on which the nucleolus organizer region is located in B. blakistoni. Each probe was labeled with Cy3-dUTP by using the Nick Translation Kit (Roche Applied Science) and ethanol-precipitated along with salmon sperm DNA and E. coli tRNA.

Single-color FISH mapping was performed for three types of probes on the same chromosome samples by triplicate hybridization and washing in the order MHCIIβ, MHCI, and 18 S/28 S. Chromosome samples were prepared on slides from cultured fi broblasts of B. blakistoni and B. bubo. After FISH, the slides were washed with 3:1 methanol/acetic acid, rinsed in distilled water, and air-dried. The chromosomal locations of the genes were examined on the same metaphase spreads for each probe by using a Nikon Eclipse E800 fl uorescence microscope.

Sequence variation and phylogeny of the MHC alleles

One MHCI allele was identifi ed from more than 10 clones in both B. blakistoni and B. bubo partial sequences (1008 bp). Some sequence variants were also detected in each species, but we could not rule out PCR and/or sequencing errors as the cause, as each variant was identifi ed from only a single clone. More than 10 MHCIIβ alleles were found in multiple clones of both B. blakistoni and B. bubo (792–1003 bp), which included part of exon 1 and complete exons 2–6. Two distinct allelic lineages for MHCIIβ, DAB1, and DAB2, based on exon 3 sequences (Burri et al., 2014), were detected among the sequences obtained by cloning. The newly identifi ed allele sequences were deposited in the DDBJ (DNA Data Bank of Japan) database under accession numbers LC169500–LC169503.

Fig. 1.

Phylogenetic relationships among MHC class IIβ exon 3 (A) and MHC class IIβ exon 2 (B) alleles from owl species, reconstructed by Bayesian analysis. Black font, DAB1 alleles; gray font, DAB2 alleles. Numbers near nodes are posterior probability values. Accession numbers for allele sequences are in parentheses.

In the Bayesian tree of MHCIIβ alleles based on exon 3 sequences (215 bp) (Fig. 1A), DAB1 and DAB2 alleles from strigid species each form a distinct clade, whereas those from the barn owl (Tyto alba; Tytonidae) group together in yet another clade. By contrast, in the Bayesian tree based on MHCIIβ exon 2 sequences (268 bp) (Fig. 1B), in the strigid clade DAB1 alleles form an unresolved polytomy with DAB2 alleles, i.e., the distinct allelic lineages seen in Fig. 1A are completely intermixed. Some clades consist of alleles from a single strigid species, but most include alleles from two or more strigid species.

Copy numbers of the MHC genes

The Ct values (Ct^MHCI, Ct^MHCIIβ, and Ct^control) obtained by qPCR varied only slightly among replicates from the same samples, and so mean Ct values were used to estimate the copy numbers of the genes. The qPCR results indicate there are two to four MHCI genes: 2-^ΔCt = 3.9 in B. blakistoni, 2.0 in B. fl avipes, 2.6 in B. bubo, 3.7 in B. virginianus, 3.7 in B. scandiacus, and 2.1 in S. uralensis (Fig. 2). The number of MHCIIβ genes was much more variable among these species: 2-^ΔCt = 11.9 in B. blakistoni, 6.0 in B. fl avipes, 4.1 in B. bubo, 5.9 in B. virginianus, 3.6 in B. scandiacus, and 4.1 in S. uralensis (Fig. 2). Bubo blakistoni had the largest number of genes for both MHCI and MHCIIβ, and in each case twice as many as in B. flavipes.

Fig. 2.

Number of gene copies in six owl species, estimated by quantitative real-time PCR. Black and white bars refer to MHC class I and class IIβ genes, respectively. The error bars show standard deviations.

Table 1.

Sample sizes, chromosome number (Rebholz et al., 1993) and numbers of MHC class IIβ alleles per individual detected by massively parallel pyrosequencing. The pyrosequencing data for B. blakistoni are from Kohyama et al. (2015).

Fig. 3.

FISH mapping of MHC class I (red), MHC class II (green) and the 18 S/28 S locus (purple) on chromosomes of B. blakistoni (A) and B. bubo (B). Solid arrows indicate the MHC loci; dashed arrows indicate the 18 S/28 S locus. Scale bars, 10 µm.

Table 1 shows the number of alleles identifi ed from each individual by massively parallel pyrosequencing of MHCIIβ exon 2. The maximum number of alleles per individual was larger in B. blakistoni (16) than in the other species: 12 B. flavipes, seven in B. bubo, four in B. virginianus, and seven in B. scandiacus (Table 1). Two sequences in MHCIIβ exon 2 were shared by Bubo fl avipes and B. blakistoni.

Location of the MHC genes on chromosomes

MHCI, MHCIIβ, and 18 S/28 S were mapped onto chromosomes for B. blakistoni and B. bubo. The FISH mapping (Fig. 3) shows that the MHCI and MHCIIβ loci are closely linked on a single pair of microchromosomes in both species, whereas the 18 S/28 S loci are on a different pair of microchromosomes from the MHC loci in both species (Fig. 3).

Duplication patterns of MHC genes in Bubo species

The MHC region contains the most diverse set of coding genes in vertebrates because variation in the MHC genes is necessary to recognize many types of antigens in immune responses. The chicken (Gallus gallus) has the minimum essential MHC complement (Kaufman et al., 1999), whereas the MHC of the quail (Coturnix japonica) and some passerine birds contains many more genes (Shiina et al., 2004; Westerdahl et al., 2000; Balakrishnan et al., 2010). Most genes in the MHC family are assumed to have evolved by gene duplication, which is mediated mainly by recombination, replication slippage, retrotransposition, or chromosome duplication. The number of retroelements in C. japonica and some passerine birds is high for birds, but lower than in humans, whereas no retroelements are present in G. gallus (Hess et al., 2002). Balakrishnan et al. (2010) reported highly duplicated MHC genes and interspersed among many chromosomes in the zebra finch (Taeniopygia guttata).

In our study, qPCR-based copy number validation revealed high variation in the numbers of MHCI and MHCIIβ genes among Bubo species. The number of MHCIIβ genes was more variable, and that of B. blakistoni (12 ± 2) was markedly larger than in the other Bubo species (fewer than six) (Fig. 2). Massively parallel pyrosequencing detected the largest number of MHCIIβ alleles in B. blakistoni, up to 16 per individual. This result indicates that B. blakistoni has more than eight copies of the MHCIIβ gene. The pyrosequencing results for the other Bubo species are congruent with the number of MHCIIβ genes estimated by qPCR (Fig. 2). These results indicate frequent duplications of MHCIIβ genes have occurred in Bubo species.

Bubo blakistoni had twice as many MHCI and MHCIIβ genes as B. flavipes. According to a phylogenetic study based on mitochondrial DNA (Omote et al., 2013), B. flavipes is closely related to B. blakistoni. The two sequences shared by the two species for the ABS in MHCIIβ exon 2 also support a close relationship between B. blakistoni and B. flavipes. This finding suggests that duplications of the MHC genes occurred in B. blakistoni after speciation; probably the entire complement MHC gene underwent duplication simultaneously.

To reveal the pattern of gene duplication, we investigated the location of the MHC gene loci on chromosomes of B. blakistoni and B. bubo, whose MHCI and MHCIIβ genes were less than half as numerous of those of B. blakistoni. Bubo blakistoni and B. bubo have 82 and 80 chromosomes, respectively, including microchromosomes (Rebholz et al., 1993) (Table 1). Although the MHC gene family is linked with 18 S/28 S on a microchromosome in chicken and turkey (Delany et al., 2009; Chaves et al., 2007), our FISH results showed that the MHCI and MHC IIβ genes were on another microchromosome than 18 S/28 S in B. blakistoni and B. bubo. The FISH results also showed that the MHCI and MHCIIβ loci are located in the same sites on a pair of chromosome, indicating that duplication of the MHC genes was not at the chromosomal level, but was the result of tandem duplications within the gene family region. To clarify the detailed duplication history of MHC genes in B. blakistoni, it will be necessary to survey the whole-genomic structure in the MHC region, using such methods as BAC clone sequencing.

Variation in MHCIIβ alleles in Bubo blakistoni

Among strigid species, the DAB1 and DAB2 alleles of MHCIIβ exon 3 formed completely separate clades (Fig. 1A), indicating that the two allele lineages diverged before radiation of the Strigidae and evolved as paralogs. In contrast, the DAB1 and DAB2 alleles with respect to MHCIIβ exon 2, including the ABS, did not form distinct clades (Fig. 1B). This might be due to gene conversion among DAB1 and DAB2 alleles, resulting in rearrangements of sequences and an increase in allelic variation. The ABS sequences of DAB1 and DAB2 alleles of strigid species formed multi-genus clades (Fig. 1B). The similarities in ABS sequences among species could have been caused by convergent or parallel evolution. Because MHC variation affects the adaptation of populations through having functional importance such as parasite resistance (Sommer, 2005), the trans-species polymorphism may be maintained by heterozygote advantage or frequency-dependent selection.

In general, the high level of variation in MHC genes often observed in wild animals is related to resistance to diseases, but it is difficult to maintain this variation in small populations because of the strong effect of genetic drift and inbreeding. Bubo blakistoni preys on freshwater fishes in the subarctic region and needs large ranges along rivers. Other fish owl species, including B. flavipes, inhabit tropical or subtropical regions in Southeast Asia. The B. blakistoni populations could have experienced severe bottlenecks during glacial periods, which could have reduced genetic variation. In addition, recent population fragmentation and decline owing to human activities has resulted in decreased genetic variation (Omote et al., 2015), and the level of MHCIIβ variation has decreased in the B. blakistoni population on Hokkaido Island during the past decades (Kohyama et al., 2015). However, the number of MHCIIβ alleles per individual remains large among owl species (Table 1), probably because of multiple loci. When loci are duplicated, the number of alleles per individual can markedly increase. Duplications of MHC genes could have prevented a decline in immunity in the small B. blakistoni population and contributed to its successful range expansion into sever subarctic habitats.