A shifting carnivore ’ s community: habitat modeling suggests increased overlap between the golden jackal and the Eurasian lynx in Europe

,


Introduction
Despite the increasingly severe effects of climate change, land use alterations, and habitat degradation, biodiversity is not spatially declining evenly: there are several cases of species that expanded their ranges (Nolet and Rosell, 1998;VerCauteren, 2003;Arnold et al., 2012;Deinet et al., 2013;Chapron et al., 2014;Halliez et al., 2015;Carvalho et al., 2018;Pacifici et al., 2020;Wehr, 2021) or managed to survive by adapting to new environmental conditions, which also involve the human presence (Basille et al., 2009;Chapron et al., 2014;Bouyer et al., 2015;Razěn et al., 2016;Fenton et al., 2021).Of course, these processes did not occur alone, but together with successful conservation actions which allowed many species to recover or survive (Hoffmann et al., 2010;Deinet et al., 2013).This is the case of Europe, where the Living Planet Index, reports a slow decline rate linked both to conservation successes (Hoffmann et al., 2010;Chapron et al., 2014) and to the natural adaptability and recolonization capacity of some species (Ledger et al., 2022).Also, a meta-analysis noted increases in species richness and abundance in many bioregions of Northern and Eastern Europe (Pilotto et al., 2020).
Although Europe is one of the least wild continents, considering the human footprint which led (and continues to lead) to the loss of ecosystems (Williams et al., 2020) and the declining trend of some species (Inger et al., 2015;Warren et al., 2021), many species appear to be able to adapt, or in general, to persist.We are witnessing the comeback of various species of mammals, including the European bison (Bison bonasus), the elk (Alces alces), the Eurasian beaver (Castor fiber), the southern chamois (Rupicapra pyrenaica), and the alpine ibex (Capra ibex) (Deinet et al., 2013;Ledger et al., 2022).Also, the main European predators, such as the brown bear (Ursus arctos), the grey wolf (Canis lupus), and Eurasian lynx (Lynx lynx) are expanding their ranges (Chapron et al., 2014).
From a conservation perspective, the recovery of these species and their range expansion is particularly relevant, as predators require large spaces to survive and can easily come into conflict with local communities.It should not be forgotten that their return is linked to a series of factors such as an increase in forest cover, a decrease in human population densities, protective legislation, acceptance by public opinion, and the use of measures that allow coexistence with humans (Chapron et al., 2014;Cimatti et al., 2021).Preserving this delicate equilibrium between apex predators and human activities is thus important since not far as a hundred years ago these have been hunted and persecuted, and their decrease has facilitated the expansion of mesopredators (Prugh et al., 2009;Ripple et al., 2013;Newsome et al., 2017).
In Europe, the golden jackal (Canis aureus), represents an excellent example of this phenomenon; it is currently an expanding species, whose growth was initially favored by the absence of the wolf (Krofel et al., 2017).The European distribution of this medium-sized canid species was once limited to the Balkan area (Krysťufek et al., 1997), and saw a decline until 1950 caused by human persecution (Spassov, 1989) which reduced it to a few populations among Former Yugoslavia, Bulgaria, Turkey, and Greece (Spassov, 1989;Giannatos, 2004;Spassov, 2007).Then, two waves of expansion, the first occurred around the 1950s, and the second from the 1980s, also favored by the legal protection (in Bulgaria from 1962) led to a north-and westward expansion (Trouwborst et al., 2015;Krofel et al., 2017); in the following years, the jackal colonized most of the countries of central and southern Europe.More recently, a small population has established in the Baltic countries (Pysǩováet al., 2016).Dispersing individuals have also been observed in northern Europe (Sørensen and Lindsø, 2021;Rykov et al., 2022), and central Europe (Hatlauf et al., 2016;Böcker et al., 2023), and its presence has been reported in 33 European countries, 10 of them in the last decade (Hatlauf et al., 2021).
The golden jackal is a generalist species with opportunistic feeding habits, with its prey spanning from ungulates to rodents (the most common prey), also comprising invertebrates, vegetables, and anthropogenic resources, including waste and carcasses (Markov and Lanszki, 2012;Aleksandra and Dusǩo, 2015;C ́irovicé t al., 2016;Lange et al., 2021).The ecological plasticity of the golden jackal also emerges in the choice of habitat: the most used one is a mixture of shrub and herbaceous vegetation, including mesophilic woods or riparian habitats (Lapini et al., 2011;S ̌aĺek et al., 2014;Spassov and Acosta-Pankov, 2019), but sometimes it may also be attracted by human settlements (Lapini et al., 2011).However, the habitat choice of the golden jackal is highly influenced by the wolf presence (Krofel et al., 2017), and climatic-based models predict its change for the next future (Cunze & Klimpel, 2022).
The changing distribution of a species, such as in the case of the golden jackal which is currently expanding its range, can lead to new biotic communities, and new interactions with other species that find themselves coexisting for the first time (Gilman et al., 2010;Pecl et al., 2017).In the case of the golden jackal, it is exceptional the recent observation of kleptoparasitism towards the Eurasian lynx (Krofel et al., 2022), another predator expanding in Europe.
The Eurasian lynx is a European top predator that was intentionally eradicated from many parts of Europe between the 19th century and the beginning of the 20th century (Breitenmoser, 1998).Only four populations survived, in Scandinavia, the Carpathians, the Baltic area, and in the southwestern Balkans (von Arx et al., 2004).The species recolonized parts of its former range thanks to several factors, from legal protection and reintroductions to an increase in prey availability and greater public acceptance (Linnell et al., 2009), and is now slowly recovering (9-10,000 individuals estimated in Europe) (Chapron et al., 2014).This species is a predator specialized on roe deer (Capreolus capreolus) (Jędrzejewski et al., 1993), but the diet varies regionally and may include other species, from hare (Lepus europaeus) and chamois (Rupicapra spp.) to red deer (Cervus elaphus), also including livestock (Okarma et al., 1997;Stahl et al., 2001;Odden et al., 2006).The Eurasian lynx is associated with forest habitats with good environmental heterogeneity for hunting (Belotti et al., 2013) and many hiding places to rest (Podgoŕski et al., 2008).However, the species can select also more open habitats, especially at night (Filla et al., 2017) adapting to them by generations (Nagl et al., 2022).Also, it can adapt to habitats with moderate human activity if prey abundance is elevated and refuges are available, seeking a compromise between human avoidance and prey hunting (Basille et al., 2009;Bouyer et al., 2015;Oeser et al., 2023).The overlap between the golden jackal and Eurasian lynx ranges was scanty until the mid-20th.Indeed, the recent expansion of the jackal, especially since 2011, caused an increase in this overlap, currently estimated at 13% (Krofel et al., 2022).Considering the expansion of the golden jackal and the recovery of the Eurasian lynx, this overlap is likely to increase, even though no information is available to date.Moreover, considering the positive trends of the wolf populations, and their influence on the golden jackal (Krofel et al., 2017;Newsome et al., 2017), it could be also important to assess whether the range changes of both golden jackal and Eurasian lynx would increase in areas where the wolves are well-established, or where is predicted an expansion in the future.
To lay the foundation for the study of this ecological issue, we infer the areas of potential current and future overlap between the two target species, as they may indicate the zones where the interactions are more likely to occur.To achieve this result, we take advantage of the most up-to-date Ecological Niche Modeling techniques and post-modeling GIS spatial analyses.The use of Ecological Niche Models for the management of expanding species and the understanding of their dynamics has increased in recent years (De Marco et al., 2008;Jimeńez-Valverde et al., 2011;Maiorano et al., 2011;Barbet-Massin et al., 2012a;Barbet-Massin et al., 2012b;Dawe and Boutin, 2016;Iannella et al., 2019;Rodrıǵuez-Rey et al., 2019;Serva et al., 2023).Also, Ecological Niche Models can be refined through a post-modeling phase, converting them into more precise Species Distribution Models (SDMs), in a process which increases the predictive performances and allows the creation of more realistic models (Iannella et al., 2021).
Given the ongoing range changes observed for both the Eurasian lynx and the golden jackal, and the recently observed interactions between them, we build SDMs for these two species under current and future conditions.Then, we use those SDMs to analyze their possible range overlap under both current and future scenarios, using for these latter the areas predicted to remain stable and the ones predicted to be gained.We finally use these zones to highlight the areas where possible interaction may occur.Moreover, given the great influence of the wolf on the golden jackal, we build environmental niche models (ENMs) to estimate the current and the future predicted suitable areas for the wolf, to understand if the predicted overlap areas of the Eurasian lynx and golden jackal would overlay with this other large European carnivore.

