Modeling Spatially Explicit Densities of Endangered Avian Species in a Heterogeneous Landscape

Abstract

Relating population density to spatially explicit habitat characteristics can inform management by directing efforts to areas with lower densities or focusing conservation and land protection on high-density areas. We conducted point-transect surveys for the endangered Golden-cheeked Warbler (Setophaga chrysoparia) and Black-capped Vireo (Vireo atricapilla) in the live-fire region of Fort Hood, Texas. We used mark—recapture distance sampling and combined a Horvitz-Thompson estimator with a habitat-based, resourceselection gradient to estimate spatially explicit density for both species. We detected Golden-cheeked Warblers at 120 locations (202, 197, and 89 detected by primary, secondary, and both observers, respectively) and Black-capped Vireos at 173 locations (241, 255, and 107 detected by primary, secondary, and both observers, respectively). For Golden-cheeked Warblers, the average (± SE) composite detection probability estimate within a 100-m point-sample radius was 0.57 ± 0.14, and for vireos it was 0.24 ± 0.02. Estimated mean density (singing males ha^-1) was 0.14 ± 0.03 (95% confidence interval [CI]: 0.08–0.23) and 0.47 ± 0.05 (95% CI: 0.38–0.60) for Golden-cheeked Warblers and Black-capped Vireos, respectively. Our analysis suggested evidence of heterogeneity in the detection process for both species, as well as imperfect detection at distance g(0), both of which would bias estimated densities if ignored. Additionally, both species exhibited spatial variability in estimated densities, with those areas that had higher occurrence probabilities typically having higher estimated density. In the absence of spatially explicit density prediction, managers must treat all losses of potential habitat for endangered species uniformly, despite likely differences in conservation value. Our approach could be used to ascertain areas of changing density in relation to changing habitat conditions over time and space.

Resumen

Relacionar la densidad poblacional de una especie con características espacialmente explícitas de su háabitat puede ayudar a su manejo dirigiendo los esfuerzos hacia áreas con bajas densidades o enfocando la conservación y la protección en las áreas con alta densidad. Hicimos conteos de Setophaga chrysoparia y Vireo atricapilla en transectos de puntos en la región de Fort Hood, Texas. Usamos muestreo de marcado-recaptura a distancia y combinamos un estimador Horvitz-Thompson con un gradiente de selección de recursos basado en el hábitat para estimar la densidad espacialmente explícita de ambas especies. Detectamos S. chrysoparia en 120 localidades (202, 197 y 89 detecciones por observadores primarios, secundarios y juntos, respectivamente), y V. atricapilla en 173 localidades (241, 255 y 107 detecciones por observadores primarios, secundarios y juntos, respectivamente). Para S. chrysoparia la probabilidad media de detección (±EE) estimada dentro de un radio de muestra de 100 m fue 0.57 ± 0.14, y para V. atricapilla fue 0.24 ± 0.02. La densidad media estimada (machos cantando por hectárea) fue 0.14 ± 0.03 (intervalo de confianza del 95% [IC]: 0.08–0.23) y 0.47 ± 0.05 (95% IC: 0.38–0.60) para S. chrysoparia y V. atricapilla, respectivamente. Nuestros analisis sugieren que hay heterogeneidad en los procesos de detección de ambas especies, así como detección imperfecta en g(0), y ambos patrones pueden sesgar las densidades estimadas si son ignorados. Además, ambas especies exhiben variabilidad espacial en las densidades estimadas, de modo que aquellas áreas que tuvieron mayores probabilidades de presencia de las aves típicamente presentaron una mayor densidad estimada. En la ausencia de predicciones de densidad espacialmente explícitas, los gestores deberían tartar todas las pérdidas de hábitat potencial para especies amenazadas de manera uniforme, a pesar de sus posibles diferencias en su valor para conservación. Nuestra aproximación puede ser usada para determinar áreas de densidad cambiante en relación con condiciones cambiantes en el hábitat en tiempo y espacio.

