A cost efficient spatially balanced hierarchical sampling design for monitoring boreal birds incorporating access costs and habitat stratification.

Predicting and mitigating impacts of climate change and development within the boreal biome requires a sound understanding of factors influencing the abundance, distribution, and population dynamics of species inhabiting this vast biome. Unfortunately, the limited accessibility of the boreal biome has resulted in sparse and spatially biased sampling, and thus our understanding of boreal bird population dynamics is limited. To implement effective conservation of boreal birds, a cost-effective approach to sampling the boreal biome will be needed. Our objective was to devise a sampling scheme for monitoring boreal birds that would improve our ability to model species-habitat relationships and monitor changes in population size and distribution. A statistically rigorous design to achieve these objectives would have to be spatially balanced and hierarchically structured with respect to ecozones, ecoregions and political jurisdictions. Therefore, we developed a multi-stage hierarchically structured sampling design known as the Boreal Optimal Sampling Strategy (BOSS) that included cost constraints, habitat stratification, and optimization to provide a cost-effective alternative to other common monitoring designs. Our design provided similar habitat and spatial representation to habitat stratification and equal-probability spatially balanced designs, respectively. Not only was our design able to achieve the desired habitat representation and spatial balance necessary to meet our objectives, it was also significantly less expensive (1.3-2.6 times less) than the alternative designs we considered. To further balance trade-offs between cost and representativeness prior to field implementation, we ran multiple iterations of the BOSS design and selected the one which minimized predicted costs while maximizing a multi-criteria evaluation of representativeness. Field implementation of the design in three vastly different regions over three field seasons showed that the approach can be implemented in a wide variety of logistical scenarios and ecological conditions. We provide worked examples and scripts to allow our approach to be implemented or adapted elsewhere. We also provide recommendations for possible future refinements to our approach, but recommend that our design now be implemented to provide unbiased information to assess the status of boreal birds and inform conservation and management actions.

Introduction

Tackling ongoing [1-3] and projected [4,5] global biodiversity losses will require difficult decisions about resource allocation and the need to consider where, when, and how to focus conservation efforts. Species conservation is often most successful and least expensive when implemented early, before species require dedicated recovery efforts [6]. Early identification of species declines should also improve the likelihood of successful conservation. In addition, developing effective conservation strategies requires identification of the locations and factors contributing to species’ declines [7]. Sparse or biased data may inaccurately identify which species require conservation actions or misidentify key drivers, resulting in misdirected or ineffective conservation actions [8]. Thus, well-designed ecological monitoring is necessary to prioritize conservation activities effectively [9-11].

Even comparatively well-monitored taxa such as terrestrial birds have significant gaps in species coverage that hinder effective status and trend assessments and associated conservation actions [12]. In North America, many species of terrestrial birds are primarily monitored using the North American Breeding Bird Survey (hereafter BBS; [12]). However, this survey, which is based on roadside surveys and mainly run by volunteers, has very limited coverage in remote northern locations such as the boreal forest [13,14]. Furthermore, much of the coverage for species that breed in the boreal forest is limited to the southern edges of their range, where population change may be quite different from elsewhere. The lack of data from most of the boreal forest may present significant risks to conservation given expanding resource development [15-17] and projected climate change impacts on boreal birds [18].

Monitoring data from northern ecoregions are needed to inform conservation actions (e.g., prioritize selection for protected areas), contribute to management decisions (e.g., adaptive management), detect range shifts, assess threats to species (e.g., species response to developments), and meet legislative requirements. In Canada, monitoring data are required to support conservation under the Migratory Birds Convention Act (1994, c. 22), which includes informing listing decisions under Canada’s Species at Risk Act (hereafter SARA) [19], and environmental impact assessments required by the Impact Assessment Act [20] or by provinces and territories. Key criteria used by the Committee on the Status of Endangered Wildlife in Canada (COSEWIC) to recommend listing a species under SARA include rate of decline in the total number of mature individuals, population size, and extent of species occurrence [21]. In addition, habitat- and region-specific estimates of species’ abundance or density are required as baseline data in environmental impact assessments [20]. To support these diverse requirements, multi-species bird monitoring should be capable of accurately capturing not only trends in species’ population sizes, as well as spatial and temporal variation in species’ status and distribution.

Given the breadth of these objectives and the large and remote areas to be monitored, a cost-effective sampling approach is required. As most of these locations have limited road networks [22] with few human settlements, there is little potential for volunteer surveys [23]. Similarly, a simple, design-based sampling strategy (e.g., a randomized design in which sample units are selected with the same probabilities) is logistically and financially impractical [9,24]. Spatially balanced sampling offers an attractive alternative as it tends to be more efficient, providing more precise estimates than alternative designs [25-28]. Despite the advantages offered by spatial sampling designs, they have not been more widely implemented in large-scale (e.g., biome- or continent-wide) studies, perhaps because they force sampling into locations that can be costly to access. Instead, we propose a design-based sampling scheme that integrates multiple approaches from the statistical sampling-theory literature, including spatially balanced sampling [25,26], stratified sampling [29-32], cost-constrained sampling [31] and optimization [33,34]. The integration of these approaches facilitates simultaneous maximization of both spatial balance (i.e., evenness of the spatial coverage) and habitat representation, while minimizing program costs.

Here, we describe a cost-effective hierarchically structured unequal-probability sampling design known as the Boreal Optimal Sampling Strategy (BOSS). Specifically, our objectives are to: (1) outline the development of the BOSS sampling design; (2) describe BOSS design implementation in three ecologically distinct regions in Canada (Newfoundland and Labrador, Saskatchewan, Yukon Territory); (3) evaluate how the BOSS design meets spatial balance, habitat representation, and cost objectives in these three distinct regions; and (4) evaluate implementation costs in field trials. We compared our design against three competing spatially balanced sampling designs: a purely spatial (equal-probability) design, a habitat stratification design, and a design that minimized access costs. We hypothesized that the BOSS design should provide the second best habitat representation after the habitat stratification design since we incorporated habitat in the inclusion probabilities. We also hypothesized that our BOSS design should rank as the second-least expensive alternative after the cost design owing to incorporation of costs within the inclusion probabilities.

Materials and methods

Sampling frame

We defined the target population as all adult individuals of terrestrial bird populations occupying an ecoregion, province/territory (hereafter jurisdiction), and Bird Conservation Region (BCR) within the boreal region of Canada [35] during the breeding season for a particular year. We defined the Primary Sampling Unit (PSU) as a 5 km diameter hexagonal grid cell (i.e., circumscribed by a circle with a radius of 2500m). We generated a hexagonal grid of PSUs covering the entire sampling frame: the boreal region of Canada, plus a 100 km northern buffer (Fig 1). We used hexagonal geometry because it introduces fewer edge effects and performs well in nearest neighbour analyses since the distances between centroids are the same in all directions [36]. We included a northern buffer to allow the design to capture possible breeding range shifts associated with climatic changes in the north [37]. Each PSU was assigned a unique number/identifier. With the exclusion of PSUs that have cumulative inclusion probabilities of zero (e.g., exclusively open water or having too little land to fit secondary sampling locations), all PSUs in the sampling frame were available to be sampled. As a result, the BOSS sampling design can make valid inferences to entire populations of species inhabiting the boreal region [38].

Fig 1

Spatial extent of the sampling frame of the Boreal Optimal Sampling Strategy (BOSS) design.

The green region represents the boundary of the boreal region as defined in [35] with a 100 km northern buffer within the 4-Northwestern Interior Forest, 6-Boreal Taiga Plains, 7-Taiga Shield and Hudson Plains, and 8-Boreal Softwood Shield Bird Conservation Regions. Reprinted from [35] under a CC BY license, with permission from NRC Research Press, original copyright 2009.

Stratification

We developed a hierarchical stratification scheme to estimate population variables at multiple spatial scales. At the first level, we stratified the sampling frame using the intersection of BCRs and jurisdictions (e.g., area of Yukon Territory within BCR 4-Northwestern Interior Forest). At the second level, we further stratified each BCR by jurisdiction and ecoregion (e.g., the area of Selwyn Mountains ecoregion within Yukon Territory within BCR 4), where ecoregions are a subdivision of an ecozone characterized by similarity in surface forms, flora, fauna, hydrology, soil, and macro climate [39]. Regional monitoring programs as part of the BOSS design are defined by the area of all boreal ecoregions falling within a given jurisdiction. We set target sample sizes for the number of PSUs to select within an ecoregion using the process outlined below (see sample size allocation).

Sample size allocation

