Introduction

The delineation of conservation units such as evolutionary significant units (ESUs, Moritz 1994) and management units (MUs) can be highly nuanced. Such units can be difficult to define in practice because of similarities in species-specific ecology and life history traits that are likely to promote cohesion among populations across a broad landscape. For example, highly vagile taxa such as pelagic organisms and species capable of flight may have high rates of gene flow, and are therefore unlikely to have many reciprocally monophyletic groups based on their mitochondrial DNA (mtDNA) sequences (Crandall et al. 2000; Medina et al. 2018).

Traditionally, studies aiming to identify conservation units using genetics have sampled a small number of variable genetic loci such as mtDNA and/or microsatellites. However, a greater number of genetic markers increases the power to detect fine scale population structure, and hence, increases the reliability of defining conservation units and their associated characteristics such as past demographic changes and effective population sizes (Allendorf et al. 2010; Funk et al. 2012; Shafer et al. 2015). Genomic methods, including restriction enzyme associated DNA sequencing (RAD-seq, Baird et al. 2008; Miller et al. 2007), produce data from hundreds or thousands rather than tens of loci distributed throughout the genome, and can reveal previously uncharacterised genetic structure (e.g. Ruegg et al. 2014; Saenz-Agudelo et al. 2015; Barth et al. 2017; Vendrami et al. 2017; Younger et al. 2017; Sadanandan et al. 2019). In doing so, genomic methods have enabled novel, high resolution identification of new management units in species of high conservation priority (Palsbøll et al. 2007; Kjeldsen et al. 2016; Peters et al. 2016; Barth et al. 2017; Younger et al. 2017). In addition, these methods have proven successful in corroborating previously proposed taxonomic or conservation units that were first identified using low-resolution genetic data (Mason and Taylor 2015; Attard et al. 2018; Doyle et al. 2018) or ecological variation (Lemay and Russello 2015; Prince et al. 2017).

Once defined, the vulnerability of different conservation units to population declines can be assessed, in part, by using genetic data to reconstruct past demographic changes (e.g. Carnaval et al. 2009; Palsbøll et al. 2013; Garrick et al. 2015; Shafer et al. 2015; Stoffel et al. 2018). Importantly, the accuracy of genetic demographic reconstructions can be enhanced with data from a larger number of loci (Shafer et al. 2015), as well as using multiple modelling methods to corroborate findings (Brüniche-Olsen et al.2014).

The snowy plover, Charadrius nivosus (Cassin, 1858; Küpper et al. 2009), is one of the least abundant shorebird species endemic to the Americas, but also one of the most thoroughly studied. This iconic, partially migratory shorebird exhibits an unusual breeding system of sequential polyandry where females desert their families to re-mate with a new partner, and inhabits salt flats of alkaline lakes, coastal lagoons and sandy beaches (Warriner et al. 1986; Page et al. 2009; Eberhart-Phillips et al. 2017). Coastal populations can be made up of migratory and resident individuals, however, inland breeding populations tend to migrate during the winter towards coastal habitat (Page et al. 2009). In addition, both sexes occasionally undergo long within and between breeding season movements (50–1140 km; Stenzel et al. 1994). The snowy plover subspecies taxonomy is debated, with some authorities (Clements et al. 2018; del Hoyo et al. 2019) recognising two subspecies (Charadius nivosus nivosus and C. n. occidentalis), whereas, a third, C. n. tenuirostris (American Ornithologists’ Union 1957; Funk et al. 2007), is currently recognised by the United States Fish and Wildlife Service (USFWS 2011) although this subspecies shows limited phenotypic distinction from C. n. nivosus (Funk et al. 2007). C. n. nivosus (nivosus herein) is found within the United States and Mexico, the core distribution of C. n. tenuirostris (tenuirostris herein) comprises the Caribbean Islands and Bermuda, whereas, C. n. occidentalis (occidentalis herein) is found along the west coast of South America from Colombia to Chile (Fig. 1).

Fig. 1
figure 1

Sampling locations of snowy plovers (Charadrius nivosus) included in each of the four genetic datasets (see legend). Deme identity is given by colour (blue: western nivosus, light blue: eastern nivosus, red: tenuirostris, yellow: occidentalis). Species range is indicated in dark grey shading. For details on sample sizes see Table 1 (all markers) and Table S2 (mtDNA only)

Full size image

The most recent genetic analysis of snowy plovers supported the recognition of three subspecies (nivosus, tenuirostris, occidentalis) and updated their geographic distributions (Funk et al. 2007). Specifically, snowy plovers breeding in Florida and along the Gulf of Mexico coast, which had been previously assigned to tenuirostris, were more genetically similar to other populations of nivosus (Funk et al. 2007) than to the island populations of tenuirostris. However, a more recent analysis based on microsatellite loci with additional sampling of Mexican populations showed that the eastern nivosus snowy plover population from Florida may be differentiated from western nivosus snowy plovers (D’Urban Jackson et al. 2017). Among the western snowy plover populations, molecular analyses did not support distinct conservation units between the Pacific coast and inland populations (Funk et al. 2007; D’Urban Jackson et al. 2017). Despite this genetic evidence, the conservation status of Pacific snowy plovers as a distinct population segment under the Endangered Species Act has been maintained because available banding data have not confirmed movement between inland and coastal sites (Haig et al. 2011). Genetic homogeneity of western populations may be the result of gene flow from long distance breeding dispersal of females (D’Urban Jackson et al. 2017), which has been observed in Pacific populations (Stenzel et al. 1994). Overall, previous molecular studies using microsatellites and mtDNA sequences indicate low genetic diversity across all snowy plover subspecies (Funk et al. 2007; D’Urban Jackson et al. 2017). Furthermore, inbreeding has been detected in a small isolated population of snowy plovers through a long-term pedigree study (Colwell and Pearson 2011).

