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, 2019) 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-occupancy (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 co-occupancy models, multispecies occupancy models do not requires 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 >0.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-occupancy. 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 determineda priori and we hypothesized 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-occupancy between lynx and wildcat, lynx and wolf, and wildcat and wolf, and compared the models using AIC and biological relevance (Table 2). Due to data limitations (small sample size), we did not implement a three-species co-occupancy parameterization.