Whole Mitochondrial DNA (mtDNA) Sequencing. For sequencing the whole mtDNA genome of B. blakistoni, we utilized blood and tissue samples from Russia (n = 6) and blood samples from Hokkaido Island (n = 10). We used Tawny Fish-Owl (Ketupa flavipes, also known as Bubo flavipes), a closely related species distributed in southeastern Asia, as an outgroup taxon. We cultured fibroblasts from a skin sample of K. flavipes obtained from the Ueno Zoological Gardens to extract the DNA. Blood and tissue samples were preserved in ethanol at −20°C until use. Total DNA was extracted from blood, tissue, or fibroblasts by using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, North Rhine-Westphalia, Germany), and from feather roots by using the QIAamp DNA Micro Kit (Qiagen).

Whole mitochondrial genomes were amplified from total DNA in six overlapping fragments (2.5–7.8 kilo base-pairs [kbp]) by polymerase chain reaction (PCR) using PrimeSTAR GXL DNA Polymerase (Takara, Kusatsu, Shiga, Japan). The PCR cycling conditions were 98°C for 3 min; 30 cycles of 98°C for 10 sec, 60°C for 15 sec, and 68°C for 15 sec to 3 min; and 72°C for 10 min. The PCR products to be used as templates for nucleotide sequencing were purified with the QIAquick PCR Purification Kit (Qiagen). Sequencing was performed with the BigDye Terminator 3.1 Cycle Sequencing Kit and an ABI 3730 DNA automated sequencer (Applied Biosystems: Waltham, MA, U.S.A.). Primers for PCR and sequencing (Table 1) were newly designed, or were used or modified from those in Sorenson et al. (1999) and Omote et al. (2013). Nucleotide sequences were aligned by using MEGA 5.0 software (Tamura et al. 2011). The whole mtDNA sequences of B. blakistoni and K. flavipes were deposited in GenBank: accession nos. LC099100-LC099107.

Table 1. 

Primers for polymerase chain reaction and sequencing used to delineate phylogeny of Blakiston's Fish-Owl (Bubo blakistoni).

Phylogenetic Analysis. To examine the phylogenetic relationships among mtDNA haplotypes and divergence times between clades, we used the sequences of all 13 protein-coding genes (11,373 base pairs [bp]) because extensive sequence data were required for the accurate estimation of branch lengths among the closely related lineages. The 13 protein-coding sequences of three samples from Russia and four samples from Hokkaido Island were used for analysis. In addition, the entire sequence of the mitochondrial genome of K. flavipes was used as an outgroup taxon. Because two duplicated genes (cytochrome b-2 and ND6-2) were each identical in sequences within individuals, we removed one set from subsequent molecular phylogenetic analyses.

The average Kimura-2-parameter (K2P) differences among and within populations were calculated by using MEGA 5.0. Phylogenetic trees were reconstructed with the maximum likelihood method implemented with Treefinder March 2011 software (Jobb 2011), and by Bayesian inference implemented with BEAST 1.6.2 software (Drummond et al. 2006; Drummond and Rambaut 2007). The best-fit model of substitution was determined independently for each of the 13 protein-coding genes by using MEGA 5.0. A Bayesian skyline model or a coalescent constant size model was used as a tree prior. A posterior set of trees was obtained through the Bayesian MCMC analysis, which was run for 10,000,000 iterations. We used the Bayesian analysis to estimate the divergence times between and within the continental and insular populations using a strict clock model. To calculate the divergence times among haplotypes, we used 1.20% (entire mtDNA genome) and 1.81% (mitochondrial cytochrome b [cytb]) per million yr as the evolutionary rate of all mitochondrial protein-coding genes according to those in Accipitridae (Eo and DeWoody 2010, Weir and Schluter 2004).

Haplotyping of Mitochondrial DNA Sequences. To reveal geographic distributions of mtDNA haplotypes, we utilized blood, tissue, feather, and skull samples from locations in Russia, including Primorye (n = 22), Amur (n = 2), and Magadan (n = 1), and feathers obtained from museum specimens from Kunashir Island in the southern Kuril Islands (n = 2) and Sakhalin Island (n = 1). On Hokkaido Island, the population structure may have changed due to habitat fragmentation and reduction caused by human activities during the 1980s (Omote et al. 2015). To reconstruct the structure of the Hokkaido Island population before this disturbance, we utilized sequence data from Omote et al. (2015) derived from museum specimens (n = 19) collected before 1980. Total DNA was extracted from blood or tissue using the DNeasy Blood and Tissue Kit (Qiagen), and from feather roots using the QIAamp DNA Micro Kit (Qiagen). Total DNA was extracted from a partial skull sample following Masuda et al. (2001).

