Piping plovers demonstrate regional differences in nesting habitat selection patterns along the U.S. Atlantic coast

,


INTRODUCTION
Habitat selection studies that encompass all or a significant portion of a species' geographic distribution are valuable for species with broad geographic ranges, such as Red-cockaded Woodpeckers (Picoides borealis) (McKellar et al. 2014), Northern Spotted Owls (Strix occidentalis) (Noon and McKelvey 1996), and sea turtles (Caretta sp.) (Liles et al. 2015).Selection patterns across a broad distribution can allow managers to identify resources that are consistent and important universally across all sites versus resources that are only used locally by specific populations (Flesch and Steidl 2010).In the United States, these differences may have important implications for protection and restoration of habitats implemented through recovery plans and regulatory mechanisms, such as section 7 of the U.S. Endangered Species Act.Such patterns can also allow managers to identify environmental thresholds that explain a species' distribution or abundance, which is particularly relevant for protecting species in the face of global climate change (Flesch and Steidl 2010).
Knowledge of differences in habitat use within a species' range is also important for crafting effective management plans.A one-size-fits-all management approach developed from selection patterns observed for only a few populations may result in plans that ignore vital resources for some populations or set unrealistic goals in others (Oliver et al. 2009, McKellar et al. 2014).In addition, management plans informed by habitat-use patterns known from only a portion of a species' range can ignore patterns at range boundaries and neglect the species' complete niche (Holt and Keitt 2005).Failure to identify and protect key habitat resources throughout the species' range, including along the boundaries, can reduce important opportunities for conserving imperiled species (Channell and Lomolino 2000).
In this study, we investigate regional habitat selection patterns by Piping Plovers (Charadrius melodus), with the objective of informing habitat management and restoration for the species.The Piping Plover is a migratory shorebird that nests along the Atlantic coast, northern Great Plains, and Great Lakes of Canada and the United States.The Atlantic coast breeding population, which is the focus of this study, nests along approximately 2000 km from Newfoundland, Canada, to North Carolina, USA.Most habitat selection studies for this species to date have focused on localized study areas in relatively close proximity (e.g., Loegering and Fraser 1995, Jones 1997, Cohen et al. 2008, 2009, Maslo et al. 2011, 2012).Understanding potential differences in habitat requirements within the Piping Plover's Atlantic breeding range is a long-standing research need with important management implications (USFWS 1996), including tailoring habitat protection and restoration actions to regional habitat selection patterns to increase their efficacy.However, because previous studies have used different methods or metrics and have focused on selection in local study areas, it has been difficult to establish whether Piping Plovers perceive and utilize habitat consistently or whether availability and selection of certain resources differ throughout the range.In this study, we identify landscape characteristics at Piping Plover nests along the Atlantic coast from Massachusetts to North Carolina, USA (Fig. 1), to explore habitat selection patterns.This study is the first to investigate patterns over a significant portion of the species' Atlantic coast breeding range, and we explicitly consider whether habitat selection patterns differ among geographic regions.Results discussed here are important for filling knowledge gaps present at the time the Piping Plover's original recovery plan (USFWS 1996) was written.This study will be important for providing consistent rangewide recommendations where possible, tailoring local site management plans to regional patterns of habitat selection, and informing how limited management resources are used.