ANIMAL ABUNDANCE AND distribution within landscapes or across vegetation types are standard state variables for ascertaining current population status, monitoring population changes over time, and predicting proximate and ultimate factors regulating populations (Williams et al. 2002, MacKenzie et al. 2006). Species distribution and density are focal variables of ecological research (Holt et al. 2002) and are commonly used for targeting regions and habitats for conservation and management efforts (Guisan and Zimmerman 2000, Austin 2007). For many avian species of conservation concern, broad-scale distribution may remain fairly constant, whereas local variation in distribution and density (Reijnen and Foppen 1994, Forman and Deblinger 2000) may result from altered habitat suitability due to fragmentation (Tewksbury et al. 2002), brood parasitism (Tewksbury et al. 1998), predator dynamics (Danielson et al. 1997, Kluza et al. 2000), or anthropogenic disturbances (Reijnen and Foppen 1994, Miller et al. 1998, Forman et al. 2002).

The Golden-cheeked Warbler (Setophaga chrysoparia; hereafter “warbler”) and Black-capped Vireo (Vireo atricapilla; hereafter “vireo”) are federally endangered species that breed across Texas. Warblers are thought to prefer mature mixed oak—Ashe Juniper (Quercus spp.—Juniperus ashei) woodlands with trees about 4–6 m in height and relatively homogeneous closed canopy (Ladd and Gass 1999). Vireos are thought to prefer midsuccessional, mixed-species shrubland, with vegetation that is about 0.5–3.0 m in height, has moderate woody cover with substantial breaks or openings, and is often qualitatively described as being structurally heterogeneous (Grzybowski 1995).

Within Texas, human population growth is projected to increase by 20–30%. Concomitant with growth will be an increased human— wildland interface across the breeding range of both the warbler and the vireo. Conservation planning for these species to date has been guided largely by regulatory restrictions (U.S. Fish and Wildlife Service [USFWS] 1992, Tear et al. 1995, Campomizzi et al. 2012) that generally assume anthropogenic factors to have deleterious effects on species habitat use, habitat quality, occurrence, or density (Reijnen and Foppen 1994, Chapin et al. 2000, Krausman and Harris 2001). Conservation planning has generally treated all direct effects on habitat (e.g., brush management and land development) the same or has relied on expert opinion to determine whether affected areas were high- or low-quality habitat. Without a spatially explicit prediction of bird density across space, the capacity for strategic decision making, such as siting disturbance or development to minimize negative impacts or to focus conservation effort, has been limited.

In the United States, large expanses of wildlife habitat exist on military installations, where species conservation is conducted in the context of other, sometimes competing, land uses requisite for meeting military training needs (Boice 2006). In central Texas, Fort Hood encompasses 87,890 ha and has been actively managed for the federally endangered Golden-cheeked Warbler and Black-capped Vireo since listing in 1990 and 1987, respectively (55 FR 53153; 52 FR 37420). In addition to habitat for the warbler, vireo, and other federally endangered and threatened species, Fort Hood provides infrastructure and training lands for two divisions and combat service support units totaling ∼50,000 soldiers. Thus, Fort Hood's landscape is host to a wide array of anthropogenic disturbances, including artillery fire, dismounted and mechanized maneuvers, aircraft gunnery, and aviation training. The live-fire region of Fort Hood where we conducted our study covers ∼24,000 ha (27% of the installation) where live artillery training is conducted ∼240 days year^-1, on all or nearly all of the 15 training ranges using small-arms fire to missile and aircraft ordnance. Flyovers of helicopters and other aircraft are frequent, and Explosive Ordnance Disposal [EOD] units periodically sweep the area to locate and detonate unexploded ordnance.

Information on how distribution and population density vary across identifiable habitat conditions (Kareiva et al. 1999) would enable managers to plan projects and disturbance events strategically to minimize deleterious impacts and identify areas where conservation actions should focus, allowing the military to meet training needs while optimizing conservation. As such, we conducted point-transect surveys at random locations within a systematic sample grid using two independent observers in the live-fire region of the Fort Hood. We used mark—recapture distance sampling (Laake 1999; Borchers et al. 2006; Laake et al. 2008, 2011) and mapped our spatially explicit density surface using a Horvitz- Thompson-like estimator (Borchers and Burnham 2004) to a habitatbased resource selection gradient (Farrell et al. 2013). In addition to detailing the potential benefits of mark—recapture distance sampling and the point-independence assumption for avian sampling, our work provides a more detailed understanding of fine-scale density— distribution relationships for two endangered species in a regularly disturbed environment.

Methods