Sequencing from the DNA samples was performed on shorter fragments (140–590 bp), amplified by using the Multiplex PCR Master Mix (Qiagen): the PCR cycling conditions were 95°C for 15 min; 40 cycles of 94°C for 30 sec, 57°C for 3 min, and 72°C for 30 sec; and 72°C for 10 min. The PCR products to be used as templates for nucleotide sequencing were purified with the QIAquick PCR Purification Kit (Qiagen). Sequencing was performed with the BigDye Terminator 1.1 Cycle Sequencing Kit and an ABI 3730 DNA automated sequencer (Applied Biosystems). On each museum and archaeological sample, we performed PCR and sequencing at least three times to eliminate the potential for contamination and sequencing errors, and thus confirmed authenticity of the data. Haplotypes of the samples were assigned based on the nucleotide sequences of mitochondrial control region 2 (CR-2), using the criteria for haplotyping specified in Omote et al. (2015). The relationships among haplotypes were estimated using statistical parsimony network analysis using TCS 1.21 software (Clement et al. 2000). Calculations of haplotype and nucleotide diversities were performed with ARLEQUIN 3.5.1.2 software (Excoffier and Lischer 2010). Migration rate among populations and growth rate for each population were estimated using LAMARC 2.1.10 software (Kuhner and Smith 2007).

Organization of Whole Mitochondrial Genes of B. blakistoni. The size of the whole mitochondrial genome of B. blakistoni was at least 21 kbp. The genome contained a tandemly duplicated gene unit from the cytb to control region (CR; Fig. 2); due to sequencing difficulties and resulting ambiguity in the large tandemly duplicated repeat region in CR-2, we were unable to determine the exact size of the genome. In the duplicated repeat region, only cytb gene sequences were partial, and the sequences of tRNAThr-2, tRNAPro-2, tRNAGlu-2, and ND6-2 were identical with those of tRNAThr-1, tRNAPro-1, tRNAGlu-1, and ND6-1. Both CR-1 and CR-2 included tandem repeats consisting of four to six copies of 49-bp units (51-bp units in K. flavipes), whereas only CR-2 had a much longer tandem repeat region consisting of more than 20 times of 78-bp or 79-bp units. There were variations in the number and/or order of tandem repeat units as reported by Omote et al. (2013). Sequence regions other than the tandem repeat regions were similar between CR-1 and CR-2, and a few single nucleotide polymorphisms (SNPs) within individuals from the Russian Far East and from Hokkaido Island were found in the first domain of CR-2 (Omote et al. 2013).

Figure 2. 

Schematic diagram of the entire mitochondrial genome of Blakiston's Fish-Owl (Bubo blakistoni). The positions of polymerase chain reaction and sequencing primers (listed in Table 1) are indicated by arrowheads.

Phylogenetic Relationships and Divergence Times between Clades. The phylogenetic trees from the maximum likelihood and Bayesian analyses were identical in topology (Fig. 3). Both analyses strongly supported high genetic differentiations between continental and insular clades. The results showed that individuals from Hokkaido, southern Kuril, and Sakhalin islands were assigned to the same clade. The average K2P difference based on entire mtDNA sequences within continental and insular populations was 0.052% and 0.017%, respectively, and that between the populations was 1.31%. The Bayesian analysis based on the mutation rate of entire mtDNA genome indicated that continental and insular clades diverged 669 kilo yr before present (kyrBP; 95% highest posterior density interval [HPDI], 560–770 kyrBP). The timing of the most recent common ancestor for the continental clade was estimated to be 9.3 kyrBP (95% HPDI: 1.6–20 kyrBP), and that for the insular clade was estimated to be 10 kyrBP (95% HPDI: 2.3–23 kyrBP; Fig. 3). The divergence time between B. blakistoni and K. flavipes was estimated to be 3461 kyrBP (95% HPDI: 3200–3700 kyrBP). In addition, the analysis based on the mutation rate of mitochondrial cytb indicated that divergence time was 506 kyrBP (95% HPDI: 424–583 kyrBP) between the clades, 7.0 kyrBP (95% HPDI: 1.2–15 kyrBP) and 7.6 kyrBP (95% HPDI: 1.7–17 kyrBP) within the clades, and 2620 kyrBP (95% HPDI: 2400–2800 kyrBP) between B. blakistoni and K. flavipes, respectively.

