DNA sampling and sequencing

A total of 57 Oriental greenfinch tissue samples were obtained from 19 locations across the entire breeding range (Fig. 1; Table 2; see  Supplementary Table S1 (10.2108.zsj.37.280.s1.pdf) online), mostly during the breeding season. We also obtained old skin samples from Chinese bird populations but were unable to generate DNA sequence data from them. Instead, we used a C. s. sinica DNA sequence from a sample collected in Szechuan, China (accession no. L76592; Arnaiz-Villena et al., 1998). Tissue samples were also obtained from the closest relative of the Oriental greenfinch—the European greenfinch C. chloris (accession nos. LC485969 and LC485970) as an outgroup species.

DNA was extracted from muscle and blood samples using a DNeasy Blood & Tissue Kit (Qiagen, Tokyo, Japan), in accordance with the manufacturer's instructions. We amplified and sequenced 846 bp of the cytochrome b (Cytb) gene, 389 bp of the NADH subunit 6 (ND6) gene, 69 bp of the tRNA-Glu gene, and 1034 bp of the control region (CR) of the mtDNA. Cytb was amplified via polymerase chain reaction (PCR) and sequenced using the primer set L14841 (Kocher et al., 1989) and H15767 (Edwards et al., 1991). Each reaction mixture included 0.5 µL of purified total DNA template, 2 µL of dNTP mixture (2.5 mM each), 2.5 µL of 10× Ex Taq Buffer (Takara Bio, Inc., Shiga, Japan), 0.3 µL of each primer (0.24 pmol/µL), and 0.36 µL of Ex Taq (1.25 units/µL; Takara Bio). The PCR conditions consisted of an initial denaturation step at 95°C for 3 min; 35 cycles of 95°C for 30 s, 53°C for 30 s, and 72°C for 45 s; and a final extension step at 72°C for 5 min. The other mtDNA regions (ND6, tRNA-Glu and CR) were amplified using the primer set CRTPPO (5′-CCATCTCCAACTCCCAAAGC-3′; P. Boag, personal communication) and H1261Cs (5′-GAAGATGTCAAGATGGCTGCC-3′; modified from Mindell, 1997). Sequencing primers used were CRTPPO, H1261Cs, MATs (5′-CCATTGTCCCCTCCAGGCGC-3′; P. Boag, personal communication), F304, and F389 (Mindell, 1997), CRTPPO and F389 for the amplification of ND6 to CR, H1261Cs, MATs and F304 for the amplification of CR. The reagent concentrations in each reaction mixture were the same as those used to amplify Cytb. The PCR conditions consisted of an initial denaturation step at 95°C for 3 min; 35 cycles of 95°C for 30 s, 67°C for 30 s, and 72°C for 1 min; and a final extension step at 72°C for 5 min. All reactions were performed in a 25 µL volume using a Veriti Thermal Cycler (Applied Biosystems, CA, USA). The PCR products were purified using an illustra ExoProStar PCR Clean-up kit (GE Healthcare, UK). The resulting products were sequenced using a BigDye Terminator v. 3.1 Cycle Sequencing kit and visualized on an ABI PRISM 3130xl Genetic Analyzer (Applied Biosystems). All sequences were deposited in GenBank (see  Supplementary Table S1 (10.2108.zsj.37.280.s1.pdf) online).

Phylogenetic and dating analyses

There were no gaps in the target regions. Therefore, sequences were aligned by eye using GENETYX-MAC ver. 19.0.1 (GENETYX Corporation, Tokyo, Japan). Substitution models were selected based on the Akaike information criterion calculated using jModelTest ver. 2.1.10 (Darriba et al., 2012). Inter- and intra-specific phylogenetic relationships were evaluated using the concatenated sequences (Cytb, ND6, tRNA-Glu, and CR, hereafter referred to as Cytb–ND6–tRNAGlu–CR). Phylogenetic trees were determined for all concatenated sequences via Bayesian inference using BEAST ver. 1.10.4 (Suchard et al., 2018) with the Hasegawa, Kishino and Yano model (HKY model, Hasegawa et al., 1985), a strict clock, and constant-size coalescent tree priors. To construct phylogenetic trees, we used the European greenfinch C. chloris as an outgroup species. Chloris chloris is the species most closely related to C. sinica, based on the topologies of the Cytb tree published by Arnaiz-Villena et al. (1998, 2008) and the ND2/ND3 tree published by Zuccon et al. (2012).

