Skip to main content

Advertisement

Genetic, phenotypic and ecological differentiation suggests incipient speciation in two Charadrius plovers along the Chinese coast

Article metrics

Abstract

Background

Speciation with gene flow is an alternative to the nascence of new taxa in strict allopatric separation. Indeed, many taxa have parapatric distributions at present. It is often unclear if these are secondary contacts, e.g. caused by past glaciation cycles or the manifestation of speciation with gene flow, which hampers our understanding of how different forces drive diversification. Here we studied genetic, phenotypic and ecological aspects of divergence in a pair of incipient shorebird species, the Kentish (Charadrius alexandrinus) and the White-faced Plovers (C. dealbatus), shorebirds with parapatric breeding ranges along the Chinese coast. We assessed divergence based on molecular markers with different modes of inheritance and quantified phenotypic and ecological divergence in aspects of morphometric, dietary and climatic niches.

Results

Our integrative analyses revealed small to moderate levels of genetic and phenotypic distinctiveness with symmetric gene flow across the contact area at the Chinese coast. The two species diverged approximately half a million years ago in dynamic isolation with secondary contact occurring due to cycling sea level changes between the Eastern and Southern China Sea in the mid-late Pleistocene. We found evidence of character displacement and ecological niche differentiation between the two species, invoking the role of selection in facilitating divergence despite gene flow.

Conclusion

These findings imply that ecology can indeed counter gene flow through divergent selection and thus contributes to incipient speciation in these plovers. Furthermore, our study highlights the importance of using integrative datasets to reveal the evolutionary history and assist the inference of mechanisms of speciation.

Background

Understanding how strongly evolutionary processes (i.e. selection, gene flow and genetic drift) shape the divergence of closely related species has been a long-standing interest in evolutionary biology [1,2,3,4,5]. Allopatric speciation is conventionally considered to be the prevalent mode of speciation in which physical barriers completely restrict gene flow between two populations, facilitating the initiation of divergence through genetic drift or selection [6,7,8]. If populations remain isolated for a long enough period of time after divergence has been established, genome-wide differentiation is expected to be accumulated even with fixation of novel mutations. Such divergence can be maintained and eventually promote reproductive isolation due to genetic incompatibilities even in the presence of gene flow after secondary contact [9, 10].

An increasing number of studies have shown that divergence can arise and be maintained due to selection imposed by heterogeneous environments despite the homogenizing effect of gene flow [11, 12]. Under this scenario, divergent selection or ecologically-mediated sexual selection operating on certain traits can lead to reproductive isolation. Incompatibilities in a few “speciation genes [13, 14]” may be enough to constrain the homogenizing effect caused by gene flow, even at an early stage of speciation [15,16,17]. For instance, different host plant preferences result in divergence in digestive, physiological and morphological traits within sympatric ecotypes of Timema walking-stick insects [18]. In such a system, adaptive loci which are responsible for color pattern were confirmed to overcome local gene flow [19]. Locally-adapted phenotypes in body size are related to swimming and foraging ecology in guppies (Poecilia reticulata), which can maintain population diversification and fixed divergence in genomic regions [20] in the face of extensive gene flow [15]. Though evidence of speciation with gene flow is accumulating [11, 12, 21, 22], to what extent a balance between selection and gene flow can drive divergence remains largely unknown.

Confirming whether genetic differentiation between two taxa is due to strict allopatry or speciation-with-gene-flow remains challenging [23]. Towards this end, ideal study systems include populations found across a variety of environments, and which have a high propensity for introgressive hybridization such as in waterfowl (e.g. [24]) and shorebirds [25, 26]. Specifically, we expect differences in environments to confer differential selective regimes across the populations that would counteract the amalgamating forces of gene flow; and thus, providing a means to understand how selection, gene flow, and genetic drift effect the genomes of these species.

Shorebirds are a group of migratory species with remarkable movement ability. Seasonal migration may increase the probability of individual dispersal between populations [27, 28] and consequently drive frequent gene flow between geographically distant populations [29, 30]. Direct evidence, obtained using tracking approaches, confirmed that extremely long-distance gene flow can be mediated through individual movements among remote breeding colonies of pectoral sandpipers (Calidris melanotos) [31]. Range-wide phylogeographic studies also revealed extensive gene flow in several migratory shorebird species (e.g. [32,33,34,35,36]). This may result in weak genetic structure across a species and consequently prevent population-level divergence [25, 26].

In this study, we test the role of ecology and gene flow on the divergence of two shorebird species, the Kentish Plover (Charadrius alexandrinus) and the White-faced Plover (C. dealbatus; Fig. 1). C. alexandrinus breeds in coastal areas and inland lakes in Europe, Asia and North Africa [37]. A previous study uncovered low genetic differentiation of C. alexandrinus across the Eurasian continent, and also between continental and island populations in East Asia [32]. C. dealbatus was formerly described as a distinctive species [38, 39]. Breeding records of C. dealbatus were reported along the coast of China from Fujian to Hainan Island, and in south-central Vietnam [40], yet its geographic range and boundaries with C. alexandrinus are uncertain [41]. Previous studies have found subtle differences in morphometric, plumage and behavioral traits between the two taxa (Fig. 1d), leading to the erection of C. dealbatus as a full species [41, 42]. But no evidence of genetic differentiation suggested by the first genetic investigation, which concluded that these species may only differ in a few genomic regions [42]. Therefore, this taxonomic treatment has been accepted by some world bird checklists [43], but is yet to be accepted by others [44].

Fig. 1
figure1