The estimated census breeding population sizes of the three subspecies are 25,869 (nivosus Thomas et al. 2012,), 1000–10,000 (occidentalis, Wetlands International 2019) and < 100 (tenuirostris, Elliot-Smith et al.2004; Brown 2012). Many coastal populations are declining because of habitat deterioration and human disturbance (Colwell et al. 2007; Küpper et al. 2011; Powell and Collier 2011; Page et al. 2009; Thomas et al. 2012; Cohen et al. 2014; Galindo-Espinosa and Palacios 2015; Cruz-López et al. 2017). As a result, many populations are protected by federal and/or state-level statutes (Haig et al. 2011; Cohen et al. 2014; Galindo-Espinosa and Palacios 2015). Population viability models suggest that the decline of coastal populations will become more severe throughout the species’ range without intensive conservation management (Aiello-Lemmens et al.2011; Eberhart-Phillips and Colwell 2014; Cruz-López et al.2017). Furthermore, climate change might exacerbate risks of local extinctions due to rising sea levels and changes in the dynamics of tropical cyclone events, which are important for creating suitable nesting habitats (Aiello-Lemmens et al.2011; Convertino et al. 2011).

Snowy plover breeding habitat includes sandy beaches, making the conservation management of this species in North America controversial, as habitat protection often conflicts with development and recreational activities (Lafferty et al. 2006). Regular reassessment of population structure in snowy plovers by implementing current methodological advances is therefore required to guide effective conservation policy. Both previous genetic analyses (Funk et al. 2007; D’Urban Jackson et al. 2017) used ≤ 15 nuclear loci, which may have been insufficient to detect fine scale population structure in this highly vagile species. Here, we comprehensively re-evaluate genetic differentiation among snowy plover populations using four different genetic data sets: mtDNA, microsatellites, autosomal and Z-linked single-nucleotide polymorphisms (SNPs), respectively. We use these data to: (1) examine genetic support for the current subspecies delineation, examine population structure within the most widespread subspecies nivosus and identify putative distinct conservation units based on genetic data; (2) compare the sensitivity of the four different types of genetic data to detect fine scale genetic structure and spatial genetic patterns among controversial population segments; and (3) reconstruct the demographic history of snowy plover populations using three modelling methods to understand the origin of low genetic diversity.

Materials and methods

Combining previously published (Funk et al. 2007; Küpper et al. 2009; D’Urban Jackson et al. 2017) and newly collected data, we gathered genetic data from all three proposed subspecies (genetic lineages herein): nivosus, occidentalis, and tenuirostris. Sampling sites represented in each genetic dataset are shown in Fig. 1 and sample information is summarised in Table 1, with a more detailed description of each data set below.

Table 1 Genetic diversity of snowy plover (Charadrius nivosus) populations using four genetic marker datasets: autosomal and Z-linked ddRAD loci, microsatellites, and mtDNA. Shaded cells indicate populations/demes with genetic data for all four markers
Full size table

mtDNA

We downloaded available snowy plover sequences of the mitochondrial D-loop (control region) from GenBank (n = 114; Table S1). In addition, we sequenced the D-loop for 67 new samples using primers and polymerase chain reaction (PCR) conditions described in Küpper et al. (2012). We aligned all sequences using the ClustalW alignment tool in BioEdit v7.2.5 (Hall 1999) and trimmed the alignment to 424 base pairs (bp) for downstream analyses. For clarity, only locations with both mtDNA and additional genetic information (microsatellites or single nucleotide polymorphism markers [SNPs]) are presented in Table 1 and Fig. 1. Details for mtDNA only samples, including GenBank accession numbers, are given in Table S1.

Microsatellites

Additional samples from occidentalis (n = 18) and tenuirostris (n = 13) were added to those described in the microsatellite dataset for nivosus published by D’Urban Jackson et al. (2017, n = 146). We genotyped all samples at 15 microsatellite loci (including seven additional loci in comparison to D’Urban Jackson et al. (2017)) and removed first order relatives from the data set based on field records, i.e. we retained only one chick per brood. Primers, PCR conditions and reagents are described in Küpper et al. (2008, 2009) and D’Urban Jackson et al. (2017).

Autosomal and Z-linked SNPs

We extracted DNA from blood samples using the phenol–chloroform extraction method (Sambrook and Russell 2006). We prepared a multiplexed fragment library for double digest restriction site associated DNA sequencing (ddRAD) for 47 samples following the protocol described by DaCosta and Sorenson (2014) using EcoRI and SbfI restriction enzymes and adaptors incorporating an individual 6 bp barcode. The multiplexed library was paired-end sequenced (100 bp reads) on one lane of an Illumina HiSeq 2000 platform in the Bauer Core Facility at Harvard University.

Bioinformatics pipeline

