Co-occurrence patterns and habitat selection of the mountain hare, European hare, and European rabbit in urban areas of Sweden

Assessing the underlying mechanisms of species co-occurrence patterns can be challenging as biotic and abiotic factors are hard to disentangle. To date, few studies have investigated co-occurrence patterns of mammals within urban areas. As urban areas are increasingly used as habitat by wildlife, there is a need for a better understanding of urban ecology to facilitate human-wildlife co-existence. Here, we investigated co-occurrence patterns and habitat selection of the European hare (Lepus europaeus), mountain hare (L. timidus), and European rabbit (Oryctolagus cuniculus) inside urban areas of Sweden, using joint species distribution models and generalized linear mixed models based on citizen science observations. All three species were observed within urban areas, but European hares and rabbits appear to be more successful urban colonizers compared to mountain hares. Overall, our findings suggested that urban occurrence by all three lagomorphs was related to suitable conditions within the distribution of each species (e.g., climate and elevation), rather than by the presence of other lagomorph species or specific land cover types within urban areas. On a finer spatial scale, European hares and rabbits generally selected for green urban areas and mountain hares for residential gardens, which likely constitute suitable foraging sites. Moreover, overlap in activity times between European hares and rabbits was mediated by land cover type and sympatry. Our findings contribute to the understanding of urban ecology and provide insights for management measures of the three lagomorphs in urban areas of Sweden.


