Dracula’s ménagerie: A multispecies occupancy analysis of lynx, wildcat, and wolf in the Romanian Carpathians

Abstract The recovery of terrestrial carnivores in Europe is a conservation success story. Initiatives focused on restoring top predators require information on how resident species may interact with the re‐introduced species as their interactions have the potential to alter food webs, yet such data are scarce for Europe. In this study, we assessed patterns of occupancy and interactions between three carnivore species in the Romanian Carpathians. Romania houses one of the few intact carnivore guilds in Europe, making it an ideal system to assess intraguild interactions and serve as a guide for reintroductions elsewhere. We used camera trap data from two seasons in Transylvanian forests to assess occupancy and co‐occurrence of carnivores using multispecies occupancy models. Mean occupancy in the study area was highest for lynx (Ψwinter = 0.76 95% CI: 0.42–0.92; Ψautumn = 0.71 CI: 0.38–0.84) and wolf (Ψwinter = 0.60 CI: 0.34–0.78; Ψautumn = 0.81 CI: 0.25–0.95) and lowest for wildcat (Ψwinter = 0.40 CI: 0.19–0.63; Ψautumn = 0.52 CI: 0.17–0.78) We found that marginal occupancy predictors for carnivores varied between seasons. We also found differences in predictors of co‐occurrence between seasons for both lynx‐wolf and wildcat‐wolf co‐occurrence. For both seasons, we found that conditional occupancy probabilities of all three species were higher when another species was present. Our results indicate that while there are seasonal differences in predictors of occupancy and co‐occurrence of the three species, co‐occurrence in our study area is high. Terrestrial carnivore recovery efforts are ongoing worldwide. Insights into interspecific relations between carnivore species are critical when considering the depauperate communities they are introduced in. Our work showcases that apex carnivore coexistence is possible, but dependent on protection afforded to forest habitats and their prey base.