Habitat use and availability data
To collect habitat-use data throughout the species' U.S. Atlantic coast breeding range, we collaborated with conservation practitioners and biologists from federal and state government agencies, universities, and private conservation organizations through the smartphone application iPlover (Sturdivant et al. 2016, 2018, Thieler et al. 2016, Zeigler et al. 2017).iPlover data used in this study contained landscape characteristics at locations where Piping Plovers established nests and at random points on 21 beaches and barrier islands (or sites) from Massachusetts to North Carolina during the 2014 and 2015 breeding seasons.We primarily collected data from federally or state-protected properties and surrounding areas, including national wildlife refuges (NWRs) and national seashores (NSs), although five sites were in slightly more urban environments (i.e., Rhode Island NWR complex, Fire Island, Rockaway Peninsula, Long Beach Island, and Cape Hatteras; Fig. 1).Study areas were selected because Piping Plover sub-populations were intensively monitored by trained professionals and because lidar and orthoimagery datasets were available for those locations for dates concurrent with Piping Plover monitoring.
iPlover users characterized landscape characteristics in the area immediately surrounding any Piping Plover nest found during the course of monitoring efforts.This procedure was repeated at the coordinates of random points, which we disseminated to iPlover users at the start of each breeding season.Random points could fall anywhere within the boundaries of the study area (e.g., the administrative boundary of Assateague Island NS), were not limited to sandy beach habitats, and were considered independent of each other and nest points.In general, the number of random points visited each breeding season for each study area was equivalent to the number of Piping Plover nests observed in that study area in the prior breeding season, and we attempted to create a dataset with roughly equal numbers Fig. 1.Study extent for the regional analysis of Piping Plover habitat selection patterns, including site locations where landscape characteristics for used and available resource units were observed.Site names are listed according to the name of the barrier island or beach, and many sites contained federally protected national wildlife refuges (NWR), seashores (NS), and recreation areas (NRA; listed in parentheses under site name).Under the U.S. Endangered Species Act, the Atlantic coast Piping Plover population is managed in three recovery units: the (A) New England, (B) New York-New Jersey, and (C) Southern recovery units.
of nest and random point observations.We follow the nomenclature of Lele et al. (2013), where an individual nest or random point is a resource unit that encompasses a 5 × 5 m area with the point at its center.We considered a resource unit to be used if it was associated with a Piping Plover nest and unused/available if associated with a random point.
We did not observe a Piping Plover nest at any random point in the same breeding season despite intensive species monitoring throughout the breeding season.A previous analysis of a model using this dataset (Zeigler et al. 2017) found a substantially higher number of false positives (i.e., habitat was predicted to be suitable but no nest was present) compared to false negatives (i.e., landcover was predicted to be unsuitable but a nest was present).This suggests that some random points within the iPlover dataset may have characteristics similar to those selected for by nesting Piping Plovers.Under Manly et al. (2002), our sampling efforts represent a special case of a Design I study, where used sites are identified but individual animals are not and where available and used resource units are independently sampled.
We associated resource units in situ with the geomorphic setting, substrate type, vegetation type, and vegetation density that best described the 5 × 5 m area comprising the used or available resource unit (detailed variable definitions in Zeigler et al. 2017Zeigler et al. , 2019a; see also Appendix S1: Table S1).We associated all resource units with the federally recognized recovery unit within which they were located: (1) New England (Maine through Connecticut), (2) New York-New Jersey, and (3) Southern (Delaware through North Carolina; Fig. 1; USFWS 1996).We also characterized resource units ex situ in terms of beach width, elevation, least-cost path distance to moist substrate on low-energy shorelines (henceforth, "distance to MOSH"), and Euclidean distance to the ocean shoreline (referenced to mean high water; henceforth, "distance to ocean") using remotely sensed lidar and orthoimagery (for variable definitions, see Appendix S1: Table S1; Zeigler et al. 2019b).Variables considered in this study were identified as known or suspected drivers of Piping Plover habitat selection based on peer-reviewed literature and expert opinion (e.g., Burger 1991, Patterson et al. 1991, Loegering and Fraser 1995, Jones 1997, Cohen et al. 2008, 2009, Maslo et al. 2011, 2012).
For distance to MOSH, we used orthoimagery to delineate low-energy shorelines along the bayor sound-side of the barrier island, interior water bodies, and ephemeral pools with moist substrates.These environments have previously been described as optimal foraging habitats for this species (Loegering 1992, Cohen and Fraser 2010, Maslo et al. 2012).We did not consider ocean shorelines to be low-energy MOSH in this study, and we may not have included some interior moist substrates if they were not visible at the time orthoimagery was captured.We measured distance to MOSH as a least-cost path distance from a resource unit to the nearest lowenergy shoreline with moist substrate, where barriers to the movement of flightless chicks (e.g., human development, water bodies, dense vegetation) increased the shortest potential chick-walking distance above the direct Euclidean distance.We assumed that a resource unit did not have access to MOSH if movement barriers completely blocked all paths between the resource unit and MOSH.See Appendix S1: Table S1 and Zeigler et al. (2019a, b) for additional information.
We used lidar and orthoimagery captured between 2013 and 2015 (sources listed in Zeigler et al. 2019b) to create raster layers (resolution 5 × 5 m) spanning an entire beach or barrier island for each of the ex situ variables and geomorphic setting.Detailed processing methods and spatial datasets are published separately (Sturdivant et al. 2019a, b, Zeigler et al. 2019b; but see Appendix S1: Table S1).Using ArcGIS Toolbox 10.4.1 (ESRI), we assigned values from the nearest pixel centroid in each raster layer to the resource unit points from the iPlover dataset.The final dataset includes only resource units for which all covariates were known, resulting in 399 resource units from the New England recovery unit (298 used, 101 available), 335 from the New York-New Jersey recovery unit (178 used, 157 available), and 835 from the Southern recovery unit (452 used, 383 available).

