Seaside Sparrows (Ammospiza maritima) are tidal salt marsh endemic passerines found along the Atlantic and Gulf coasts of North America. Currently there are 7 described subspecies and "MacGillivray's" Seaside Sparrow (A. m. macgillivraii) is the name given to the Atlantic coast subspecies breeding from North Carolina to northern Florida. In 2019 the US Fish and Wildlife Service received a petition to list this subspecies under the Endangered Species Act due to shrinking populations and loss of marsh habitat, which necessitated a Species Status Assessment. As part of the Species Status Assessment, we analyzed genetic (microsatellite and mitochondria) data from 464 Seaside Sparrows from Connecticut to Florida, USA, to infer population connectivity (gene flow) among Atlantic coast populations, and to assess the concordance of population genetic data with the putative ranges of the subspecies. Bayesian cluster analysis (program Structure) indicates three genetically distinct population segments that were recovered: (1) Florida + Georgia, (2) Charleston, South Carolina, and (3) North Carolina to Connecticut. The population in Charleston, South Carolina was the most strongly differentiated based on microsatellite data, and harbored a unique mitochondrial haplotype not shared by other sampling locations, possibly reflecting long-standing isolation. Our results indicate population genetic discordance with currently described ranges of the subspecies of Seaside Sparrow and provide grounds for the consideration of separate management plans for the three populations.

Sampling

We analyzed data from Seaside Sparrows (N = 464) captured in mist nets at 14 sites along the Atlantic coast (Figure 1, Table 1). We banded each bird with a uniquely numbered USGS aluminum band and collected blood samples from the brachial vein which we stored in Queen’s lysis buffer (Seutin et al. 1991). We chose sampling localities throughout the entire range of Seaside Sparrows on the Atlantic coast, with more intensive sampling across the currently described boundary (Dare County, North Carolina) between A. m. maritima and A. m. macgillivraii, as well as the type locality for MacGillivray’s Seaside Sparrow (Charleston, South Carolina) (see Figure 1, Table 1). There are no diagnostic plumage criteria to distinguish the two subspecies, so we used only samples of putative macgillivraii taken between 20 May and 30 July to avoid including wintering maritima. All putative macgillivraii sampled in 2017-2018 were either adults in breeding condition (based on the presence of either a cloacal protuberance or brood patch) or birds hatched in the years of capture, indicating that these were birds from locally breeding populations.

Laboratory Methods

Microsatellites. We extracted DNA using Qiagen DNeasy or ThermoScientific Tissue and Blood kits, and quantified DNA concentration using a Nanodrop 2000 (ThermoScientific, Delaware, USA). We amplified the following 15 microsatellite loci for each individual (N = 464) using polymerase chain reaction: Aca01, Aca11 (Hill et al. 2008), Am02, Am12, Am14, Am18, Am20, Am32 (Lehmicke et al. 2012), Sosp13, Sosp14 (Sardell et al. 2010), ZoleC06, ZoleC11, ZoleE11, ZoleF11, and ZoleG03 (Poesel et al. 2009). We modified all primers with a 19 bp M13 tag (Boutin-Ganache et al. 2001). We ran polymerase chain reactions (PCR) in 10 µl reactions with 3.00 mM MgCl₂, 0.16 mM dNTPs, 1X buffer, 0.06 µM forward primer, 0.36 µM reverse primer, 0.60 µM dye-labeled m13 primer, 0.10 units Taq polymerase, and 10-20 ng DNA. We improved amplification of some loci by adding dimethyl sulfoxide (DMSO; 0.30 μL): Am32, Am20, Aca01, Am18, Sosp13, Sosp14, ZoleH02, ZoleC11, ZoleC06, ZoleE11, ZoleF11; and/or betaine (1.00 μL): Aca01, Am12, Am14, Am18, Am20, Am32, Sosp13, Sosp14, ZoleC06, ZoleC11, ZoleE11. Thermocycler protocols were as follows: 94˚C for 1 minute, 33 (Am14, Sosp14, ZC06, ZG03) or 35 (Aca01, Aca11, Am02, Am12, Am18, Am20, Am32, ZC11, ZE11, ZF11, ZG03) cycles of 94 ˚C for 30 seconds, an annealing temperature of 60 ˚C for most loci (55 ˚C for Aca01, Aca11, Am32, and ZG03, and 57 ˚C for Am12 and Am20) for 30 seconds, 72 ˚C for 30 seconds, and a final extension of 72 ˚C for 5 minutes. We pooled PCR products (post-PCR) into four panels (dyes: PET, VIC, NED, and 6-FAM) and ran them on an Applied Biosystems 3730 capillary sequencer. We used GeneMarker (v. 2.7.0; SoftGenetics, LLC.; State College, Pennsylvania) for genotyping.