Sampling localities, morphology and trophic level differentiation of the Kentish Plover Charadrius alexandrinus (blue) and the White-faced Plover C. dealbatus (yellow). a Samples were from an inland site (1, Qinghai Lake), 16 Chinese mainland coastal localities (2, Tangshan; 3, Cangzhou; 4, Weifang; 5, Lianyungang; 6, Nantong; 7, Ningbo; 8, Zhoushan; 9, Wenzhou; 12, Fuzhou; 13, Xiamen; 14, Jinmen; 15, Shanwei; 16, Yangjiang; 17, Zhanjiang; 18, Beihai and 19, Dongfang), and two localities on Taiwan Island (10, Xinbei and 11, Zhanghua). Map source: National Geomatics Center of China (http://www.ngcc.cn). Differences observed between two Charadrius species in several characters. C. alexandrinus has on average shorter bill length (except for Taiwan populations) (b), wing length (c) and lower body mass (except for Taiwan populations) (d) than C. dealbatus. Localities with less than five measured adults were excluded in these analyses. e C. alexandrinus has darker plumage than C. dealbatus in general and males of the latter species have less black tinges on their face and neck during breeding season. Shaochong Peng drew the artwork of the two species of plovers and these images used with permission. f Plots of the first two principal components and their associated variance showed the subtle morphometric differences between the two species. g Pairwise morphological difference QST plotted against geographical distance between coastal populations of C. alexandrinus and C. dealbatus. QST values within species are marked by blue or yellow dots. QST between species (grey circles) declined as the distance increased (regression showed as the dashed line)

Here, we provide comprehensive data from multiple sources to explore divergence patterns between these two species of plovers. We carried out extensive sampling of breeding populations covering the potential area of contact along the Chinese coast. We then obtained genetic polymorphism data from multiple markers with various mutation rates, i.e. mitochondrial DNA, exons and autosomal microsatellites. Using these new datasets, it is possible to estimate the intensity and direction of gene flow [45, 46] as well as other important demographic parameters such as the effective population size (Ne) and the timing of divergence [47] between C. alexandrinus and C. dealbatus. Because it is difficult to directly quantify divergent selection, it is necessary to make inference from the comparison of traits, which are potential targets of selection in order to offer hints on selective mechanisms in maintaining divergence [18, 48, 49]. By collecting data on morphology, diet, and environmental niche, we attempted to characterize geographic variation in genetic and phenotypic traits, estimate demographic history and intensity of gene flow, and infer the role of multiple factors, i.e. geographical premating isolation and gene flow during divergence, between these two taxa. Finally, because we applied multiple lines of evidence to demonstrate divergence patterns, we attempt to re-evaluate species delimitation between the two plovers. Taken together, our investigation provides new insights into the evolutionary history of two incipient bird species.

Results

Morphometric differentiation

Both plover species showed subtle but significant differences for most morphological traits (ANCOVA, p < 0.001, Fig. 1) even though there was considerable overlap in morphology between the species at the level of the individual (Fig. 1f). On average, C. dealbatus had a longer bill (17.91 ± 0.94 mm vs. 16.69 ± 1.06 mm in C. alexandrinus, p < 0.001, Fig. 1b), longer wings (118.58 ± 2.90 mm vs. 115.59 ± 3.02 mm in C. alexandrinus, p < 0.001, Fig. 1c), and larger body mass (48.99 ± 3.26 g vs. 47.56 ± 3.73 g in C. alexandrinus, p = 0.001, Fig. 1d, Additional file 1: Table S5) than those of C. alexandrinus. There was no difference in tarsus length between the two species (p = 0.962). Individuals of C. alexandrinus from Taiwan Island were heavier than continental populations of C. alexandrinus and C. dealbatus (both p < 0.001). QST between coastal populations of C. alexandrinus and C. dealbatus was negatively correlated with geographic distance (R2 = 0.293, p < 0.001, Fig. 1g).

Genetic diversity and population structure

We sequenced 357 individuals at all three mtDNA loci (224 C. alexandrinus and 133 C. dealbatus, GenBank Accession No. MK830738-MK830815). For each site, the sample size ranged from 11 to 30 individuals. For each individual, we obtained in total 1729 base pairs (bp) of mtDNA sequence, including 846 bp ATPase6/8, 505 bp D-loop and 378 bp ND3. C. alexandrinus showed higher intraspecific genetic diversity than C. dealbatus (Table 1). Haplotype networks show that at all loci, most individuals were sorted into two major haplogroups, corresponding to the two species of plovers. One non-synonymous substitution separated the two haplogroups at both ATPase6/8 and ND3 loci (Fig. 2a). Moreover, a subset of samples containing 20 individuals of each species were sequenced at 16 loci (range 440–902 bp for each locus; Table 1 and Additional file 1: Table S1, GenBank Accession No. MK830816-MK830957) for a total of 11,209 bp of nuclear DNA sequence. The haplotype networks from autosomal and Z-linked loci did not show strong patterns of lineage sorting like the mtDNA. The most common haplotypes were shared by both species of plover (Fig. 2b and Additional file 1: Figure S1). Moreover, both species showed signs of recent demographic expansion as detected by significant Tajima’s D values (Table 1).

Table 1 Sampling localities and genetic polymorphism of C. alexandrinus and C. dealbatus. The number of individuals (n) analyzed for mtDNA, autosomal microsatellites and nuclear exonic loci (nuDNA) are given. Site number corresponds to the numbers in Fig. 1. Estimates of h, number of haplotypes; Hd, haplotype diversity; π, nucleotide diversity; Tajima’s D value, Ho, observed heterozygosity; He, expected heterozygosity were calculated for each locality and species
Fig. 2
figure2

Haplotype networks based on mitochondrial and nuclear DNA sequences, and population genetic structures based on microsatellite loci. C. alexandrinus is marked in blue and C. dealbatus in yellow. a Haplotype networks of three mitochondrial loci. For ATPase6/8 and ND3, over 95% of the individuals are sorted into two major haplogroups. b Examples of haplotype networks based on four exonic nuclear loci (NFIL3 and S1PR3 are Z-linked loci) of 55 to 80 individuals. c Genetic clustering inferred from microsatellite genotypes using Geneland. Blue and yellow bars show the assignment probability of a location for alternative genetic clusters. Location numbers are consistent with Fig. 1. d Genetic clustering inferred from microsatellites using STRUCTURE. Map source: National Geomatics Center of China (http://www.ngcc.cn)

Overall genetic differentiation was significant and high in mtDNA data (ΦST = 0.506, p < 0.001) and low at microsatellite loci (FST = 0.036, p < 0.001). For nuclear sequence data, genetic differentiation between species was also significant at autosomal loci (ΦST = 0.100, p < 0.001) and particularly high at the Z-linked loci NFIL3 (Φ ST = 0.726, p < 0.001) and S1PR3 (ΦST = 0.309, p < 0.001). In the AMOVA analysis, we observed the largest difference between groups when the data were partitioned by species (i.e. C. alexandrinus and C. dealbatus), with this grouping explaining 49.7% of the variance in mtDNA and 2.4% in the microsatellites (p < 0.001; Table 2). These were significantly higher than the values of genetic variation observed when we partitioned samples as coastal vs. island populations or coastal population vs. Qinghai vs. Taiwan populations. Within C. alexandrinus, minor genetic differentiation was found between the inland population (Qinghai Lake) and coastal populations (mtDNA ΦSC = 0.042, p = 0.022; mst FSC = 0.021, p < 0.001). C. alexandrinus populations on Taiwan Island also shared haplotypes with other coastal populations (mtDNA ΦSC = 0.021, p = 0.146), but were significantly differentiated in microsatellites (microsatellite FSC = 0.016, p < 0.001).

Table 2 Hierarchical analyses of molecular variance (AMOVA) based on concatenated mtDNA data and 13 microsatellite loci for C. alexandrinus and C. dealbatus. Samples were partitioned in three groupings: 2 groups (C. alexandrinus and C. dealbatus); 3 groups (C. alexandrinus continental populations and Taiwan island populations, C. dealbatus); 4 groups (C. alexandrinus inland (Qinghai Lake), C. alexandrinus coastal, C. alexandrinus Taiwan island, and C. dealbatus). Va, genetic variation among groups; Vb, variation among populations within groups; Vc, variation within populations. Between-group genetic differentiation was highest when populations were partitioned into two groups (C. alexandrinus and C. dealbatus)

For microsatellite loci, 18 out of 22 markers were successfully genotyped. However, four markers (Calex-04, 08, 19, C204) showed a large proportion of missing data, while another four loci (Calex-11, 24, 26, 43) showed an estimated frequency of null alleles over 10%. Moreover one locus (Calex-35) showed H-W disequilibrium (Additional file 1: Table S2). Consequently, the aforementioned nine loci were removed from the dataset, so the final microsatellite dataset used for further analysis consisted of the genotypes of 402 individuals (271 C. alexandrinus and 131 C. dealbatus) at 13 different loci.

The results of the STRUCTURE analysis clearly showed two genetic clusters representing each species (Fig. 2d) and no obvious gradual transition along the coastline that would be expected from a hybrid zone. The average Delta K value when K = 2 was much higher than that of other options (Additional file 1: Figure S2). Using georeferenced data in GENELAND, our results corroborated the genetic clustering patterns inferred from STRUCTURE. The posterior probability was 0.70 when K = 2, which contrasted with a value of 0.25 when K = 3. We visualized the GENELAND results on a map with the probability distribution of posterior mode of class membership, which further supports the separation between C. alexandrinus and C. dealbatus along the China coast (Fig. 2c). The divide between these two species was located between Wenzhou and Fuzhou, according to the GENELAND results, but it is unclear if these two species come into contact or form a putative hybrid zone in this region. Again, the individuals from the sites in Taiwan were assigned to the cluster of C. alexandrinus (Fig. 2c).

Inferred demographic history

Isolation with migration analyses suggested that C. alexandrinus and C. dealbatus diverged 0.56 (0.41–5.19) million years ago. Estimated migration rates in both directions were significant (p < 0.001) with slightly higher gene flow from C. alexandrinus into C. dealbatus (2.69 migrants per generation) than vice versa (ca. two migrants per generation). The estimated effective population size of C. alexandrinus (Ne ≈ 1.59–4.44 million) was about 8–10 times higher than for C. dealbatus (Ne ≈ 0.15–0.52 million).

The estimation of recent gene flow between the two species, carried out with BayesAss using microsatellites, also suggested bidirectional gene flow. Gene flow from C. dealbatus to C. alexandrinus was slightly higher (0.013, p = 0.028) than in the other direction (0.010, p = 0.028). In C. alexandrinus populations, one likely migrant from C. dealbatus and two hybrid F1 individuals were identified (probability higher than 0.5). With the same threshold value, two migrants and one hybrid F1 individual were identified in C. dealbatus populations.

The detection of hybrids

Based on the STRUCTURE results and simulations, the optimal threshold for distinguishing hybrids from non-hybrids was q = 0.836. Based on this threshold, 81.9% of the individuals (204 out of 246) collected from sites at the northern Chinese coastline, Qinghai Lake and Taiwan Island were assigned to C. alexandrinus (Fig. 2c, d). Individuals with intermediate q-values, which were possibly hybrids, were found at most northern Chinese sampling sites (Fig. 2d). For the southern coastline in China, 145 out of 156 Individuals (92.9%) belonged to C. dealbatus with a probability larger than 0.836 (Fig. 2c, d). Only one individual of each species was assigned to the a migrant from the other species with high probability (Fig. 2d).

Ecological niche profiles

Our ecological niche models effectively captured the current distribution of both C. alexandrinus and C. dealbatus (Fig. 3) with a high discrimination capacity (AUC values > 0.88 for training and test data). Jackknife tests on variable importance for C. alexandrinus revealed that isothermality, precipitation seasonality and mean temperature of the warmest quarter were the three highest ranked variables when used in isolation. For C. dealbatus, mean diurnal range and annual precipitation were the most important variables. The simulation of the three periods, i.e. Last Interglacial (LIG, 120-100Ka), Last Glacial Maximum (LGM, 21Ka, MIROC model) and current times, respectively, showed range shifts in both species. In particular, these results suggest that the ranges of both species shrank during the LIG in the Chinese coastal area accompanied by an increase in climatic suitability in the inland region for C. alexandrinus. In contrast, suitable habitats expanded for both species during the LGM in the coastal area, the East China Sea and the northern part of the South China Sea (sea between Hainan and Taiwan), probably due to a fall in sea level.

Fig. 3
figure3

Predicted environmental suitability for C. alexandrinus (left) and C. dealbatus (right), via ecological niche modeling (ENM). ENM results are shown for the Last Interglacial (LIG, 120-100Ka), Last Glacial Maximum (LGM, 21Ka, MIROC model) and current time periods, respectively. Maps were drawn based on the projected distributions using ArcGIS 9.3 (ESRI, Redland, CA. URL http://www.esri.com/)

The ordination approach, using PCA-env, suggested that the overlap of the current climatic niches of the two species is relatively low (Additional file 1: Figure S3a, b). The two species can be separated based on the first two PCs with an accumulative 81.5% of the total variance explained (Additional file 1: Figure S3c). The niche equivalency test rejected the null hypothesis that the species pair is distributed in identical environmental space (p = 0.019; Additional file 1: Figure S3d).

Stable-isotope profiles

Our stable-isotope analysis showed that C. dealbatus exhibited significantly higher δ15N values than C. alexandrinus (p < 0.001, Fig. 4a, Additional file 1: Figure S4). In contrast, we did not find significant differences in δ13C between the two species (p = 0.161). Further isotope space overlap analysis showed that C. alexandrinus individuals had a high probability of being found within the isotopic niche space of C. dealbatus (95.3%), while C. dealbatus individuals showed a relatively low probability of being found in the isotopic niche space of C. alexandrinus (46.5%, Fig. 4b). Moreover, we found that C. alexandrinus showed higher variability in δ15N profile across breeding populations than C. dealbatus.

Fig. 4
figure4

a Raw data (bottom-left), density distribution (top-left and bottom-right), and stable isotopes related to dietary niche regions (top right) of δ13C and δ15N profiles of C. alexandrinus and C. dealbatus feathers generated by “nicheROVER” (Swanson et al. 2015). The niche region of C. alexandrinus (blue) contained that of C. dealbatus (yellow) and was the broader of the two regions. b Distribution of the posterior probability that an individual from one species is found within the niche region of the other species. Vertical lines for mean and 95% credible intervals are included in the histogram of each overlap metric

Discussion

Our genetic data show that C. alexandrinus and C. dealbatus have diverged to a level of advanced allele sorting between the two species, particularly in mtDNA and Z-linked genes with lower effective population sizes (Fig. 2a, b). For autosomal microsatellite data, we also found that genetic differentiation is low at the intraspecific level but substantially high between the species (Fig. 2c, d, Table 2). Our results suggest the two plover species diverged approximately 0.6 mya, and both have large effective population sizes. Though we found no or just a narrow hybrid zone exists in the contact area on the Chinese coast, a considerable level of symmetric gene flow occurred between the two plovers (Table 3). Further, we found significant differences in morphometric traits and ecological characters between the two plovers along the coast (Fig. 1e-g, 3 and 4).

Table 3 Posterior mode, mean and range of 95% highest probability distribution (HPD) of six demographic parameters inferred with IMa2p between C. alexandrinus and C. dealbatus. Divergence times (T) are given in million years ago (Ma). The effective population size (Ne) of C. alexandrinus was about eight times that of C. dealbatus. Migration rate (2NM) into each species was about the same. K represents Kentish Plover C. alexandrinus; W represents White-faced Plover C. dealbatus; A is the most recent common ancestor of two species

Phylogeographic patterns

Two genetic lineages were found among breeding Charadrius plovers along the Chinese coast, the northern lineage corresponding to C. alexandrinus and the southern one to C. dealbatus (Fig. 1a, Table 2). The sharp genetic break between the two lineages lies between Wenzhou and Fuzhou (latitude 26–27 °N), north of the Taiwan Strait (Fig. 2c, d). Furthermore, samples from Taiwan Island belong to C. alexandrinus but the population in Jinmen Island (site 14 in Fig. 1a), off Fujian coast is affiliated with C. dealbatus (Fig. 2c, d). Our IMa analysis estimated that the divergence time between the two plovers was back to approximately 0.6 mya (Table 3, Fig. 5), which falls into Marine Isotope Stage 16 of the mid-late Pleistocene climatic fluctuation periods [50]. For ocean marginal and coastline taxa, this implies that fluctuations in sea level can lead to alternation between population isolation and contact throughout their evolutionary histories. Because Charadrius plovers originated in the Northern hemisphere and then radiated to the Southern hemisphere [51], it is thereby likely that C. dealbatus evolved from a recent common ancestor with C. alexandrinus and could have been diverging along the east China coast since the mid-late Pleistocene. During this period, vicariance occurred because the coastline was separated by a land-bridge raised due to the shallow sills between the East and the South China Sea [52]. As shown in our niche modeling analysis, both species had decreased distribution ranges during the LIG followed by enlarged suitable habitats and overlapping ranges between the East and South China Sea during the Last Glacial Maximum (Fig. 3). While our niche modeling demonstrates a cycle of range dynamics caused by sea-level changes, one should bear in mind that the last million years witnessed multiple glaciation cycles [53, 54]. It is possible that the divergence between the two plovers proceeded through cycles of allopatric stages interspersed by secondary contacts due to sea-level fluctuation.

Fig. 5
figure5

Posterior densities of population demographic parameters estimated using the isolation-with-migration model (IM) implemented in IMa2p. These analyses used the combined dataset of three mitochondrial (1729 bp) and 15 nuclear exonic loci (11,209 bp). K represents the Kentish Plover C. alexandrinus; W represents the white-faced Plover C. dealbatus. a Population divergence times (T) of C. alexandrinus and C. dealbatus. b Population migration rates (2NM) represent the average number of individuals in each group that had previously migrated from the other group. Gene flow in both directions was significant. c Effective population sizes (Ne) of the two species and their most recent common ancestor

To the best of our knowledge, our study shows the first documented phylogeographic break in a bird species along the coastline of China. Phylogeographic patterns have been relatively well characterized in shorebirds breeding at high latitudes [32,33,34]. Panmixia or weak genetic differentiation was usually suggested, most likely driven by extensive gene flow. However, it seems that islands can act as a barrier to natal dispersal for the continental Kentish Plovers [32]. In contrast, population structures in temperate and subtropical shorebird populations are poorly documented, probably due to a low level of species diversity. Interestingly, concordant phylogeographic patterns have been described in several coastal marine taxa, such as plants [55], fishes [56, 57], shellfishes [58], and crustaceans [59]. Though the exact splitting times are not congruent, the observed split line is at approximately 25°N latitude between the East and the South China Sea (reviewed in [52]). This pattern resulted in the hypothesis that historical factors, i.e. sea level fluctuation during the Pleistocene, caused a convergent phylogeographic pattern in multiple coastal marine fauna in the marginal northwestern Pacific Ocean [52, 56]. Thus, our study contributes a vertebrate case to the accumulating literature about this hotspot of species divergence. Apart from the major role of physical barriers, comparative phylogeographic studies also revealed that other abiotic factors, like ocean currents and hydrothermal conditions, as well as the ecological characters of species, i.e. dispersal ability, habitat preference, life-history, and population demography, can also play a role in contributing to the divergence of coastline fauna [52, 57].

Phenotypic and ecological differentiation along a latitudinal gradient

The two plovers, C. alexandrinus and C. dealbatus, are distributed along the Chinese coastline across latitudinal and associated environmental gradients from temperate to tropical zones. Within each species, we found a general trend that northern populations have larger morphometric values when compared with their southern counterparts (Fig. 1). This pattern may be related to Bergmann’s rule which states that the body size of homeothermic animals is larger in colder climates than in warmer ones [60]. Populations of both plover species may benefit from optimal control due to their distribution across such a north-south gradient [61]. A mutually non-exclusive explanation is that the difference in morphological traits, i.e. body mass and bill length, link to potential differences in resource exploitation and life history. Our data show that C. dealbatus has a larger average body mass than C. alexandrinus. Body mass is a comprehensive trait reflecting nutrition assimilation, energy reservation and expenditure [62]. The difference in body mass might be related to the difference in migratory behavior where a lighter body mass in C. alexandrinus is favored by decreased transport cost of fuel storage during migration [63, 64]. In contrast, larger body mass in C. dealbatus, a short-distance migrant or resident, may be beneficial for multiple reproductions within a single breeding season while the lighter C. alexandrinus usually only produces a single clutch (Lin et al. in prep.). The difference in bill length may be driven by the difference in the use of food resources and also foraging strategies [65] but recently, the function of the bill as a temperature regulator has also started to attract scientific attention [66]. An increasing number of studies have demonstrated that birds at higher latitude and in cooler environments have shorter bills (European sparrows, [67], Australian shorebirds, [68]), consistent with Allen’s Rule [69]. Nevertheless, we found that populations of C. alexandrinus in Taiwan had similarly large body mass and bill lengths when compared to C. dealbatus populations (Fig. 1b-d, g), which are significantly larger than mainland conspecific populations. This is probably caused by phenotypic plasticity of the Taiwan population towards sub-tropical ecomophology and a resident life history, similar to the populations of C. dealbatus in the South China Sea coast.

The niche-equivalency test in environmental (E)-space rejects the possibility that the two species are distributed in strictly equivalent environmental space (Additional file 1: Figure S3). Corresponding to this, in geographical (G)-space, C. dealbatus is restricted to breeding sites close to the coast, particularly where there is a warmer tropical climate. In contrast, C. alexandrinus has a wider climatic niche, as represented by a broader climatic zone (Fig. 3). This species can breed not only on temperate coasts [70] but also on inland saline lake shores [71], such as Qinghai Lake. In addition, our isotope analysis revealed that C. alexandrinus covered a wider range of isotope ratios than C. dealbatus, but C. alexandrinus exploited a lower trophic range (δ15N). This indicates that C. dealbatus probably feeds on a higher energy diet than its sister species [72]. However, it is unclear whether such differences result from divergence in diet preference or food resource availability [73]. Taken together, these results suggest that the ecological niches of the two plovers are significantly different in several aspects, and support a role for ecology in constraining the range limits, and perhaps promote reproductive barriers between the two shorebird species.

Molecular signatures of early species divergences in two plover species

After the two nascent plover species split, how was their divergence maintained in the presence of gene flow? Speciation theory predicts that divergence is initiated either by genetic drift or divergent selection [22, 74, 75]. However, genetic drift is unlikely to have been the only force to initiate divergence in the plovers. Given the large effective population sizes (Ne) and recent divergence (Table 3), it is difficult that genetic drift could have led to fixed differences between populations [76, 77]. Gene flow can directly counteract the diverging effects of drift [78, 79] but sufficient selection could help overcome the homogenizing effects of gene flow [1,2,3,4,5]. We found some hints of selection by comparing the level interspecific genetic differentiation among genetic markers. For sequence data, the ratio of ΦST was five times larger at the concatenated mtDNA dataset, and 3~7 times larger at the Z-linked genes than at the autosomal loci, respectively. Under migration-drift equilibrium, one would expect that divergence should be four and 1.33 times larger than the extent of genetic diversity at autosomal markers, respectively [80, 81], and thus the observed elevated extent of genetic divergence for mtDNA, and especially for Z-linked markers cannot be explained by genetic drift alone. For maternally inherited mtDNA, both female philopatry and selective sweeps in the mitochondrial genome could cause the deviation from the neutrality model [24]. However, a high level of female philopatry is less likely due to evidence of female-biased dispersal suggested by field observation and genetic inference in C. alexandrinus [24]. Thus selection on mitochondrial DNA could explain the observed pattern, and this might be further supported by the significantly negative values of Tajima’s D (Table 1) throughout all study populations.

In addition, elevated levels of genetic divergence are often found at Z-linked markers compared with autosomal markers, especially at the nascent stage of speciation, termed as the faster-Z effect [82]. Indeed this pattern has been shown in several empirical studies in birds [reviewed in 125]. Beside neutral explanations, such as the lower Ne (3/4 that of the autosome) and higher mutation rate than autosomes, the Z chromosome has been demonstrated to contribute to reproductive isolation by playing a large role in promoting both prezygotic and postzygotic isolation, because some regions in the Z chromosome can be responsible for encoding sex-related or mating preference traits [83]. Hence divergent selection or ecologically-mediated sexual selection can reduce genetic introgression on the Z chromosome and play a disproportional role in the divergence in birds, but the current data do not have enough power to resolve this link in the plover case.

The loci analyzed here represent just a small proportion of the genome and further investigations using whole genome data hold great potential for providing insights on of the genetic basis of divergence in the two plovers. Particularly, it would be important to identify specific genomic regions (“genomic islands of speciation”) characterized by a high level of divergence and reduced level of gene flow [2,3,4,5]. Data from several organisms suggest that incipient speciation can be maintained with divergence at a small number of genomic regions [11, 12, 84,85,86,87]. Because several morphological and ecological differences between the two plovers have been characterized in the present study, it is worth to understand whether genomic islands of high differentiation included relevant candidate genes that are involved in reproductive isolation.

The geographical boundary between the two plovers was ambiguously defined in previous studies [41, 42]. Our results show that the discontinuity in genetic structuring (Fig. 2c) and morphometric values (Fig. 1b-d) between the two species is situated at the coastline between Wenzhou and Fuzhou, indicating a contact zone (Fig. 2c). For incipient species, individuals may show a clinal pattern of allele frequencies and morphology at their contact area and form a hybrid zone [53, 88, 89]. Our STRUCTURE results revealed no obvious signs of a hybrid zone (Fig. 2d), and rather indicated sporadic migrants in the respective range of the other species. NewHybrid 1.1 [90] was also used to identify potential hybrid individuals and the result was highly concordant with STRUCTURE. However, NewHybrids failed to identify simulated hybrid individuals and recognized them as migrants (Additional file 1: Figure S5), probably due to incomplete divergence between two species. With the possibility that sampling was not fine-grained enough, the potential hybrid zone should be narrower than the 200 km distance. Another possibility for an apparent lack of a hybrid zone is assortative mating between the two plovers at the contact area [91]. The latter explanation is very possible because populations of the two species close to this region were much more differentiated in morphological traits than the ones farther apart (Fig. 1b-d) and showed ecological segregation (see discussion above). In this case, reproductive [92, 93] and ecological [94, 95] character displacement are called as mechanism to contribute reproductive isolation between the two species. Even with the aforementioned approach to disentangle the heterogeneous landscape of genomes, genomic data alone might not reveal the form of reproductive isolation between the two plovers [23]. It is essential to carry out detailed experiments to test assortative mating in the contact zone between C. alexandrinus and C. dealbatus.

Taxonomic implications

Kennerley et al. 2008 [41] recommended that the tropical breeding population, previously defined as the subspecies dealbatus, warranted species status, based on differentiation in morphology, behavior and distribution from C. alexandrinus [41]. Uncertainty about the taxonomic status arose because the first genetic evaluation found no evidence of genetic differentiation between the two plovers [42]. At odds with this, the present dataset is comprised of a larger suite of genetic markers and well- characterized ecological traits (Fig. 1), obtained by systematic sampling along the Chinese coast, while Rheindt et al. 2011 [42] analyzed samples collected outside the breeding season. Using the same mtDNA markers we found interspecies divergence not detected in Rheindt et al. [42]. Such underestimated genetic differentiation may be due to the erroneous species assignment of non-breeding birds, because the distinguishing features in non-breeding plumage are relatively subtle [41,42,43] and non-breeding C. alexandrinus and C. dealbatus mix frequently on the coast of the South China sea [41]. We further found a very small number of hybrids, and a lack of broad hybrid zone, implying some mechanisms to retain species boundary. Furthermore, we provide new data to support divergent morphometric characters and distinct ecological niches, which are likely key factors for maintaining species limits. Although we did not apply explicit Bayesian species delimitation analyses (e.g. [96]), our results offer compelling evidence that C. dealbatus is a distinct lineage from C. alexandrinus, and deserves a full species status under the General Lineage Concept of Species (GLC hereafter) [97, 98], nowadays widely adopted by taxonomists. The key essence of GLC is that species are defined as independently evolving metapopulation lineages, and can be assessed from an integrative taxonomic approach [99]. Beyond resolving the taxonomic puzzles between the two plovers, this study forms a basis to allow C. dealbatus to be considered as a target for conservation. Because it has a restricted range and is one of few breeding shorebirds in subtropical China coastline, where habitat loss may cause a risk of population decline [100].

Conclusions

Resolving the balance between diverging selection and gene flow is of fundamental importance to understand speciation processes. Here, we show that C. alexandrinus and C. dealbatus represent a case of incipient species in which divergent selection associated with ecological differences likely works as an efficient mechanism for the maintenance of divergence in the face of gene flow. While the full species status of C. dealbatus may be justified under the General Lineage Concept of Species, it remains untested whether a strong level of reproductive isolation has been established between C. alexandrinus and C. dealbatus. Furthermore, we have shown low genetic divergence between the two plovers, so it would be of particular importance to explore the patterns of divergence at the genome level and determine whether specific regions are related to reproductive isolation and adaptation. In addition, the IM analysis presented here estimated average historical gene flow, but it would be interesting to evaluate more realistic (and complex) demographic models in order to better estimate changes in Ne and the span of isolation and secondary contact. Genome-wide population genetic approaches promise the power to resolve the evolutionary history of the two Charadrius plovers but also open an avenue to characterize the genetic architecture associated with phenotypic trait divergence and local adaptation [84, 101].

Methods

Sampling

We mainly collected samples along the eastern coastal area of China (Fig. 1a), and also obtained samples from the two biggest continental islands, Taiwan and Hainan, as well as two small islands that are close to the Chinese coastal line (Zhoushan and Jinmen). Sampling sites were separated by approximately 200–250 km along a 2300 km transect, spanning almost the entire Chinese coastline. Furthermore, one high altitude population (breeding at 3350 m above sea level) near Qinghai Lake was sampled as an inland outgroup. Breeding individuals were captured using the funnel-trap method described in [102] between March and July in 2014–2015. For each nest in each locality, blood samples of the breeding pair in the most cases or a single chick were collected. Tissue samples were also obtained from dead individuals found in the field. Species identification performed during sampling was based on the summary of plumage characters as well as morphometric and ecological differences between C. alexandrinus and C. dealbatus, as described by [41]. Birds were released after morphometric measurement and blood collection. Overall, we collected samples from 454 individuals from 19 breeding sites.

Morphometrics

For each adult individual, four basic morphometric measurements were taken: body mass, wing length, tarsus length and bill length. Body mass was measured with an electronic scale (± 0.1 g). Wing length (flattened) was measured with a wing ruler (± 1 mm). Bill length to skull and tarsus lengths were measured using vernier callipers (± 0.1 mm). Measurements were taken by Q.H., P.J.Q. and C.Y.C. following the standard described in [103]. To avoid potential biases, some individuals were measured twice or three times by different authors between 20 and 30 trials to make sure there was no significant measurement difference (p < 0.05,) at the beginning and end of the fieldwork in the year 2014–2015.

We carried out principal component analysis (PCA) for the four traits to visualize the variation of morphometric data. Analyses of covariance (ANCOVAs) were conducted based on the PCA results to test statistical significance. We carried out t-tests for each measurement to assess the difference between species. Pairwise morphological difference QST based on PC1 scores from PCA were calculated following the method used by JF Storz [104] and plotted against linear distances between coastal populations. All aforementioned analyses were performed in PAST 3.12 [105].

Molecular genetic methods

We extracted genomic DNA using Tiangen Blood & Tissue Genome DNA Kits, following the manufacturer’s protocols (Tiangen, Beijing, China). DNA quality was measured with a NanoDrop 2000 (Thermo Scientific, USA). Three mtDNA loci, namely partial ATPase subunit 6/8, partial D-Loop of the mitochondrial control region (CR) and NADH dehydrogenase subunit 3 fragment (ND3), were amplified from all samples using primers from [42]. PCR reactions for mtDNA amplification were carried out in 20 μl volumes containing 1X PCR buffer (Takara Shuzo, Japan), 10–50 ng DNA template, 0.5 U Taq DNA Polymerase (Takara Shuzo, Japan), 1.0 μM of each primer, 2.0 mM of each dNTP and 1.0 μM MgCl2. We also sequenced 14 autosomal and two Z-linked exonic loci [106] for a subset of 40 individuals (Table 1) representing both species and several populations within species. PCR for nuclear loci was performed using a higher concentration of dNTP (4.0 mM) and MgCl2 (2.5 μM). The PCR cycling profiles used are listed in Additional file 1: Table S3. Each PCR product was checked on a 1% agarose gel prior to being sequenced on an ABI3730XL (Applied Biosystems, USA) by MajorBio, Shanghai, China.

Additionally, all samples were genotyped for 22 autosomal microsatellite markers (mst), mostly those genotyped by Küpper et al. [107], but also the markers C204 [108] and Hru2 [109], respectively. Microsatellite loci were amplified using three multiplex PCRs with respective cycling profiles (Additional file 1: Table S3). Each PCR was in a total volume of 15 μl containing 1X HotStart buffer (Tiangen), 10 ng DNA template, 1 unit Multi HotStart DNA Polymerase (Tiangen), 0.4 μM of each fluorescently labeled primer, 2.0 mM of each dNTP and 2.0 μM MgCl2. Multiplex PCR products and associated genotypes were isolated on an ABI3730XL (Applied Biosystems, USA) by Invitrogen, Shanghai, China and their length was determined using GeneMapper software v.3.7 (Applied Biosystems) against an internal size standard (GeneScan-500LIZ; Applied Biosystems).

Population genetic analyses

For DNA sequence data, we aligned each mitochondrial or nuclear locus using the CLUSTAL W algorithm in MEGA v.6.06 [110], and the alignment was checked by eye and manually edited if needed. For nuclear DNA sequence data, we first used PHASE 2.1.1 [111] to reconstruct haplotypes of nuclear sequences with heterozygous sites. Each run was set to 10,000 iterations, 100 burn-in and 10 thinning intervals. For both mtDNA and nuclear loci, basic genetic polymorphism statistics, such as haplotype number (h), haplotype diversity (Hd), number of segregating sites (S), nucleotide diversity (π) and Tajima’s D [112], of each locus and each population were calculated in DnaSP 5.10.1 [113]. Haplotype networks of each locus were constructed using a median-joining algorithm [114] in PopART 1.7.2 [115].

We used FreeNA [116] to check for the frequency of null alleles at each microsatellite locus. Further tests for Hardy-Weinberg equilibrium and pairwise linkage disequilibrium (LD) were carried out with Arlequin 3.1.1 [117]. Hardy-Weinberg equilibrium tests were run with 10,000 permutations. LD tests were run for 100,000 steps of the Markov chain. To obtain genetic diversity estimates, we calculated observed heterozygosity (Ho) and expected heterozygosity (HE) in GenAlEx 6.5.1 [118].

We estimated population structure between species and among sampling sites within species with several approaches. First, for the concatenated mtDNA sequences and microsatellites, we performed analyses of molecular variance (AMOVAs) in Arlequin to assess the proportion of genetic variance explained by the different partition settings: the two species; C. alexandrinus continental populations and Taiwan island populations; C. dealbatus, and coastal populations (including Hainan Island), Qinghai and Taiwan. Second, we also calculated pairwise ΦST and FST between breeding sites for mtDNA and mst, respectively, and we derived significance levels using 10,000 permutations in Arlequin.

For microsatellite genotypes only, we carried out assignment analyses with two model-based Bayesian approaches. First, we performed a Bayesian clustering analysis using the admixture model with correlated allele frequencies in STRUCTURE 2.3.4 [119, 120]. Ten independent analyses were run from K = 1 to K = 8 for 500,000 Markov chain Monte Carlo (MCMC) generations with 100,000 burn-in. Replicate runs were combined using STRUCTURE Harvester 0.6.94 [121] and CLUMPP 1.1.2 [122]. The most likely number of genetic clusters was also determined using Structure Harvester using the criteria described in Evanno et al. [123]. Second, we used a Bayesian clustering algorithm, that takes the geographical coordinates of each sampling site into account, using the R package, GENELAND 4.0 [124]. One million MCMC iterations were performed using a thinning setting of 100, from 1 to 10 populations, and a maximum Poisson process rate of 100. The uncertainty of special coordinates was set to zero and the maximum number of nuclei in the Poisson-Voronoi tessellation was 300. An independent Dirichlet distribution model for allele frequencies was also used. With the same package, the most likely number of clusters was determined based on their posterior density. To confirm the consistency of the results, we repeated the MCMC simulation 10 times.

Furthermore, to determine the probability that an individual was a hybrid or a migrant, we estimated the posterior probability of each individual based on microsatellite genotypes. First, we used HYBRIDLAB 1.0 [125] to simulate 100 individuals of each species, F1, F2 and back-crosses in each direction. Microsatellite genotypes of 260 individuals with a q-value higher than 95% in STRUCTURE analysis were used as the genotype pool where the genotypes of simulated parents were drawn from. One hundred simulated parents from each species and 100 F1 were used to calculate the threshold for individual hybrid assignment in STRUCTURE following the method used by J Vähä and C Primmer [126].

Historical demographic analyses

To infer the demographic histories of the two plover species, we applied Isolation with Migration model (IM) [82] analyses using the combined sequence data set of 16 nuclear loci and the concatenated mtDNA sequences based on 20 individuals of each species. IM analysis allows the inference of genealogies under different demographic scenarios and the estimation of population genetic parameters, such as divergence times, effective population sizes and migration rates between species since their divergence from the common ancestor. The homologous Killdeer Charadrius vociferus sequences from GenBank were used as an outgroup. The generation time was set to 2.5 years. The substitution rate of each locus was calculated using the method in Li et al. [127]. The ratio of net genetic distance of each locus across ingroup–outgroup was calculated, compared with net distance of mitochondrial cytochrome b (cytb) and then multiplied by the substitution rate for cytb (0.0105 ± 0.0005 substitution/site/mya) [128]. We used the Phi test [129] for recombination within nuclear loci and no recombination was detected. We implemented models in IMa2p [130], the parallel version of IMa2 [131]. For each analysis, we ran 48 MCMC chains for 2,500,000 steps of burn-in, after which 500,000 genealogies were saved, each recorded after 100 steps. Because IM analysis only provides estimates of average gene flow since the divergence of the two species, we used BayesAss 3.0.4 to perform a Bayesian analysis based on [132] based on microsatellite genotypes to characterize the level of gene flow within recent generations. We performed 10,000,000 iterations and 1,000,000 steps of burn-in.

Ecological niche analyses

To infer potential past range shifts induced by climatic changes, we carried out ecological niche modeling (ENM) using Maxent 3.3.3 k [133]. The geographically independent occurrence records of the two species of plovers (93 for C. alexandrinus and 24 for C. dealbatus) were obtained from our sampling expeditions during 2002–2016, China Bird Report (http://www.birdreport.cn), and the Global Biodiversity Information Facility (GBIF, http://www.gbif.org/). Furthermore, bioclimatic variables were obtained from the WorldClim database v.1.4 [134]. For highly correlated temperature/precipitation variable pairs (|r| ≥ 0.8) [135], we retained the variable with larger contribution to model development and putative greater biological importance (Additional file 1: Table S6). As a result, climatic conditions were measured as a function of eight bioclimatic variables (i.e., the mean diurnal range, isothermality, minimum temperature of the coldest month, mean temperature of the warmest quarter, annual precipitation, precipitation of the wettest month, precipitation of the driest month and precipitation seasonality, see Additional file 1: Appendix S1).

To explore niche similarity between the two species, we performed an ordination null test of PCA-env in environmental (E)-space [136, 137]. We used the methodology previously described by Broennimann and co-workers [136]. This method calculates the occurrence density and environmental factor density along environmental (principal component) axes for each cell using a kernel smoothing method and then uses the density of both occurrences and environmental variables to measure niche overlap along these axes. An unbiased estimate of Schoener’s D metric was calculated for our data, using smoothed densities from a kernel density function to measure niche overlap between the two species to ensure independence from the resolution of the grid. Statistical confidence in niche overlap values was then tested through a one-sided niche-similarity test. All statistical analyses were performed in R 3.0.2 [138]. Details of ENM construction and niche-similarity tests are available in Additional file 1: Appendix S1.

Stable-isotope analysis

Because the stable isotopic compositions of consumed tissues can be used to estimate the relative contribution of assimilated dietary sources [139], stable-carbon (C) and nitrogen (N) isotope analysis is widely used as a tool to study avian diet and migratory patterns [72, 140]. Carbon isotope ratios differ between C3, C4 and CAM plants due to differences in the photosynthetic pathways, and these differences are incorporated into an animal when the plants are consumed and so can be used to infer information about dietary niches [72] and geographic origins [141]. N isotopes are useful for identifying species/individuals which occupy different trophic positions (high δ15N implies higher trophic level [142]). In order to compare dietary differences based on differences in δ15N and δ13C between the two species, we collected the outer pair of rectrices from seven adults per site at eight sites: Qinghai Lake, Tangshan, Lianyungang and Rudong for C. alexandrinus; Fuzhou, Shanwei, Zhanjiang, Dongfang for C. dealbatus within our sample set. Since both species perform a complete post-breeding molt within their breeding grounds ([143], personal observation), isotope ratios represent trophic level and habitat preferences during the breeding period. We estimated niche width and overlap per species using an isotopic Bayesian approach based on δ13C and δ15N profile. Detailed information on lab procedures and statistical analyses can be found in Additional file 1: Appendix S2.

Availability of data and materials

Sequences deposited at GenBank: accession number MK830738- MK830957.

Morphometric and geo-referred occurrence data, stable isotope profiles and microsatellite genotypes available at: Dryad Doi: doi:https://doi.org/10.5061/dryad.3r01b8d

Abbreviations

AMOVAs:

Analyses of molecular variance

ANCOVA:

Analyses of covariance

CR:

The mitochondrial control region

cytb :

Cytochrome b

ENM:

Ecological niche modeling

h :

Haplotype number

Hd :

Haplotype

H E :

Expected heterozygosity diversity

Ho:

Observed heterozygosity

IM:

Isolation with Migration model

LD:

Linkage disequilibrium

LGM:

Last Glacial Maximum

LIG:

Last Interglacial

MCMC:

Markov chain Monte Carlo

mst:

Microsatellites

mtDNA:

Mitochondrial DNA

ND3:

NADH dehydrogenase subunit 3 fragment

Ne:

Effective population size

PCA:

Principal Component Analysis

S :

Number of segregating sites

π :

Nucleotide diversity

References

  1. 1.

    Coyne JA, Orr HA. Speciation, vol. 37. Sunderland: Sinauer Associates; 2004.

  2. 2.

    Wolf JB, Ellegren H. Making sense of genomic islands of differentiation in light of speciation. Nat Rev Genet. 2017;18(2):87.

  3. 3.

    Wagner CE, Harmon LJ, Seehausen O. Ecological opportunity and sexual selection together predict adaptive radiation. Nature. 2012;487(7407):366–9.

  4. 4.

    Nosil P, Feder JL. Genomic divergence during speciation: causes and consequences. In: Proceedings of the Royal Society of London; 2012.

  5. 5.

    Feder JL, Flaxman SM, Egan SP, Comeault AA, Nosil P. Geographic mode of speciation and genomic divergence. Annu Rev Ecol Evol Syst. 2013;44:73–97.

  6. 6.

    Mayr E. Animal species and evolution, vol. 797. Cambridge: Belknap Press of Harvard University Press; 1963.

  7. 7.

    Carson H, Clague D. Geology and biogeography of the Hawaiian Islands, Hawaiian Biogeography: Evolution on a Hot Spot Archipelago; 1995. p. 14–29.

  8. 8.

    Bernardi G, Alva-Campbell YR, Gasparini JL, Floeter SR. Molecular ecology, speciation, and evolution of the reef fish genus Anisotremus. Mol Phylogenet Evol. 2008;48(3):929–35.

  9. 9.

    Zhou Y, Duvaux L, Ren G, Zhang L, Savolainen O, Liu J. Importance of incomplete lineage sorting and introgression in the origin of shared genetic variation between two closely related pines with overlapping distributions. Heredity. 2016;118:211–20.

  10. 10.

    Le Moan A, Gagnaire PA, Bonhomme F. Parallel genetic divergence among coastal–marine ecotype pairs of European anchovy explained by differential introgression after secondary contact. Mol Ecol. 2016;25:3187–202.

  11. 11.

    Morales AE, Jackson ND, Dewey TA, O’Meara BC, Carstens BC. Speciation with gene flow in north American Myotis bats. Syst Biol. 2016;66(3):440–52.

  12. 12.

    Martin SH, Dasmahapatra KK, Nadeau NJ, Salazar C, Walters JR, Simpson F, Blaxter M, Manica A, Mallet J, Jiggins CD. Genome-wide evidence for speciation with gene flow in Heliconius butterflies. Genome Res. 2013;23(11):1817–28.

  13. 13.

    Rieseberg LH, Blackman BK. Speciation genes in plants. Ann Bot. 2010;106(3):439–55.

  14. 14.

    Meier JI, Marques DA, Mwaiko S, Wagner CE, Excoffier L, Seehausen O. Ancient hybridization fuels rapid cichlid fish adaptive radiations. Nat Commun. 2017;8:14363.

  15. 15.

    Fitzpatrick SW, Gerberich JC, Kronenberger JA, Angeloni LM, Funk WC. Locally adapted traits maintained in the face of high gene flow. Ecol Lett. 2015;18(1):37–47.

  16. 16.

    Nosil P, Schluter D. The genes underlying the process of speciation. Trends Ecol Evol. 2011;26(4):160–7.

  17. 17.

    Slatkin M. Gene flow and selection in a cline. Genetics. 1973;75(4):733–56.

  18. 18.

    Nosil P. Divergent host plant adaptation and reproductive isolation between ecotypes of Timema cristinae walking sticks. Am Nat. 2007;169(2):151–62.

  19. 19.

    Riesch R, Muschick M, Lindtke D, Villoutreix R, Comeault AA, Farkas TE, Lucek K, Hellen E, Soria-Carrasco V, Dennis SR. Transitions between phases of genomic differentiation during stick-insect speciation. Nat Ecol Evol. 2017;1(4):82.

  20. 20.

    Fraser BA, Axel K, Reznick DN, Christine D, Detlef W. Population genomics of natural and experimental populations of guppies (Poecilia reticulata). Mol Ecol. 2015;24(2):389–408.

  21. 21.

    Schluter D. Evidence for ecological speciation and its alternative. Science. 2009;323(5915):737–41.

  22. 22.

    Seehausen O, Butlin RK, Keller I, Wagner CE, Boughman JW, Hohenlohe PA, Peichel CL, Saetre G-P, Bank C, Brännström Å. Genomics and the origin of species. Nat Rev Genet. 2014;15(3):176–92.

  23. 23.

    Yang M, He Z, Shi S, Wu CI. Can genomic data alone tell us whether speciation happened with gene flow? Mol Ecol. 2017;26(11):2845–9.

  24. 24.

    Liu Y, Keller I, Heckel G. Breeding site fidelity and winter admixture in a long-distance migrant, the tufted duck (Aythya fuligula). Heredity. 2012;108-116(2):108.

  25. 25.

    Jackson DU, dos Remedios N, Maher KH, Zefania S, Haig S, Oyler-McCance S, Blomqvist D, Burke T, Bruford MW, Székely T. Polygamy slows down population divergence in shorebirds. Evolution. 2017;71(5):1313–26.

  26. 26.

    Eberhart-Phillips LJ, Hoffman JI, Brede EG, Zefania S, Kamrad MJ, Székely T, Bruford MW. Contrasting genetic diversity and population structure among three sympatric Madagascan shorebirds: parallels with rarity, endemism, and dispersal. Ecol Evol. 2015;5(5):997–1010.

  27. 27.

    Arguedas N, Parker PG. Seasonal migration and genetic population structure in house wrens. Condor. 2000;102(3):517–28.

  28. 28.

    Procházka P, Stokke BG, Jensen H, Fainová D, Bellinvia E, Fossøy F, Vikan JR, Bryja J, Soler M. Low genetic differentiation among reed warbler Acrocephalus scirpaceus populations across Europe. J Avian Biol. 2011;42(2):103–13.

  29. 29.

    Peters JL, McCRACKEN KG, Pruett CL, Rohwer S, Drovetski SV, Zhuravlev YN, Kulikova I, Gibson DD, Winker K. A parapatric propensity for breeding precludes the completion of speciation in common teal (Anas crecca, sensu lato). Mol Ecol. 2012;21(18):4563–77.

  30. 30.

    Winker K, McCracken KG, Gibson DD, Peters JL. Heteropatric speciation in a duck, Anas crecca. Mol Ecol. 2013;22(23):5922–35.

  31. 31.

    Kempenaers B, Valcu M. Breeding site sampling across the Arctic by individual males of a polygynous shorebird. Nature. 2017;541(7638):528–31.

  32. 32.

    Küpper C, Edwards SV, Kosztolanyi A, Alrashidi M, Burke T, Herrmann P, Arguelles-Tico A, Amat JA, Amezian M, Rocha A, et al. High gene flow on a continental scale in the polyandrous Kentish plover Charadrius alexandrinus. Mol Ecol. 2012;21(23):5864–79.

  33. 33.

    Trimbos KB, Musters C, Verkuil YI, Kentie R, Piersma T, de Snoo GR. No evident spatial genetic structuring in the rapidly declining black-tailed godwit Limosa limosa limosa in the Netherlands. Conserv Genet. 2011;12(3):629–36.

  34. 34.

    Miller MP, Gratto-Trevor C, Haig SM, Mizrahi DS, Mitchell MM, Mullins TD. Population genetics and evaluation of genetic evidence for subspecies in the Semipalmated sandpiper (Calidris pusilla). Waterbirds. 2013;36(2):166–78.

  35. 35.

    Conklin JR, Reneerkens J, Verkuil YI, Tomkovich PS, Palsbøll PJ, Piersma T. Low genetic differentiation between Greenlandic and Siberian sanderling populations implies a different phylogeographic history than found in red knots. J Ornithol. 2016;157(1):325–32.

  36. 36.

    Rönkä N, Kvist L, Pakanen VM, Rönkä A, Degtyaryev V, Tomkovich P, Tracy D, Koivula K. Phylogeography of the Temminck’s stint (Calidris temminckii): historical vicariance but little present genetic structure in a regionally endangered Palearctic wader. Divers Distrib. 2012;18(7):704–16.

  37. 37.

    del Hoyo J, Collar N, Kirwan GM, Sharpe CJ. White-faced Plover (Charadrius dealbatus). In: del Hoyo J, Elliott A, Sargatal J, Christie DA, de Juana E (editors). Handbook of the Birds of the World Alive. Lynx Edicions, Barcelona. 2019. Retrieved from https://www.hbw.com/node/467300.

  38. 38.

    Swinhoe R. On the plovers of the genus Aegialites found in China; 1870.

  39. 39.

    Hartert E, Jackson AC. Notes on some waders. Ibis. 1915;57(3):526–34.

  40. 40.

    IUCN Red List for birds. Downloaded from http://www.birdlife.org Accessed 8 Sept 2014.

  41. 41.

    Kennerley PR, Bakewell DN, Round PD. Rediscovery of a long-lost Charadrius plover from South-East Asia. Forktail. 2008;24:63–79.

  42. 42.

    Rheindt FE, Szekely T, Edwards SV, Lee PL, Burke T, Kennerley PR, Bakewell DN, Alrashidi M, Kosztolanyi A, Weston MA, et al. Conflict between genetic and phenotypic differentiation: the evolutionary history of a ‘lost and rediscovered’ shorebird. PLoS One. 2011;6(11):e26995.

  43. 43.

    de Hoyo J, Elliot A, Sargatal J, Christie D, de Juana E. Handbook of birds of the World Alive. Barcelona: Lynx Editions; 2019. http://www.hbw.com/ (Accessed Feb 2019)

  44. 44.

    Gill F, Donsker D: IOC World Bird List (v 9.1). 2019. Avaialble online at https://www.worldbirdnames.org/ioc-lists/crossref

  45. 45.

    Hey J, Nielsen R. Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics. 2004;167(2):747–60.

  46. 46.

    Won Y-J, Hey J. Divergence population genetics of chimpanzees. Mol Biol Evol. 2005;22(2):297–307.

  47. 47.

    Hey J. Recent advances in assessing gene flow between diverging populations and species. Curr Opin Genet Dev. 2006;16(6):592–6.

  48. 48.

    Lande R. Rapid origin of sexual isolation and character divergence in a cline. Evolution. 1982;36(2):213–23.

  49. 49.

    Hoskin CJ, Higgie M. Speciation via species interactions: the divergence of mating traits within species. Ecol Lett. 2010;13(4):409–20.

  50. 50.

    Railsback LB, Gibbard PL, Head MJ, Voarintsoa NRG, Toucanne S. An optimized scheme of lettered marine isotope substages for the last 1.0 million years, and the climatostratigraphic nature of isotope stages and substages. Quat Sci Rev. 2015;111:94–106.

  51. 51.

    Dos Remedios N, Lee PL, Burke T, Székely T, Küpper C. North or south? Phylogenetic and biogeographic origins of a globally distributed avian clade. Mol Phylogenet Evol. 2015;89:151–9.

  52. 52.

    Ni G, Li Q, Kong L, Yu H. Comparative phylogeography in marginal seas of the northwestern Pacific. Mol Ecol. 2014;23(3):534–48.

  53. 53.

    Taylor SA, White TA, Hochachka WM, Ferretti V, Curry RL, Lovette I. Climate-mediated movement of an avian hybrid zone. Curr Biol. 2014;24(6):671–6.

  54. 54.

    Taylor SA, Larson EL, Harrison RG. Hybrid zones: windows on climate change. Trends Ecol Evol. 2015;30(7):398–406.

  55. 55.

    Wang J, Ling MT, Dong YW. Causations of phylogeographic barrier of some rocky shore species along the Chinese coastline. BMC Evol Biol. 2015;15(1):114.

  56. 56.

    Ding S, Mishra M, Wu H, Liang S, Miyamoto MM. Characterization of hybridization within a secondary contact region of the inshore fish, Bostrychus sinensis, in the East China Sea. Heredity. 2017;120(1):51–62.

  57. 57.

    Liu JX, Gao TX, Wu SF, Zhang YP. Pleistocene isolation in the northwestern Pacific marginal seas and limited dispersal in a marine fish, Chelon haematocheilus (Temminck & Schlegel, 1845). Mol Ecol. 2010;16(2):275–88.

  58. 58.

    Ni G, Li Q, Kong L, Zheng X. Phylogeography of bivalve Cyclina sinensis: testing the historical glaciations and Changjiang River outflow hypotheses in northwestern Pacific. PLoS One. 2012;7(11):e49487.

  59. 59.

    Wang CH, Li CH, Li S. Mitochondrial DNA-inferred population structure and demographic history of the mitten crab (Eriocheir sensu stricto) found along the coast of mainland China. Mol Ecol. 2008;17(15):3515–27.

  60. 60.

    Bergmann C. Über die Verhältnisse der Wärmeökonomie der Thiere zu ihrer Grösse; 1848.

  61. 61.

    Salewski V, Watt C. Bergmann’s rule: a biophysiological rule examined in birds. Oikos. 2017;126(2):161–72.

  62. 62.

    McNab BK. Ecological factors affect the level and scaling of avian BMR. Comp Biochem Physiol A Mol Integr Physiol. 2009;152(1):22–45.

  63. 63.

    McWilliams SR, Karasov WH. Phenotypic flexibility in digestive system structure and function in migratory birds and its ecological significance. Comp Biochem Physiol A Mol Integr Physiol. 2001;128(3):577–91.

  64. 64.

    Alerstam T, Hedenström A, Åkesson S. Long-distance migration: evolution and determinants. Oikos. 2003;103(2):247–60.

  65. 65.

    Badyaev AV, Young RL, Oh KP, Addison C. Evolution on a local scale: developmental, functional, and genetic bases of divergence in bill form and associated changes in song structure between adjacent habitats. Evolution. 2008;62(8):1951–64.

  66. 66.

    Tattersall GJ, Arnaout B, Symonds MR. The evolution of the avian bill as a thermoregulatory organ. Biol Rev. 2016;92:1630–52.

  67. 67.

    Johnston RF. Character variation and adaptation in European sparrows. Syst Zool. 1969;18(2):206.

  68. 68.

    Nebel S, Rogers KG, Minton CD, Rogers DI. Is geographical variation in the size of Australian shorebirds consistent with hypotheses on differential migration? Emu. 2013;113(2):99–111.

  69. 69.

    Allen JA. The influence of physical conditions on the genesis of species. Radic Rev. 1877;1:108–40.

  70. 70.

    Que P, Chang Y, Eberhart-Phillips L, Liu Y, Székely T, Zhang Z. Low nest survival of a breeding shorebird in Bohai Bay, China. J Ornithol. 2015;156(1):297–307.

  71. 71.

    Cramp S, Simmons KL, editors BD, Collar N, Dunn E, Gillmor R, Hollom P, Hudson R, Nicholson E, Ogilvie M. Handbook of the birds of Europe, the Middle East and North Africa. The birds of the Western Palearctic: 3. Waders to gulls; 1983.

  72. 72.

    Hobson KA, Clark RG. Assessing avian diets using stable isotopes I: turnover of 13C in tissues. Condor. 1992;94:181–8.

  73. 73.

    McCormack JE, Zellmer AJ, Knowles LL. Does niche divergence accompany allopatric divergence in Aphelocoma jays as predicted under ecological speciation? Insights from tests with niche models. Evolution. 2010;64(5):1231–44.

  74. 74.

    Rundle HD, Nosil P. Ecological speciation. Ecol Lett. 2005;8(3):336–52.

  75. 75.

    Lynch M, Ackerman MS, Gout J-F, Long H, Sung W, Thomas WK, Foster PL. Genetic drift, selection and the evolution of the mutation rate. Nat Rev Genet. 2016;17(11):704–14.

  76. 76.

    Ellstrand NC, Elam DR. Population genetic consequences of small population size: implications for plant conservation. Annu Rev Ecol Syst. 1993;24(1):217–42.

  77. 77.

    Lanfear R, Kokko H, Eyre-Walker A. Population size and the rate of evolution. Trends Ecol Evol. 2014;29(1):33–41.

  78. 78.

    Poelstra JW, Vijay N, Bossu CM, Lantz H, Ryll B, Müller I, Baglione V, Unneberg P, Wikelski M, Grabherr MG. The genomic landscape underlying phenotypic integrity in the face of gene flow in crows. Science. 2014;344(6190):1410–4.

  79. 79.

    García-Navas V, Bonnet T, Waldvogel D, Wandeler P, Camenisch G, Postma E. Gene flow counteracts the effect of drift in a Swiss population of snow voles fluctuating in size. Biol Conserv. 2015;191:168–77.

  80. 80.

    Zink RM, Barrowclough GF. Mitochondrial DNA under siege in avian phylogeography. Mol Ecol. 2010;17(9):2107–21.

  81. 81.

    Irwin DE. Sex chromosomes and speciation in birds and other ZW systems. Mol Ecol. 2018;27(19):3831–51.

  82. 82.

    Hey J, Nielsen R. Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proc Natl Acad Sci U S A. 2007;104(8):2785–90.

  83. 83.

    Qvarnström A, Bailey RI. Speciation through evolution of sex-linked genes. Heredity. 2009;102(1):4–15.

  84. 84.

    Nadachowska-Brzyska K, Burri R, Olason PI, Kawakami T, Smeds L, Ellegren H. Demographic divergence history of pied flycatcher and collared flycatcher inferred from whole-genome re-sequencing data. PLoS Genet. 2013;9(11):e1003942.

  85. 85.

    Toews DP, Taylor SA, Vallender R, Brelsford A, Butcher BG, Messer PW, Lovette IJ. Plumage genes and little Else distinguish the genomes of hybridizing warblers. Curr Biol. 2016;26(17):2313–8.

  86. 86.

    Cruickshank TE, Hahn MW. Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow. Mol Ecol. 2014;23(13):3133–57.

  87. 87.

    Burri R. Linked selection, demography and the evolution of correlated genomic landscapes in birds and beyond. Mol Ecol. 2017;26(15):3853–6.

  88. 88.

    Harrison RG. Pattern and process in a narrow hybrid zone. Heredity. 1986;56(3):337–49.

  89. 89.

    Bastos-Silveira C, Santos SM, Monarca R, Mathias ML, Heckel G. Deep mitochondrial introgression and hybridization among ecologically divergent vole species. Mol Ecol. 2012;21(21):5309–23.

  90. 90.

    Anderson EC, Thompson EA. A model-based method for identifying species hybrids using multilocus genetic data. Genetics. 2002;160(3):1217–29.

  91. 91.

    Bearhop S, Fiedler W, Furness RW, Votier SC, Waldron S, Newton J, Bowen GJ, Berthold P, Farnsworth K. Assortative mating as a mechanism for rapid evolution of a migratory divide. Science. 2005;310(5747):502–4.

  92. 92.

    Brown WL, Wilson EO. Character displacement. Syst Zool. 1956;5(2):49–64.

  93. 93.

    Robinson BW, Wilson DS. Character release and displacement in fishes: a neglected literature. Am Nat. 1994;144:596–627.

  94. 94.

    Rybinski J, Sirkiä PM, McFarlane SE, Vallin N, Wheatcroft D, Ålund M, Qvarnström A. Competition-driven build-up of habitat isolation and selection favoring modified dispersal patterns in a young avian hybrid zone. Evolution. 2016;70(10):2226–38.

  95. 95.

    Reifová R, Reif J, Antczak M, Nachman MW. Ecological character displacement in the face of gene flow: Evidence from two species of nightingales. BMC Evol Biol. 2011;11(1):138.

  96. 96.

    Yang Z, Rannala B. Unguided species delimitation using DNA sequence data from multiple loci. Mol Biol Evol. 2014;31(12):3125–35.

  97. 97.

    De Queiroz K. The general lineage concept of species, species criteria, and the process of speciation. In: Howard DJ, Berlocher SH, editors. Endless forms: species and speciation. New York: Oxford University Press; 1998. p. 57–75.

  98. 98.

    De QK. Species concepts and species delimitation. Syst Biol. 2007;56(6):879–86.

  99. 99.

    Fujita MK, Al E. Coalescent-based species delimitation in an integrative taxonomy. Trends Ecol Evol. 2012;27(9):480–8.

  100. 100.

    Ma Z, Melville DS, Liu J, Chen Y, Yang H, Ren W, Zhang Z, Piersma T, Li B. Rethinking China's new great wall. Science. 2014;346(6212):912–4.

  101. 101.

    Freedman AH, Gronau I, Schweizer RM, Ortega-Del Vecchyo D, Han E, Silva PM, Galaverni M, Fan Z, Marx P, Lorente-Galdos B. Genome sequencing highlights the dynamic early history of dogs. PLoS Genet. 2014;10(1):e1004016.

  102. 102.

    Székely T, Kosztolányi A, Küpper C. Practical guide for investigating breeding ecology of Kentish plover. 2008. Available at https://www.researchgate.net/publication/228494424_Practical_guide_for_investigating_breeding_ecology_of_Kentish_plover_Charadrius_alexandinus.

  103. 103.

    Redfern CP, Clark JA, Wilson A, Gough S, Robertson D. Ringers manual: British trust for ornitology; 2001.

  104. 104.

    Storz JF. Contrasting patterns of divergence in quantitative traits and neutral DNA markers: analysis of clinal variation. Mol Ecol. 2002;11(12):2537–51.

  105. 105.

    Hammer Ø, Harper D, Ryan P: PAST: Paleontological Statistics Software Package for education and data analysis. Palaeontolia Electronica 4. 2001.

  106. 106.

    Liu Y, Liu S, Yeh C-F, Zhang N, Chen G, Que P, Dong L, Li S-H. The first set of universal nuclear protein-coding loci markers for avian phylogenetic and population genetic studies. Sci Rep. 2018;8(1):15723.

  107. 107.

    Küpper C, Horsburgh GJ, Dawson DA, Ffrench-Constant R, Szekely T, Burke T. Characterization of 36 polymorphic microsatellite loci in the Kentish plover (Charadrius alexandrinus) including two sex-linked loci and their amplification in four other Charadrius species. Mol Ecol Notes. 2007;7(1):35–9.

  108. 108.

    Funk WC, Mullins TD, Haig SM. Conservation genetics of snowy plovers (Charadrius alexandrinus) in the Western hemisphere: population genetic structure and delineation of subspecies. Conserv Genet. 2007;8(6):1287–309.

  109. 109.

    Primmer C, Møller A, Ellegren H. Resolving genetic relationships with microsatellite markers: a parentage testing system for the swallow Hirundo rustica. Mol Ecol. 1995;4(4):493–8.

  110. 110.

    Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9.

  111. 111.

    Stephens M, Smith NJ, Donnelly P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001;68(4):978–89.

  112. 112.

    Tajima FV, Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123(3):585–95.

  113. 113.

    Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25(11):1451–2.

  114. 114.

    Bandelt H-J, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16(1):37–48.

  115. 115.

    Leigh JW, Bryant D. Popart: full-feature software for haplotype network construction. Methods Ecol Evol. 2015;6(9):1110–6.

  116. 116.

    Chapuis M-P, Estoup A. Microsatellite null alleles and estimation of population differentiation. Mol Biol Evol. 2007;24(3):621–31.

  117. 117.

    Excoffier L, Laval G, Schneider S. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinforma. 2005;1:47–50.

  118. 118.

    Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in excel. population genetic software for teaching and research--an update. Bioinformatics. 2012;28(19):2537–9.

  119. 119.

    Pritchard JK, Stephens MJ, Donnelly PJ. STRUCTURE, version 2.3.3, inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–59.

  120. 120.

    Hubisz MJ, Falush D, Stephens M, Pritchard JK. Inferring weak population structure with the assistance of sample group information. Mol Ecol Resour. 2009;9(5):1322–32.

  121. 121.

    Earl DA. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4(2):359–61.

  122. 122.

    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(14):1801–6.

  123. 123.

    Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14(8):2611–20.

  124. 124.

    Guillot G, Mortier F, Estoup A. GENELAND: a computer package for landscape genetics. Mol Ecol Notes. 2005;5(3):712–5.

  125. 125.

    Nielsen EE, Bach LA, Kotlicki P. HYBRIDLAB (version 1.0): a program for generating simulated hybrids from population samples. Mol Ecol Notes. 2006;6(4):971–3.

  126. 126.

    Vähä J, Primmer C. Efficiency of model-based Bayesian methods for detecting hybrid individuals under different hybridization scenarios and with different numbers of loci. Mol Ecol. 2006;15(1):63–72.

  127. 127.

    Li JW, Yeung CK, Tsai P-W, Lin RC, Yeh C-F, Yao CT, Han L, Hung L, Ding P, Wang Q. Rejecting strictly allopatric speciation on a continental island: prolonged postdivergence gene flow between Taiwan (Leucodioptron taewanus, Passeriformes Timaliidae) and Chinese (L. canorum canorum) hwameis. Mol Ecol. 2010;19(3):494–507.

  128. 128.

    Weir J, Schluter D. Calibrating the avian molecular clock. Mol Ecol. 2008;17(10):2321–8.

  129. 129.

    Bruen TC, Philippe H, Bryant D. A simple and robust statistical test for detecting the presence of recombination. Genetics. 2006;172(4):2665–81.

  130. 130.

    Sethuraman A, Hey J. IMa2p–parallel MCMC and inference of ancient demography under the isolation with migration (IM) model. Mol Ecol Resour. 2016;16(1):206–15.

  131. 131.

    Hey J. Isolation with migration models for more than two populations. Mol Biol Evol. 2010;27(4):905–20.

  132. 132.

    Stephens M, Donnelly P. A comparison of bayesian methods for haplotype reconstruction from population genotype data. Am J Hum Genet. 2003;73(5):1162–9.

  133. 133.

    Phillips SJ, Dudík M. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography. 2008;31(2):161–75.

  134. 134.

    Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25(15):1965–78.

  135. 135.

    Qiao H, Peterson AT, Ji L, Hu J. Using data from related species to overcome spatial sampling bias and associated limitations in ecological niche modelling. Methods Ecol Evol. 2017;8(12):1804–12.

  136. 136.

    Broennimann O, Fitzpatrick MC, Pearman PB, Petitpierre B, Pellissier L, Yoccoz NG, Thuiller W, Fortin M-J, Randin C, Zimmermann NE, et al. Measuring ecological niche overlap from occurrence and spatial environmental data. Glob Ecol Biogeogr. 2012;21(4):481–97.

  137. 137.

    Hu J, Broennimann O, Guisan A, Wang B, Huang Y, Jiang J. Niche conservatism in Gynandropaa frogs on the southeastern Qinghai-Tibetan plateau. Sci Rep. 2016;6:32624.

  138. 138.

    R: a language and environment for statistical computing. https://www.r-project.org.

  139. 139.

    DeNiro MJ, Epstein S. Carbon isotopic evidence for different feeding patterns in two hyrax species occupying the same habitat. Science. 1978;201(4359):906–8.

  140. 140.

    Pagani-Núñez E, Renom M, Mateos-Gonzalez F, Cotín J, Senar JC. The diet of great tit nestlings: comparing observation records and stable isotope analyses. Basic Appl Ecol. 2017;18:57–66.

  141. 141.

    Hobson KA. Tracing origins and migration of wildlife using stable isotopes: a review. Oecologia. 1999;120(3):314–26.

  142. 142.

    Chen Q, Liu Y, Ho W-T, Chan SK, Li Q-h, Huang J-R. Use of stable isotopes to understand food webs in Macao wetlands. Wetl Ecol Manag. 2017;25(1):59–66.

  143. 143.

    Ginn H, Melville D: Moult in birds. BTO guide 19. British Trust for Ornithology, Tring 1983.

Download references

Acknowledgements

We thank Qiaoyi Liang, Xiaoyan Long, Derong Meng, Demeng Jiang, Li Tian and Hebo Peng for their help with field sampling and Guoling Chen and Yangkai Zhou for their assistance with laboratory work. We are also grateful to Dr. Jennifer McDowall for editing the English text of a draft of this manuscript. Special thanks to Shaochong Peng for preparing sketches of plovers in Fig. 1. We are very grateful to the three anonymous reviewers and Dr. Robert Kraus for their constructive comments and suggestions.

Funding

This study was supported by the National Natural Science Foundation of China, Grant/Award Number: 31301875 and 31572251 to Y.L., 31600297 to P.J.Q. and 31572288 to Z.W.Z., the Open Grant of the State Key Laboratory of Biocontrol of Sun Yat-sen University to Y.L. and T.Z. (SKLBC13KF03) and the Youth Innovation Promotion Association CAS (2015304) to J.H.H.. Computational machinery work was granted by the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase) under Grant No. U1501501 to Y.L. We declare that these funding bodies played no role in the design of the study, the collection, analysis, and interpretation of data and in writing the manuscript.

Author information

YL, ZWZ and TS conceived and designed the study. PJQ, C-YC, QH, JM, NZ, XJW, and TS collected data in the field. XJW, SML, and ZN conducted molecular genetic laboratory work. XJW and PJQ analyzed genetic data with input from YL and GH. JHH performed niche modeling analysis, XCZ carried out stable isotope analysis with input from EP-G, YYL (Yu Yan Leung) and CD. XJW, PJQ and YL wrote the manuscript with input from JHH, GH, XCZ and CD. All authors gave final approval for publication.

Correspondence to Yang Liu.

Ethics declarations

Ethics approval and consent to participate

This study did not involve endangered or protected species. The protocol of our study, including blood collection procedures, measurements and stable isotope analysis was reviewed and approved by the Institutional Ethical Committee of Animal Experimentation of Sun Yat-sen University (2005DKA21403-JK). We strictly complied with the ethical conditions by the Chinese Animal Welfare Act (20090606). Also, bird captures, banding, measurement and sample collection were performed with permission from the respective authorities (Mainland China: Beijing Normal University to P.J.Q., which covered the study site 2–9 in Fig. 1, and Sun Yat-sen University to Y.L., which covered the study site 1, 12–19 Fig. 1; Taiwan: Changhua County Government, New Taipei City Government and Kinmen County Government to C.Y.C., which covered the study site 10–11 Fig. 1).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1:

Table S1. Estimates of h, number of haplotypes; s, number of segregating sites; π, nucleotide diversity; Tajima’s D value, substitution rate of each locus estimated based on the substitution rate for cytb and GenBank accession number for each mtDNA and nuclear locus of C. alexandrinus and C. dealbatus were provided. The substitution rate of each locus was calculated using the method in Li et al. (2010): the ratio of net genetic distance of each locus across ingroup–outgroup was calculated, compared with net distance of mitochondrial cytochrome b (cytb) and then multiplied by the substitution rate for cytb (0.0105 ± 0.0005 substitution/site/mya, Weir & Schluter 2008). Table S2. Estimates of pairwise FST for each microsatellite locus between C. alexandrinus and C. dealbatus, LD test result (same for both species), Na, number of alleles, Ho, observed heterozygosity and He, expected heterozygosity in each species. Table S3. PCR amplification protocols for three mtDNA, 16 nuclear exons and 22 microsatellite loci genotyped for C. alexandrinus and C. dealbatus. Table S4. Genetic differentiation between each pair of sampling localities of C. alexandrinus and C. dealbatus. Estimates were based on 1729 bp mtDNA (lower diagonal, ΦST) and 13 microsatellite loci (upper diagonal, FST). Values highlighted in red represent significant value after Bonferroni correction. Table S5 The morphological measurements of C. alexandrinus and C. dealbatus in different breeding sites along the Chinese coast, and Taiwan and Hainan Islands. Site reference number corresponds to numbers in Fig. 1 and Table 1. The number of sample size for each site (n), the mean value of measurements including bill and wing length, body mass and their respective standard deviations (SD) are given. Only DNA samples from a breeding site in Cangzhou were collected but not the measurement, the corresponding data is missing for this population. Mean and SD of bill length, wing length and body mass for each population. Measurement data of the site 3-Cangzhou was not collected. Table S6. Pearson’s correlation coefficients between pairwise of bioclimatic variables with the ranges of C. alexandrinus and C. dealbatus. Figure S1. Haplotype networks based on 80 individuals of 12 nuclear loci not shown in Fig. 2. C. alexandrinus (blue) and C. dealbatus (yellow). Figure S2. The Bayesian clustering analysis with STRUCTURE clearly suggested two genetic clusters corresponding to C. alexandrinus and C. dealbatus. Shown is the maximum value of the Delta K (Δ K) in posterior likelihood Ln P (X/K) over 10 runs per K of STRUCTURE. Figure S3. Niche of C. alexandrinus and C. dealbatus in climatic space from a principal component analysis (PCA-env). a) and b) show the niche characteristics of C. alexandrinus and C. dealbatus, respectively, along the first two axes of the PCA. Grey shading shows the density of the occurrences of the species by cell. The solid and dashed contour lines illustrate, respectively, 100 and 50% of the available (background) environment. c) The contribution of the variables on the first two axes of the PCA and the percentage of inertia explained by the two axes. d) Observed niche overlap D between the two ranges (bars with a diamond) and simulated niche overlaps (grey bars) on which tests of niche equivalency were calculated with 100 iterations. Figure S4. Differentiation in the stable isotope ratios among breeding sites (top-left: δ13C and bottom-left: δ15N) and between the two plovers (top-right: δ13C and bottom-right: δ15N). In all representation, C. alexandrinus marked in blue and C. dealbatus in yellow. Figure S5. (a) Newhybrid simulations failed to reliably detect simulated hybrid individuals. Real data of individuals with posterior density higher than 95% in STRUCTURE, and data of 5 simulated F1 and F2 individuals, and 10 back-cross individuals on each direction was used to imitate a situation when hybrids were the minority in the population. (b) Newhybrid results from real data are highly consistent with the STRUCTURE results. KP represents C. alexandrinus, WFP represents C. dealbatus, bx represents back-crosses. Appendix S1. Reconstruction of potential range shifts induced by climatic changes and inference of environmental niche overlap between the two-plover species, Charadrius alexandrinus and C. dealbatus. Appendix S2. Inference of interspecific diet overlap using stable-isotope analysis between the two plover species, Charadrius alexandrinus and C. dealbatus. (DOCX 530 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Parapatry
  • Character displacement
  • Gene flow
  • Hybridization
  • Stable isotope analysis
  • Ecological niche