A decade of systematic camera trapping in two strictly protected areas reveals the demography of a Eurasian lynx (Lynx lynx) population in Central Europe

Large carnivores are crucial for ecosystems but are increasingly threatened by human persecution and habitat destruction. Successful conservation of this guild requires information on long-term population dynamics through the demographic surveys. We camera trapped Eurasian lynx between 2009 and 2018 in two strictly protected areas in the Bohemian Forest Ecosystem, estimating sex-specic demographic parameters using spatial capture-recapture (SCR) models as well as the relative abundance index of its prey species and a mesopredator seeking potential drivers of lynx population dynamics. Over 48 677 trap nights, we detected 65 lynx individuals. Density increased to 1.31 and 2.39 individuals/100 km 2 for open and closed population SCR models, respectively, with positive population growth rates. Survival rates were high (females 83%, males 88%) and per capita recruitment was low (females 13%, males 9%), indicating a low yearly population turnover. Reproductive parameters showed successful reproduction. Our results reveal the importance of the study area as host to a saturated population and key source for the potential lynx metapopulation. The observed increasing lynx number is most likely represented by dispersing individuals due to reduced mortality outside the protected areas as the number of reproducing females inside remained constant. In what is the rst open population SCR study on lynx, we provide demographic parameters contributing to the development of model-based population viability forecasting and inform lynx management in the study area as well as in similar contexts.


Introduction
Large carnivores shape ecosystems through top-down control of herbivores and intraguild predation of mesocarnivores, triggering trophic cascades 1 . Their role is therefore crucial for the functioning of ecosystems. However, large carnivore conservation is particularly challenging. Carnivore lifestyles can generate con icts with human activities 2 and their wide-ranging behaviour necessitates complicated transboundary management 3 . Moreover, low densities and reproduction rates make these species vulnerable to the effects of human persecution 4

