Phylogeography

Samples and gene regions - this study aims to provide more detailed demographic analysis of the Eurasian green woodpecker based on previously published multi-locus data. It includes mitochondrial DNA (mtDNA) haplotypes of the cytochrome b (cyt-b) and the NADH dehydrogenase 2 (ND2) genes. Sample locations and haplotypes, together with their frequencies, were given in Pons et al. (2010) for the cyt-b gene, and in Perktas et al. (2011) for the ND2 gene. The total sample size of the Eurasian green woodpecker for the phylogenetic analysis was 123 for mtDNA loci. Sample locations for mitochondrial loci were shown in Fig. 1.

Fig. 1.

Distribution of the Eurasian green woodpecker (Picus viridis) in Europe, Anatolia, and Caucasus. Genetic structure of populations is indicated by pie charts (A, the cyt-b gene; B, the ND2 gene); the area of each pie chart is proportional to the sampling size for each population. Both figures include haplotype networks. In figure B, scatter dot plot shows nucleotide (pi) and haplotype diversity (hd) based on the ND2 gene for each sampling areas.

Genetic variation and phylogenetics — standard measures of genetic variation, such as number of haplotypes, haplotype diversity, and nucleotide diversity, were estimated in Pons et al. (2010) and Perktas et al. (2011) for mtDNA loci. Thus, they were not estimated once again in this paper. Haplotype networks were constructed using PopART (Population Analysis with Reticulate Trees) version 1.7.2 ( www.popart.otago.ac.nz) with median joining network criteria implemented in NETWORK (Bandelt et al. 1999). The median joining network combines the minimum-spanning trees within a single haplotype network and yields the best gene genealogy. Therefore, it provides a good methodological approach to investigate phylogenetic relationships between mtDNA haplotypes and also to describe genetic population structure of the Eurasian green woodpecker. On the basis of the results of haplotype network of mtDNA loci, we concluded that there was no phylogeographic structure among the Eurasian green woodpecker populations. In addition, Pons et al. (2010) used nuclear BRM intron in their phylogenetic analyses. They did not find significant variation within distribution range of the Eurasian green woodpecker. Thus, all sequences from all populations were combined to explore the historical demography of the Eurasian green woodpecker. This approach made the results of phylogeography more comparable with those of ecological niche modelling (e.g. Perktaş et al. 2015).

Historical demography — to infer the demographic history of the Eurasian green woodpecker, we used the extended Bayesian skyline plot analysis (Heled & Drummond 2008) implemented in BEAST version 1.7.5 (Drummond et al. 2012). This analysis uses coalescent approach to estimate Bayesian inference of the effective population size changes through the time based on the information of multiple mtDNA loci. Before the extended Bayesian skyline plot runs, the best-fit substitution models were identified for each locus in MEGA version 6.06 (Tamura et al. 2013). These were the Tamura & Nei (1993) with a discrete Gamma distribution (TN93 + G, AIC_c = 3462.72) for the ND2 gene, Hasegawa, Kishino and Yano with a discrete Gamma distribution (HKY + G, AICc = 1466.53) for the cyt-b gene. Because we focused on the different mtDNA data sets of the Eurasian green woodpecker, we first performed independent separate runs for each dataset (the cyt-b gene and the ND2 gene) using extended Bayesian skyline plot (Heled & Drummond 2008) implemented in BEAST version 1.7.5. Because of the fact that results of separate extended Bayesian skyline plots were congruent with each other, we used both mtDNA loci (the cyt-b and the ND2) in one extended Bayesian skyline plot analysis. Multiple independent extended Bayesian skyline plot runs were performed using the following parameters: linear models, 100 million steps, parameters sampled every 10000 steps, and a burn in of 10 %. We used strict clock model with default mutation rate under uniform prior distribution, and allowed the analysis to estimate the rates, relative to the cyt-b gene, for the ND2 gene. Then, we calculated the expansion time using a wide range of mutation rates of mtDNA, i.e. 1.65 % substitutions/site/million years as the mutation rate of the mtDNA cyt-b gene of Piciformes (Weir & Schluter 2008), and the range of 1–4 % (Brito 2005). Effective sample size values of the parameters were over 200 for each run.

