Density dependence influences competition and hybridization at an invasion front

Landscape and climatic change are promoting range shifts, potentially leading to competition and hybridization between formerly isolated species. However, density‐dependent interactions can impede the timely identification of associated conservation problems. The barred owl's expansion into the spotted owl's range provides a natural experiment to test for density dependence in niche overlap and hybridization in the early versus late stages of a biological invasion, thus illuminating an important biogeographical process.


| INTRODUC TI ON
Species distributions are becoming increasingly fluid in response to anthropogenic transformations of the biosphere (Dornelas et al., 2014), which may result in increasing sympatry between closely related species or species with similar ecological niches. Theory predicts that competitive exclusion will limit the coexistence of ecologically similar species (Levin, 1970). Ecological similarity manifests as niche overlap, and as niche overlap increases, the resources to which each species has exclusive access will decrease and competition will increase, potentially leading to the decline or extinction of the competitively inferior species (Ritchie, 2002;Weber & Strauss, 2016).
Hybridization may also occur and constitutes both competition for reproductive opportunities and a potential threat to the genetic integrity of either or both species. Yet the influence of population density on competition and hybridization between species in novel sympatry is poorly understood. Studying these cases will increase our understanding of the evolutionary forces that shape biodiversity and potentially serve as an early warning of conservation crises arising from biological invasions.
The intensity of interspecific competition may increase with population density when resources are limiting (Figure 1, solid line) (Muñoz & Cavieres, 2008). Early in an invasion when the colonizing species has low population density, competition will occur heterogeneously and infrequently because there are not enough individuals to occupy all preferred habitat (Anderson et al., 2013). As the density of the colonizing species increases, intraspecific competition may drive competitively inferior individuals out of preferred habitat and into more marginal areas (Holmgren, 1995). That process will eliminate low-density differences in landscape-scale habitat selection between the native and colonizing species, if any exist. As a result, distributional patterns may be dynamic during invasions, and reflect species' realized, or post-interactive, niches (Hutchinson, 1957). On the other hand, facets of the niche manifested at the level of the individual, as opposed to the population, may be less sensitive to density. Animals may face behavioural and physiological constraints in foraging habitat selection and diet from factors like perceptual range and body size (Fisher et al., 2011;Gutiérrez et al., 2007). This would result in interspecific competition whose intensity at an individual level would be fixed but whose prevalence at the population level would increase with density.
In contrast to competition, hybridization can be negatively density-dependent ( Figure 1, dashed line), because colonizing individuals may fail to find conspecific mates when their population density is low. When hybridization is possible, it may therefore facilitate invasions by allowing colonizing species to circumvent Allee effects (Yamaguchi et al., 2019). Thus, higher rates of hybridization and individual-level niche overlap during the early phases of an invasion would be indicative of species with relatively little evolutionary and ecological divergence. Such similarity would also likely be manifested as high niche overlap, and the invasion habitat selection may influence the population processes that can determine the outcome of those invasions.

K E Y W O R D S
biogeography, conservation, evolution, landscape transformation, niche overlap, secondary contact, speciation F I G U R E 1 The intensity of interspecific competition and prevalence of hybridization between related species in secondary contact may be densitydependent. If so, different niche overlap, a proxy for competition, and hybridization rates should differ between the leading edge of a range expansion and areas where the expanding species has reached high densities in novel sympatry with its congener could therefore threaten the persistence of the competitively inferior species. Yet density-dependent competition, particularly when manifested as landscape-scale niche overlap, could obscure looming conservation challenges until they become crises that are difficult to manage.
Over the last century, barred owls (BOs; Strix varia) have expanded from eastern North America into much of western North America, including the range of the northern spotted owl (NSO; S. occidentalis caurina), a process that was likely facilitated by anthropogenic landscape change (Livezey, 2009). The two species can hybridize, and the competitively superior BO poses an existential threat to the NSO (Dugger et al., 2015;Hanna et al., 2018;Wiens et al., 2014;Yackulic et al., 2019). Barred owl density in the Sierra Nevada, the core range of the CSO (CSO; S. o. occidentalis), is low but increasing , indicating that the invasion of that region is still in its early stages. California and northern spotted owls are very similar (e.g. both respond positively to old forest, Dugger et al., 2015;Jones et al., 2018;they have similar diets, Hobart et al., 2019;Wiens et al., 2014;and they interbreed Miller et al., 2017), and the competitive and reproductive dynamics of NSO and BOs have been well-documented when the latter species has reached high densities. These factors make the Sierra Nevada invasion front an ideal natural experiment with which to test the hypothesis that competition and hybridization between related species are inversely related density-dependent processes (Weber & Strauss, 2016).
For CSOs and BOs, we assessed (a) landscape-scale habitat associations using multi-scale dynamic occupancy models parameterized with data from passive acoustic surveys covering >6,000 km 2 , (b) fine-scale habitat selection using resource selection functions parameterized with GPS tag data, (c) dietary niche space with carbon and nitrogen stable isotopes, and (d) hybridization with expert phenotypic evaluation. For each metric, we also quantified overlap. We contextualized those findings with literature on niche overlap between NSOs and BOs when the latter were at high densities. We predicted that (a) landscape-scale habitat associations would be different at low but not high BO densities, indicative of low realized niche overlap early in the invasion; that (b) foraging habitat selection would be similar between species and density-independent; that (c) the BO's dietary niche would overlap and exceed that of the spotted owl and be density-independent; and that (d) hybridization would be negative density-dependent. Support for these predictions would indicate that individual-level aspects of the niche, as opposed to population-level attributes like landscape-scale habitat selection, are independent of density and therefore predict competition early in an invasion.