Study area.—Fort Hood is located in north-central Texas (Coryell and Bell counties) in the Cross Timbers and Texas Blackland Prairies level III ecoregions. Vegetation types include pasture, grassland, mixed woodland—shrubland, and mature oak—juniper woodland. Woodland edges, shrubland, and grassland matrices often include young Ashe Juniper, Bastard Oak (Q. sinuata), Texas Redbud (Cercis canadensis L. var. texensis), Rooseveltweed (Baccharis neglecta), Yaupon (Ilex vomitoria), and other small woody species, forbs, and grasses. Our study area was the live-fire region of Fort Hood, covering ∼24,000 ha of the 87,890-ha installation (27%), where live artillery fire, aircraft gunnery, and aviation training are conducted on a near daily basis (Fig. 1).

Survey design and analysis.—We generated a point grid across the survey area that was initiated at a random starting point (Thompson 2002) within the live-fire region. The point grid was composed of 1,341 locations on a 300 × 300 m grid, and we randomly selected 453 points from this grid as the sample locations based on standard-occupancy- model sample size estimators (MacKenzie and Royle 2005). Detection estimates (Collier et al. 2010) targeted variance of the occupancy estimate at ≤6% to determine appropriate sample size (Collier et al. 2010). We used this sampling frame rather than restricting sampling a priori to areas considered habitat for either species (Wilsey et al. 2012) because recent research has expanded our understanding of what constitutes habitat for both species (Farrell et al. 2012, Klassen et al. 2012, Smith et al. 2012). To develop our selection gradient (detailed in Farrell et al. 2013), we used a single-season occupancy model with occurrence (ψ) and detection (p) parameters and a suite of species- specific predictive models using habitat metrics we hypothesized would explain warbler and vireo occurrence. We used 100-m fixed-radius point-transect surveys (Laake et al. 2011) because territory sizes of 2 to 4 ha are commonly reported for both species, making the 3.14-ha sample area biologically appropriate for modeling (but for a discussion of occupancy modeling in continuous habitats, see Efford and Dawson 2012). For each survey, we had two independent observers conduct 5-min surveys of singing males at the same time and location. For each survey occasion, detection—nondetection histories were denoted as 10, 11, 01, or 00 (detected at the point by first observer only, by both observers, by second observer only, or not detected, respectively; MacKenzie et al. 2006, Laake et al. 2011). We collapsed bird-encounter histories for each observer into detection—nondetection data at each point for occupancy estimation (Farrell et al. 2013). We surveyed each point ≤3 times over the spring breeding season of 2011, representing a maximum of 3 pairs of detection—nondetection distance surveys (MacKenzie et al. 2006, Collier et al. 2012). Previous research and monitoring has shown that once territories are established early in the breeding season, territories remain fairly stable over the course of a breeding season (Lackey et al. 2011, Campomizzi et al. 2012); thus, we assumed closure within the season although it is guaranteed that there is some restructuring among territories within a season (Betts et al. 2008). Each observer independently searched for singing males and categorized each unique detection into distance bins (0–50, 50–100 m for warblers and 0–25, 25–50, 50–75, 75–100 m for vireos). Additionally, to ensure that we did not double-count singing males, if there was any doubt that a bird might have moved during our survey period, observers were instructed to not include those observations, hence introducing some conservatism in our data collection. Observers were randomly allocated to sample locations and to primary or secondary observer order, and observers compared observations after each survey to ensure correct paired classification of all detections. Survey points were sampled from 0600 to 1200 hours, from 19 March to 16 June 2011, and all repeated surveys were separated by >9 days.

Density estimation methods that correct for bias from incomplete detection in point-transect surveys have typically used mark—recapture approaches with multiple observers or periods, or single-observer distance sampling (see review in Nichols et al. 2009). Single-observer distance sampling assumes that detection at the point or line is 1, an assumption that is seldom realized and generally results in negatively biased density estimates (Kissling and Garton 2006, Bächler and Liechti 2007). Even with mark— recapture approaches that incorporate distance as a covariate, unmodeled heterogeneity in the detection process may cause significant positive bias in detection and, thus, negative bias in density estimates (Laake et al. 2008, 2011). For example, birds that are singing loudly are more likely to be detected (“ captured”) at farther distances by multiple observers than birds that are singing quietly, and this unmodeled heterogeneity will induce a positive covariance in detection probability. Thus, for density estimation, we used a composite approach wherein the mark—recapture and distance-sampling approaches were integrated into a single analysis, which copes with the deficits of the two methods when used alone (Laake et al. 2011). Composite methods can either assume independence for all distances (full independence) or assume independence only at zero distance (point independence). The fullindependence model requires the restrictive assumption that the detection functions for the mark—recapture and distance data have identical shapes, whereas the point-independence model allows the shapes to differ. In the full-independence model the mark—recapture detection function can be distorted by unmodeled and often unknown sources of heterogeneity in the detection process, which may then induce dependence that increases with distance from the observer (Borchers et al. 2006, Laake et al. 2011).