Isolation-by-distance - a matrix of genetic distances between all pairs of populations was estimated from standard fixation index (F_st) for two mtDNA loci (the cyt-b and the ND2) by DnaSP version 5.10 (Librado & Rozas 2009). Because each gene region was evaluated independently in different studies (the cyt-b gene in Pons et al. 2010, the ND2 gene in Perktas et al. 2011), two different matrixes of geographic distances (km) between all pairs of populations from each study were estimated using Geographic Distance Matrix Generator (version 1.2.3). A Mantel test with 10000 random permutations was performed between matrices of genetic (F_st) and geographic (log) distances (Slatkin 1993, Rousset 1997).

Coalescent simulations - approximate Bayesian computation implemented in DIYABC version 2.0 beta (Cornuet et al. 2008) was used to infer most likely coalescent scenarios on the demographic history of the Eurasian green woodpecker. In approximate Bayesian computation methods, different population scenarios were modeled by simulating datasets fitting each defined scenario. Then, statistical tests were used to define which scenario best fits the observed dataset. For this study, two different scenarios were tested using the ND2 gene, which was geographically the most complete data set in the species distribution range. First scenario simulated decreasing ancestral population effective size (N_ef) before the Last Glacial Maximum (22000 years ago) to the refuge area, then after the Last Glacial Maximum N_ef started to increase. On the other hand, second scenario simulated stable effective population size through the time. So, no refuge area was taken into account in determining the second scenario. In both scenarios, the prior distributions for all historical parameters and mutation rate were set as uniform. Mutation rate was set between 5.00E-09 and 2.00E-08, which corresponded to 1 % and 4 % sequence divergence per Myr. The substitution model was Hasegawa-Kishino-Yano. Six different summary statistics were taken into account: number of haplotypes, number of segregating sites, mean of pairwise differences, variance of pairwise differences, Tajima's D, and private segregating sites. Two million simulated datasets were generated (1 million for each scenario). Then, posterior probability of each scenario was estimated using logistic regression on the simulated datasets.

Ecological Niche Modelling

Modelling - species occurrence data were compiled from the following sources: GBIF ( www.gbif.org), KuşBank ( www.kusbank.org), and ORNIS ( www.ornisnet.org). The occurrence records showed more intense sampling (e.g. France, Sweden) in the northwest of and much less (e.g. Poland, Romania) or no sampling (e.g. Russia, Ukraine) in the east of the Western Palearctic. To correct for sampling bias and to reduce spatial autocorrelation (Boria et al. 2014, Fourcade et al. 2014), the occurrence records were removed by reducing multiple records to a single record within 50 km distance, resulting 454 records used for ecological niche modelling. This reduced the spatial aggregation of the occurrence records, but could not correct the lack of occurrence records in the east of the Western Palearctic (see Fig. 5). Thus, a Gaussian kernel density of spatially filtered occurrence records was also created by a sampling bias distance of two decimal degrees, producing a bias grid used for ecological niche modelling (Fig. S1; see below; Elith et al. 2010, Fourcade et al. 2014).