We discarded the paired read to simplify downstream data processing and checked the quality of the raw single-end reads using FastQC (v.0.11.5; Babraham Bioinformatics; https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Using the ddRAD software pipeline created by DaCosta and Sorenson (2014; https://github.com/BU-RAD-seq/), we first de-multiplexed sequences based on their unique barcodes. For each individual, we condensed identical sequences while maintaining a count of the sequences, and then filtered low quality, singleton sequences using the UCLUST function in USEARCH v5 (Edgar 2010). Condensed and filtered reads across all samples were then clustered at an 85% identity level using UCLUST. We compared the highest quality sequence from each cluster to the closest available reference genome (killdeer Charadrius vociferous Zhang et al. 2014,) using BLASTN (Altschul et al. 1990) and then combined clusters with identical or nearly identical hits. We aligned sequences within each cluster (i.e. putative locus) using MUSCLE v3 (Edgar 2004a, b). Finally, we scored genotypes for each individual at each locus using the python script RADgenotypes.py.

We flagged individual genotypes that had: (i) fewer than three reads, (ii) an uneven ratio between two alleles, with the minor allele having a proportion less than 0.3 (if unique across all samples) or 0.2, and (iii) “extra” reads appearing with a proportion ≥ 0.1, implying the presence of more than two alleles. We removed clusters if they: (i) had more than five samples with flagged genotypes, (ii) had a length < 19 bp, (iii) had a median per-sample depth less than 9.5, or (iv) had more than five samples with missing data. Clusters with potential alignment errors (see manual at https://github.com/BU-RAD-seq/ddRAD-seq-Pipeline) were also removed, as were clusters with > 4 SNPs/indels, indels of length > 4 bp, and very high read depths (> 800 reads total across all individuals), as these characteristics may be indicative of multiple genomic loci being clustered together.

Z-linked RAD loci

We identified SNPs on the Z chromosome using a two-step method. First, we compared mean read depths for males and females for each quality filtered RAD locus, with the underlying expectation that sequencing depth for female samples should be about half that of male samples for loci on the Z chromosome, as females are hemizygous (ZW) and males are ZZ in birds. RAD loci with a female to male average depth ratio between 0.4 and 0.65 were then checked for alignment to the Z chromosome of the killdeer genome (Zhou et al. 2014) using our BLASTN results. RAD loci with both read depth and alignment support for Z-linkage were assigned to the Z chromosome dataset. RAD loci that were only “Z-linked” by one of the two methods, were removed from all datasets. RAD loci only present in females (i.e., indicative of being located on the W chromosome) were excluded by the above filter (i.e., missing data in five or more samples). We assigned all other clusters/loci to an “autosomal dataset.”

SNP filtering

For SNP based analyses we retained one randomly selected biallelic SNP per RAD locus to minimise linkage. Finally, to reduce potential biases in Bayesian clustering solutions (Linck and Battey 2017), singleton SNPs were excluded from population structure analyses. For all other analyses, we included singleton SNPs. For diversity, correlation tests between geographic and genetic distances, and PCA analyses of the Z-linked SNP dataset we included only (diploid) males, whereas all other analysis included all individuals.

Genetic diversity and differentiation

For mtDNA we used DNAsp v6 (Librado and Rozas 2009) to calculate haplotype diversity. We calculated microsatellite allelic richness in MSA v4.05 (Dieringer and Schlötterer 2003) and expected heterozygosity in Arlequin v3.5 (Excoffier and Lischer 2010). For microsatellite and SNP datasets we used the adegenet package (Jombart 2008) in R to calculate inbreeding coefficient (FIS) and the PopGenome (Pfeifer et al. 2014) package to calculate per site Watterson’s theta (ϴW). Where referred to, we used R statistical software v3.6.0 (R Development Core Team 2019).

We calculated pairwise FST estimates between demes (nuclear genetic clusters identified with Structure Fig. 2) for autosomal SNPs, Z chromosome SNPs, microsatellites and mtDNA datasets in Arlequin. We calculated diversity statistics separately for all sampling sites, demes, and genetic lineages.

Fig. 2
figure 2

Results of Bayesian clustering analyses of snowy plovers using three genetic datasets (1089 autosomal SNPs; 15 autosomal microsatellites and 52 Z-linked SNPs) estimated by Structure. The most likely models (K = 3 and K = 4) are shown. See Table 1 for population codes and respective genetic lineages are indicated below plots. Eastern nivosus are represented by population FLO and western nivosus by populations UTA, MT, CEU, SAQ, NAY and TEX

Full size image

Population structure

For the mtDNA D-loop sequences, we constructed a haplotype network using the Templeton, Crandall and Sing (TCS, Templeton et al. 1992) method implemented in PopArt (Leigh and Bryant 2015) to visually inspect the distribution of haplotypes across demes and genetic lineages.

To determine genetic structure with SNPs and microsatellites we performed Bayesian cluster analysis using the program Structure v2.3.4 (Pritchard et al. 2000) automated for command line usage with StrAuto (Chhatre and Emerson 2017). We tested K values between one and five for 1 million repetitions following a burn-in of 250,000 repetitions with 20 replicates for each K value using an admixture model. Due to an a priori expectation of limited population differentiation, we additionally ran the program with the LOCPRIOR model (Hubisz et al. 2009). We merged runs for each K value using Clumpak (Kopelman et al. 2015) implemented through Pophelper (Francis 2017). We determined the most biologically relevant K value by inspecting Delta K values (Evanno et al. 2005), Ln Pr (X|K) values, as recommended by Janes et al. (2017), and the bar plots of ancestry proportions grouped according to geographic distance between sampling sites and populations.

As an alternative approach to visualizing population structure, we also ran a principal components analysis (PCA) based on the same set of SNPs used in the Structure analysis. Following the approach of Novembre and Stephens (2008), we scored each genotype as having 0, 1 or 2 copies of the reference allele. In the case of missing data (~ 1.5% of the data), individuals were assigned a value equal to two times the global frequency of the alternative allele at that SNP. In the case of Z-linked SNPs, females were scored as having 0 or 1 copies of the reference allele and one allele missing, such that input genotypes for females were calculated as 0 or 1 plus the global frequency of the alternative allele at each SNP. PCAs were computed in R using the function prcomp.

Finally, we analysed the full set of variable autosomal loci in fineRADstructure (Malinsky et al. 2018), which leverages the information available in the linkage of SNPs within each RAD locus to infer patterns of recent shared ancestry. More specifically, each individual’s ancestry at a given locus is allocated to other individuals that carry either an identical haplotype (or most closely related haplotype if the focal individual has a unique allele). Rare haplotypes, which are on average of more recent origin (Kimura and Ohta 1973), make the greatest contribution to the resulting pairwise co-ancestry coefficients. Thus, the fineRADstructure analysis uses all of the variable sites at each locus and focuses on a different kind of information in the data than does Structure or PCA.

Correlation between geographic and genetic distances

To test for a correlation between geographic and genetic distances we used Mantel tests with 9999 permutations in the R package adegenet (Jombart 2008). For mtDNA sequences, we calculated genetic distances using the Maximum Composite Likelihood method implemented in MEGAx (Tamura et al. 2004; Kumar et al. 2018). We calculated individual pairwise genetic distances from microsatellites and SNPs as 1 − kinship coefficient (Loiselle et al. 1995) estimated using SPAGeDi v1.5 (Hardy and Vekemans 2002). We then created Euclidean distance matrices for geographic distances using GenALEx 6.501 (Peakall and Smouse 2012). We performed these analyses across all populations and separately within nivosus, which was sampled at multiple breeding sites (Fig. 1).

Demography

We used the microsatellite dataset to investigate past effective population size (Ne) changes using Bottleneck v1.2.02 (Piry et al. 1999) and Msvar v1.3 (Beaumont 1999; Storz and Beaumont 2002). Bottleneck tests for heterozygosity excess from the observed number of alleles under the assumption that after a recent bottleneck there is higher heterozygosity under Hardy–Weinberg equilibrium than under a mutation drift model (Cornuet and Luikart 1996). We used the Wilcoxon signed-rank test implemented in Bottleneck with a two-phase mutation model to assess the significance of the excess. Using Msvar we inferred the demographic history of each deme (western nivosus, eastern nivosus, tenuirostris and occidentalis) by comparing three models: bottleneck, expansion and stable population size based on the full likelihood, coalescent simulation method implemented in the software. Populations within western nivosus were also tested independently because of the controversial status of the Pacific coast population segment, but these analyses produced nearly indistinguishable results (Fig. S1). Therefore, we combined these populations to provide deme-wide estimates. The stable population model assumed the same Ne values for N0 (current) and N1 (ancestral). For the bottleneck, we specified the priors for N1 to be wider than at N0 (bottleneck), whereas, we assumed the opposite (priors for N0 wider than N1) for the expansion model. We set the microsatellite mutation rate (μ) to 5 × 10–4 substitutions per site per generation (Peery et al. 2012), with a range of 10–3 to 10–5. We ran each model with 1 × 1010 MCMC iterations, thinning every 100,000 steps to end with 100,000 thinned updates. The initial 20% were discarded as burn-in. We performed multiple runs with broad priors (Table S3) to influence posterior distributions as little as possible. We assessed convergence of the MCMC simulations for each model with the Gelman and Rubin’s diagnostic (Gelman and Rubin 1992; Brooks and Gelman 1998) calculated using the Boa package v1.1.7 (Smith 2007) in R.

To examine population demography with the autosomal SNPs we applied Approximate Bayesian Computation (ABC, Beaumont et al. 2002) in ABCtoolbox v2 (Wegmann et al. 2010) to compare our empirical data to simulated results from three different demographic models (stable, bottleneck and expansion). ABC was conducted for single demes independently: western nivosus (California (USA), Utah (USA) and Ceuta (Mexico) combined), eastern nivosus (Florida (USA)), occidentalis (Peru) and tenuirostris (Puerto Rico).

ABC uses the following logic (reviewed by Beaumont 2010). First, summary statistics are calculated from the empirical dataset in addition to a set of simulated datasets created based on pre-defined demographic models (involving, for example, population divergence, bottlenecks, expansions, migration). Second, summary statistics from the empirical dataset are compared with those from the simulated datasets to reveal the best fitting demographic model. Finally, differences between the empirical and best fitting model are used to predict the posterior distributions of model parameters (e.g. past/current effective population size and time of divergence; Beaumont et al. 2002).

We specified wide priors (Table S4) with uniform distributions, and modelled effective population sizes with log10 transformed values. Our dataset contained 1288 individual independent autosomal SNPs (see “Results”: SNP filtering). To simulate 1 million datasets for model optimisation we used fastsimcoal v2.6 (Excoffier and Foll 2011; Excoffier et al. 2013). We computed summary statistics (total heterozygosity, mean and standard deviation (SD) for the heterozygosity over SNP loci, number of alleles over all SNP loci, mean and s.d. for the number of alleles over SNP loci) for each dataset (simulated and observed SNPs) with the console version of Arlequin v3.5, arlsumstat (Excoffier and Lischer 2010).

To infer the demographic model providing the best fit to our observed data we calculated the Euclidean distance between observed and simulated datasets, and carried forward the closest 5000 simulations for model fitting and posterior estimations. For each model we used the general linear post-sampling adjustment step within ABCtoolbox to calculate posterior probabilities and marginal densities. We evaluated the best fitting demographic scenario by pairwise comparisons of the marginal densities and p values of each model (Wegmann et al. 2010). Spearman’s Rho statistics to assess pairwise correlation between summary statistics were calculated, and highly correlated statistics were removed before model comparisons. The uncorrelated summary statistics used in model comparisons were the total number of alleles over loci and the mean heterozygosity. We considered models with p-values greater than 0.05 as possible and the best fitting demographic scenario was required to have a Bayes Factor (BF) ≥ 3. We tested for stable, expansion, or bottleneck models within the last 1000 years as this focused our analysis on the period of intensive human induced change to coastal habitats (Lotze et al. 2006; Spatz et al. 2017). The generation time for ABCtoolbox and Msvar was set to one year based on monitoring data from natural populations (e.g. Eberhart-Phillips et al. 2017).

In addition, we generated recent migration estimates between demes using the 1288 autosomal SNP dataset with a modified version of BayesAss v3 (Wilson and Rannala 2003; available from https://github.com/smussmann82/BayesAss3-SNPs), which allows the analysis of large SNP datasets. To optimise acceptance rates, we first explored migration mixing parameter combinations (migration rates m ranging from 0.1 to 0.6, allele frequencies a ranging from 0.1 to 0.6; inbreeding coefficients f ranging from 0.1 to 0.5) with short default numbers of MCMC steps, burn-in, and sampling values as recommended by Rannala (2015). Once the acceptance rates fell between 20 and 60% (Rannala 2015) we conducted three longer runs with different seeds, with the following parameter combination: m = 0.4, a = 0.6, f = 0.4, for 10 million MCMC steps, 1 million burn-in and sampling every 500 iterations. We visualised the results and determined the convergence of the three runs using TRACER v1.7 (Rambaut et al. 2018) and averaged the results of the runs to provide migration rate estimates.

Results

SNP generation with ddRAD

We obtained a total of 62.5 million sequence reads from 47 individuals. After initial filtering and de-multiplexing, read coverage ranged from 0.46 to 1.4 million reads per individual with a mean (± SD) of 1.1 (± 0.15) million reads. From these we identified 23,026 putative RAD loci, including 4859 loci that passed quality filters (4582 putative autosomal and 277 putative Z-linked). Of these, 1338 autosomal and 65 Z-linked RAD loci included between 1–4 SNPs and/or indels. Selection of one bi-allelic SNP per locus resulted in datasets comprising 1288 autosomal, and 65 Z-linked SNPs. Excluding singletons (for Structure runs and PCA), the dataset were reduced to 1089 autosomal SNPs and 52 Z-linked SNPs.

Genetic diversity

Overall, we identified 25 unique haplotypes for a 424 bp section of the mtDNA control region (D-loop). MtDNA haplotype diversity was highest in occidentalis, followed by nivosus (eastern and western combined), and tenuirostris (Table 1). Microsatellite mean allelic richness was low and ranged from 1.8 in tenuirostris up to 2.5 in western nivosus populations (Utah and Ceuta; Table 1). Similarly, per site genetic diversity (Watterson’s theta) across all quality filtered autosomal RAD loci (n = 4582 loci) was low overall but highest in nivosus (eastern and western combined, 0.0008), followed by tenuirostris (0.00048) and occidentalis (0.00045) (Table 1). The highest diversity of Z-linked loci was found in nivosus (eastern and western combined), followed by occidentalis and lastly tenuirostris (Table 1). No substantial level of inbreeding was detected, however, comparing between populations, eastern nivosus consistently had the greatest mean inbreeding coefficient ranging from 0.02 (± 0.32) based on autosomal SNPs to 0.08 (± 0.12) based on microsatellites (Table 1).

Population structure

mtDNA

Based on mtDNA haplotype data, the three genetic lineages were not strictly distinct. We observed multiple shared haplotypes between nivosus and tenuirostris, whereas, only one haplotype was shared between occidentalis and tenuirostris (Fig. 3). On the other hand, haplotypes from occidentalis and tenuirostris appear to be more closely related (Fig. 3). Eastern nivosus had one exclusive haplotype, whereas, we detected twelve exclusive haplotypes in western nivosus, three in occidentalis, and four in tenuirostris (Fig. 3).

Fig. 3
figure 3

Haplotype network of snowy plover D-loop mtDNA sequences constructed using the statistical parsimony method (TCS). Samples were collected from localities representing all previously identified genetic lineages. Colour coding follows Fig. 1

Full size image

Microsatellites and SNPs

The Bayesian structure analysis conducted with SNP and microsatellite datasets all split occidentalis from the remaining two genetic lineages at K = 2 (data not shown). At K = 3, the SNP datasets, but not the microsatellite dataset, suggest three distinct clusters that are consistent with the updated delineation of the three subspecies proposed by Funk et al. (2007) (Fig. 2). Without location prior information, the microsatellite dataset showed high levels of estimated admixture between tenuirostris the nivosus at K = 3 (Fig. 2). However, with location prior information the three genetic lineages could be clearly discriminated with microsatellites (Fig. S2). At K = 4, Structure analysis further split the nivosus lineage into eastern nivosus (Florida) and western nivosus (all other nivosus sampling sites) for both SNP datasets (Fig. 2). Furthermore, the autosomal SNP dataset showed minimal apparent admixture between tenuirostris and eastern nivosus. Increasing K to 5 did not result in further meaningful clustering consistent with geographic information. The eastern and western nivosus distinction was most pronounced with Z-linked SNPs compared to autosomal SNPs, and the microsatellite dataset only showed weak distinction between eastern and western nivosus (Fig. 2) but improved when location information was added as a prior (Fig. S2).

The three genetic lineages were separated in the PCAs using autosomal and Z-linked SNPs (Fig. S3). Separation of eastern and western nivosus became clear when plotting PC1 against PC3 (Fig. S3B). Fine scale structuring was less pronounced within Z-linked SNP dataset (Fig. S3C).

The co-ancestry matrix of individuals using the autosomal RAD loci further supported the clear distinction of the three genetic lineages and additional sub-structuring within nivosus to recover eastern and western demes (Fig. 4). Occidentalis was consistently the most distinct lineage from Structure, PCA (Fig. S3), and fineRADstructure (Fig. 4) analyses compared to tenuirostris and nivosus.

Fig. 4
figure 4

Pairwise co-ancestry matrix and dendrogram from fineRADstructure analysis of snowy plover demes. The co-ancestry coefficients, ranging from low (yellow) to high (blue), provide a relative index of recent shared ancestry between the focal individual (along the vertical axis) and other all other individuals in the analysis (along the horizontal axis). The co-ancestry matrix was created from 1338 autosomal RAD loci

Full size image

Pairwise differentiation

Between genetic lineages we found high and significant pairwise population differentiation with all markers (Table 2). The occidentalis lineage was most differentiated from all other demes/lineages (Table 2). Pairwise differentiation between eastern and western nivosus demes showed significant differentiation (FST = 0.04–0.11) in microsatellites as well as autosomal and Z-linked SNPs, but not in mtDNA (ΦST = 0.03).

Table 2 Pairwise FST (and ΦST mtDNA) comparisons between demes of snowy plover (Charadrius nivosus), see Fig. 2 for population clustering
Full size table

Correlation between geographic and genetic distances

We found significant correlation between geographic and genetic distances for autosomal and Z chromosome datasets across all lineages and within nivosus individuals only (all lineages/nivosus only: autosomal SNPs r (mantel test statistic) = 0.70/0.61, Z chromosome SNPs r = 0.75/0.45, microsatellites r = 0.38/0.11, all p < 0.05; Fig. 5). For mtDNA, the correlation between geographic and genetic distances was only significant for the dataset that included all genetic lineages (all lineages r = 0.52, p < 0.05, Fig. 5a; nivosus only r = -0.02, p = 0.73, Fig. 5b). For nuclear markers, the correlation between geographic and genetic distances was weakest for microsatellites and strongest for Z chromosome SNPs across all subspecies (Fig. 5a). Within nivosus only, we detected the strongest correlation between geographic and genetic distances with autosomal SNPs, followed by Z chromosome SNPs, and a weak but significant correlation with microsatellites (Fig. 5b).

Fig. 5
figure 5

Pairwise individual-based correlation between geographic and genetic distances across snowy plover breeding populations (a) across all lineages, and (b) within C. n. nivosus only, assessed with four genetic datasets. Regression lines are provided where the correlation is significant (p < 0.05)

Full size image

Demography

Overall our demographic analyses provided strong evidence for population bottlenecks in all genetic lineages during the past 1000 years. Specifically, Msvar predicted bottlenecks occurring between 56 and 445 years ago (YA) when the effective population size decreased from an average of 5983 (± 581) to 33 (± 19) (Table 3). We found the largest effective population decline in western nivosus and tenuirostris and the least decline in occidentalis (Fig. 6 and Table 3; Table S5). According to the ABC analysis, the population bottleneck model fitted our autosomal SNP datasets best for eastern nivosus, tenuirostris and occidentalis with Bayes Factors of > 5. However, none of the three tested scenarios received sufficient support in western nivosus (Table S6). ABC estimated that the bottlenecks occurred most recently in tenuirostris and occidentalis (221 YA and 302 YA respectively) than in eastern nivosus (784 YA) (Table 3). Posterior and prior distributions of the ABCtoolbox bottleneck scenario are presented in Supplementary Material (Fig. S4). Only tenuirostris showed statistical support for a bottleneck based on heterozygosity excess with the program Bottleneck (Wilcoxon test: p = 0.04).

Table 3 Effective population size (Ne) and time of bottleneck of snowy plover (Charadrius nivosus) demes estimated with Msvar and ABCtoolbox (shaded cells) with 95% highest posterior densities (HPD)
Full size table
Fig. 6
figure 6

Effective population size (current and ancestral) of Eastern nivosus, tenuirostris, and occidentalis snowy plover genetic lineages estimated by Msvar and ABCtoolbox. Density is normalised for comparison of the methods. ABCtoolbox posteriors are restricted by prior parameter values. Effective population size is reported on log10 scale. Western nivosus did not fit any of the modelled demographic scenarios tested with ABCtoolbox and is therefore not shown

Full size image

We found low but significant levels of recent migration (m = 0.02–0.03) between lineages in the BayesAss analyses (Table S7). Among demes, there was one notable exception of high migration (m = 0.26, s.d. 0.04) indicating gene flow from western to eastern nivosus. However, given the low sample sizes and the weak population differentiation between eastern and western nivosusFST ~ 0.05, these results should be interpreted with caution (Meirmans 2014).

Discussion

Examining population structure based on a comprehensive assessment that included four types of genetic data, we find genetic support for the three genetic lineages delineated previously as subspecies (Funk et al. 2007). These genetic lineages largely correspond to the proposed geographic ranges of the subspecies: nivosus (Mexico and USA), tenuirostris (Bermuda and Caribbean), and occidentalis (South America). The level of genetic divergence among the snowy plover genetic lineages/demes in microsatellite (FST = 0.04–0.47), autosomal SNPs (FST = 0.05–0.58), Z-linked SNPs (FST = 0.12–0.71) and mtDNA data (ΦST = 0.03–0.79) is similar to that observed among other plover subspecies, such as the piping plover (Charadrius melodus, Miller et al. 2010) and the chestnut banded plover (Charadrius pallidus dos Remedios et al. 2017). Genetically, the three snowy plover lineages are more differentiated than the subspecies of some other shorebirds such as dunlin, Calidris alpina, Miller et al. 2015) or redshank, Tringa tetanus, Ottvall et al. 2005). The differences at microsatellites and mtDNA are also more pronounced than between white-faced plover (Charadrius alexandrinus dealbatus) and Kentish plover (Charadrius alexandrinus), FST = 0.01, Rheindt et al. 2011). The latter two taxa are considered separate species by some authorities based on phenotypic differences (del Hoyo et al. 2018; but see Clements et al. 2018; Gill and Donsker 2018) and recent analyses have described further genomic differences between both taxa (Sadanandan et al. 2019; Wang et al. 2019). We found that occidentalis is the most distinct lineage of snowy plovers with consistently the greatest level of differentiation across all analyses. We suggest the uniqueness of this lineage should be considered when prioritising conservation efforts for South American shorebirds.

By combining new data with those from two studies (Funk et al. 2007; D’Urban Jackson et al. 2017) we have substantially increased the number of genetic markers, individuals and sampling locations assessed. The SNP results corroborate previous findings (Funk et al. 2007) that tenuirostris lineage is distinct from nivosus and should be considered as a separate conservation unit. At the same time, individuals sampled from Florida (eastern nivosus) are more closely related to other nivosus populations, rather than the earlier delineation which grouped it with Caribbean tenuirostris (American Ornithologist’s Union 1957). This finding supports the current state-wide protection of Florida snowy plovers as a separate management unit. Avian subspecies are defined by phenotypic differences between geographically distinct populations. So far, such differences between Caribbean tenuirostris and eastern nivosus have not been reported, preventing the delineation of a separate subspecies for either tenuirostris or eastern nivosus. The exact geographic ranges of the two nivosus demes (eastern and western) are currently unknown. Our study lacked samples from existing intermediate breeding populations such as Kansas, Texas, Alabama and Lousiana, which are needed to fully characterise the boundaries of eastern and western nivosus demes. Previously, the American Ornithologist’s Union assigned all mainland snowy plovers east of Louisiana to tenuirostris and all western snowy plovers to nivosus. This delineation may serve as a useful primer for further work that aims to characterize the geographic ranges of the demes in detail. Our estimates of recent migration suggest that eastern snowy plover populations receive gene flow from western populations, whereas gene flow in the opposite direction is minimal. This directional bias maybe the result of wintering western migrants occasionally remaining to breed within the breeding distribution of the eastern snowy plover.

Snowy plovers have been the subject of several movement and demographic studies, though none so far have encompassed the entire range (Warriner et al. 1986; Stenzel et al. 1994, 2007, 2011; Eberhart-Phillips and Colwell 2014; Aiello-Lammens and Akcakaya 2016; Cruz-López et al. 2017; Eberhart-Phillips et al. 2017). Therefore, the long distance dispersal behaviour of this species remains poorly described and anecdotal. We currently lack sighting data to support or reject the genetic inference of unidirectional gene flow from western nivosus populations in Florida. Such gene flow from larger western to smaller eastern nivosus populations may provide a crucial bolstering effect against population declines in the more vulnerable eastern deme (Lamonte et al. 2002).

Despite the significant differentiation between the eastern and western nivosus demes, we observed only a weak correlation between geographic and genetic distances within the entire nivosus subspecies. At first, this appears surprising given the large distance between sampled sites (> 3000 km). However, weak clinal structure or genetic homogeneity is common in shorebird populations, including Temminck’s stint Calidris temminckii (Rönkä et al. 2012), ruff, Philomachus pugnax (Verkuil et al. 2012), mountain plover Charadrius montanus (Oyler-McCance et al. 2008) and the closely related Kentish plover, C. alexandrinus (Küpper et al. 2012). Interestingly, across the entire range of snowy plovers, we find that Z-linked markers show the strongest correlation between geographic and genetic distances compared to autosomal SNPs, microsatellites, and maternally inherited mtDNA. This result is consistent with male philopatry and female biased dispersal. Female biased dispersal is generally common across avian families regardless of mating system (Greenwood 1980; Clarke et al. 1997). Field observations suggest strong female biased dispersal in snowy plovers (Stenzel et al. 1994, 2007; Paton and Edwards 1996; Colwell et al. 2007; Pearson and Colwell 2014), a pattern that may be explained by the “dispersal-to-mate” hypothesis. The “dispersal-to-mate’ hypothesis predicts that polygamous females disperse to access more mates and consequently achieve higher breeding success (Küpper et al. 2012; D’Urban Jackson et al. 2017).

Autosomal polymorphic SNPs showed strong support for four conservation units in population structure analyses (Figs. 2, 4, S3). This highlights the advantages of using as many independent genetic markers as possible when investigating potentially controversial fine scale population structure, especially in high gene flow species such as the snowy plover (e.g. Ruegg et al. 2014; Vendrami et al. 2017). Differences in effective population size between mitochondrial, autosomal and sex-linked loci might also contribute to differences in level of genetic structure observed between genetic data sets. Assuming similar variance in male and female reproductive success, autosomal loci have the highest effective population size, followed by Z-linked loci and then mtDNA. Relatively high Z chromosome divergence is commonly observed in birds (Irwin 2018), which may be enhanced by “Fast-Z” evolution (Mank et al. 2007; Meisel and Connallon 2013) in addition to differences in effective population size. In snowy plovers, Z chromosome differentiation based on FST values was greater than for autosomal SNPs. The ability of a relatively small number of Z-linked SNPs to distinguish between eastern and western snowy plovers in Bayesian cluster analysis was higher than that of the autosomal SNPs, although the Z-linked SNP PCA indicated weak distinction between the two nivosus demes. Furthermore, Z-linked SNPs and autosomal SNPs produced similar correlations between geographic and genetic distances. The similarity of general spatial patterns produced by the two SNP datasets is striking given the difference in number of polymorphic markers (Z-linked: n = 52/65 excluding/including singletons; autosomal: n = 1089/1288 excluding/including singletons).

Consistent with previous findings (Funk et al. 2007), we find no evidence of genetic differentiation between coastal and interior western nivosus populations and no differentiation among coastal Pacific populations. This is not surprising given the high dispersal observed in this species (Stenzel et al. 1994). To detect putative fine scale differentiation within this deme, substantially greater sampling of the genome might be required, with the most fruitful approach being an analysis of whole genome sequences (WGS). Combining a WGS approach with improved sampling strategy to encompass the entire range of the species would enable not only the delineation of deme boundaries but also allow for an investigation of adaptive variation (Kjeldsen et al. 2016).

Also supporting previous studies, we found that the current effective population sizes of snowy plover subspecies and demes are low (Funk et al. 2007). Our demographic models suggest that this is most likely the result of strong population bottlenecks in all genetic lineages/demes that have occurred within the last 1000 years. As a result of these bottlenecks, the effective population sizes of the genetic lineages/demes have been reduced by at least 98%. Overall there was remarkable similarity between the full likelihood (Msvar, microsatellites) and ABC (ABCtoolbox, autosomal SNPs) methods, with the exception of a lack of fit between any of the tested ABC demographic scenarios and data for the western nivosus deme. This may have been caused by conflicting demographic trends across the vast breeding range. The western nivosus range spans large latitudinal and longitudinal gradients (more than 30 and 40 degrees, respectively), and encompasses populations that have likely experienced episodes of both growth and contraction. For example, western nivosus populations at the Pacific coast declined in the 1980s (Page et al. 1991) but due to intensive conservation effort they have stabilised or even increased over the last 20 years based on census results (Mullin et al. 2010; Colwell et al. 2017; Feucht et al. 2017). In contrast, Pacific populations in Mexico appear to be rapidly declining (Galindo-Espinosa and Palacios 2015; Cruz-López et al. 2017).

Both ABC and Msvar models inferred effective population sizes far smaller than the census sizes of breeding snowy plover populations (Thomas et al. 2012). For example, the census size of eastern snowy plovers has been estimated between 683 and 1022 (Thomas et al. 2012), whereas, our effective population size estimate was between 47 (Msvar) and 59 (ABC) albeit with wide confidence intervals (CI) of < 1–1455 (Msvar) and 10–5941 (ABC). Even more surprisingly, the western nivosus population size was estimated at 17 (Msvar, CI: < 1–576), whereas, the census size is over 8000 (Thomas et al. 2012). These discrepancies may either have technical or biological explanations. Shafer et al. (2015) found that a larger number of genetic markers/loci yielded more power in estimating the parameters of demographic models using ABC, with over 50,000 loci being needed for accurate effective population size estimations. The accuracy of ABC compared to other methods in estimating effective population size is largely unknown (Wang et al. 2016), and ascertainment bias in SNP datasets can lead to erroneous temporal estimates (Quinto-Cortés et al. 2018). Furthermore, the large confidence intervals of all the estimated parameters should be considered when interpreting these results. Snowy plover reproductive behaviour and population ecology may further contribute to the discrepancy between effective and census population sizes and reduce genetic diversity. These populations are highly male-biased (Carmona-Isunza et al. 2017; Eberhart-Phillips et al. 2017), deviating from the basic model of unbiased sex ratios. The polyandrous mating system will further reduce the effective population size in comparison to the census population (Nunney 1993), especially in the nivosus subspecies, which shows the highest rates of observed polygamy (Warriner et al. 1986; Eberhart-Phillips et al. 2017).

Overall, our results indicate a recent and substantial reduction in population size across snowy plover populations. Two of the demographic estimation methods we implemented (Msvar and Bottleneck) indicated that population size reduction has been pronounced in tenuirostris. This finding is supported by current population extirpations of tenuirostris on Caribbean islands (Brown 2012). Bottleneck did not find evidence of population reduction in occidentalis or nivosus subspecies, which could be because the program is best suited for detecting low magnitude population bottlenecks (Peery et al. 2012). Again, the power to detect bottlenecks depends on the number of loci analysed (Hoban et al. 2013).

During the past 1000 years, humans have dominated and profoundly altered global ecosystems (Vitousek et al. 1997; Barnosky et al. 2011; Ceballos et al. 2015) especially in coastal and island environments (Lotze et al. 2006; Spatz et al. 2017). The introduction of predators such as the non-native red fox (Vulpes vulpes) and feral cats (Felis catus) plus other forms of human disturbance in coastal environments have led to reduced snowy plover breeding productivity and survival (Lafferty 2001; Ruhlen et al. 2003; Colwell et al. 2007; Webber et al. 2013; Dinsmore et al. 2017). Therefore, we hypothesise that reductions in snowy plover population size may have resulted from the combined impacts of increasing anthropogenic modification of habitats. Estimating demographic histories is becoming more common in conservation genetic studies (Carnaval et al. 2009; Palsbøll et al. 2013; Shafer et al. 2015; Stoffel et al. 2018) and can be used among other to prioritize conservation units (Stockwell et al. 2013).

Implications for snowy plovers and future work

Evidence of low and declining effective population size, combined with continued threats to snowy plovers highlight the vulnerable status of this species. Our understanding of this species is strongly reliant on studies from the nivosus subspecies, with little dedicated monitoring programs for tenuirostris and occidentalis populations. The demographic findings we present here and recent population extirpations occurring in tenuirostris (Brown 2012), provide strong justification for increased efforts at the Caribbean subspecies to halt the loss of genetic diversity and adaptive potential (Willi et al. 2006). In addition, basic monitoring of occidentalis breeding populations, such as those conducted by Wetlands International (Wetlands International 2019), are required to obtain reliable information on population trends.

Due to the high vagility and potential for gene flow among snowy plover populations, future work may benefit from population genomics based on whole genome sequence information with comprehensive sampling to focus on adaptive diversity in combination with actual movement data (e.g. GPS tagging) and range-wide demographic studies. This would allow us to fully characterise population connectivity and thus, infer conservation management units (Crandall et al. 2000; Fraser and Bernatchez 2001; Palsbøll et al. 2007; Lowe and Allendorf 2010; Funk et al. 2012; Flanagan et al. 2018). Our results suggest that eastern snowy plovers deserve a continuation of distinct conservation management, whereas we caution against immediate changes in management actions of the western nivosus deme such as delisting of the threatened Pacific distinct population segment. The high vagility and challenges provided by the unusual mating system in this species call for further work that should be based on an integrative approach using genomic, demographic, and movement, data to model population connectivity and viability.