Modeling the occurrence of the Yellow Rail (Coturnicops noveboracensis) in the context of ongoing resource development in the oil sands region of Alberta

Yellow Rails (Coturnicops noveboracensis) are among the most secretive bird species in North America. They are poorly sampled by common survey protocols, and as a result their occurrence across much of their range is uncertain. We compiled occurrence records of the species and used resource selection functions to classify habitats as selected, neutral, or avoided using four different land cover maps in the oil sands region of northeastern Alberta. We assessed the accuracy of these maps using 279 previously unsurveyed locations and showed that a consensus-based ensemble classifier predicted occurrence more accurately than any single map. We combined the four maps into one map that rated habitat on a scale from 0 (consensus avoided) to 8 (consensus selected). Occupancy analysis showed increasing occupancy rates in areas with higher habitat suitability classes, with maximum occupancy rates of 0.18 (95% CI: 0.07-0.32) in class 8 habitat. We combined detections of 169 male Yellow Rails at surveyed locations with model predictions for unsurveyed locations to produce two population estimates for our study area, based on two estimates of the detection radius of the species. The estimate assuming a 150-m detection radius was 2747 males (95% CI: 588-5563), and the estimate assuming a 250-m detection radius was 1650 males (95% CI: 416-3266). Although estimates contained substantial uncertainty, our results suggest a larger number of Yellow Rails in the region than previously thought, which alters the current understanding of the distribution of this species. We estimated that about 17% of the population in our study area resides on oil sands leases that cover 14% of the study area, in habitats facing ongoing and future industrial development. The availability of a habitat map based on empirical evidence and detailed analyses for this species of conservation concern will improve targeted monitoring and promote mitigation of potential effects of development. Modélisation de l'occurrence du râle jaune (Coturnicops noveboracensis) dans le contexte du développement continu des ressources dans la région des sables bitumineux de l'Alberta RÉSUMÉ. Le râle jaune (Coturnicops noveboracensis) est l'une des espèces d'oiseau les plus secrètes d'Amérique du Nord. Les protocoles d'enquête courants ne permettant pas de les échantillonner correctement, leur occurrence sur une grande partie de leur territoire demeure incertaine. Nous avons compilé les enregistrements d'occurrence de cette espèce et utilisé des fonctions de sélection des ressources pour classer les habitats comme suit : sélectionnés, neutres ou évités. Nous avons pour ce faire utilisé quatre cartes différentes de couverture du territoire dans la région des sables bitumineux du nord-est de l'Alberta. Nous avons évalué l'efficacité de ces cartes en utilisant 279 sites qui n'avaient encore jamais été étudiés et démontré qu'un classificateur global reposant sur un consensus permet de prévoir l'occurrence avec plus de précision que toute carte individuelle. Nous avons combiné les quatre cartes en une seule, qui classait les habitats sur une échelle de 0 (globalement évités) à 8 (globalement sélectionnés). L'analyse d'occupation a démontré des taux d'occupation croissants dans les zones présentant des classes supérieures d'adéquation à l'habitat de 0,18 (IC de 95 % : 0,07 0,32) dans l'habitat de classe 8. Nous avons combiné les détections de 169 râles jaunes mâles sur les sites observés avec les prévisions modélisées des sites non étudiés afin de produire deux estimations de population pour notre zone d'étude, reposant sur deux estimations du rayon de détection de l'espèce. L'estimation supposant un rayon de détection de 150 m était de 2 747 mâles (IC de 95 % : 588-5 563), et l'estimation supposant un rayon de détection de 250 m était de 1 650 mâles (IC de 95 % : 416-3 266). Bien que ces estimations présentent une incertitude substantielle, nos résultats suggèrent que les râles jaunes sont plus nombreux dans la région qu'on ne le croyait précédemment, ce qui modifie la compréhension actuelle de la répartition de cette espèce. Nous avons estimé qu'environ 17 % de la population dans notre zone d'étude réside dans des territoires riches en sables bitumineux qui couvrent 14 % de la zone d'étude. Ces habitats sont confrontés à un développement industriel actuel et futur. La disponibilité d'une carte d'habitat basée sur des preuves empiriques et des analyses détaillées de cette espèce dont la conservation est préoccupante améliorera la surveillance ciblée et favorisera l'atténuation des effets potentiels de ce développement.