Mitochondrial DNA. We sequenced the protein coding NADH dehydrogenase subunit 2 (ND2; 1042bp) from 9-20 samples from each site (Table 1). We ran PCRs in 25 μL volumes with 1.5 mM MgCl₂, 0.8 mM dNTPs, 1.25 μM forward primer (L5215; Hackett 1996), 1.25 μM reverse primer (H6313; Johnson and Sorenson 1998), 1X buffer, 2.5 units of Taq polymerase, and ~10 ng DNA. The thermocycling protocol was as follows: 94˚C for 30 s, 34 cycles of 94 ˚C for 30 s, 52˚C for 30 s, and 72˚C for 1 min, with a final extension of 72˚C for 7 minutes. Eurofins Genomics (Louisville, KY, USA) performed SimpleSeq Sanger DNA sequencing and we aligned sequences using Sequencher (v. 5.4.6, Genecodes Corporation, Ann Arbor, Michigan).

Data Analysis

Summary Statistics and relatedness. We used MicroChecker (van Oosterhout et al. 2004) to test for evidence of scoring issues (e.g., stutter, large allele dropout) and null alleles within the entire dataset. To test for evidence of within-population deviation from Hardy Weinberg Equilibrium (HWE) and Linkage Equilibrium (LE) for each sampling location we used GenePop v. 4.2 (Rousset 2008). The presence of close relatives (full-siblings and parent-offspring) can bias population genetic analyses (Sánchez-Montes et al. 2017) by overestimating population structure when using unsupervised Bayesian clustering algorithms and by causing deviations from HWE and LD (Rodríguez-Ramilo and Wang 2012). We identified 35 dyads of close relatives within our dataset using coancestry (Wang 2010), and one of each pair of full-siblings or parent-offspring (r_xy, estimate of 0.5 or higher across all relatedness estimators) were removed prior to downstream analyses, resulting in a sample size of 429 individuals. The intention of this cut-off was to remove full-sibling and parent-offspring pairs from the dataset; however, relatedness estimates can be imprecise (Taylor 2015), and shifting the threshold from a relatedness estimate of 0.5 to 0.4 yielded an addtional 3 dyads, resulting in a sample size of 426 individuals. Most dyads, regardless of the threshold chosen, were restricted to our northern-most sampling locations (sites 9 and 10). These samples were collected from demographic sites where netting occurred frequently throughout the season, and thus the likelihood of sampling relatives (i.e. young birds that have not completed natal dispersal) is higher than the sites in the south that we visited once or twice. We reran all analyses with this further reduced dataset and our findings did not change. With both datasets we identified the same number of genetic clusters, saw the same patterns of pairwise differentiation, and found no deviations from HWE or LD suggesting that relatives were not biasing our estimates (Rodriguez-Ramilo and Wang 2012); we chose to present the results using the more conservative dataset (N = 426, r_xy, estimates < 0.4).