Multivariate analysis using Bayesian networks
We used Bayesian networks for multivariate analysis and comparison among recovery units.This methodology allowed us (1) to consider variables that were hierarchical, continuous, and categorical and (2) to incorporate previous knowledge about the relationships among variables.Landscape variables for geomorphic setting, substrate type, vegetation type, vegetation density, beach width, elevation, distance to ocean, and distance to MOSH were represented as nodes on a graph.We considered a single binary output node for the predicted resource unit type (available or used).Variables could be connected to each other or to the output node through edges, which convey correlations or dependencies among variables.Because we had a mix of categorical and continuous variables, we discretized continuous variables into bins for subsequent analyses (Appendix S1: Table S1).
For each recovery unit, we created a database of all used and available resource units and used this dataset to fit a Bayesian network using the R package bnlearn 4.5 (Scutari 2010).We determined the optimal network structure using the score-based hill-climbing learning algorithm and the Bayesian information criterion (BIC) score.This greedy algorithm iteratively adds, removes, and reverses the direction of edges between variables to find the optimal network structure with the lowest BIC score given the data (Scutari 2010).We partially constrained network inference based on previous knowledge of barrierisland dynamics and Piping Plover habitat selection to enhance the structural learning process (Chen and Pollino 2012).We prevented resource unit type (output node) from influencing the environmental variables (input nodes; Appendix S1: Table S2).We forced connections from the substrate type, vegetation density, and distance to MOSH nodes to the resource type node and from the geomorphic setting to the distance to MOSH node because prior studies have found prominent relationships among similar variables (Patterson 1988, MacIvor 1990, Patterson et al. 1991, Jones 1997, Cohen et al. 2008, Maslo et al. 2011).We also forced connections among elevation, distance to ocean, and geomorphic setting because of known inherent correlations among these variables in coastal systems (Young et al. 2011, Halls et al. 2018, Enwright et al. 2019, Lentz et al. 2019).
A variety of network configurations may equally maximize relationships among variables based on the data, so we used a two-step process to determine the final optimal network structure for each recovery unit dataset (Montesinos-Navarro et al. 2018).We first applied nonparametric bootstrap to the data (500 replicates) and estimated the relative frequency that a given edge, or connection, appeared in the graph.We then determined the best-averaged network, assuming a threshold of 0.5 (i.e., a given edge appeared in at least 50% of the bootstrap replicates).The proportion of times that an edge was present in the 500 bootstrapped iterations served as a measure of edge strength in the averaged network, with values approaching 1 associated with high strength and values approaching 0 associated with low strength (Scutari and Denis 2014).
After network structure was determined, we performed two different sets of 10-fold cross-validation: first to validate the accuracy of habitatuse predictions (Marcot 2012) and then to determine network sensitivity to individual nodes.Finally, we trained networks on data from one recovery unit and then tested predictive accuracy for data from another recovery unit to determine the interchangeability of selection patterns between recovery units.We conducted all analyses in bnlearn using the method Bayesian parameter estimation to generate underlying probability distribution tables in the network (Scutari and Denis 2014).

Comparisons among individual variables
To further explore how individual variables may be driving differences among recovery unitspecific Bayesian networks, we assessed habitat selection for individual variables using resource selection ratios (Manly et al. 2002) for categorical variables and Welch's unequal variances t-test (Welch 1947) for continuous variables.For these analyses, individual used and available resource units were combined by recovery unit.The value of a resource selection ratio (W i ) is proportional to the probability of that resource type being utilized, given that the selecting animal has unrestricted access to the entire distribution of available resource types: where O i is the proportion of used resource units in category i, and A i is the proportion of resource in category i in the recovery unit.For geomorphic setting, we were able to reliably quantify the areas of each setting using remotely sensed orthoimagery.Therefore, A i for this variable is considered known and based on area (km 2 ).To provide adequate sample sizes for geomorphic setting categories and to reduce common misclassification errors present in the iPlover dataset (Thieler et al. 2016), we evaluated resource selection ratios for the following combined geomorphic setting categories: (1) beach and backshore, (2) dune and ridge/swale complex, (3) washover, (4) barrier interior, and (5) marsh.For substrate type, vegetation type, and vegetation density, we were unable to reliably differentiate categories and therefore quantify areas at the resolution of available orthoimagery.Therefore, we estimated A i for these variables using a random sample (i.e., the random points in the iPlover dataset).
When considering resource selection ratios, a Wi > 1 indicates selection, Wi ~1 indicates opportunistic use (or use in proportion with availability), and Wi < 1 indicates low use in proportion with availability (Manly et al. 2002).We calculated a standardized index (Β i ) for each variable within each recovery unit where In Β i , ratios are standardized such that they sum to 1 within a recovery unit and give the estimated probability that a randomly selected used resource unit will be in category i if all settings are equally available (Manly et al. 2002).Standardization further facilitates comparison within and among recovery units.Statistical analyses were conducted in the adehabitatHS 0.3.13package (Calenge 2006) in the R statistical environment 3.6.1 (R Development Core Team 2018).
For continuous variables, we report mean AE standard deviation.We compared means for each variable between used and available resource units within a single recovery unit as well as between used resource units among recovery units using Welch's unequal variances t-test.Null hypotheses assumed no difference between used and available resource units or between used resource units in different recovery units.Null hypotheses were rejected when P < 0.05.Statistical analyses were conducted in the R statistical computing environment.

Bayesian networks
In the development of recovery unit-based Bayesian networks, the hill-climbing algorithm detected no additional variables that directly influenced resource unit type for any recovery unit beyond the relationships we forced (Figs.2-4).Therefore, substrate type, vegetation density, and distance to MOSH were the only variables that directly influenced Piping Plover habitat selection throughout the species' U.S. Atlantic coast distribution in this analysis.The remaining variables indirectly influenced selection and, with two exceptions, were also identically linked among recovery units.In the Southern recovery unit, edge strength between beach width and distance to ocean was 0.02, below the threshold for inclusion in the final model, compared to 0.92 and 0.52 in the New England and New York-New Jersey Bayesian networks, respectively (Figs. 2-3).The algorithm also detected a relationship between substrate type and vegetation density in the Southern recovery unit but not in the other two recovery units (Fig. 4).
The resulting averaged Bayesian networks had relatively low error rates ranging from 0.05 (Southern recovery unit) to 0.18 (New York-New Jersey; Table 1).The largest net change in error rate for the New England network occurred with the removal of the vegetation density node, which reduced accuracy by 4% (Table 1).The removal of the distance to MOSH node slightly increased accuracy (Table 1).For the New York-New Jersey Bayesian network, the largest losses in accuracy occurred with the removal of nodes for substrate type, vegetation density, and distance to MOSH (Table 1).The removal of the beach width node slightly increased accuracy in this recovery unit (Table 1).Finally, the removal of the node for substrate type resulted in the largest loss in accuracy for the Bayesian network trained on resource units from the Southern recovery unit (Table 1), and accuracy did not improve with the removal of any one node.
We also observed that the accuracy of Bayesian networks trained on resource units from one recovery unit and tested with resource units from another resulted in substantially higher error rates in a manner that suggests a latitudinal v www.esajournals.orggradient (Table 2).The highest error rates occurred when data from New England were used to test networks trained on New York-New Jersey and Southern data or when data from the Southern recovery unit were used to test networks trained on New England or New York-New Jersey data.However, the error rate when New York-New Jersey data were used to test networks trained on either of the other two networks was not substantially different from that of the New York-New Jersey network (Table 2).This suggests that patterns observed in either the Southern or New England recovery units were not shared by other recovery units; however, some patterns observed in New York-New Jersey were observed in both New England and the Southern recovery units.