Figure 3. 

The molecular phylogeny based on the whole mitochondrial genome sequences of Bubo blakistoni and K. flavipes as an outgroup. Numbers above the tree branches indicate the Bayesian posterior probabilities / ML bootstrap values. Numbers below branches show the estimated divergence times, and bars show mean 95% highest posterior density intervals; kyrBP = kilo yr before present.

Geographical Distribution of CR-2 Haplotypes. Seven CR-2 haplotypes based on nucleotide substitutions were found: three haplotypes from the Russian Far East, and four haplotypes from Hokkaido Island (Fig. 4). In northeastern continental Asia, we determined new haplotypes for 13 samples from Primorye, two from Amur, and one from Magadan (Fig. 4A, B). The most common continental haplotype occurred in 12 of the 13 samples from Primorye and in the sample from Magadan, the northern boundary of the range of B. blakistoni. The other two haplotypes were minor and found in inland areas. On Hokkaido Island, three of the four haplotypes were widespread (Fig. 4C). The most frequent haplotype on Hokkaido Island was also found on Sakhalin and Kunashir islands. The parsimony network of the seven haplotypes clearly showed two clusters, which consist of continental and insular haplotypes, respectively (Fig. 4D). Haplotype and nucleotide diversities were 0.34 and 0.001 in the continental population, 0.69 and 0.004 in the insular population, respectively. The likelihood estimation showed positive growth rate for the continental and insular populations (3338 and 771, respectively) and low migration rate between the populations.

Figure 4. 

Estimated geographical distribution of mtDNA haplotypes of Bubo blakistoni in: (A) northeast Asia; (B) Primorye; and (C) Hokkaido and southern Kuril islands. (D) Parsimony network of the haplotype based on mitochondrial CR-2 sequences. Each symbol refers to one haplotype and one bar means one substitution.

Phylogeography of B. blakistoni Populations. Blakiston's Fish-Owl is classified into two subspecies, B. b. blakistoni on islands and B. b. doerriesi on the continental areas (Slaght and Surmach 2008); however, the genetic relationships between the two subspecies have not been well studied to date. Our results showed that haplotypes based on mtDNA CR-2 were not shared between insular and continental populations (Fig. 4D). The phylogenetic analysis using the whole mitochondrial genome showed two clearly separated clades: one clade consisted of individuals from Hokkaido, southern Kuril, and Sakhalin islands; and the other of individuals from the Eurasian continent (Fig. 3). The results suggest that the geographically isolated populations have been separated for a long time.

Generally, 98% of sister-species pairs of vertebrates were observed to have K2P mtDNA divergences more than 2%; however, avian taxa had less genetic divergence than same-rank taxa in other vertebrate groups (Johns and Avise 1998). According to analyses on birds using mitochondrial cytochrome c oxidase I sequences, the average K2P difference between sister-species pairs varied from 0.78% to 11.77% (Tavares and Baker 2008). Because differences within species were much lower than those among sister-species, the suggested threshold of 10 times the mean intraspecific variation to screen for splits referred to as putative species has also been criticized (Hebert et al. 2004). For B. blakistoni, the migration rate between the continental and insular populations was estimated to be very low. The average K2P mitochondrial difference between the populations was 1.310%, and that within the populations was <0.052%; therefore, the ratio of between-populations to within-populations difference was >25. The results suggest that the level of difference between the populations is enough to recognize them as separate species. Because the two subspecies differ in plumage markings (e.g., a white spot on the head in B. b. doerriesi) and call patterns (Slaght and Surmach 2008), the geographical populations may be behaviorally isolated as well. Blakiston's Fish-Owl is listed as endangered on the IUCN Red List, and the size of insular population was much smaller than that of continental population. The continental population size may be 800 pairs, based on extrapolation, whereas the population was only an estimated 140 individuals on Hokkaido Island and 70–85 individuals on Kunashir Island, and there was no confirmation of breeding on Sakhalin Island (Slaght and Surmach 2008, Omote et al. 2015). As the present study suggests, both populations are unique and irreplaceable. It means that reintroduction of individuals across the populations is inappropriate, and indicates the greater necessity for conservation of the tiny and divided population on the islands.