| INTRODUC TI ON
Terrestrial carnivores are some of the most imperiled species today due to their large home range requirements, high metabolic demands, sensitivity to habitat fragmentation, and persecution by humans (Crooks, 2002;Palomares & Caro, 1999;Ripple et al., 2014;Woodroffe & Ginsberg, 1998). Carnivores can also be important top-down regulators in ecological communities (Beschta & Ripple, 2009;Ripple & Beschta, 2006;Ripple & Beschta, 2012). The loss of key carnivore species can have devastating ecosystem effects (Effiom et al., 2013;Ripple et al., 2014) and changes in abundance or occurrence of carnivores can trigger trophic cascades (Ripple & Beschta, 2012). As such, the recovery of apex predators as a conservation tool to restore ecosystem functions (termed trophic rewilding) has become increasingly popular (Jørgensen, 2015;Seddon et al., 2014). Trophic rewilding is an ecological restoration strategy used to promote self-regulating ecosystems (Svenning et al., 2016).
Rewilding efforts in the context of apex predators requires not only an understanding of their ecological interactions within the carnivore guild but also the broader context of these interactions including sources of anthropogenic impacts. Many apex predators readily reestablish in human-dominated landscapes and exhibit potential coexistence with humans (Chapron et al., 2014;Lamb et al., 2020). Although the effects of apex predator recovery in natural landscapes are relatively well understood, there are significant knowledge gaps regarding the effects of their recovery in shaping species interactions (both intraguild and across trophic levels) in human-dominated landscapes (Dorresteijn et al., 2015;Kuijper et al., 2016). Interactions between carnivores are complex in nature and are integral to shaping the ecology and structure of wildlife communities. Therefore, examining such interactions in landscapes that harbor viable carnivore populations may provide important insights into the effects of carnivore recovery on the mesocarnivore communities that often dominate landscapes where apex predators have been eliminated.
Grey wolf (Canis lupus) and Eurasian lynx (Lynx lynx) are top predators in many temperate ecosystems in Europe and Asia, but their co-occurrence has been severely limited by extirpation of one species (most often wolf). This is particularly the case for most of Western and Central Europe due to a long history of human habitation and persecution of carnivore species. Both wolves and Eurasian lynx are recovering in Europe's landscapes (Chapron et al., 2014;Kaczensky et al., 2013) either through natural range expansion (wolf) or reintroductions and population augmentation (lynx). The European wildcat (Felis silvestris) is a mesocarnivore that was once common in Europe and has also been extirpated and currently at the core of reintroduction programs in some European Union states. In this context, the Romanian Carpathians represent one of the few natural areas in Europe that still harbor intact viable populations of all three species and serve as a stronghold for carnivore populations in Europe, despite anthropogenic influences common (hunting, forestry, farming, and livestock production) (Popescu et al., 2016;Salvatori et al., 2002).
While no work has been conducted on understanding the spatial relations and interactions between these three species simultaneously, research exists on pairwise interactions between species, particularly for lynx and wolf. Lynx and wolf are sympatric across most of their range and there is some diet overlap between them.
Research addressing coexistence between these species differ in their findings, but recent studies looking at spatial interactions between these species in Europe found that these two apex predators coexist and competition between them is low Wikenros et al., 2010). In Poland, lynx and wolf territories overlap and researchers concluded that the co-occurrence of these two species was facilitated by heterogeneous habitat and specialization on different prey . These predictors, habitat heterogeneity and diet, are also explaining competitive interactions between canids and felids in North America, with a lack of interference competition in heterogeneous habitat . Therefore, we expect to observe similar co-existence (high co-occurrence) and little evidence of interference competition (neutral or positive conditional occupancy values) between lynx and wolf in our study area.
Additionally, we expect to observe differences in co-occurrence based on seasonal changes in these species' behaviors. For example, the daily movement distances of male lynx are greater during the mating season (January-March) and for female lynx are greater during periods of extensive kitten care (May-August) (Jedrzejewski et al., 2002), which could cause increased interactions with wolves as lynx cover a larger geographic area during these periods. Research on wildcats is scarce, but a study conducted in the Jura Mountains of central Europe found no evidence of avoidance between lynx and wildcat (Hercé, 2011). No published research examines interactions between wildcats and wolf. Given the size difference between wolf and wildcats and their different diets, it is likely that the relationship between wildcats and wolf will be similar to that of wildcats and lynx.
In this study, we aimed to address these knowledge gaps by studying the intraguild interactions of two apex carnivores, the Eurasian lynx and the grey wolf, and a mesocarnivore, the wildcat in the Romanian Carpathians using multispecies occupancy models (Rota et al., 2016). Unlike traditional occupancy modeling, multispecies occupancy models allow for the estimation of K E Y W O R D S carnivores, coexistence, human-dominated, interactions, landscapes, multispecies occupancy

T A X O N O M Y C L A S S I F I C A T I O N
Conservation ecology; Population ecology; Spatial ecology; Trophic interactions co-occurrence probabilities for more than two species and do not assume asymmetric interactions (i.e., dominant and subordinate species). This is useful for estimating co-occurrence probabilities between species for which there is not a priori knowledge about interspecific relationships or for which there is not an obvious dominant or subordinate species. Multispecies occupancy models also allow for the estimation of marginal occupancy (occupancy of a single species irrelative of other species) and conditional occupancy (occupancy of a single species based on the presence or absence of another species) probabilities in relation to variables of interest (e.g., altitude). This approach has been used effectively to assess habitat use, interspecific interactions of carnivores in a variety of landscapes (Dechner et al., 2018;Lombardi et al., 2020;Van der Weyde et al., 2018). Previous research on lynx-wolf and lynx-wildcat interactions suggests a high capacity for coexistence, low interspecific competition, and little to no intraguild killing.
However, this research is limited and there has been no work on wolf-wildcat dynamics or interactions of lynx, wildcat, and wolf in the same region. Additionally, none of the published literature in Europe has been conducted in an area with a fully intact carnivore guild, whereas the Romanian Carpathians have viable, reproducing populations of many large carnivores and meso-carnivores that have not been extirpated (see study area). This information is crucial to understanding the effects of apex predators on mesocarnivores and the carnivore guild. By using a multispecies occupancy approach, we can analyze complex intraguild interactions and better understand competition and coexistence patterns. Results can elucidate variables and thresholds important for occurrence and coexistence of elusive species and help inform management or reintroduction efforts. Our specific objectives were as follows: (1) evaluate seasonal predictors for occupancy of each species, (2) characterize the spatial relationships (co-occurrence) of each species in winter and autumn, and (3) identify predictors that facilitate co-occurrence. Specifically, we analyzed the effects of potentially dominant apex carnivores on the occupancy and detection of a mesocarnivore to understand potential impacts reintroductions of apex predators may have on smaller carnivores. We also evaluated seasonal changes is marginal and co-occurrence probabilities to better understand how species persist and interact under different environmental conditions.