For the dating analysis, we used only the Cytb dataset. The molecular clock for Cytb has been calibrated to approximately 2.1% sequence divergence per million years, across most orders of birds (Weir and Schluter, 2008). Therefore, we used the GTR + Γ model for Cytb, a strict clock with a mean rate of 0.0105 substitutions/site/ million years, based on the calibration method (Weir and Schluter, 2008), and constant-size coalescent tree priors. To calculate divergence times among taxa, we used sequences of the following outgroup species: C. chloris, the Black-headed greenfinch C. ambigua (accession no. U78322), and the Yellow-breasted greenfinch C. spinoides (accession no. U79018; Arnaiz-Villena et al., 1998).

All analyses were run for 150 million generations and sampled every 1000 generations. The first 25% of the trees generated were discarded as “burn-in”. Good mixing of the Markov chain Monte Carlo parameters and reproducibility were ensured via multiple runs from independent starting points. Trees were summarized using TreeAnnotator ver. 1.10.4 (included in the BEAST package) with the “Maximum clade credibility tree” and “Mean heights” settings, and then displayed using FigTree ver. 1.4.4 (Rambaut, 2012).

In addition to phylogenetic analysis using BEAST software, we used MrBayes ver. 3.1 (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003) to estimate a posterior probability (PP) at each tree node. We used the HKY model and default priors for the concatenated sequences (Cytb–ND6–tRNA-Glu–CR). Four Metropolis-coupled Markov chain Monte Carlo chains were run for 10⁷ generations and sampled every 1000 generations. The average standard deviation of split frequencies passing below 0.01 and the potential scale reduction factor were close to 1.00 for all parameters.

