Recent fragmentation of the endangered Blakiston’s fish owl (Bubo blakistoni) population on Hokkaido Island, Northern Japan, Revealed by Mitochondrial DNA and Microsatellite Analyses
© Omote et al.; licensee BioMed Central. 2015
Received: 10 October 2014
Accepted: 13 January 2015
Published: 29 April 2015
Blakiston’s fish owl (Bubo blakistoni) was previously widespread on Hokkaido Island, Japan, but is now distributed only in limited forest areas. The population size on Hokkaido decreased during the 20th century due to reduction and fragmentation of the owl’s habitat. To elucidate temporal and spatial changes in population structure and genetic diversity, we analyzed 439 individuals collected over the last 100 years.
We detected a population bottleneck and fragmentation event indicated by mitochondrial DNA (mtDNA) haplotype and microsatellite analyses. The lowest value for effective population size, which was estimated by moment and temporal methods from microsatellite data, occurred in the 1980s. Five haplotypes were found in the mtDNA control region; most haplotypes were previously widespread across Hokkaido, but have become fixed in separate areas after the bottleneck period. Genetic differentiation among local populations, as indicated by both mtDNA and microsatellite data, likely arose through population fragmentation.
The owl population may have been divided into limited areas due to loss of habitats via human activities, and have lost genetic variability within the local populations through inbreeding. Our mtDNA and microsatellite data show that genetic diversity decreased in local populations, indicating the importance of individuals moving between areas for conservation of this species on Hokkaido.
Population genetic theory predicts that bottlenecks will lead to decreases in genetic variability, individual fitness, and the decreased capacity of populations to adapt to environmental change . Fragmented populations may also face reduced gene flow, potentially leading to lower genetic diversity, inbreeding depression, fixation of deleterious mutations, and a higher risk of extinction . A microsatellite analysis of the Hokkaido population of Blakiston’s fish owl (120 individuals) in the intermittent intervals (1986–1993, 1997–1999, and 2009–2010) suggested diminished genetic diversity and isolation of local populations . The cause of the low level of genetic variation and gene flow was assumed to be a recent bottleneck; however, no genetic data has demonstrated the changes in population size and structure before and after the bottleneck event assumed from field observations. In the present study, to clarify the temporal changes of population structure and genetic diversity in the Hokkaido population over longer period, we used microsatellite and mitochondrial DNA (mtDNA) analyses on more current, museum, and archaeological samples collected for the last century.
Materials and Methods
Sampling and DNA extraction
Numbers of Blakiston’s fish owls examined, by sampling locations and periods
Analysis of mtDNA sequences
PCR primers newly designed in the present study for amplification and/or sequencing for the Blakiston’s fish owl mtDNA control region
Nucleotide sequence (5′–3′)
Omote et al. 
Omote et al. 
Nucleotide sequences were aligned by using MEGA 5.0 software . Haplotype diversity (h) and nucleotide diversity (π) were calculated using ARLEQUIN 126.96.36.199 software . Genetic differences (Fst) among local populations (Ishikari, Daisetu, Hidaka, Akan, Konsen and Shiretoko) on Hokkaido Island (Figure 1) were calculated by analysis of molecular variance (AMOVA) implemented, and tested with 10,000 permutations by using ARLEQUIN.
Analysis of microsatellite genotypes
The 287 samples of blood, tissues, cultured fibroblasts, and feathers were genotyped for seven autosomal microsatellite loci. In addition, data from 120 samples (62 samples collected from 1986 to 1993; 37 from 1997 to 1999; 21 from 2009 to 2010) from Omote et al.  were also included, and in total, 407 individual were analyzed. We cited microsatellite markers designed for other owl species: Oe058, Oe128, and Oe129 [15,16]; 13D8 and 4E10.2 ; and FEPO5 and FEPO43 . Forward primers were fluorescently labeled with 6-FAM, NED, PET, or VIC. PCR conditions were as previously reported by Omote et al. . The molecular size of PCR products was determined on an ABI 3730 DNA automated sequencer using the GS600 LIZ size standard and GENEMAPPER 4.0 software (Applied Biosystems).
Departures from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) were tested for each microsatellite locus by using ARLEQUIN. The observed (Ho) and expected (He) heterozygosities were calculated using ARLEQUIN. Allelic richness (Ar) and inbreeding coefficient (Fis) were calculated with FSTAT 2.9.3 software . Changes among sampling periods of Ho, He, Ar and Fis values were tested by pairwise non-parametric Wilcoxon signed-rank test by using STATISTICA 10 software (StatSoft). Because the sample size in older periods were small, Fst values among local populations were calculated for each locus and combined by jackknifing, providing confidence intervals via bootstrapping, using FSTAT, and tested with 10,000 permutations by using ARLEQUIN. Effective population size (Ne) was estimated by moment method, the bias-corrected version of the method based on LD , and temporal method using moment biased F-statistics  by using NeEstimater 2.01 software . In the moment method, Ne was estimated on each sampling period with monogamous mating model. In the temporal method, Ne was estimated between continuous two sampling periods, assuming that generation time was five to ten years and that census size was 100 to 150. Populations that have undergone a bottleneck often exhibit a reduction in allele number and Ho, with allele number decreasing faster than Ho . Thus, after a bottleneck, Ho is larger than He estimated from allele frequencies assuming mutation-drift equilibrium . Based on this effect, known as heterozygosity excess, whether a population bottleneck had occurred was tested with BOTTLENECK 188.8.131.52 software [24,25]. To identify heterozygosity excesses, the Wilcoxon signed-rank test was used in default settings under three mutation models: infinite alleles (IAM), two-phased (TPM), and stepwise mutation (SMM).
To determine the number of genetic clusters that best fit the microsatellite data, a Bayesian assignment analysis was performed with STRUCTURE 2.3 software [26,27]. The STRUCTURE program groups multilocus data into K clusters without considering population origin. The STRUCTURE analysis was run for 100,000 iterations after a burn-in of 50,000 iterations, using the admixture model, and posterior probabilities were estimated for K = 1 through K = 7. For each value of K, 10 independent runs were conducted to quantify the amount of variation in the likelihood value for K. To determine the probable number of clusters, ΔK was calculated . Individual sets of 10 replicate STRUCTURE runs were aligned using CRUMPP software . Changes in Ho for the next 50 years were simulated based on the latest genotype data in the 2005–2012 interval by using BOTTLESIM 2.6.1 software . We assumed a single population without barrier against random mating (Ne = 40) or completely isolated five local populations in equal sizes (Ne = 8). The Ne values were based on estimation by moment and temporal methods in the latest interval. For estimates, 1000 iterations were performed with the following constant parameters: mean life span = 16 years; age at maturity = 2 years; completely overlapping generations; each female mates with a single male each year; and sex ratio = 1:1.
Variation and distribution of mtDNA haplotypes
Polymorphic sites of the mtDNA control region haplotypes identified from the Blakiston’s fish owls
Genetic differentiations and diversities caluculated by mtDNA haplotype and microsatellite data for the periods indicated
Diversity and population structure indicated by microsatellite analysis
Temporal changes in the population structure
A previous microsatellite analysis of Blakiston’s fish owls on Hokkaido in the intermittent intervals  suggested that the local populations had been genetically isolated in each area. But it was not directly shown whether the local populations have been originally isolated or recently divided. To address this question, the present study analyzed a larger number of samples (442 individuals), including samples collected before the population decline, by using mtDNA haplotype and microsatellite markers.
Before 1964, most mtDNA haplotypes were distributed over large areas of Hokkaido (Figure 1a), but subsequently became restricted in distribution (Figure 1b,c and d). AMOVA for both the mtDNA haplotype and microsatellite data detected no genetic differentiation among populations in different areas before 1980 (Table 4). These results indicate that gene flow occurred across large areas of Hokkaido before 1980. Although the detailed population structure was unclear due to the smaller amount of microsatellite data before the 1980s, there may have been some genetic variation and clines between areas in historical periods. Among-population Fst values from both the mtDNA haplotype and microsatellite data markedly increased after the 1980s (Table 4), suggesting major changes in the haplotype and allele frequencies. The separation and isolation of habitats could have rapidly grown the genetic differentiation among local populations. Departures from HWE were detected for the overall Hokkaido population in the 1981–1996 and 1997–2012 intervals; this phenomenon is known as the Wahlund effect and suggests the separation and isolation of the local populations. Because marker pairs detecting LD were not common, biased mating rather than links of loci could have caused it. Heterozygosity excess in some local populations after 1981 also indicated a recent bottleneck in local areas. The results suggest that the fragmentation of the owl population and subsequent reduction in population size have occurred recently on Hokkaido, probably around 1980 as judged from the Fst values (Table 4). Field researches observed that the population size on Hokkaido Island was minimized in the 1970s to 1980s [2,6]. The change of Ne estimated from microsatellite data also supported a bottleneck event around 1980s (Figure 2). Therefore, the population fragmentation occurred at the nearly the same time as the recent bottleneck. Human activities such as cultivation, housing, and building roads and dams damaged forest and river environments, and may have caused population fragmentation and bottleneck. The Blakiston’s fish owl population would have been decreased and fragmented due to loss and division of habitats via vigorous human activities around 1980s on Hokkaido.
The STRUCTURE analysis showed that the four genetic clusters mainly correspond to the local populations. Although individuals from the Hidaka and Akan areas belonged to the same cluster, the pairwise Fst value between the two local populations was significantly higher. The Hidaka population would have been isolated rather than interacted with the Akan, and the small sample size in the Hidaka probably caused the STRUCTURE result. The 1981–1996 interval was unique in that the Daisetu, Akan, and Konsen populations were each genetically uniform. In the 1997–2012 interval, some individuals were assigned to clusters that had not previously been detected in their local populations, especially in Daisetsu (Figure 3), suggesting gene flow among local populations, and in fact a few owls were observed in the field to have moved from Akan and Shiretoko to Daisetsu and to have bred there in 1990 and 2002 (T. Takenaka, personal communication). In addition, the mtDNA haplotype diversity and microsatellite diversity have increased in Daisetsu from 1981–1996 to 1997–2012 intervals, and significant departure from HWE in Daisetsu might have been caused by the migration. Due to the small population size and long life span of Blakiston’s fish owl, these movements of a few individuals among local populations may have strongly affected the population structure and genetic diversity.
Temporal changes in genetic diversity
Populations that undergo bottlenecks often lose genetic diversity . Based on research on many threatened birds, Heber et al.  concluded that population bottlenecks increase hatching failure due to inbreeding and the consequent loss of genetic diversity, particularly when a population experiences a severe bottleneck to fewer than 100–150 individuals. Changes in genetic diversity have been reported in several avian species, for which samples were analyzed before and after a bottleneck. Microsatellite analysis of the endangered black-capped vireo (Vireo atricapilla), which apparently experienced a bottleneck early in the 20th century, showed lower genetic diversity and increased differentiation in recent (2005–2008) samples than in historical (1899–1915) samples . The great prairie-chicken (Tympanuchus cupido) in Wisconsin, which decreased in number by 91% via loss and fragmentation of the habitat, lost genetic diversity for both mtDNA and microsatellite markers over a period of 50 years . Other studies, however, have documented no changes of genetic diversity in populations that experienced bottleneck events. For example, white-tail eagle (Haliaeetus albicilla) populations in Europe experienced dramatic declines during the 20th century, but have retained high levels of genetic diversity, probably because of long generation time (average life span is around 17 years) . The Spanish imperial eagle (Aquila adalberti), the population of which fragmented into small patches, exhibits reduced mtDNA h and π, but a microsatellite analysis showed undiminished level of nuclear diversity .
Blakiston’s fish owl also has long generation time, and the present study indicates that the Hokkaido population experienced bottleneck and fragmentation event. Although species with long generation times show resistance to loss of nuclear diversity, as reported by Hailer et al.  and Martinez et al. , we observed the loss of a few alleles after the recent bottleneck and detected significant decrease in Ar and Ho (Table 4). This indicates severe bottleneck, and Ne value of the Hokkaido population around 1980 (Figure 2) was actually lower than Ne values reported by Athrey et al.  and Johnson et al.  in populations during bottlenecks. Generally, bottleneck events reduce heterozygosity via loss of alleles, and make He smaller than Ho . Heterozygosity excesses were observed in some local populations, indicating decline of the population size in each local area. However, Ho was decreased more rapidly than He in the overall Hokkaido population. Population fragmentation could explain the results of the present study, as it increases biased mating and cuts population into tinier local patches. The Fis values increasing after 1981 (Table 4) indicate a higher level of inbreeding. The mean of haplotype diversity on local populations was significantly decreased much faster than haplotype diversity on the overall Hokkaido population (Table 4). Our results also support recent fragmentation of the Blakiston’s fish owl population. Both studies on white-tail eagle and Spanish imperial eagle reported undiminished levels of nuclear diversity and declines of mtDNA diversity [35,36]. Because Ne of the mitochondrial genome is 4× lower, mtDNA haplotype diversity may be more sensitive for demographic events than nuclear diversity. The present study showed the level of genetic diversity similar to the two eagle species, but Ho and Fis were significantly changed. Therefore, it indicates that the Blakiston’s fish owl population on Hokkaido would have undergone more serious fragmentation and inbreeding.
Inbreeding often significantly affects birth weight, survival, reproduction, and resistance to factors such as disease, predation and environmental stress . For conservation of species, we need to consider genetic factors as well as ecological factors. Our study demonstrates that the genetic diversity in local populations decreased for the last 30 years due to population fragmentation, and that even a few movements among local populations may counter this decline in genetic diversity in Blakiston’s fish owls. Simulations of changes in Ho for the next 50 years indicated a more rapid decrease in genetic diversity when the local populations remained isolated. Therefore, conserving appropriate habitats for reproduction and promoting movements among local populations will be critical for the survival of populations and the species. Because this species has the ability to move among areas on Hokkaido, as shown in the present study, their movement should be supported and promoted through conservation projects such as corridors connecting isolated habitats.
Hokkaido population of endangered Blakiston’s fish owls has been fragmented into local areas. In the present study, we clarified that most mtDNA haplotypes of old museum samples were distributed over large areas, and that the local populations have been isolated since 1980s. Microsatellite analyses and field observations showed the smallest population size on Hokkaido in around 1980s when forest and river environments were severely damaged by human activities. These indicate that the Blakiston’s fish owl population has recently decreased and fragmented due to loss and division of habitats on Hokkaido. The genetic diversity has been decreased for the last 30 years due to genetic drift and/or inbreeding in the tiny local populations. Our study demonstrates that movement of a small number of individuals between local areas could restore genetic diversity, and that promoting movements among local populations will be effective for conservation of the species. For long-term conservation management of Blakiston’s fish owls on Hokkaido, it is necessary to monitor changes in genetic diversity and study the movement of individual owls in the field.
We thank Sumio Yamamoto (Nemuro City), Yuko Hayashi (Sapporo University) and the Kushiro Nature Conservation Office of the Ministry of the Environment for providing modern samples and invaluable suggestions. The samples were originally obtained through activities to conserve Blakiston’s fish owls by the Ministry of the Environment, Japan. Feather samples from stuffed specimens were provided by the Yamashina Institute for Ornithology (Takeshi Yamasaki, Yasuko Iwami, and Takema Saitou), the Botanical Garden of Hokkaido University (Fumihito Takaya), Tomakomai Experimental Forest Hokkaido University (Hiroshi Sugiyama), The University of Tokyo Hokkaido Forest, the National Museum of Nature and Science Tokyo (Isao Nishiumi), Abiko City Museum of Birds (Ken-ichi Tokita), Betsukai Town Local Museum (Kazuto Ishiwata), Erimo Town Local Museum (Toshiyasu Nakaoka), Kamishihoro Town Higashitaisetsu Museum, Kushiro City Museum, Mombetsu City Museum, Nakashibetsu Town Local Museum, Nemuro City Museum of History and Nature, Rausu Town Local Museum, Shari Town Shiretoko Museum, Sounkyo Visitor Center (Toru Katayama), Tomakomai City Museum (Aiko Kodama) and Urakawa Town Local Museum. We are also grateful to Hiroshi Utagawa and Toshiro Kumaki (University of Tokyo) for providing archaeological bone samples. We deeply thank Prof. Matthew H. Dick (Hokkaido University) for invaluable comments and editing the manuscript.
This work was supported by a Grant-in-Aid for Scientific Research (No. 23310164) from the Japan Society for the Promotion of Science (JSPS); a JSPS fellowship (No. 13 J01622) to K. Omote; and the Environment Research and Technology Development Fund of the Ministry of the Environment, Japan (No. D-1201).
- Slaght JC, Surmach SG. Biology and conservation of Blakiston’s fish-owls (Ketupa blakistoni) in Russia: a review of the primary literature and an assessment of the secondary literature. J Raptor Res. 2008;42:29–37.View ArticleGoogle Scholar
- Takenaka T. Habitat environments, and reasons for reduction of the endangered Blakiston’s fish owl in Hokkaido, Japan. Doctoral Dissertation. Sapporo: Hokkaido University; 1998.Google Scholar
- Yamamoto S. Simahukurou. Sapporo: the Hokkaido Shimbun Press (in Japanese); 1999.Google Scholar
- Hayashi Y. Home range, habitat use and natal dispersal of Blakiston’s Fish-owls. J Raptor Res. 1997;31:283–5.Google Scholar
- 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.View ArticleGoogle Scholar
- Brazil M, Yamamoto S. The status and distribution of owls in Japan. In: Meyburg BU, Chancellor RD, editors. Raptors in the modern world. Berlin, London and Paris: WWGBP; 1989. p. 389–402.Google Scholar
- Keller LF, Waller DM. Inbreeding effects in wild populations. Trends Ecol Evol. 2002;17:230–41.View ArticleGoogle Scholar
- Couvet D. Deleterious effect of restricted gene flow in fragmented populations. Conserv Biol. 2002;16:369–76.View ArticleGoogle Scholar
- 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. Zool Sci. 2012;29:299–304.View ArticlePubMedGoogle Scholar
- Masuda R, Amano T, Ono H. Ancient DNA analysis of brown bear (Ursus arctos) remains from the archaeological site of Rebun Island, Hokkaido, Japan. Zool Sci. 2001;18:741–51.View ArticleGoogle Scholar
- Sato T. Faunal remains excavated from the Nijibetsu Shuwan Kumaokuriba site: a zoo-archaeological study of a Nusa site at an Ainu kotan. Bull Natl Mus Jpn Hist. 2003;107:119–65 (in Japanese with an English abstract).Google Scholar
- Omote K, Nishida C, Dick HM, Masuda R. Limited phylogenetic distribution of a long tandem-repeat cluster in the mitochondrial control region in Bubo (Aves, Strigidae) and cluster variation in Blakiston’s fish owl (Bubo blakistoni). Mol Phylogenet Evol. 2013;66:889–97.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMed CentralPubMedGoogle Scholar
- Excoffier L, Laval G, Schneider S. Arlequin ver. 3.0: an integrated software package for population genetics data analysis. J Evol Bioinform Online. 2005;1:47–50.Google Scholar
- Hsu Y-C, Severinghaus LL, Lin Y-S, Li S-H. Isolation and characterization of microsatellite DNA markers from the Lanyu scops owl (Otus elegans botelensis). Mol Ecol Notes. 2003;3:595–7.View ArticleGoogle Scholar
- Hsu Y-C, Li S-H, Lin Y-S, Severinghaus LL. Microsatellite loci from Lanyu scops owl (Otus elegans botelensis) and their cross-species application in four species of Strigidae. Conserv Genet. 2006;7:161–5.View ArticleGoogle Scholar
- Thode AB, Maltbie M, Hansen LA, Green LD, Longmire JL. Microsatellite markers for the Mexican spotted owl (Strix occidentalis lucida). Mol Ecol Notes. 2002;2:446–8.View ArticleGoogle Scholar
- Proudfoot G, Honeycutt R, Slack RD. Development and characterization of microsatellite DNA primers for ferruginous pygmy-owls (Glaucidium brasilianum). Mol Ecol Notes. 2005;5:90–2.View ArticleGoogle Scholar
- Goudet J. FSTAT, a program to estimate and test gene diversities and fixation indices version 2.9.3. 2001 [http://www2.unil.ch/popgen/softwares/fstat.htm]Google Scholar
- Waples RS, Do C. Linkage disequilibrium estimates of contemporary Ne using highly variable genetic markers: a largely untapped resource for applied conservation and evolution. Evol Applications. 2010;3:244–62.View ArticleGoogle Scholar
- Jorde PE, Ryman N. Unbiased estimator for genetic drift and effective population size. Genetics. 2007;177:927–35.View ArticlePubMed CentralPubMedGoogle Scholar
- Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR. NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol Ecol Res. 2014;14:209–14.View ArticleGoogle Scholar
- Luikart G, Cornuet JM. Estimating the effective number of breeders from heterozygote excess in progeny. Genetics. 1999;151:1211–6.PubMed CentralPubMedGoogle Scholar
- Cornuet JM, Luikart G. Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics. 1996;144:2001–14.PubMed CentralPubMedGoogle Scholar
- Piry S, Luikart G, Cornuet JM. BOTTLENECK: a computer program for detecting recent reductions in the effective population size using allele frequency data. J Hered. 1999;90:502–3.View ArticleGoogle Scholar
- Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.PubMed CentralPubMedGoogle Scholar
- Falush D, Stephens M, Pritchard JK. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics. 2003;164:1567–87.PubMed CentralPubMedGoogle Scholar
- Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.View ArticlePubMedGoogle Scholar
- Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23:1801–6.View ArticlePubMedGoogle Scholar
- Kuo CH, Janzen FJ. BOTTLESIM: a bottleneck simulation program for long-lived species with overlapping generations. Mol Ecol Notes. 2003;3:669–73.View ArticleGoogle Scholar
- Bellinger MR, Johnson JA, Toepfer J, Dunn P. Loss of genetic variation in greater prairie chickens following a population bottleneck in Wisconsin, U.S.A. Conserv Biol. 2003;17:717–24.View ArticleGoogle Scholar
- Heber S, Briskie JV. Population bottlenecks and increased hatching failure in endangered birds. Conserv Biol. 2010;24:1674–8.View ArticlePubMedGoogle Scholar
- Athrey G, Barr KR, Lance RF, Leberg PL. Birds in space and time: genetic changes accompanying anthropogenic habitat fragmentation in the endangered black-capped vireo (Vireo atricapilla). Evol Appl. 2012;5:540–52.View ArticlePubMed CentralPubMedGoogle Scholar
- Johnson JA, Bellinger MR, Toepfer JE, Dunn P. Temporal changes in allele frequencies and low effective population size in greater prairie-chickens. Mol Ecol. 2004;13:2617–30.View ArticlePubMedGoogle Scholar
- Hailer F, Helander B, Folkestad AO, Ganusevich SA, Garstad S, Hauff P, et al. Bottleneck but long-lived: high genetic diversity retained in white-tailed eagles upon recovery from population decline. Biol Letters. 2006;2:316–9.View ArticleGoogle Scholar
- Martinez BC, Godoy JA, Negro JJ. Population fragmentation leads to spatial and temporal genetic structure in the endangered Spanish imperial eagle. Mol Ecol. 2007;16:477–86.View ArticleGoogle Scholar
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.