Univariate analysis
We observed consistent use of resources throughout the study area for several variables (center portion of overlapping ovals in Fig. 5).Nesting Piping Plovers showed strong selection for resource units with substrates that are a mix of sand and shell, gravel, or cobble (henceforth, "mixed substrates") but low selection of predominantly sandy substrates compared to availability (Table 3).When vegetation type but not density Fig. 2. Averaged Bayesian network for the New England recovery unit, with structure derived using a hillclimbing algorithm and data on Piping Plover habitat use in that recovery unit.Landscape variables are represented as nodes (boxes) with probability distributions (shown here, distributions associated with used resource units only).Correlative relationships among variables, or edges, are represented as black or gray arrows.Black arrows signify edges that were determined by the learning algorithm directly, and the numeric value associated with these edges shows the proportion of the 500 bootstrapped replications in which that connection was identified.For the final averaged networks, we only retained edges that were present in at least half of 500 bootstrapped replications of the network structure learning process.Gray arrows are edges that were forced into the network learning process based on our previous knowledge of the species' habitat selection patterns.Networks were built for display purposes here in the program Netica 6.05.
was considered, birds either selected (New England, New York-New Jersey) or opportunistically used (Southern) resource units with no vegetation over units with herbaceous or shrub vegetation (Table 3).In addition, nesting Piping Plovers exhibited strong selection for resource units with a sparse (<20%) cover of vegetation (Table 3).Birds never used resource units in marsh or containing mud/peat, water, forest, or human development.We also observed that nesting Piping Plovers used resource units at significantly higher elevations and at shorter distances to the ocean shoreline compared to available resource units in all recovery units (P < 0.05; Table 4; Appendix S1: Fig. S1).Beach widths for used and available resource units did not significantly differ in New York-New Jersey and New England.However, used resource units were found on significantly wider beaches in the Southern recovery unit compared to available resource units (P = 0.002).Among recovery units, the only significant variable difference for used resource units was between New England and New York-New Jersey, where resource units Fig. 3. Averaged Bayesian network for the New York-New Jersey recovery unit, with structure derived using a hill-climbing algorithm and data on Piping Plover habitat use in that recovery unit.Landscape variables are represented as nodes (boxes) with probability distributions (shown here, distributions associated with used resource units only).Correlative relationships among variables, or edges, are represented as black or gray arrows.Black arrows signify edges that were determined by the learning algorithm directly, and the numeric value associated with these edges shows the proportion of the 500 bootstrapped replications in which that connection was identified.For the final averaged networks, we only retained edges that were present in at least half of 500 bootstrapped replications of the network structure learning process.Gray arrows are edges that were forced into the network learning process based on our previous knowledge of the species' habitat selection patterns.Networks were built for display purposes here in the program Netica 6.05.
We noted variability in the characteristics of used resource units for all remaining variables (non-overlapping regions in Fig. 5).First, proportions of used resource units in beach/backshore areas, dune complexes, washover, and barrier interior geomorphic settings varied by recovery unit (Table 3).In New England, although breeding Piping Plovers selected both beach/backshore areas and washovers, they were two times more likely to select beach/backshore according to standardized resource selection ratios (Bi; Table 3).In New York-New Jersey, although birds also selected beach/backshore areas, they were two to three times more likely to select washover or dune complexes over beach/backshore when scaled by availability in this recovery unit (Table 3).In the Southern recovery unit, birds selected resource units in the beach/backshore and in washovers, with breeding pairs four times more likely to nest in washover than any other geomorphic setting (Table 3).
We also observed subtle variability in vegetation density in used resource units.Although Correlative relationships among variables, or edges, are represented as black or gray arrows.Black arrows signify edges that were determined by the learning algorithm directly, and the numeric value associated with these edges shows the proportion of the 500 bootstrapped replications in which that connection was identified.For the final averaged networks, we only retained edges that were present in at least half of 500 bootstrapped replications of the network structure learning process, and a link from beach width to another variable did not meet that threshold in this network.Gray arrows are edges that were forced into the network learning process based on our previous knowledge of the species' habitat selection patterns.Networks were built for display purposes here in the program Netica 6.05.birds consistently selected for sparse vegetation in all recovery units (Table 3), breeding pairs in New England exhibited a higher tolerance for resource units with moderate (20-90% cover) vegetation.One percent of resource units even contained dense vegetation (Table 3).Selection for sparse vegetation over no vegetation or moderate vegetation became stronger with decreasing latitude, and no resource units contained dense vegetation in the New York-New Jersey or Southern recovery units (Table 3).
We observed significantly greater elevations associated with used resource units compared to available resource units among all recovery units (P < 0.05).Furthermore, we observed significant differences in the mean elevation of used resource units between all pairs of recovery units (P < 0.05), with elevation decreasing with latitude (Table 4; Appendix S1: Fig. S1).Used resource units were also significantly closer to the ocean shoreline in New England compared to New York-New Jersey (P < 0.001) and the Southern recovery units (P < 0.001).We did not observe significant differences in distance to ocean shoreline for used resource units in New York-New Jersey and the Southern recovery units (P = 0.264; Table 4).
In a comparison of used versus available resource units within a single recovery unit, used resource units were significantly farther from MOSH in New England (P = 0.01) but significantly closer to MOSH in the Southern recovery unit (P < 0.001) compared to available resource units.There was not a significant difference in distance to MOSH between used and available resource units in New York-New Jersey (P = 0.94; Table 4).In a comparison of used resource units among recovery units, resource units with access to MOSH in New England were significantly farther from MOSH than those in New York-New Jersey (P < 0.001) and Southern recovery units (P < 0.001; Table 4; Fig. 6).We did not observe a significant difference in distance to MOSH for used resource units in the New York-New Jersey and Southern recovery units (P = 0.05; Table 4; Fig. 6).Although large proportions of used resource units were located within 1 km of MOSH in all three recovery units (Fig. 6), more used resource units were found at distances >1 km to MOSH in New England (59% of used resource units) compared to the New York-New Jersey (11%) and Southern (7%) recovery units.We also noted that a higher percentage of used resource units had no clear path to MOSH in New England (11%) compared to the New York-New Jersey (8%) and Southern (1%) recovery units (Table 4).Notes: We iteratively removed variables and recalculated the network's predictive error rate through 10-fold cross-validation to determine the network's sensitivity to a given variable †.Dashed lines indicate scenarios for which an error rate was not calculated ‡.
‡ Averaged Bayesian network structures included beach width in the NE and NY-NJ recovery units, and we tested how the removal of this variable affected error rates.However, this variable did not meet sufficient thresholds to be included in the final averaged network for the South recovery unit.Therefore, we tested the change in error rate with the inclusion of this variable in this recovery unit's network.