From sampling theory, a stratified sample is optimal when it provides the greatest precision for the lowest cost [40]. Optimal stratified sampling accounts for both the size of the stratum (e.g., land area for geographic regions), and the variance within the stratum [40]. By placing greater sampling effort in strata that are larger and more variable, the overall precision of the estimate(s) is optimized [29,31]. Here, we calculated the sample size of PSUs (n) for a given jurisdiction using the following equation to allocate sampling effort between all (k) ecoregions within the jurisdiction: Where n is the target sample size for a given ecoregion within a given jurisdiction, N is the total target sample size to sample within the jurisdiction, Area is the area of ecoregion e, σ is a measure of observed or predicted variation within ecoregion e for the variable(s) of interest, and r is the described species richness of the bird community in ecoregion e (i.e., the number of all bird species with ranges overlapping the ecoregion derived from overlays of species range maps with the ecodistricts of Canada [39] by GeoInsight Corporation [41]. Since our goal is to monitor multiple species, employing a univariate measure of variance would provide an inadequate representation of variance [31]. We included species richness as a weighting factor on within-stratum variation (σ), to increase sampling in ecoregions where more species occur and improve (1) the likelihood of adequately monitoring a greater number of species, and (2) the ability to monitor and predict species distributions.

Many of the ecoregions in our sampling frame have little or no historic sampling from which we could derive estimates of spatial or temporal variance in the bird community. As an alternative, we used proxy variables that correspond with the spatial and/or temporal variance in bird abundance within an ecoregion to set relative PSU sample size targets for the ecoregion (i.e., as proxies for σ). We used these proxy variables [31] to calculate a measure of multivariate dispersion (variance estimate) for each ecoregion following [42,43], representing anticipated spatial and temporal variation in bird abundance and distribution. We used several well-described links between species abundance and distribution and spatial and temporal climatic variation [44], vegetation [44,45] and variation in forest age [45,46]. Specifically, we used the following proxy variables: (1) variance of mean temperatures for May-July across years (1981–2014); (2) variance of total precipitation for May-July across years (1981–2014); (3) standard deviation of annual burn rates (% area burned by wildfire per year [47]) between 1980 and 2015 within an ecodistrict; (4) richness of land cover classes [48] within an ecodistrict; (5) standard deviation of elevation; and (6) percent area of open water since riparian edges significantly increase avian species richness [49]. We calculated variance of mean May-July temperatures and variance of total precipitation from monthly gridded (0.5 x 0.5 degrees) values from the climate datasets of Harris et al. [50]. We calculated the between year standard deviation of burn rates from data obtained from [47]. Since burn rates from [47] were calculated in sample units of 10,000 km2 hexagons, we first calculated the between-year standard deviation of burn rates within hexagons, and derived the mean of those values for ecodistricts. Land cover of Canada 2010 [48] is a 250 m resolution land cover layer of 15 land cover classes from 13,350 Landsat-7 satellite images taken between 2009 and 2011. Ecodistricts (subdivisions of ecoregions characterized by distinctive landforms, relief, geology, soils, water bodies, flora and fauna [39]) were treated as sampling units to allow us to calculate dispersion within ecoregions. Proxy variable data are provided in S1 Data.

Sampling design

We developed a two-stage sampling design in which the first stage of the design involved selection of PSUs from ecoregions within the sampling frame and the second stage involved selection of Secondary Sampling Units (SSU) within the selected PSUs (Fig 2). The SSUs represent on-the-ground locations where point count surveys [51] will be conducted using trained observers and/or autonomous recording units (hereafter ARUs) [52].

Fig 2

The two-stage sampling design of the Boreal Optimal Sampling Strategy (BOSS) design.

This figure illustrates the selection of Primary Sampling Units (PSU, 5 km diameter hexagonal grid cell) from ecoregions within the sampling frame and selection of Secondary Sampling Units (SSU, point count survey locations) within the selected PSUs. Reprinted from [35] under a CC BY license, with permission from NRC Research Press, original copyright 2009.

Selection of primary sampling units

To select PSUs within each jurisdiction (Stage 1), we first developed an approach to combine the following three design components: access costs, habitat representation, and spatial balance. We followed the steps outlined below in each of our three test jurisdictions (Newfoundland and Labrador, Saskatchewan, and Yukon Territory) to integrate all components of our sampling design across our sampling frame (Fig 3). First, we created separate cost surfaces for logistically feasible access methods within each of our jurisdictions. Specifically, we developed spatially explicit models predicting the cost to access a given PSU assuming access from roads (trucks), helicopters, float planes, canoe, motorized boat, snowmobile and/or all-terrain vehicles as appropriate for a given jurisdiction (see examples in S2 Data). The resulting cost models were represented as pixel-based (250 x 250m) estimates of cost, from which we subsequently calculated the average cost of access for all pixels within a PSU separately for each access type (e.g., trucks, helicopters, float planes). We assumed that the lowest-cost access method will generally be used to access a given PSU and calculated the minimum cost across all access cost surfaces. We subsequently used this value as the estimated access cost for each PSU.

Fig 3

Components of the Boreal Optimal Sampling Strategy (BOSS) design depicted for Newfoundland and Labrador, Saskatchewan, and Yukon Territory: (a) access cost; (b) weighted habitat inclusion probability; (c) cost and weighted habitat inclusion probability. Reprinted from [35] under a CC BY license, with permission from NRC Research Press, original copyright 2009.

Second, we calculated inclusion probabilities based on land cover classes (hereafter habitat classes) from the North American Land Change Monitoring System [48]. This land cover product consisted of nineteen land cover classes in total, of which only fifteen occurred in our sampling frame. We calculated habitat-based inclusion probabilities in an attempt to achieve balanced representation of habitat classes within an ecoregion. Inclusion probabilities for habitat classes used a weighted sampling approach to balance sampling across common (high occurrence/high proportion) and uncommon (low occurrence/low proportion) habitats within an ecoregion. Weighted habitat inclusion probabilities were calculated as follows: where the inclusion probability for a given pixel p of habitat class i within a given ecoregion (p) was a function of the inverse of the number of habitat classes within the ecoregion divided by the area of habitat i within ecoregion e. Our weighted habitat inclusion probability thus results in the cumulative probability of sampling a given habitat within an ecoregion being equal among all habitat classes despite differences in habitat area, to help ensure that we can draw inferences for rare habitat types.

Third, we combined access costs with weighted habitat inclusion probabilities to derive our final summed value for inclusion probabilities for each PSU as follows: Where I = the sum of all inclusion probabilities within a given primary sample unit (PSU) 1 through k, P = area-weighted selection probabilities for habitat class i for pixels (p) 1 through j in ecoregion e per Eq 2, and cost for a given PSU was the minimum access cost for sample unit k (as described above). Cost was included as the inverse square root of costs to preferentially (all else being equal) favour sampling at PSUs that are less expensive to sample in accordance with optimal allocation theory [53]. Across all sample units within the given study area, the inclusion probabilities sum to unity based on Eq 3.

Fourth, we used unequal-probability sampling based on the calculated inclusion probabilities (P) to draw a spatially balanced sample within each ecoregion using Generalized Random Tessellation Stratified sampling (GRTS; [25]). In brief, the GRTS approach maps two-dimensional space into a single dimension (that maintains the spatial properties of the original space) and takes a systematic sample along that representation, yielding a sample that is neither over- nor under-represented across space. For a detailed description, we refer the reader to [25] and [30]. The random draw to select PSUs within ecoregions in each jurisdiction uses PSU sample size targets for each ecoregion (where ecoregion boundary was clipped to the jurisdiction boundary). We included a 20% oversample of PSUs for each ecoregion to allow substitution of inaccessible or unsafe PSUs (e.g., cliffs or other features not visible in the available input layers). All PSUs were labeled with their draw order and oversample PSUs were additionally labeled as oversamples. Should a PSU be inaccessible or unsafe to survey, the oversample and draw order can be used to select a replacement PSU that retains the spatial balance of the sample [25].

Finally, we derived an approximately optimal draw of PSUs using a combination of repeated random sampling and multi-criteria evaluation. Specifically, we propose that multiple random unequal-probability GRTS samples can be obtained for which the predicted sampling costs and metrics of habitat representation and spatial balance can be derived. The iteration resulting in the highest combined representativeness of habitat sampling combined with spatial balance for the lowest predicted sampling cost should provide an approximately optimal design. As a metric of habitat representation, we calculated the sum of squared differences (SSD) between the area of each habitat class within all PSUs in a randomized draw versus the proportional area of each habitat class. As a metric of spatial balance, we use Pielou evenness (PE in Eq 4 below) to represent the even spread of PSUs across the sample frame [30]. Possible values of PE range from 0–1. In order to evaluate which draw represents the best trade-off of cost versus combined spatial and habitat representation across multiple (1…j) iterations, we propose evaluating which iteration minimizes predicted costs while maximizing the weighted sum [34] of the two metrics: Where SSD for sample draw i is relativized to the minimum SSD across draws 1 through j, by subtracting the resulting value from 1 so that the metric scales from 0–1 and thus a value of 1 would represent idealized representation of habitat composition within the draw. We applied equal weights to both the habitat (w) and spatial balance (w) components as we did not want either component to be more heavily valued than the other in our final design (i.e., w = w = 0.5). However, this approach can be generalized by allowing weights to vary between habitat representation and spatial representation if there is a desire to emphasize habitat representation (e.g., for environmental impact assessments) or spatial representation (e.g., for distribution modeling).

Selection of secondary sampling units

To select secondary sample units (SSU) within each PSU in each jurisdiction (Stage 2), we first created a uniform (systematic) grid of potential secondary sample unit locations at 300 m spacing over the extent of each ecoregion. Each SSU was separated by 300 m since the effective detection radius for most boreal songbirds is less than 150 m [51] and therefore double counting should be infrequent for most species. We used the grid of SSUs to query the habitat-based inclusion probabilities. Within each PSU, we selected a stratified random sample of four SSUs using inclusion probabilities as calculated in Eq 2. We used the four SSUs selected during this draw as plot centroids and the eight SSUs surrounding the plot centroid to construct plots of nine SSUs for sampling (see Fig 2). We set a sampling target of two plots of nine SSUs per PSU so that a team of two staff can complete one PSU per day in most circumstances. Furthermore, having two survey plots within a PSU allows staff to have a partner(s) to provide support in an emergency situation. Our PSU size of 5 km diameter was specifically selected to ensure field staff would be working within a reasonable walking or flight distance of a partner if assistance or emergency response is required. We drew a sample of four plots of SSUs per PSU to create an oversample in case of inaccessibility (e.g., due to large water barrier or cliff) or safety concerns (e.g., hazardous terrain or aggressive wildlife).

Data collection

We tested the BOSS design within three jurisdictions across the boreal biome of Canada (Newfoundland and Labrador [NL], Saskatchewan [SK], and Yukon Territory [YT]) during 2017, 2018, and 2019. Newfoundland and Labrador is located at the eastern edge of the boreal region within BCR 7-Taiga Shield and Hudson Plains and BCR 8-Boreal Softwood Shield (Fig 1). Saskatchewan is located in the center within BCR 6-Boreal Taiga Plains, BCR 8, and BCR 7 (see Fig 1); while the Yukon Territory is located at the northwestern edge within BCR 4-Northwestern Interior Forest, BCR 6, and BCR 3-Arctic Plains and Mountains. Jurisdictions were also ecologically distinct, differing in their geology, topography, vegetation, terrestrial bird diversity, and climate as indicated by the number of ecozones (NL: 2; SK: 3; YT: 3), ecoregions (NL: 23; SK: 8; YT: 22), and elevation range (NL: 0−1,652 m; SK: 204−823 m; YT: 0−5,900 m). Differences in the location of mountain ranges, rivers and lakes, human settlements, and primary and secondary road networks within each jurisdiction present distinct logistical challenges and therefore useful comparative test cases. To conduct point count surveys at SSUs within selected PSUs, we used10-minute point counts with trained observers following standard protocols recommended by Ralph et al. [54] and Matsuoka et al. [51] and conducted during suitable weather conditions from 30 minutes prior to sunrise until 4–5 hours after sunrise, from the last week of May to the first week of July. Where we conducted surveys using ARUs, we analysed 10-minute recordings from [52] to match the count duration conducted by human point counts. ARUs were deployed on a variety of different schedules and using multiple access methods. In some cases, ARUs were deployed in February by snowmobile, and were retrieved in the summer. Wherever possible, we attempted to get recordings conducted over ≥ 4 mornings with good survey conditions; however, ARUs carried on canoe trips were deployed for as little as a single night. We programmed ARUs to record a minimum of six 10-minute intervals over the course of a morning during the same time-period as point counts, as well as additional times for other objectives. An analysis of the relative merits of different programming schedules and a comparison of field observers with ARUs is beyond the scope of this paper.

We sampled 295 PSUs selected using the BOSS design between 2017−2019 in Newfoundland and Labrador (n = 69), Saskatchewan (n = 103) and Yukon Territory (n = 123). Three PSUs in Saskatchewan had to be replaced with PSUs from the oversample owing to access and safety issues associated with two open-pit uranium mines and one on a bombing range within Department of National Defence lands. No oversamples were required in either Newfoundland and Labrador or Yukon Territory. In Newfoundland and Labrador, we sampled 5 out of 9 ecoregions on the island of Newfoundland, with all ecoregions being in BCR 8. PSUs were sampled in the Avalon Forest ecoregion (n = 13), Central Newfoundland (n = 1), Maritime Barrens (n = 8), South Avalon-Burin Oceanic Barrens (n = 23), and Southwestern Newfoundland (n = 24). We sampled seven out of eight ecoregions within Saskatchewan: Selwyn Lake Upland (BCR 7, n = 12), Tazin Lake Upland (BCR 7, n = 9), Athabasca Plain (BCR 8, n = 19), Churchill River Upland (BCR 8, n = 22), Mid-boreal Uplands (BCR 6, n = 28), Mid-boreal Lowland (BCR 6, n = 5), and the Boreal Transition (BCR 6, n = 8). We sampled ten out of 22 ecoregions within Yukon Territory, sampling in the British-Richardson Mountains (BCR 3, n = 1), Old Crow Basin (BCR 4, n = 9), Old Crow Flats (BCR 4, n = 2), Eagle Plains (BCR 4, n = 3), North Ogilvie Mountains (BCR 4, n = 5), Mackenzie Mountains (BCR 4, n = 8), Klondike Plateau (BCR 4, n = 72), Yukon Plateau Central (BCR 4, n = 10), Yukon Plateau North (BCR 4, n = 10), and Ruby Ranges (BCR 4, n = 3).

Ethics statement

We obtained the required permits to conduct avian monitoring to test the BOSS design from provinces/territories, Provincial and Territorial Parks, Parks Canada Agency, and Yukon First Nations with Traditional Territories and Settlement Areas. Field observers spend <15 minutes at SSU locations during observer surveys and <30 minutes (time for set up and subsequent retrieval) at SSU locations during ARU surveys, resulting in minimal disturbance to breeding birds.

Statistical analyses

Prior to analysis of our multivariate proxy data, we centered and standardized all variables to zero means and unit variance and then calculated Euclidean distances between all points (ecodistricts) in the dataset. We derived the multivariate estimate of dispersion by calculating mean distances to ecoregion centroids using the ‘betadisper’ function in the vegan package [55] in the R statistical computing environment [56]. The resulting estimates of dispersion for all ecodistricts within the boreal region, as per Brandt [35], are available in S3 Data.

In order to assess our BOSS design, we used the ‘spsurvey package’ [57] within the R statistical computing environment [56] to draw samples under four alternative spatially balanced sampling designs. Specifically, we drew samples under one equal-probability spatially balanced design and three unequal-probability designs: (a) a habitat stratification design, calculated as per Eq (2) above; (b) a cost design in which inclusion probabilities were based solely on access costs as per Eq (3); and (c) our BOSS design in which inclusion probabilities included both habitat stratification and access costs as per Eq (4) above. For each scenario, we ran 100 iterations in which we used the GRTS algorithm [57] to draw a random sample of n = 400 PSUs from each of the three jurisdictions considered here. Within each iteration, our script calculated the access costs (CDN $), spatial balance (Pielou evenness [30,57]), and habitat representation as SSD (see above). All implementations are available within scripts in the (S2 and S6 Datas) along with geospatial data examples (S4, S5 and S7 Datas).

We used linear models to test for differences (for each jurisdiction separately) in each of our response variables (cost, SSD, and Pielou evenness) between the alternative designs described above (i.e., spatial, habitat, cost, and BOSS). Study design was included as a factor, and in each model we specified an a priori reference factor level to test our hypotheses. We specified the cost design as the reference factor level for models examining access costs, as we anticipated this design should have the lowest access costs on average. Similarly, we specified the habitat stratification design as the reference level for models examining variation in SSD as we predicted this design should have the lowest SSD.

All four designs we considered employed a GRTS algorithm to generate spatially balanced samples; however, the use of unequal inclusion probabilities based on spatially structured variables could alter the ability of the design to achieve spatial balance. We therefore tested for differences in spatial balance between our designs and we treated the spatial design (equal-probability spatial sampling design) as the reference factor level when comparing spatial balance between sampling designs. We calculated robust standard errors (‘sandwich’ estimators) to overcome problems with heteroscedastic errors using the ‘lmtest’ [58] and ‘sandwich’ [59] packages within the R statistical computing environment [56]. All data are available in the (S8 Data).

Results

Sample size allocation

Sample size allocation to geographic strata (ecoregions × jurisdiction) under the BOSS design resulted in a similar sample size on average as allocating samples based on stratum area, because data are centered along a 1:1 correspondence line (Fig 4). However, allocation to individual strata differed markedly, with the BOSS design suggesting significantly lower sample sizes for some strata than an area-based allocation, and allocating disproportionately larger sample sizes in other strata (Fig 4). For example, sampling proportional to area would allocate 151 PSUs to the New Québec Central Plateau ecoregion of Québec, compared to only 50 based on the BOSS allocation, because the stratum has lower than average dispersion (i.e., 0.36 vs average of 1.04). In contrast, 103 PSUs would be allocated to the Abitibi Plains ecoregion of Ontario under sampling proportional to stratum area, but 173 PSUs were allocated to the same stratum under the BOSS design owing to higher than average dispersion (i.e., 20.1 vs. average of 1.04). Several of the larger and more variable strata (sample sizes > 100) require more sampling under the BOSS design than predicted solely based on stratum area (Fig 4). Overall, the BOSS sample size allocation generally results in sampling intensities being greatest in the south and decreasing northward (Fig 5).

Fig 4

Allocation of sampling effort between boreal ecoregions of Canada under the Boreal Optimal Sampling Strategy (BOSS; see Methods) versus traditional allocation based on Stratum areas to achieve a total sample size of 4980 primary sample units (i.e. 2% sample).

Dashed line indicates the 1:1 correspondence line.

Fig 5

Geographic distribution of sampling intensity (number of allocated primary sampling units per 1000 km2) across the boreal ecoregions of Canada under the Boreal Optimal Sampling Strategy (BOSS) design (see Methods).

Reprinted from [35] under a CC BY license, with permission from NRC Research Press, original copyright 2009.

Predicted sampling costs

Predicted costs to access 400 primary sample units varied by the jurisdiction and sampling design considered (Fig 6). Across all three jurisdictions, designs in which the inclusion probabilities were derived solely based on access costs (cost design) were the least expensive designs (Fig 6) and all the alternative designs were more expensive regardless of jurisdiction (Table 1). The BOSS design had the next lowest cost (Fig 6), with costs predicted to be between 1.02 (Yukon Territory) and 1.2 (Saskatchewan) times more expensive on average than the cost design for a given jurisdiction (Table 1). In both Newfoundland and Saskatchewan the equal-probability spatial design had the third-lowest cost followed by the habitat stratification design, whereas it was the most expensive design for Yukon Territory (Fig 6); it was predicted to be 1.4−2.6 times more expensive on average (~$362,000−562,000; Table 1) than the cost design (maximum of 1.5−3.0 times; Fig 6). The habitat stratification design was predicted to be 1.4−2.0 times more expensive on average (~$122,000−1,020,000; Table 1), or a maximum of 1.6−2.3 times more expensive across draws and jurisdictions (Fig 6).

Fig 6

Variation in predicted costs (in Canadian dollars x 106) of accessing 400 primary sampling units in each of Newfoundland and Labrador, Saskatchewan and Yukon Territory.

Data are from 100 separate random draws of 400 primary sampling units in each jurisdiction. Note differing y-axis scales due to variation in overall costs of operating in each jurisdiction. Dark solid line represents the median, box indicates the inter-quartile range, whiskers are 1.5 times the inter-quartile range, and dots indicate extreme values.

Table 1

Linear model results comparing variation in access costs ($) between alternative monitoring designs.

The cost design was the reference category (Intercept).

Design	β	SE	t-value
Newfoundland
Intercept	1,006,006.0	2870.2	350.51
BOSS design	121,535.7	3951.3	30.76
Spatial design	362,012.3	4072.3	88.90
Habitat design	482,229.8	4023.7	119.85
Saskatchewan
Intercept	1,053,370.5	4927.6	213.77
BOSS design	175,634.7	7288.8	24.10
Spatial design	561,151.8	7587.2	73.96
Habitat design	1,020,042.6	8439.4	120.87
Yukon
Intercept	331,627.1	990.8	334.70
BOSS design	8700.2	1615.1	5.39
Spatial design	541,671.3	1658.7	326.57
Habitat design	121,812.7	1507.5	80.80

Habitat representation

As predicted, the habitat stratification design had the lowest sum of squared differences between the habitat areas represented within sample draws and the desired equal habitat representation (Fig 7). The BOSS design provided the second-most representative sampling of habitat classes on average across all three jurisdictions (Fig 7 and Table 2). In contrast, both the equal-probability spatial design and the cost design (Fig 7) had substantially lower habitat representation relative to the habitat stratification design in all three jurisdictions (Table 2).

Fig 7

Variation in habitat class representation (measured as the sum of square differences between the area of each habitat class within all PSUs in a randomized draw versus the proportional area of each habitat class, see Materials and Methods) between alternative sampling designs in Newfoundland and Labrador, Saskatchewan and Yukon Territory.

Data are from 100 separate random draws of 400 primary sampling units in each jurisdiction. Dark solid line represents the median, box indicates the inter-quartile range, whiskers are 1.5 times the inter-quartile range, and dots indicate extreme values.

Table 2

Linear model results comparing variation in habitat representation (i.e., sum of squared differences between area of land cover classes represented within sample draws against a balanced sample (see Methods)) between alternative sampling designs.

The habitat design was the reference category (Intercept).

Design	β	SE	t-value
Newfoundland
Intercept	0.06	3.84 x 10−4	158.59
BOSS design	4.90 x 10−3	5.60 x 10−4	8.75
Spatial design	0.11	9.97 x 10−4	108.35
Cost design	0.10	8.97 x 10−4	114.37
Saskatchewan
Intercept	0.04	2.06 x 10−4	215.01
BOSS design	9.36 x 10−4	3.10 x 10−4	3.02
Spatial design	0.01	4.02 x 10−4	32.31
Cost design	0.02	3.93 x 10−4	41.56
Yukon
Intercept	0.10	4.78 x 10−4	211.13
BOSS design	1.53 x 10−3	7.17 x 10−4	2.14
Spatial design	0.04	6.47 x 10−4	61.51
Cost design	0.04	8.40 x 10−4	53.36

Spatial representation

All sampling designs provided relatively balanced spatial sampling of all three jurisdictions (range of Pielou evenness across all designs: 0.96−1.01; Fig 8). Interestingly, the spatial balance of the cost design provided the closest approximation to the equal-probability spatial design (Fig 8) and with substantial overlap of the data distribution compared to the equal-probability spatial design (Table 3). The BOSS design showed small differences in spatial balance relative to a purely spatial design in both Newfoundland and Yukon Territory (Fig 8 and Table 3), but not Saskatchewan (Table 3). The habitat-stratification design had the largest overall differences in spatial balance relative to the equal-probability spatial design, with the largest difference occurring in Yukon Territory (Fig 8 and Table 3).

Fig 8

Variation in spatial balance (Pielou evenness) between alternative sampling designs in Newfoundland, Saskatchewan and Yukon Territory.

Data are from 100 separate random draws of 400 primary sampling units in each jurisdiction. Dark solid line represents the median, box indicates the inter-quartile range, whiskers are 1.5 times the inter-quartile range, and dots indicate extreme values.

Table 3

Linear model results examining variation in spatial balance (Pielou evenness) between alternative sampling designs.

The equal-probability spatially balanced sampling design was the reference category (Intercept).

Design	β	SE	t-value
Newfoundland
Intercept	0.99	6.52 x 10−4	1514.78
BOSS design	8.70 x 10−3	8.80 x 10−4	9.87
Cost design	1.56 x 10−3	9.03 x 10−4	1.72
Habitat design	4.18 x 10−3	8.92 x 10−4	4.68
Saskatchewan
Intercept	0.99	3.32 x 10−4	2977.07
BOSS design	-1.10 x 10−3	4.23 x 10−4	-2.61
Cost design	9.22 x 10−4	4.52 x 10−4	2.04
Habitat design	-1.39 x 10−2	4.41 x 10−4	-31.61
Yukon
Intercept	1.00	3.12 x 10−4	3204.27
BOSS design	-1.59 x 10−2	5.81 x 10−4	27.41
Cost design	-6.15 x 10−3	5.69 x 10−4	-10.82
Habitat design	-1.21 x 10−2	5.89 x 10−4	-20.68

Optimal design selection

Across 100 random sample draws, the BOSS design provided the best combination of cost savings and both spatial and habitat representation amongst the four designs considered. Regardless of jurisdiction, the BOSS design had costs competitive with the cost design (Figs 6 and 9), with no loss of habitat or spatial representation. For all three jurisdictions combined, selecting the draw that simultaneously maximized representation while minimizing cost (see red-filled squares in Fig 9) only results in an estimated savings of approximately $2,000 relative to the median of all BOSS sample draws. However, access costs for 400 PSUs varied substantially between draws (Newfoundland: $1,076,897−1,199,693; Saskatchewan: $1,089,811 −1,348,550; and Yukon Territory: $309,497−371,174) and thus relative to the most expensive draws for each jurisdiction, the combined savings could be up to an estimated $225,000 for these three jurisdictions when using the draw providing the best trade-off based on our multi-criteria evaluation.

Fig 9

Between draw variation in combined spatial and habitat (weighted sum, see Methods) representativeness of sampling relative to access costs in Newfoundland and Labrador, Saskatchewan and Yukon Territory for four alternative sampling designs: Cost based sampling design (open red triangles), the BOSS design (open black squares), the habitat design (open green triangles), and an equal probability spatial sampling (blue filled triangles).

The red-filled square within the cluster of points representing the BOSS design represents the lowest cost solution that maximizes combined spatial and habitat representation. Note differing x-axis scales due to variation in overall costs of operating in each jurisdiction.

Field implementation

Of the 295 PSUs sampled across all three jurisdictions, we sampled 85 PSUs in remote areas that we predicted would require charter flights to access. In Newfoundland, cost models accurately predicted an average cost of $1,721 (SD = $86) per PSU to access 16 remote PSUs, compared to actual costs which averaged $1,721 (SD = $441) per PSU (Fig 10). In comparison, the cost models for Saskatchewan predicted an average access cost of $3,741 (SD = $1,928) for the 43 remote PSUs sampled, but the actual access costs were substantially lower (mean = $1,573, SD = $1,562). In Yukon, the mean predicted access cost to reach the 26 remote PSUs sampled was $1,574 (SD = $82; Fig 10) per PSU, but the mean cost to access these PSUs was $3,016 (SD = $132). Weighted least squares regression (weights proportional to sample size) suggests across all three jurisdictions the access actual costs were 57% of the model predicted cost on average (β = 0.57, SE = 0.10; Fig 10).

Fig 10

Predicted versus realized access costs (in Canadian dollars) for each of Newfoundland and Labrador, Saskatchewan and Yukon Territory.

The solid line represents weighted (weights proportional to sample size) least squares fit and dashed lines represent 95% confidence intervals. Points represent average costs to access multiple (n = 1–17) primary sampling units, in line with how charter services were billed. Numbers beside symbols indicate the number of primary sampling units included with the associated mean cost estimate.

Discussion

We integrated key concepts from sampling theory to develop a hierarchically structured survey design that can provide stratified sampling of predefined strata (e.g., ecoregions) and achieve the desired spatial balance and habitat representation while minimizing access costs. The probability-based BOSS design allows for valid statistical inference from field samples across the sampling frame and reduces operational costs. We found that our integrated design was significantly less expensive than all but a design based only on access cost. Predicted access costs were 1.3−1.7 times more expensive for the habitat stratification design and 1.6−2.7 times more expensive on average for the spatial design compared to our BOSS design. Additionally, the distribution of habitat representation across sample draws from the BOSS design overlapped almost entirely with that from the habitat stratification design (see Fig 7). While we did find minor differences in the distribution of spatial balance between the designs considered (Fig 8), spatial balance was always close to a value of one, indicating the samples were well spread over jurisdictions irrespective of the type of design. Across jurisdictions, even the spatial design did not always produce spatial balance metrics that overlapped one, suggesting that the inclusion of unequal sampling probabilities in the BOSS design did not introduce systematic biases in spatial representation. Combined with our randomizations to balance trade-offs, we therefore simultaneously achieved habitat and spatial representation goals while minimizing access costs.

Our pilot field seasons suggest that logistic considerations might result in access costs differing from those suggested by our cost models. While cost models for Newfoundland have thus far provided accurate prediction of average access costs per PSU, our cost models underestimated access costs in Yukon. Overall, regression analysis suggested that field implementation was 43% cheaper than suggested by our cost models. Many reasons likely contributed to the cost differences we observed, including external logistical support. One key difference is that our cost models largely portray the cost of access to a sample unit as independent of access to other nearby sampling units. In practice, we were able to combine logistics and access methods for many sampling units and thereby substantially lowered the average cost per sample unit. For example, in the summer of 2019 we chartered a Twin Otter to drop off two crews on or near opposite sides of Cree Lake, Saskatchewan (57°23′57″N 106°40′10″W) to access six PSUs via canoe. During crew retrieval, an additional crew was dropped off to canoe down river to a roadside pickup location, allowing access to two additional PSUs. We were therefore able to access eight PSUs for the cost of two charter flights (the base cost for accessing a single PSU), substantially reducing the average cost per PSU. We therefore recommend that programs revise their cost models as field implementation improves knowledge of local logistics. As a result, final program costs can be lowered by refinement of preliminary cost models.

Our sampling design has several key differences from other large-scale terrestrial bird monitoring programs that improve its utility for remote areas such as the boreal forest. Other programs such as the Integrated Monitoring in Bird Conservation Regions (IMBCR) program [9], the UK Breeding Bird Survey [60], the North American BBS [12] and the US National Park Service Vital Signs Monitoring program [61] also have defined sampling frames and randomly selected sample units that facilitate valid statistical inference from field samples [62]. Importantly, our design has fewer restrictions on the sampling frame while minimizing sampling costs, and implicitly accommodates health and safety considerations. While other programs restrict the sampling frame to roadsides or trail networks to facilitate access and alleviate safety concerns, we instead dealt with access by incorporating access costs into sample unit selection, and designing surveys to include multiple staff within primary sample units. Thus, our design increases safety and reduces costs without limiting statistical inference about bird populations to accessible areas [62-64]. In addition, our hierarchical stratification with static strata (BCRs, jurisdictions and ecoregions) allows for assessment of the population-level responses of birds to natural and human disturbance and climate change [9,65]. Indeed, Pavlacky et al. [9] recently demonstrated the ability of their hierarchical spatial design to allow estimation of both strata-specific population sizes and species-habitat relationships. Our approach of integrating both cost and weighted habitat-inclusion probabilities within a similar spatially balanced design will allow us to draw similar inferences to Pavlacky et al. [9] at multiple spatial scales, but at reduced costs compared to alternative approaches.

Our design directly incorporated several factors into the stratification and layout of primary and secondary sampling units to address the need to minimize costs and provide for the health and safety of field staff. Field implementation of the optimal designs selected here during pilot field seasons across all three jurisdictions showed that the BOSS design was achievable on the ground. We were able to apply the design across a wide variety of challenging conditions from mountainous terrain (e.g., >2000 m elevation), wetlands (floating bogs, fens and marshes), forests with blow-down (recent burns, tornados and insect outbreaks), as well as steep, rocky, and densely vegetated slopes across all three jurisdictions. Furthermore, drawing an oversample using the GRTS algorithm allowed us to replace unsafe PSUs without altering the representativeness of the sample [25]. Using pairs of observers to conduct surveys within each PSU means that staff can work alone to each survey a plot of nine SSUs within a PSU or they can work together to complete one plot of nine SSUs before travelling together to complete the second plot. The latter option would be appropriate if a PSU contains dangerous terrain, wildlife or limits to communication between staff (radio, satellite transmitter).

One complication imposed by our design is that most analyses assume that sampling is proportional to habitat availability and, by extension, population sizes. We used habitat stratification to ensure representation across the range of available habitat classes and thereby provide reasonable sample sizes to facilitate modeling abundance and distributions of species associated with rare habitats. Habitat stratification should also improve our future capacity to detect trends of species using rare habitats. In particular, our approach to habitat stratification should increase model precision for species using rare habitats at the cost of reducing precision for common species, a trade-off we feel is worthwhile because estimating trends, distribution and habitat relationships for rare species is always more limiting. However, extrapolating from the sample to estimate stratum-specific population sizes and/or population trends is more complicated because strata boundaries and sizes can change with shifts in vegetation through time [10]. This problem can be addressed analytically via post-hoc stratification [66], weighting via the inclusion probabilities [68,69], including random-effects for the strata indexed by the inclusion probabilities [70], or inclusion of interval-specific covariates reflecting habitat supply [69]. Recent spatially explicit approaches to avian population size estimation [71] provide one potential approach to build trend estimation models that could implicitly incorporate changes in strata sizes through time.

While the BOSS design presented here incorporates several key improvements over current large-scale bird monitoring designs, we foresee several areas where further improvements could occur. First, we stratified habitat using coarse satellite-derived products due to their availability and consistency across the entire sampling frame. Incorporating higher-resolution products capable of further stratifying habitat, and potentially better quantifying the amount and distribution of rare habitats, should improve our capacity to monitor species associated with these habitats. One possible approach worth investigation would be to create an avian habitat classification (sensu [45]) from higher-resolution data, such as forest inventory data, for the portion of the sample frame for which they may be available. Higher-resolution data could be integrated into the current framework by (1) using the coarse-resolution satellite data to select PSUs and using the high-resolution habitat classification to stratify SSUs within the selected PSUs, or (2) splitting the selection of PSUs between areas with versus without access to these high-resolution data products.

Second, once sufficient on-the-ground data are available, it would be desirable to re-examine the stratification by using the design-based sample to estimate spatial and temporal variance in bird community composition directly for optimal allocation, rather than the proxy variables used in the initial sample. We anticipate that sampling efficiency will improve by incorporating these direct estimates of variability into the design, but overall costs could be inflated if these data are not incorporated until they are available for all regions. Using data from preliminary rounds of sampling under the sampling design to validate allocation based on multivariate dispersion of proxy variables will ensure allocation between strata is cost effective, and could be used for second-phase sampling to improve the precision of estimates further [28].

Future simulation studies are needed to examine the likely trend precision that can be achieved under the BOSS design under alternative PSU revisit schedules (e.g., every 5, 10, or 20 years). Evaluating whether the BOSS emphasis on greater sampling in southern regions (Fig 5) could impede our ability to detect trends for northern-distributed species or potentially improve our capacity to detect anthropogenic impacts on avian populations due to the concentration of anthropogenic disturbance in southern regions would be an important addition to these simulation studies. Finally, optimizing the use and allocation of humans versus ARUs within our design is also a priority for future simulation modelling. Recent studies suggest that biases in species detection between human and ARU surveys can largely be controlled [72,73], but potential biases between single versus multiple visits would ideally be controlled by spatially balanced allocations of survey effort between humans and ARUs. However, logistical constraints and/or the aggregation of hard to detect species (e.g., secretive marsh birds) into a subset of habitat may make other approaches more efficient and is the subject of ongoing work.

To avoid inefficiencies, there is a clear need for improvements in cost models or habitat layers to be reflected in revised sample draws; however, revised sample draws should account for existing sampling under the BOSS design to avoid inducing spatial imbalance. Foster et al. [26] proposed using a squared-loss distance metric to alter the inclusion probabilities around “legacy” monitoring sites (i.e., pre-existing sampling locations that are known to be a randomly selected representative sample from a known sample frame [26]). This approach would provide a clear way to improve future sampling under the BOSS design without introducing biases and inefficiencies due to spatial imbalances. In addition to data from legacy sites, a large collection of historic data exist within the Boreal Avian Modelling (hereafter BAM) project database [74] for a large portion of boreal Canada. The difficulty with incorporating data from the BAM database is that these largely represent “iconic sites” (sensu [26]) which are either known to be non-randomly selected and/or there is insufficient documentation to know whether the sampling frame included the full suite of habitat classes. Including iconic sites in a randomized draw could introduce bias into the design and result in biased status and trend estimates that may not match those in the population [26]. Since these historic datasets could potentially add valuable ecological insights into factors affecting populations of boreal birds, it would be fruitful to consider when and how they could be included within a long-term monitoring design.

It is worth nothing that our approach to “optimizing” which randomized draw to use derives from the concept of Pareto optimality [34], but does not strictly meet these criteria. In optimization, a Pareto optimal solution is a solution in which the results cannot be modified without resulting in the objective function becoming worse. In our approach, this could only be accomplished by considering all permutations of a given sample size. Considering all permutations would become computationally difficult as well as result in a solution that is no longer truly random. As such, we feel our approach provides a reasonable approximation of optimality while maintaining the desirable feature of representing a randomized draw.

We have demonstrated that including access costs and habitat within a hierarchically structured sampling scheme was the most cost-efficient design that maintained the desired spatial and habitat representation. Although we have focused on a design for monitoring boreal birds, our approach should be broadly applicable and adaptable as a template for other regions where remote locations form a large proportion of the sampling frame. We recommend that our approach be adapted for use in other large-scale studies such as Breeding Bird Atlases [75] if access and sampling costs are limiting for much, if not all, of the study areas. For example, the Saskatchewan Breeding Bird Atlas (https://sk.birdatlas.ca/) is currently using systematic sampling within accessible areas, and the BOSS design within the boreal portion of the province. To facilitate adaptation and implementation of our approach and/or further refinements to our method, we provide worked examples and scripts within the Supporting Information. Using design-based approaches to inform the design of large-scale studies would add to their statistical rigor and efficiency and facilitate their integration (i.e., treated as legacy sites) into long-term monitoring to provide more responsive conservation efforts in light of widespread declines in biodiversity [1-3]. Our hierarchical sampling design should be widely implemented to monitor boreal birds because it can provide an unbiased representation of when and where conservation and management should be targeted [2, 76]. Using monitoring data to inform conservation decisions is crucial, given ongoing and projected ecological changes in the boreal biome [18,37] and population declines already observed for boreal birds [3].

Data used as proxy variables for assumed spatial and temporal variation in boreal forest bird communities.

(CSV)

Click here for additional data file.

R script implementing primary and secondary sample unit selection, including an example calculation of access costs for Northwest Territories, Canada.

(R)

Click here for additional data file.

Allocation of 4980 primary sampling units amongst geographic strata (intersection between ecoregions and jurisdictions) across Canada.

(CSV)

Click here for additional data file.

zip file containing land cover data for habitat stratification derived from Latifovic et al.

(Latifovic, et al., 2012) used in association with R scripts in S2.

(ZIP)

Click here for additional data file.

Shapefiles required to run example analyses in script S2.

(ZIP)

Click here for additional data file.

R script implementing parallel computing for primary and secondary sample unit selection for Manitoba, Canada, where access-cost models were created using ArcGIS using the Spatial Analyst™ extension.

(R)

Click here for additional data file.

Zip file containing example geospatial data, metadata describing of cost modeling, and supporting R scripts used in association with R scripts in S6.

(ZIP)

Click here for additional data file.

Comma-separated values file of data from 100 iterations of drawing 400 primary sample units for three Canadian jurisdictions (Newfoundland & Labrador, Saskatchewan and Yukon Territory).

(CSV)

Click here for additional data file.

6 Mar 2020

PONE-D-19-35579

A cost efficient spatially balanced hierarchical sampling design for monitoring boreal birds incorporating access costs and habitat stratification

PLOS ONE

Dear Mr Van Wilgenburg,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

**From the Academic Editor: I apologize for the delay in my initial decision, but I wanted to get comments from multiple reviewers regarding your manuscript, which covers a highly relevant topic that should be of interest to many PLOS ONE readers. The reviewers agree that the manuscript is very well written. Each reviewer has noted just a few aspects that could be improved. I ask that you consider their comments and revise your manuscript accordingly. I don't think this will require too much of your time, and should result in an even stronger paper.**

We would appreciate receiving your revised manuscript by Apr 20 2020 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Frank H. Koch, PhD

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements:

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at http://www.plosone.org/attachments/PLOSOne_formatting_sample_main_body.pdf and http://www.plosone.org/attachments/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. We note that you have indicated that data from this study are available upon request. PLOS only allows data to be available upon request if there are legal or ethical restrictions on sharing data publicly. For more information on unacceptable data access restrictions, please see http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions.

In your revised cover letter, please address the following prompts:

a) If there are ethical or legal restrictions on sharing a de-identified data set, please explain them in detail (e.g., data contain potentially sensitive information, data are owned by a third-party organization, etc.) and who has imposed them (e.g., an ethics committee). Please also provide contact information for a data access committee, ethics committee, or other institutional body to which data requests may be sent.

b) If there are no restrictions, please upload the minimal anonymized data set necessary to replicate your study findings as either Supporting Information files or to a stable, public repository and provide us with the relevant URLs, DOIs, or accession numbers. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories.