For each species, we compared a suite of a priori candidate models for the distance and mark—recapture components of the detection functions, including sets that assumed full or point independence. For the mark—recapture component of the likelihood in both full- and point-independence models, we evaluated variability in detection probability due to distance, distance plus an additive effect of observer, and distance plus interactive effects with observer, day, or both. For the distance-sampling component of the likelihood for the point-independence models, we evaluated uniform (warbler only, with a cosine adjustment term), half-normal, and hazard-rate (vireo only, no adjustment terms) key functions with no covariates.

Estimates of mean density across the entire survey area, while useful for monitoring population trends over time (Williams et al. 2002), may not be informative for more detailed conservation planning because habitat loss or degradation often occurs at a scale of <5 ha. Thus, we estimated spatially explicit densities of both species across our study area by combining our density estimates from point-independence models with occupancy-model-based resourceselection functions previously developed for each species on our study area (Farrell et al. 2013). The best-fitting occupancy models for both species included high-resolution metrics of vegetation structure (height and canopy cover). We apportioned the resultant occupancy probabilities into 10 equal intervals (hereafter “occupancy strata”). For each occupancy stratum, we used a Horvitz-Thompsonlike abundance estimator (Borchers and Burham 2004) to estimate mean (± SE) density derived from the point-transect data within each stratum, weighted by survey effort. We then mapped these mean density estimates to the corresponding occupancy stratum to develop our density surface. Occupancy modeling by Farrell et al. (2013) was conducted using Program MARK (White and Burnham 1999) via RMARK, version 2.1.4 (Laake 2013), and mark—recapture distance-sampling analyses were conducted using R package mrds 2.1.2 (Laake et al. 2012) in R, version 3.0.0 (R Core Team 2013). An R package (Txmrds) containing the point-transect density data and analysis are available from the primary author (B.A.C.). All data and code have also been incorporated into R package mrds (Laake et al. 2012) as an example of point-based mark—recapture distancesampling analysis with independent observers and binned distances.

Results

We surveyed 453 randomly selected point-transect locations within our study area and detected singing male warblers at 120 locations, singing male vireos at 173 locations, and both species at 42 locations (Fig. 1). We detected 310 warblers (202, 197, and 89 detected by primary, secondary, and both observers, respectively) and 389 vireos (241, 255, and 107 detected by primary, secondary, and both observers, respectively).

There was little separation between the full- and pointindependence candidate models for the warbler. The best fullindependence model used a mark—recapture portion that included the variables distance to detection, survey day, and their interaction (Table 1), found no difference between observers (p₁ = p₂ = 0.47 ± 0.08), and had a combined average mark—recapture observer detection probability at the point (p[0]) of 0.69 ± 0.09. The best point-independence candidate model for the warbler used a half-normal key function with constant scale, whereas the mark—recapture model included the variables distance to detection, survey day, and their interaction (Table 1 and Fig. 2A). The mark—recapture model also did not include a difference between primary and secondary observers (e.g., p₁ = p₂ = 0.44 ± 0.11), with a combined average mark-recapture observer detection probability at the point (p[0]) of 0.66 ± 0.14. Detection probability from the distance sampling portion of the point-independence model was 0.85 ± 0.11 for an average estimate of the composite pointindependence model's detection probability of 0.57 ± 0.14 within the 100-m point-sample radius. For warblers, estimated mean density (singing males ha^-1) based on the point-independence model was 0.14 ± 0.03 (95% confidence interval [CI]: 0.08–0.23).