Range-wide habitat selection patterns
In all three recovery units, Piping Plovers selected for mixed substrates and sparse vegetation, avoided dense vegetation, and used areas that were significantly closer to the ocean and at higher elevations than available resource units.Furthermore, vegetation density and substrate type were consistently among the top-ranked variables in Bayesian network sensitivity analyses for all three recovery units, suggesting that habitat selection was largely driven by these characteristics.These universal habitat selection patterns throughout the study area are also consistent with more localized studies.On Long Island, New York, Cohen et al. (2008) found that all 1-m 2 plots containing Piping Plover nests had <47% vegetative cover (mean 7.5 AE 1.7%), but models indicated that nest plots were more likely to be vegetated than random plots.Wilcox (1959) also reported that Piping Plovers abandoned his Long Island site as vegetation encroached.Twom 2 plots with Piping Plover nests in New Jersey had a mean vegetative cover of 12.9 AE 1.4% in dune complexes and 2 AE 6% in all other geomorphic settings in New Jersey (Maslo et al. 2011).Farther south, mean vegetative cover within 1 m of nests in Maryland and Virginia ranged from 9.0 to 16.1% (Patterson 1988), while Loegering (1992) observed that 81% of nests occurred in unvegetated areas in Maryland.The slightly higher tolerance for areas with moderate vegetation that we observed in New England may also be reflected by MacIvor (1990), who found that 0.5-m 2 plots had a mean vegetative cover of 15.2 AE 23.6% in Massachusetts.In addition, the particularly high importance of coarse-grained or shelly substrates for Piping Plovers is also evident in localized studies in Massachusetts  ------------------------Notes: Ratios were calculated as the geomorphic setting, substrate type, vegetation type, and vegetation density of resource units proportionately used by nesting Piping Plovers (Prop Used) given the availability of those resources (Prop Avail) †.Habitat selection ratios in bold indicate selection (or use above what is expected given availability) of a particular resource type.Because no nests were found in dense vegetation in New York-New Jersey or Southern recovery units, habitat selection ratios were not calculated for this variable in these recovery units (dashed lines).
† W i > 1 indicates selection (in bold), W i = 1 indicates opportunistic use (or use in proportion with availability), and W i < 1 indicates use in low proportion given availability.Used resource units were never located in marsh, mud/peat, water, forest, or development, and these variables were not included in analysis.Availability for geomorphic settings was known, as interpreted through remotely sensed spatial datasets.Availability for the remaining variables was estimated using a random sample of available resource units.(MacIvor 1990, Jones 1997), New Jersey (Burger 1987, Maslo et al. 2011, Grant et al. 2019), and Maryland (Patterson 1988).In contrast with our results, a study on Long Island (Cohen et al. 2008) suggested that predominantly sandy substrates were preferred.However, the authors found model support for selection of coarsegrained sands, which were less available than at other sites in the literature, perhaps due to recent renourishment in the study site with predominantly fine-grained sediment.These consistently selected attributes throughout the range demonstrate a fundamental component of the species' ecological niche, which optimizes fitness under the constraints of interand intraspecific competition (Morris 2003).For birds nesting in coastal environments, predation and nest inundation are among the primary causes of nest failure (Patterson et al. 1991, USFWS 1996, Cohen et al. 2009, Maslo et al. 2016), and habitat selection could reflect mechanisms for minimizing these threats (Espie et al. 1996, Lauro and Tanacredi 2002, Greenberg et al. 2006, Doherty and Heath 2011).Sandy substrates, particularly those mixed with shell or rock fragments, provide camouflage for mottled Piping Plover eggs and allow adults and chicks to "hide in plain sight" (Fraser and Catlin 2019).Preference for nest sites with little to no vegetation also affords Piping Plovers visibility to detect approaching predators and deploy common nest defense behaviors-including false brooding, running, feigning injury, or staying motionless-to maximize benefits of their cryptic coloration (Cairns 1982), while sparse vegetation may reduce detection of nests by avian predators (Maslo et al. 2016).Some studies have suggested that these characteristics increase nest success for plovers and other ground-nesting species (Prindiville Gaines and Ryan 1988, Colwell et al. 2011, Troscianko et al. 2016); however, others have found contrary relationships (Patterson et al. 1991, Darrah et al. 2018).

