Accounting for detection unveils the intricacy of wild boar and rabbit co-occurrence patterns in a Mediterranean landscape

The patterns of species co-occurrence have long served as a primary approach to explore concepts of interspecific interaction. However, the interpretation of such patterns is difficult as they can result from several complex ecological processes, in a scale-dependent manner. Here, we aim to investigate the co-occurrence pattern between European rabbit and wild boar in an estate in Central Portugal, using two-species occupancy modelling. With this framework, we tested species interaction for occupancy and detection, but also the interdependencies between both parameters. According to our results, the wild boar and European rabbit occurred independently in the study area. However, model averaging of the detection parameters revealed a potential positive effect of wild boar’s presence on rabbit’s detection probability. Upon further analysis of the parameter interdependencies, our results suggested that failing to account for a positive effect on rabbit’s detection could lead to potentially biased interpretations of the co-occurrence pattern. Our study, in spite of preliminary, highlights the need to understand these different pathways of species interaction to avoid erroneous inferences.


Materials and methods
Study area. The study site was the "Charneca", NE section (100 km 2 ) of Companhia das Lezírias S.A., the largest agroforestry farmstead in Portugal, with nearly 180 km 2 . This area was primarily used for silviculture and pastoral practices. Over the years, different management options shaped a complex landscape, where cork oak Montado was the primary land-cover (~66 km 2 ). The Montado occurred in pure or mixed patches and with variable composition and density of understory, depending on grazing pressures and/or shrub clearance activities. This landscape was interspersed by maritime pine stands, scrublands, olive groves, rice fields and irrigation plots 42 (Fig. 1).
The study site was also explored as a hunting estate where game species such as the European rabbit and wild boar were targeted. In Companhia das Lezírias' Forest Management reports 43 , rabbit's local population has been recognised as declining. Since 2008 several conservation measures have been implemented (e.g. scrubland management, food and water supply network, population monitoring) and hunting ceased in 2013. Contrasting, the wild boar has been growing steadily in population numbers and was widespread. As a consequence, the population was controlled through intensive hunting practices 43  Field sampling protocol to assess species presence. We surveyed rabbit and wild boar signs of presence between December 2015 and March 2016. We used a regular sampling strategy by dividing the area into 1 × 1 km plots and prospecting a 25 meters radius buffer around the centre of each plot (i.e. sites, N = 73; Fig. 1). Presence was confirmed by signs such as pellets, latrines and burrows in the case of European rabbit, and droppings, footprints and rooted area in the case of wild boar. Each buffer was prospected by the same two observers, who actively searched for these signs of presence. A species was considered present at a site if at least one of these signs was detected.
We conducted three sampling occasions and during each occasion, all sites (n = 73) were surveyed once, over a period of one to one and a half weeks (X = 9 days, SE = 1.528), depending on logistic constraints. Sampling occasions were spaced by a period of two to three weeks. All signs were registered so that in the next sampling occasion we could exclude highly persistent signs such as rabbit burrows, ensuring temporal independence between occasions. Rabbit and wild boar detections were coded as 1 for detected and 0 for non-detected, which allowed us to build site-and species-specific detection histories. For example, the detection history X l A = 001, means the location i was surveyed on three occasions, with species A being detected only in the last occasion. Overall, we had four sites with missing observations in the last two sampling occasions.
Occupancy modelling. Occupancy models make use of spatial-temporal replicated data to generate a likelihood based estimate using probabilistic arguments that account for false absences 44,45 . Site-and survey-level covariates can be incorporated when estimating both detection and occupancy parameters via a logit link 45 . Two species occupancy modelling 12,44 extends this approach allowing the probability of occupancy and detection of a subordinate species to be modelled as a function of the occupancy and detection status of the dominant species. A species interaction factor (SIF) is calculated as a ratio of how likely the two species are to co-occur compared to what would be expected under a hypothesis of independent occupancy 15 . A value <1 suggests avoidance, while values >1 indicate co-occurrence. This parameterization assumes: i) sites are closed to changes in the occupancy status of the species, ii) sites and occasions are, respectively, spatially and temporally independent, and iii) there is no un-modelled heterogeneity in both detection and occupancy. Here, we relaxed these assumptions by assuming that, during our sampling period: rabbit populations were stable, as the survey would have ended before the reproductive season 46 , and wild boar's occupancy was interpreted as probability of site use 12 , given it is a highly mobile species with large home ranges and likely random movement patterns 47 within the estate (assumption i,); furthermore, the time interval between sampling occasions (>2 weeks) mitigated the temporal dependence (assumption ii); and the inclusion of detection and occupancy covariates achieved a compromise between model complexity and interpretability (assumption iii).
We used a two-stage modelling approach 14 : firstly, we used single-species occupancy models to identify relevant covariates for both species and to avoid overparameterization in the next stage 15 ; secondly, we used two-species models, to assess interspecific interactions. For single-species models, we first tested the effect of the covariates on species detection (p) while keeping occupancy (ψ) constant; then carried the covariate from the best-fitting model and built a second set of candidate models explaining species occupancy probabilities 3 . For detectability we considered habitat features within the surveyed buffers, but also the rabbit's marking behaviour. Since in our study area European rabbit occurs in low density, large latrines are uncommon 48 and scattered pellets were more often used to ascertain species presence. We considered vegetation influenced the detection of such signs and therefore measured the following variables: shrub (S), herbaceous (H), and litter cover (L), shrub (SH) and herbaceous height (HH), and percentage of bare ground (B). For each site (N = 73), the same two observers visually estimated an average value of cover percentage and height, during the first sampling occasion. We also used Julian day (JD) as a sampling covariate as the interval between sampling occasions could have influenced species detection.
For occupancy, we selected covariates related to habitat, food availability and disturbance. The wild boar has been described as a generalist species, with preference for forested areas with dense shrub cover and high food availability (e.g. acorns) 49 . Studies on European rabbit established the importance of landscape mosaics of scrublands and pastures 50 . However, the presence of livestock poses a threat to rabbit populations, and ultimately to wild boar, as high grazing pressure leads to shrub clearance and landscape simplification 51 . The covariates tested are described in Table 1. These were extracted for a 100 meters radius buffer around each site (according, approximately, to rabbit's largest home-range size in Doñana, Spain 52 ) in software Quantum GIS (QGIS Development Team, 2016), using a GIS database available for the study area 3 . We mostly tested univariate models, but also a few with covariate combinations representing specific ecological hypothesis (e.g. food and habitat: Mont+Sparse+Cult). We assessed covariate correlation and to avoid multicollinearity excluded the less ecologically meaningful covariate when the coefficient surpassed 0.7 53 . Herbaceous cover and shrub height covariates were excluded from the analysis (see Supplementary Table S1 online). Also, the dense scrubland covariate led to convergence errors when modelling wild boar's occupancy and was thus removed for this species. Prior to the analysis, all covariates were standardized to z-scores.
In the second stage, we tested for species interaction using the best supported covariates previously identified and following the parameterization of Richmond et al. (2010) (see Table 2 for parameter abbreviations). Wild boar ("WB" indicates present; "wb" absent) was considered the dominant species and European rabbit (ER) the subordinate species. To test the hypotheses outlined in the Introduction we built four types of models: (i) both occupancy and detection probabilities are unconditional (ΨWB ΨER pWB pER); (ii) occupancy is unconditional but detection is conditional (ΨWB ΨER pWB pER rER/WB); (iii) occupancy is conditional but detection is unconditional (ΨWB ΨER/WB ΨER/wb pWB pER); (iv) both occupancy and detection are conditional (ΨWB ΨER/WB ΨER/wb pWB pER rER/WB).

Model selection.
Prior to model ranking we determined the goodness-of-fit of each species global model, using the Pearson chi-square statistic (1000 parametric bootstrap samples). Following the methodology described in MacKenzie and   54 , we used the Akaike Information Criterion corrected for small sample size (AICc) for model ranking. This approach assumes that, within the candidate set, models with ΔAICc ≤ 2 comparatively to the best model are strongly supported and their covariates are good predictors of the dependent variable. However, the use of an apparently arbitrary cut-off value has been previously criticized, and other values of ΔAICc of 4, 6 or even 8 may also include relevant models 55 . Thus, for the selection of single-species models we maintained a stricter rule of ΔAICc ≤ 2 to avoid overparameterization in the next stage. Then, for the two-species model analysis we resorted to the Akaike's weight (AICw), which indicates the weight of the evidence in favour of a certain model being the most parsimonious in the set. According to Burnham and Anderson (2002) 56 , unless a single model has a AICw >0.9, then other models are likely to be relevant and therefore should also be considered for inference. In the two-species model set, we compared AICw between models with conditional occupancy and/ or detection probabilities to compare the evidences supporting our hypothesis.
We used model averaging across the two-species model set to obtain parameter estimates of occupancy and detection probabilities. This also allowed us to compare the different hypothesis, since: if ψER/WB < ψER/wb, rabbit's occupancy probability is lower at sites where the wild boar is present, which would indicate spatial displacement (SIF < 1). If p ER > r ER , rabbit's detection probability is higher when wild boar is absent, and therefore detection is negatively influenced by wild boar's presence. The covariate effect was also ascertained by model averaging. If a covariate was not included in a given model the value was set to zero, and we calculated the weighted average using the respective AICw 57 . For the standard error (SE) we used the formula from Anderson (2008) 58 and estimated 90% unconditional c onfidence intervals (CI). We considered a well-supported effect when this interval did not overlap zero. Single-species single-season modelling was done in R statistical software (version 3.4.4) using the "unmarked" package 59 , and two-species single-season modelling was done in software PRESENCE 11.3 60 .

Results
Overall, European rabbit's naïve occupancy (i.e. proportion of sites where species was detected) was 0.36 (26 of the 73 sites), which contrasts with the wild boar with 0.80 naïve occupancy (58 out of 73 sites; Fig. 1). European rabbit and wild boar were detected together in 32% of sites, meaning rabbit was detected alone in only three of the 26 sites it was detected at.
Model averaging of the top-ranking single-species models indicated a well-supported negative effect of the proportion of pine stands on wild boar's occupancy probability and a negative effect of the distance to cultivated areas on European rabbit's occupancy (see Supplementary Tables S2 and S3). Also, litter cover had a negative effect on wild boar's detection probability and shrub cover had a positive effect on rabbit's detection, both well-supported (see Supplementary Tables S3 and S4).
To build the single season two-species candidate model set we combined the four types of models described in the "Occupancy modelling" sub-section, with these four well-supported detection and occupancy covariates. This model set comprised 16 models (see Supplementary Table S5), however, for final inference, we maintained (2020) 10:6651 | https://doi.org/10.1038/s41598-020-63492-9 www.nature.com/scientificreports www.nature.com/scientificreports/ only the top better fitted 12 models (Table 3). A few of the least supported models had convergence errors and therefore we discarded them, and all inferences were done for this subset that accounts for 90% of the explained variance. From this set, models that support some type of interaction, either for detection, occupancy or both, had a cumulative AICw of 0.569.
Our first hypothesis of a negative effect on European rabbit's occupancy by wild boar's presence was not strongly supported given independence models for occupancy probability (ψ ER/WB = ψ ER/wb ; AICw = 0.564) out-performed models that indicated spatial displacement (ψ ER/WB < ψ ER/wb ; AICw = 0.116). Furthermore, the average SIF across the model set supports a scenario of independent occupancy (SIF = 1.003 ± 0.114). Regarding our second hypothesis of wild boar's presence negatively influencing rabbit's detection probability, the evidences supported species interaction (pER <rER/WB; AICw= 0.247). For these models, rabbit's detection probability was higher when wild boar was present and averaging across the model set further supported this interaction with estimates of pER = 0.383 ± 0.068 and rER/WB = 0.489 ± 0.087 (Fig. 2).
SIF estimates for models considering rabbit's occupancy conditional upon wild boar's presence exhibited opposite signs depending on the parameterization of rabbit's detection probability. When rabbit's occupancy was conditional upon wild boar's presence, but detection was independent (AICw = 0.322), the models suggested slight co-occurrence (SIF=1.067 ± 0.276). Models that assumed simultaneously conditional occupancy and detection for rabbit (AICw=0.116), estimated an avoidance pattern instead (SIF=0.845 ± 0.217) ( Table 3).
In terms of covariate effect, model average of the beta coefficients indicated a positive effect of shrub cover on European rabbit's detection probability (β = 0.158 ± 0.757), but the distance to cultivated areas had a negative effect on species occupancy (β = −0.362 ± 0.484). For wild boar, detection probability was negatively influenced by litter cover (β = −0.413 ± 0.43) and the proportion of pine stands had a negative influence on species occupancy (β = −0.535 ± 0.639). However, the large SEs hindered the accurate assessment of the predictors' influence (see Supplementary Table S6), although they had a well-supported effect in the preliminary stage.

Discussion
Our study explored the different pathways of European rabbit and wild boar interaction in an agroforestry farmstead. This fine-scale analysis demonstrated the importance of considering interdependencies between conditional detection and occupancy, an issue often overlooked when unveiling ecological patterns.
The results from occupancy modelling supported our previous knowledge of the study area 43 : European rabbit had a low occupancy while the wild boar was widespread. Furthermore, the evidences indicated an independent occurrence pattern of European rabbit and wild boar in the study area, which did not support our first hypothesis. Wild boar's wide distribution most likely hindered predictions regarding species co-occurrence, since sites were wild boar was absent were too scarce to draw inferences. Also, the small study area and depleted rabbit population further limited the assessment of species co-occurrence.
Nonetheless, there was some support for interaction in species detection, which indicated that wild boar's presence positively influenced rabbit's detection probability. The model averaged estimates showed that rabbit detection was higher when the wild boar was present, a result that supported our second hypothesis but contrariwise to what was expected. This could have resulted from the combined effect of wild boar's influence on rabbit's behaviour and surveying success. Wild boar's rooting creates open areas which rabbits could have used more intensively for marking 61  where LSU is the number of livestock units in a plot of a given area (ha) during a certain number of days (n). Weighted average as a function of the grazing plot area within the 100 meters radius buffer. Table 1. Covariates selected to model European rabbit (ER) and wild boar (WB) occupancy probability (ψ) in Charneca (CL), according to species' ecological requirements. For each covariate an abbreviation code and brief description are presented. (2020) 10:6651 | https://doi.org/10.1038/s41598-020-63492-9 www.nature.com/scientificreports www.nature.com/scientificreports/ observer's visibility (and thus detection) in areas of dense vegetation. This effect could be greater in low density populations where scattered pellets are more common than big latrines 48 . However, further research into the factors mediating such an interaction, either through behavioural or methodological studies, are necessary.

ER/WB
Assuming this positive effect on rabbit's detection, we further analysed how this parameter influenced interpretations of species co-occurrence. For models where occupancy was conditioned but detection was not, the pattern was of slight co-occurrence; but for models where both occupancy and detection were conditioned the resulting pattern was of spatial displacement. Although these interaction patterns did not have a strong support, this suggests that failing to account for an interaction in species detection probability could lead to erroneous interpretations of wild boar's influence on rabbit's occupancy pattern. Ultimately, if species interaction for occupancy and detection probabilities have opposite signs (positive effect on detection, and negative on occupancy), not accounting for wild boar's positive effect on detection could lead to interpretations of independent occurrence or even co-occurrence. Furthermore, these interdependencies could contribute to mask similar spatial interactions between species' pairs 62 , especially for species that affect the habitat's physical structure 63,64 and hence potentially influence detectability indirectly. Several authors have demonstrated the importance of considering imperfect detections 65 , but only a few have verified an interaction in detection patterns 4,17,66 .
When modelling single-species occupancy and detection most of the covariates included in top-ranked models were well-supported. The positive effect of percentage of shrub cover on rabbit's detection can potentially reflect increased intensity of fine-scale habitat use in more vegetated patches 50 , increasing detectability. Such effect may surpass an anticipated negative influence by reduced visibility of signs of presence. For occupancy, an association to cultivated areas is well documented, since this is an important food resource for rabbit 67 , especially in Mediterranean ecosystems, characterized by strong fluctuations in annual primary production 68 . Wild boar's detection probability was lower when litter cover was high, probably because during the sampling period, dead leaf cover, especially in Montado, was abundant and could have reduced visibility of signs of presence (i.e. footprints and droppings). A lower occupancy of wild boar in sites with higher proportion of pine stands could be due to the low shrub cover which reduced suitability as a refuge. However, the influence of the landscape variables seems to be unclear, since in the final two-species model framework none of the covariates included had a   Table 3. Model ranking of single season two-species occupancy models testing the different scenarios for species interaction and according to the hypothesis presented. The models with convergence errors were discarded maintaining the sub-set that accounted for a cumulative AICw of 0.9. Model parameters are described in Table 2 and covariate abbreviations in Table 1. "WB" stands for wild boar, "ER" for European rabbit, "un" for unconditional and "cond" for conditional. (2020) 10:6651 | https://doi.org/10.1038/s41598-020-63492-9 www.nature.com/scientificreports www.nature.com/scientificreports/ well-supported effect. The most likely explanation for this outcome is that the covariates were required because the influence of the other species presence on detection and occupancy was being ignored. Once it was accounted for these covariates' significance was no longer supported.
The small-scale of our study and the current population scenario of European rabbit and wild boar in the area, limited our assessment of a co-occurrence pattern. Therefore, we recommend that field sampling design considers a landscape approach to capture a gradient of wild ungulates' and European rabbits' densities. Also, the use of species abundance instead of binary data could make this pattern more evident at larger scales. Nonetheless, our study was a preliminary approach and resorting to occupancy modelling allowed us to evaluate potential interactions in detection and occupancy, but also the interdependencies between them. The above-mentioned limitations also hindered a definitive assessment of wild boar's effect on rabbit's detection. However, if the suggested positive effect on detection probability is confirmed, then failing to account for it could falsely indicate a co-occurrence pattern. Understanding the reasons behind this interaction requires more in-depth knowledge into European rabbit's fine-scale habitat use. Besides, even if assessing co-occurrence is not of interest per se, this means one could overestimate rabbit's presence in areas of wild boar presence. Therefore, despite the limitations, this complexity of detection and occupancy interdependencies should be acknowledged in future research. Especially, for other wild and domestic ungulates with similar foraging behaviour that creates open areas, potentially influencing rabbit's detection. Furthermore, questions may also arise for any other setting in which the influence of a species promotes the detection of a second one, leading to potential bias if ignored. We reinforce the urgency of research on the potential threat by ungulates and hope that our results and recommendations will better guide future research on this topic. This is particularly relevant since the European rabbit is central to many species conservation plans in Mediterranean ecosystems 69 and has just recently been reclassified as "Endangered" worldwide by the IUCN 40 .

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.  Table 2.