Introduction
Studying the underlying mechanisms of species co-occurrence and interactions can be challenging, because disentangling abiotic and biotic factors affecting the occurrence and abundance of species is difficult, especially in heterogeneous environments. However, studies on co-occurrence and new methods, which can untwine abiotic and biotic factors, have received more attention in recent years. For example, niche differences, distinct habitat preferences, competitive exclusion, environmental filtering, or a combination of these factors were proposed as ecological explanations for species segregation or co-occurrence (Bar-Massada 2015; Estevo et al. 2017;Kohli et al. 2018;Pollock et al. 2014;Ulrich et al. 2018).
Urban areas might lead to altered species co-occurrence and interactions because they restructure animal communities (Brown 2001;Grimm et al. 2008). Both the expansion of urban areas into animal habitats as well as active urban colonization leads to the increasing occurrence of wildlife in these novel habitats (Luniak 2004). The causes of urban colonization are often not well understood. For example, urban colonization might be driven by poorer habitat conditions or increased hunting pressure outside urban areas (Mayer and Sunde 2020;Rutz 2008). Thus, urban areas can constitute an advantageous habitat, e.g., due to relaxed predation (Møller 2012) or increased resource availability (Contesse et al. 2004). While some species proliferate in urban areas, others are not able to adapt and become locally extinct (McKinney 2006;Shochat et al. 2006). Consequently, an increased understanding of habitat preferences by urban wildlife can be a valuable tool to aid conservation actions, ensuring suitable habitats for urban colonizers. At the same time, it is important to consider biotic interactions, as competition between species might cause exclusion from otherwise suitable habitats (Thulin 2003). While habitat selection within urban areas has been previously addressed in numerous mammal species (Bozek et al. 2007;Chambers and Dickman 2002;Duduś et al. 2014;Mayer and Sunde 2020), research on urban community ecology and species interactions are scarce (Carrete et al. 2010;Magle et al. 2012;Ramírez-Cruz et al. 2019).
Using citizen science observations, we investigated the factors affecting co-occurrence patterns and habitat selection of the European hare (Lepus europaeus), mountain hare (L. timidus), and European rabbit (Oryctolagus cuniculus, hereafter rabbit) in urban areas of Sweden. Sweden was selected as a case study due to its high quantity of citizen science data, and for harboring three lagomorph species of similar ecology, providing an optimal model to investigate species interactions (Leach et al. 2015).
Both European hares and rabbits now occur in urban areas (Mayer and Sunde 2020;Ziege et al. 2020) that, under certain conditions, appear to constitute advantageous habitat. For example, rabbits became more diurnal, spent less energy on anti-predator behaviors, and reduced their home range size, possibly due to increased resource availability (Ziege et al. 2016(Ziege et al. , 2020. There is little information regarding urban colonization by mountain hares, although some urban and suburban observations exist (Haigh and Lawton 2007;Levänen et al. 2019).
Mountain hares are native to Sweden, typically associated with tundra, open forest, and heathland in upland areas (Flux and Angermann 1990;Thulin 2003). Moreover, their occurrence is positively associated with deep and lasting snow cover, and negatively with human influence (Jansson and Pehrson 2007;Leach et al. 2016). Mountain hares are declining and categorized as near threatened in Sweden (Artdatabanken 2020), with milder winters and competitive exclusion by European hares expanding their distribution northwards proposed to be responsible for this decline (Jansson and Pehrson 2007;Thulin 2003). Both European hares and rabbits were introduced to Sweden (Artdatabanken 2020). European hares are typically associated with agricultural lowland (Mori et al. 2022), and their densities have been found to be positively correlated with higher temperatures and lower precipitation (Leach et al. 2016;Smith et al. 2005). While European hare populations have been declining in large parts of Europe since 1960 due to agricultural intensification (Smith et al. 2005), they might still be expanding their distribution in Sweden (Jansson and Pehrson 2007). The rabbit is categorized as 'near threatened' in its native range in the Iberian Peninsula, but introduced populations, e.g., in Australia and New Zealand, have proliferated before the introduction of biological control measures (Cooke and Fenner 2002;Lees and Bell 2008). Although flexible in their habitat preferences, rabbits are predominantly found in grassland, pastures or arable land bordering scrubland, providing cover from predators (Calvete et al. 2004;Tapia et al. 2014). They prefer sandy soil that allows them to dig burrows, and their distribution is positively correlated with temperature and negatively with precipitation and slope (Calvete et al. 2004;Leach et al. 2016). All three species are game species in Sweden, with European hares and rabbits being regulated in areas where they might cause damage (https:// jagar eforb undet. se/).
Previous studies are not in compliance on European hare and rabbit interactions. Most studies have found no or limited evidence for competition between the two species (Flux 2008;Katona et al. 2004;Stott 2003), with one study suggesting facilitation (Leach et al. 2017). However, an assessment of the effectiveness of different measures to eradicate rabbits from islands showed that European hares were markedly more effective than both cats and myxomatosis in removing rabbits due to competitive exclusion (Flux 1993). The distribution of the mountain hare, apart from abiotic factors (Leach et al. 2016), might be affected via competitive exclusion by the European hares' northward expansion (Bedson et al. 2021;Thulin 2003).
Here, we first described patterns of urban occurrence by the three lagomorphs and compared them to non-urban population trends obtained from hunting bag statistics to assess their urban colonization success. We then used joint distribution models to investigate the underlying mechanisms (i.e., environmental filtering or biotic interactions) of the three species' co-occurrence patterns. Moreover, we investigated species occurrence and habitat selection within urban areas, assessing the role of urban area size, climate, elevation (occurrence analysis only), urban land cover types and occurrence of the other lagomorph species. We predicted that urbanization might increase competition for resources between European hares and rabbits, which should lead to segregation of the two species within urban areas. We further predicted that mountain hares segregate from both European hares and rabbits, due to environmental filtering, given the mountain hares' distinct habitat and climatic preferences (occurring in colder and more forested areas compared to the other two species), and potentially due to competitive exclusion. Regarding habitat selection, we predicted that the lagomorphs selected land covers that resemble those of their preferred habitats outside urban areas, i.e., European hares and rabbits selecting open herbaceous vegetated areas, e.g., green urban areas and residential lawns (and rabbits additionally for sandy soils), and mountain hares selecting forested areas.

Study areas and preparation of spatial data
Our study area comprised Urban Morphological Zones (UMZ) of the CORINE Land Cover 2000 version 16, defined as areas within 200 m of each other considered to contribute to the urban tissue and function (https:// www. eea. europa. eu/ data-and-maps/ data/ urban-morph ologi cal-zones-2000-2), within Sweden, obtained from The European Environment Agency (EEA) (http:// ftp. eea. europa. eu/ www/ umz/ v4f0/ UMZ20 00. zip). Because higher human population densities increase sampling effort, thereby reducing the number of pseudo-absences and the effect of spatially biased sampling effort in point occurrence data (Geldmann et al. 2016;Mair and Ruete 2016), we only considered UMZ's > 10 km 2 (hereafter urban areas) for our analysis, leaving 97 urban areas (Fig. 1A). Moreover, we created 1 × 1 km grid cells within urban areas using ArcGISPro 2.8.3 (Esri Inc. 2020), resulting in 4915 grid cells, to analyze species associations on a finer spatial scale (see below).

Citizen science observations and hunting bag data
Point occurrence data (human observations of live and dead animals) for the three lagomorph species within Sweden were derived from the Global Biodiversity Information Faculty (GBIF) (GBIF.org (21 March 2022) GBIF Occurrence Download https:// doi. org/ 10. 15468/ dl. du6h5m) for the years . There were very few observations before 2007 (likely due to limited citizen science engagement before that time), which is why we used this year as cut-off. GBIF is a database providing institutions from all over the world with common standards and open-source tools for sharing information on when and where species have been recorded (https:// www. gbif. org/ what-is-gbif). The bulk of data in this study (98.5%) came from Artportalen (https:// www. artpo rtalen. se/ Home/ About), a website for reporting species in Sweden. To reduce variation in data quality and ensure a certain degree of location precision, observations with a spatial uncertainty of > 1000 m were excluded. We intersected all observations with urban areas and grid cells to assign them to environmental variables (see below), using the R package 'raster' (Hijmans et al. 2015).
Moreover, we used hunting bag data, derived from the Swedish Hunters Associations databank Viltdata (https:// rappo rt. viltd ata. se/ stati stik/), as a relative measure of nonurban population trends for the three species. We compared hunting bag data with the proportion of urban observations (compared to all observations) for each species separately for each year, to assess concurrence and temporal patterns in urban colonization. Changes in human urban population were obtained from the World Bank database (https:// data. world bank. org/ indic ator/ SP. URB. TOTL. IN. ZS? locat ions= SE).

Environmental data
We obtained environmental data known to and/or suspected to affect the occurrence of the three species (Schai-Braun et al. 2021;Smith et al. 2005). We obtained climate data (i.e., annual mean temperature, mean temperature of the coldest quarter, and annual precipitation) from 1970 to 2000 at 2.5 arc-minute resolution from WorldClim version 2.1. (https:// www. world clim. org), and mean soil sand content and bulk density at 250 m resolution from the International Soil Reference and Information Centre version 2.0.1. (https:// maps. isric. org/). We selected a depth of 15-30 cm, because rabbit burrow depth is 20 cm on average (Serrano and Hidalgo de Trucios 2011). Elevation data at 90 m resolution were obtained from the Shuttle Radar Topography Mission (http:// srtm. csi. cgiar. org/). To extract climate, soil, and elevation information for the urban areas, we created 100 random points per urban area using ArcGIS Pro2.8.3. We extracted the environmental values from each raster layer using the R package 'raster' (Hijmans et al. 2015) and values were averaged for each urban area. CORINE land cover vector data were intersected with urban areas and urban grid cells to calculate the area of each land cover patch. Moreover, to describe the land cover surrounding urban areas, we buffered each urban area by 1000 m, and then intersected this buffer with the land cover vector. We reclassified the CORINE land cover classes into 8 categories within urban areas: (1) continuous urban fabric (e.g., city centers, > 80% of the ground covered by artificial surfaces, i.e., soil sealed), (2) discontinuous urban fabric (e.g., suburbs, > 30% scattered urban fabric without sealed soil), (3) industry (including airports, railways, etc.), (4) green urban areas (e.g., parks), (5) agriculture, (6) forest (including other (semi)natural areas including heathland), (7) water, and (8) other areas (beaches, bare rock, etc.) ( Table S1). The land covers surrounding urban areas were categorized into (1) agriculture, (2) forest, (3) urban areas (merging the abovementioned urban categories due to little urban land cover surrounding urban areas), and (4) other areas. We then calculated the proportion of each land cover category per urban area, grid cell, and surrounding urban areas.

Defining species occurrence
We categorized a species occurring in an urban area when there were ≥ 7 observations within the urban area, i.e., at least one observation per year when most observations were recorded (from 2015 to 2021; see results). We chose this categorization to minimize defining species occurrence based on misidentifications, observations of escaped/released pet hares and rabbits, or dispersing individuals that had not established in the area. Defining occurrence based on a single observation did not markedly change the results (not shown). On grid cell level, we defined species presence in grid cells with ≥ 1 observation, because on this fine scale point occurrence data likely was more prone to false-negatives rather than false-positives (Crall et al. 2011).

Data analyses
First, to assess whether there was evidence of biotic interactions between the three species and whether the three species shared or had distinct environmental affiliations, we used the joint species distribution model (JSDM) provided by Pollock et al. (2014), which accounts for co-occurrence patterns of multiple species. We modeled predicted probabilities of occurrences, using a binary response variable (i.e., presence/absence), using a multivariate probit regression model (Pollock et al. 2014). The model provides environmental correlations for species pairs indicating shared or differing environmental responses, while residual correlations suggest biotic interactions, such as competition or facilitation (Pollock et al. 2014). We did two model runs, one on an urban area level and one on a grid cell level, as environmental effects and competitive interactions are known to appear at different scales (Leach et al. 2017). For the JSDMs, collinearity among environmental variables was assessed building a correlation matrix, defining correlation as Pearson's coefficient > 0.6 (Zuur et al. 2010). We removed mean annual temperature (positively correlated with temperature of coldest quarter), soil bulk density (positively correlated with soil sand content) and the proportion of surrounding forest (negatively correlated with surrounding agriculture). Moreover, we removed the proportion of forest (positively correlated with surrounding forest and negatively with surrounding agriculture) for the urban area level analysis and the proportion of discontinuous urban fabric (negatively correlated with industry) for the analysis on grid cell level. Consequently, we included temperature of coldest quarter, soil sand content, precipitation, elevation, the proportion of continuous urban fabric, discontinuous urban fabric (for the urban area level analysis only), forest (for the grid cell level analysis only), industry, green urban areas, agriculture, surrounding agriculture and surrounding urban fabric. We centered and scaled covariates prior to analyses to obtain comparable estimates (Grueber et al. 2011).
Second, we investigated the factors affecting species occurrence within urban areas, separately for the three lagomorphs within their distribution, using generalized linear models with a log link and a binomial response distribution (present = 1 versus absent = 0). To estimate the European hares' and rabbits' distribution in Sweden (the mountain hares' range covers the whole of Sweden), we created 100% minimum convex polygons based on citizen science observations (excluding obvious outliers). We initially built 4 candidate models based on biological hypotheses: (1) land covers within urban areas affect the presence of a species, including the proportion of continuous urban fabric, discontinuous urban fabric, green urban areas, industry, forest, and agriculture; (2) land covers surrounding urban areas affect species occurrence, including the proportion of surrounding agriculture, forest, and urban areas; (3) climate and the size of an urban area (as proxy for the number of observers) affect occurrence, including mean temperature of coldest quarter, mean annual precipitation, elevation, soil sand content (for rabbits only); (4) competition with or facilitation by other lagomorphs affects occurrence, including the presence of other lagomorph species (estimated as above).
The proportion of surrounding forest and agriculture were highly correlated (Pearson's correlation coefficient > 0.6 and variance inflation factor > 3 (Zuur et al. 2010)) so we only included agriculture (European hares and rabbits) or forest (mountain hares) in our analysis. Additionally, as measure of relative abundance, we analyzed the number of observations per urban area separately for each species (again including all urban areas within the species' distribution), using generalized linear models of the R package 'glmmTMB' (Magnusson et al. 2017) with a log link function and negative binomial distribution to account for overdispersion and zero-inflation (O'hara and Kotze 2010). We again built the same candidate models as for the species presence analyses. We scaled all numeric variables (mean = 0; standard deviation = 1) to obtain comparable estimates. We initially compared the 4 models based on biological hypotheses for both the analyses of species presence and relative abundance using Akaike's Information Criterion (AIC). To obtain the most parsimonious (hereafter best) model, we performed a stepwise backward selection, starting from the full model including all variables, and removed variables that lead to an increase in AIC, selecting the model with the lowest AIC (Wagenmakers and Farrell 2004).
Third, to analyze habitat selection within urban areas, we selected urban areas that had at least 10 observations of a given species. We excluded observations with spatial uncertainty of > 500 m, because this analysis was conducted at a finer spatial scale. To get a measure of resource availability, we created 5 × the number of random positions than we had obtained from citizen observations within each urban area (Mayer et al. 2019). We then assigned each random and used (observed) position to the land cover type (as defined above) and the soil sand content (for rabbits only). To analyze habitat selection (observed location = 1 versus random location = 0, dependent variable), we used generalized linear mixed models with a binomial distribution and a logit link, using the R package 'lme4' (Bates et al. 2015). We included the land cover type (excluding 'water' and 'other' land cover), soil sand content (for rabbits only), the presence of other lagomorph species, and the interaction of land cover type with lagomorph presence (to investigate if habitat selection differs in the presence of other lagomorphs) as fixed effects, and urban ID as random intercept to control for non-independence of the data. We again performed a stepwise backward selection, starting from the full model including all variables, selecting the model with the lowest AIC. Parameters that included zero within their 95% confidence interval were considered uninformative (Arnold 2010).
Finally, we compared temporal overlap in observations between European hares and rabbits in urban areas where they were sympatric versus allopatric, and separately for the main land cover types where they were observed (discontinuous urban fabric, green urban areas, and industrial areas). This latter comparison was only possible in urban areas where European hares and rabbits were sympatric (sample sizes were too low in urban areas where the two species were allopatric). To do so, we converted the time of day into radians (Wang and Fisher 2012) and calculated observation density overlap Δ between European hare and rabbit observations, separately for urban areas where they occurred sympatric and allopatric (as defined above) (Meredith and Ridout 2014; Ridout and Linkie 2009). We used the estimator "Dhat4" to calculate the coefficient of overlap (Δ4), because all groups had > 50 observations (Meredith and Ridout 2014). Coefficient of overlap values range from 0 to 1, with 0 representing no overlap and 1 being total overlap. We obtained 95% confidence intervals (CI) from 500 bootstrap samples. All analyses were carried out in R4.0.3.

Patterns of urban occurrence
Out of a total of 20,931 observations (both within and outside urban areas), European hares constituted 12,492 observations (60%), mountain hares 5727 observations (27%), and rabbits 2712 observations (13%). The number of observations increased over time, with few observations before 2015 ( Fig. 2A). Of the three species, rabbits had the highest proportion of urban observations, accounting for 39% (1049 observations) of all rabbit observations. Most observations were recorded in the morning (around 07:00) and evening (around 19:00), though this pattern was stronger for the two hare species compared to rabbits, and the fewest observations were recorded between 23:00 and 04:00 (Fig. 2B). The proportion of urban rabbit observations fluctuated between years, with noticeable decreases in -2009-2015. For European hares, 22% (2769) were urban observations, with the proportion of urban observations being relatively stable over time (Fig. 2B). For mountain hares, urban observations accounted for 12% (714 observations), and the proportion of urban observations increased over the years, being 2.4% in 2007 and 21% in 2021 (Fig. 2B). Hunting bag numbers decreased for European hares and mountain hares, and fluctuated for rabbits, with pronounced increases in 2009 and 2015 (Fig. 2D). The percentage of people living in urban areas increased from 85% in 2007 to 88% in 2020.

Environmental and residual correlations
European hares and rabbits shared environmental responses (mean temperature of the coldest quarter, annual precipitation, soil sand content, elevation, and land cover proportions), while mountain hares had distinct environmental responses from the other two species (Table 1). These responses were more pronounced on a grid cell level. All three species pairs had positive residual correlations (especially European and mountain hares, and rabbits and mountain hares), suggesting that all species pairs cooccurred more than expected, due to unmodelled factors (Table 1). For European and mountain hares, and rabbits and mountain hares, this pattern was more pronounced on  urban level compared to grid cell level, but for European hares and rabbits residual correlation was stronger on grid cell level (Table 1).

Species occurrence and observations per urban area
Within their respective ranges (Fig. 1), 63 of 77 urban areas (82%) contained European hare observations, 38 of 97 urban areas (39%) contained mountain hare observations, and 40 of 69 urban areas (58%) contained rabbit observations. When defining species presence as at least 7 observations within an urban area, European hares occurred in 45% of urban areas within their distribution, mountain hares in 10%, and rabbits in 26% of urban areas. For all three species, the probability of occurrence within an urban area was best explained by the model including climate variables, elevation, and the size of the urban area, followed by the model including the surrounding land cover (European and mountain hare analysis) or other species (rabbit analysis), and finally the model including land cover within urban areas (Table 2). After model selection, the best model explaining the probability of urban European hare occurrence included urban area size (positive correlation), elevation (positive correlation; Fig. 3A), mean annual precipitation (negative correlation; Fig. 3B), and temperature of the coldest quarter (uninformative positive correlation; Table 3). The probability of urban mountain hare occurrence also increased with urban area size and declined with increased temperature of the coldest quarter (Fig. 3C, Table 3). Elevation was included in the best model but was uninformative (positive correlation). The probability of urban rabbit occurrence also increased with urban area size, and with the proportion of green urban areas (Fig. 3D), though this effect was uninformative (Table 3).
The number of observations per urban area ranged from 0 to 846 (mean ± SD = 36 ± 116, median = 5) for European hares, from 0 to 461 (mean ± SD; 8 ± 48, median = 0) for mountain hares, and from 0 to 226 (mean ± SD; 15 ± 40, median = 1) for rabbits. The number of European hare observations per urban area was positively correlated with the size of the urban area, the proportion of forest, continuous urban fabric, surrounding agriculture, and rabbit presence, and negatively correlated with increasing precipitation and temperature of the coldest quarter (Table S2, Table S3, Fig.  S1). Urban mountain hare observations were positively correlated with urban area size, elevation, the proportion of surrounding urban areas, and European hare presence, and negatively with the proportion of discontinuous urban fabric and temperature of the coldest quarter (Table S2, Table S3, Fig. S2). Proportion of agriculture was included in the best model (positive correlation) but was uninformative. Urban rabbit observations were positively correlated with increasing urban area size, soil sand content, proportion of discontinuous urban fabric, green urban areas, industry, and the proportion of surrounding urban areas, and negatively with increasing elevation and proportion of urban forest (Table S2, Table S3, Fig. S3).

Habitat use and selection within urban areas
Based on random positions (located within urban areas where lagomorphs were present), urban areas were dominated by discontinuous urban fabric (61%), followed by industrial areas (20%), green urban areas (13%), forest (3%), continuous urban fabric (< 2%), and agriculture (< 2%). All three species were mostly observed in discontinuous urban fabric (especially mountain hares), followed by green urban and industrial areas (Fig. 4A). For any species, < 5% of observations came from continuous urban fabric, forest, and agriculture combined. Further, there were more European hare and rabbit observations in discontinuous urban fabric when they occurred in absence of the other lagomorph species (Fig. 4A). The opposite was the case in green urban and industrial areas, i.e., more European hares were observed in green urban areas when rabbits were present, and more rabbits were observed in industrial areas when European hares were present (Fig. 4A).
Habitat selection by European hares differed between urban areas where rabbits were absent versus present (Table 4). European hares selected for green urban areas and avoided discontinuous urban fabric when rabbits were present but avoided green areas and selected for discontinuous urban fabric when rabbits were absent (Fig. 4B). Moreover, they showed no clear selection or avoidance of continuous urban fabric and agriculture when rabbits were present but avoided these land covers when rabbits were absent (Table 4). They consistently avoided industrial areas and forests independent of rabbit presence (Fig. 4B). Only 6 urban areas had at least 10 mountain hare observations, all located outside the distribution of European hares and rabbits. Mountain hares selected for discontinuous urban fabric, and avoided green urban and industrial areas, and forests (Fig. 4C, Table 4). We removed the continuous urban fabric and agriculture from this analysis because there were no mountain hare observations in these areas, and they constituted a negligible portion of the area (< 1%). Habitat selection by rabbits was not affected by European hare presence. Rabbits selected for green urban areas, showed no clear selection or avoidance of continuous urban fabric, and avoided discontinuous urban fabric, industrial areas, forests, and agriculture within urban areas (Fig. 4D, Table 4).

Fig. 3
The predicted probability of urban occurrence by A European hares in relation to mean annual precipitation and B elevation, and C by mountain hares in relation to the temperature of the coldest quar-ter, and D by rabbits in relation to green urban areas (this effect was uninformative, i.e., 95% confidence intervals overlapped zero). 95% confidence intervals are shown as shading

Diurnal observation patterns and overlap between European hare and rabbit observations
The proportion of observations by land cover differed depending on the time of the day, with more observations recorded in industrial and green urban areas during the day than at night, dusk, and dawn (Table 5). There was considerable temporal overlap between urban European hare and rabbit observations (Fig. 5). When the two species were sympatric, temporal overlap of observations across all land covers was 0.88 (95% confidence interval (CI): 0.86-0.91), with both species mostly observed during dawn and dusk. Overlap was greater in discontinuous urban fabric (mean: 0.89, 95% CI: 0.82-0.89) compared to green urban areas (mean: 0.81, 95% CI: 0.77-0.86) and industrial areas (mean: 0.81, 95% CI: 0.77-0.92; Fig. 5A-C). Compared to observations in discontinuous urban fabric, European hare observations were shifted toward midday in green urban areas, as were rabbit observations in industrial areas (Fig. 5B-C). In urban areas where the two species were allopatric, the observation density of European hares was narrower, showing a clear peak in the morning, and rabbit observations showed a different temporal distribution, being more spread throughout daytime (Fig. 5D). Observation overlap (estimated across all land covers) was 0.74 (95% CI: 0.68-0.78).

Discussion
Citizen observations were useful in describing urban occurrence and habitat selection of the three lagomorphs in Sweden. The data suggest that European hares and rabbits are successful urban colonizers, and mountain hares are beginning to establish populations in some urban areas in the northern part of Sweden. Urban occurrence by all species was generally better explained by climatic conditions, elevation, and urban area size, rather than by the proportion of land cover types within urban areas or the presence of other lagomorph species. Thus, urban colonization was likely driven by suitable conditions within the distribution of each species. In contrast to our prediction, the JSDM and habitat selection analyses indicated no direct competition among the three species but indicated a positive relationship between European hares and rabbits.

Trends in urban observations
Both citizen observations and hunting bag data suggest that European hares were the most abundant of the three lagomorphs in Sweden. Hunting bag reports indicated that European hare populations are declining, a trend seen throughout Europe (Smith et al. 2005). However, this pattern must be taken cautiously, because hunting bags could have declined for other reasons (Imperio et al. 2010), such as declining hunting effort on hares. Moreover, Jansson and Pehrson (2007) reported a northward range expansion of European hares in Sweden, likely due to milder winters. The relatively stable proportion of urban European hare observations over time suggests that European hare populations have established in urban areas of Sweden, similar to urban areas in Denmark (Mayer and Sunde 2020). Like European hares, rabbits appeared to be strong urban colonizers, with nearly 40% of all observations coming from urban areas, consistent with previous findings showing that rabbits are successful urban colonizers (Ziege et al. 2016(Ziege et al. , 2015(Ziege et al. , 2020. Assuming hunting bag data to be a measure of population trends, rabbit populations fluctuated over the years. Hunting bag reports and the proportion of urban observations mirrored each other well, i.e., increases in hunting bag were accompanied by decreases in the proportion of urban observations. A potential explanation might be that rabbits were culled in urban areas to prevent damage to city parks. For example, in Stockholm 6000 rabbits were culled in 2008 and 3000 in 2009 (https:// abcne ws. go. com/ Inter natio nal/ rabbi ts-burned-fuel-sweden/ story? id= 88245 40), which coincided with the decrease in proportion of urban rabbit observations. Additionally, fluctuations in rabbit numbers could be related to fluctuations in climate conditions and/or disease outbreaks (Calvete et al. 2002;Rödel and Dekker 2012). Based on hunting bag data, mountain hare numbers were intermediate compared to European hares and rabbits, and were declining, consistent with the species' red list status in Sweden (Artdatabanken 2020). This decline has been attributed to climate warming and competitive exclusion by and hybridization with the European hare (Thulin 2003). The proportion of urban mountain hare observations increased in recent years, indicating that urban areas are increasingly colonized by mountain hares. However, this increase might also be partly related to an increased proportion of humans living in urban areas. Overall, urban mountain hare observations were less common compared to the other two species, consistent with findings showing that mountain hares select for areas of low human influence (Leach et al. 2016). Conversely, the proportion of urban observations for both European hares and rabbits might be biased in relation to Fig. 4 A The proportion of urban observations in the different land cover categories separately for the three species. For European hares and European rabbits, observations are further separated by the presence or absence of European rabbits/European hares (mountain hares were only observed in urban areas without the other two species).
Moreover, the relative probability of use by European hares (B), mountain hares (C), and European rabbits (D). For European hares, European rabbit presence affected habitat selection, but not for mountain hares. Values > 0.2 indicate selection, whereas values < 0.2 indicate avoidance. The 95% confidence intervals are given as bars mountain hare observations, because their range covered the more densely populated south of the country, potentially leading to a comparatively greater sampling effort inside urban areas (Geldmann et al. 2016).

Biotic interactions and environmental filtering
European hares and rabbits shared environmental responses, while mountain hares had distinct environmental responses, consistent with previous findings (Leach et al. 2017). This was likely related to the distribution of the three species, with the European hares and rabbits' southern distribution characterized by higher temperatures and lower elevations compared to northern Sweden, where only mountain hares occurred. Both European hares and rabbits are generally associated with comparatively warm and dry climate, and lowland areas (Calvete et al. 2004;Leach et al. 2016;Smith et al. 2005;Tapia et al. 2014), whereas mountain hares typically occupy colder areas at higher elevations (Jansson and Pehrson 2007;Thulin 2003). For all species pairs, environmental correlations were stronger on a grid cell level, probably because this finer spatial scale captured more detailed environmental differences. The positive residual correlations (both on urban area and grid cell level) between European hares and rabbits suggest that the two species co-occurred more than expected from their shared environmental responses, consistent with a previous study (Leach et al. 2017). Co-existence between the two species has been proposed to be mediated by the larger home range of the European hare, enabling localscale avoidance and diet partitioning (Lush et al. 2017;Stott 2003). Although there is evidence that European hares and rabbits are not in competition (Katona et al. 2004;Stott 2003), the study by Leach et al. (2017) and this study, to our knowledge, are the only implying a potential facilitative interaction between European hares and rabbits (also see discussion regarding habitat selection). However, positive residual correlations might represent unmodelled shared environmental preferences from variables not included in Table 4 Estimate, standard error (SE), lower 95% confidence interval (LCI) and upper 95% confidence interval (UCI) of explanatory variables for the analyses of urban habitat selection separately for European hares, mountain hares and rabbits, using generalized linear mixed models The land cover 'discontinuous urban fabric' was used as reference category, with positive estimates indicating a higher relative probability of use (selection) and negative values indicating a lower relative probability of use (avoidance) in comparison to this land cover. Sample sizes (number of observations) are given in parenthesis. Informative parameters are in bold the models, or they represent biases in citizen observations (also see discussion of habitat selection below).

Species occurrence, relative abundance, and habitat selection
Urban area size was the most important factor explaining the occurrence of all three lagomorphs. This might indicate that urban areas must be large to allow enough individuals to adjust (either via selection of bold individuals or behavioral adaptations) to the novel conditions (e.g., high level of human disturbance), and consequently establish a population. Alternatively, there might not be sufficient observers in smaller urban areas to reliably detect the presence of a species, cautioning against interpreting this finding too much in the absence of a true measure of observation effort (Kelling et al. 2015). Apart from urban area size, the probability of urban European hare occurrence decreased with higher precipitation (which included both rain and snow) and tended to increase with higher temperatures (95% confidence intervals slightly overlapped zero), suggesting that warmer and drier areas generally favor European hare occurrence (Leach et al. 2016;Smith et al. 2005), which in parts might be responsible for recent range expansions (Schai-Braun et al. 2021). Moreover, the probability of European hare occurrence increased with elevation; a counterintuitive finding, as this species is typically associated with lowland. However, the average elevation of urban areas within the European hares' distribution was 62 m, and only a single urban area was located > 210 m asl (at 300 m), i.e., all urban areas were located at comparatively low elevations. The probability of mountain hare occurrence markedly decreased when temperatures were higher, in line with this species' preference for colder climates (Jansson and Pehrson 2007). Rabbit occurrence, apart from urban area size, tended to increase when more green urban areas were present, suggesting that parks and other green areas constitute important habitat (providing forage and areas for denning) for this species, also reported in other (peri)urban areas (Sogliani et al. 2021). The general lack of urban land cover in the best models explaining the probability of urban occurrence suggests that factors explaining the general distribution of the species (climate and elevation) are better at predicting urban occurrence, especially for European and mountain hares. We have no evidence that species competition affected urban occurrence by any of the three species. The analyses of the number of citizen observations per urban area yielded different results compared to the urban occurrence and habitat selection analyses. For example, the number of mountain hare observations decreased with the proportion of discontinuous urban fabric, whereas the habitat selection analysis indicated that mountain hares selected for this land cover type. Similar contrasting results were found for European hares in relation to forest and for rabbits  (47) concerning discontinuous urban fabric. We deem the analyses of relative abundance less reliable, because the number of observations was likely more biased (based on observer distribution) compared to a presence/absence measure and compared to accounting for availability in the habitat selection analysis, though the latter might have also resulted in biases due to creating random positions in areas where no observers went. This highlights that using different analytical approaches can be useful to test the generality of findings, especially when using heterogeneous citizen science data. Inside urban areas, European hares selected for green areas (parks, sport facilities, cemeteries, etc.) in the presence of rabbits, but avoided them when rabbits were absent. General selection of green urban areas is consistent with previous findings of urban habitat selection by European hares in Denmark (Mayer and Sunde 2020), likely because these areas resemble the hares' preferred habitat, characterized by low vegetation height, providing high-quality forage (Lush et al. 2017;Mayer et al. 2018). Similarly, hares selected for discontinuous urban fabric (often consisting of residential areas) in the absence of rabbits but avoided them when rabbits co-occurred. Residential gardens, which have been found to constitute important habitats for other urban wildlife (Van Helden et al. 2020), might also constitute foraging sites for European hares. It is harder to explain the difference in habitat selection depending on the presence of rabbits that seemingly facilitated the use of green urban areas by European hares (also selected for by rabbits) at the expense of discontinuous urban fabric. One explanation could be that the presence of rabbits increased overall grazing intensity and fertilization via defecation on lawns, leading to increased grass growth, benefitting European hares. This potential facilitation of European hares by rabbits might be mitigated by dietary differences between the two species (Lush et al. 2017). Similarly, it has been shown that megaherbivore trampling and feeding stimulates highquality grass regrowth, making it more accessible for smaller ungulates (Wegge et al. 2006).
We found no evidence that the presence of European hares affected habitat selection by rabbits, indicating that rabbit space use and occurrence was unaffected by hares, as suggested in previous studies (Flux 2008;Katona et al. 2004;Leach et al. 2017;Stott 2003). Rabbits generally selected for green urban areas that likely provided good forage opportunities (Bakker et al. 2005). They showed no selection or avoidance for continuous urban fabric, and avoided the other land cover types, including forest. An avoidance of areas that likely provided cover (such as forest and discontinuous urban fabric via hedgerows) might indicate that urban rabbits experienced relaxed predation pressure, as previously proposed, reducing the need for cover (Ziege et al. 2016), in combination with these areas probably providing less forage (Lombardi et al. 2003). However, as most observations likely came from active rabbits, our results might not apply to inactive rabbits that might select for areas with more cover, leading to a reduced detection probability (Geldmann et al. 2016; also see discussion of study limitations).
The finding that the proportion of observations differed between day and night depending on land cover likely represents an observer bias. People are usually at home between dusk and dawn, potentially leading to increased observations in residential areas during this time and fewer observations in industrial areas and (to a lesser degree) parks. However, we assume that this bias varied little among urban areas, allowing us to compare temporal observation overlaps. Our findings indicate that overlap in activity times between European hares and rabbits was mediated by land cover type and sympatry. In sympatry, European hare and rabbit observations overlapped more in structured habitats (discontinuous urban fabric) compared to more open habitats (green urban areas and industrial areas). This suggests that the two species temporally avoid each other more in open areas, potentially because these more homogenous areas offer greater potential for exploitative competition due to a smaller availability of foraging species. More generally, resource availability can be a driver of inter-specific competition, thereby affecting spatio-temporal overlap between species (Karanth et al. 2017). In allopatry, the observation density distribution of European hares was narrower than in sympatry, suggesting that co-existence in areas where the two species co-occur was partly mediated by temporal avoidance of rabbits by European hares. Conversely, it appears that rabbits adjusted their temporal activity to resemble that of hares in sympatry, suggesting interactions between the two species. These findings must be taken cautiously due to the comparatively small sample size of rabbit observations, and because these patterns might have been caused by alternative factors not measured, such as predator densities and human disturbance (Moll et al. 2018).
As urban mountain hare observations were located outside the current distribution of the other two lagomorphs, we could not investigate habitat selection depending on species co-occurrence. Mountain hares selected for discontinuous urban fabric, potentially providing both forage and cover, and avoided green urban areas, industry, and forest. The apparent avoidance of forest might be related to observer biases (see below). The avoidance of green urban areas might be related to the absence of cover, as mountain hares are typically associated with habitats providing some cover, typically tundra and open forest (Flux and Angermann 1990;Thulin 2003).

Study limitations, future considerations, and conclusions
Citizen science data are susceptible to spatial biases with regards to infrastructure and human population density (Geldmann et al. 2016). Consequently, citizen observations might have measured human-lagomorph encounters rather than actual habitat preferences, e.g., shown for canids (Mueller et al. 2019). Urban areas, while generally having high levels of infrastructure and human population densities, yielding a high sampling effort overall, might still be prone to varying sampling efforts due to being highly heterogeneous (Crall et al. 2011;Dickinson et al. 2010). For example, it is plausible that citizens rather recorded animal observations in their own gardens and in parks compared to city centers and industrial areas. Moreover, detectability also differs between land cover types, accessibility, and depending on group size and animal activity (Mair and Ruete 2016;Pereira-Ribeiro et al. 2019). As most observations likely came from active lagomorphs, our results probably represent occurrence and habitat selection of active individuals and from areas that were easily accessible to observers. However, habitat selection by active and inactive lagomorphs differs (Mayer et al. 2018;Neumann et al. 2012), implying that we might have underestimated the importance of certain land cover types that are predominantly used by resting individuals (e.g., forest patches). Avoiding such biases in citizen observations will be hard. One potential solution would be to select larger spatial scales, as scaling up generally decreases spatial bias and reduces pseudo-absences (Rondinini et al. 2006), and to define species occurrence rather than relative abundance. Finally, species might have been misclassified in some cases, resulting in false-positives (Dickinson et al. 2010). GPS tagging individuals would enable us to obtain more detailed information on habitat selection and movements by lagomorphs in urban areas, shedding more light on their adaptations to this novel environment. To quantify urban population densities, transect counts could be used (Mayer and Sunde 2020), potentially conducted by citizen scientists if incentivized correctly, like for example the Great Backyard Bird Count (https:// www. birdc ount. org/).
Our study contributes to the understanding of species co-occurrence patterns and habitat preferences within urban areas, while highlighting the benefits and challenges of citizen science data. We generally found little evidence for competition between the three lagomorphs, though we cannot exclude that urban mountain hare occurrence is inhibited by interspecific competition. Future studies should also investigate how the presence of predators, in this case predominantly red foxes (Vulpes vulpes), affects the occurrence and habitat selection of lagomorphs within urban areas. Moreover, it would be of interest to shed more light on the drivers of urban colonization by wildlife, to be able to predict urban species occurrence. Insights into species habitat associations within urban areas and depending on co-occurrence with other species can help in targeting urban management plans, which will be useful to identify suitable habitats for desired species and efficient management of conflict species (Apfelbeck et al. 2020;Gaertner et al. 2017).