| Study area
The study area is situated in the Southern Carpathians, Romania, covering 1200 km 2 in the eastern part of the Făgăraș Mountains, Piatra Craiului, and parts of Leaota Mountains (Figure 1). The altitude of the study area ranges from 600 to 2400 m; forests cover most of the area (62%), along with a mosaic of urban-rural landscape and agriculture with significant areas of natural vegetation (22%), and alpine grasslands and subalpine shrubs (16%) (Iosif et al., 2022).

| Camera trapping and environmental variables
We divided the study area into a grid of 2.7 × 2.7 km cells ( Figure 1) and removed cells with more than ⅔ of their area exceeding 1800 m altitude and cells more than ½ of their area covered by urban landscape features. From the remaining cells, we sampled every other cell, when it was not possible to reach a selected cell, we used an adjacent cell. Each sampled cell contained a trap station, randomly located within the cell. We conducted two seasons of monitoring: (1) December 17th, 2018, to March 31st, 2019 (winter) and (2) October 9th, 2019, to January 15th, 2020 (autumn). We installed 64 camera trap stations during winter, and 76 during autumn, with high spatial overlap between seasons ( Figure 1). Each trap station had two opposite cameras installed at a height of 40 to 60 cm positioned toward animal paths. We used two camera models per trap station, a CuddeBack C1 Model 1279 with white flash for highquality color pictures in night conditions, and a Bushnell Trophy infrared camera. Camera traps were installed on animal trails along mountain ridges, mid-slopes, upper valleys, and bottom of slopes to detect carnivores at various altitudes/habitats. Camera traps were installed 1-2 weeks prior to the start of monitoring to account for additional anthropogenic disturbance from the camera installation process. We checked camera trap stations every two weeks to replace batteries and SD cards.
At each camera trap location, we recorded the presence or absences of anthropogenic disturbance (i.e., logging or settlements) as a binary variable for species detection and occurrence. We also recorded altitude (m) via GPS and extracted distance to stream (m), distance to settlement (m), and distance to roads (m) from the camera trap location using Geographic Information Systems (ArcGIS 10.7, ESRI, Redlands CA). Within a 500-meter buffer around each camera trap location (Lombardi et al., 2020), we calculated the density of local roads (km/km 2 ), the proportion of forested area and a terrain ruggedness index (TRI) (Riley et al., 1999). Full covariate descriptions and summaries are available in Table 1.