We will update your Data Availability statement on your behalf to reflect the information you provide.

3. Thank you for stating the following in the Acknowledgments Section of your manuscript:

"This research was supported by Environment and Climate Change Canada."

We note that you have provided funding information that is not currently declared in your Funding Statement. However, funding information should not appear in the Acknowledgments section or other areas of your manuscript. We will only publish funding information present in the Funding Statement section of the online submission form.

Please remove any funding-related text from the manuscript and let us know how you would like to update your Funding Statement. Currently, your Funding Statement reads as follows:

"The author(s) received no specific funding for this work."

4. We note that Figures #1-3 and 5 in your submission contain map images which may be copyrighted. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.

We require you to either (1) present written permission from the copyright holder to publish these figures specifically under the CC BY 4.0 license, or (2) remove the figures from your submission:

1.    You may seek permission from the original copyright holder of Figures #1-3 and 5 to publish the content specifically under the CC BY 4.0 license.

We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:

“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”

Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission.

In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”

2.    If you are unable to obtain permission from the original copyright holder to publish these figures under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only.

The following resources for replacing copyrighted map figures may be helpful:

USGS National Map Viewer (public domain): http://viewer.nationalmap.gov/viewer/

The Gateway to Astronaut Photography of Earth (public domain): http://eol.jsc.nasa.gov/sseop/clickmap/