Present (1950–2000) and Last Glacial Maximum (22000 years ago) bioclimatic data were downloaded from the WorldClim database (Hijmans et al. 2005,  www.worldclim.org) at a resolution of 2.5 arc-min (4.63 km at the equator). Last Glacial Maximum bioclimatic data are based on three general circulation model (GCM) simulations: CCSM4, MIROC-ESM, and MPI-ESM-P (for detailed information for downscaled paleoclimate data, see CMIP5 - Coupled Model Intercomparison Project Phase 5; cmip-pcmdi.llnl. gov/cmip5). Bioclimatic data include 19 bioclimatic variables derived from monthly temperature and precipitation values (for detailed descriptions of these bioclimatic variables, see  www.worldclim.org/bioclim). Since the geographic distribution of the Eurasian green woodpecker covers a large range in the Western Palearctic, these bioclimatic variables were masked to include only - 12° to 62° E and 29° to 72° N (representing our hypothesized M area, as in Soberón & Peterson 2005). To reduce multicollinearity among bioclimatic variables (Dormann et al. 2013), some highly intercorrelated (r > 0.9 or < -0.9) variables were removed, leaving Bio 2–4, 8–11, and 15–19 (see  www.worldclim.org/bioclim) as input variables used for ecological niche modelling. Correlation structure among these bioclimatic variables appeared to be stable and consistent across time periods (Mantel test: r ≥ 0.857, P < 0.0001), suggesting that the model may be transferable to Last Glacial Maximum bioclimatic conditions (Jiménez-Valverde et al. 2009).

The maximum entropy machine learning algorithm in the software MaxEnt version 3.3.3k (Phillips et al. 2006, Phillips & Dudík 2008, Elith et al. 2011;  www.cs.princeton.edu/∼schapire/maxent) was used to predict the geographic distribution of the Eurasian green woodpecker under present (1950–2000) and reconstructed past (Last Glacial Maximum, 22000 years ago) bioclimatic conditions. MaxEnt, which is among the most effective methods of ecological niche modelling (Elith et al. 2006), was run with the default settings, except for the cases noted below. To optimize the model performance (Shcheglovitova & Anderson 2013), 40 models were tested, using all combinations of five sets of feature types (1. linear, 2. linear and quadratic, 3. hinge, 4. linear, quadratic, and hinge, and 5. linear, quadratic, hinge, product, and threshold, as in Brown 2014) and eight values of regularization multiplier (0.05, 0.25, 0.50, 0.75, 1, 2, 5, and 10, as in Richmond et al. 2010, Gür 2013). The best model was chosen in the following order of preference: (1) model with the lowest omission rate, (2) model with the highest AUC (see below), and (3) model with the simplest set of feature types (Brown 2014) and included all feature types (linear, quadratic, hinge, product, and threshold) and a regularization multiplier of 0.05. To test whether the model performance is robust to variation in the selection of the occurrence records for training and testing, a ten-fold cross-validation was performed in which a different 90 % of the occurrence records were used to train the model and 10 % were used to test it for each of 10 runs (Gür 2013). To up-weight the occurrence records with fewer neighbors in the geographic landscape, the bias grid (Fig. S1; see above) was implemented in the “bias file” option (Elith et al. 2010). To remove the degree of clamping (negative effects caused by novel bioclimatic conditions) from the predictions under Last Glacial Maximum bioclimatic conditions, fade by clamping was performed (Habel et al. 2011), because the MESS (multivariate environmental similarity surface, Elith et al. 2010) maps showed that the model was required to extrapolate into novel Last Glacial Maximum bioclimatic conditions. However, novel bioclimatic conditions encountered during projecting were observed mainly in areas covered by ice during the Last Glacial Maximum (Fig. S2). Thus, extrapolation should not be a critical issue.

Fig. 2.

Isolation-by-distance of populations of the Eurasian green woodpecker (Picus viridis) based on two loci (the cyt-b and the ND2 genes).

Fig. 3.

The extended Bayesian skyline plot for the Eurasian green woodpecker (Picus viridis). The bar under the coalescent time axis shows the period of the Last Glacial Maximum based on different mutation rates used in this study.

Fig. 4.

Alternative biogeographic scenarios for the Eurasian green woodpecker (Picus viridis). t indicates the time scale and N indicates the effective population size. The prior distributions for the time scale and the effective population size are given in the brackets.

Fig. 5.

(A) Present and (B, C, D, E) Last Glacial Maximum predictions of the geographic distribution of the Eurasian green woodpecker (Picus viridis). The visible area in maps is - 12° to 62° E and 29° to 72° N, representing our hypothesized M area. The red circles indicate the occurrence records and warmer colors (more red) indicate areas of relatively high suitability. The white shading shows the ice extent during the LGM (Svendsen et al. 2004). Note that Last Glacial Maximum coastlines differ from present coastlines, because the sea level was lower in the Last Glacial Maximum than in the present.