Latitudinal gradients and differences across recovery units
In addition to universal drivers of habitat selection noted above, we also found evidence for latitudinal gradients in other attributes.First, although network structure was largely the same, networks trained on one dataset had a reduced capacity to predict habitat selection in another recovery unit.Analyses of individual variables indicate that differences in selection patterns are likely driven by latitudinal gradients in selection of geomorphic setting and distance to MOSH, which appear to be related.The typical Atlantic coast barrier-island cross section starts with beach/backshore areas closest to the ocean shoreline, followed by dune complexes, upland areas with thickets and forests, and, finally, bay-or sound-side beaches or marsh flats along the back-barrier (Godfrey 1976, Hayes 1979, Davis Jr 1994, Neuendorf et al. 2011).Washovers, which occur frequently but can either be intermittently spaced or widespread on low-elevation barrier islands, are typically formed by storm-induced high water conditions and waves overtopping dunes and, in extreme cases, marshes (Morton and Sallenger Jr 2003).These features tend to be flat, low-elevation fans of sandy sediment that form deeper into the barrier island's interior, sometimes stretching the entirety of the barrier-island cross section (Hayes 1979, Morton and Sallenger Jr 2003, Neuendorf et al. 2011).Piping Plovers in New England most strongly selected for areas closest to the oceanbeach/backshore areas-while those in the Southern recovery unit predominantly selected washovers located deeper into the barrier island.Selection patterns in New York-New Jersey were intermediate between those of the New England and Southern recovery units, with selection occurring in relatively equal frequency in beach/ backshore, dune complexes, and washover.This pattern is also evident in the probability distributions associated with geomorphic setting nodes in the Bayesian networks.
The tendency to nest deeper into the island's interior along the gradient from north to south echoes the shorter average distances to MOSH for used resource units and lower percentages of nests without access to MOSH as one moves south.Previous studies have indicated that Piping Plovers spent more time foraging in bay-or sound-side intertidal zones, ephemeral pools, and ponds-features that we categorize as MOSH and that tend to occur on the back-barrier or barrier interior-than in other areas (e.g., ocean intertidal, dunes; Fraser et al. 2005, Cohen and Fraser 2010, Maslo et al. 2012).Washovers tend to provide chick access unimpeded by v www.esajournals.orgdense vegetation to the ephemeral pools and bay-side shorelines that constitute MOSH.
This increasing reliance on washover nesting habitats and access to MOSH as one moves south is reflected in literature from more localized studies.At the Cape Cod National Seashore in Massachusetts, nests occurred on beaches with access to bay-side feeding habitat in significantly greater proportion to availability, but 63% of nests occurred on beaches without access to bayside foraging habitats (Jones 1997).On Long Island, New York, Piping Plovers nested and reared broods on all 1-km beach segments with ephemeral pools or bay tidal flats but nested in less than half of segments without those foraging habitats (Elias et al. 2000).Cohen et al. (2009) documented rapid Piping Plover colonization of a Long Island site after a series of storms created new nesting habitat adjacent to bay-side intertidal flats, while others (Maslo et al. 2019, Zeigler et al. 2019a) observed rapid colonization in washovers by nesting piping plovers throughout the Mid-Atlantic region following Hurricane Sandy in 2012.In New Jersey, Piping Plover nests were split evenly between sites with and without nonocean foraging habitats, but birds initiated nests in washovers if those features were present (Maslo et al. 2011).Nesting habitat selection studies related to MOSH are limited in the south; however, chick fledging success was substantially higher (69%) for Piping Plover broods that had access to bay flats or tidal pools compared to chicks that only foraged on the ocean beach (19%) in Maryland and Virginia (Patterson et al. 1991).Likewise, Piping Plover chicks reared on the bay beach and island interior in Maryland had significantly higher daily survival rates and foraging rates than chicks reared on the ocean beach (Loegering and Fraser 1995).This contrasts with a similar Massachusetts study, which found that fledging success did not differ significantly with the availability of bay access (Jones 1997).Broods with access to mudflats surrounding an interior pond in Rhode Island had higher survivorship and fledging success compared to those foraging on ocean intertidal, but higher human disturbance along the ocean shoreline was noted as a confounding variable (Goldin and Regosin 1998).
Although our study design does not support conclusions as to why access to MOSH and use of washovers becomes increasingly important with decreasing latitude, we hypothesize that it is related to changes in ocean wrack biomass, species composition, species richness/diversity, or community composition as one moves south.High energy beach shorelines tend to have limited productivity (Liebowitz et al. 2016), and back-barrier shorelines are typically cited as ideal foraging habitats for Piping Plovers (Loegering andFraser 1995, USFWS 1996).However, sea wrack deposits-defined as dead, shore-cast seaweeds and grasses-can provide a rich nutritional supply on ocean-facing beaches for insects and other invertebrates, which support higher trophic levels (e.g., crabs, birds; Liebowitz et al. 2016).Although Liebowitz et al. (2016) did not find systematic latitudinal patterns on the North American Pacific coast, the authors did find strong regional signals in wrack composition and abundance that were driven by local factors, particularly physical shoreline characteristics (e.g., slope, substrate, wave exposure, etc.) and donor marine ecosystem habitat.Beaches with coarse substrates (e.g., cobble) and that were sheltered from wave energy-morphological conditions more common on New England pocket beaches compared to the New York-New Jersey and Southern barrier islands (Hapke et al. 2010)-retain more wrack (Orr et al. 2005, Barreiro et al. 2011, Wickham et al. 2020).We hypothesize that New England may support higher ocean wrack availability with richer invertebrate prey base compared to more southern latitudes.A foraging study in Cape Cod, Massachusetts, indicated high arthropod abundances on fresh wrack and suggested Piping Plover chicks preferred this substrate (Hoopes 1993).However, similar studies in New York found significantly lower arthropod abundances in ocean intertidal zones and fresh wrack compared to those on the bay shoreline (Cohen et al. 2009) or interior ephemeral pools (Elias et al. 2000).In Maryland, terrestrial arthropod abundance was significantly higher and chicks had higher foraging rates on bay beaches and on the island interior compared to ocean beaches (Loegering and Fraser 1995).
Our results and the broader literature preliminarily suggest that Piping Plovers may become more reliant on or limited to foraging habitats on the back-barrier as one moves from north to south.This would also explain the higher selective pressure for more interior habitats, like washovers that are in closer proximity to backbarrier foraging habitats, in the south.Further study is necessary to support or refute this hypothesis.We caution, however, that the value of low-energy shorelines with moist substrates in New England should not be dismissed.Our study found a large proportion of nests within 1 km of MOSH, and Fraser et al. (2005) observed a large proportion (87%) of pre-nesting Piping Plovers foraging along sound and tidal pond habitats on South Monomoy Island, Massachusetts.
Latitudinal differences in the distribution of geomorphic settings and reliance on low-energy shorelines with moist substrates may help explain greater overall growth in abundance of breeding Piping Plover pairs in New England (445%) than in the New York-New Jersey and the Southern recovery units (122% and 148%, respectively) between 1989 and 2018 (USFWS 2019).We argue that birds in New England may be less habitat-limited given that they are much less restricted to washovers and areas with easy access to MOSH than in the Southern recovery unit.While beach/backshore is a ubiquitously available geomorphic setting and comprises 13-14% of available area in all 3 recovery units compared to 4-10% available area of washovers, beach/backshore is much more heavily used in New England.This could be a major factor explaining the 1.40 breeding pairs of piping plovers per mile of sandy ocean beach in New England in 2015, compared with 0.81 and 0.75 breeding pairs per mile, respectively, in New York-New Jersey and the Southern recovery units (Appendix S1: Table S3).Furthermore, variable trends in the number of breeding pairs in these two recovery units (USFWS 2019) appear to respond strongly to availability of washover and accessibility of nest sites to MOSH (Schupp et al. 2013, Zeigler et al. 2019a).Cohen et al. (2009) found that peak density of Piping Plovers in the portion of their Long Island, New York study area where chicks had access to both ocean-and bay-side foraging habitat was more than double the density of Piping Plovers in adjacent habitat without accessible MOSH.All six Atlantic Coast examples of Piping Plover population irruptions following storm-induced habitat creation that are presented by Robinson et al. (2019) occurred in New York-New Jersey and the Southern recovery units.Hence, carrying capacity and resulting population trends in more southerly latitudes-where plovers use more consistently available beach/backshore habitats less often-may be especially sensitive to natural processes that create washover habitat and to anthropogenic activities that accelerate its loss.
Finally, we observed strong selection for dunes in New York-New Jersey, opportunistic use in New England, and low use relative to availability in the Southern recovery unit.We found that all used resource units in dune complexes were located in bare or sparsely vegetated areas with mean elevations of 2.3 AE 1.1 m (New England), 1.5 AE 0.8 m (New York-New Jersey), and 1.5 AE 0.6 m (Southern).Similar use of low, bare dunes was observed in other studies.For example, Maslo et al. (2011) found that dunes containing nests in New Jersey averaged 1.1 AE 0.1 m in height from the apex to the seaward toe, were significantly lower than dunes within the surrounding landscape, and were never higher than 3.1 m.In Maryland, Loegering (1992) found that 87% of nests were located on flat areas or on shallow dunes <0.75 m high.We hypothesize that strong selection for dunes in New York-New Jersey observed in our study is due in part to use of small remnant dunes close to MOSH on beaches created during widespread restoration of engineered dunes in New York and New Jersey beaches following Hurricane Sandy in 2013 and 2014.These constructed dunes may have provided sparsely vegetated habitat and access to MOSH at the time of our study, but further analysis is needed to determine if these dunes fostered rapid vegetation regrowth or impeded overwash that would accelerate future habitat loss.