Maps at the CIA (public domain): https://www.cia.gov/library/publications/the-world-factbook/index.html and https://www.cia.gov/library/publications/cia-maps-publications/index.html

NASA Earth Observatory (public domain): http://earthobservatory.nasa.gov/

Landsat: http://landsat.visibleearth.nasa.gov/

USGS EROS (Earth Resources Observatory and Science (EROS) Center) (public domain): http://eros.usgs.gov/#

Natural Earth (public domain): http://www.naturalearthdata.com/

5. Please include captions for your Supporting Information files at the end of your manuscript, and update any in-text citations to match accordingly. Please see our Supporting Information guidelines for more information: http://journals.plos.org/plosone/s/supporting-information.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: Your well-written paper provides a superb blueprint for addressing a very important need, how best to monitor terrestrial breeding bird populations across the boreal region of Canada. You proposed four scientifically rigorous sampling designs and used simulation modeling to test their efficiencies relative to cost, habitat-representativeness, and spatial balance, all of which are important for implementing a successful and robust monitoring program. You also tested the sampling design empirically in three large study areas, comparing actual vs. projected costs for implementation, a pragmatic consideration for successful adoption and funding of such a large-scale program. Overall, your paper is based on sound analytical procedures, knowledge of how birds are distributed across the landscape, thoughtful reasoning about how populations are likely to change over time, and an understanding of the logistical features needed to implement such an ambitious program in remote areas. Including your R code further illustrates the complexity of the task you have undertaken and will be extremely helpful to those who want to understand the details about how best to design such a program.