Bayesian analysis shows that the insular and continental populations diverged at least 500 kyrBP. Considering the data from Antarctic ice cores (Lüthi et al. 2008), the two populations have experienced five or six global glacial-interglacial cycles since their divergence. Northern Sakhalin and the adjacent continental areas have non-forested areas that may act as a barrier to dispersal of this owl. The divergence times within the continental and insular populations were estimated to be approximately 9.3 and 10 kyrBP, respectively (Fig. 3), which were later than the end of LGM in 19–20 kyrBP (Clark et al. 2009). Populations of this owl would have decreased in size in the LGM due to limited forested habitats and drastic changes in the environment. The low differentiation in the mitochondrial genome within the populations indicates that the mtDNA haplotypes diverged from single mitochondrial lineages that survived the LGM. In continental northeastern Asia, forested areas in the LGM were much smaller and more fragmented than those at present (Fig. 1; Harrison et al. 2001), although Primorye and southern areas were not glaciated during the Pleistocene (Krestov 2003). Blakiston's Fish-Owls need forests and non-frozen rivers or lakes during winter. This means that most of the continental areas where the owls currently occur were unsuitable for them in the LGM. Refugia for the continental populations were likely located in northeastern China or Korean Peninsula, and owls likely dispersed to more northern areas after the LGM. The lower haplotype and nucleotide diversities, higher growth rate, and the distribution pattern of the mtDNA haplotypes in the continental population suggest severe bottleneck and rapid dispersion through the last glacial period. In contrast, the insular population retains a higher haplotype variation. Although environments on northeastern Hokkaido, Sakhalin, and Kuril islands could have been unsuitable for the fish-owl in the LGM, Southern Hokkaido Island was covered with forests (Harrison et al. 2001). The insular population could have survived on southern Hokkaido or Honshu Island, and dispersed to northeastern Hokkaido, Sakhalin, and the Kuril islands after the LGM. Since then, more variations might have been accumulated in the mtDNA and been maintained in the Hokkaido Island population. There are no paleontological records of Blakiston's Fish-Owl from the Japanese islands, and future study on its historical distribution is required for clarification of its dispersion history.

Because all of the three fish-owl species phylogenetically closest to B. blakistoni occur in tropical forests in southern Asia (Omote et al. 2013), B. blakistoni is assumed to have originated in southern Asia. The divergence time between the B. blakistoni and K. flavipes was estimated to be 3500 kyrBP, or during the late Pliocene. During the following Pleistocene, the climate and vegetation were repetitively changed via glacial-interglacial cycles (Lüthi et al. 2008). Glacial-interglacial cycles are thought to have stimulated speciation due to geographical isolation of populations and high selective pressure for adaptation (Weir and Schluter 2004). Ancestors of B. blakistoni could have adapted to colder environments during glacial periods and dispersed to northeastern Asia during interglacial periods.

Duplication and Variation in the Mitochondrial Genomes. In the present study, we found duplication of several genes in mtDNA sequences of B. blakistoni and K. flavipes. This is the first report of mtDNA rearrangements in the order Strigiformes. Variations of gene rearrangements in the mitochondrial genome have been reported in several avian taxa (Gibb et al. 2007). The pattern of the gene duplication in B. blakistoniand K. flavipes was similar to that in albatrosses (Abbott et al. 2005) in that there was an additional incomplete cytb gene and duplicates of ND6, tRNAThr, tRNAPro, tRNAGlu, and CR (Fig. 2). The rearrangements in mitochondrial genomes have occurred several times in various lineages of birds. For example, Schirtzinger et al. (2012) identified six independent origins of mitochondrial CR duplications within the order Psittaciformes. Because the same regions in the mitochondrial genome between both B. blakistoni and K. flavipes are duplicated, this gene rearrangement could have occurred before their speciation. In the duplicated region, the differentiation between paralogous sequences within each species was smaller than that between orthologous sequences between the species, indicating the result of gene conversions. A similar pattern has been reported in some other avian lineages, such as Amazona sp. parrots (Eberhard et al. 2001), Thalassarchesp. albatrosses (Abbott et al. 2005), and Black-faced Spoonbill (Platalea minor; Cho et al. 2009). In these avian lineages, the paralogous sequences, which duplicated before speciation, kept a higher similarity in each species rather than orthologous sequences, indicating that there would be some system synchronizing sequences in an individual.