Assessing patterns of genetic structure. To quantify pairwise population differentiation, we used the diffCalc function in the diveRsity package (Keenan et al. 2013) in program R (R Core Team 2019). To determine significance for pairwise comparisons that differ from zero (P < 0.001), and to account for multiple comparisons, we used the resampling approach built into the diffCalc function (Keenan et al. 2013). The function calculates confidence intervals with a bias-corrected bootstrapping method, then tests for significance using permutation (Keenan et al. 2013). We used 20,000 permutations of individuals within samples to generate our sampling distribution and chose a conservative significance threshold of P < 0.001 to assure 99.9% confidence that the comparative test statistic was different than zero. We report the results of two methods for calculating differentiation: G_ST (Nei and Chesser 1983) and G´_ST (Hedrick 2005). We report G_ST for familiarity and ease of comparison with other literature (Alcala and Rosenberg 2019), but suggest caution in its interpretation because we are comparing levels of differentiation at multiple alleles that exhibit high heterozygosity within our populations, which can constrain estimates of G_ST (Hedrick 2005). We therefore also report G´_ST, as this standardized measure of differentiation is more robust in situations of high heterozygosity (Hedrick 2005).

To explore the distribution of genetic structure throughout our sampling range we used a Bayesian clustering method in the program Structure (v 2.3.4; Pritchard et al. 2000). We ran admixture models with correlated allele frequencies using sampling locations as priors (loc prior) to test for genetic clusters (K) (Pritchard et al. 2000, Hubisz et al. 2009). Run lengths ranged between 5 x 10⁵ – 1 x 10⁶ MCMC iterations after a burn in of 5 x 10⁴ – 1 x 10⁵ iterations. We used loc prior and correlated allele frequencies to account for the non-independence of the data based upon geographic location. To ensure these priors were not overestimating signals of structure we also explored several variations of the model, including loc prior vs no loc prior, and correlated vs non-correlated allele frequencies. Initial runs explored up to 16 possible populations (K) that were subsequently reduced to 6 K in the final runs. We first processed Structure outputs using Structure Harvester (Earl and vonHoldt 2012) and determined the best support for K number of clusters based upon mean estimated log-normal probabilities and the Evanno method (Evanno et al. 2005). We used clumpp (Jakobsson and Rosenberg 2007) followed by distruct (Rosenberg 2004) to visualize Structure results. To confirm Structure findings with multivariate statistics, which do not rely on traditional population genetic assumptions like HWE or LE, we conducted discriminant analysis of principal components (DAPC; Jombart et al. 2010) with the package adegenet v2.1.2  (Jombart 2008) using the dapc function in R (R Core Team 2019); missing data (<1% of total genotypes) were replaced with mean allele frequencies (Jombart 2008).

Population-level diversity. We estimated descriptive population genetic characteristics (mean number of alleles, allelic richness, and observed and expected heterozygosity) for each sampling location with GenAlEx (v 6.503, Peakall and Smouse 2012). Given the variation in temporal sampling within our dataset we used several methods to estimate effective population size (N_e). For populations visited within a single year we used the sibship assignment method, full-likelihood option, in Colony (v 2.0.6.3, Jones and Wang 2009) and the Linkage-disequilibrium (LD) method in NEestimator (v 2.1, Do et al 2014). We used the full dataset (N = 464) before removing the close relatives (full-siblings and parent-offspring pairs) across methods because Colony uses estimated levels of relatedness (probabilities of individuals being full-siblings, half-siblings, or completely unrelated) to estimate N_e (Jones and Wang 2009, Ackerman et al. 2017) and because an increase in upward bias is correlated with aggressive sibling purging when using the LD method (Waples and Anderson 2017). Because of the aforementioned uncertainty in relatedness estimates we did not use our relatives identified in coancestry to inform the sibship method in Colony and instead allowed the method to determine those relationships independently; using a single run for each sampling locality. Our reasoning to employ two methods of estimation for single year visits is because both the sibship and LD method were developed for use with discrete, non-overlapping generations (Jones and Wang 2010, Waples et al. 2013). When applied to iteroparous species with overlapping generations both methods estimate the number of breeders (N_b) that produced the randomly sampled cohort which can be biased by ± 14% (Jones and Wang 2010, Waples et al. 2013). Fortunately, the downward bias of N_b estimates from single samples of overlapping generations, calculated with the LD method, can be reduced by applying corrections based on life history data (adult lifespan and age of maturity; Waples et al. 2013). To that end, we generated  estimates of N_b using Colony and NeEstimator to compare estimates by method and subsequently applied two-trait life history corrections from table 3 of Waples et al. (2014) to the LD estimates to generate corrected-N_b (N_bAdj2) and corrected-N_e (N_eAdj2).