There are a few aspects that you could address more fully. First, you use maximization of the variance of diversity of species (with habitat as proxy) rather than minimization of the variance in estimates of population trends or population size as one of the key criteria for allocating effort in terms of primary sampling units (PSUs). I realize that optimization relative to estimating population size or trend would be difficult to achieve, particularly for an omnibus survey that encompasses a multitude of species with varying levels of population change and abundance. I surmise, however, that the proposed allocation would tend to weight southern areas of the boreal more heavily (where species richness is higher) and that there may be low power to estimate trends for populations that are declining rapidly (and have lower population sizes) or those that may be changing rapidly in northern ecotones (perhaps even increasing). It would be interesting if you could address some of these limitations or tradeoffs in the discussion.

Secondly, you mention that the monitoring program will use point count surveys with trained observers and/or autonomous recording units (ARUs). I assume that the decision on which of these methods to use would be related to cost of access for any given PSU. If so, I would conclude that more remote, northern areas could conceivably be monitored exclusively through the use of ARUs alone, because they can be deployed and retrieved least expensively via snow machine or ice roads in winter. It would be useful if the authors could expound at least minimally upon some of the implications of this aspect of the design. Use of ARUs alone without validation from comparable surveys by human observers can result in biased estimates of species’ occurrence and density, particularly for those species that are detected more often by visual than aural cues. ARUs can provide more detailed information on temporal and spatial occurrence of species in a given area if they are set to record samples throughout the season. Analysis of recordings, however, incurs a significant cost, particularly for monitoring of all species, which at this point cannot be automated efficiently. Also, I don’t believe that the cost of the ARUs was factored into the cost of access, but perhaps it should be. In the end, it might be worth a point of discussion to address some of the tradeoffs of using a combination of human observers and ARUs and to suggest that this is one other aspect of efficiency that remains to be designed. I would advocate for a spatially balanced use of human observers and ARUs so that any habitat- or temporal-related biases can be estimated and corrected for. A couple of other recent papers on ARUs in sparse northern habitats you might examine include Thompson et al. (2017, J. Wildl. Manage. 81:1228-1241) and Vold et al. (2017, Wildl. Soc. Bull. 10.1002/wsb.785).