Fig. 6.

Creation of present and Last Glacial Maximum dispersal networks for the Eurasian green woodpecker (Picus viridis). (A) Present and (B) Last Glacial Maximum friction layers. (C) Present and (D) Last Glacial Maximum dispersal networks. The red circles indicate the occurrence records and warmer colors (more blue in A and B, more black in C and D) indicate areas of low dispersal cost. The white shading shows the ice extent during the LGM (Svendsen et al. 2004). Note that Last Glacial Maximum coastlines differ from present coastlines, because the sea level was lower in the Last Glacial Maximum than in the present.

The area under the receiver operating characteristic (ROC) curve (AUC) was used to evaluate the model performance (Fielding & Bell 1997). An AUC > 0.5 indicates that the model performs better than a random prediction, i.e. AUC ≥ 0.9 = very good; 0.9 > AUC ≥ 0.8 = good; and AUC < 0.8 = poor (Gassó et al. 2012). Landscape connectivity - first, the occurrence records used for ecological niche modelling were further removed by reducing multiple records to a single record within 200 km distance, resulting 73 records (see Fig. 6), to reduce time necessary for analysis and three predictions (i.e. logistic values) under Last Glacial Maximum bioclimatic conditions (i.e. three GCM simulations) were averaged to generate a final summary prediction. Then, to create dispersal networks for once again spatially filtered occurrence records in the present and the Last Glacial Maximum, the predictions under present and Last Glacial Maximum bioclimatic conditions were inverted for use as friction layers, i.e. areas of high suitability were converted to areas of low dispersal cost. The least-cost corridors (LCC) were calculated in three classes (the lowest 1 % LCPs, lowest 2 % LCPs, and lowest 5 % LCPs) by percentage of “the least-cost path (LCP) values” method. The weighed values were applied to each class of LCCs (values of 5, 2, and 1 to low, mid, and high classes, respectively), because all LCCs were summed to create dispersal networks (Chan et al. 2011, Brown 2014).

Except that the predictions under present and Last Glacial Maximum bioclimatic conditions were obtained by the software MaxEnt, all analyses were conducted using the software SDMtoolbox version 1.1b (Brown 2014,  www.sdmtoolbox.org) implemented in the software ArcGIS version 10.1 ( www.arcgis.com).

Phylogeography

In this study, we focused on previously published multilocus mtDNA dataset of the Eurasian green woodpecker. In total, we detected six unique haplotypes for the cyt-b gene and 24 unique haplotypes for the ND2 gene. Seventeen variable sites occurred in the cyt-b gene, six of which were parsimony informative. Eight of transitions were A/G, eight were C/T, and only one transversion was detected, which was A/T. The ND2 gene contained more variable sites than the cyt-b gene. 35 variable sites were found, 11 of which were parsimony informative. 13 of transitions were A/G, 21 were C/T, and only one transversion was detected, which was C/A.

Table 1.

Relative posterior probabilities with 95 % credibility intervals for each scenario.

Haplotypes were not geographically structured in neither the cyt-b gene nor the ND2 gene (Fig. 1). The ND2 gene had higher genetic diversity than the cyt-b gene. Estimated mutation rate of the ND2 gene was almost 4.5 times faster than the cyt-b gene (for details see discussion). According to the variation within the ND2 gene, genetic diversity of southern populations was higher than northern populations. Interestingly, however, a population from United Kingdom had relatively higher nucleotide diversity (see Table 1 in Perktas et al. 2011). Four populations in Europe (Austria, Balkans, Italy and United Kingdom) contained the most common haplotype of the ND2 gene. However, northern Europe (Sweden), Caucasus, and Turkey did not contain this haplotype. Geographic origin of the most common haplotype of the ND2 gene was south and almost mid-latitudes of Europe.