The best competing candidate model for the vireo was a pointindependence model with a half-normal key function with constant scale, and the mark—recapture model included observer (primary or secondary), distance, and the interaction of observer and distance (Table 2 and Fig. 2B). Observer detection probability at the point (p[0]) was 0.61 ± 0.06 for primary observers, 0.51 ± 0.06 for secondary observers, and 0.81 ± 0.05 for both observers combined, whereas detection probability from the distance-sampling portion of the model was 0.30 ± 0.02 for a composite detection probability estimate of 0.24 ± 0.02 within the 100-m point-sample radius. For vireos, estimated mean density (singing males ha^-1) was 0.47 (0.05; 95% CI: 0.38–0.60).

For the warbler, the fit of the full-independence model was roughly equivalent (ΔAIC = -1.81) to that of the best-fitting point-independence model; for the vireo, however, the full-independence model had a very poor fit (ΔAIC = 167.89) in relation to the less restrictive point-independence model, suggesting significant heterogeneity in detection of vireos in relation to distance. Thus, the choice of model will influence estimates of abundance for the two species (Tables 1 and 2, respectively). For example, estimated abundance within the covered area for the full- and point-independence models for the warbler were 464 ± 31 and 541 ± 138, respectively, with a model-averaged covered-area abundance estimate of 478 ± 64.7. Estimated abundance within the covered area for the full- and point-independence models for the vireo were 637 ± 53 and 1,586 ± 157, respectively, with a model-averaged covered- area abundance estimate of 1,577 ± 172. If we had assumed p(0) = 1 and used standard distance sampling, abundance (warbler = 364, vireo = 1,285) would have been underestimated in relation to estimates from the point-independence models as a result of overestimation of p(0) (364/541 = 0.67; 1,285/1,586 = 0.81). Had we required the fullindependence assumption for vireos, the abundance estimate allowing p(0)<1 would have been ∼50% (637/1,285) of the standard distance sampling estimate with p(0) =1 because of unmodeled heterogeneity (Laake 1999). The interaction of survey day and distance for the warbler conditional detection function implies that heterogeneity of detection probability in relation to distance increased throughout the season, as shown for two specific periods in the survey (day = 12 and 65) in which the distance slope changed from negative to positive (Fig. 3).

Occurrence models for both the warbler and the vireo were developed using LiDAR-derived metrics for vegetation height and canopy cover (Farrell et al. 2013). Thus, our resource selection function describes how spatial variation in avian distribution was related to underlying vegetation conditions across the live-fire region (Fig. 4A, C). When categorized according to the resource selection surface, both the warbler and the vireo exhibited similar variability in estimated density (Fig. 4B, D), dependent on the underlying selection gradient, with those areas that had higher occurrence probabilities typically having higher estimated density (Fig. 5A and B, respectively).

Discussion

For avian point-transect surveys, capture—recapture methods have dominated methodological development (Nichols et al. 2009). Although potential impacts of heterogeneous detection probabilities in capture—recapture studies have been well documented (Otis et al. 1978, Pollock et al. 1990), most mitigation of heterogeneity using mark—recapture methods has focused on use of identifiable covariates to reduce bias (Huggins 1989). However, as detailed by Borchers et al. (2006), observer independence does not guarantee statistical independence, and unmodeled heterogeneity induces positive covariance and, thence, positive bias in estimated detection probability. Of significant concern for avian surveys is that unmodeled heterogeneity becomes more likely when auditory cues underlie the observation process (Efford and Dawson 2009, Laake et al. 2011). As described by Buckland et al. (2010), relaxing the assumption of observer independence provides a flexible framework to assess whether full- or point-independence assumptions are reasonable for a given data set. Our results for the vireo validate concerns that population estimates that do not address unmodeled heterogeneity will be negatively biased, whereas our results for the warbler found similar levels of support for models that assumed full and point independence. In both cases, however, if we had not collected the distance data as part of our study, we would have been unable to evaluate either the full- or point-independence assumptions, and abundance would have been biased low. However, if we had collected and used distance only as a covariate under a mark—recapture design (and not also separately as its own detection function), consistent with pointtransect methods currently recommended for avian surveys (detailed in Nichols et al. 2009), the conditional detection functions would have been similar for both species and we would have concluded that distance has little effect on detection probability, an obviously incorrect conclusion (Laake et al. 2011). Additionally, for both species, our results support the contention that the standard distance-sampling assumption of p(0) =1 was untenable for avian point-count surveys (Bachler and Liechti 2007). Thus, application of a composite method, mark-recapture distance sampling, and evaluation of the full- and point-independence assumptions should allow ecologists conducting auditory point-transect surveys to better address issues associated with both g(0) and unmodeled heterogeneity in a wide variety of study designs (Laake 1999, Laake et al. 2011).