Thirdly, I was a bit confused about how secondary sampling units (SSUs) were selected relative to the presence of water. It seems that any PSUs that encompassed only water were excluded from the sampling frame, which is reasonable. However, it was unclear how the mini-grids of 3 x 3 points were selected with respect to water features within a given PSU. From a cursory look at the R code, it appears that all points that fell within water were also removed from the sampling frame. I’m not certain what this did in terms of influencing which grids might be selected, but you might want to ensure that you are not biasing (low) your sample of points adjacent to water bodies, which are extremely important in terms of species richness and density in the boreal region.

Fourthly, given your proposed, very complicated sampling design, I anticipate that statistical analyses will be extremely complex, which you acknowledge in the discussion. The complexity will depend, of course, on what parameter you are estimating (e.g., habitat relationships, density, changes in distribution, or population trends). Once you add in the temporal dimension, with replication across years, the complexity will be even more challenging, particularly if your sampling intensity or inclusion probabilities change through time. I caution against changing the boundaries of your strata through time for long-term monitoring of population trends, particularly relative to habitat. Not only will vegetation be changing through time, but community composition will likely shift as well, with individual species responding differentially to changes in vegetation, temperature, precipitation, predator abundance, and other ecological factors. Thus, if you retain habitat diversity as a major determinant of inclusion probability, you should consider the ramifications of altering that through time. There are many powerful techniques for modeling trends through time, so this might not be a problem. There could be consequences, however, for the resulting precision of your estimates.

Pertinent to this point is the temporal replication of sampling. It would be helpful to mention what you have in mind for this aspect of your sampling design, particularly for estimating changes in distribution or population size (will there be annual estimates of population trend?). You should also address the possible ramifications of missing data from repeated samples, which will inevitably occur and will prove no small headache, particularly in remote areas where access is so challenging and unpredictable. Typically, repeated samples at the same sites will be more efficient statistically in estimating population trends, but an alternative would be to set up different sampling frames at various intervals through time, with changes in inclusion probabilities. In the boreal zone, sampling individual sites on a biennial basis may be more efficient than sampling them every year because of high interannual correlation (Handel and Sauer 2017). I was glad to see you consider the important issue of how best to incorporate legacy sites. Such data sets can provide key information on long-term changes in distribution and abundance but provide their own challenges when trying to account for potential bias in terms of selection of samples.

Finally, I note a few other minor points that would be helpful to clarify or correct:

P. 11, line 8. What landcover map did you use and how many cover classes were there?

P. 15, line 16. Ralph et al. (1993) is missing from Lit Cited and Matsuoka et al. (2014) should be cited as numbered reference. I also noted that references 71-80 don’t seem to be cited in the text.

P. 19, line 14. Final sentence of the paragraph seems to be missing text after ‘where.’

P. 21, line 17. I would change this to ‘red-filled squares’ so that it is clear what you are referencing. It took me a bit to figure this out (I had to go back to methods to understand what you were saying, especially since red is also used for the triangles representing cost).

P. 22, lines 6-10. I think you meant to reference Fig. 10 instead of Fig. 8. I was confused, however, by the number of PSUs listed. I counted only 2 PSUs in Newfoundland and Labrador (not 16), 2 in the Yukon (not 26), and 12 in Saskatchewan (not 43).