| Bioacoustics and occupancy modelling
We deployed autonomous recording units ( Survey sites were 400-ha hexagonal grid cells, the approximate size of CSO and BO territories in that region (Tempel et al., 2016;. These sites were randomly selected such that 20% of the study area was sampled and were non-contiguous to promote independence among sites. Sites ridges rather than gullies) without knowledge of owl occupancy history; surveys were separated by at least 4 weeks.
We derived detection/non-detection data for both species from the resulting audio data using Raven Pro 2.0 (Bioacoustics Research Program, Cornell Lab of Ornithology 2017). We developed a sliding window template detector, which calculates the correlation between audio data and a user-defined target signal, for both species using libraries of vocalizations recorded across a range of conditions in our study areas (see Wood, Popescu, et al., 2019 for further detail). We applied the template for each species to all data recorded 20:00-06:00, approximately 49,800 hr in 2017 and >145,600 hr in 2018, and manually confirmed all potential owl vocalizations. Confirmed vocalizations were considered detections, and if any of the ARUs deployed at a given site generated a detection, we considered that site to be occupied. Importantly, the randomized survey design and naïve ARU deployments meant that the occupancy of a grid cell had no a priori biological significance; sampling sites were the unit of occupancy, rather than historically active owl territories.
To allow direct comparisons between our study and those conducted in the range of the NSO where habitat conditions (and BO densities) differ (e.g. Yackulic et al., 2019), we calculated landscapescale niche overlap as the proportion of CSO sites at which a BO was also present. Earlier work indicated that CSO site occupancy in our study area was 0.43-0.53 (Wood et al., 2020;Wood, Popescu, et al., 2019), meaning that the BO population could grow to the same size as the CSO population without any overlap in site occupancy. For this reason, equating site overlap with niche overlap in the context of an increasing population  is not a tautology. We then contextualized those findings with the results of multi-scale dynamic occupancy models (MacKenzie et al., 2017;Nichols et al., 2008). In our analyses and almost all the pertinent literature, owl occupancy was measured, not density per se. However, occupancy is a reliable proxy for abundance-and thus density on the landscape-for territorial owl species (Tempel & Gutiérrez, 2013).
First, we tested whether variation in sampling effort (total ARUnights) affected, whether an owl was detected during a given secondary sampling period (θ; MacKenzie et al., 2017, pp. 276) and then whether site occupancy (ψ; resolution = 400 ha) varied by year. Next, we sequentially tested whether site occupancy, detection at a given ARU (p; resolution = 15.6 ha), site colonization (γ) and site extinction (ε) varied with habitat (see below). In each step, we compared univariate models while holding untested parameters constant, considered multivariate models if biotic and abiotic covariates (e.g. amount of open forest and elevation) both had substantial support (ΔAIC c < 2), and we then advanced the model structure with the lowest Akaike Information Criterion with a correction for small sample size (AIC c ) F I G U R E 2 Acoustic survey coverage with 2018 results (naïve occupancy = 0.47 and 0.13 for spotted and barred owls, respectively) (a), and the minimum convex polygons of the GPS data and the locations at which feather samples were collected (b) [Colour figure can be viewed at wileyonlinelibrary.com] (b) (a) (Burnham & Anderson, 2010). We used packages RMark v2.2.8 (Laake, 2013) and xlsx (Dragulescu & Arendt, 2018) in program R (R Core Development Team, 2014) for these analyses.

| GPS tagging and resource selection functions
We captured and GPS-tagged 11 CSO and 10 BOs across the north- We conducted two analyses of the GPS data. First, we tested whether foraging locations of the two species differed in their habitat composition by comparing spotted and BO foraging locations with logistic regression. Second, we developed Resource Selection Functions (RSFs) for each species following an SP-A, Design III approach (Manly et al., 2002). We then compared used and available points with logistic regression for each species. For all three sets of models (direct comparison, CSO RSF, BO RSF), we treated individual habitat variables as covariates, included a distance to territory centre (i.e. the mean centre point of all locations) term because owls are central place foragers, and treated individual owls as random intercepts to stratify comparisons by territory and account for unbalanced sample sizes (Gillies et al., 2006). We compared the seven univariate models in each of the three sets with ΔAIC c . Finally, we used the β estimates from each of those models and the R package EcoSimR (Gotelli et al., 2015) to calculate Pianka's niche overlap (Pianka, 1974).

| Stable isotope analyses
We collected body feathers from CSO (n = 29) and BO (n = 34) when they were captured as part of ongoing mark-recapture studies, or, for some BOs, when they were collected in a removal study ( Figure 2b). We used those feathers to conduct carbon and nitrogen stable isotope analyses. The stable isotope composition of a consumer's tissue reflects the average stable isotope composition of its food integrated over the time that tissue was grown. However, discerning dietary differences can be challenging for at least two reasons. Individuals consuming a diverse range of isotopically indistinct prey will have a relatively small isotopic niche. Conversely, individuals that consume isotopically distinct prey will have an average isotopic signature; if that foraging pattern is uniform across the population, there will again be a relatively small isotopic niche.
Indistinct prey and isotopic averaging can thus bias the estimated niche space low (see below), making the niche overlap between species more important than their absolute niche space. While stable isotopes lack the taxonomic resolution of pellet analyses, they can be measured efficiently over larger areas and reflect assimilated diet over much longer time periods.
Feather samples were washed three times with 2:1 chloroform:methanol, homogenized, and dried at 55°C for ≥72 hr. δ 15 N and δ 13 C analysis was conducted with a Thermo Scientific Delta V mass spectrometer at the University of New Mexico Center for Stable Isotopes (see Hobart et al., 2019 for further detail). Because all sampled owls occupied a subset of the same ecoregion, we assumed that baseline prey isotopic values were consistent across individuals.
We quantified the dietary niche for each species as the 95% prediction interval of the maximum likelihood estimated ellipse area in two-dimensional isotopic space (Jackson et al., 2012). This metric allowed us to quantify the degree to which the two species' dietary niches overlapped and, of lesser importance, determine which species' dietary niche was larger. We used the R package SIBER (Jackson & Parnel, 2017) for these analyses. We then used the mean stable isotope values for each species to calculate Pianka's niche overlap.

| Phenotypic assessments of hybridization
In 2019, we expanded our acoustic survey coverage to include >85% of the ~6,000 km 2 study area. We supplemented these surveys with vocal lure surveys for spotted owls, yielding comprehensive coverage of the study area. Following all BO detections, we conducted vocal lure surveys during which we classified birds as pure BOs or hybrids based on their plumage and territorial vocalizations. Individuals with vertical bars on their breast and horizontal barring on the nape that produced either two-phrased hoots ("who cooks for you") or inspection calls ("hoo-aw") were considered BOs; individuals with bars and spots and that exclusively produced vocalizations uncharacteristic of either parent species ( Figure S1) were considered hybrids. Birds were not classified without both visual and aural identification. We estimated reproductive overlap by dividing the number of hybrids by the number of BOs.

| Landscape-scale occupancy patterns
BOs were detected at 12% and 21% of sites at which CSO was also detected in 2017 and 2018, respectively (Table 1); co-occupied sites accounted for 50% and 73% of BO sites. The most-supported occupancy models indicated that CSO site occupancy was 0.44 and 0.48 in 2017 and 2018, respectively, while BO site occupancy increased from 0.07 to 0.13. By the late stages of the BO invasion, BOs occupied 40%-95% of historically occupied NSO territories (Yackulic et al., 2019). Across the range of the NSO from 1995 to 2013, spotted owl territory occupancy across 11 long-term study areas decreased from 0.6-0.9 to 0.1-0.5 while BO territory occupancy increased from 0.0-0.7 to 0.4-0.95 (Wiens et al., 2014;Yackulic et al., 2019).
Though the estimates from the two regions were derived with different methodology (thus yielding site vs. territory occupancy), the results corroborate the premise of our experimental design: that the relative density of BOs was high in the range of the NSO and low in the range of the CSO.  Table S1).
Overall, the two species displayed some differences in realized niche space as measured by landscape-scale habitat selection. While they shared a strong aversion to open forest, SOs and BOs preferred steeper versus shallower slopes, respectively. These findings indicate that the increasing number of co-occupied survey sites was not merely a function of increasing BO density, but rather a reflection of niche overlap. Across the range of the NSO, spotted owls were positively associated with older forests, while BOs tended to colonize sites with (montane) riparian forest and older forest (Dugger et al., 2015;Yackulic et al., 2012). NSO had a stronger preference for old forests and for steeper slopes than BOs (Pearson & Livezey, 2003).