INTRODUCTION
Assessing the status of most wildlife populations requires dedicated survey effort across large geographic areas. For many North American bird species, standardized survey programs such as the Breeding Bird Survey (BBS) provide information used to estimate population sizes and trends (Sauer et al. 2017). However, birds breeding in the boreal forest region are poorly sampled by the BBS, owing to the scarcity of roads in the region and lack of volunteers to conduct surveys (Cumming et al. 2010). Wetlandassociated and nocturnal species are particularly poorly sampled (Conway 2011), because surveys typically take place during daylight hours, and wetlands are under-represented in road-based survey routes.
Yellow Rails (Coturnicops noveboracensis) reside in wetlands, are nocturnal, and have a breeding distribution that includes much of the North American boreal forest (Leston and Bookhout 2015). These traits make them one of the most secretive and least known bird species in North America. The lack of information on the distribution and population size of Yellow Rails poses challenges for conservation. For example, losses of wetland habitats have occurred across the southern portion of their breeding range and in wintering areas (Eddleman et al. 1988), but in only a few cases have predisturbance observations been available to link habitat loss with local population declines (Alvo and Robert 1999). Concerns about ongoing habitat loss, coupled with the high likelihood that declines might go undetected, motivated their listing as a species of special concern by the Committee on the Status of Endangered Wildlife in Canada (COSEWIC) in 1999 and as a Schedule 1 species under the Species at Risk Act (S.C. 2002, c.29;SARA) in 2005 (COSEWIC 2009). A recent management plan for the species called for researchers to "conduct regular surveys throughout and beyond the species' range within suitable habitat to find additional key sites" (Environment Canada 2013:16). A first step toward that goal is to identify suitable habitat for surveys. Minimal information is available regarding the location of Yellow Rail habitat for most parts of their range. Moreover, in many areas, only sporadic records of the species exist (Smith 1996, Federation of Alberta Naturalists 2007 and few systematic survey efforts have been carried out (Prescott et al. 2002, Robert et al. 2004). Thus, whether Yellow Rails occur at all remains unknown for many regions.
An efficient way to fill survey gaps and facilitate the discovery of breeding sites is through the use of land cover maps, which can be combined with empirical observations to model habitat suitability across large geographic areas. A common way of analyzing habitat suitability using land cover maps is through resource selection functions (RSFs). In an RSF, the occurrence of animals is compared to the availability of different vegetation types or land cover categories (Boyce and McDonald 1999). Each land cover category can then be classified as either selected, avoided, or neutral. Resulting maps can be used in risk assessment, impact mitigation, and conservation planning. This is straightforward when applied to a single land cover map, but when several land cover maps are available for a region, researchers must decide which map to use or whether to combine their predictions. Evaluating these options is important because these choices will influence our understanding of habitat suitability. Differences between land cover maps can arise from different definitions of land cover categories, different numbers of land cover categories, different classification rules, and classification errors. In a review of research on global forest cover, for example, Sexton et al. (2016) noted that the definition of forest can range from > 10% tree cover to > 30% tree cover, with implications for land cover mapping and global forest area calculations. In a study of Philippine forest cover, eight land cover maps had just 30% consensus agreement on the classification of forest vs. nonforest, with differences attributed to a combination of different definitions of land cover categories and classification errors (Estoque et al. 2018). Similar problems exist with respect to the identification of wetlands. In a comparison of four land cover maps in Siberia with ground-truth observations, the maps classified wetland vs. upland with between 2% and 56% accuracy (Frey and Smith 2007). These findings highlight the potential for large discrepancies among maps.
Our study area covers much of the oil sands region of northeastern Alberta. At least four land cover maps are available for our study area based on four different sensor technologies including MODIS, Landsat, Radarsat, and aerial imagery. These maps differ in their spatial resolution as well as the number and identity of land cover categories mapped. Each map also contains classification errors. We sought to explore the consequences of these differences for predicting the occurrence of Yellow Rail.
Field-based descriptions of Yellow Rail breeding habitat show that the species uses wet areas with water depths typically < 30 cm and vegetation cover dominated by graminoid vegetation such as sedges and rushes (Bookhout and Stenzel 1987, Martin 2012, Austin and Buhl 2013. This relatively narrow habitat niche may give the impression that modeling and mapping this species should be simple, because occurrence may be concentrated in one or a few land cover classes. However, of the four land cover maps of our study area, only Ducks Unlimited Canada's Enhanced Wetland Classification map includes a category specific to graminoid wetlands (Ducks Unlimited Canada 2011). In another, the Alberta Biodiversity Monitoring Institute's 2010 land cover classification, low classification accuracy of wetland classes during the creation of the map led wetlands and open upland habitats to be collapsed into categories labeled "shrubland" and "grassland," to increase the overall accuracy of the map (Castilla et al. 2014). Collapsing wetland categories in this way might reduce the utility of a map for predicting wetland bird occurrence because Boreal wetland classes harbor distinct bird assemblages (Morissette et al. 2013). On the other hand, finer subdivisions between wetland classes are likely useful only if they are accurately classified; graminoid wetlands may be difficult to classify accurately from imagery because they share similar spectral characteristics with other wetland and upland land cover classes (Ozesmi and Bauer 2002). In light of these considerations, it is not obvious a priori which map or maps might perform best for modeling this species. In addition, it is unknown whether it is better to use "the best fitting map" or to combine them to look for areas where the different products agree or disagree.
Habitats in the oil sands region face threats from ongoing resource development. A recent analysis found that the proportion of land area that has been directly disturbed by residential, recreational, or industrial developments nearly doubled from 4.8% to 8.4% from 1999 to 2015 (Alberta Biodiversity Monitoring Institute 2017). These changes have centered around the development of oil leases. Two types of leases exist, with different disturbance profiles: in situ leases involve extracting oil from deep underground using steam injection. Such activities produce a large number of roads, linear features, well pads, and other infrastructure but most of the landscape remains vegetated (Alberta Chamber of Resources 2004). In contrast, mining leases employ methods that remove most vegetation and can dewater wetlands for the operational life of the mine with reclamation occurring after mining is completed (Alberta Chamber of Resources 2004).
The importance of this region for Yellow Rails has historically been poorly known, with records of the species including approximately 10 records in our study area during the two Alberta Breeding Bird Atlases (Semenchuk 1992 Resources Ltd. 2011). Five additional historical records from our study area were also described by Prescott et al. (2002). Taken together, evidence prior to this study pointed to Yellow Rails using the region, but the degree of use and the abundance of the species has been unclear. The potential for resource extraction to affect wetland habitats in this region has amplified the need to map and survey critical habitat.
In this study, we present the results of bioacoustic surveys for Yellow Rails in the oil sands regions of northeastern Alberta, Canada. Our goals were threefold. First, we sought to produce habitat selection maps for the species using four existing land cover maps and compare the performance of the four individual maps to the performance of an ensemble approach for predicting Yellow Rail occurrence. Second, we aimed to provide a quantitative assessment of the population numbers of Yellow Rails in our study area. Third, we examined the potential risk that proposed oil sands developments represent for the species.

Study area
The study took place in the eastern portion of the Athabasca Oil Sands Area in northeastern Alberta, Canada (Fig. 1). The study area is in the Boreal Plains ecozone and is 48,715 km² in size, with an irregular shape owing to the patchy availability of remote sensing products. The study area covers much of Alberta's oil sands region, so many of its habitats face pressing conservation challenges. Habitats in the region include deciduous, evergreen, and mixed upland forests, and several types of wetlands. A recent examination of the boreal region of Alberta estimated that wetlands cover 45% of the region (DeLancey et al. 2019). Wetlands broadly fall into five main classes, bog, fen, marsh, swamp, and shallow open water, with finer classifications possible within each of these (Ducks Unlimited Canada 2018). Of these classes, some types of fens and marshes might be suitable for Yellow Rails because of their dense cover of sedge (Bookhout and Stenzel 1987). Swamps and bogs are expected to be less suitable because of their high cover of shrubby vegetation and trees (Martin 2012, Austin andBuhl 2013), and, similarly, open water is not expected to be suitable because Yellow Rails require dense herbaceous vegetation cover.

Data collection
We used three sets of Yellow Rail survey data: (1) to summarize survey effort and known detections of Yellow Rails; (2) to produce habitat selection maps; and (3) to assess the predictive accuracy of the selection maps and produce population estimates. To produce the habitat selection maps, we compiled a set of existing Yellow Rail detections as inputs (the "historical detections dataset," detection-only data). In addition to the historical data, we also compiled another dataset to summarize more recent detection/nondetection effort (the "recent detection/nondetection dataset," detection/nondetection data). To assess the performance of the selection map, we conducted additional acoustic surveys in 2018 and 2019 to produce an independent dataset for testing (the "test dataset," detection/nondetection data). The test dataset is distinct from the other two datasets, with surveys being conducted at previously unsurveyed locations. The historical detections dataset includes both incidental human-based detections and autonomous recording unit (ARU) acoustic detections, while the latter two datasets were exclusively derived from ARU data. Additional details on each dataset are provided below. Summaries of the three data sets, and the analyses each contributed to, are provided in Table 1.

Historical detections dataset
We used several sources of data to compile detections of Yellow Rails for the purposes of creating habitat selection maps. We sourced observations from three online databases:

Recent detection/nondetection dataset
We retrieved additional recent detection/nondetection data from the Bioacoustic Unit database in 2018. Our intent was twofold: to extract additional Yellow Rail detections from data processed after the creation of the habitat selection maps, and to establish which locations had been surveyed without detecting a Yellow Rail (nondetections). These data comprise bioacoustic surveys analogous to point counts in which a trained surveyor listened to recordings and visually scanned spectrograms for vocalizing animals (H. E. Lankau, A. MacPhail, M. Knaggs et al., unpublished data). The surveyor identified detected species, estimated abundance, and noted the time of first detection within the recording for each individual. In this dataset, abundance estimation was carried out via the listener's impression of the number of rails in the recording. In some cases, the listener noted there were "too many to count" and provided a range of possible values (e.g., three to five), in which case we used the lower bound as a conservative estimate. If there were too many individuals to count, but no subjective estimate was available, we coded abundance as two. Automated detectors were not used because of the lack of reliable and accurate software options. Sound recordings were made using Songmeter SM2, SM3, or SM4 recording units (Wildlife Acoustics, Inc.).
We considered a location to be surveyed, and therefore worthy of inclusion in this dataset, if at least four one-minute acoustic surveys had been conducted within a single year at times of day when Yellow Rails are known to be most vocal (11:00 pm to 3:00 am), and during the period of the year when they are known to be present and vocal in this region (1 June to 15 July). We have detected Yellow Rails in this region as early as 15 May and as late as 29 July (E. Bayne, unpublished data), suggesting that our chosen date range is appropriate. Our experience in this regard is supported by the findings of Martin et al. (2014), who found no evidence for temporal changes in detection probability of Yellow Rails on surveys conducted from 23 May to 5 July in Manitoba. When multiple years of data were available from a location, we discarded all but the most recent year. We did this in order to avoid bias on our counts, and to ensure that our estimates incorporated the most up-to-date information from each station. There was some overlap between this dataset and the historical detections dataset, specifically at locations in the historical detections dataset that had sufficient survey effort to meet the above inclusion criteria for this dataset. In some instances, however, locations included in the historical detections dataset were converted from detections to nondetections if the location was surveyed in multiple years and a rail was not detected in the most recent year. The recent detection/nondetection dataset included surveys from 2012 to 2017, with a total of 2474 surveys at 564 locations within our study area.

Test dataset
To test the predictive accuracy of the habitat selection maps, one observer (RWH) conducted acoustic surveys at a total of 279 previously unsurveyed locations from 22 May to 15 July of 2018 (N = 236 locations) and 2019 (N = 43 locations). A completely held-out test dataset was deemed desirable because of concerns related to overfitting when a model (in this case habitat selection maps) is tested using data that was also used as input data (James et al. 2013). In this study, the overlap between the historical detections dataset and the recent detection/nondetections led us to use an entirely distinct set of data for testing.
Ideally, a test dataset would include a random sample of locations covering each of the habitat classes in each habitat selection map. Because of the remoteness of much of our study area, producing a truly random test set was not possible. Instead, we used acoustic recordings collected in the course of various Bioacoustic Unit projects to estimate the occupancy rate of each habitat class in each map. We chose some locations (N = 17) specifically to assess the occupancy rate of habitat deemed selected by all four habitat selection maps. We targeted these consensus-selected habitats because they were rare, covering less than 1% of the study area, but likely important for rails. The remainder (N = 262) were selected for other projects without consideration of any habitat map: some survey locations were placed systematically (e.g., in a grid formation spaced by 600 m, N = 92 locations), and some were placed with other specific research questions in mind (N = 170 locations, e.g., placed in forest to assess the effect of forest regeneration on the bird community).
Recording units deployed at these locations recorded three minutes of audio at the start of each hour during the middle of the night. Units were deployed during the day, so human presence would not influence bird vocal activity. We extracted the first minute of recordings to conduct one-minute surveys at each location on four randomly selected nights between 1 June and 15 July for 267 locations. The other 12 locations did not have enough recordings in this date range; for those locations we used recordings as early as 22 May. We conducted surveys at 2:00 am, except where recordings from that time were not available (N = 14 locations); in such cases we analyzed recordings from other times between 11:00 pm and 3:00 am. Surveys involved visually scanning spectrograms in Audacity (Audacity Team, https://www. audacityteam.org/) to detect Yellow Rails, with listening used to confirm ID. Visual scans have been shown to be an accurate and fast way of processing data for this species (Sidie-Slettedahl et al. 2015). The surveyor estimated the number of calling Yellow Rails using the methods described by Drake et al. (2016), which involved counting the number of ticks within windows of approximately 170 milliseconds. Counts made via this method are highly accurate up to approximately six individuals; the highest count on any survey in this study was four. We used one-minute surveys because we have found that Yellow Rails have a high detection probability at 2:00 am (E. Bayne, unpublished data). Most often, when a rail was present, it was detected within the first few seconds of a recording, making surveys longer than a minute unlikely to provide additional information. This process produced data from 1116 surveys at 279 locations.

Creation of habitat selection maps
The four map products used to produce habitat selection maps were (1) Alberta Biodiversity Monitoring Institute (ABMI) 2010 Land cover classification (Castilla et al. 2014). This is an Albertawide land cover map based on digital classification of 30-metre resolution Landsat satellite images and enhanced using ancillary datasets. It consists of 15 classes, including water, shrubland, grassland, agriculture, exposed land, developed land, and several forest types; (2) Ducks Unlimited (DU) Enhanced Wetland Classification (Ducks Unlimited Canada 2011). This is a map based on a combination of Landsat and Radarsat, with some additional data collected via helicopter surveys. Using objectbased supervised classification methods, this product classifies 30-metre resolution cells into one of 19 boreal plains wetland classes and 10 other land cover classes. Images used to build this map were acquired after 2010; (3) Land cover classification of Canada (LCC; Latifovic et al. 2008). Based on MODIS satellite data from 2005, the product has 250 m spatial resolution with 39 classes. Reported accuracy of land cover classification in this map was about 70% (Latifovic et al. 2008). We converted it to a 30 m by 30 m raster; (4) Alberta Vegetation Inventory (AVI; Alberta Sustainable Resource Development 2005). The Alberta Vegetation Inventory is an aerial photo-based inventory developed to identify the type, extent, and conditions of forest vegetation. The aerial imagery is updated on a continuous basis, and in our case may have included photographs from as early as the 1980s to 2013. Photo interpreters digitized polygons and assigned attributes that were then used to classify them into 37 land cover classes including nine wetland types. The accuracy of land cover classification is not known for this map. Additional details on the creation of this land cover map can be found in Mahon et al. (2016). We converted land cover polygons on this map to a 30 m by 30 m raster.
We analyzed the data in the historical detections dataset with resource selection functions based on the logistic discriminant method to determine whether land cover classes within the four mapping products were avoided, neutral, or selected (Manly et al. 2002). Avoided land cover classes were those with a negative beta coefficient and a 95% confidence interval that did not include zero. Selected land cover classes ("selected habitat") were those with a positive beta coefficient and a 95% confidence interval that did not include zero. All other land cover classes were classified as neutral. Results are shown as habitat selection ratios where a ratio > 1 indicates selection and < 1 indicates avoidance.

Assessment of habitat map accuracy
In assessing predictive performance, we sought to answer two main questions: first, how accurate are the four individual maps at predicting Yellow Rail occurrence? Second, do the four maps combined outperform the individual maps? To assess accuracy, we determined whether each location in the test dataset was within selected habitat for each of the four individual maps. Stations in neutral or avoided habitat were considered "not selected habitat," producing a binary classification scheme. To combine the four maps, we used a simple consensus-based ensemble approach: a location that was deemed selected by all four maps was coded as selected, and all other locations coded as not selected. Next, we calculated the number of true positives, false positives, true negatives, and false negatives for each map and for the combined map. A true positive was a station in selected habitat where a Yellow Rail was detected. A false positive was a station in selected habitat where a Yellow Rail was not detected. A true negative was a station in nonselected habitat where no Yellow Rails were detected. A false negative was a station in nonselected habitat where a Yellow Rail was detected. The accuracy of each map was calculated as the proportion of the 279 stations in the test dataset that were correctly classified, i.e., accuracy was true positives plus true negatives divided by 279. To assess uncertainty in these values, we used bootstrapping to resample the test data set 10,000 times.

Production of a composite map
Because the ensemble approach was more accurate than any of the individual maps, we combined the four maps to produce a habitat suitability map for the entire study area. We used this composite map for all subsequent analyses, and have made it publicly available (see results). To produce the composite map, we coded avoided, neutral, and selected habitat as 0, 1, and 2, respectively, then took the sum of the four RSF maps for each 30 m x 30 m cell across the study area. Thus, the composite map rated each cell on a scale from 0 to 8. Cells with a value of 8 corresponded to locations where the four individual maps agreed were selected by Yellow Rails, while cells with a value of 0 corresponded to locations where the four individual maps agreed were avoided by Yellow Rails. Mapped values between 1 and 7 included locations where the individual maps agreed were neutral (class 4), or where at least one of the four maps disagreed with the others in their assessment of whether the location was selected.
Most intermediate values could arise through multiple possible combinations of selected or avoided. For example, class 6 habitat could result from areas selected under three maps and avoided in the fourth, or selected in two and neutral in the other two. In general, however, we expected cells with higher values to be more likely to be occupied by Yellow Rails.

Occupancy analysis for population estimation
We used occupancy analysis to assess the rate of occupancy in each of the nine habitat suitability classes (0-8) of the composite map for the purposes of estimating the number of Yellow Rail males in our study area. We restricted our analysis to males because males are the only sex known to produce the "tic-tic, tictic-tic" vocalization that we detected on surveys (Stalheim 1975).
Occupancy analysis provides a simultaneous estimate of occupancy rate (the proportion of locations with a Yellow Rail) and detection probability (the probability of detecting a rail on each survey at an occupied location), while correcting for the possibility of missing the target species at an occupied location (Mackenzie et al. 2002).
We calculated occupancy rates using the 279 locations surveyed in the test dataset. Locations spanned the nine habitat suitability classes, with a range from 17 to 45 locations per class. We included the habitat suitability class from the composite map as the lone categorical occupancy predictor. To account for possible differences in vocal activity throughout the breeding season, we used AIC to compare three models that differed in having no detection covariates, a linear function of survey date as a predictor, or a quadratic function of survey date as a predictor.
The model with the lowest AIC was considered the best model, and models with AIC scores within two of each other were considered similar. We also estimated the expected occupancy rate of a randomly selected point by calculating a weighted average occupancy of the nine habitat classes in the composite map, where the predicted occupancy rates were weighted by the proportion of the study area in each habitat class.
We conducted occupancy analysis in the package "unmarked" (Fiske and Chandler 2011) in R (version 3.4.2, R Core Team 2019). We used bootstrapping to estimate the precision of parameter estimates. To do so, we drew 10,000 bootstrap samples of the 279 locations in the test dataset, sampled with replacement (James et al. 2013). To avoid instances when bootstrapping left few or no data points in a particular habitat class, leading to artificially high uncertainty, we conducted bootstrapping in a stratified way such that the within-class sample size remained constant across bootstrap samples.

Population estimation
We used the occupancy model to extrapolate occupancy rates to the broader study area, by creating a hexagonal grid across the study area and applying predicted occupancy rates and abundance estimates to each hexagon centroid. We used a hexagon grid because, of the three types of regular tessellations of a plane (triangles, squares, and hexagons), the hexagonal grid is the nearest in shape to a circle, resulting in the most even coverage (Birch et al. 2007). The hexagon simulation methodology is depicted in Figure 2.
For a hexagonal grid, the only parameter to be adjusted is the area of each hexagon or, equivalently, the distance between adjacent hexagon centroids. With a small spacing parameter,

Fig. 2.
A representation of the hexagonal grid methodology used to produce a population estimate for the study area. A grid was simulated across the study area (dotted lines) and points generated at the centroids (gray dots). Two grid sizes were used (only one is shown here), arrived at via simulation under different assumptions regarding the detection radius of autonomous recording units (ARUs), Yellow Rail (Coturnicops noveboracensis) territory size, and territory shape (see text). The habitat class at each point was extracted (background colors represent habitat classes from zero to eight; darker colors of blue represent higher estimated habitat suitability). This allowed estimation of occupancy rates across the study area based on the occupancy model. Hexagon tiles were excluded if they contained a Yellow Rail detection (black square) or at least four nocturnal surveys without a detection (white square). These previously surveyed locations composed the confirmed population of Yellow Rails, while simulated points were used to produce an estimated population; added together, the confirmed and estimated populations composed the final population estimate. more individuals will be double-counted and few will be missed. In contrast, if the spacing parameter is large, many individuals will be missed and few or none will be double-counted. If doublecounting or missed birds predominate, occupancy rate will overor underestimate the true population size, respectively (Efford and Dawson 2012). We propose, however, that occupancy will be closely correlated with population size with minimal bias if the grid size is chosen such that the numbers of double-counted birds and missed birds are equal. We determined the optimal grid size via simulation, by simulating 10,000 Yellow Rail territories in a hypothetical 20 km x 20 km area. An individual was counted if its territory overlapped the survey radius of a hexagon centroid, and missed if it did not.
Assumptions were required regarding three separate parameters in these simulations: survey radius (plot size), Yellow Rail home range size, and home range shape. We set our survey radius to be the effective detection radius of Yellow Rails recorded on ARUs. Uncertainty exists in this parameter, having been estimated using call broadcast experiments as low as 150 m (Drake et al. 2016) and as high as 250 m (Yip et al. 2017). Accordingly, we ran two separate simulations (one assuming a 150 m detection radius and one assuming a 250 m detection radius) to produce two hexagonal grids of different tile size and two population estimates. For Yellow Rail home range size, we randomly sampled each home range from a range of published values provided by Bookhout and Stenzel (1987), who found that males had territories from 5.8 to 10.5 ha. We randomly drew home range sizes from a uniform distribution between these two values. No information was available regarding home range shape, so we assumed territories were elliptical with a major-to-minor axis ratio varying randomly from one, i.e., circular, to three. Home ranges were placed at random within the simulation area, then we simulated hexagonal grids, calculated centroids, and buffered those centroids by the appropriate detection radius. We determined whether each territory overlapped (occupied) exactly one survey location, two survey locations (double-counted) or no survey locations (missed). The optimal grid size was determined to be the grid size at which the number of double-counted individuals was equal to the number that were missed in the simulation.
A source of bias not considered in the above process is the tendency for territories to show clumped distributions, which can also decouple occupancy rates from population size (Efford and Dawson 2012). We accounted for this by estimating the average cluster size (Buckland et al. 2001), defined as the mean high count of individuals at survey locations with at least one detection. We multiplied this number by the occupancy rate to correct for clustered Yellow Rail territories. We expect that this method is conservative for two reasons: first because human listeners are trained to be conservative when counting individuals calling at any one time, and second because all individuals present at a site do not always vocalize simultaneously, so high counts during short surveys tend to underestimate the true number present at a location. As described above, when surveyors in the recent detection/nondetection dataset found that there were "too many to count," we coded abundance conservatively by taking the lower bound of their subjective estimate when possible, or coded abundance as two if no subjective estimate was made.
To produce the final population estimate, we generated a hexagonal grid across the entire study area with the two tile sizes arrived at in the above simulations. The grid was generated using the generate tessellation tool in ArcMap 10.6.1 (ESRI 2017). We classified each hexagon as either surveyed or unsurveyed. Surveyed hexagons were those in which we had conducted at least four nighttime surveys in a year, in either the recent detection/ nondetection dataset or the test dataset. We did not include the historical detections dataset in this analysis because of concerns that detections dating back as far as 1998 may not provide an accurate assessment of the current population. Yellow Rails detected in surveyed hexagons were considered a confirmed population (N Confirmed ), and surveyed hexagons were excluded from subsequent population estimation. For each unsurveyed hexagon, we retrieved the habitat class at the centroid. The population size for the study area, E(N), was estimated by adding the confirmed population to an estimated population with the following formula: Where Ψ i was the occupancy rate for habitat class i; X i was the number of unsurveyed simulated points that fell within habitat class i; and E(s̅ ) was the cluster size, or mean number of detected individuals at occupied locations. We used the 10,000 bootstrapresamples of the occupancy rates to assess the uncertainty in the estimate of E(N).

Assessing threats of oil sands development to Yellow Rails
Concurrent with the above population estimate, we assessed whether each survey location fell within the boundaries of an oil sands lease. For the surveyed portions of the study area, this produced a simple count of the number of Yellow Rail individuals detected on oil leases. For the unsurveyed areas, hexagon centroids that fell within an oil lease were used to estimate the number of individuals that may be found on lease, as above. We produced separate estimates for in situ oil leases and surface mining leases.

Survey results
Considering only the recent detection/nondetection dataset and the test dataset, a total of 169 male Yellow Rails were detected. Of these, 146 were from the recent detection/nondetection dataset (ARU-based detections with ≥4 nocturnal surveys per location), and 23 from the test dataset (ARU-based surveys in 2018 and 2019 with 4 nocturnal surveys per location).

Creation of habitat selection maps
The four original land cover maps differed substantially in which habitat types, and how many habitats, were deemed selected (Fig.  3). This was expected, because the original maps differed with respect to the number and types of habitat categories. However, in some cases, habitat types with the same or similar names were assigned to different selection categories. For example, water was classified as selected by the ABMI and LCC maps and avoided by the DU and AVI maps. Similarly, recently burned areas were classified as selected by the AVI map and avoided by the other three maps. These disparate results highlight discrepancies between maps, which could propagate in subsequent analyses.

Comparison of habitat selection map accuracy
Surveys at 279 locations in the test dataset found a total of 23 Yellow Rails at 12 survey locations. The accuracy of the four individual habitat selection maps at predicting the occupancy of these 279 stations varied from 0.58 to 0.77. In contrast, the accuracy of the ensemble classifier that classified by consensus was 0.86 (

Occupancy analysis
In our occupancy analysis based on the test dataset and the composite map, the three potential models, with no detection probability covariates, a linear date covariate, and a quadratic date covariate, respectively, were similar when compared via AIC (Δ = 0.03, 0 and 0.59, respectively). This suggested minimal variation in detection probability across the range of survey dates.
In addition, the model including habitat classes had lower AIC than the null model without them (ΔAIC = 9.5). We therefore selected the model with no detection covariates and with habitat class as the occupancy covariate for subsequent analysis. Using this model, detection probability was estimated to be 0.63 per 1minute survey (95% CI: 0.41 to 0.81). With this detection probability, the estimated false negative rate across four oneminute surveys at an occupied location is 1.9% (95% CI: 0.1% to 12%).
We detected 16 Yellow Rails at 8 of 45 survey locations in class 8 habitat in the test dataset. Occupancy analysis showed that locations in this class had occupancy rates of 0.18 (Fig. 4a, 95% CI: 0.070 to 0.315). We found four Yellow Rails at two locations in class 7 habitat, which had the second highest occupancy rate (0.069, 95% CI: 0 to 0.174). This implies that this habitat class was also predictive of breeding habitat, albeit less so than class 8 habitat. We detected no individuals in class 6 habitat, and two individuals at one location in class 5 habitat. The lack of detections in class 6 habitat may have been due to relatively low survey effort in that habitat class (N = 20 locations). Our results generally support the notion that higher habitat suitability classes in the composite map tend to be occupied at higher rates. As a point of comparison, the weighted average occupancy rate across the study area was estimated to be 0.010 for a random location (95% CI: 0.0018 to 0.021). We estimated, therefore, that by restricting survey effort to the most suitable habitat class, Yellow Rails can be encountered roughly 18 times more efficiently than if searches were conducted at random (0.18 vs 0.010 occupancy rates).

Population estimate
The mean number of detected males per occupied survey location in our data was 1.69 (N = 100 occupied locations in the recent detection/nondetection dataset and the test dataset), which was used as an estimate of cluster size, E(s̅ ), in our population estimates. Assuming a 250 m effective detection radius, our simulations suggested that the optimal distance between adjacent hexagon centroids was 802 m. This spacing produced 87,982 hexagon centroids in our study area, after removing 605 hexagons that had already been surveyed at least once. Using this hexagonal grid, the population of our study area was estimated to be 1650 males (95% CI: 416 to 3266 males; Fig. 4b, Table 3).
The model assuming a 150 m effective detection radius led to a higher estimated population size, as expected. The optimal distance between adjacent hexagon centroids was 605 m in this model, which produced 153,681 hexagon centroids in our study area, not including 718 hexagons that had already been surveyed. The population of our study area was estimated from this model to be 2747 males (95% CI: 588 to 5563 males; Fig. 4c, Table 3).

Threats posed by oil sands developments
Summaries of the estimated numbers of male Yellow Rails on and off lease, and by lease type, are provided in Table 3. For off lease areas and in situ areas, there was a close concordance between the percent of estimated Yellow Rails in a lease type category and the percent of the total study area. Off lease portions of the study area covered 86.2% of the total area and contained roughly 84% of the estimated Yellow Rails; in situ leases covered 9.2% of the study area and contained about 9% of the estimated rails. In contrast, mining leases covered 4.7% of the area, and contained about 7% of the estimated rails. Hence, mining leases were estimated to contain a disproportionate number of individuals relative to their area.

Production of a Yellow Rail habitat suitability map
We mapped the selected habitat of Yellow Rails in the Lower Athabasca region of Alberta, Canada using four available land cover maps, Yellow Rail surveys, and resource selection functions. We showed that an ensemble approach based on the consensus of the four maps produced more accurate predictions than any of the individual maps. Relying on consensus locations reduced false positives and increased true negatives relative to predictions of any of the four individual maps.
Substantial differences existed between the four individual maps in terms of which habitats were deemed selected or avoided. We expect that when researchers rely on a single map, problems may arise during the classification of habitats into discrete selection categories. For example, we found that the AVI map was the only map to identify recently burned areas as suitable for Yellow Rails. Despite evidence that fires can improve the suitability of wetland habitats for Yellow Rails Buhl 2013, Morris et al. 2017), we maintain that the classification of burned areas as "selected" is misleading in the context of our study. This is because the "burn" category in the four land cover maps used here supersedes other habitat classifications, resulting in lumping of burned forest and burned wetland habitats. Therefore, although fire may improve the suitability of wetland habitats, presumably via clearing of shrubs and woody vegetation, there is no evidence that it renders upland habitats suitable for Yellow Rails. In light of this, we believe that the blanket classification of burned areas as selected will lead to overestimation of suitable habitat. Combining the four maps into a composite map appeared to ameliorate this issue because burnt areas in the AVI map were only classified as class 8 if the other three maps similarly classified them as selected, i.e., if the underlying habitat was high quality. This example illustrates the potential benefits that can be attained by merging multiple land cover maps, which may each contribute potentially unique information to the final output.
We combined the four maps using a simple summation of the habitat selection class of the four individual maps. In the resulting composite map, the highest suitability class represented consensus high-quality habitat, but the map also had additional categories reflecting different degrees of agreement among the four individual maps. In general, this method takes N input maps with categories of 0 (avoided), 1 (neutral), and 2 (selected), and produces a map with 2N + 1 classes labeled from 0 to 2N. In our map with nine classes, occupancy analysis confirmed that the highest suitability class in this map had the highest occupancy rate (0.18), and that there was a general tendency toward higher occupancy rates in higher suitability classes (Fig. 4a).
One limitation of combining multiple maps is that doing so inhibits a detailed examination of habitat preferences. In our composite map, for example, locations were classified on a scale from 0 to 8, but this map alone lacked further information on the characteristics of the habitat on the ground. In contrast, when working with individual input maps, the relationship between land cover categories and selection classes was more explicit, as can be seen in Fig. 3. From that figure, it is evident that the land cover types with the highest selection coefficients tended to be open wetland habitats, which broadly aligns with previous descriptions of the habitat preferences of this species (Bookhout and Stenzel 1987, Martin 2012, Austin and Buhl 2013. Some of this interpretability is lost when multiple maps are combined. We propose that combining maps might be most appropriate when the primary goal of research is prediction (as in this study), rather than ecological understanding (Elith and Leathwick 2009).
A general factor that may have limited our ability to accurately predict Yellow Rail occurrence was the timeliness of the input maps. Although our test dataset used to assess the predictive accuracy of the maps was collected in 2018 and 2019, the four land cover maps were based on imagery from 2013 and earlier.
One map (LCC) was based on 2005 imagery, while another (AVI) used some aerial imagery dating back as far as the 1980s. The other two maps (DU and ABMI) used 2010 imagery. Outdated imagery may lead to a mismatch between mapped land cover types and current land cover, which likely has the effect of reducing predictive accuracy. A broad mismatch between mapped and current conditions would have resulted from the Richardson fire in 2011 and the Horse River fire in 2016, which together burned roughly 15% of our study area (Pinno et al. 2013, Wilkinson et al. 2018). In addition, oil sands developments have continued (Alberta Biodiversity Monitoring Institute 2017), which may also influence Yellow Rail habitat use and are not represented in older land cover maps. Recent innovations such as Google Earth Engine present a possible solution by allowing land cover and surface water to be mapped and updated annually (Gorelick et al. 2017). Yellow Rail occurrence is known to fluctuate from year to year, apparently in response to changing water levels (Kehoe et al. 2000), so models that can take interannual variation in conditions into account may improve predictive mapping.
In spite of the limitations described above, the success of our approach implies that when the goal is to predict species occurrence, researchers may benefit from using all available maps and combining them to produce a composite map where the highest selection category represents consensus high quality habitat. Methods such as these are expected to become more important as the number of available remote-sensing platforms and land cover classifications continues to increase.

Population estimate for Yellow Rails in the Lower Athabasca region
We applied the composite map along with our survey data and occupancy estimates to produce a population estimate for our study area. Two models with different assumptions regarding the survey radius of ARUs produced point estimates of 1650 and 2747 males, respectively. The models included as much as an order of magnitude of uncertainty, with 95% confidence intervals ranging from as low as 416 to as high as 5563 males (Table 3). This uncertainty was due to the relatively low level of survey effort relative to the size of the study area (surveys covered less than 1% of the study area), leading to a reliance on extrapolation to unsurveyed areas. Uncertainty over the occupancy rates of lowsuitability habitat was problematic. For example, we detected a single bird during surveys in class 2 habitat, which led to the occupancy rate of this class to be estimated between approximately 0 and 0.11 (Fig. 4a). This habitat class covers 18.5% of the study area, compared with just 0.8% for class 8 habitat, so uncertainty in this category contributed disproportionately to overall uncertainty in the regional population estimate.
Detections in low suitability habitat classes can occur for three reasons. First, an ARU placed in poor quality habitat might detect rails calling from high quality habitat nearby; second, high quality rail habitat may be misclassified as low-quality habitat; third, rails may occur in low quality habitat. All three of the above types of detections are plausible, if rare. A more precise population estimate could be derived in the future by visiting a larger number of previously unsurveyed locations to improve the precision of estimated occupancy rates in all habitat classes.
Even given the uncertainty in our estimate, these results significantly alter the current understanding of the distribution of Yellow Rails. The most authoritative existing population estimate suggested there were "500 or more" pairs in Alberta, based on expert opinion (Alvo andRobert 1999, COSEWIC 2009). By our estimation, there are likely more than 500 pairs in our study area alone, which covers just 7% of Alberta's total area. This is not to say that Yellow Rails are common in our study area: our estimates place the overall density of rails at no more than one male per 10 km². In contrast, Robert et al. (2004) estimated about 400 males in approximately 71 km² of wetlands in Ontario, for an average density 56 times greater than we report. The key difference between their study and ours is that their study area encompassed three patches of mostly suitable habitat, while our 48,715 km² study area was a mosaic of suitable and unsuitable habitat. The higher total population estimate in our study can be attributed to the size of our study area, rather than high average density.
Although we caution against extrapolation from our study area to the rest of the province, it is noteworthy that some known and presumed hotspots for Yellow Rails in Alberta lie outside of our study area. High densities have been reported in the Hay-Zama wildland (Prescott et al. 2002), 300 km west of our study area. Other wetland complexes, including the Peace-Athabasca delta that covers 3212 km² to the immediate north of our study area, have large amounts of apparently suitable habitat and scattered records of Yellow Rails, but have not been systematically surveyed. We conclude that the province of Alberta should be considered to contain a significant proportion of the global population of this species, which has been estimated at 10,000-25,000 individuals (Wetlands International 2019).
Our results also inform the distribution of Yellow Rails in Alberta, which has historically been thought to be concentrated south of the boreal zone. Yellow Rails were found in just seven out of 687 10 km x 10 km atlas squares (1%) that were surveyed in the Boreal forest during Alberta's first provincial breeding bird atlas (Semenchuk 1992). In the second atlas, the species was found in only 20 out of 758 surveyed squares (2.6%) in the Boreal region (Federation of Alberta Naturalists 2007). The situation is similar in Saskatchewan, where Yellow Rail records are almost entirely restricted to the southern half of the province (Smith 1996). The scarcity of records in the Boreal likely reflect a lack of targeted nocturnal surveys, rather than distributional patterns. Prescott et al. (2002) conducted standardized surveys throughout Alberta, and most of the 42 Yellow Rails detected were in the northwest of the province, contrary to the conventional wisdom at the time. Additional surveys at boreal wetlands will almost certainly reveal additional breeding locations for this species, with important implications for its conservation.

Conservation threats
We estimated that roughly 17% of the Yellow Rails in our study area occur within the boundaries of oil sands leases. These birds and their habitats likely face a higher degree of threat than those in comparable off lease areas. Among lease types, birds residing on in situ leases face less imminent threat than those on surface mining leases because in situ extraction practices leave surface habitats functionally intact, albeit fragmented by linear features used for seismic exploration (Bayne et al. 2005). In contrast, surface mining involves the complete removal of the top layer of the earth (R. B. Dunbar, unpublished data). The result is an area devoid of vegetation that is unsuitable for Yellow Rails during the period of mining. Although reclamation may recover suitable habitat in the future, the suitability of regenerated wetlands for Yellow Rails is unknown. An analysis of the reclamation plans submitted for existing surface mining leases in the oil sands region showed that after reclamation, a 36% decline in wetland habitats over predisturbance levels may occur, in addition to changes in the composition of wetland classes (Rooney et al. 2012).
Surface mining leases appear to harbor a disproportionately high number of Yellow Rails relative to their area. One likely reason for this is the existence of several surface mining leases abutting or overlapping the McLelland Lake wetland complex, which covers several square kilometers in the northern portion of our study area. We detected 78 Yellow Rails on surveys in this fen, making it the most significant known wetland for this species in our study area, with the most confirmed individuals for any single wetland complex in Alberta. Extensive proposed surface mining developments make it highly threatened. The availability of our habitat selection map should allow decision makers to prioritize important areas for mitigation or avoidance. A large portion of the McLelland complex has plans for mining to occur which will result in an absolute loss of habitat. One mine has already expanded to their lease boundary, which is directly adjacent to high quality Yellow Rail habitat. Monitoring of this area is ongoing to determine if changes in hydrology occur that impact Yellow Rail population dynamics.
Our analyses only loosely approximate the threats facing this species in the study area. As mentioned, not all areas within lease boundaries are necessarily slated for development. On the other hand, areas outside lease boundaries are also not necessarily free from human impacts. For example, oil sands extraction can affect hydrology beyond the mine footprint (Hackbarth 1980). Yellow Rails prefer sites with standing water within the narrow range from 0-30 cm in depth (Bookhout and Stenzel 1987), and their occurrence at any given wetland may fluctuate annually, a pattern that has been linked to changes in water depth (Kehoe et al. 2000;K. Drake, and L. Latremouille, unpublished data). If fluctuations in water level drive changes in site occupancy, even relatively small hydrological changes caused by mining could influence habitat quality well beyond the boundaries of a lease. Future research should prioritize understanding the effects of hydrology on this species, monitoring population trends over time, and conducting additional surveys in potentially suitable habitat, with a goal of better characterizing population parameters and mitigating threats to the long-term persistence of Yellow Rail at local and regional scales.
Responses to this article can be read online at: http://www.ace-eco.org/issues/responses.php/1538