CC-BY 4.0
Introduction
The Eurasian lynx (Lynx lynx) was eradicated in most of Europe between the 17th and 20th centuries. The main reasons for its disappearance were habitat degradation, human persecution and a decrease in prey availability (Breitenmoser et al., 2000). The species has recently recolonized parts of its historical range in Central and Western Europe thanks to different reintroduction programs which started in the 1970s. Nowadays, there are eleven identified lynx populations in Europe (von Arx, 2020), and the species benefits from a conservation status across its whole range area. The species is considered as “least concerned” at the European level of the IUCN Red list. However, its status greatly differs from one population to another, notably because of their difference in population sizes, even though these populations share similar threats, mostly habitat destruction and fragmentation, illegal killings and collisions with vehicles (von Arx, 2020). Long-term persistence remains notably uncertain for the Alpine population (France and Switzerland) and for the Upper Rhine meta-population, which encompasses the Jura population (France and Switzerland), the Vosges-Palatinian population (France and Germany) and few individuals located in the Black Forest and its surroundings (e.g., the Swabian Alb, Germany) (Drouet-Hoguet et al., 2021; Germain & Schwoerer, 2021; Herdtfelder et al., 2021; Idelberger et al., 2021; Molinari-Jobin et al., 2021). These populations are currently defined as endangered (Jura and Alpine) or critically endangered (Vosges-Palatinian) (von Arx, 2020). Indeed, the populations forming the Upper Rhine meta-population remain currently small and isolated. Only a few exchanges between them and with the Alpine population are documented, most likely because of habitat fragmentation and little functional connectivity (Zimmermann & Breitenmoser, 2007; Morand, 2016). Low female dispersal rates and movements most likely slow population expansion due to their conservative dispersal behavior compared to male lynx (Port et al., 2021). This situation, added to the common origin from the Carpathian population (Breitenmoser et al., 1998; Vandel et al., 2006) of the few reintroduced individuals which did reproduce, resulting in a bottleneck and genetic drift, may lead to high inbreeding within populations (Breitenmoser-Würsten & Obexer-Ruff, 2003; Bull et al., 2016; Premier et al., 2020). A national action plan has been defined in France with several objectives to restore the lynx to a favorable conservation status in France (Gatti, 2022).
In this context, wildlife conservationists, scientists and policy-makers face significant challenges for lynx conservation when individuals inhabit human-dominated landscapes. Several studies have used individual-based models to inform decision-making processes for population conservation (IBMs; Railsback & Grimm, 2012). These bottom-up models flexibly integrate species demography (dispersal, territory establishment, reproduction, mortality) and how animals interact with their inhomogeneous environment (e.g., habitat selection, collisions, illegal killings) and other individuals, while accounting for individual characteristics (e.g., sex, age, dispersal status). Population-level results emerge from the individual-level simulations (Railsback & Grimm, 2012). These models are used for the management and conservation of large and small carnivores (e.g., Kramer-Schadt et al., 2004; 2005; 2011; Marucco & McIntire, 2010; Hradsky et al., 2019). Recent applications of IBMs for the lynx species allowed estimating illegal killing (Heurich et al., 2018), linking movements and genetic diversity (Premier et al., 2020) and evaluating reintroduction scenarios (Ovenden et al., 2019; Philips, 2020).
Lynx are large carnivores with extensive spatial requirements (Breitenmoser-Würsten et al., 2007). They are territorial and live in vast connected forested areas (Vandel & Stahl, 2005; Zimmermann & Breitenmoser, 2007; Ripari et al., 2022), which sustain important prey populations (Basille et al., 2009). Natural barriers (e.g., broad waterbodies) as well as anthropogenic barriers (e.g., urban areas, road network) limit lynx movements (Kramer-Schadt et al., 2004). Roads are complex infrastructures acting as semi-permeable movement barriers (Klar et al., 2006; Zimmermann & Breitenmoser, 2007; Bastianelli et al., 2024) in addition to being an important source of mortality (Schmidt-Posthaus et al., 2002; Ceia-Hasse et al., 2017; Bencin et al., 2019). Therefore, a major modelling challenge when simulating lynx populations is to account for movements over large areas and consequently, for the impact of landscape characteristics on different demographic processes such as dispersal, territory establishment and survival rates (Premier et al., 2025). Simulation rules in IBMs are defined at the individual level and result in different behavior of each simulated individual due to their characteristics. Therefore, these models appear to be a great tool to integrate the complex interactions between landscape characteristics and both lynx demographic and spatial processes, to simulate population changes on a large scale.
Our study aimed to assess the long-term persistence of the Eurasian lynx in the large range area delimited by the Upper Rhine meta-population and the Alpine population, using the available population status in 2018. We build a spatially-explicit IBM that integrates up-to-date information on lynx ecology while accounting for habitat preferences and the risk of collisions with vehicles. We implement our model in the R programming language and share it publicly with open-source code allowing full reproducibility and its adoption by ecologists and conservationists. We use this model to forecast the fate of lynx populations over a 50-year time horizon. Finally, we provide several population metrics that help in diagnosing the long-term species persistence in its most western range area.
Methods
Study area and populations
We conducted the study on the Eurasian lynx populations located in France, Germany and Switzerland (Figure 1). The three main populations inhabiting the vast mountainous and forested areas of these countries are the Vosges-Palatinian (France-Germany), the Jura (France-Switzerland) and the Alpine (France-Switzerland) populations (von Arx, 2020). Some individuals were also observed in the Black Forest and its surroundings (e.g., Swabian Alb) in Germany (Figure 1), but were not considered as a lynx population per se (von Arx, 2020). Therefore, we should refer to “mountain ranges” when speaking of the Vosges-Palatinate, Jura, Alps and Black Forest with its surroundings, but for simplicity and clarity, we use “populations” throughout our paper.
The Jura and Alpine populations originated from individuals reintroduced in Switzerland in the 1970s (Breitenmoser et al., 1998; Vandel & Stahl, 2005), followed by a natural recolonization of the territories. After the complete decline of the Vosges-Palatinian population in the 18th century, a reintroduction of 21 lynx occurred in the southern part of the Vosges Mountains (France) between 1983 and 1993 (Vandel et al., 2006). Only 10 individuals contributed to the local establishment of the species (Vandel et al., 2006), without conclusive stabilization (Charbonnel & Germain, 2020). Since 2016, the lynx has been back in the Palatinate Forest (Germany), thanks to a new reintroduction program finalized in 2020, with a few individuals that already moved to the Vosges Mountains (Scheid et al., 2021; Schwoerer, 2021). Finally, only a few male individuals have been observed in the Black Forest area since 2013, most of them coming from the Swiss Jura Mountains (Stiftung KORA, 2018; Ministerium für Ländlichen Raum und Verbraucherschutz, 2019; Wölfl et al., 2021).
Figure 1 - Eurasian lynx presence as available in 2017-2018 in France (Fr.), Germany (Ger.) and Switzerland (Sw.) over the study area (black rectangle). Data for France cover the period from 01/04/2013 to 31/03/2017 (OFB Réseau Loup-Lynx), for Germany from 01/05/2017 to 30/04/2018 (Bundesamt für Naturschutz) and for Switzerland from 01/01/2015 to 31/12/2017 (KORA). We used a standardized 10x10 km grid from the European Environment Agency for France and Germany, and a grid derived from the 1:25,000 map for Switzerland). The four colors are for the four different populations: the Vosges-Palatinian population (blue; 1,800 km2) with cells in France (Vosges Mountains) and Germany (Palatinate Forest), the Black Forest population (purple; 900 km2) in Germany, and the Jura (green; 12,057 km2) and the Alpine (red; 11,190 km2) populations with cells in France and Switzerland. Top right corner inset: Europe with the location of the study area.
Lynx population model
Based on previous works by Kramer-Schadt et al. (2004; 2005; 2011), we built a spatially-explicit individual-based model (SE-IBM) to simulate lynx populations dynamics and dispersal, while accounting for the risk of lynx-vehicle collisions and lynx habitat preferences. An SE-IBM is an IBM where individual responses to behavioral rules are constrained by environmental characteristics. Our lynx SE-IBM is made of four components (Figure 2) for which the complete description is available in Supplementary Information 1. The first component represents the impact of the road network on lynx survival. Hereby, the “collision layer” (Supplementary Information 1, Figure S1.1) represents lynx-vehicle lethal collision probabilities and is based on a collision model. The second component represents the impact of land cover on lynx space use with a habitat model yielding a “habitat layer” (Supplementary Information 1, Figure S1.2) and showing the distribution of the different lynx habitat types. These two spatial layers influence the behavioral rules followed by the simulated lynx individuals. The third component represents the initial lynx populations made of lynx individual’s location (Supplementary Information 1, Figure S1.3) and characteristics (e.g., age, sex). The fourth component details all SE-IBM rules including lynx demography, dispersal movement and territory establishment. The four components were combined to create an SE-IBM and simulate lynx populations. The initial lynx populations are the starting individuals created at the beginning of the simulation. Simulated individuals are either resident (i.e., established with a territory) or dispersers (Supplementary Information 1, Figure S1.4). Dispersing lynx move through the landscape, driven by a preference towards certain habitat types defined by the “habitat layer”, a constant direction and their memory (Supplementary Information 1, Figure S1.5). Dispersers can die with a fixed annual mortality probability or from vehicle collisions (additive spatial mortality), defined using the “collision layer”. When arriving at a new location, dispersers look for a territory (Supplementary Information 1, Figure S1.6). Females need a large area of good quality habitat to establish, whereas male establishment is driven by the presence of available resident females. A female resident may reproduce if there is a resident male available nearby. Residents may also die with a fixed mortality or by vehicle collision. This sequence represents one year and is repeated with the updated individuals (e.g., newborns added, dead individuals removed, individuals at different locations, new residents) as long as desired. A complete description of the SE-IBM structure and processes following the Overview, Design concepts, and Details (ODD) protocol (Grimm et al., 2006, 2010) is provided in Supplementary Information 2. The full list of model parameters with their values is also provided in Supplementary Information 2 (Supplementary Information 2, Table S2.1). We calibrated parameters to fit the model to our studied populations better. We also performed a sensitivity analysis to evaluate the changes in the main results following the changes in parameter values. Full details on model calibration, sensitivity analysis and validation are available in Supplementary Information 3.
Figure 2 - The four components of the lynx SE-IBM. The first component represents the impact of the road network on lynx survival through vehicle collision. A generalized linear model predicts collision probabilities (“Collision layer”, Supplementary Information 1, Figure S1.1). The second component represents the impact of land cover on the lynx populations. A site-occupancy model predicts lynx habitats influencing lynx movement (“Habitat layer”, Supplementary Information 1, Figure S1.2). The third component represents lynx populations with individuals’ locations and characteristics (Supplementary Information 1, Figure S1.3). The fourth component encompasses all ecological rules and SE-IBM parameters. All four components are included in the SE-IBM to simulate lynx populations and assess their persistence with different metrics.
Population persistence
We ran 100 replicates of the lynx SE-IBM to forecast the populations over 50 years. We used a different initial population (Supplementary Information 1) at each replicate to avoid bias due to the initial locations of simulated individuals. Initial populations were created based on real monitoring data (Supplementary Information 1). We defined a burn-in phase of 2 years after the start of the simulation (i.e., only simulation outputs from the 3rd year of simulation were included in the conclusions of the analyses) to let the population settle down as all individuals from the initial population were defined as dispersers. We did not simulate any mortality during the first year of simulation to allow all individuals to find a territory, if they can, without dying. To evaluate lynx persistence, we analyzed a) population growth rates, b) number of movements between populations looking at establishments outside of the lynx native population, c) female territory occupancy, and d) lynx density.
a) We extracted the number of individuals in each population, for every year and each replicate. Then, we calculated the growth rate for each replicate and each population over the simulated years as the number of individuals at time t divided by the number of individuals at time t-1. Finally, we calculated the mean growth rate per year and population, and the 95% confidence interval of the mean over the 100 replicates.
b) We extracted the number of times individuals established their territory in a population area not corresponding to the one they were born in (“Population layer”, Supplementary Information 1) for every year and for each replicate. We did not account for individuals which moved to another population area during their dispersal but finally came back to their native population area to establish, nor those which died in another population area while dispersing. We calculated the number of lynx establishing outside of their natal population per year for each replicate and direction (e.g., from the Alpine to the Jura population and vice versa, from the Jura to the Vosges-Palatinian population and vice versa). Finally, we calculated the mean and 95% confidence interval of the mean across all replicates.
c) We extracted female territory locations at the last year of simulation for each replicate. We focused on female territories as male territories are based on those of the females. We assigned a “1” to each cell of the gridded study area (1 km2) that was occupied by a female territory and a “0” otherwise. Finally, we estimated territory occupancy by calculating the mean value per cell from these rescaled maps over the 100 replicates. We also visually validated our model predictions using GPS and VHF tracks from several collared female residents. To this end, we overlaid telemetry locations on the resulting territory occupancy map to see if they mostly fall in cells often selected by simulated females, to establish their territory (details in Supplementary Information 4).
d) We extracted the locations of all adult residents in the last year of simulation for each replicate. As density measure, we computed the mean number of lynx on each 1 km2 cell over all replicates. Finally, we summed these means at the 100 km2 resolution cells to obtain lynx density values at this scale.
Results
Population growth rates
Simulations predicted similar growth rate patterns for all lynx populations, with growth rates above one (i.e., growing phase) slowly decreasing towards reaching one (i.e., stabilization phase at saturation) along the 50 years of simulation (Figure 3). The Alpine and Jura populations had similar patterns with their maximum growth rates equal to 1.09 (sd = 0.04) and 1.05 (sd = 0.06) at the 10th and 13th year of simulation, respectively. The Jura population quickly reached a stabilization phase around the 20th year and fluctuated a little over the final period of the simulation. The Alpine population growth rate similarly decreased rapidly, reaching a low value, but still significantly over 1. The Vosges-Palatinian and Black Forest populations had more fluctuating growth rates over the simulation. It reached a maximum of 1.18 (sd = 0.51) and 1.45 (sd = 0.93) at the 23rd and 14th year of simulation, respectively. Confidence intervals around mean growth rates were much larger than those of the Alpine and Jura populations. The Vosges-Palatinian and Black Forest populations also had their growth rates decreasing after reaching their maximum, until reaching values close to 1 but only towards the end of the simulation period.
Figure 3 - Annual rates of increase over the simulated years for each population. The grey area represents the 2-year burn-in phase at the beginning of the simulation. Points are means over 100 replicates and colored envelopes are 95% confidence intervals from 100 replicates. The dashed line represents a stabilization of the population (growth rate equals to 1).
Lynx establishments outside their native populations
The number of lynx establishments in a population area different from their native one showed that movements between populations were rare for the four populations. Indeed, the highest values were 2.5 lynx per year on average, concerning individuals born in the Alps and successfully established in the Jura when reaching 25 years of simulation (Figure 4.a). The Jura population seemed to be at the center of lynx movements among populations in Western Europe (Figure 4). However, the number of individuals born in Jura and establishing elsewhere decreased around the 20th year of simulation, matching when the population started to stabilize (Figures 3 and 4.b). There was almost no exchange of individuals between the Alpine and Vosges-Palatinian populations, in any direction (Figure 4.a and 4.c); as these are the most distant populations, the other ones must serve as stepping stones. Surprisingly, the Black Forest showed an exchange of individuals mainly with the Vosges-Palatinian population, even though the Black Forest seemed evenly distant to all three other populations (Figure 4).
Figure 4 - Number of lynx establishments outside their native populations per year when born in the a) Alpine, b) Jura, c) Vosges-Palatinian or d) Black Forest population. Colors represent different directions from the native populations (each figure) to those where individuals established (different colored lines). The grey area represents the 2-year burn-in phase. Points are mean over 100 replicates and colored envelopes are 95% confidence intervals around the mean from 100 replicates.
Territory occupancy
At the end of the simulation period, female territories occupied most of the study area (Figure 5). In the core areas of the four populations, female occupancy was high, reaching values above 0.85. The overlaid VHF and GPS paths of collared resident females (Figure 5 and Supplementary Information 4) matched well the limits of the most frequently occupied cells. Most patches far from these core areas were never occupied over the 100 replicate simulations (Figure 5). Nevertheless, there were distant patches, especially the bigger ones, occupied by female territories with low probabilities (Figure 5).
Figure 5 - Occupancy by female territories over the study area in the last year of simulation. Values, ranging from zero (excluded) to one represented by the color scale blue-yellow-red, are mean occupancy probability per cell of 1 km2 over 100 replicates (e.g., a cell equal to 1 means the cell was included in a female territory, at the 50th year of simulation, in all 100 simulation replicates). Grey cells are “breeding habitats” never included in a female territory across all replicates (equal to zero). White cells are cells of a habitat type different than “breeding habitat” and therefore could never be included in a female territory. GPS and VHF recorded paths for female residents are overlaid as thin black lines. A zoom of the different areas with telemetry data is presented in Supplementary Information 4. Axis (x and y) are in meters from the WGS 84 / Pseudo-Mercator coordinate system.
Lynx density
Adult resident density ranged from 0 to 1.69 lynx per 100 km2 over the study area, in the last year of simulation (Figure 6). The highest density value was in the Alpine population and this population also had the highest mean density over the four populations (0.55, sd = 0.4). The Jura population had the lowest density values (max = 0.97, mean = 0.26, sd = .024). The Vosges-Palatinian and Black Forest had intermediate values (max = 1.13 and 1.16, mean = 0.33 and 0.47, sd = 0.29 and 0.34, respectively).
Figure 6 - Adult resident lynx density per 100 km2 in the last year of simulation. Lynx density (color scale blue-yellow-red) are mean values over the 100 simulation replicates. White cells had a density equals to 0. Axis (x and y) are in meters from the WGS 84 / Pseudo-Mercator coordinate system.
Discussion
The well-established Alpine and Jura populations
The model predicted a steady growth for the Alpine and Jura populations, quickly reaching a more stable phase, indicating that they may soon be at saturation, especially for the Jura population. However, saturation could be higher than what we found in our simulations, depending on the strength of density-dependence processes (Zimmermann et al., 2007). Although our model assumed that the number of individuals an area can support is mainly defined by female territory size, it does not include a relationship between lynx behavior and density. Yet, studies found that lynx density also influences home range sizes (Pesenti & Zimmermann, 2013), and differently for females and males (Aronsson et al., 2016). For instance, in the North-Western Swiss Alps, lynx territories were much larger in the 80s when lynx density was lower compared to the 90s (Breitenmoser-Würsten et al., 2001). Simulated lynx density in 50 years was a little smaller than those estimated in the field for both populations (Gimenez et al., 2019; Pesenti & Zimmermann, 2013). Our model did not seem to allow a significant increase in lynx density. The Alpine and Jura areas could therefore support more individuals if such density-dependent mechanisms would occur.
The Jura population as the crossroad of Western Europe lynx movements
The Jura population was found to be the only one connected by lynx exchanges with all the other populations which is coherent with observed data, also in terms of the number of individuals moving per year (Drouet-Hoguet et al., 2021). Individual exchanges simulated were low, with only a couple of individuals moving from their native population to establish their territory in another one on average per year. Exchanges were predicted to be the highest for the individuals born in the Alps and moving to settle down in the Jura Mountains, with two to three individuals per year after 25 years of simulation. In Switzerland, only a few lynx movements between the Alpine and Jura populations have been observed despite a large camera-trapping effort over the years. But even though camera trapping can detect movements between populations, it cannot guarantee an exhaustive sampling nor ensure that individuals successfully established in their destination, which contrasts with the computed results. However, Alpine and Jura populations also differ genetically, suggesting very few exchanges between them (Breitenmoser-Würsten & Obexer-Ruff, 2003). Our model may slightly over-estimate these dispersal movements due to the absence of movement barriers in the “habitat layer”.
When inspecting model output maps (Figures 5 and 6), individuals are contained in the Swiss Alps, restricted to an area of suitable habitat almost totally surrounded by less favorable habitats (Supplementary Information 1, Figure S1.2). However, individuals have recently started to settle permanently on the Swiss Plateau (between the Alps and the Jura Mountains) predicted with less favorable habitats. Some have even reproduced successfully (F. Zimmerman, pers. comm.), indicating potential connections between the two populations in the future. On the other hand, our model suggests connectivity between the Jura and Alpine populations on the French side with continuity of a few forested corridors until the Chartreuse Mountains in the Alps (Zimmermann & Breitenmoser, 2007). Many observations were made via camera traps of lynx moving between the southern part of the Jura population and the Chartreuse Mountains (Bailly, 2021). However, movement barriers (e.g., highways) may prevent connections with the rest of the French Alps and the Alpine populations in Switzerland. Illegal killings, notably in Switzerland in movement corridors as well as a slow demographic dynamic of the Swiss population may also have prevented a larger development of the lynx population in the French Alps (Arpin et al., 2024).
Lynx movements simulated from the Jura to establish within the Vosges-Palatinian population were possible but rare. This is coherent with field monitoring that highlighted the presence of one male in the south of the Vosges Mountains, coming from the Jura Mountains during the winter 2014-2015 (Chenesseau & Briaudet, 2016; Hurstel & Laurent, 2016). Connectivity between these two populations remains far from optimal because of multiple barriers (e.g., highways, railroads and rivers) and restricted habitat that impede lynx dispersal between the two mountain ranges and increase collision risk (Zimmermann & Breitenmoser, 2007; Morand, 2016).
Regarding the Black Forest (including the adjoint Swabian Alb area), six males are known to have immigrated from the Swiss side of the Jura population since 2013 (unpublished data, Forest Research Institute Baden-Wuerttemberg; Stiftung KORA, 2018; Drouet-Hoguet et al., 2021), as well as two males from the Northeastern Alps (Herdtfelder et al., 2021). Simulated lynx movements towards the Black Forest seem therefore coherent with field data. After 25 years of simulation, it is the Vosges-Palatinian population which seemed the most connected to the Black Forest. This is concurrent with the growth of the Vosges-Palatinian population, which could act as a source population only after several years, as it started only with a few individuals. It lasted several years until the first lynx immigrated to the Black Forest from the other populations in our simulation compared to the others. This is due to male establishment which is driven by female presence in our model, and because the Black Forest population did not have any female at first, new males arriving could not establish and were doomed to die while searching for females. A more realistic rule for male establishments would be to let them search for females but still allow them to settle after a defined period, even without females. It happens quite regularly that males show territorial behavior even without females (M. Herdtfelder, pers. comm.).
The increasing Vosges-Palatinian and Black Forest populations
Projections for the Vosges-Palatinian and Black Forest populations showed positive growth over the next 50 years. However, there was a difference in their initial populations. The lynx reintroduction program in the Palatinate Forest conducted since 2016 constituted a new larger initial population for the Vosges-Palatinian one (Scheid et al., 2021; Schwoerer, 2021) compared to the few individuals, only males (Drouet-Hoguet et al., 2021), located in the Black Forest. Growth rates reached higher values for the Black Forest population, but growth rates were very heterogeneous along the simulated period and confidence intervals around the mean were very large. Demographic stochasticity greatly impacted the Black Forest population because of its small initial population. Precision, as well as robustness, to interpret the results for this population were indeed lower than for the Vosges-Palatinian population. In the field, one female was reintroduced in the Black Forest by the end of 2023 to support this only-male population, but it died in Summer 2024 (www.baden-wuerttemberg.de), highlighting the high stochasticity of the recolonization of this habitat patch.
Exchanges of individuals between the Vosges-Palatinian population and the Black Forest one increased only at the end of the simulation, probably when both populations had more individuals. The two also seemed connected with the Jura population, highlighting the Upper Rhine meta-population structure. However, functional connectivity within the meta-population may not be at its best and may be altered by anthropogenic barriers. For example, in the Vosges-Palatinian population, at the Col de Saverne, this forest bottleneck is fragmented by both a highway and high-speed railway (Klar et al., 2006; Morand, 2016; Scheid et al., 2021). However, movements are possible as lynx from the Palatinate Forest are known to have crossed this pass, some even going back and forth (Idelberger et al., 2021; Scheid et al., 2021). Moreover, no female lynx has been observed until now dispersing from the Jura Mountains to the Black Forest crossing the Rhine valley that separates the Black Forest population from the Jura one, probably due to their conservative dispersal behavior compared to males (Port et al., 2021). Monitoring data from 2021 indicate that one female from the north-eastern Swiss population crossed this barrier for the first time (M. Herdtfelder, pers. comm.).
Model limitations
Although predicted growth rates did not seem very sensitive to model calibration and model predictions were well validated (Supplementary Information 3 and 4), several aspects of our IBM could still be improved.
The “habitat layer” is defined in categories with strict “barriers” that the simulated lynx cannot cross, whereas the species may be tolerant of human activities (Basille et al., 2009; Bouyer et al., 2015). Lynx may live near urban areas, cross small lakes (F. Zimmermann, pers. comm.), and sometimes large rivers. In that context, the “habitat layer” could be improved by being defined as a continuous variable representing a degree of selection or permeability of the landscape to avoid prohibiting the movement through certain landscape elements. The “habitat layer” could also be improved by accounting for roads and their associated structures as movement barriers (e.g., highways and other types of roads difficult to cross, over-pass facilitating the movement) instead of only mortality sources (Klar et al., 2006; Marchand et al., 2017).
In our model, lynx movement is not impeded by roads. Including permeability of these linear barriers to the lynx movements may help to refine the behavior rules (Marchand et al., 2017). We could also redefine the movement behavior rules to be sex or age-dependent. Indeed, females seem more conservative (i.e., disperse close to their natal range) compared to males and some males can disperse over long distances (Port et al., 2021). Moreover, lynx movement and therefore the impact of habitats and collision risk were modeled only for dispersing individuals. Movements inside territories were not simulated for resident lynx and we assumed an even use of the space. Model predictions could be different if resident individuals would move across their territory in an uneven way, increasing or decreasing the averaged collision risk of their home range. However, mortality risk for residents, thanks to their knowledge about their territory where and how to cross linear infrastructures, would probably be lower than for dispersing individuals.
Moreover, our model predictions were made in a static environment, as the landscape did not change over the simulated period. Modifying the landscape (e.g., agricultural intensification, forest cover increase) would give different lynx population predictions. However, a landscape model would be necessary to predict a realistic future landscape, or different scenarios could be defined (Bauduin et al., 2018). We used a medium-term simulation period of 50 years in which a static landscape was a reasonable assumption.
Finally, we did not include any genetic aspect in the demographic rules of our SE-IBM (Mueller et al., 2022). A collective expertise recently done in France on lynx viability showed that inbreeding depression could greatly impact population size and extinction risk (Arpin et al., 2024). We restricted our predictions over a 50-year period to limit the (non-included) effects of inbreeding depression on population dynamics. However, even within this time frame, population trends could shift from increasing when inbreeding depression is excluded, to potentially decreasing when included. Therefore, our results need to be interpreted carefully and more in lights of population relative differences than raw values for each population per se.
Perspectives for assessing lynx conservation strategies
Our model used to evaluate lynx population persistence could be applied to better understand aspects of lynx conservation that we did not include yet, such as illegal killings (Heurich et al., 2018). Moreover, thanks to its individual-based structure, a genetic component could easily be included to track relatedness between individuals and allow studying inbreeding risk and Allee effects (Arpin et al., 2024; Premier et al., 2020). Our model could also be used to test the effect of different scenarios, either by modifying the lynx population sizes (e.g., due to illegal killings or reintroduction programs) or the landscape at different scales (e.g., green bridges, roads construction, habitat destruction or restoration) and assess the potential benefits or negative effects on lynx population persistence. Modifications of road networks to improve connectivity, such as the removal of road segments or the addition of overpasses or underpasses could be tested. Then, their effect could be included through an updated layer of collision probabilities and the population persistence calculated accordingly. Similarly, other modifications of land cover (e.g., restoring forest areas) could also be tested. This model could be very useful for stakeholders working on corridors and the reduction of lynx-vehicle collisions as well as reintroduction programs and species acceptance. Helping reducing road mortality using model simulations is one of the objectives of the French national action plan (Gatti, 2022).
Using IBMs, Herdtfelder (2014) showed that Black Forest population reinforcement with lynx females might be one solution considering that habitat is of good quality in this area. The collective expertise on lynx viability also used a model similar to the SE-IBM presented here to make recommendations on lynx translocations in the Vosges area to support the population viability (Arpin et al., 2024). IBMs are great tools to evaluate the impact of such management actions. However, because our model does not include human attitudes towards the species, we warn against its blind use to assess the effect of reinforcement on lynx long-term persistence without accounting for the human dimension component. Illegal killings, as they occurred after the 90’s reintroduction program in the southern part of the Vosges Mountains (Vandel et al., 2006), as well as more recent ones (Germain, 2020), may be an additional mortality that the model does not account for. The extent of acceptance towards species reinforcement from some local stakeholders is an essential element (Charbonnel & Germain, 2020). In this case, we recommend extending our model to include the dynamics of the whole socio-ecosystem (Behr et al., 2017; Guerrero et al., 2018).
Conclusion
In this paper, we built and analyzed a spatially-explicit individual-based model to forecast the fate of lynx populations over the next 50 years. Our results suggest that exchanges of individuals between Vosges-Palatinian, Black Forest, Jura and Alpine populations to establish new territories were limited, and emphasize that the Jura population plays the role of a Western Europe crossroad. Overall, lynx persistence in the Upper Rhine meta-population and the Alpine population over the next 50 years seems likely on a large scale by considering road mortality and habitat selection.
Acknowledgements
We thank OFB and Réseau Loup-Lynx for the lynx presence and lynx-vehicle collision data to build the habitat and collision models; OFB, Réseau Loup-Lynx, KORA and Bundesamt für Naturschutz for population data to create the initial populations; OFB, Réseau Loup-Lynx and KORA for the collision data for model validation; and KORA and Stiftung Natur und Umwelt Rheinland-Pfalz for the female telemetry data for model validation. We thank Elodie Vercken, Hector Ruiz and Henrik Andren for their very constructive reviews. Preprint version 4 of this article has been peer-reviewed and recommended by PCIEcology (https://doi.org/10.24072/pci.ecology.100408; Vercken 2025).
Funding
SB and OG were funded by French National Research Agency (ANR‐16‐CE02‐0007). SB was funded as well by OFB and OG was funded by CNRS and the “Mission pour l’Interdisciplinarité” through the “Osez l’Interdisciplinarité” initiative. CEFE, Cerema and CROC were funded by CILB, MTES (ITTECOP) and FRB through the research program ERC-Lynx. Cerema received support from the “Direction générale des infrastructures de transports et de la mer (Ministère de la transition écologique)”. CROC was funded in 2019/2020 by the European Union within the framework of the Operational Program FEDER-FSE “Lorraine et Massif des Vosges 2014–2020”, the “Commissariat à l’Aménagement du Massif des Vosges” for the FNADT (“Fonds National d'Aménagement et de Développement du Territoire”), the DREAL Grand Est (“Direction Régionale pour l’Environnement, l’Aménagement et le Logement”), the “Région Grand Est” and the “Fondation d’entreprise UEM”. This research is partly a product of the DISCAR group funded by the synthesis centre CESAB of the French Foundation for Research on Biodiversity.
Conflict of interest disclosure
The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.
Olivier Gimenez is a PCI Ecology recommender.
Data, scripts, code, and supplementary information availability
Data, scripts, code and supplementary information are available online: https://doi.org/10.5281/zenodo.15102511 (Bauduin et al., 2025).