.
A deep understanding of demographic parameters such as abundance, density, survival and growth rate, is crucial for large carnivore conservation (e.g. Breitenmoser-Würsten et al., 2007), providing information on population dynamics and forming a basis for population forecasting using model-based predictions. However, these parameters can only be calculated over substantial time intervals, hence long-term data are required for meaningful studies 6 .
The use of camera traps as a wildlife monitoring tool has gained popularity in recent years (O'Connell et al., 2011). These non-invasive devices monitor different species simultaneously and avoid stressful animal immobilisation 8 . Camera traps can provide high-quality pictures enabling the identi cation of naturally marked animals 9 . This has led to their extensive use in combination with capture-recapture (CR) methods in demographic studies. Recently developed spatial capture-recapture (SCR) models present improvements over conventional non-spatial CR models in estimating demographic parameters because they incorporate spatial information (i.e. location of individuals and traps) 10 . Closed population SCR models assume demographic closure (i.e. constant population size) and are normally used to estimate abundance and density within one session 10 . More recent open population SCR models provide the advantage that they can be applied across multiple sessions providing further parameters such as survival, per capita recruitment and population growth rate (e.g. Chandler and Clark, 2014).
In the last century, the Eurasian lynx (Lynx lynx), was eradicated across Central Europe but, following legal protection, the establishment of protected areas and population reintroductions, have recolonised parts of its former range 12 . However, most of the reintroduced populations in Central Europe remain isolated and small in size 12 . As a typical example, the Bohemian-Bavarian-Austrian (BBA) population is classi ed as "endangered" and exhibits reduced genetic variability 13 . Müller et al. (2014) showed that the Bohemian-Bavarian Forest (BBF) could serve as an important source area within a Central European lynx metapopulation although subadults cannot connect with neighbouring populations (i.e. German Harz, Czech Carpathian, and Austrian Alpine). The illegal killing has been implicated as the main mortality cause driving this population and the protected areas in the region as the crucial factor for its persistence 4 . The protection of source populations is proposed as a strategy for recovering species such as the tiger (Panthera tigris) 15 . However, whether protected areas in Central Europe, with their limited size, can act as sources is poorly understood.
Here we conducted a demographic study on a Eurasian lynx (hereafter lynx) population in one of Central Europe's largest strictly protected areas, the Bohemian Forest Ecosystem (BFE). Studying this reintroduced population helps to improve the understanding of lynx ecology across its range since protected areas provide the opportunity to investigate population dynamics under low human disturbance. We systematically camera trapped lynx from 2009 to 2018 in the BFE. To our knowledge, this is the rst demographic study on lynx comprising annual camera trapping over a decade and using open population SCR models. We expected higher densities of females since males have larger home ranges 16 . As lynx is a K-selected species, we expected generally high survival and low per capita recruitment rates. Speci cally, we expected lower male survival rate because they generally take higher risks by getting closer to human activities due to high prey density 17 and patrolling their larger home ranges. We expected an average litter size of about 2 kittens, in line with the optimum number shown for this species 18 . Finally, we used the relative abundance index (RAI) 19 of lynx's main prey species and a mesopredator to explore potential in uences of the prey base on lynx.

Study area
The BBF is situated in Central Europe at the border between Austria, Czechia, and Germany, and includes two adjacent national parks: the Bavarian Forest National Park (BFNP) (240 km²) in Germany and the Šumava National Park (SNP) (690 km²) in Czechia (Fig. 1). The rst is surrounded by the Bavarian Forest Natural Park (3007 km 2 ) and the second by the Šumava Protected Landscape Area (1000 km 2 ), together forming the BFE. Elevation ranges from 600 to 1456 m.a.s.l. and snow cover can persist for 5-8 months with greatest depths from January to March. The area is covered by a mixed mountainous forest composed mainly of Norway Spruce (Picea abies), followed by European Beech (Fagus sylvatica) and Silver Fir (Abies alba) 20 and hosts ungulate species such as roe deer (Capreolus capreolus), red deer (Cervus elaphus), wild boar (Sus scrofa) and moose (Alces alces) 21 . In the BFNP roe deer density ranges from 1.1 to 5 animals/km² 22 and that of red deer was estimated as 1.56 animals/km² via coordinated counting at winter feeding stations. Densities of both species are higher in the SNP 21 . Wildlife control within both national parks is conducted by trained staff and is limited to red deer and wild boar, while roe deer, red deer and wild boar are killed by hunters outside their boundaries 21 .
The BBA lynx population originates from 5-10 individuals illegally reintroduced in the 1970s followed by an o cial reintroduction of 17 individuals (11 males and 6 females) in the 1980s (Wöl et al., 2001).
Population size and density were recently estimated as 97-143 24 and 0.4-0.9 individuals/100 km² 25 , respectively. This predator is legally protected in all three countries but illegal killing still occurs 4 . The main prey species of lynx is roe deer followed by red deer (80% and 17% of kills, respectively) 26 .

Camera trap monitoring
Data were collected between 2009 and 2018. The monitoring design was developed for lynx but also provided data on other species. The spatial organisation of camera traps was designed by Weingarth et al. (2012). This consisted of a 2.7 x 2.7 km grid ( Fig. 1) maximising detection probability and avoiding gaps that might include female home ranges, i.e. a minimum of 122 km 2 estimated in the study area with radiotelemetry 27 , where camera trap sites were located in every second cell along forest paths, roads and trails. The grid occupied the entire BFNP and two-thirds of the SNP. Since individual identi cation requires highquality pictures of both animal anks, most sites included two opposing cameras (Cuddeback®) with white ash. In the BFNP, camera traps functioned almost constantly for 10 "lynx years". A lynx year is de ned as the period from 01/05 to 30/04 since kittens are mostly born in May and leave the mother in April of the following year 28 . In the SNP, most of the camera traps were only active from mid-September to the end of December. Camera trapping sites were moved slightly at the beginning of each session and changed in number.
We identi ed lynx individuals from images by comparing their unique coat patterns. Sex was determined by observing females with kittens or the genital area of the animal. Age was di cult to assess except for individuals rst photographed as kittens. Therefore, we followed Weingarth et al. (2012) and assigned animals into two categories referring to their status: "juvenile" and "independent". The juvenile category included individuals < 1-year-old, i.e. kittens detected with their mother. This category was discarded from the SCR analysis because of the high mortality of kittens (e.g. Andrén et al., 2002). The independent category included all individuals > 1-year-old (i.e. subadults and adults) and individuals of unknown age but with proof of independence. All individuals of unknown status were discarded from analyses.

Demography of lynx
Demography was investigated using both open and closed population SCR model frameworks (Supplementary Table S1). Although the former method allows estimation of the full range of demographic parameters, we included the latter for backward compatibility with previous studies conducted on lynx in Europe using closed population SCR models (e.g. Gimenez et al., 2019), as well as for comparison of both methods within one study.
Following Weingarth et al. (2015), we selected a time frame of 100 days from 15/09 to 24/12 of every year ensuring demographic closure as a primary period. We combined all primary periods according to a classical "robust design" 33 , and thus performed open population SCR analysis. We de ned one secondary period as one day and restricted the number of detections to at most one per site in any secondary period in line with a Bernoulli distribution, thereby reducing temporal autocorrelation. Our study included ten primary periods with 100 secondary periods each. The SCR method assumes the baseline detection probability g 0 of any individual declines with the distance from its theoretical home range centre, the detection function scale σ, in the state space S, which should be at least 3σ, 10 . We, therefore, created a state-space mask with a continuous buffer of 18,000 km 2 around camera traps based on results of preliminary closed SCR analyses (i.e. σ ~ 3.5-4 km). We used the density-independent population growth model in a Bayesian framework developed by 11 and an adaptation to sex-speci city of the demographic parameters by Satter et al. (2019). Using the R package "OpenPopSCR" 35 , we ran 10 pooled chains for 25,000 Markov chain Monte Carlo iterations each with an augmented observed population size of 400 individuals, i.e. much greater than the overall number of independent individuals detected across all primary periods. We estimated sexspeci c density and used a Markov activity centre relocation type 36 between primary periods to estimate male and female yearly survival and per capita recruitment separately from emigration and immigration, respectively 37 . The yearly movement between primary periods was estimated for combined males and females. Per capita recruitment estimates indicated the number of individuals per sex class added to total abundance each year. Sex-speci c population growth rate and sex ratio were estimated as derived parameters. The estimates of the latter indicated the probability of any individual to be a female. All point estimates were obtained using posterior modes and interval estimates were calculated through 95% highest posterior density (HPD) intervals using the 2.5% and 97.5% quantiles of the posterior distribution 10 .
Closed population SCR models in a maximum likelihood framework were tted in the R package "secr" (Efford, 2019). We de ned one sampling occasion as ve days 9 and used detector type proximity (Efford, 2019) according to a Bernoulli distribution to estimate the combined density of males and females. We used a buffer of 13'000 km 2 according to secr buffer predictor to create state-space masks with continuous and discrete habitat. For the latter, we followed Rovero and Zimmermann (2016) and created a habitat suitability mask excluding highly unsuitable patches 27 . Sex was tested as a covariate on g 0 and σ. We followed Akaike's Information Criterion corrected for small sample sizes (AICc) to determine the best model for our data and kept those with ΔAICc < 2 for model-averaging 40 . Statistical signi cances were evaluated using the 95% con dence intervals (CI).
We calculated generation time as the mean age of reproducing females (i.e. females with kittens) at their rst recorded litter including only those of known age (assessed via all recaptures), while all known reproducing females were included for calculating average litter size. Regarding the latter, we included all available data referring to complete lynx years, since some kittens are detected later in January/February. Finally, we included all the independent individuals of known age and sex detected during the primary periods for investigating the age distribution.

Camera trap monitoring
We detected 65 unique independent individuals (25 males, 28 females and 12 individuals of unknown sex) in a total of 48,677 trap nights ( Table 1). The standardised session length of 100 days was not achieved in all primary periods, especially in the earlier ones, due to technical problems such as camera trap failure, risk of theft or snowfall. We observed increasing juveniles and independent individuals including males, females and individuals of unknown sex while the number of reproducing females (i.e. females with kittens) remained quite stable. The age status was determined in most cases except for few individuals in the rst and third sessions. We were able to con rm mortality causes of some detected individuals in the vicinity of the study area, which showed an increase in tra c accidents and occasional illegal killing. In particular, tra c accidents killed ve independent individuals and two juveniles in the BFNP, while illegal killing affected two independent individuals and one juvenile of unknown sex outside its boundaries, in Germany.

Demography of lynx
The combined density of independent males and females estimated for each primary period using open population SCR models ranged from 0.68 (95% HPD intervals 0.48-1.07) to 1.31 (95% HPD intervals 1.03-1.76) individuals/100 km² while, as expected, sex-speci c density estimates were on average higher for females (Fig. 2). However, the overlapping 95% HPD intervals suggest there was no signi cant difference between sexes.
The parameter posterior modes (Table 2) of the yearly baseline detection probability g 0 were equal between sexes. The yearly capita recruitment was higher for females. Contrary to our expectation, the yearly survival rate was higher for males. As expected, males showed a higher yearly detection function scale σ, due to larger home ranges. However, considering the overlap of 95% HPD intervals, none of the parameters indicated a signi cant difference between sexes.  Table S4). As expected, the latter was close to the lynx optimum number.
The age distribution pyramid included 24 independent individuals of known age and sex and was slightly female-biased in size (15 females vs 9 males) (Supplementary Figure S1). The oldest individual of known age was a 10-year-old male individual. The number of females at each age was equal to or higher than that of males except for ages 8-10.