For multi-year sampling efforts, we calculated temporal estimates of N_e for all sites with data that extended beyond a single generation length (Sites 3b, 4a, and 10; Table 1, Table 2) using the Pollak method (Pollak 1983) in the program NeEstimator (v 2.1, Do et al 2014). This method is thought to be particularly robust when sampling efforts are separated by multiple generations when the appropriate demographic data and sampling design for the Jorde-Ryman modification are not available (Waples and Yokota 2007). However, temporal estimates have the most power when sample sizes are large and are susceptible to bias (Waples and Yokota 2007) associated with small sample sizes; all of our sample sizes, by year, are ≤ 40 individuals so we also calculated point estimates (as described above) for each sampling year. We chose to calculate separate generation lengths (the mean age of parents of newborn cohorts) based upon migratory strategy as demographic studies of Seaside Sparrows suggest that annual survival rate is lower for migratory populations than sedentary ones and that the mean number of offspring produced annually differs between northern and southern populations (Shriver and Gibbs 2004, Roberts et al. 2017, Post and Greenlaw 2020, U.S. Fish and Wildlife Service 2018). We calculated generation length (T) using the formula in Waples and Yakota (2007) - the sum of survivorship at each age times the mean number of offspring produced per time interval by individuals in each age group (Waples and Yokota 2007). Estimates of T for sedentary populations (all populations south of Dare County, NC based on current range delineations; AOU 1957) were calculated with life span and survivorship estimates from mark-recapture data from a population in South Carolina (Post and Greenlaw 2020, Post Unpublished Data); fecundity estimates were based upon data collected for the Species Status Assessment in South Carolina, Florida, and Georgia (U.S. Fish and Wildlife Service 2018). Estimates of T for migratory populations were based upon average estimates of Seaside Sparrow survival and fecundity over 4 and 5-year periods collected in New York and Massachusetts (Post and Greenlaw 1982, Post and Greenlaw 2020); though not used in our calculation of T survival estimates from New Jersey were similar to those of New York and Massachusetts and confirmed survival rates trend higher in migratory versus non-migratory Seaside Sparrows (Roberts et al. 2017).

We used the adegenet package (Jombart 2008, Jombart and Ahmed 2011) and ade4 package (v1.7.15; Dray and Dufour 2007) in R (R Core Team, 2019) to estimate the degree of isolation by distance (IBD) across all sampling localities by implementing a Mantel test between matrices of Nei’s genetic distance (Nei 1972, 1978) and straight-line geographic distance (km) using 10,000 permutations of the data. We used the program Bottleneck (V 1.2.02; Luikart and Cornuet 1999) to test for evidence of recent bottlenecks in each sampling locality. Drastic population reductions, or bottlenecks, are accompanied by similar reductions in allelic frequencies, but this reduction occurs faster than a reduction in population heterozygosity causing a quantifiable deviation from mutation-drift equilibrium (MDE; Luikart and Cornuet 1999). Identifying deviations from MDE, and the subsequent allelic frequency changes, is important for several reasons. First, equilibrium is a required assumption for many population genetic analyses and testing for deviations gives one confidence, or caution, in the interpretation of results (Luikart and Cornuet 1999). Second, as one of the goals of this study is to identify potential areas of conservation concern, looking for evidence of bottlenecks potentially allows us to identify populations with recent reductions of effective population size (N_e) that may be more susceptible to the effects of genetic drift. 

Mitochondrial analyses. Mitochondrial ND2 sequences were analyzed in DnaSP (v 6, Rozas et al. 2003) and 95% confidence interval haplotype networks were made with TCS (v 1.21, Clement et al. 2000).We conducted an analysis of molecular variance (AMOVA) to hierarchically test for population differentiation, with stratification levels set as the full dataset and sampling locations, in our mitochondrial dataset using poppr v2.8.5 (Kamvar et al. 2014). To test for significance, we implemented a randomization test with the ade4 v1.7.15 package using the randtest function in R (R Core Team 2019) with 1000 iterations (Dray and Dufour 2007).