Implications for management
Our study shows evidence of both range-wide consistency in the importance of some habitat characteristics as well as latitudinal patterns in others for nesting Piping Plovers along the U.S. Atlantic coast.Resilient Piping Plover populations distributed across their U.S. Atlantic coast breeding range require sufficient nesting habitats that provide all of the characteristics required for successful breeding-particularly mixed substrates and sparse vegetation.Caution should, therefore, be exercised when implementing projects (e.g., building artificial dunes, vegetation planting, constructing sea walls) that accelerate vegetation succession or impede natural coastal processes that would otherwise maintain habitat.In addition, beach renourishment projects that utilize borrow materials of similar color and texture to native substrates and contain coarser grains, including shell fragments, are more likely to create or maintain nesting habitat for Piping Plovers.
Our results also highlight particular areas of consideration for each recovery unit when designing and implementing beach management programs.Coastal management practices that allow washover processes and provide connectivity between ocean-facing beaches and backbarrier foraging areas are particularly critical to habitat suitability for this species in the New York-New Jersey and Southern recovery units.Built features (e.g., buildings, sand fencing, sea walls, constructed dunes) designed to stabilize coastal areas, prevent overwash, or otherwise obstruct the cross-island movement of sediments (Kratzmann and Hapke 2012, Lentz et al. 2013, Schupp et al. 2013, Rogers et al. 2015) could have a particularly detrimental impact on habitat availability and ultimately population dynamics in this part of the range.Piping Plover chicks are vulnerable to crushing by off-road vehicles range-wide, but activities such as beach renourishment and sand raking that degrade backshore habitats and invertebrate populations may be especially detrimental to nesting Piping Plovers in New England and at those New York-New Jersey beaches that lack chick access to MOSH (Elias et al. 2000, Peterson et al. 2000, Peterson et al. 2006, Kluft 2009).Our results suggest that considering regional habitat preferences in the design and evaluation of proposed habitat modifications, including restoration projects, could improve the likelihood of project success.Furthermore, caution should be taken when extrapolating from studies conducted at sites in distant parts of the range, as Piping Plover selection patterns differed by regional recovery unit.Results and conclusions presented here rely on nest presence as a proxy for habitat availability instead of, for example, egg fate or fledging rate.Therefore, we cannot make conclusions regarding whether selected habitat was optimal for nesting or acted as an ecological trap.Continued investigation into which habitats are not only selected but also provide superior fledgling rates is an important next step for managing this and other federally listed species.