To construct a network diagram for the haplotypes of the four concatenated sequences, we used Network ver. 5.0.1.1 (Bandelt et al., 1999;  http://www.fluxus-engineering.com) with the Median Joining algorithm, and the post-processing calculation was set to the Maximal Parsimony option. These network diagrams can often reveal intraspecies affinities and homoplasies.

For the Cytb, ND6, and CR data, the number of haplotypes (h), haplotype diversity (hd), and nucleotide diversity (π) were calculated, and neutrality tests were performed, using DnaSP ver. 5.10 (Librado and Rozas, 2009). We did not use tRNA-Glu data for genetic diversity analysis, because it was revealed to have only two haplotypes.

Morphometric analyses

A total of 114 museum specimens of the Oriental greenfinch had been collected between 1884 and 2014 across the eastern part of the species' breeding range (see  Supplementary Table S2 (10.2108.zsj.37.280.s1.pdf) online). This set of specimens was not linked to the DNA samples. All of the measurements recorded from the museum specimens kept at the Yamashina Institute for Ornithology (YIO; Abiko, Japan) and the Natural History Museum (NHM; Tring, UK), including type specimens, are listed in  Supplementary Table S2 (10.2108.zsj.37.280.s1.pdf) online. The morphometric analyses were based entirely on adult male specimens, because few female and juvenile specimens were available. Using digital calipers, we measured natural wing length (NW), tail length (TAIL), tarsus length (TAR), bill height (BH) and width (BW) at the posterior edge of the nostrils, gape width (GW), total culmen (bill length from the skull, TC), total head length (TH), and the distance between the longest tertiary and the longest primary feather on the folded wing (MP–T; Fig. 2; Mitutoyo Corporation, Kanagawa, Japan). In addition, we obtained some body weight measurements from specimen labels and converted those recorded using old Japanese weight units (Monme) into grams.

Fig. 2.

Museum specimen measurements. The code for each measurement is shown in Table 5.

To identify external morphological differences among populations, we tested for statistically significant differences among mean values for each of the measurements using one-way analysis of variance (ANOVA) and Scheffe's post-hoc test. To perform the calculations for these tests, morphometric values were transformed using natural logarithms. These analyses were performed using Mac Statistical Analysis ver. 2.0 (Esumi, Tokyo, Japan). In addition, we performed principal component analysis (PCA) using Mac Excel Multivariate Analysis ver. 1.0b (Esumi). The nine characteristics included in the PCA were NW, TAIL, TAR, BH, BW, GW, TC, TH, and MP–T. The principal component scores were extracted from the correlation matrix.

Phylogeny and divergence times

Sequencing Cytb, ND6, tRNA-Glu, and CR markers from 57 individual Oriental greenfinches (see  Supplementary Table S1 (10.2108.zsj.37.280.s1.pdf) online) produced 15, 5, 2, and 12 haplotypes, respectively (Table 3). Cytb was sequenced from all of the 57 samples, whereas the full Cytb–ND6–tRNA-Glu–CR was determined only for 33 specimens (see  Supplementary Table S1 (10.2108.zsj.37.280.s1.pdf) online). These sequences were deposited in the DDBJ nucleotide database (accession nos. LC485304–LC487318; Table 3).

Inter- and intra-specific phylogenetic relationships were evaluated with BEAST and MrBayes software using the concatenated Cytb–ND6–tRNA-Glu–CR sequences from 33 samples shown in Fig. 3 and in  Supplementary Figure S1 (10.2108.zsj.37.280.s1.pdf) online. The HKY model was selected for the concatenated sequences using jModelTest software. The Oriental greenfinch formed a strongly supported clade that was separated from C. chloris. The Oriental greenfinch clade was split into two strongly supported clades, with the population from the Ogasawara Islands (Clade B) diverging from the other populations (Clade A) (Fig. 3). Clade A was also strongly supported (BEAST/MrBayes PP ≥ 0.95) and was split into two subclades: one comprised the southern Primorsky populations (PRI_RYA334 and PRI_RYA336) and the other comprised the remaining populations. The latter subclade was only supported by PP in the BEAST analysis (MrBayes PP < 0.95) and was split into two subclades: one consisted of the Ulleung-do Island populations (ULL_2802 and ULL_3546) and the other consisted of the remaining populations. However, the Ulleung-do subclade was placed in another clade by the MrBayes analysis (see  Supplementary Figure S1 (10.2108.zsj.37.280.s1.pdf) online). The remaining subclade in Clade A was only supported by the BEAST analysis, and was further subdivided into at least three subclades. However, these subclades did not correspond to particular geographical regions (Fig. 3, see  Supplementary Figure S1 (10.2108.zsj.37.280.s1.pdf) online).

Table 3.

Variable sites in the mtDNA haplotypes of the Oriental greenfinch. The dots indicate bases that are identical to those in haplotype h1 (2003-5103), and the frequency of each haplotype within each population is indicated. Site numbers are shown for each cytochrome b, NADH subunit 6, tRNA-Glu, and control region sequence. The area codes for DNA samples in each sample are listed in Table 2. DNA sequences have been deposited in the DDBJ nucleotide database (accession nos. LC485304–LC487318).

Fig. 3.

Phylogeny of 33 Oriental greenfinches based on the mitochondrial cytochrome b, NADH subunit 6, tRNA-Glu, and the control region inferred by BEAST software. The European greenfinch Chloris chloris was used as an outgroup species. Molecular phylogenies were determined using Bayesian inference according to the Hasegawa, Kishino and Yano model using concatenated DNA encoding mitochondrial Cytb, ND6, tRNA-Glu and CR. Support values at nodes indicate posterior probabilities (PPs) as determined by BEAST/MrBayes. * PP ≥ 0.95. Area codes are listed in Table 2.

The ages of the major clades were estimated with BEAST software using the Cytb sequences from 58 samples (57 sequences generated in the present study in addition to one sequence of C. s. sinica reported by Arnaiz-Villena et al., 1998, Fig. 4). We were unable to determine the relationships among the outgroup taxa and C. sinica due to a low PP value (0.57) at base node of the C. chloris and outgroup taxa. Therefore, we conclude that the outgroup species and C. sinica diverged either 1.38 or 1.75 million years ago (MYA), with 95% highest posterior density (HPD) intervals of 0.85–1.94 or 1.25–2.36 MYA, respectively. This divergence would have occurred within the Calabrian (0.78–1.80 MYA) or Gelasian (1.80–2.58 MYA) age, during the Pleistocene epoch (Cohen et al., 2013). The deepest split between Clade A and Clade B within C. sinica occurred approximately 1.06 MYA, with a 95% HPD interval of 0.63–1.59 MYA, which falls between the Chibanian (0.13–0.78 MYA) and the Calabrian ages. This estimate is much older than the timing of the split between C. ambigua and C. spinoides (0.59 MYA; 95% HPD interval, 0.29–0.96 MYA; Fig. 4). The divergence times of the subclades in Clade A are less than 0.28 MYA (95% HPD interval, 0.14–0.47 MYA), which fits to the Chibanian Pleistocene age.

A median joining network diagram was constructed to show the optimal connections between the haplotypes of C. sinica and C. chloris (Fig. 5). The topology of this network diagram was almost identical to that of the concatenated tree (Fig. 3), except at the base. There were at least 48 substitutions between the Primorsky and Ogasawara bird populations, via several unobserved haplotypes and at least 66 substitutions between C. sinica and C. chloris.

Because the sample size in Clade A was much larger than that in Clade B, the genetic diversity in these two clades was different: the former had higher diversity values than the latter. Few haplotypes were identified in Clade B (representing subspecies kittlitzi): only one haplotype was detected for Cytb and ND6, and three haplotypes were detected for the CR (Table 3). Therefore, the haplotype diversity (hd) and nucleotide diversity (π) of Cytb and ND6 were 0, whereas those of the CR were 0.60 (Table 4). In Clade A (which included all subspecies except kittlitzi), the haplotype diversity of Cytb and the CR was high (0.622 and 0.801, respectively), whereas that of ND6 was low (0.276). For all genes, nucleotide diversity within both clades was low, supporting the possibility of a sudden population expansion. Tajima's D, Fu's Fs, and Fu and Li's D values of CR for both clades were negative, which also would indicate a sudden population expansion; however, the values were only significant for Cytb within Clade A (Table 4). At the subspecific level, we could not compare the degree of genetic diversity among subspecies, because the sample size of each subspecies was too small.

Fig. 4.

Chronogram for 58 Oriental greenfinches inferred by BEAST software and a 2.1%/million year molecular clock for cytochrome b. Values above the nodes indicate PPs (* PP ≥ 0.95). Bayesian PP values > 0.5 are shown above the nodes. Values in parentheses below the nodes indicate divergence times from the present in million years. All branch lengths are scaled according to time, and 95% highest posterior density age estimate intervals are shown as shaded bars on the nodes. Area codes for the DNA samples are listed in Table 2.

Fig. 5.

Median-joining network diagram for Chloris sinica and C. chloris. The network diagram is based on haplotypes of concatenated sequences (Cytb, ND6, tRNA-Glu, and CR). Circle size denotes sample size, and numbers on the midline denote the number of base substitutions.

Morphology

The measurements recorded for the 114 adult male Oriental greenfinch museum specimens are presented in Table 5. Significant differences among the mean measurements from each population are shown in Fig. 6A–G. The longest and shortest mean NW measurements were recorded for birds from the Kuril Islands (87.36 ± 1.74 mm, n = 11) and Ogasawara Islands (74.35 ± 1.52 mm, n = 13), respectively (Table 5). There was a significant difference of NW between these populations (Scheffe's post-hoc test, P < 0.01; Fig. 6A). The mean NW of the population from the Kuril Islands was significantly longer (F_(9,111) = 31.5, P < 0.01; Fig. 6A), and that of the population from the Ogasawara Islands significantly shorter (Scheffe's post-hoc test, P < 0.01; Fig. 6A), than those of each of the other populations. Furthermore, the northern populations exhibited a trend for a longer NW compared with the southern populations (Fig. 6A, Table 5).

Table 4.

Indices of genetic diversity, Tajima's D, Fu's Fs and Fu and Li's D tests of neutrality. Abbreviations: n, number of samples; h, number of haplotypes; hd, haplotype diversity with variance; π, nucleotide diversity with variance.

Although the bird population from the Ogasawara Islands had the shortest NW, it had the longest mean TC (13.25 ± 0.84 mm, n = 13; F_(9,111) = 9.6, P < 0.01; Fig. 6B). This was significantly longer than the TC measurements recorded for birds from northeastern China, mainland Korea, Honshu (Scheffe's post-hoc test, P < 0.01), and Cheju Island (Scheffe's post-hoc test, P < 0.05). The bird population from the Kuril Islands had a greater mean BH than the other populations (F_(9,111) = 22.1, P < 0.01; Fig. 6C), with the exception of bird populations from southeastern China (Scheffe's post-hoc test, P < 0.01; Fig. 6C). By contrast, bird populations from northeastern China had a shorter BH compared to populations from Sakhalin Island, Hokkaido, southeastern China, and the Ogasawara Islands (Scheffe's post-hoc test, P < 0.01). Bird populations from the Kuril Islands had wider BW (F_(9,111) = 8.9, P < 0.01; Fig. 6D) and GW (F_(9,111) = 9.0, P < 0.01; Fig. 6E) than those of other populations. Bird populations from northeastern China had a shorter TH than those of most other populations (F_(9,111) = 9.6, P < 0.01; Fig. 6F).

Thus, the bird populations from the Kuril Islands and Sakhalin Island had deep and wide bills, whereas the birds from northeastern China had short and narrow bills, and those from the Ogasawara Islands had long bills (Fig. 6A–F).

The bird population from the Ogasawara Islands had significantly shorter MP–T lengths than the other bird populations (F_(9,111) = 6.8, P < 0.01; Fig. 6G). The trend of MP–T length differences among populations was similar to the trend shown by NW, it was also associated with size indicators, as was NW. We observed some significant differences among populations but no remarkable trends with regard to TAIL (F_(9,111) = 4.8, P < 0.01; not shown) or TAR (F_(9,111) = 4.5, P < 0.01; not shown) measurements.

Table 5.

Measurements of Oriental greenfinch males recorded from the Yamashina Institute for Ornithology museum specimens. Sample sizes, means, and standard deviations are shown. Lower values in each line indicate the range of measurements.

Fig. 6.

Comparison of measurements from specimens of Oriental greenfinch males. For each parameter, a mean value was calculated from nine measurements (mm): (A), natural wing length; (B), bill length; (C), bill height; (D), bill width; (E), gape width; (F), total head length; and (G), distance between the longest tertiary and the midpoint of the primaries on the folded wing. Medians are indicated by filled circles. Maximum, minimum, lower, and upper quartile values, as well as outliers are indicated by filled squares. The tail length and tarsus length data are not shown. ** P < 0.01, * P < 0.05, one-way analysis of variance or Scheffe's test. Area codes for the measurements are shown at the bottom of the figure and are listed in Table 5.

Fig. 7.

Principal Component Analysis (PCA) of Oriental greenfinch males. The PCA was based on measurements of nine characteristics recorded from museum specimens of Oriental greenfinch males, including type specimens of Chloris sinica kittlitzi (specimen no. 98.11.1.60, BNHM), C. s. sitchitoensis (no. YIO-00110) and C. s. clarki (no. YIO-00132). Each type specimen is shown by a symbol and an arrow. Area codes with the symbols are listed in Table 5.

Table 6.

Principal component scores for the three principal component axes based on nine characteristics.

For body weight, we were only able to read a few specimen labels from each population and could not obtain data for the bird populations from Hokkaido or Ogasawara Islands (Table 5). The birds from the Kuril Islands population were the heaviest (26.3 ± 2.22 g, n = 7), while those from northeastern China were the lightest (18.8 ± 2.01 g, n = 7).

The PCA supported the conclusions described above (Fig. 7). The contributions of the principal components PC1, PC2, and PC3 were 41.28%, 18.21%, and 9.49%, respectively (Table 6). Although the total cumulative variance among the three components was low (68.98%), we can use PCA to understand the major morphological characteristics of each population. For PC1, the top four factor loadings were related to bill morphology, and larger positive y-axis values indicated taller and wider bills. By contrast, larger positive x-axis values for PC2 indicated larger wings, and more negative x-axis values indicated longer culmina (Fig. 7; Table 6).

A scatter plot of measurements from a type specimen of C. s. kittlitzi (no. 98.11.1.60, BMNH) fell within the range of measurements obtained for birds from the Ogasawara Islands (Fig. 7). This result suggests that the morphology of C. s. kittlitzi is consistent with it being a subspecies of C. sinica from the Ogasawara Islands. By contrast, a scatter plot of measurements from a type specimen of C. s. sitchitoensis (YIO-00110), which originated from a wintering ground on Hachijojima Island, fell within the range of measurements obtained for birds from the Kuril Islands and Sakhalin Island (Fig. 7). Thus, the breeding population that produced this specimen probably originated in the Kuril Islands or Sakhalin Island, not Hachijojima Island in the Izu Islands.

Divergence times and relationships within the species

We found a highly diverse and well-supported clade within the Oriental greenfinch species, corresponding to a population of birds from the Ogasawara Islands. This clade diverged from the other Oriental greenfinch populations approximately 1.06 MYA, during the Calabrian age of the Pleistocene epoch (Fig. 4). This divergence date is 1.8-fold more ancient than the split between C. ambigua and C. spinoides (0.59 MYA; Fig. 4), which fits to a species-level divergence. This very old divergence date for the Oriental greenfinch clade indicates that the populations on the Ogasawara Islands were isolated for a long time, accumulating genetic differences and acquiring unique morphological characteristics. Chloris sinica diverged from its closest continental relatives 1.38–1.75 MYA. Avian species from oceanic islands that diverged from their continental ancestors in the distant past also reportedly include Darwin's finches on the Galapagos Islands (2–3 MYA; Grant and Grant, 2008), Hawaiian Honeycreepers on the Hawaiian Islands (7.2 MYA; Heather et al., 2011), and the dark-eyed junco Junco hyemalis on Guadalupe Island (0.6 MYA; Aleixandre et al., 2013).

Understanding this ancient divergence within the Oriental greenfinch species would be easier if we knew the estimated divergence times between other extant species (or subspecies) and their close relatives on the Ogasawara Islands. This would enable us to reconstruct the phylogeographical history of the islands' avifauna. Unfortunately, there are no records for many endemic species and subspecies from the Ogasawara Islands for the last 50 years, and these populations may be currently extinct on the islands. These include the Bonin wood pigeon Columba versicolor, the Bonin thrush Cichlopasser terrestris, the Bonin grosbeak Chaunoproctus ferreorostris, the Nankeen night heron Nycticorax caledonicus crassirostris, the white-browed crake Pozana cinerea brevipes, the peregrine falcon Falco peregrinus furuitii, and the Bonin white-eye Apalopteron familiare familiare (The Ornithological Society of Japan, 2012). Therefore, we cannot examine these taxa using live birds. However, DNA extracted from museum specimens may be useful to study the molecular phylogeny of these species, and such research is ongoing.

Several studies have revealed the molecular phylogeny of some endemic extant subspecies from the Ogasawara Islands using mtDNA. Genetic differences or divergence times between the subspecies living on the Ogasawara Islands and other regions are as follows: the Japanese wood pigeon Columba janthina nitens and other subspecies (CR, 3.56%; Seki et al., 2007), the Japanese bush warbler Cettia diphone diphone and other subspecies (COI, 0.2 MYA; Emura et al., 2013) and the brown-eared bulbul Hypsipetes amaurotis squamiceps and H. a. ogawae, H. a. pryeri and H. a. borodinonis (COI, 2.53–3.19%; Saitoh et al., 2015). These studies have shown that there are fewer genetic differences or shorter divergence times among the subspecies living beyond the Ogasawara Islands than we observed between Clades A and B of the Oriental greenfinch (1.06 MYA; Fig. 4), except in the brown-eared bulbul, which shows large genetic differences within the species. It is unclear why intraspecific divergence occurred in the Oriental greenfinch so long ago, but we know it happened during the Calabrian age in the Pleistocene epoch (1.80–0.78 MYA), when many Palearctic bird species also emerged (Finlayson, 2011). During the Pliocene (5.3–2.6 MYA), the climate of the Palearctic became more arid. Recently, species have often emerged due to habitat fragmentation caused by the opening up of forested regions and the expansion of desert regions (Finlayson, 2011). After the Pliocene, the Pleistocene (2.58–0.01 MYA) climate cooled, and lineages of closely related species of Palearctic birds continued to split (Finlayson, 2011). The ancestral Oriental greenfinch may have been affected by such changes in the continental climate, leading to a split in its lineage.

Clade A is strongly supported and subdiverged into three minor clades, one of these clades consisting of two individuals from southern Primorsky (PRI_RYA334 and PRI_RYA336), another clade consisting of two individuals from Ulleung-do Island (ULL_2802 and ULL_3546), and the other consisting of the remaining populations (Fig. 3). However, the Ulleung-do individuals were completely nested in the last clade in the MrBayes analysis (see  Supplementary Figure S1 (10.2108.zsj.37.280.s1.pdf) online). Therefore, we were unable to find any associations between genetic structure and distribution range for the subspecies. This may be due to a sudden expansion in population from a few founder individuals after 0.28 MYA, during the Chibanian Pleistocene (Fig. 4; Table 4). The existing haplotypes may have originated in the Upper–Chibanian Pleistocene (Cohen et al., 2013), because this was when several extant populations of Far East Asian songbirds diverged, including the brown-eared bulbul H. amaurotis (0.09 MYA), the fan-tailed warbler Cisticola juncidis (0.17 MYA), the jungle crow Corvus macrorhynchos (0.10 MYA; Nishiumi et al., 2006), and the Siberian rubythroat Luscinia calliope (0.76 MYA; Spiridonova et al., 2013).

Origin of the Ogasawara bird populations

The Ogasawara Islands originated as a result of volcanic activity approximately 48 MYA and appeared above sea level during the Eocene epoch (Umino et al., 2016). Therefore, it is difficult to know when these islands were colonized by different bird species, with the exception of those species that have only appeared there for a few decades, such as the common moorhen Gallinula chloropus, the bull-headed shrike Lanius bucephalus, and the scaly thrush Zoothera dauma (Kawakami, 2019).

Nakamura (2002) suggested that the Ogasawara Oriental greenfinch population originated as the result of an accidental invasion from a northern population located in Kamchatka, the Kuril Islands, or Hokkaido. A deviation from the normal migration route may have resulted in birds straying onto the Izu Islands during autumn and colonizing these islands, before finally reaching the Ogasawara Islands. However, this hypothesis is not supported by our genetic data for two reasons. First, the Ogasawara population diverged from the common ancestor of C. s. kittlitzi and the other subspecies, but did not originate from extant continental taxa (Figs. 3 and 4). Second, the Miyake Island (Izu Islands) population belongs to a different clade than the Ogasawara population (Figs. 3 and 4). Although we do not know where the two Oriental greenfinch clades diverged, the sequence data strongly suggest that a common ancestor of the Oriental greenfinch and its closest relatives further diverged into two clades, with Clade B (Ogasawara Islands) dispersing to the remote oceanic islands and Clade A dispersing to the remaining range inhabited by several extant subspecies. The phylogenetic network (Fig. 5) clearly shows close affinity of Primorsky mainland population for that of Ogasawara Islands. This however cannot be directly interpreted in terms of origination, but this would be considered in further investigations with more data.

Seki et al. (2012) reported similar results as this study. They identified two clades within the Japanese robin Luscinia akahige, using mtDNA sampled across its breeding range in the Japanese Archipelago, that had diverged in the distant past (1.1 MYA). One of these clades is endemic to the Izu Islands, whereas the other nests on mainland Japan. Those two mtDNA lineages probably originated in different refugia during the mid-Pleistocene glacial period, one in southwestern Japan and the other on the Izu Islands. Similarly, the two distinct Oriental greenfinch clades may have evolved in different locations during the Pleistocene epoch, with one clade emerging from the continental Far East or the Japanese archipelago and the other emerging from the Ogasawara Islands. Throughout the glacial periods, the distribution of these finches would have been restricted. Thereafter, populations from the northern refugium may have expanded north, whereas those from the Ogasawara refugium may have remained at their present location, unable to spread north due to their local adaptations.

Morphological differences among Oriental greenfinch populations

Our morphological analyses showed that Oriental greenfinches from the Kuril and Ogasawara Islands had the longest and shortest mean wing lengths, respectively (Fig. 6A; Table 5). In general, wing length is the most reliable indicator of body size for intraspecific comparisons in birds (Lanyon, 1960; Selander and Johnson, 1967). If so, bird populations from the Kuril Islands (subspecies kawarahiba) and Ogasawara Islands (subspecies kittlitzi) should have the largest and smallest bodies, respectively. However, other study argued that body weight is a more reliable indicator of size than wing length (James and Somers, 1989). We were unable to obtain sufficient body weight data from the museum specimens. Therefore, we compared body weights from two populations only. The mean body weight of the Kamchatka population was 29.1 ± 1.7 g (four adult males; Sugawa et al., 2010), whereas that of the Ogasawara population was 18.09 ± 0.99 g (53 adult males; K. Kawakami, unpublished data). Therefore, the body sizes of these two populations were different.

Across the greenfinch populations, we found a moderate trend in size that was associated with latitude. The northern populations had longer wings than the southern populations (Fig. 6A, G; Table 5). Although only Japanese populations were analyzed, a similar latitudinal trend in body weight, wing and tail length was reported for Oriental greenfinch males ranging from the northern populations of Hokkaido (43°50′N) to the southern Ogasawara Islands (26°35′N; Nakamura, 1997). He observed that the mean body weights of birds inhabiting more northerly latitudes increased gradually from the Ogasawara Islands (17.7 g, n = 12) to Nagano Prefecture (Honshu; 20.5 g, n = 30), and Hokkaido (23.0 g, n = 44). If the mean body weight of C. s. kawarahiba (29.1 g) is included in the Japanese finch dataset, then the latitudinal trend for increasing body size in the Oriental greenfinch progresses along the Pacific island chain, through the Japanese archipelago and the Kuril Islands, to the Kamchatka Peninsula. This trend in body size is consistent with Bergmann's rule (Bergmann, 1847). Several studies have re-examined Bergmann's rule using bird populations and found that it is supported by intraspecific (Ashton, 2008) and interspecific data (Olson et al., 2009). This means, higher latitudes and cooler climates tend to favor larger bodies.

However, the continental Oriental greenfinch populations of Korea (C. s. ussuriensis), northeastern China (C. s. chaborovi), and southeastern China (C. s. sinica) are relatively small compared to the insular subspecies (Fig. 6; Table 5). In particular, the northeastern China bird population clearly had the smallest bill (Fig. 6B–F). The mean body weight of the birds from northeastern China (18.8 g, n = 7) was also low (Table 5), although these birds were slightly heavier than those from the Ogasawara Islands (17.7 g, n = 12, Nakamura, 1997; 18.09 g, n = 53, K. Kawakami, unpublished data).

Our morphological analysis suggested that C. s. kawarahiba (Kuril Islands), C. s. kittlitzi (Ogasawara Islands), and C. s. chabarovi (northeastern China) can be identified based on their anatomical measurements, which were significantly different from those recorded from other populations (Fig. 6; Table 5). However, it was difficult to identify the remaining subspecies based on their anatomical measurements because the morphological differences between neighboring subspecies were often small (Fig. 7).

Morphological characteristics of the insular subspecies C. s. kittlitzi

Chloris sinica kittlitzi, which inhabits the Ogasawara Islands, has a long TC, despite its small body (Figs. 6 and 7; Table 5). Nakamura (1997) reported that the C. s. kittlitzi bill was longer, deeper, and thicker than that of C. s. minor. Probably, the C. s. kittlitzi bill is specifically adapted for feeding on the Ogasawara Islands. Nakamura (1997) studied feeding behavior and surveyed the crop contents of finches on Hahajima Island (part of the Ogasawara Islands) and showed that the finches frequently fed their chicks with seeds from trees, especially Wikstroemia pseudoretusa. Chloris s. kittlitzi inhabits mountain ridges and feeds on arid shrubs, including W. pseudoretusa. The seeds of W. pseudoretusa are much bigger than grass seeds, and C. s. kittlitzi probably developed its large bill as an adaptation for feeding on these seeds (Nakamura, 1997).

A similar morphological adaptation was reported by Aleixandre et al. (2013). The dark-eyed junco J. hyemalis is a widespread songbird species from North America. It comprises several subspecies, including the insular subspecies J. h. insularis, which is restricted to Guadalupe Island, 257 km from Baja California. This island population consists of diverse genetic lineage that differ morphologically from mainland juncos. In particular, these birds have a large bill and small body. The island juncos have adapted to feeding on large seeds within hard cones of the endemic cypress Cupressus guadalupensis. Increased beak size has been observed in other insular passerines and is commonly associated with a wider range of food resources and fewer competing species (Illera et al., 2007). The Japanese bush warbler Cettia diphone is another insular passerine with adaptations of its bill size and wing length. The insular subspecies Ce. d. diphone breeds on the Ogasawara Islands (The Ornithological Society of Japan, 2012) and has a shorter tail and wings, and a longer and narrower beak, than its conspecifics on mainland Japan and the Izu Islands (Morioka, 1977).

Taxonomic recommendations

Our study has shown that the subspecies C. s. kittlitzi breeding on the Ogasawara Islands has diverged significantly by its mtDNA from that of all other Oriental greenfinch populations. Moreover, it has distinctive morphological characteristics, including short wings, a long bill, and a small body. In addition, it has a smaller mean clutch size and larger eggs than other Japanese Oriental greenfinch populations (Nakamura, 1997). This subspecies has become adapted to its subtropical island environment.

Therefore, we recommend that the subspecies C. s. kittlitzi should be classified as a distinct species, Chloris kittlitzi (originally Fringilla kittlitzi Seebohm, 1890). In addition, we propose that C. kittlitzi is named ‘the Ogasawara greenfinch.'

By contrast, we found no divergent lineage within other mtDNA clades that could be considered as a separate species. Therefore, all other subspecies should retain their present status (Fig. 3). Moreover, we were unable to identify clear relationships between differences in mtDNA sequences and geographical locations of the various subspecies within Clade A (Fig. 3). Different molecular markers, such as microsatellite sequences within the nuclear genome, may be used to improve our understanding of these relationships.

Conservation of the Ogasawara greenfinch

The Ogasawara greenfinch C. kittlitzi is now critically endangered on the Ogasawara Islands. The estimated bird population size on the Hahajima and Volcano Islands is approximately 200–400 birds. Although these estimates vary, it is clear that the risk of extinction is increasing year by year.

The black rat Rattus rattus, is not found within the breeding range of the Ogasawara greenfinch. This is fortunate, because populations of finches would be unable to survive on islands inhabited by the black rat (Kawakami, 2019). Predation by the black rat is probably the major factor in the local extinction of finch populations. The Norway rat, R. norvegicus, has invaded the satellite islands of Hahajima Island and likely preys on finch nests (Kawakami, 2019). Another threat to the finches is predation by the domestic cat, Felis catus, finch feathers have been found in cat droppings at Hahajima Island during non-breeding season (Kawakami and Higuchi, 2002).

In addition to these threats, eradication of the alien tree Casuarina equisetifolia may also affect the Ogasawara greenfinch population. Finches frequently nest in these trees and feed on their seeds (Nakamura, 1997; Kanto Regional Forest Office, 2009). However, C. equisetifolia has been eradicated from the Ogasawara Islands because it forms dense monospecific stands and reduces biodiversity (Kawakami, 2019).

To reduce the risk of the Ogasawara greenfinch becoming extinct, urgent action is necessary to conserve the species' present habitat and provide rearing facilities (e.g., in a zoo).