Relative abundance index
The RAI of red fox ranged from 0.07 (SD 0.01) to 0.14 (SD 0.02), that of red deer from 0.01 (SD 0.00) to 0.07 (SD 0.01) and that of roe deer from 0.00 (SD 0.00) to 0.04 (SD 0.01) (Fig. 3). Fox abundance increased from 2009 to 2014 and decreased in the remaining years. Red deer and roe deer slightly increased until 2014, after which the former almost doubled while the latter decreased to earlier values.

Discussion
Our results prove the number of independent lynx in the BBF has slightly increased over the past decade with a positive population growth rate. Survival rates were high and per capita recruitment was low indicating a low yearly population turnover. Reproductive parameters showed successful reproduction occurring every year, con rming the area represents a key source for the surrounding landscape. The yearly baseline detection probability g 0 was equal between sexes ( Table 2). Although males are generally more active patrolling larger home ranges, females with kittens hunt at a higher rate resulting in augmented activity 43 . We found a higher yearly detection function scale σ for males (Table 2), corresponding to larger home ranges, which in line with, e.g., Schmidt et al. (1997) and common in felids Open population SCR models in a Bayesian framework allowed us to separate survival from emigration and recruitment from immigration. Although this approach is still under development due to large camera trapping array requirements 37 , it is reasonable to assume that true survival is at least as high as the value we estimated when the underestimation due to the inclusion of emigration is considered. Our combined yearly survival rates for subadults and adults (population mean = 85%) are among the highest observed for lynx in Central Europe. This is because our estimates come from strictly protected areas where animals have a higher chance to survive. However, the only values available for comparison are from telemetry studies, which are generally conducted with smaller sample sizes because of high-cost devices but allow to track animals' fates. In Switzerland, Breitenmoser-Würsten et al. (2007) found an overall survival rate of 76% for adults and 53.3% for subadults. In Poland, the combined survival rate of subadults and adults was 63% 44 . In three Scandinavian study sites, Andrén et al. (2002) found higher survival rates for subadults (70, 77 and 71%) and adults (87, 91 and 84%), likely due to lower human-related mortality throughout adulthood. Against our expectations, the survival rate was higher for males. This was due to a higher overall number of male apparent survival events, which consist of the number of consecutive detections over years including gaps during which the animal was alive but eluded detection. However, we did not nd any strong difference in survival rates between sexes. A similar open population SCR study, conducted on a low-density ocelot population in Belize, also showed no signi cant differences in sex-speci c survival, despite a survival rate of 0.86 for females and 0.78 for males (Satter et al. 2019). The authors suggested the statistical power was not enough to detect signi cant differences between sexes in these parameters, although they were able to determine sex for a large number of adult ocelots (n = 322).
Camera trapping does not allow assessment of the fates of all disappearing individuals. While our ndings show increasing vehicle collisions in recent years (Table 1), we could not assess the actual impact of illegal killing since carcasses are seldom found and it only occurs outside the protected area, as con rmed by the high survival rates found inside it. Illegal killing still represents one of the main lynx mortality causes Europewide, responsible for 32% of known lynx death in Switzerland 5 , 46% in Scandinavia 45 and 71% in Poland 44 . In the vicinity of our study area, a high poaching rate has been suggested by a modelling approach 4 .
We did not detect a signi cant difference in per capita recruitment between sexes (Table 2) 47 . In Switzerland, the average litter size was 2.00 5 . We were not able to assess the proportion of reproductive females because sex was not determined for all detected individuals.
The maximum age observed in the study area was 10 years, which was constrained by the study duration.
Age cannot always be determined through camera trapping. As such, individuals already classi ed as independent in the rst session reached higher values than those we observed, even if their exact age remained unknown.
According to the RAI results, roe deer abundance on the German side remained stable despite the ban on population control of this species in 2012. Contrastingly, the red deer population increased (Fig. 3). This is consistent with the top-down limitation of roe deer abundance by lynx in the study site (e.g. Heurich et al., 2012), with little in uence on that of red deer. This nding might indicate that a higher consumption rate of roe deer was not possible in the area. However, especially male lynx might take advantage of increasing red deer numbers as they hunt calves and yearlings of that species 26 . The observed abundance decrease in red fox could be due to intraguild predation as it coincides with increasing in lynx numbers, as reported in Scandinavia (e.g. Elmhagen and Rushton, 2007). However, considering that in the study area the red fox is only 1% of lynx kills 26 , other factors such as food availability and diseases are likely driving red fox dynamics. Although RAI can be affected by biases attributed to changes in detection probability 50 , we believe this method well represented the abundance trend of the species in question since the study design was uniform throughout the monitoring period.
For lynx, we found positive population growth rates, high survival rates and low per capita recruitment resulting in a low yearly population turnover. A stable number of reproducing females and stable prey base indicate carrying capacity for lynx in the study area has already been reached. Therefore, the increasing number of independent lynx most likely comprises dispersing individuals crossing the study area, a few of which settle replacing disappearing resident individuals. This indicates that lynx distribution has changed rather outside the study area.
Contemporary assessments of the BBA lynx population, based on coordinated transboundary camera trapping within an area of 13,000 km2, show a slightly positive population trend and con rm lynx presence and reproduction in the surroundings of the study area 24 . In conclusion, the changes we observed could be related to locally decreasing poaching possibly as a result of both law enforcement and intensive long-term campaigns aiming to raise acceptance in the wake of high-pro le poaching incidents.
Our study stresses open population SCR models can provide a wider range of demographic parameters useful for lynx monitoring than closed population SCR models. Data gaps affecting the earlier sessions possibly resulted in biased estimates, especially regarding those of closed population SCR models in a maximum likelihood form (Supplementary Table S2). Bayesian methods can better deal with the problem of incomplete detection because they produce posterior distributions of the demographic parameters incorporating the uncertainty resulting from data gaps by using the information deriving from other primary periods 11 . This highlights the robustness of this method, which we recommend in future studies and monitoring, especially with incomplete detections. Nonetheless, closed population SCR models still represent a reliable tool for abundance and density assessments (e.g. Gimenez et al., 2019).

Conclusions
The protected areas in focus appear to host an important population core with high density and survival rate and low yearly turnover. As proposed for other large carnivores 15 , protected areas should be considered as strategic for lynx recovery across its range. However, our results clearly show that the suitability of the surrounding landscapes is fundamental for long term population viability, as the constant number of reproducing females throughout the monitoring period indicates lynx reached carrying capacity in the study area of two adjacent National Parks. This view is also supported by a stable prey base. Additionally, conservation of species with large spatial requirements cannot rely on protected areas alone especially when considering long term genetic viability 51 .
Our study revealed important demographic parameters of lynx that can be used for population forecasting across broad spatial contexts, thus improving conservation and management plans. Finally, we stress the importance of long-term systematic monitoring as a basis for the understanding of population dynamics of this and other large carnivore populations.

Con ict of interests
The authors declare that they have no known competing nancial interests or personal relationships that could have appeared to in uence the work reported in this paper.
Authors contribution