| Occupancy modeling
We implemented a multispecies occupancy model of two or more interacting species (Rota et al., 2016) in program R 3.5.1 (R Core Team, 2021) via package unmarked (Fiske & Chandler, 2011) to explore how environmental and anthropogenic variables affect the marginal occupancy (occupancy without accounting for interactions with other species), co-occurrence (overlap in marginal occupancy between species), and conditional occupancy (effects of each species presence on other species detection and occupancy) of lynx, wildcat, and wolf in the Romanian Carpathians. Unlike traditional cooccurrence models, multispecies occupancy models do not require a priori assumptions of asymmetric interactions; therefore, species were not considered dominant or subordinate to one another (Rota et al., 2016). Data from the two seasons were analyzed separately, and sessions were divided into 14-day sampling occasions, with the winter and autumn seasons having eight and seven sampling occasions respectively. Camera trap photos were cataloged by FCC staff and volunteers, and the date, time, location, and species identification were recorded for each animal detection (Iosif et al., 2022).
Covariates were checked for correlation using Pearson's correlation tests and Pearson's chi-squared test (for numerical and factors respectively), those with high correlations r > .7 were not included in the same models for the same parameter. We first explored combinations of five detection covariates for species-specific detection probabilities (Table 1) by comparing models with the same marginal occupancy parameterization for each species. Detection covariates were kept the same for all three species as we did not have a biological reason to vary them between species. We also included the latent presence/absence of every other species as species-specific detection covariates (e.g., lynx detection predicted by the presence/ absence of wildcat and wolf). Although multispecies occupancy models do not assume asymmetric interactions between species, we wanted to explore the possibility that dominant species could exist in our system and affect the presence of other species. Therefore, we also included species-specific detections of lynx as a function of the latent presence/absence of potentially dominant wolf, and wildcat as a function of lynx and wolf.  From these models, we determined a best model for each season based on Akaike information criterion (AIC), using R package MuMIn (Bartoń, 2020). We included the top detection covariates in the models exploring marginal occupancy and co-occurrence. We then ran a series of models to assess the marginal occupancy of our three species using environmental and anthropogenic variables ( Table 1) that were determined a priori and we hypothesized it would affect the marginal occupancy of each species. The candidate set of marginal occupancy models was similar for both seasons, models were only removed if variation in covariates was not great enough to allow estimation (i.e., models produced NAs or unreasonable estimates and standard errors). We compared the marginal occupancy models for each season using AIC to identify the best covariates explaining occupancy of each individual species. Using the top covariates from the marginal occupancy analysis, we ran a series of additional candidate models that reflected a priori hypotheses regarding pairwise co-occurrence between lynx and wildcat, lynx and wolf, and wildcat and wolf, and compared the models using AIC and biological relevance (Table S1). Due to data limitations (small sample size), we did not implement a three-species co-occurrence parameterization.