Fig. 4 .
Fig. 4. Averaged Bayesian network for the Southern recovery unit, with structure derived using a hill-climbing algorithm and data on Piping Plover habitat use in that recovery unit.Landscape variables are represented as nodes (boxes) with probability distributions (shown here, distributions associated with used resource units only).Correlative relationships among variables, or edges, are represented as black or gray arrows.Black arrows signify edges that were determined by the learning algorithm directly, and the numeric value associated with these edges shows the proportion of the 500 bootstrapped replications in which that connection was identified.For the final averaged networks, we only retained edges that were present in at least half of 500 bootstrapped replications of the network structure learning process, and a link from beach width to another variable did not meet that threshold in this network.Gray arrows are edges that were forced into the network learning process based on our previous knowledge of the species' habitat selection patterns.Networks were built for display purposes here in the program Netica 6.05.

Fig. 5 .
Fig. 5. Summary of habitat selection patterns across the Piping Plover's U.S. Atlantic coast breeding range.Each oval represents selection patterns observed in a given recovery unit (RU).Variables in areas of overlap indicate consistent patterns observed in two or more RUs, while variables in one oval indicate selection patterns observed in only that RU.Categorical variables were compared using resource selection ratios and indicate whether a variable was selected for (+), used in proportion with availability (0), or used in low proportion given availability (−).Continuous variables were compared using Welch's unequal variances t-tests (P < 0.05).Abbreviations: vegetation (veg); significant (sig.);low-energy shorelines with moist substrate (MOSH).

Table 1 .
Error rates for averaged Bayesian networks derived from data on Piping Plover habitat use in the species' New England (NE), New York-New Jersey (NY-NJ), and Southern (South) recovery units.

Table 2 .
Bayesian network error rate for networks trained on one recovery unit-based dataset and tested with data from another.Where training and testing data came from the same recovery unit, we used 10-fold cross-validation, and error rate represents the mean ensemble error rate.

Table 3 .
Raw (W i ) and standardized (B i ) habitat selection ratios.

Table 4 .
Mean (AE standard deviation) values for landscape characteristics in resource units used by nesting piping plovers compared to available/unused resource units, evaluated by federally designated recovery unit †.