P. 23, lines 9-12. I did not understand the sentence beginning “In addition, the distribution of spatial balance metrics…”

P. 28, line 6. Two l’s in Boreal Avian ‘Modelling.’

Figure 1. You might consider outlining jurisdictions that you sampled in bold lines.

Figure 2. Consider adding a distance scale to each of these submaps or else noting in the figure heading that hexagons are 5 km in diameter and points are spaced 300 m apart.

Figure 4. ‘Proportional’ is misspelled on x-axis label.

I commend you all for this impressive and well-executed effort.

Colleen M. Handel

USGS Alaska Science Center

Reviewer #2: This is a nicely written paper. I would suggest that the authors provide a stronger background on spatial sampling, and sampling optimization:

Delmelle, E. M., & Goovaerts, P. (2009). Second-phase sampling designs for non-stationary spatial variables. Geoderma, 153(1-2), 205-216.

Van Groenigen, J. W., Stein, A., & Zuurbier, R. (1997). Optimization of environmental sampling using interactive GIS. Soil Technology, 10(2), 83-97.

Delmelle, E. (2009). Spatial sampling. The SAGE handbook of spatial analysis, 183, 206.

Reviewer #3: This paper by Van Wilgenburg et al. is well-written, technically sound, and of broad interest to the research community. There are too few papers on the topic of sampling design when it comes to broad scale monitoring programs supported by public funds; and it is critical that we develop strategic methods in areas like the Boreal where many countries have an international responsibility to maintain biodiversity in the face of multiple stressors.

I selected minor revisions because most of my recommended changes can be made fairly easily and no re-analysis or major reworking of the MS are necessary. I made comments directly on the pdf if that helps. More substantive comments can be found from p12 on. My most major concern is that the authors need to carefully consider their use of the term optimal. The optimization literature uses this term in a very specific and mathematical way and it is not clear from this work that the authors performed an optimization. I think a better way of describing what they did was to perform a spatial benefits costs analysis of various broad scale monitoring strategies. They do balance trade offs, but it is not clear that their preferred design is 'optimal' per se. I also think that they might want to call their approach something other than Boreal Optimal Sampling Strategy (BOSS).

The discussion is on the long side. It talks a lot about the advantages of their approach. I think the authors can shorten it. But I also think they should work in a discussion of other literature that examines large scale monitoring programs and trade offs between spatial coverage, costs, and target species representation.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Colleen M Handel

Reviewer #2: No

Reviewer #3: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: PONE-D-19-35579_reviewer_s.pdf

Click here for additional data file.

23 Apr 2020

We would like to thank the reviewers for a very collegial and constructive set of reviews. You will find our responses to your queries below. We hope that you find our responses have adequately answered your questions.

Comments from Reviewer 1:

Reviewer #1: Your well-written paper provides a superb blueprint for addressing a very important need, how best to monitor terrestrial breeding bird populations across the boreal region of Canada. You proposed four scientifically rigorous sampling designs and used simulation modeling to test their efficiencies relative to cost, habitat-representativeness, and spatial balance, all of which are important for implementing a successful and robust monitoring program. You also tested the sampling design empirically in three large study areas, comparing actual vs. projected costs for implementation, a pragmatic consideration for successful adoption and funding of such a large-scale program. Overall, your paper is based on sound analytical procedures, knowledge of how birds are distributed across the landscape, thoughtful reasoning about how populations are likely to change over time, and an understanding of the logistical features needed to implement such an ambitious program in remote areas. Including your R code further illustrates the complexity of the task you have undertaken and will be extremely helpful to those who want to understand the details about how best to design such a program.

There are a few aspects that you could address more fully. First, you use maximization of the variance of diversity of species (with habitat as proxy) rather than minimization of the variance in estimates of population trends or population size as one of the key criteria for allocating effort in terms of primary sampling units (PSUs). I realize that optimization relative to estimating population size or trend would be difficult to achieve, particularly for an omnibus survey that encompasses a multitude of species with varying levels of population change and abundance. I surmise, however, that the proposed allocation would tend to weight southern areas of the boreal more heavily (where species richness is higher) and that there may be low power to estimate trends for populations that are declining rapidly (and have lower population sizes) or those that may be changing rapidly in northern ecotones (perhaps even increasing). It would be interesting if you could address some of these limitations or tradeoffs in the discussion.

Response: Thank you very much, this highlights a need for us to further elaborate and should therefore be a substantial improvement since this is a key component of our approach. You are correct that our stratification will emphasize increased sampling effort in regions with the greatest variance. It is important to note however that optimal sampling theory suggests that increasing sampling in strata with the greatest variance will result in increased precision of the estimate(s) of annual abundance, and therefore trend. We now elaborate on this in the text (1st paragraph of the “Sample size allocation” subsection of the Methods). The addition of a weighting for avian species richness does tend to emphasize sampling in southern regions (Figure 5); however, we see this as advantageous owing to the number of species involved, the access costs, and the spatial overlap with anthropogenic threats. It is true that low densities of northern distributed species may magnify this issue, but an initial sample collected under our design will help to further evaluate the trade-offs involved (e.g. in future simulation modeling). We elaborate on this in the discussion (8th paragraph), and now further discuss how we will re-evaluate this stratification using bird data collected under our design i.e. instead of environmental proxies.

Secondly, you mention that the monitoring program will use point count surveys with trained observers and/or autonomous recording units (ARUs). I assume that the decision on which of these methods to use would be related to cost of access for any given PSU. If so, I would conclude that more remote, northern areas could conceivably be monitored exclusively through the use of ARUs alone, because they can be deployed and retrieved least expensively via snow machine or ice roads in winter. It would be useful if the authors could expound at least minimally upon some of the implications of this aspect of the design. Use of ARUs alone without validation from comparable surveys by human observers can result in biased estimates of species’ occurrence and density, particularly for those species that are detected more often by visual than aural cues. ARUs can provide more detailed information on temporal and spatial occurrence of species in a given area if they are set to record samples throughout the season. Analysis of recordings, however, incurs a significant cost, particularly for monitoring of all species, which at this point cannot be automated efficiently. Also, I don’t believe that the cost of the ARUs was factored into the cost of access, but perhaps it should be. In the end, it might be worth a point of discussion to address some of the tradeoffs of using a combination of human observers and ARUs and to suggest that this is one other aspect of efficiency that remains to be designed. I would advocate for a spatially balanced use of human observers and ARUs so that any habitat- or temporal-related biases can be estimated and corrected for. A couple of other recent papers on ARUs in sparse northern habitats you might examine include Thompson et al. (2017, J. Wildl. Manage. 81:1228-1241) and Vold et al. (2017, Wildl. Soc. Bull. 10.1002/wsb.785).

Response: Here we have chosen to focus on the implications of alternative approaches to stratification and selection of the spatial locations for sampling, but see the points you raise here as key next steps in our program development. Decisions on where and when to allocate ARUs can vary substantially based on access logistics, potential partnerships, etc. that add further complications to the decision. We agree with your suggestion of spatially balanced use ARUs versus human observations would be the optimal approach to avoid spatial, habitat and temporal biases. We briefly (paragraph 8) discuss this in the discussion, but largely point to ongoing and future work that will help us optimize this component of our program.

Thirdly, I was a bit confused about how secondary sampling units (SSUs) were selected relative to the presence of water. It seems that any PSUs that encompassed only water were excluded from the sampling frame, which is reasonable. However, it was unclear how the mini-grids of 3 x 3 points were selected with respect to water features within a given PSU. From a cursory look at the R code, it appears that all points that fell within water were also removed from the sampling frame. I’m not certain what this did in terms of influencing which grids might be selected, but you might want to ensure that you are not biasing (low) your sample of points adjacent to water bodies, which are extremely important in terms of species richness and density in the boreal region.

Response: You are correct that we did indeed eliminate grid points that fell in open water. We also set a minimum target of points to ensure we can meet our within PSU sampling goals. This along with the systematic nature of the grids that are created by the algorithm tends to ensure sampling across a gradient of distances from water (even if the centroid does not fall in or immediately adjacent to water). For example, samples visited under this design thus far in Saskatchewan have represented a range of 0m – 5791 m from large waterbodies, with 5% of the sampling occurring less than 50 m from a waterbody, 12.5% less than 124m, 25% within 218m of water. Below are quantiles of that distribution.

0% 5% 12.5% 25% 50% 75% 87.5% 95% 100%

0m 46m 124m 218m 457m 819m 1231m 1788m 5791m

Fourthly, given your proposed, very complicated sampling design, I anticipate that statistical analyses will be extremely complex, which you acknowledge in the discussion. The complexity will depend, of course, on what parameter you are estimating (e.g., habitat relationships, density, changes in distribution, or population trends). Once you add in the temporal dimension, with replication across years, the complexity will be even more challenging, particularly if your sampling intensity or inclusion probabilities change through time. I caution against changing the boundaries of your strata through time for long-term monitoring of population trends, particularly relative to habitat. Not only will vegetation be changing through time, but community composition will likely shift as well, with individual species responding differentially to changes in vegetation, temperature, precipitation, predator abundance, and other ecological factors. Thus, if you retain habitat diversity as a major determinant of inclusion probability, you should consider the ramifications of altering that through time. There are many powerful techniques for modeling trends through time, so this might not be a problem. There could be consequences, however, for the resulting precision of your estimates.

RESPONSE: We appreciate that the analyses will be more complicated by our choice of design. We have however considered the implications for resulting estimates of trend precision and we anticipate that our choices will result in decreased precision for common species, but should improve precision and accuracy for species associated with rare habitats. As such, we feel that this trade-off will be worthwhile because it is the rare species for which we always have the greatest difficulty estimating trends. We have added brief text to the discussion (paragraph 5) to address this.