Densities of both the warbler and vireo increased in approximate concordance with occupancy probability strata, which generally supports the long-standing positive occupancy-abundance relationship common to field ecology studies (Holt et al. 2002, He and Gaston 2003). Monitoring programs using occupancy-based methods are well suited for species that are rare or widely distributed or for studies that focus on distribution rather than abundance (Bried and Pellet 2012). However, our results represent a case where occupancy models at the fine scale of our evaluation may be of limited use for population trend monitoring. Occupancy, when measured at any scale (Efford and Dawson 2012), maximizes at 1, whereas density often can and will exceed 1 on the basis of bird—habitat relationships, territoriality, or a myriad of other factors. Thus, in areas of moderate to high bird density, occupancy estimates, although necessary for distribution modeling and mapping, may not provide the complete picture necessary for population monitoring, because densities can and often will vary significantly at high levels of occupancy probability, as shown by our results.

The live-fire region on Fort Hood in which we surveyed is likely one of the most highly disturbed areas of wildlife habitat in the United States. Despite this, we observed warbler and vireo densities that were consistent with those reported outside this impact area and elsewhere in the species range (Wilkins et al. 2006, Groce et al. 2010, Cimprich and Heimbuch 2011, Peak 2011). We suggest that the differences in distribution and density across our study area are likely a response to variability in underlying habitat conditions, mainly vegetation height and structure, both of which are known to affect warblers and vireos (Ladd and Gass 1999). As such, our density surface represents a categorization of the density distribution dependent on an underlying resource selection model describing habitat complexity. Note that the north-central region of our study area was affected by a wildfire before our sampling, so although singing males of both species were detected in that area during our surveys (Fig. 1), we did not predict either occurrence or density for that region.

There are a suite of other potential assumptions that we did not address, such as assumptions regarding availability over time (Diefenbach et al. 2007, Stanislav et al. 2010), accuracy of distance measurement from auditory surveys (Alldredge et al. 2007), or count duration (Peak 2011, Bonthoux and Balent 2012). For our approach, the above approaches, and other avian point-count methods currently in use, the primary objectives are to address potential biases inherent on how sample data are collected. However, it is likely that effects of heterogeneity from these sources will be minimal in relation to the effect of unmodeled heterogeneity that we have identified here. However, we suggest that future research efforts focus on quantifying the relative contribution of different types of samplebased bias in relation to density or abundance estimation.

Although resources available to individuals can vary between occupied locations and can cause differential species-specific responses, managers addressing issues of habitat loss for endangered species must often treat all losses uniformly, despite likely differences in conservation value of different habitat locations (Kareiva et al. 1999, USFWS 2005). Our results show how densities for two endangered species with overlapping habitat requirements can vary within small spatial units across a landscape. As such, our work enables quantitative, reproducible identification of local areas of high, moderate, or low densities for two endangered species in Texas. If monitoring data under an appropriate sample design were available over time, our approach could be used to ascertain areas of changing density in relation to changing habitat conditions (MacKenzie et al. 2011), thus providing insights into habitat drivers and allowing for conservation actions, including measures of habitat loss and gain, to be evaluated in a rigorous and efficient manner. Spatially explicit approaches that incorporate local conditions into line-transect sampling are becoming more frequent in the literature (Hedley and Buckland 2004, Royle et al. 2004, Johnson et al. 2010), and methods such as ours that link avian pointsampling information to local conditions are likely to follow (Webster et al. 2008), better supporting conservation decisions.

Acknowledgments

The U.S. Army Integrated Training Area Management Program, Office of the Secretary of Defense, Department of Defense, provided funding support. S. Manning, D. Petty, J. Tatum, and M. Marshall provided additional indispensable logistical and field support. We gratefully acknowledge the Geospatial and Information Technology Research group in the Institute of Renewable Natural Resources for GIS support. We thank field technicians M. Gasner, M. Boone, A. Stuckert, I. Walker, A. Hill, S. Quasius, H. Pruett, A. Connor, J. Hennig, and S. Walters for their efforts and contribution to data collection and management. We greatly appreciate the excellent review comments from C. Handel, D. Diefenbach, and an anonymous reviewer that helped clarify several topics in our manuscript.

Literature Cited

Author notes