Target species and study area
The target species are represented by the Eurasian lynx, a top predator, and the golden jackal, a mesopredator.Secondly, we focused on the grey wolf, another European apex predator.
The Eurasian lynx is classified as Least Concern by the IUCN and is included both in the Bern Convention (Appendix III except for the Balkan lynx balcanicus, listed in Appendix II) and in the EU Habitats Directive (Annexes II and IV, except for some countries where it is included in Annex V and excluded from Annex II).The golden jackal is also classified as Least Concern by the IUCN, and it is included in Annex V of the Habitats Directive.However, the protection of the golden jackal in Europe has been under discussion, since some countries reported it as an alien species, and the Habitats Directive imposes distinct limitations on national policies and national management options (Trouwborst et al., 2015).Indeed, currently, both populations of the target species are increasing (Arnold et al., 2012;Deinet et al., 2013;Chapron et al., 2014;Spassov and Acosta-Pankov, 2019).Similarly, the grey wolf is increasing in Europe, and is classified as Least Concern by the IUCN.Like the Eurasian lynx, the grey wolf is protected both by the Bern Convention (Appendix II) and the EU Habitats Directive (Annex II and IV), even if there are some local derogations.
The area selected for our study is the western Palearctic, a region where the golden jackal is expanding and co-occurring with the lynx.

