Integrating dynamic processes into waterfowl conservation prioritization tools

Traditional approaches for including species' distributions in conservation planning have presented them as long‐term averages of variation. Like these approaches, the main waterfowl conservation targeting tool in the United States Prairie Pothole Region (US PPR) is based primarily on long‐term averaged distributions of breeding pairs. While this tool has supported valuable conservation, it does not explicitly consider spatiotemporal changes in spring wetland availability and does not assess wetland availability during the brood rearing period. We sought to develop a modelling approach and targeting tool that incorporated these types of dynamics for breeding waterfowl pairs and broods. This goal also presented an opportunity for us to compare predictions from a traditional targeting tool based on long‐term averages to predictions from spatiotemporal models. Such a comparison facilitated tests of the underlying assumption that this traditional targeting tool could provide an effective surrogate measure for conservation objectives such as brood abundance and climate refugia.


| INTRODUC TI ON
The traditional approach to including species' distributions in conservation planning has been to pool spatiotemporal variation and create a static snapshot of conditions (Pressey et al., 2007). However, species' distributions and the processes on which they depend are not static, and conservation plans require consideration of the dynamic and highly complex ecological processes that change and maintain the biodiversity within an ecosystem (Pressey et al., 2003(Pressey et al., , 2007Soule et al., 2004;Van Teeffelen et al., 2012;Wilson et al., 2009). A highly variable climate, for example, might cause changes in species' habitat use (Groves et al., 2012). Alternatively, natural disturbances can increase the overall habitat needed to support viable populations (Allison et al., 2003). Highly mobile species pose additional challenges for conservation planners, because their natural intraand inter-annual movements also require consideration (Gilmore et al., 2007;Johnston et al., 2020;Runge et al., 2014;Schuster et al., 2019).
North American waterfowl are perhaps one of the best studied highly mobile groups in the literature with a long history of management and conservation planning (North American Waterfowl Management Plan Committee, 2012, 2018. Species distribution models for waterfowl have helped to support this history of conservation, particularly in the Prairie Pothole Region (PPR: Figure 1), where a disproportionately large number of North American waterfowl breed each year. Most waterfowl modelling efforts have focused on describing patterns of breeding pair abundance and distribution (Barker et al., 2014;Doherty et al., 2015;Feldman et al., 2016;Janke et al., 2017). More recently, there have been efforts to model waterfowl brood abundance and distribution in the PPR as well (Carrlson et al., 2018;Kemink et al., 2019;Walker, Rotella, Schmidt, et al., 2013). Both avenues of investigation have highlighted spatial and temporal trends in both pair and brood distributions (Doherty et al., 2015;Janke et al., 2017;Kemink et al., 2019). However, we know of no studies that have contrasted pair and brood distributions during these different stages of reproduction. Further, the prevailing trend for conservation planning in the PPR still focuses on pooling variation to create a static distribution for targeting purposes (Prairie Habitat Joint Venture, 2014; Barker et al., 2014;Prairie Pothole Joint Venture, 2017, but see Humphreys et al., 2019;Adde et al., 2020).
In the US PPR, the leading tool for supporting decisions about breeding waterfowl conservation is developed through methods that parallel the traditional use of static distributions. The Waterfowl Breeding Pair Accessibility Map, colloquially known as the thunderstorm map ( Figure 2; Reynolds et al., 2006Reynolds et al., , 2007, is used to display categorical ranges of duck pair numbers (mallard  (Niemuth et al., 2010). These pair abundance values are scaled to a 0.152 km 2 resolution grid and were collected through an annual regional survey known as the "Four Square Mile Survey" . To produce the map of "accessibility," they are adjusted by species-specific constant values of waterfowl hen travel distances from core breeding wetlands to upland nest sites during the breeding season (Reynolds et al., 2006 [ personal communication, Chuck Loesch, USFWS).
While these pair abundance values and their derivatives have provided support for decades of valuable conservation work, they preclude the explicit consideration of wetlands' inter-annual wet-dry cycles and ignore any intra-annual changes in wetland ponding across the region. Historically, the US Fish and Wildlife Service (FWS) conducted brood count surveys in the late summer to complement the May breeding population and habitat surveys.
However, due to funding cuts and concern about methodology, this data collection was curtailed in the early 2000s. Conservation planners in the PPR might consequently be overlooking areas that have conservation value to waterfowl during periods of extreme weather variation (e.g. drought or deluge: Doherty et al., 2015;Wilson et al., 2009) or during the brood rearing period (Carrlson et al., 2018).
Periods of drought and deluge are a well-known characteristic of the PPR (Johnson et al., 2010;Karl & Riebsame, 1984;Larson, 1995;Niemuth et al., 2010;Woodhouse & Overpeck, 1998). These weather patterns are the primary drivers of the region's wetland hydrology F I G U R E 1 Prairie Pothole Region and the major North American Level III Ecoregions that it encompasses and thus of aquatic invertebrate abundance and diversity (Euliss & Mushet, 2004;Euliss et al., 1999), which fulfil dietary requirements for breeding duck pairs, nesting hens and growing waterfowl recruits (Cox et al., 1998;Stafford et al., 2016). While both the adults and broods of wetland obligate birds often depend on resources provided by wetlands for survival and growth during the breeding season, the amount and type of habitat available to and used by each group can be quite different (Carrlson et al., 2018;Johnson et al., 2010).
Breeding dabbling duck pairs arrive in the early spring (April-May) to establish territory in the PPR prior to nesting. It is widely accepted that densely ponded areas attract the highest number of pairs. At more local extents, small, seasonal (sensu Stewart & Kantrud, 1971) wetlands tend to provide the best habitat for breeding dabblers (Bartzen et al., 2017;Fields, 2011;Reynolds et al., 2006). These ponds receive most of their water as spring snowmelt running over frozen ground (Hayashi et al., 2016) and thus are available earlier in the spring than their deeper semipermanent counterparts. Dabbling duck pairs feed along the edges of these ponds, concealing themselves from predators and conspecifics (Bartzen et al., 2017;Kantrud & Stewart, 1977;Reynolds et al., 2006). Many of the temporary ponds used by dabbling duck pairs settling in the PPR are dry in the late summer (July-August) by the time waterfowl hens are raising broods (Johnson et al., 2010).
Greater numbers of broods are often found on the deeper seasonal or semipermanent ponds (Kemink et al., 2019;Talent et al., 1982). As a result, conservation targeting for successful reproduction requires a diverse mix of wetland types, or hydrologic regimes, ranging from temporary, shallow ponds able to thaw early in the year, to deeper semipermanent wetlands that will remain inundated through hot, dry summers.
In this paper, we develop spatiotemporal models of waterfowl pair and brood abundance that incorporate layers describing water and land use changes on the landscape. Specifically, we seek to use these models to evaluate: (a) whether the pair abundance values scaled to a 0.152 km 2 resolution grid (hereafter averaged pair abundance) that are used to develop the thunderstorm map are a good surrogate measure for other conservation objectives including brood abundance and climate refugia and (b) whether spatiotemporal predictions of pair abundance provide a surrogate measure for brood abundance.

| Study area
The PPR is a 700,000 km 2 landscape dominated by small, shallow wetlands and historically covered in perennial grasslands (van der Valk, 1989). The region's major land uses, physiography, geography and climate have been described in detail elsewhere (Cowardin & Golet, 1995;Johnson et al., 1994;Reynolds et al., 2006). The PPR covers five states and three Canadian provinces. However, independently collected brood data and the averaged FWS pair abundance data were available only for the PPR in North Dakota, South Dakota and part of the Montana PPR. Similarly, the annual pair data we used for this analysis were not available for the Iowa and Minnesota portions of the PPR. Consequently, any spatial comparisons made between distributions were limited to the PPR of North Dakota, South Dakota and eastern Montana. The time period for which we modelled pair and brood abundance (2008-2017) is described as one of the wetter periods of the PPR's climatic history since the mid-1900s.
However, as was typical for the region, precipitation and temperature varied spatially within and between years (NOAA, 2020).

| Spatiotemporal breeding pair data
We used data from the publicly accessible Waterfowl Breeding  (Smith, 1995: Figure 3a). During the annual survey, the transects are flown by a fixed-wing aircraft 30-45 m above the ground. An observer and the pilot count ducks and ponds 200 m on both sides of the segments (Smith, 1995). Ground counts are also completed simultaneously to allow estimation of detection rates (see Smith, 1995).
The dependent variable in our duck pair analysis was the total number of dabbling duck pairs counted within a segment. We included the dabbling duck species considered in the averaged pair abundance data, which are the five most common dabbling duck species in the PPR: mallard, gadwall, Northern pintail, Northern shoveler and blue-winged teal. These species are the most targeted in wetland and waterfowl management plans in the region (Prairie Pothole Joint Venture, 2017). We calculated the total number of pairs per segment from raw counts such that: where P represents a duck pair (male and female), LM (isolated lone drake) represents an indicated pair, and VIF represents the detection adjustment factor specific to the strata relevant to that segment, year and species (Smith, 1995). We included only counts for segments that were completely within the US or Canadian PPR ( Figure 3b).

| Spatiotemporal pair models: predictor variables
The predictor variables we tested were supported by previous studies and tied ecological and anthropogenic processes together. They included two variables describing wetlands and moisture and variables describing our hypotheses about human-driven processes ( Table 1).
The variables describing wetlands included the number of wet wetlands counted per segment in the survey (pond count) and climate moisture index, which is the difference between annual precipitation and potential evapotranspiration on a vegetated landscape. Landscapes with more wet area and higher wetland densities overall generally provide more habitat for breeding duck pairs (Johnson & Grier, 1988). As most wetlands used by pairs in the spring are filled through rainfall and snowmelt, we expected areas with more ponded wetland counts and higher climate moisture indices to coincide with higher pair counts each year (Doherty et al., 2015;Johnson et al., 2010;Zimpfer et al., 2009). Total = (P + LM) × VIF TA B L E 1 Description of fixed effects incorporated in pair and brood abundance models with brief justifications for their inclusion as well as the sources of raw data Areas with more growing degree days are more conducive to cropping and will be less likely to have large expanses of perennial cover available for nesting ducks Human-driven processes like agriculture that alter the landscape might also impact pair abundance. Perennial cover surrounding wetlands has been shown to increase nest success and productivity and thus is believed to be the preferred habitat of pairs (Greenwood et al., 1995;Reynolds et al., 2001;Stephens et al., 2005, but see . We included a variable to represent the amount of perennial cover surrounding a survey segment as well as the amount of growing degree days (degree days >5°C; Doherty et al., 2015). We expected that perennial cover would demonstrate a positive relationship with breeding pair abundance, while areas with higher growing degree days would be more conducive to cropping and thus have less habitat suitable for breeding duck pairs. Like Doherty et al. (2015), we summarized the climate moisture index, perennial cover and degree day variables using a moving window analysis in ArcMap 10.6 with an area equivalent to the average area of a survey segment (11.52 km 2 ). We extracted the value of the resulting layers to the centroid of each survey segment within the PPR.

| Spatiotemporal pair models: analysis
Preliminary analyses indicated that the Poisson distribution provided the best fit for dabbling duck pair abundance between 2008 and 2017 and that residuals contained spatial and temporal correlation (Zuur et al., 2007). We used Bayesian hierarchical models to examine the data. The hierarchical approach allowed us to test several hypotheses about the structure of spatial and temporal correlation.
We binned the data by year and randomly selected 80% of the data for the analysis and withheld 20% of the dataset to test model fit.
The remaining analysis contained two stages. We first compared support for different global model structures with regards to the presence or absence of spatial and/or temporal correlation. Global models contained all four fixed effects: pond count, climate moisture index, perennial cover and growing degree days. We assessed support for the fixed effects within the most supported model structure in the second stage of analysis.
In the first stage of our analysis, we considered six model structures to test different hypotheses about how the spatial random field changed over time. The first model contained no spatial or temporal correlation and was an ordinary Poisson model (M1). We modelled spatial correlation in M2-M6 using the stochastic partial differential equation (SPDE: Lindgren et al., 2011).
The SPDE approach models spatial autocorrelation across a triangular mesh rather than a grid or polygons and has been used to model spatial autocorrelation in a similar manner on waterfowl data from eBird (Humphreys et al., 2019) and Eurasian crane data (Soriano-Redondo et al., 2019) as well as on processes such as tornadoes (Gómez-Rubio et al., 2015) and pollution spread (Cameletti et al., 2013). More recently, a study has also applied the SPDE approach to Canadian WBPHS data to predict the abundance of 15 waterfowl species (Adde et al., 2020). We used a low-resolution mesh (fewer and larger triangles) in the first stage of analysis to speed processing time as recommended by Krainski et al. (2018) and Bakka (2019).
In models M3-M6, spatiotemporal correlation was represented using SPDE in combination with an autoregressive structure AR1 process for residuals (Zuur et al., 2017). Because we used a Bayesian analysis, the models required priors as starting values. For all fixed effects but the intercept, we used normal priors provided by the INLA package (Rue et al., 2009). For the intercept, we provided a prior with a mean of 0 and precision of 0.001 (Kifle et al., 2017).
We used penalized complexity (PC) priors for the latent effects in our models as recommended by both Simpson et al. (2017) and Fuglstad and Beguin (2018). These priors penalize departure from a base model and encourage parsimony in model selection. We also used information from the early stages of analysis to inform the prior nominal range of the SPDE mesh in final models. The nominal range is the distance at which residual autocorrelation declines to 0.1 (Krainski et al., 2018). We fitted all models using the INLA package (Rue et al., 2009) in the r statistical environment (R Core Team, 2019). We compared the six described model structures using our hold-out dataset and Spearman's correlation test (Humphreys et al., 2019).
The model that provided the highest R-squared values was then used for the second stage of the analysis, in which we applied a remove-one approach to test support for our predictor variables (Chambers, 1992;Walker, Rotella, Schmidt, et al., 2013). In this approach, a variable was removed from the global stage-one model, its Watanabe-Akaike's information criterion recorded, and then the variable put back into the model (WAIC: Gelman et al., 2014;Vehtari et al., 2017). When the removal of a variable decreased the WAIC score of a model by any amount, that variable was not included in the final reduced model. After we applied the remove-one approach to all variables in the model, we ran the reduced model with a high resolution SPDE mesh to acquire parameter estimates.
We assessed the fit of the most supported model from stage 2 using the hold-out data. We compared model-based predictions to actual pair counts using Spearman's correlation test. R-squared values over 0.7 with p-values below .01 were considered to support correlation and model predictive ability.

| Spatiotemporal brood count data
We used data from several previous studies conducted from 2008 to 2010 (Walker, Rotella, Schmidt, et al., 2013), 2012 to 2013 (Carrlson et al., 2018) and from 2014 to 2017 (Kemink et al., 2019) to develop spatially explicit brood abundance models (Figure 3c). Data were not collected during 2011. The data collection for these surveys was conducted at individual wetland basins. Observers surveyed basins either from a vehicle on the roadside or on foot from the edge of the basin. Each basin was visited two to three times in a 36-hr period. Because the models we intended to use did not permit missing response data, and most of our data were collected via two visits per basin, we selected only two visits from surveys with three visits (Walker, Rotella, Schmidt, et al., 2013). We then had early morning (sunrise-12:00) and late afternoon surveys (15:00-sunset) for comparison. More details on data collection can be found in previously published literature (Carrlson et al., 2018;Kemink et al., 2019;Walker, Rotella, Schmidt, et al., 2013).

| Spatiotemporal brood models: predictor variables
We tested the explanatory strength of a suite of covariates that had significant influence on brood abundance in previous analyses (Table 1) Kantrud, 1971). This covariate differentiated between wetlands that were permanent (lakes), experienced strong summer drawdowns (semipermanent), were ponded only through July or August (seasonal) and those that were ponded for only 1-2 months early in the breeding season (temporary: Johnson et al., 2010). We also incorporated several wetland-level variables in the brood detection models.
The detection models were, however, not the focus of the analysis, and we included them largely so that we could ensure abundance estimates were being adjusted for imperfect detection rates (Pagano & Arnold, 2009;Royle, 2004).
Two of the landscape covariates we included in our models we expected to have positive relationships with brood abundance.
Here, we define landscape as a 10.36 km 2 plot on which brood data were collected during the survey. As described previously, we expected higher May pond counts to lead to more breeding duck pairs (Johnson & Grier, 1988) and subsequently greater numbers of broods on surveyed basins. Similarly, we predicted a positive relationship between perennial cover and brood abundance.
Hypotheses regarding this relationship have typically stemmed from the relationship of covariates with pair nesting success (Carrlson et al., 2018;Kemink et al., 2019;Stephens et al., 2005;Walker, Rotella, Schmidt, et al., 2013). In contrast, we predicted that higher amounts of wet area on the landscape in July would provide greater opportunity for birds to spread out, fewer detection opportunities and lower basin-level abundance (Carrlson et al., 2018;Kemink et al., 2019).
During all brood surveys used in our modelling, concurrent flights were used to acquire ponding data on surveyed wetlands and the surrounding landscape. Both technicians and automated software techniques were used in combination to classify the resulting imagery. Specific methodologies can be viewed in previous publications (Carrlson et al., 2018;Kemink et al., 2019;Walker, Rotella, Schmidt, et al., 2013). We used these shapefiles in addition to data collected by observers during the surveys to parameterize the models.

| Spatiotemporal brood models: analysis
We analysed brood count data (2008)(2009)(2010)(2012)(2013)(2014)(2015)(2016)(2017) in two stages. Our main impetus was to minimize processing time because the final models we used would have been temporally prohibitive to run through model selection criteria. Prior to any modelling, we stratified the data by year and randomly split them into training (80%) and testing (20%) datasets.
In the first stage of the analysis, we tested the explanatory strength of our selected predictor variables on the training dataset, modelling data within a maximum likelihood framework using Nmixture models in the r package unmarked (Fiske & Chandler, 2011).
Applying a remove-one approach, we identified variables that increased the model AIC value and earmarked those to be removed from the final reduced model. We used the reduced model in the second stage of analysis.
We modelled brood abundance in the second stage using Bayesian N-mixture intrinsic conditional autoregressive models (iCAR: Besag, 1974), which allowed us to account for both imperfect detection and spatial autocorrelation (Guélat et al., 2018;Latimer et al., 2006;Vielledent et al., 2015). This model combines an ecological process dealing with the abundance of duck broods due to habitat suitability and an observation process that accounts for the probability of detection being less than one (Pagano & Arnold, 2009). Others have used this modelling approach in a similar manner on shorebirds and pintails (Specht, 2018) and on cetaceans (Vilela et al., 2016).
These models treated the true wetland-level abundance (N) as a latent variable with a Poisson distribution and estimated N via a simple reflective random walk algorithm (Hastings, 1970;Vielledent et al., 2015). The observed counts of broods (y) on site i during visit j followed a binomial distribution with index parameter N i and success parameter p ij . The ecological process (Abundance: λ i ,) was modelled through a log link as a function of U covariates and the observation process (detection probabilities) through a logit link as a function of V covariates. The ecological process contained an additional term rho ji to account for the spatial autocorrelation between observations wherein the abundance of broods on one wetland depends on the abundance of the broods on neighbouring wetlands. Here, u j is the mean of ρ j in the neighbourhood of j, V ρ is the variance of the spatial random effects, and n j is the number of neighbours for the spatial entity j. The models were parameterized with flat priors and fitted using the "hSDM" package (Vielledent, 2019) in the r statistical environment (R Core Team, 2019).
We assessed model fit in the second stage by conducting Spearman correlation tests between predicted and actual count values for the hold-out dataset (Humphreys et al., 2019;Kendall, 1938).
We conducted these tests at both the basin and the plot (10.36 km 2 ) resolution because previous analyses have advised that the plot is the best grain for planning with these data and models (Carrlson et al., 2018). Model fit was considered sufficient if correlation values were over 0.70 with p-values less than .01.

| Spatiotemporal model-based predictions
Developing predictions for each year within the time period 2008-2017 required annual PPR-wide layers describing spring and summer ponding as well as overall wetland seasonality. We developed these layers using the Global Surface Water Layer (Pekel et al., 2016). We used layers describing the monthly maximum ponding extent (April-May and July-August) to describe May pond counts (pair and brood models), July wet areas (brood models), basin regime (brood models) and ponded wetland hectares (brood models). We assessed these input variables for accuracy and excluded outliers and data points with missing or invalid predictor data. Other input variables for the pair and brood predictions were obtained from layers used in the original modelling process. In the brood models, the exception to this was the emergent cover variable. Because it was not feasible to obtain region-wide information on the status of this variable, we devel- were scaled by this amount to obtain per km values. This process was completed for each year of the analysis (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017), to obtain 10 raster layers.

| Comparison of distributions
To facilitate comparison to the brood data, we applied similar methods to both the averaged pair data and our pair prediction raster layers. We aggregated the averaged pair abundance data to a 1 km × 1 km raster layer. Then, we applied focal statistics using a 10.36 km 2 neighbourhood to both the averaged pair data layer and the 10 modelled pair abundance layers. We clipped the spatiotemporal pair distributions to the extent of the spatiotemporal brood and averaged pair abundance data. Next, we used the Spearman correlation statistic to test for similarities between the averaged pair, and spatiotemporal pair, and brood data distributions in these areas. All raster comparisons were completed using the stats package in program r (cor: R Core Team, 2019).
For all correlation results, we considered values greater than 0.70 to be significant, indicative of highly similar distributions and to suggest the potential for surrogacy as a conservation measure.
Finally, we examined the overlap among our predicted pair and brood distributions and the averaged pair abundance data (Reynolds et al., 2006). We considered larger proportional areas of overlap to be more indicative of similar distributions and to suggest the potential for surrogacy as a conservation measure. We examined only

| Spatiotemporal pair models
The two stages of our modelling process for the breeding duck pair data provided support for a reduced model with spatial and temporal autocorrelation. In the first stage of our analysis, we found most support for a model structure demonstrating an additive relationship between spatial and temporal autocorrelation (Appendix S1: Table S1). In the second stage, the remove-one analysis showed support for the removal of all variables except adjusted pond count (Appendix S1: Table S2).  (Table 2) and a high spatial autocorrelation with a median nominal range of 78 km (CI: 70-88 km: Appendix S1: Figure   S1).

| Spatiotemporal brood models
Our initial remove-one analysis did not support the removal of any predictor variables in the brood abundance or detection models. Thus, we used a global model in the Bayesian analysis to obtain parameter estimates. Results supported the major conclusions of previous studies, indicating that wetland area is a strong driver of duck brood abundance in the PPR. Further, variables at a larger spatial resolution had both positive (May pond count, perennial cover) and negative (July wet area) associations with brood abundance. However, the credible intervals for the perennial cover relationship crossed zero, suggesting some ambiguity in this effect (Table 3). We also saw support for inclusion of variables describing the seasonality of ponds. The largest difference was between the "Lake" category and the more ephemeral pond types.
Finally, model parameter estimates suggested that abundance varied significantly across years and that spatial correlation was relatively high (Appendix S1: Figure S2).

| Spatiotemporal pair and brood predictions
Using the top models from each analysis, we provided yearspecific pair and brood predictions of abundance for all years 2008-2017, except for year 2011 when no data were collected on broods. Median bootstrapped estimates of pair abundance for the traditional WBPHS area within the US and Canadian PPR varied annually and ranged from 12,188,879 (2008: 95% CI 10,753,552 -13,861,340) to 16,898,488 (2014; Appendix S1: Figure S3a). Predicted distributions at the 10.36 km 2 resolution reflected these temporal changes but did not change dramatically across the study period, with the highest densities of pairs remaining concentrated in the western PPR each year (Appendix S1: Figure S4).

| Distribution comparisons
The strongest correlations occurred between the averaged pair distribution (Appendix S1: Figure S6) and our predicted pair distri-

| D ISCUSS I ON
Our results underscore the contributions that current conservation targeting tools have made to waterfowl conservation to date but also suggest that conservation plans in the PPR would benefit from the additional consideration of intra-and inter-annual dynam- In the PPR, the cyclic weather patterns of drought and deluge drive many of the changes in annual carrying capacity for waterfowl. Evidence of these dynamics has been displayed in other studies (Doherty et al., 2015;Janke et al., 2017;Johnson & Grier, 1988) and which parallels reports of improved pond conditions in that area and period (USFWS, 2008(USFWS, , 2009(USFWS, , 2010(USFWS, , 2011. The averaged pair data identified the Northwest Glaciated Plains as an area that was consistently important for breeding pairs. However, in portions of the Northern Glaciated Plains there were areas where our distributions predicted higher densities than the averaged pair data because of changes in wetland numbers. Model-based predictions of brood abundance also suggested disagreement with the averaged pair data, supporting our hypothesis and the research of others (Carrlson et al., 2018;Talent et al., 1982;Walker, Rotella, Schmidt, et al., 2013)  that targeted areas will be highly dependent on whether an organization's conservation goal is minimizing poor, increasing average or facilitating excellent production in good years. If the latter is true, a conservation strategy that targeted areas with consistently high brood numbers would be most appropriate. However, if the goal was to minimize poor production, areas used less often but during drought years might be equally if not more important because of their value as refugia (Bino et al., 2015;Murray et al., 2012;Stralberg et al., 2020).
Even with the addition of pair and brood spatiotemporal distributions, the efficacy of a conservation prioritization tool for the PPR would depend, in part, on the uncertainty and error accompanying the predictions. The noisy nature of the input data and the questions we asked resulted in uncertainty in our predictions, particularly for the pair data which were modelled at a coarser spatial resolution than the brood data (Hermoso & Kennard, 2012). While we feel the results presented herein are robust given the spatial and temporal resolution of the data used, we also note that the datasets incorporated for developing annual predictive surfaces could be improved. The Global Surface Water layer we used had a 30 × 30 m resolution and was not developed for identifying wetlands obscured by vegetation (Pekel et al., 2016). As a result, we expect that abundance was underestimated in some areas, most likely the abundance of pairs because of their preference for small, temporary and seasonal wetlands Johnson & Grier, 1988;Reynolds et al., 2006). However, preliminary correlation analyses indicated that the layers developed from the data were positively correlated with both May pond counts from the WBPHS and brood survey wetland data. Thus, we were comfortable using these data for predictions and maintain that, until an easily accessible data source at a comparable spatiotemporal scale

F I G U R E 5
is made publicly available, the Global Surface Water data might represent the best option for regional geospatial wetland data in the PPR (Davidson, 2014;Guo et al., 2017). Further, we emphasize the importance of addressing uncertainty in any conservation planning strategy (Langford et al., 2009).

| CON CLUS ION
Waterfowl conservation is perhaps one of the oldest fields of conservation management but has yet to adopt many of the new conservation practices such as the integration of spatiotemporal processes addressed in this analysis (Prairie Pothole Joint Venture, 2017). Future studies will need to improve upon our work here by incorporating better remote sensing data for predictions (more geared towards the PPR), brood data from Canada so that predictions can be expanded to that area and sensitivity analyses regarding uncertainty that include cost data. Foundation. We appreciate the efforts of the numerous seasonal and permanent personnel who collected the brood and pair count data as well as the private landowners who allowed access to their property.

ACK N OWLED G EM ENTS
The authors would also like to thank Gwen Iacona, two anonymous reviewers, Johann Walker, Mason Sieges, Adam Janke, Kyle Kuechle, Sean Fields, and Chuck Loesch for their comments and review.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13218.