Haplotype network based on the cyt-b gene had low frequency haplotypes differed from the most common haplotype by one, two or three base pairs. A similar network pattern was found for the ND2 gene. However, the only difference between both haplotype networks was the number of divergent haplotypes in the ND2 gene. These haplotypes differed from the most common haplotype by up to six base pairs, and located in Italy, Turkey, Caucasus, and Russia; they made genetic diversity higher in southern part of distribution range of the Eurasian green woodpecker in Europe, Anatolia and Caucasus (Fig. 1).

Table 2.

Properties of the approximate posterior distribution of parameters under the scenario 1.

Despite mostly low-level genetic variation between populations of the Eurasian green woodpecker in Europe, Anatolia and Caucasus (maximum divergence was 0.007 for the cyt-b gene and 0.0077 for the ND2 gene), genetic distances were positively correlated with geographic distances for both mtDNA genes (Fig. 2), which indicate an evidence of isolation-bydistance (the cyt-b gene, r = 0.568, P < 0.0001; the ND2 gene, r = 0.374, P = 0.001). All of these results led us to justify combining all sequences of two different mtDNA genes from all populations while assessing historical demography of the Eurasian green woodpecker. Extended Bayesian skyline plot showed a good resolution of effective population size fluctuations over the history of the Eurasian green woodpecker. Therefore, the result of the extended Bayesian skyline plot detected a clear population expansion after the Last Glacial Maximum (Fig. 3). Finally, statistical phylogeography (i.e. coalescent simulations) supported the result of the extended Bayesian skyline plot and suggested that the best explanation of population expansion was scenario 1 (Fig. 4), which indicated a population expansion after the Last Glacial Maximum. The posterior probability was 0.9984 for this scenario (Table 1). Based on the approximate posterior distribution of parameters under scenario 1, population expansion began about 14500 years before present and the effective population size increased from 215000 to 375000 (Table 2).

Ecological Niche Modelling

The model performed better than a random prediction (AUC for training data (mean ± standard deviation (SD), range) = 0.928 ± 0.001, 0.926–0.930; AUC for test data (mean ± SD, range) = 0.842 ± 0.013, 0.822–0.861, based on the 10-fold cross-validation runs). The small SDs for the mean AUCs suggested that the model performance was robust to variation in the selection of the occurrence records for training and testing.

Under present bioclimatic conditions, areas of relatively high suitability were predicted across most of the known distribution of the Eurasian green woodpecker. However, the model underpredicted the distribution in the east of the Western Palearctic (e.g. Russia, Ukraine) and overpredicted the distribution in Iberia (Fig. 5). The former is most likely to be an effect of sampling bias (e.g. no sampling in the east of the Western Palearctic), although correction for sampling bias was applied (Fig. 5, Fig. S1), and/or of that the distribution is poorly described in that region. The latter is most likely to be an effect of biotic interaction, because a sister species (Iberian green woodpecker) is distributed in Iberia, i.e. two species show parapatric distributions suggestive of competitive interaction.

Under the Last Glacial Maximum bioclimatic conditions, areas of relatively high suitability were predicted across southern Europe (including France), Anatolia, and the Caucasus/Caspian region. Three predictions under three GCM simulations were largely similar in this respect (Fig. 5).

During present bioclimatic conditions, a relatively high degree of population connectivity was observed across most of the Western Palearctic, with the higher degree of connectivity from the northern Anatolia through the Balkans and central Europe to England and Scandinavia. During the Last Glacial Maximum bioclimatic conditions, a relatively high degree of population connectivity was observed across southern Europe (including France), Anatolia, and the Caucasus/Caspian region, with the higher degree of connectivity from the northwestern Anatolia through the Balkans to Italy (Fig. 6).

All these results suggested that the Eurasian green woodpecker survived the Last Glacial Maximum in a single large glacial refugium, encompassing southern Europe (including France), Anatolia, and the Caucasus/Caspian region, and, afterwards, substantially expanded its range from this LGM refugium to its present range (Fig. 5, 6).