Environmental predictors
To encompass as much as possible the ecological needs and the factors influencing the two target species' distribution, we used climatic, topographic, habitat, hydrographic, human disturbance and biotic predictors, each of which is described below.
The climatic aspect was assessed through the 19 bioclimatic variables available on WorldClim (version 2.1, 2.5 arc-min, Supplementary Document 1) (Fick and Hijmans, 2017) downloaded for current and future (2050 and 2070) conditions.For each future scenario, we selected the shared socio-economic pathways (SSP) 2.45, 3.70, and 5.85 to involve all the different possible trajectories except the most optimistic one (SSP 1.26) (Riahi et al., 2017).Since the use of a single General Circulation Model (GCM) in future projections is shown to produce variable predictions (Stralberg et al., 2015), we selected and managed three different GCMs, namely the BCC-CSM2-MR (Wu et al., 2019), the IPSL-CM6A-LR (Boucher et al., 2020) and the MIROC6 (Tatebe et al., 2019) for all future models' forecasts.
We collected elevation data from the ASTER NASA database (https://asterweb.jpl.nasa.gov/gdem.asp),deriving the slope raster map at the same spatial resolution (2.5 arc-min) ('Slope' function, Spatial Analyst in ArcGIS Pro), because the slope is considered as a relevant variable for the ecology of both the Eurasian lynx (acting as a proxy for resting sites) and the golden jackal (Filla et al., 2017;Spassov and Acosta-Pankov, 2019;Hocěvar et al., 2021).
About the habitat predictors, we took advantage of the 100-m resolution Global Corine Land Cover map of the Copernicus repository (https://land.copernicus.eu/global/products/lc),containing a discrete classification with 23 classes according to the UN-FAO Land Cover Classification System (Buchhorn et al., 2020).
To incorporate the human influence over the two target species, we considered the influence of roads, since they represent a major limit to the dispersal of lynx and jackal, and are capable to alter movement behavior, increase direct mortality and lower patches' connectivity (Kramer-Schadt et al., 2004;Lapini et al., 2011;Heurich et al., 2018).For these reasons, we downloaded the roads data from the "World Roads" Layer Package managed by ESRI (ESRI Inc, 2022) and divided them into major roads and highways.To use this information as a modeling predictor, we calculated the distances from those, classifying them into four intervals (0-500 m, 500-1,000 m, 1,000-5,000 m, 5,000-20,000 m).To further account for human disturbance, especially important for the Eurasian lynx (Ripari et al., 2022), we also downloaded the GHSL SMOD data (Global Human Settlement Layer, available at https:// ghsl.jrc.ec.europa.eu/) at ~1 km resolution: this map takes into account both population densities and built-up area densities, in a unique spatial data (Pesaresi et al., 2019).
Also, we used three other predictors for golden jackal modeling only: the average snow cover map (since snow cover was used in previous studies and could partially limit its suitable habitat) (Ranc et al., 2015;Spassov and Acosta-Pankov, 2019;Männil and Ranc, 2022), the river network (as riverine habitats are known to be important dispersal corridors; Arnold et al., 2012;Spassov and Acosta-Pankov, 2019) and the grey wolf current range (since this species is known to be a major limiting factor for the jackal, especially through intraguild predation; Krofel et al., 2017;Ranc et al., 2018;Spassov and Acosta-Pankov, 2019).
Specifically, for the snow cover information we used the MODIS data (Terra Snow Cover Monthly L3 Global, ~5 km resolution, available at https://nsidc.org/data/MOD10CM; Hall and Riggs, 2021) averaging each month from 2010 to 2020.
The rivers' network information were first obtained from the HydroRIVERS dataset (https://www.hydrosheds.org/products/hydrorivers; Lehner and Grill, 2013) and then processed through the 'Euclidean distance' tool in ArcGIS Pro 10.3 (ESRI Inc, 2022), to obtain distance intervals form each watercourse, with the following intervals: 0-1 km, 1-2 km, 2-3 km, 3-10 km, > 10 km.Finally, the grey wolf current range was obtained at a resolution of 10 km from the work of Kaczensky et al. (2021), giving a different weight to the cells in which the presence of the wolf is permanent, compared to the cells in which the presence of the wolf is sporadic.

Ecological niche modeling
The occurrence datasets were obtained from the GBIF repositories for both golden jackal (GBIF, 2022a) and Eurasian lynx (GBIF, 2022b).Two filters were applied, namely a temporal (selecting the occurrences 1970-current only) and a spatial (coordinate uncertainty < 5 km) one.Moreover, we added occurrence data from recent reports in literature of both the golden jackal and the Eurasian lynx (Supplementary Tables), comparing the final datasets with the known distributions of the species.For the purposes of the modeling process, we processed the GBIF-filtered records so as to obtain a comparable spatial resolution between the occurrences and environmental predictors (Sillero and Barbosa, 2021), moreover reducing the spatial autocorrelation.In order to do that, we used the 'Spatial rarefy' tool of SDMtoolbox (Brown et al., 2017).
To better explore the potential expansion of the golden jackal in the future, we decided to keep the records from dispersing individuals (~19% of the total occurrences), given that these individuals are likely to form stable population, considering his ongoing expansion pattern.The records belonging to dispersing individuals can be noted in Supplementary Figure 1, together with the core areas for each population.However, to assess possible biases caused by the inclusion of dispersing individuals in the species distribution models or related to the nature of the dataset (mainly composed of citizen-science-derived data), we repeated the analysis (building ENMs) using only records from the dataset of Ranc et al. (2022) selecting only the cells where golden jackal is reported as "permanent" or "high confidence presence" (expertbased ENMs).
To model the grey wolf predicted distribution for the current and future projections, we relied on to the dataset of Kaczensky et al. (2021), considering only the cells where the species is reported as "permanent".Then, we assessed the multicollinearity of the 19 bioclimatic variables by measuring the Variance Inflation Factor ('vifstep' algorithm of the 'usdm' R package (Naimi, 2015), threshold = 10); we included the ones both not exceeding this threshold and being important for the species' ecology (Brandt et al., 2017).To minimize the variation caused by the use of individual General Circulation Models (GCMs) in future projections (Stralberg et al., 2015), we used during the modeling process the three different GCMs, namely the BCC-CSM2-MR (Wu et al., 2019), the IPSL-CM6A-LR (Boucher et al., 2020) and the MIROC6 (Tatebe et al., 2019), and the three SSPs mentioned before for future projections (2.45, 3.70, and 5.85).
To build the ENMs we used the 'biomod2' package (Thuiller et al., 2016) in R (R Core Team, 2016).The advantage of this package is that allows to create ensemble models, which are the result of the combination of individual algorithms' outcomes, corresponding to different modeling techniques.We generated five sets of 1,000 pseudoabsences each (for the wolf, five sets of 10,000 pseudo-absences) using the Surface Range Envelope algorithm (quantile = 0.05) (Barbet-Massin et al., 2012a).We parametrized the single models through the 'BIOMOD_ModelingOptions' as follows: GLM: type = "quadratic", interaction.level= 3; GBM: n.trees = 10,000, interaction.depth= 3, cv.folds = 10; GAM: type="s_smoother", interaction.level=3;RF: n.trees=10,000, bag.fraction = 0.75, shrinkage = 0.001.Then, the 'BIOMOD_Modeling' function was used to calibrate the single models, and the 'BIOMOD_EnsembleModeling' algorithm was finally run to obtain the Ensemble models, using as input only single models with high performances (AUC >0.8 and TSS >0.7).To merge the single models, we use the 'wmean' algorithm, that averages models into the Ensemble one based on a weighted mean function, in which weights are the discrimination scores (Thuiller et al., 2016).Finally, we projected to future scenarios with the 'BIOMOD_EnsembleForecasting' function.
For each projection of the same future scenario (i.e., keeping as fixed the year/SSP combination, for instance the 2050 SSP 3.70), we reduced the uncertainties deriving from the calibration performed on each GCM as follows.We first checked for model extrapolation (i.e., the dissimilarity from the calibration conditions) by assessing the Multivariate Environmental Surface Similarity (MESS) (Elith et al., 2010), computed through the function 'mess' of the 'dismo' package (Hijmans et al., 2016).The MESS maps so obtained were then included in the Multivariate Environmental Dissimilarity Index (MEDI) algorithm (Iannella et al., 2017), a form of weighted raster average which down-weights extrapolation based on projections and MESS, thus reducing the uncertainties in future projections.We finally used those MEDI-corrected future projections as baseline models for all the following analyses.

Weighted overlay process and postmodeling analyses
To obtain a more informative and robust models, we applied the 'Weighted overlay' framework of Iannella et al. (2021) to both target species.This procedure permits to refine Ecological Niche Models into more precise Species Distribution Models, incorporating, in a post-modeling GIS environment, the information not included during the Ecological Niche Modeling process, also increasing predictive performances (Iannella et al., 2021).Thus, we included in this process all the predictors mentioned above except the climatic ones, which were included in the ENM calibration and projection.The information about slope, habitat, population density, snow cover, and wolf range was extracted in occurrence localities for both species through the 'Extract multi values to points' tool in ArcGIS.Then, relative frequencies were calculated for each predictor/occurrence localities and converted into absolute frequencies, on a 1-to-10 scale, to make them comparable with the ENMs.The only exception we made is the conversion into the 1-to-10 scale of the roads and rivers distance categories, deduced from literature (roads information for Eurasian lynx model referring to Dobbert et al. (2021); rivers data for the golden jackal modeling referring to Spassov and Acosta-Pankov (2019) and thus not obtained through the occurrences' frequencies).
Specifically, we used as input into the weighted process on current conditions the ENMs, habitat, population/built areas density, slope, distance from roads and rivers, snow cover, and wolf range predictors to model the distribution of the golden jackal; about the Eurasian lynx, we used the ENMs, habitat, population/ built areas density, slope, and distance from roads ones.
For the projections to future scenarios, we kept all the predictors as fixed, except the ENMs; when needed, weighted models were binarized using as a threshold the 10 th percentile value (Iannella et al., 2021).Subsequently, the 'BIOMOD_RangeSize' function of the 'biomod2' package (Thuiller et al., 2016) was applied to calculate the areas predicted to be lost, remaining stable, or gained from the two target species in the future scenarios.Considering the presenceonly nature of the data, we assess the predictive performance of the weighted models (as well as for the ENMs presented above) through the Boyce index (Hirzel et al., 2006;Leroy et al., 2018).
To evaluate areas where sympatry between the two species may occur, we assessed the overlap of their respective "gain + stable" areas for both the present and the future (considering all the SSPs for this latter case).
Wolf ENMs were binarized through the 10 th percentile; we overlap the predicted suitable areas with the zones where future sympatry between golden jackal and Eurasian lynx is predicted to occur.

Ecological Niche Models
After the presence localities filtering processes, we obtained 2,164 occurrences for the Eurasian lynx, and 564 for the golden jackal within the study area (Figures 1A, B).These datasets are roughly comparable to current distributions of the target species, even if they show some gaps, especially for the golden jackal in the core distribution areas.
Considering the results of the VIF analysis and the ecological characteristics of the species, we kept ten bioclimatic variables for the modeling of the The ENMs for the Eurasian lynx, the golden jackal, and the wolf resulted in high predictive performances, scoring a Boyce index value of B = 0.897, B = 0.869, and B = 0.921 respectively.The models calibrated on current climatic conditions reveal a wide, highly suitable territory in central Europe, between Germany and Poland, in the Pannonian basin, and the Po valley in Italy, for the golden jackal (Supplementary Figure 2A).Also, some minor and highly suitable patches appear in southern France and eastern Iberian peninsula, in the Northwestern coasts of the Black Sea, and Caucasus (Supplementary Figure 2A).Regarding the future predictions, a great expansion of the climatically-suitable areas is inferred (Supplementary Figures 3A-F).In general, a northeastern suitability gain is observed, as well as an increase in Germany and France, for all the future scenarios (Supplementary Figures 3A-F).The expert-based ENM for the current condition concerning the golden jackal shows similar habitat suitability compared with the prediction obtained with the first ENMs, exception made for a minor suitability between Poland and Germany, and southern Sweden (Supplementary Figure 4).The future habitat suitability maps obtained from expert-based ENMs indicate an analogous pattern of increase in habitat suitability towards the northernmost part of the study area (Supplementary Figures 5A-F).However, a slight difference is observed in Finland and among Ukraine and Russia (Supplementary Figures 3G-L, 5A-F).Regarding the Eurasian lynx, two highly suitable areas emerge, namely the main European mountain systems (Pyrenees, Alps, Apennines, Dinarides, Carpathians, and Caucasus) and the central and northern European territories, with the latter being vaster than the former (Supplementary Figure 2B).In future scenarios, the climatic suitability is generally predicted to massively decrease from all the aforementioned territories, especially in the central and northern European ones (but excluding the Scandinavian peninsula) (Supplementary Figures 3G-L).

Species' frequencies in nonclimatic predictors
The extraction of the information about slope, habitat, population density, snow cover, and wolf range performed in occurrence localities for both species highlighted different trends between the species.About the habitat information (Supplementary Figure 6A), 40% of the occurrences of the golden jackal fall into the "Cropland" habitats, with the lynx reaching only 10% of its occurrences in this category.On the other side, the "Closed forests" resulted to be the most preferred habitat for the Eurasian lynx (~52%), while the golden jackal occurred in 27% of its total occurrences.The "Open forest" category had similar percentages for the two species: 11.4% for the Eurasian lynx and 11% for the golden jackal.Finally, the "Herbaceous vegetation" category is more used by the lynx (21.4%) compared with the golden jackal (6.4%).
Regarding slope, the golden jackal highly prefers low values of it (53.9% of occurrences in the 2% slope category, followed by 17.8% in the 3% slope one), while the Eurasian lynx uses territories both with low (32.5% for the 2% category, 16.6% for the 3%) and medium (10% of slope used in the ~5% of the occurrences) slopes (Supplementary Figure 6B).
Another result from the frequency of occurrence emerges from the ones falling into territories where human presence occurs.The golden jackal is mainly found in human settlements classified as low-density rural or rural clusters (78% and 13%, respectively), while the Eurasian lynx, when occurring in those settlements, is mainly found in lowdensity rural patches (89%) (Supplementary Figure 6C).

Weighted suitability models and overlap analyses
The weighted models obtained starting from the ENMs and from the frequency of occurrence of the two target species into the non-climatic predictors are reported for current (Figures 1C, D) and future (Figures 2A-L) conditions.The predictive performance resulted in B = 0.988 for the Eurasian lynx and B = 0.943 for the golden jackal.
The current weighted map resulted for the Eurasian lynx (Figure 1D) emphasized the high-suitability areas which are compliant with the current distribution of the species (Figure 1B), except for the suitable territories in the Apennines.In contrast, the golden jackal current weighted suitability map (Figure 1C) shows high values in areas from which the jackal is expanding, such as southern and eastern Europe.On the other side, areas where the Future suitability predictions (low to high on a blue-to-red scale) for (A-F) the Eurasian lynx and (G-L) the golden jackal, for each different Shared Socio-economic Pathway and projection year considered.Serva et al. 10.3389/fevo.2023.1165968Frontiers in Ecology and Evolution frontiersin.orgjackal is currently extending its range (e.g., Baltic countries, Finland, and Norway) show lower environmental suitability.About future weighted scenarios, the Eurasian lynx is predicted to face a reduction in its environmental suitability.The greatest reduction is observed for the 2070_SSP8.5scenario, where only some clusters in central and southern Europe persist with mediumhigh values of environmental suitability (Figures 2A-F).Similarly, all the other projections in other year scenarios report similar trends, indicating that all the southern populations of the Eurasian lynx (the ones occurring in the Alps, Dinarides, Balkans, and Carpathians) could face a decrease in habitat suitability in the near future.Specifically, a clear decrease in suitability is observable in all future projections, and it is highlighted the fragmentation effect of the road network, especially in the Alps and the Balkans, resulting in a patchy mosaic of suitable areas separated by these infrastructures.By contrast, Scandinavia, the Baltic area, and northern Russia's suitability remains stable at high values, even though the road network heavily influences those areas.
The future scenarios for the golden jackal, on the other hand, report an increase in its environmental suitability in the entire study area, except for the Atlantic coast of Scandinavia and northern Russia (Figures 2G-L).This can be observed across central Europe and in the Baltic States.Also, Italy, France, Balkan countries, and southern Sweden are predicted to be more suitable, with the greater increment represented by scenario 2070_8.5,where also Finland appears to be environmentally suitable for the golden jackal.
When binarizing the weighted predictions (10 th percentile threshold = 7.6 for the Eurasian lynx, Figure 3A, and threshold = 6.4 for the golden jackal, Figure 3B) and subsequently assessing the current potential overlapping territories for the two species, we find that this is mainly limited to areas where mountain systems are not present (Figure 3C).

Future sympatry areas and overlap with grey wolf distribution
In the future, an evident trend emerges: the overlap between the two species gradually shifts towards northeastern Europe, concurrently decreasing in southern Europe (Supplementary Figure 7).When comparing these trends with the range shifts analyses (Supplementary Figure 7), we find that this pattern is driven both by the expansion of the environmental suitability of the golden jackal towards the north-east, and by the reduction of the environmental suitability of the Eurasian lynx in the southwest.
In particular, sympatric areas gradually increase in Finland, southern and eastern Sweden, and Russia, especially in the scenario 2070_8.5.
A slight increase in overlap arises on the slopes of the Carpathians Mountains, in the Balkan area, and in the Dinaric range, where both species occur.Also, along the Apennines, it emerges an increase in the overlapping ranges, where however the Eurasian lynx is not present and the golden jackal sporadically occurs.
Comparing these results with the current wolf-predicted suitable areas, most of the current sympatric areas appear to be overlapped (Figure 4).These mainly occur in the Carpathians Mountains, the Dinaric range, and various zones between Bulgaria, Albania, Macedonia, and Greece (Figure 4).Moreover, the current overlapping areas from the Baltic area to those in central Europe are also suitable for the grey wolf (Figure 4), which occur in some of these countries with stable populations.Similarly, between Italy and France, the predicted sympatric areas for the golden jackal and the Eurasian lynx are eligible for the wolf, which is recently spreading in these countries.Thus, considering the current wolf suitable areas, the potential sympatric areas for the golden jackal Binarized predictions and corresponding overlap.The binarized weighted suitability model for (A) the Eurasian lynx, (B) the golden jackal, and (C) the shared areas between the two models, obtained by overlapping them.Serva et al. 10.3389/fevo.2023.1165968Frontiers in Ecology and Evolution frontiersin.organd Eurasian lynx (lacking wolf suitable areas) are only a few, located in the Pannonian Basin, Germany, Slovenia, and Russia (Figure 4).Concerning the future predicted sympatric areas between the golden jackal and the Eurasian lynx, it is shown a gradual increase towards the northernmost part of the study area (Figure 5).Sympatric areas expand across the Baltic area, Russia, and Scandinavia.On the contrary, sympatric areas remain stable or decrease elsewhere in future scenarios (Figure 5).However, considering the future predictions for the wolf distribution, the great majoring of the predicted sympatric areas in Scandinavia and the Baltic area are covered by suitable areas (Figure 5).Thus, considering the predicted expansion of the wolf in these countries, it is unlikely that the golden jackal could spread significantly.A different situation could occur in central and southern Europe.In fact, the future predictions for the wolf show that its suitable areas will increase between France and Italy, and between France and Germany, but reduce in eastern Europe, with suitable areas kept only in the Carpathians, Balkans, and Dinaric range (Figure 5).Then, sympatric areas between the golden jackal and the Eurasian lynx in these areas could expand, such as in some zones in central Europe, between Germany and Poland (Figure 5).

Discussion
Europe is seeing a great return of some wildlife species in recent years (Deinet et al., 2013;Ledger et al., 2022), despite the effect of climate and land use changes.While these may favor the expansion of some species, for others they represent a threat that makes their future uncertain (Pacifici et al., 2020;Ledger et al., 2022), as also our results suggest (Figures 2A-L).
Indeed, the expansion of the jackal is favored by several factors: deforestation and land use changes (S ̌aĺek et al., 2014), additional anthropogenic trophic resources, the absence of competitors (such as the wolf) (Pysǩováet al., 2016;Krofel et al., 2017;Spassov and Acosta-Pankov, 2019).In addition to these elements, climate change may also favor the expansion of the jackal.While in the past it was believed that deep snow could act as a limiting factor for the golden jackal, recent findings proved that the species is able to adapt even in areas with stable snow cover, such as in the boreal ecoregion.Although, in the Baltic states, it was reported that the golden jackal tends to avoid zones with permanent snow cover, preferring coastal areas without snow (Männil and Ranc, 2022).In our results, the greatest expansions are observed in the scenarios with the greatest impact of climate change, especially in the northern part of the study area (Figure 2G-L).The same expansion process affects not only higher latitudes but altitudes as well.In fact, the Range Shift analysis shows a growing rise of "gains" in the main mountain ranges of the study area.The strength of our results is highlighted by the fact that they are in line with those recently obtained (Cunze and Klimpel, 2022) for a similar area, with the important difference that in our analysis we used not only climatic predictors but also non-climatic parameters.
The redistribution of a species can lead to the appearance of new interactions with the fauna already present (Pecl et al., 2017).Changes in the dynamics of predation, herbivory, competition, and mutualism can all have substantial impacts at the community level (Gilman et al., 2010).The expansion of the golden jackal may change these dynamics, particularly interspecies competition, as Predicted current overlap between the Eurasian lynx and the golden jackal (dark grey areas) and the grey wolf predicted suitable areas as inferred from the ENMs built on current climatic conditions.observed from interactions with the red fox (Vulpes vulpes), where the golden jackal is dominant (Giannatos, 2004;Stoyanov, 2012).On the contrary, with the grey wolf, the cases of competition are few, and the golden jackal is usually defeated (Krofel et al., 2017).Competition is also possible with the European badger (Meles meles), pine marten (Martes martes), and stone marten (Martes foina) (Stratford, 2015).Data about the interactions between golden jackals and other large European predators are scant, but recently it has been observed a case of kleptoparasitism toward the Eurasian lynx (Krofel et al., 2022).
The recovery of the Eurasian lynx has been linked to a series of factors: various measures of legal protection, greater acceptance of this species by the population, reintroductions and recent population reinforcements in many European countries, the return of ungulate species that have expanded the availability of prey (Breitenmoser, 1998;Linnell et al., 2009).The modern expansion of this species is also associated with reduced human population density and increased forest cover (Cimatti et al., 2021).The current suitability map shows that there are large suitable zones in the study area where the species does not occur even if it could survive (e.g., in the Alpine range).The current environmental suitability is comparable with that obtained in a previous study (Dobbert et al., 2021).Our spatially explicit results are in line with those obtained from previous studies which explore lynx habitat selection at European scale (Ripari et al., 2022;Oeser et al., 2023).In fact, our maps of environmental suitability show that at the landscape scale, lynx tend to avoid human-modified areas and human infrastructures network, with habitat suitability remaining clustered in forested areas and/or zones associated with moderate slopes.But the future of this iconic predator is likely to be worse: our results of environmental suitability for the future show clear reductions, particularly around the Alps, in the Carpathian Mountains, and the Balkan area, with the suitability remaining high only in Scandinavia and northwestern Russia (Figures 2A-F; Supplementary Figure 4).Despite these reductions, the Eurasian lynx is a plastic species, especially regarding climatic conditions, and is able to survive even in warmer and drier habitats (Flezǎr et al., 2023).
In addition to these reductions, there are also some areas with increases in habitat suitability.In particular, beyond the Scandinavia and the Baltic area, good environmental suitability will be maintained in the Alps.Considering these predictions, it is essential to encourage conservation measures for the populations of Eastern Europe (such as the Balkan population and the Dinaric population), but it is also important to improve connectivity along the Alpine range to favor the expansion of the Alpine population.
Our overlap analysis shows that in the future the two species could have much more than the current 13% overlap between the ranges, as currently found by Krofel et al. (2022): this increase is mainly linked to the predicted expansion of the golden jackal throughout Europe, especially in the area between Russia, Baltics States and Scandinavia (Figures 2G-L; Supplementary Figure 4).In the same countries, Eurasian lynx suitability is expected to increase, so it could become an area of sympatry in the future.However, considering that these areas host stable populations of wolves, and that habitat suitability is predicted to increase in the future (Figure 5), it is unlikely that the golden jackal could expand significantly in the future, or at least his expansion will be restrained only in areas where the wolf is less abundant.
From a conservation point of view are interesting the overlap dynamics in the Carpathians, the Balkan area, and the Alpine area.Although our maps show future losses for the Eurasian lynx, some Predicted future overlap (A-F) between the Eurasian lynx and the golden jackal (dark grey) and the binarized predicted suitable range, inferred for future climatic scenarios, of the grey wolf.
suitable patches will be maintained, and for the Alps, some increases are predicted.However, "gains" are also expected in the same areas for the golden jackal, therefore it is probable that in the future the slopes of these mountain ranges may also become areas of sympatry between the two species.But in the same mountain ranges the wolf is also expanding, and habitat suitability is predicted to increase (Figure 5).Thus, even in this case, the predicted sympatric areas between the golden jackal and the Eurasian lynx could be limited by the wolf's presence.However, even considering the wolf's presence, interactions between these species can't be excluded a priori.And the consequences of this coexistence, not only between the golden jackal and Eurasian lynx but involving the grey wolf and the brown bear as well, must be deepened, especially considering that the lynx populations of Alpine, Balkan, and Dinaric areas are the focus of several conservation efforts, such as translocations and reintroductions (Toplicǎnec et al., 2022).
A previous sympatric occurrence between the golden jackal, grey wolf, Eurasian lynx, and brown bear, was reported in a dense forest in Slovakia (KysuckéBeskydy Mountains), demonstrating the capacity of the golden jackal to coexist even in the presence of other large carnivores (Guimãraes et al., 2021).Thus, the adaptability of the golden jackal could led to different interactions with these large carnivores, beyond the ones reported.
Observing future scenarios, the greatest overlap among the ranges is observed in the 2070 scenario with SSP 8.5, a projection where is observed the greatest expansion for the golden jackal, especially in the northernmost areas of the study area.Despite the reduction in habitat suitability for the Eurasian lynx in future scenarios (especially 2070), the overlap between the ranges remains high because the "losses" of the lynx are balanced by the higher "gains" of the golden jackal (Supplementary Figure 4).
The results of the Point Sampling tool underline the known differences in the habitat preference between the Eurasian lynx and the golden jackal: while the jackal tends to prefer suburban rural environments, such as cropland, the lynx tends to select preferentially dense forests.Although, this does not exclude interactions, as both species may select the other's preferred habitat.In fact, the golden jackal can select woodland habitats to search for burrows and reproduce, and its ecological plasticity allows it to use even more closed forests, especially when any apical predators are absent (Krofel et al., 2017).Still, Eurasian lynx select more open habitats (Belotti et al., 2013), even if close to humans, if they have high prey densities (Basille et al., 2009;Filla et al., 2017).Ungulates often are abundant near human infrastructures because they find protection from predators due to the proximity of man, a process known as the "human shield effect" (Berger, 2007).However, given that the golden jackal habitat selection is strongly influenced by the presence of wolves (Krofel et al., 2017), where sympatry may occur between the species (Figures 4, 5), the golden jackal habitat selection would probably change significantly.
After being the most widely distributed large terrestrial predator in the northern hemisphere, wolves went to a severe decline mainly due to human persecution (Ripple et al., 2014).Recently the species is recolonizing Europe: expanding both from its permanent range in Italy and Spain, and from the Dinaric-Balkan, Carpathian, and Polish lowland populations, recolonizing central Europe and Scandinavia, and expanding in eastern Europe (Boitani and Linnell, 2015;Hindrikson et al., 2017).
Several studies prove the wolf's strong influence on the golden jackal through intraguild predation (Mohammadi et al., 2017), and the wolf's absence is considered one of the main causes of golden jackal expansion (Krofel et al., 2017).
Our predicted sympatric areas for the golden jackal and the Eurasian lynx show that most of them overlap with current wolfpredicted suitable areas.In these areas it is more probable that the golden jackal could exploit different resources, like anthropogenic food sources, using more human-modified habitats and more open habitats.Consequently, the interactions with the Eurasian lynx would be significantly reduced.
Moreover, it should not be forgotten that Europe hosts other populations of a large predator (i.e., brown bear) and mesopredators, like the red fox.For example, the Eurasian lynx suffers kleptoparasitism from the red fox, which feeds on lynx prey (Krofel et al., 2019), and similarly, there are also cases of brown bears feeding on the prey hunted by the lynx (Krofel et al., 2012).
Future studies will be needed to fully understand the outcomes of the interaction between the golden jackal and the Eurasian lynx, an interaction already real, but which our study shows would be much more frequent in the future.Furthermore, future research will also have to better consider possible interactions, apart from kleptoparasitism, with other species of predators.
A possible limitation of our study could be linked to the dataset used, especially for the golden jackal, because we retained the records from dispersing individuals (i.e., out of the current core areas), and the main part of the dataset is obtained from the GBIF repository, known to be possibly spatially biased due to differences in datasharing between countries and the consequent unequal effort of sampling (Beck et al., 2014).This could have led to an inflated effect on the environmental suitability in the northernmost part of the study area, even though the robust design of the modeling methods we used, together with the limited number of the dispersing individuals' occurrences, buffer this issue.To check possible biases in the models' output, we repeated the ecological niche modeling process by using only records reported as "permanent" or "highconfidence presence" from an expert-based dataset.Comparing these modeling outputs (Supplementary Figures 2-5), only some subtle differences are noted, resulting in similar habitat suitability values, both for current and future conditions.Thus, this interesting finding proves that for some specific contexts (i.e., modeling habitat suitability of a generalist species at the landscape scale), the use of citizen-science data could be as useful as data from expert-based datasets.However, considering the small differences in the northernmost part of the study area between the two different ENMs (the ones built with our dataset and the one obtained from the expert-based dataset), we recommend considering this aspect when exploring our maps.However, the future availability of more robust and expert-based datasets would surely improve the quality of ecological modeling studies, and for species in rapid expansion, these datasets are also urgently needed.
In this complex context of large carnivores expansion through the European landscapes, is important to improve measures that could bring the conservation of the entire biocenosis: for example, regulating hunting activities, especially of wolves and ungulates, in order to reduce possible competition between large predators.Other measures could include the coverage increase of protected areas, which have favored the expansion of large carnivores in the past, together with an increasing public acceptance (Dressel et al., 2015), especially in areas where these species are making a comeback after some time.The use of economic compensation for any damage by predators can also help improve the public acceptance of these species (Dickman et al., 2011).
In conclusion, our results can support conservationists by identifying which areas these two carnivores could coexist in the future.Both in areas already occupied by both species and in those where sympatry is expected in the future, further studies will be needed to better understand the effects of this coexistence and any interactions.
FIGURE 1 Current and predicted distribution.Current European distribution of (A) the golden jackal (Canis aureus) and (B) the Eurasian lynx (Lynx lynx).Predicted species distribution model for (C) the golden jackal and (D) the Eurasian lynx, ranging from low (blue) to high (red) suitability.