TA B L E 1
Niche overlap between California spotted owls and barred owls, where the latter species is at low but increasing densities, and between northern spotted owls and barred owls, where the latter species has reached high densities

| Fine-scale foraging habitat selection
We  (Table S1). Pianka's overlap between the two species based on the coefficients of the RSFs was 0.809 (Table 1).
Movement data derived from VHF-tagged NSOs and BOs indicated that, when foraging, both species selected for old forest and for montane riparian forest, and selected against young forest; NSOs selected for steeper slopes and BOs selected for shallower slopes, stands near streams, and hardwoods (Irwin et al., 2012(Irwin et al., , 2018Wiens et al., 2014). Pianka's overlap between the two species based on proportional use of five forest cover types was 0.802 (Wiens et al., 2014) (Table 1).

| Stable isotope-based diet analysis
The 95% maximum likelihood estimate ellipse of stable isotope niche space of BOs overlapped almost entirely (99%) with that of CSOs, while CSOs had 52% overlap with BOs. Contrary to our prediction, the CSO had a larger maximum likelihood estimated isotope ellipse than the BO and thus occupied greater isotopic niche space (area = 26.51 and 13.83, respectively; Figure 4). Analyses of pellets indicated that although NSOs and BOs share the same key prey items F I G U R E 3 Parameter estimates (β) and model selection tables for logistic regression differentiating spotted and barred owl foraging locations (a) and resource selection functions for spotted owls (b) and barred owls (c). All covariates were standardized. All models included a distance to territory centre term and treated individual owls as random intercepts (flying squirrels [Glaucomys sabrinus], woodrats [Neotoma spp.], and lagomorphs), BOs relied less heavily on these items and consumed a wider range of species (Wiens et al., 2014). Thus, dietary overlap between spotted and barred owls was similar when BOs were at low and high densities (Table 1).

| Prevalence of hybridization
Comprehensive acoustic and vocal lure surveys of the northern Sierra Nevada study area yielded detections of 38 barred or hybrid owls. On the basis of plumage and vocalizations, the number of hybrids (10) was 36% the number of BOs (28). In contrast, whole-genome sequencing of tissue collected in regions with high BO densities indicated that hybridization between NSOs and BOs is very rare, with just 6% as many hybrids (2) as pure BOs (33) (Hanna et al., 2018).

| D ISCUSS I ON
The differences in BO densities between the range of the CSO and NSO and the ecological similarities between the two subspecies of spotted owl allowed us to test whether niche overlap and hybridization between these species in novel sympatry are densitydependent (Figure 1) These findings suggested that overlap of landscape-scale habitat selection and rates of hybridization were density-dependent, with the former increasing with BO density and the latter decreasing with density. In contrast, niche overlap was consistently high at individual-level aspects of the niche, namely for foraging habitat selection and diet (Table 1).
Consistent with our prediction, there was modest differentiation between spotted and BO's landscape-scale habitat selection in the Sierra Nevada but these differences in the species' realized niches were not apparent when BOs reached high densities.
Barred owls' preference for shallower slopes and older forest based on landscape-scale site occupancy data is consistent with patterns observed in the PNW (Pearson & Livezey, 2003;Yackulic et al., 2012). Mark-recapture studies have shown that older forest is important for CSO territory survival on comparable landscapes in the Sierra Nevada (Jones et al., 2018), suggesting that competitive interactions between the two species will increase as the BO population grows .
Behavioural data corroborate the prediction that competitive interactions across the landscape are density-dependent: the presence of BOs reduces NSO response rates (Bailey et al., 2009), but BO presence actually leads to increased rates of CSO territorial vocal activity, possibly because the BO remains a novel competitor in that region (Wood et al., 2020). The prediction that the prevalence of competition will increase as the invasion proceeds is further supported by the finding that niche overlap of both foraging habitat selection and diet was high regardless of BO density.
Given that CSO and BO foraging locations were indistinguishable within the limits of the habitat data we used (Figure 3a), it is not surprising that isotopic measurements of the diets overlapped by 0.52-0.99 (Figure 4). Yet the principle of competitive exclusion suggests that overlap will not persist: as the competitively inferior F I G U R E 4 Stable isotope signatures of spotted owls and barred owls and the 95% prediction interval of the maximum likelihood estimated ellipse area (Jackson et al., 2012) of their isotopic niches species is increasingly displaced into suboptimal habitat, individuals may be forced to adapt their foraging behaviour. Indeed, this process has already begun for the NSO (Irwin et al., 2019), and dietary changes can substantially influence spotted owl persistence (Hobart et al., 2019).
Access to reproductive opportunities may be the ultimate exclusive resource, and hybridization indicates that it is subject to competition (Ritchie, 2002). We found that BOs will more frequently hybridize with spotted owls when their own population density is low rather than high, suggesting that hybridization is the result ofand an escape from-Allee effects. When hybridization is negatively density-dependent, it is unlikely to threaten the genetic integrity of either species; it can lead to genetic extinction when hybridization is common and density-independent (Mank et al., 2004). Hybridization is indicative of relatively minor divergence of species' reproductive biology, and similarity in both foraging habitat selection and isotopic niche space suggests that ecological divergence can be correspondingly minimal.
Our results yielded some unexpected findings that warrant explanation. Apparent inconsistencies within the occupancy modelling results and between them and the RSFs may be explained by the fact that passive acoustic surveys rely on unprompted vocalizations which individuals produce in the face of potentially complex trade-offs. California spotted owls displayed a negative association (p) with medium forest at fine scales (16.5 ha), yet they selected for stands with more medium forest when foraging (Figure 3) and at a territory scale (400 ha) sites with more medium forest were less likely to go extinct. This suggests that CSO is less vocally active when foraging. This is corroborated by the finding that bouts of CSO territorial calling are shorter and contain fewer vocalizations in montane riparian forest, which is the key habitat of flying squirrels, one of their primary prey species (Wood, Schmidt, et al., 2019). Foraging habitat selection itself was also surprising: CSOs and BOs selected for young forest, rather than against it ( Figure 3). However, the amount of young forest in an owl territory is positively associated with the proportion of CSO's diet that is composed of woodrats and pocket gophers (Thomomys spp.), as opposed to the energetically less profitable flying squirrel, which, in turn, is associated with lower extinction probabilities (Hobart et al., 2019). The occupancy analyses indicated that, at the scale of entire territories (400 ha), owls were selecting for areas with less open forest, so individuals could be targeting edges between younger and medium or older stands, as has been reported for CSO elsewhere (Atuo et al., 2019), or the shrubby, closed-canopy stands that NSO has also been observed to select (Irwin et al., 2012

| CON CLUS ION
Landscape-scale monitoring programmes capable of yielding systematic data on multiple species of interest can offer a crucial early warning of range expansions ). Yet in the absence of prior information on the interactions between two species-as we had in this case-such programmes may fail to predict the intensity of competition because the landscape-scale habitat associations they reveal may be density-dependent. By the time, such programmes indicate that the decrease in one species is attributable to the increase in another, it can be too late for efficient management.
Our finding that foraging habitat selection and isotopic dietary niche breadth were density-independent suggests that comparing such individual-level aspects of species niches may better predict the potential level of competition. These ecological similarities combined with the ability for species to hybridize are indicative of fundamental similarities between species. High niche overlap could be driven by a lack of divergence between closely related species, or by convergence of more disparate taxa in the face of similar ecological constraints. Hybridization requires a range of physiological and behavioural congruencies that can persist for over a million years after divergence (e.g. Garroway et al., 2010). In these ways, the evolutionary histories of species pairs may inform their potential for conflict upon novel sympatry.
Instances of secondary contact and subsequent competition and hybridization are increasingly common (Acevedo et al., 2007;Garroway et al., 2010;Scott et al., 2019;Wood et al., 2016). their congener the Steller's jay (C. stelleri) (Smith, 1978). Even human recreation is eroding historical barriers (West et al., 2016). The vast anthropogenic transformation unfolding across the globe will undoubtedly result in many more cases of invasion, hybridization and competitive displacement. We have shown that it may be possible to predict which cases will threaten the persistence of native species by studying the individual-level traits, like foraging habitat selection and diet that influence the population processes that can determine the outcome of biological invasions.