Pertinent to this point is the temporal replication of sampling. It would be helpful to mention what you have in mind for this aspect of your sampling design, particularly for estimating changes in distribution or population size (will there be annual estimates of population trend?). You should also address the possible ramifications of missing data from repeated samples, which will inevitably occur and will prove no small headache, particularly in remote areas where access is so challenging and unpredictable. Typically, repeated samples at the same sites will be more efficient statistically in estimating population trends, but an alternative would be to set up different sampling frames at various intervals through time, with changes in inclusion probabilities. In the boreal zone, sampling individual sites on a biennial basis may be more efficient than sampling them every year because of high interannual correlation (Handel and Sauer 2017). I was glad to see you consider the important issue of how best to incorporate legacy sites. Such data sets can provide key information on long-term changes in distribution and abundance but provide their own challenges when trying to account for potential bias in terms of selection of samples.

RESPONSE: We have added brief text to the discussion detailing our longer-term plans. In brief, we are currently focusing on using our design to collect data to fill key gaps in our knowledge on abundance and distribution. These data will subsequently inform simulation analyses to compare and contrast alternative rotating panel designs. We envision revisiting the majority of sites on a five or ten year revisit schedule and have a subset of sites with greater temporal replication (annually) to allow improve precision of trend estimates and estimation of inter-annual variance. Pragmatically, this will involve cost benefit trade-offs and arguments for further resourcing, and thus will not be immediately settled. We attempt to briefly outline future work that will address these issues in the discussion (paragraph 8).

Finally, I note a few other minor points that would be helpful to clarify or correct:

P. 11, line 8. What landcover map did you use and how many cover classes were there?

RESPONSE: We have now added reference to the data source here and discuss the number of landcover classes involved

P. 15, line 16. Ralph et al. (1993) is missing from Lit Cited and Matsuoka et al. (2014) should be cited as numbered reference. I also noted that references 71-80 don’t seem to be cited in the text.

RESPONSE: Added reference to Matsuoka et al and Ralph et al manuscripts and have removed the other references which were included by accident from the reference management software. Thanks for catching this.

P. 19, line 14. Final sentence of the paragraph seems to be missing text after ‘where.’

RESPONSE: Deleted ‘where’.

P. 21, line 17. I would change this to ‘red-filled squares’ so that it is clear what you are referencing. It took me a bit to figure this out (I had to go back to methods to understand what you were saying, especially since red is also used for the triangles representing cost).

RESPONSE: Done.

P. 22, lines 6-10. I think you meant to reference Fig. 10 instead of Fig. 8. I was confused, however, by the number of PSUs listed. I counted only 2 PSUs in Newfoundland and Labrador (not 16), 2 in the Yukon (not 26), and 12 in Saskatchewan (not 43).

RESPONSE: We have corrected the figure reference in the text. We understand the confusion regarding sample sizes. This stems from the figure depicting average access costs to multiple PSUs within the same air charter contract(s) since billing is done on a contract by contract basis and not PSU by PSU. Thus, a single data point may represent access to many PSUs. For example, one of the Yukon data points represents access of 17 PSUs while the data point represents access of 9 PSUs. We have added text to the figure heading to clarify this.

P. 23, lines 9-12. I did not understand the sentence beginning “In addition, the distribution of spatial balance metrics…”

RESPONSE: Reworded to “Across jurisdictions, spatial balance metrics of the spatial design did not always overlap one; suggesting that the inclusion of unequal sampling probabilities did not introduce systematic biases in spatial representation”

P. 28, line 6. Two l’s in Boreal Avian ‘Modelling.’

RESPONSE: Revised

Figure 1. You might consider outlining jurisdictions that you sampled in bold lines.

RESPONSE: We have added the outline as suggested

Figure 2. Consider adding a distance scale to each of these submaps or else noting in the figure heading that hexagons are 5 km in diameter and points are spaced 300 m apart.

RESPONSE: Distance scale added as requested

Figure 4. ‘Proportional’ is misspelled on x-axis label.

RESPONSE: Corrected

Comments from Reviewer #2: This is a nicely written paper. I would suggest that the authors provide a stronger background on spatial sampling, and sampling optimization:

Delmelle, E. M., & Goovaerts, P. (2009). Second-phase sampling designs for non-stationary spatial variables. Geoderma, 153(1-2), 205-216.

Van Groenigen, J. W., Stein, A., & Zuurbier, R. (1997). Optimization of environmental sampling using interactive GIS. Soil Technology, 10(2), 83-97.

Delmelle, E. (2009). Spatial sampling. The SAGE handbook of spatial analysis, 183, 206.

RESPONSE: Thank you for theses helpful references. We have used these and added text providing slightly more reference to and explanation of spatial sampling in the second last paragraph of the introduction and paragraph seven of the discussion.

Comments from Reviewer #3: This paper by Van Wilgenburg et al. is well-written, technically sound, and of broad interest to the research community. There are too few papers on the topic of sampling design when it comes to broad scale monitoring programs supported by public funds; and it is critical that we develop strategic methods in areas like the Boreal where many countries have an international responsibility to maintain biodiversity in the face of multiple stressors.

I selected minor revisions because most of my recommended changes can be made fairly easily and no re-analysis or major reworking of the MS are necessary. I made comments directly on the pdf if that helps. More substantive comments can be found from p12 on. My most major concern is that the authors need to carefully consider their use of the term optimal. The optimization literature uses this term in a very specific and mathematical way and it is not clear from this work that the authors performed an optimization. I think a better way of describing what they did was to perform a spatial benefits costs analysis of various broad scale monitoring strategies. They do balance trade offs, but it is not clear that their preferred design is 'optimal' per se. I also think that they might want to call their approach something other than Boreal Optimal Sampling Strategy (BOSS).

RESPONSE: Thank you for the comments and we can fully appreciate the confusion here. Our use of the term optimal actually derived from the sampling theory literature and not the optimization literature, but the confusion is understandable given our further attempts apply concepts from the optimization literature to select from amongst multiple randomized draws. We now clarify that our use of the term optimal refers to optimal allocation as per the sampling theory literature, and have reworded text around our use of multiple randomized draws to avoid confusion. We further cite the optimization literature related to the weighted sum approach we used in combining the multiple objectives we are attempting to maximize while minimizing costs. In addition, we have added brief text in the discussing Pareto optimality and discussing the pros and cons of our approach to “optimizing” the design.

We also incorporated most of the suggestions made directly within the manuscript. We only provide a detailed response to the more complicated questions/comments or ones where we have chosen not to incorporate below:

In regards to Equation 3, we have modified the equation to note that all pixels (1…j) within the primary sampling unit.

With regard to the note on the 1:1 correspondence line (in the Sample Size Allocation subsection of Results), we have not highlighted accuracy as mentioned, but add further emphasis on values above the 1:1 correspondence line later in the paragraph to highlight increased sampling of more variable strata.

With regard to the comment in the discussion re: “This could be modelled spatially using a more complicated spatial optimization algorithms that factor in reduced costs with neighbouring samples. (e.g., spatial simulated annealing algorithms)”. We have not added text in this regard as we do not yet feel we can sufficiently model these as yet hard to predict logistical efficiencies, though we hope to use our increasing experience on the ground to build improved models in the future. We envision using these refined models to improve our cost models within the same framework to maintain randomized sampling rather than using simulated annealing to seek a global optimum since environmental conditions in the boreal forest change rapidly and thus a global optimum may be a fleeting concept.

With respect to the comment: “Discussion is a bit on the long side and it focuses a lot on advantages of BOSS method. Can it be shortened, but also include reference to other large scale monitoring program designs that consider cost, coverage, and strata? ”. We have endeavored to shorten the text where we could. We note however that many of the requested revisions made by reviewer one have resulted in a slightly longer discussion. We feel that the remaining text is relevant to our results and thus would suggest it remain. We would however welcome suggestions from the Academic editor as to whether there remain any sections that they feel are tangential to the main thrust of the MS. With respect to comparisons against other designs, we are unaware of any that have incorporated all the concepts we have applied, and have kept comparisons to other progams within paragraph 3 of the discussion.

Submitted filename: BOSS Response to the reviewers_Apr 20_2020 .docx

Click here for additional data file.

29 Apr 2020

PONE-D-19-35579R1

A cost efficient spatially balanced hierarchical sampling design for monitoring boreal birds incorporating access costs and habitat stratification

PLOS ONE

Dear Mr Van Wilgenburg,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

I appreciate the time and care you put into your responses to the reviewers. I believe you have addressed their concerns fully, and the manuscript is nearly suitable for publication. In reading through revision 1, I noticed some minor editorial errors (omitted words, punctuation, or the like). Rather than list them here, I've attached a tracked changes version of the manuscript for your reference. Once you address these and resubmit, I will move quickly to accept your manuscript. You do not need to submit point-by-point responses to any of my edits or comments.

We would appreciate receiving your revised manuscript by Jun 13 2020 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A marked-up copy of your manuscript that highlights changes made to the original (R1) version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Frank H. Koch, PhD

Academic Editor

PLOS ONE

Additional Editor Comments (if provided): Please see attachment

[Note: HTML markup is below. Please do not edit.]

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: BOSS_PLoS_ONE_Revision1_fhkoch.docx

Click here for additional data file.

22 May 2020

All changes made by the Academic Editor within the Manuscript Word file have been accepted. We made some additional changes to ensure clarity, accuracy, and adherence to the submission guidelines. We hope everything is now suitably formatted for publication.

Submitted filename: Response to the reviewers.docx

Click here for additional data file.

28 May 2020

A cost efficient spatially balanced hierarchical sampling design for monitoring boreal birds incorporating access costs and habitat stratification

PONE-D-19-35579R2

Dear Dr. Van Wilgenburg,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Frank H. Koch, PhD

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Thank you for addressing the latest set of comments. The manuscript reads well and should interest many in the PLOS ONE audience.

Reviewers' comments:

1 Jun 2020

PONE-D-19-35579R2

A cost efficient spatially balanced hierarchical sampling design for monitoring boreal birds incorporating access costs and habitat stratification

Dear Dr. Van Wilgenburg:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Frank H. Koch

Academic Editor

PLOS ONE