| Marginal occupancy
Mean occupancy for both seasons was highest for lynx (winter However, in autumn, marginal occupancy of wolf decreased with terrain ruggedness (Figure 2f), and lynx occupancy increased with forest cover (Figure 2b) while wildcat occupancy decreased with forest cover (Figure 2d).

| Co-occurrence
We also found differences in predictors of co-occurrence between seasons for both lynx-wolf and wildcat-wolf co-occupancies.
In winter, lynx-wolf and wildcat-wolf co-occurrence were predicted by forest cover (Figure 3b,c), but in autumn, co-occurrence for both pairs were predicted by terrain ruggedness (Figure 3e,f).

Lynx-wildcat co-occurrence was predicted by terrain ruggedness
for both winter and autumn seasons and was positively associated with terrain ruggedness in both winter and autumn (Figure 3a,d), but in autumn, the relationship was less linear (Figure 3d). In contrast, both lynx-wolf and wildcat-wolf co-occurrence were negatively associated with terrain ruggedness in autumn (Figure 3e,f). In winter, wildcat-wolf co-occurrence was negatively associated with forest cover while lynx-wolf co-occurrence was positively associated with forest cover, but only at >75% forest cover (Figure 3e,f).

| Conditional occupancy
In the winter season, we found that occupancy probabilities of all three species were higher when another species was present, regardless of the species (Figure 4). However, the occupancy probability of wildcat, decreased with increasing forest cover when either lynx or wolf were present (Figure 4), potentially a signal for mesopredator exclusion by apex predators in area of higher suitability. Similarly, in autumn, all species tended to co-occur, but this relationship was dependent on terrain ruggedness. Occupancy probabilities for both felids, lynx and wildcat, increased with terrain ruggedness when the other felid species was present and decreased when the other species was absent ( Figure 5). We observed the inverse relationship for both felids when considering the presence/absences of wolf, such that occupancy probabilities for lynx and wildcat decreased with increased terrain ruggedness when wolf were present and showed a positive relationship with terrain ruggedness when wolf were absent ( Figure 5). The presence of lynx and wildcat appeared to have no effect on wolf occupancy.

| Detection probabilities
For both seasons, the models that included that latent presence/absence of a potentially dominant species as a detection covariate performed significantly better than those that did not (∆AIC > 5). The top models for each season did not vary in their detection covariates; both models included distance to stream and the latent presence/absence of all species as species-specific detection covariates.
For both seasons, lynx, wildcat, and wolf detections were positively associated with the presence of the other two species (Table S1).

| DISCUSS ION
The results from our multispecies occupancy model of lynx, wildcat, and wolves in the Romanian Carpathians indicate that while there are seasonal differences in predictors of occupancy and co-occurrence of the three species, co-occurrence of the three species in our study area is high during both seasons. We identified useful predictors of marginal occupancy for each species; in winter were local road density (lynx and wolf) and altitude (wildcat). While in autumn, the best predictors of marginal occupancy were, forest cover (lynx and wildcat) and terrain ruggedness (wolf). We found that co-occurrence was influenced by environmental variables, forest cover, and terrain roughness, for both winter and autumn. Overall in this heavily forested landscape results from our study indicate that these species coexist but shift patterns of habitat use and co-occurrence seasonally.

| Determinants of occupancy
In winter, local road density was the most important predictor of occupancy for wolf, with higher road density associated with a lower probability of wolf occupancy (Figure 2e). Higher local road density in our study area is associated with higher human disturbance (e.g., limited logging) and habitat fragmentation; this corroborates findings from Jedrzejewski et al. (2004)  forests. In our study area, the proportion of forest was not an important predictor of wolf occupancy in either season, even though multiple studies have found it to be an important habitat characteristic for wolf (Jedrzejewski et al., 2004;Zlatanova & Popova, 2013) This may be due to the characteristics of our study area which is heavily forested (mean proportion forest =0.78 and 0.75 for winter and autumn monitoring sessions, respectively); thus, forest cover is not a limitation to wolf occurrence. In autumn, terrain ruggedness was the most important predictor of wolf occupancy; when terrain ruggedness index was >200 (moderately to highly rugged areas) the probability of wolf occupancy declined steeply ( Figure 2f). This can be explained by the fact that wolf's main prey source in Romania, wild boar (Sin et al., 2019), was documented to prefer less fragmented areas with large beech forest stands in autumn and early winter (Fonseca, 2008). Additionally, red and roe deer, which are also important prey for wolves, are known to move  corridors and for hunting and movement within their home range (Bailey, 1993;Bragin, 1986;Gordon & Stewart, 2007;Kerley et al., 2002;Matyushkin, 1977;Rabinowitz et al., 1987). Our results suggest that, in winter, Eurasian lynx are more likely to occupy areas with higher densities of local logging roads; these roads, which in our area are mostly unpaved, dirt roads, may provide easier access to resources within lynx home ranges due to decreased complexity of terrain and decreased snow depth/harder snowpack from vehicle travel. We did not observe this relationship with wildcat, however.
Rather, there was a negative relationship between density of local roads and wildcat occupancy in autumn (Figure 2d), which could be an artifact of body size; most documented examples of felids utilizing roads for movement within their home ranges was with larger bodied species (>11 kg). We also did not observe this relationship in winter; however, this is likely an outcome of the importance of altitude for wildcat occupancy, which has a negative relationship ( Figure 2c). Higher altitudes are associated with greater snow depth, and while lynx are well adapted to move in deep snow and altitude was not important for lynx occupancy, wildcats have physical limitations that make travel through deep snow more difficult. A study in Switzerland had similar findings whereby wildcats moved to areas free of snow in winter and spring and moved back to high elevations in summer (Mermod & Liberek, 2002). Similarly, in North America, the relationship between Canadian lynx (Lynx canadensis) and bobcat (Lynx rufus) is mediated by snowpack, with the distribution of the less snow-adapted, the bobcat, being limited by snow depth at the northern edge of its range (Morin et al., 2020;Reed et al., 2017). Our results for marginal occupancy of lynx, wildcats, and wolf provide insights into both habitat selection and spatial relations for these elusive carnivores in Romania. Our results suggest lynx may use roads for movement, a practice common for other felids of similar body size, but not described in this species. Additionally, we provide further support for previous findings on habitat selection and occupancy for these three European terrestrial predators.

| Determinants of co-occurrence
In winter and autumn, co-occurrence for lynx and wolf was high indicating that both species have similar habitat requirements. In winter, we found an effect of forest cover on the co-occurrence of lynx and wolf; co-occurrence increased with proportion of forest cover >0.75.  prey of lynx and wolf. Roe deer abundance is also lower in areas with high forest cover (Melis et al., 2009). Higher lynx-wolf co-occurrence in areas expected to have lower roe deer abundance indicates that lynx and wolf are likely partitioning prey resources which would reduce competition. In our study area, wolf also prey on wild boar and red deer (Sin et al., 2019). In autumn, terrain ruggedness was a negative predictor of co-occurrence for lynx and wolf, such that predicted co-occurrence was ~0 for the highest values of terrain ruggedness.
This relationship is driven by the negative relationship between marginal occupancy for wolf and terrain ruggedness, which is related to prey movements and availability as explained above (Fonseca, 2008;Sin et al., 2019) (Figure 2c). Because marginal occupancy for wolf is ~0 at high terrain ruggedness, co-occurrence for lynx and wolf is low as well. Additionally, co-occurrence between wolf and wildcat decreased with terrain ruggedness in autumn (Figure 3f) due to the low marginal occupancy for wolf at high terrain ruggedness. In winter however, co-occurrence of wolf and wildcat was predicted by proportion of forest such that increasing forest cover resulted in lower co-occurrence (Figure 3c). In both seasons, the co-occurrence of lynx and wildcat increased with terrain ruggedness, but the relationship was stronger in winter (Figure 3a,d). This relationship also provides further evidence that the negative relationship observed for lynx and wolf co-occurrence and terrain ruggedness was driven by wolf marginal occupancy.

| Management and conservation implications
The positive effect of wolf and lynx presences on detection of one another, high levels of co-occurrence in winter, and high levels of conditional occupancy in both seasons (higher occupancy probability when other species is present), for lynx and wolf provide little evidence of interference competition between these apex predators. This suggest that carnivore species may aggregate in certain habitats during winter, potentially driven by prey availability. This corroborates findings from other studies assessing interactions between co-occurring felids and canids that overlap in resource use. not affect the introduction efforts given that prey base can support both species, and releases occur in highly forested but less topographically fragmented areas. Additionally, our findings also suggest that apex predators have little negative effects on the mesocarnivore, wildcat. This information is useful for management given that wolves are recolonizing their former range in Europe (Chapron et al., 2014). Our findings suggest that wolf would not have negative impacts on wildcat given enough suitable habitat is available. In summary, studying intraguild interactions in an intact system has enabled us to observe and quantify intraspecific interactions among carnivores where they have co-existed and co-evolved for centuries. This provides insight into their potential long-term dynamics for areas where they are recovering naturally or recovering through rewilding efforts. While our study did not include the summer season, our results from two separate and partially overlapping autumn and winter seasons suggest that competition between lynx, wildcat, and wolf is low. However, additional information on the richness and abundance of the prey base, and the spatial and temporal relations between predators and their prey can augment these findings and provide additional management insights in the context of rewilding.

ACK N OWLED G M ENTS
We thank Piatra Craiului National Park Administration and the Hunting Associations Bârsa, Jderul, and GTS Muntenia, for permissions to undertake fieldwork. We thank Liviu Ungureanu, Călin Șerban, and rangers of the Foundation Conservation Carpathia for help with camera deployment and checking. We thank Ken Kellner for continued support with R code for the multispecies occupancy  Travel for MD to Romania was provided by the Ohio University College of Arts and Sciences.