Hare's affairs: Lessons learnt from a noninvasive genetic monitoring for tracking mountain hare individuals

Abstract Systematic monitoring of individuals and their abundance over time has become an important tool to provide information for conservation. For genetic monitoring studies, noninvasive sampling has emerged as a valuable approach, particularly so for elusive or rare animals. Here, we present the 5‐year results of an ongoing noninvasive genetic monitoring of mountain hares (Lepus timidus) in a protected area in the Swiss Alps. We used nuclear microsatellites and a sex marker to identify individuals and assign species to noninvasively collected feces samples. Through including a marker for sex identification, we were able to assess sex ratio changes and sex‐specific demographic parameters over time. Male abundance in the area showed high fluctuations and apparent survival for males was lower than for females. Generally, males and females showed only little temporary migration into and out of the study area. Additionally, using genotyped tissue samples from mountain hares, European hares (Lepus europaeus) and their hybrids, we were able to provide evidence for the first occurrence of a European hare in the study area at an elevation of 2,300 m a.s.l. in spring 2016. For future monitoring studies, we suggest to include complementary analysis methods to reliably infer species identities of the individuals analyzed and, thus, not only monitor mountain hare individual abundance, but also assess the potential threats given through competitive exclusion by and hybridization with the European hare.

& Leberg, 2000). The multi-tube approach has been widely accepted to minimize these errors in NiG studies, but it also highlights the importance of study-specific error estimations as well as species-specific optimization of sampling strategies (Taberlet et al., 1996).
Here, we present a 5-year NiG sampling dataset of the mountain hare (Lepus timidus) as study species, used to assess demographic parameters in a strictly protected area which allows their estimation under assumedly natural conditions (Rehnus, Marconi, Hackländer, & Filli, 2013). This enables us to draw conclusions about the conservation status and types of threat present in an exemplary alpine habitat. The mountain hare is classified as "Least Concern" by the International Union for the Conservation of Nature (Smith & Johnston, 2008), signifying sufficiently large global population sizes and no detectable or foreseeable significant decline in population abundance (IUCN, 2012). However, an increase in elevation and shift poleward has recently been predicted for mountain hare populations due to climate change (Leach, Kelly, Cameron, Montgomery, & Reid, 2015). For populations in the Alps, Rehnus, Bollmann, Schmatz, Hackländer, and Braunisch (2018) predicted a decline of suitable habitat and an increase in fragmentation due to changing climatic conditions. Besides having direct consequences for the mountain hare, climate change may also have indirect effects on mountain hare populations due to increasing competition with the European hare (Lepus europaeus).
An extension of the range of the European hare into higher areas in the Alps can be expected from the positive effect of temperature increase (Thulin, 2003). As a consequence, the smaller mountain hare could be competitively excluded by the European hare (Acevedo, Jimenez-Valverde, Melo-Ferreira, Real, & Alves, 2012).
Furthermore, hybridization between the two species could lead to introgression and threaten the genetic integrity of the mountain hare (Thulin, 2003;Thulin, Fang, & Averianov, 2006). Recent hybridization has been well documented in Scandinavia (e.g., Levänen, Thulin, Spong, & Pohjoismaki, 2018), whereas the current patterns of hybridization in the Alps are largely unknown (Beugin et al., 2017;Zachos, Slimen, Hacklander, Giacometti, & Suchentrunk, 2010). To understand the threats mountain hare populations are facing in the Alps, it is thus essential to monitor the occurrence of European hares and hybrids, in particular where species' elevational ranges overlap.
In this study, we investigate how NiG can be used to track population changes and to quantify threats for mountain hares by considering sex-specific demographic parameters and tracking possible occurrence of European hares. In this way, we used the samples of an ongoing monitoring from the years 2014-2018 and asked the following research questions: (a) Do apparent survival rates differ between sexes? As activity rates are sex-specific during mating season and have been closely linked to survival (Murray, 2002), we expect to observe sex-specific apparent survival rates. (b) Are temporary migration rates equal between sexes? Both sexes have been shown to display high site fidelity (Bisi et al., 2011); therefore, we expect to observe low temporary migration rates for both sexes. (c) Are European hares already present in the study area? Through genotyping tissue samples of known European and mountain hares and comparing the thereby obtained individual genotypes to individual genotypes obtained for the monitoring, we expected to be able to assign species identities (Beugin et al., 2017).

| Study area
Our study area is located in the Swiss National Park (46°39′N, 10°11′E), spanning an area of 3.5 km 2 ; selected both to study mountain hare under natural conditions and to provide accessibility for sampling in rugged, hazardous alpine terrain in late winter (Figure 1, Rehnus & Bollmann, 2016). The Swiss National Park (IUCN Category 1a; IUCN, 2019) is inaccessible for the public in winter, but a popular area for recreational activities in summer, whereby visitors are restricted to use marked paths and areas.

| Sample collection
We collected fresh fecal pellets during the mating period (end of March until first half of April) to detect individuals that potentially contributed to the next reproductive cycle (Luikart et al., 2010;Thulin, 2003) and during the postreproductive period (October) for recording offspring born in the summer of the same year in addition to adults (Pehrson & Lindlöf, 1984). Only fresh fecal pellets, which were identified based on surface characteristics, were collected, as Rehnus and Bollmann (2016) found significantly lower amplification success rates for pellets older than 5 days. Samples were collected both systematically and opportunistically, following the methods of Rehnus and Bollmann (2016). Systematic samples were collected at 91 sampling points on a 200 m square grid (Figure 1). Opportunistic samples were collected whenever fresh fecal pellets were found between systematic sampling points. To ensure that only fresh fecal pellets were collected, all present fecal pellets were cleared from every systematic sampling plot during a first survey and any fecal pellets detected during a subsequent survey (3-5 days later) were collected. It was assumed that feces with different shades of surface color found at the same location originated from different individuals or different dates, as previous studies have shown a considerable overlap in home ranges of mountain hare (Bisi et al., 2011;. In these cases, we collected one feces per color type. In cases where two feces from the same color type were available, a backup fecal pellet was collected from each feces location to have a potential replicate for samples with low DNA quality (see DNA extraction and PCR amplification). All cases for which both the replicate and original sample were genotyped showed identical genotypes. Both systematic and opportunistic sampling were conducted within a short time span (maximally 11 days) in each sampling period, during which the population could be assumed as closed (Rehnus & Bollmann, 2016). Samples were collected and stored in separate plastic tubes without touching by hand to minimize DNA contamination (Sloan, Sunnucks, Alpers, Behregaray, & Taylor, 2000). After collection in the field, samples were frozen and stored until analysis in the laboratory.
For the possible identification of European hares, 90 tissue samples from known European and mountain hares and their hybrids were collected by hunters, game keepers, taxidermists, and other researchers at different locations across Europe ( Figure 1) and frozen at −20°C until further analysis. Samples included different types of tissue, such as parts of the ear, paws, muscle, and bones. The species of each sample was determined by its collector based on morphological characteristics and classified-if possible-as either mountain hare (N Lt = 51), European hare (N Le = 36), hybrid of both hare species (N Hb = 2), or else was classified as unknown species identity (N NA = 1).

F I G U R E 1
Origins of tissue samples (left, dots) and location and extent of the study area (right, red) in the Swiss National Park, with systematic sampling points (right, orange circles). On the left, origins of European hare (Lepus europaeus, Le) samples are displayed in blue, mountain hare (Lepus timidus, Lt) samples in dark gray, hybrids (Hb) in black, and samples of unknown species identity in light gray

| Noninvasive genetic samples
DNA from fecal pellets collected in 2014 and 2015 was extracted after every sampling period following the protocol described by Rehnus and Bollmann (2016). DNA extraction from fecal pellets collected in 2016-2018 was performed with reagents from a customized sbeadex livestock kit (LGC Genomics) on a King Fisher Flex (Thermo Fisher Scientific). DNA samples were amplified as described in Rehnus and Bollmann (2016) in two multiplex PCRs with three independent replicates, following a modified multi-tube approach (Taberlet et al., 1997(Taberlet et al., , 1999. For amplification of samples from 2016 to 2018, concentrations of primers were lowered to 0.2-0.3 µM. If initial PCR results were negative, the backup fecal sample was used as replacement. Fragment length analysis of samples from 2016 to 2018 was performed on an ABI3130 genetic analyzer using GeneScan LIZ 500 Size Standard (Thermo Fisher Scientific), and electropherograms were visually analyzed using GeneMapper v5.0 (Thermo Fisher Scientific) after each sampling session.

| Tissue samples
DNA extraction for tissue samples was done using the DNeasy Blood & Tissue Kit (Qiagen) following the protocol for purification of total DNA from animal tissues (Spin-Column Protocol, Qiagen).
Pieces of 25 mg consisting of muscle, tendon, bone, or some hair were cut off from the samples and stored in 2 ml tubes at −20°C until further analysis. DNA concentration was assessed using Quantus Based on the concentrations assessed using Quantus (Promega), DNA samples were diluted to 2.5 ng/μl. Genotyping of tissue samples followed the protocol described above for NiG samples.

| Genotyping
For tissue samples, PCRs not showing clear peaks in GeneMapper were repeated. For NiG samples for which multiple PCRs were run, consensus genotypes were created following the rules described in Table 1. The sex of tissue and NiG samples was assigned as male, if one or more samples amplified at the SRY locus, and as female, if no sample amplified.
The three loci Sat2, Sat12, and Lsa2 could not be reliably scored as biallelic markers. Thus, they were not included in the automized multilocus genotype analysis. However, loci Sat2 and Sat12, but not Lsa2, showed consistent patterns of multiple allele peaks and stutter peaks among replicates and were thus scored qualitatively for NiG samples and used to distinguish among replicate genotypes with mismatching loci (see below). Multilocus consensus genotypes were accepted if at least six of the seven remaining loci could be scored successfully.
The raw genotype table output of the remaining seven loci was analyzed to find consensus genotypes for each replicated sample.
Using a custom script in R v3.5.1 (R Core Team, 2018), we considered each replicate for which a genotype could be scored (positive PCR; Broquet & Petit, 2004) and applied the predefined conditions (Table 1): (a) Consensus homozygote genotypes were accepted if all three replicates were consistent, and (b) consensus heterozygote genotypes were accepted if at least two replicates were consistent and no more than two alleles were found across all three replicates.
Even though we applied stringent conditions, multilocus genotypes cannot be considered error free (Broquet, Menard, & Petit, 2007). Here, we quantified error rates as (a) missing replicate genotypes and (b) false homozygote replicates in consensus heterozygote genotypes. The rate of missing replicates was calculated in R as the proportion of missing replicate genotypes per sample in the Further types of quantified error included errors in sex determination. Sex determination using SRY is based on the amplification of the marker on the Y chromosome (Wallner et al., 2001). However, the nonamplification could either indicate a female individual or a dropout of the locus, and thus, certain errors are expected especially for female individuals.

| Genotype data
For both types of errors (false homozygote replicates and missing replicate genotypes), means for each season, year, and marker were calculated and compared using the package aGricolae in R, which applies a pairwise comparison of group means including standard errors (SE) based on Tukey's honestly significant difference (Mendiburu, 2019). For the proportion of false homozygotes, we additionally assessed the amount of missing data in relation to the sex of the individual/genotype. Because the assessment of missing data was estimated based on raw genotype data and the sex could not be reliably determined for samples with large amounts of missing data, the amount of missing data could not be checked for differences between sexes.
Unique multilocus genotypes and unique genotype groups (individuals) in the NiG data set were identified based on successfully scored multilocus consensus genotypes using the alleleMatch package in R, which considers genotyping errors and missing data during the assignment of individuals (Galpern, Manseau, Hettinga, Smith, & Wilson, 2012). Subsequently, genotypes showing mismatches, female genotypes in male genotype groups, and multi-match samples were double-checked using the two previously excluded loci (Sat2, Sat12). Checking the additional loci was done qualitatively through comparing peak patterns in GeneMapper. Samples showing identical peak patterns (fragment lengths, patterns of stutter peaks) in GeneMapper were considered to originate from the same individual. Some examples of such peak patterns are given in Appendix S1. Samples were considered to belong to the same unique group of multilocus genotypes if they did not show more than two mismatched alleles including the additional loci. Females showing an exact match to male genotype groups were reclassified as males and marked as false females.

| NiG monitoring
The calculation of population genetic statistics and testing for Hardy-Weinberg equilibrium (Hardy, 1908;Weinberg, 1908) for the NiG genotype data was done in cervus v3.0.7 (Kalinowski, Taper, & Marshall, 2007). We calculated observed and expected heterozygosity, polymorphic information content (PIC), null allele frequency estimates, and the probability of identity under the assumption that individuals are unrelated (P ID ) as well as under the assumption of siblings being present in the data (P IDsib ). Loci for which potential null alleles were estimated to be present were still included in the analysis, as no true proof of null alleles can be given and estimates of null allele frequencies in cervus are based on homozygote frequencies (Kalinowski et al., 2007).
CMR methods seek to estimate animal abundance via capture of individuals (e.g., collection of genetic material) and marking (e.g., genetic identification) for future identifications (Williams, Nichols, & Conroy, 2002). The robust design (Pollock, 1982) distinguishes between primary and secondary sampling occasions, whereby primary sampling occasions (PSOs) are separated by enough time for population change to happen, and consist of one or more secondary sampling occasions (SSO), during which the population can be assumed as closed (Kendall, Nichols, & Hines, 1997). We applied the robust design model to estimate individual abundance, accounting for misidentifications in the data, under the assumption that genotyping errors could be present (Lukacs & Burnham, 2005a;Otis, Burnham, White, & Anderson, 1978), as implemented in Mark v9.0 (White & Burnham, 1999). Each sampling day was considered as SSO and each sampling session (e.g., spring 2014) was considered as PSO.
Using individual detection histories based on (

| Species assignment
Species assignment of individuals detected in the Swiss National Park based on NiG sampling were assessed through including a tissue data set consisting of high-quality DNA samples of predefined species identity. To do so, a first analysis was conducted using only the tissue data set to identify whether differences between species could be detected with the loci used in the monitoring study.
Thereafter, the second analysis was done using a combined data set consisting of both the tissue genotypes and the individual genotypes detected in the Swiss National Park from NiG sampling. First, a principal component analysis (PCA) was conducted in the R environment on a genind object in aDeGenet (Jombart, 2008), which is based on the package aDe4 (Dray & Dufour, 2007), and the results were visualized using Factoextra (Kassambara & Mundt, 2017). Assignments made by the PCA were then analyzed to look for individuals from the NiG dataset with higher similarity to European than to mountain hares. Second, to get further confirmation for species affiliations, individual assignments were assessed using structure v2.3 (Pritchard, Stephens, & Donnelly, 2000). The admixture model was applied for both data inputs, using a burn-in period of 100'000 and 1'000'000 Markov Chain Monte Carlo (MCMC) repetitions after burn-in. Additionally, the correlated allele frequency model was used, which provides greater power to detect distinct but closely related populations (Porras-Hurtado et al., 2013). Values for K were set to 1-5, using ten iterations per K. The value of K best explaining the data was assessed using structure harvester v0.6.94 (Earl & vonHoldt, 2012). After estimating the value for K with the highest likelihood, the results for the runs of the respective K were combined using the Full Search Algorithm in cluMpp v1.1.2 (Jakobsson & Rosenberg, 2007) and plotted in R. Cluster assignments were then visually compared to groupings observed using PCA to find correspondences or discrepancies.

| Genotyping and distribution of errors
Over the five study years, 1,540 fecal samples were genotyped, while the remaining 97 samples failed due to negative initial PCR results.
Of the former, 316 samples (20.5%) had to be removed due to too many missing loci in the consensus multilocus genotypes. The amount of missing data was significantly larger for samples collected in fall (P NA = 0.24) than for samples collected in spring (P NA = 0.07, p < .001), On average, 4.4 alleles per locus were detected for individuals in the study area (Table 2). For loci Lsa3 and Sat5, the estimated null allele frequency indicates the possible presence of null alleles and an excess of homozygotes, respectively (Table 2). At locus Lsa3, allele 210 was by far the most common (P 210 = 0.78) and most genotypes containing this allele were homozygous (87.5%).

| Individual identification
From the 1,224 retained multilocus genotypes, 96 unique individuals were identified, of which 37 were classified as female and 59 as male.
More samples originated from opportunistic sampling (N opp = 992) than from systematic sampling (N sys = 232  (Figure 3). Consequently, both sampling methods revealed the same relative differences between sexes and season, but systematic sampling identified fewer individuals and individuals were identified based on fewer samples compared to opportunistic sampling.

| Model comparison
Of the 31 CMR models run, 21 models were unsupported, with model weights of less than 0.01. The ranking of the models did not show clear tendencies, as the four best models all showed a model likelihood of more than 0.25. Additionally, the model weight was mainly distributed among the 15 best fitting models, whereby the ten highest-ranking models accounted for almost all the model weight (w 1 -w 10 = 0.96, Table 3). For parameter estimation, we considered Models 1-5, as these models showed a model likelihood of more than 0.2, an AICc weight of more than 0.05, and a cumulative model weight of 0.85.
The model with the best fit (Model 1) described unequal capture and recapture probabilities, variation for φ by sex, variation for α, f0, p, and c by season, and constant values for γ″ and γ′. However, the model with the second best fit (Model 2) estimated α, p, and c with variation by sex and season, φ and f0 with variation by sex and constant values for γ″ and γ′ (Table 3). Models considering equal capture and recapture probabilities generally showed a lower ranking than

| Parameter estimation
Models 1-5 described survival to be constant through time; thus, no seasonal or annual survival rates were estimated. Based on these models, females showed a higher apparent survival probability than males (Table 4).
In spring, recapture probabilities were higher than capture probabilities and both were higher than in fall for both sexes. For males, the difference between seasons was more pronounced than for females ( Figure 4), as both capture and recapture probabilities were higher in spring than in fall ( Figure 4, Table 4). Consequently, more male genotypes were observed in spring compared to fall for both newly identified and previously identified individuals. For females, capture probabilities were equal between seasons and recapture probabilities showed a less pronounced difference between seasons than for males ( Figure 4, Table 4). Additionally, recapture probabilities were higher than capture probabilities in spring and equal to capture probabilities in fall (Table 4). Therefore, the probability to observe a new female was approximately equal in both seasons, but the probability to observe a new male was higher in spring than in fall.
The probability of temporary emigration between two subsequent PSOs (γ″) was estimated to be 0.05 ± 0.04, and the return rate of temporary emigrants (γ′) was estimated to be 0.33 ± 0.30 for both sexes (Table 4). Additionally, the no-movement model was ranked relatively high (Model 3, Table 3). Consequently, temporary immigration and emigration (migration) into and out of the study area are assumed to occur at a similarly low rate for females and males.
In both seasons, the probability of correctly identifying an individual (α) was larger for females than males. Females showed a higher α in spring (0.95 ± 0.05) than in fall (0.81 ± 0.10), and α was the same for males in both seasons (0.77 ± 0.11). Model 1 estimated f0 to be variable only by season and estimated more individuals to be undetected in fall than in spring. Model 2 estimated f0 to be variable by sex and estimated more males to have remained undetected per primary sampling occasion than females (Table 4).

| Systematic sampling
The models set up using capture histories based on only systematically collected data had relatively large standard errors for capture probabilities (Figure 4a), but showed similar patterns for recapture probabilities (Figure 4b). Capture and recapture probabilities were lower compared to the probabilities calculated by the models using capture histories based on both systematic and opportunistic samples ( Figure 4). Thus, seasonal and sex-specific patterns of recapture probabilities were confirmed using systematic sampling.

TA B L E 4
Parameter values as weighted means across capture-markrecapture Models 1-5 for mountain hare (Lepus timidus), standard error and 95% confidence intervals (95% CI) as weighted mean across models (capture rate (p), recapture rate (c), rate of misidentification (α), number of missed individuals (f0), rate of temporary emigration (γ′) and immigration (γ″), and apparent survival rate (φ)) or as given by Mark (White & Burnham, 1999) for the specific model (f0) Additionally, fall 2017 was the only sampling session for which more females than males were estimated to be present in the study area ( Figure 5).

| Individual abundance
On average, a density of 6.4 hares per km 2 was observed. The maximum density in the study area was 8.3 hares per km 2 and occured in spring 2017. In fall 2017, the density was minimal with 4.9 hares per km 2 . One individual (LtNP02, female) was detected in every PSO from spring 2014 to fall 2018.

| Species assignments
The PCA based on tissue sample genotypes with predefined species identities showed that most samples of the data set were genetically similar to each other at the loci assessed (Appendix S2).
Nevertheless, the representation of the PCA depicted a lower left shift for most L. timidus samples and an upper right shift for most L. europaeus samples, also for the combined dataset consisting of NiG and tissue samples. Hereby, the two most informative axes accounted for 19.3% of the total variation ( Figure 6). Some of the samples contained in the tissue dataset showed different species assignments to the ones defined by their collectors (e.g., Lt107, Figure 6). The presumed hybrid samples contained in the F I G U R E 4 Overview of the weighted mean of capture (a) and recapture probabilities (b) calculated for the models using capture histories based on systematically collected samples and the whole dataset (systematic and opportunistic), given for male and female mountain hare (Lepus timidus) individuals. Error bars indicate weighted means of standard errors calculated by Mark (White & Burnham, 1999) tissue dataset were shown to be either more similar to mountain hare (Hb597) or to European hare (Hb604). The combined dataset further showed two samples from the NiG monitoring dataset (LtNP60, LtNP95) to be more similar to European hare than mountain hare samples ( Figure 6). structure harvester (Earl & vonHoldt, 2012) revealed the most likely number of clusters for the combined and for only the tissue genotype data to be K = 2 (Figure 7), which is expected in a two-species case. Generally, species assignments obtained for the tissue genotypes were mostly in accordance with the results obtained by the PCA. However, individual LtNP95 only showed slight signs of admixture, even though its position in the PCA plot was in-between samples of the two hare species. Unexpectedly, individual Lt107 from the tissue dataset, morphologically identified as L. timidus, was clearly confirmed to belong to the L. europaeus cluster by both structure (P LE = 0.995) and PCA (Figures 6 and 7).
Further, the assignment of LtNP60 to the L. europaeus cluster was confirmed by structure at a high probability (P Le = 0.98). Notably, individual LtNP60 was genotyped based on ten samples collected in spring 2016 at a maximum elevation of 2,300 m a.s.l. in the study area.

| D ISCUSS I ON
In this study, we show the results of a multi-year NiG monitoring on the elusive mountain hare. Through genotyping nuclear microsatellite loci and a sex marker based on DNA extracted from fecal pellets, we were able to show that male mountain hare individuals not only had lower apparent survival rates than females, but also higher fluctuations in activity rates throughout the year, presumably based on their mating behavior. We also found equally low temporary migration rates for males and females, which can be explained by the high site fidelity of both sexes. Further, we were able to show that a European hare had been present in the Alpine habitat at a maximal elevation of 2,300 m a.s.l., which marks the first observation of this species in the study area.

| Density of individuals
We found that density of individuals in the Swiss National Park was on average nearly twice as large (6.4 hares per km 2 ) as previously estimated (Rehnus & Bollmann, 2016). Rehnus and Bollmann (2016) F I G U R E 5 Individual abundance of mountain hare (Lepus timidus) in the study area in each primary sampling period, as given by each of the five best fitting models in Mark (White & Burnham, 1999) and as their weighted average. Error bars indicate standard errors as estimated by Mark for the model used samples from one sampling session (spring 2014) and derived spatial capture-recapture estimates, which consider a buffer around the study area (Efford & Fewster, 2013). However, the density values previously estimated in this study area and in a close-by location in Italy (5-11 hares per km 2 ; Nodari, 2006) are also higher than those observed in a National Park in Austria (0.4-0.7 hares per km 2 ; Slotta-Bachmayr, 1998) and in Ticino, Switzerland (3.5 hares per km 2 ; Gamboni et al., 2008). Besides each of these studies applying different methods and being conducted during different seasons, climatic conditions differ throughout the biogeographic regions of the Alps (Rehnus et al., 2018). Hence, comparisons of such density estimates should be taken with caution.

| Fluctuations in census estimates
Seasonal fluctuations of the estimated census size were mostly due to fluctuating abundance of male individuals, as more individuals were present during the mating season than in the postreproductive period. Apparent survival was lower for males than females, male individuals mostly disappeared between spring and fall, and new male individuals appeared between fall and spring. Equally low amounts of temporary migration for both sexes confirmed the high site fidelity of mountain hares (Bisi et al., 2011). We assumed that this pattern can be explained by the range use of individuals, which often depends on the mating system (Ostfeld, 1990). The distribution of males is mainly determined by the availability of potential mates, whereas range use of females is mostly affected by the availability of food and shelter for nursing offspring (Bisi et al., 2011). Our spring sampling sessions coincided with the beginning of the mating season (April; Thulin & Flux, 2003), during which males are searching for mates. During this time, home ranges have been found to be enlarged compared to fall (September-November) for both sexes and to be larger for males than females (Sweden, Dahl & Willebrand, 2005).
This was also found to be the case in the Alps, and as a consequence, home ranges also overlapped more strongly during this time, leading to more individuals being present (Bisi et al., 2011;Gamboni et al., 2008). We therefore conclude that males increase their home ranges during mating season, leading to more males being detected in the study area during this time of the year. Further, we conclude F I G U R E 6 Results of the principal component analysis for genotypes of noninvasive genetic samples (light gray circles) and of tissue samples (triangles). Some samples morphologically assigned to Lepus timidus (Lt, dark gray) showed signs of being assigned to L. europaeus (Le, blue) and vice versa. One hybrid (Hb, black) sample was assigned to L. europaeus, and one hybrid sample was assigned to L. timidus. The sample of unknown origin (Unknown, light gray) was assigned to L. timidus. The large symbols represent the centroids of the defined groups that the survival of males is lower than of females, leading to sexspecific fluctuations in abundance.

| Effects of sex and season on recapture rates
Our results show that individuals of both sexes had higher recapture probabilities during mating season than during the postreproductive period, whereby recapture probabilities were higher for males than females. From an ecological perspective, we explain higher recapture rates by higher activity rates. For example, hares are assumed to be more active during the mating season (Bisi et al., 2011), and male snowshoe hares (L. americanus) have also been found to increase their activity during this time of year (Murray, 2002). Higher activity of hares has been closely linked to lower survival probabilities via increased predation risk (Murray, 2002). Bisi et al. (2011) hypothesized that the reduction in home range sizes in fall, could be an antipredatory strategy, that is, to reduce predation risk while being midmolt.
From a methodical perspective, differences in genotyping success rates may cause differences in capture and recapture probabilities, as only successfully genotyped samples are recorded as capture occasions (Lukacs & Burnham, 2005a). We found genotyping success rates to be constant between sexes, but higher in spring than in fall. To control for a biased effect due to detectability and genotyping success caused by season-specific sampling conditions, individual histories based on systematically detected samples were analyzed. Due to the sampling scheme, detectability is thought to only minimally influence the number of samples collected using systematic sampling. Based on these models, differences in capture and recapture probabilities between sexes and seasons were confirmed. This indicates that seasonal and sex-specific variance is largely due to individual activities and reflects ecological conditions with higher activity rates and larger home ranges of individuals-especially of males-during mating season. This difference leads to higher individual abundance and a higher number of samples collected per individual, whereas a reduction in home range sizes in fall leads to lower recapture rates.

| Monitoring European hare occurrence with implications for hybridization
One individual detected in our study area was identified as European hare. The individual (male, LtNP60) was found to be present in the park only during spring 2016 at a maximum elevation of 2,300 m a. s. l. Species assignment of this individual was confirmed using mtDNA sequencing (CoI; S. Brodbeck, unpublished data). mtDNA sequences differ between mountain and European hares and can thus be used for reliable species identification (Thulin, 2003;Thulin, Stone, Tegelstrom, & Walker, 2006).
With climate change, habitat of European hares is thought to experience an upward shift, which may cause European and mountain hare habitat to increase in overlap (Acevedo et al., 2012;Rehnus et al., 2018), intensifying competition between the two species and causing both competitive exclusion of the mountain hare and potential hybridization (Thulin, 2003). Here, we detected one individual of European hare at its upper end of the current elevational distribution range (Bauer, 2001;Bisi, Wauters, Preatoni, & Martinoli, 2015). While at least one further individual showed signs of admixture, we consider a better resolution of markers to be necessary to clearly detect admixture between both species. With our marker set, structure and PCA did not always show consistent results, but also morphologically based identification proved potentially inconclusive. To address this phenomenon of hybridization more reliably, genomic methods are thought to provide better resolution (Allendorf, Leary, Spruell, & Wenburg, 2001;Carroll et al., 2018). We recommend future work to develop a diagnostic panel of single-nucleotide polymorphisms (SNPs), including mtDNA sequences for species diagnostic purposes (Thulin, Fang, et al., 2006), for integration into the standardized future F I G U R E 7 Results of the analysis run in structure (Pritchard et al., 2000): Cluster assignment probabilities are given in gray for Lepus timidus and in blue for L. europaeus. Individual genotypes identified based on noninvasive samples are named with the prefix "NP" and depicted on the left side. The morphological species assignments of tissue samples are given as "Le" for European hares, "Lt" for mountain hares, "Hb" for hybrids, and "NA" for unknown species identity monitoring to quickly assess species associations of the individuals detected.

| Caveats and limitations
While the value of 0.036 for the probability of identity under the assumption of the individuals in the sample being related (Waits et al., 2001) is not particularly low, it still lies below the recommendation of the upper threshold given by Woods et al. (1999).
The probability of identity (Waits & Leberg, 2000) is a proxy for the power of a set of markers and depends on the number of markers, the allele numbers and their frequency distributions, and the typing error rate of each marker (Wang, 2006). Even without error in the data, there is a probability that random individuals share a common genotype, which may lead to misidentifications of individuals (Mills, Citta, Lair, Schwartz, & Tallmon, 2000). However, the identification of unique genotypes was done with two additional loci than those included in the calculation of the P IDsib value, for which electrophoresis peak patterns were consistent but did not allow biallelic scoring. Integrating this information in the calculation would lower the overall P IDsib value (Wang, 2006). Including the two additional loci showed to be highly relevant. A few genotypes were identical at the applied seven loci, but showed distinctively different peak patterns at the additional loci (see Appendix S1 for examples). Therefore, interpreting banding patterns of two additional loci allowed us to resolve certain genotype groups into distinct individuals.
In our study, error rates differed between seasons and declined throughout the study years, but were constant between sexes.
Throughout all years, we reduced errors by applying a modified multi-tube approach based on Taberlet et al. (1996), performing three PCRs per sample. Further, we applied a strict rule in detecting homozygotes to ensure a minimum allelic dropout rate and increase the probability of consensus homozygous genotypes to be true homozygotes (Taberlet et al., 1996). The dropout rate across loci and years of approximately 20% is in accordance with observed dropout rates in other NiG studies (e.g., Ebert, Sandrini, Spielberger, Thiele, & Hohmann, 2012;Hansen, Ben-David, & McDonald, 2008).While such restrictive genotype calling limited the overall genotyping success rate, it allowed us to keep the error rate at a reasonably low level.
Error reduction strategies and error quantifications were done consistently across years. We assume that error variation between seasons originates mostly from sampling conditions. In spring, samples were preserved by snow and low temperatures, whereas in fall, samples were less well preserved due to higher temperatures and sampling without snow. Further, Murphy, Waits, and Kendall (2003) observed that nutritional composition in bear feces causes variability in observed DNA qualities and quantities and therefore also in observed error rates. Seasonal differences in genotyping errors have also been observed for ungulates (Maudet, Luikart, Dubray, Von Hardenberg, & Taberlet, 2004). As food composition for mountain hares is different in fall than in spring (Rehnus et al., 2013), this aspect could be an additional factor causing the observed seasonal variability in error rates. We therefore conclude that genotypes of samples collected in fall, and consequently individual identifications, are less reliable than those of spring samples. Some errors were further observed during the determination of sex. Even after removing samples with too many missing loci in consensus genotypes, some samples classified as females showed matches to groups of male genotypes and were consequently marked as false females (0.9%). Further support for this classification is given by the fact that the amount of missing genotypes was higher for false female samples than for all other samples (P (NA, false females) = 0.06, P (NA, females) = 0.01). A certain error rate in sex determination using SRY was also detected by Wallner et al. (2001).
This finding highlights the importance of replicating noninvasively collected samples and of using an adequate number of samples for individual identifications.

| CON CLUS IONS
In this study, we show how sex-specific demographic parameters of an elusive and difficult-to-spot boreo-alpine species can effectively be monitored using noninvasive genetic methods and we highlight potential pitfalls thereof. By detecting a European hare at a high elevation, we emphasize the importance of integrating hybrid identification into future mountain hare monitoring projects, especially in areas with sympatric occurrence of both species. As climate change is expected to increase the ratio of sympatric occurrence, we consider the mountain hare to be an excellent model species for assessing how increasing ambient temperatures and a possible upward shift of a competitive, phylogenetically related lowland species may jeopardize rare high-elevation species that are pushed toward their upper range limits with restricted habitat availability.

ACK N OWLED G M ENTS
We thank the Swiss National Park for granting permission to carry out this study, all involved hunters, game keepers, taxidermists, and researchers for providing us tissue samples for the analysis and several laboratory apprentices who helped with preparing samples and generating data. We further acknowledge helpful comments from two anonymous reviewers and the Associate Editor. We are grateful to the Federal Office for the Environment FOEN, the Margarethe und Rudolf Gsell Foundation, Migros, and the Bristol Foundation for financial support.

CO N FLI C T O F I NTE R E S T S
We have no conflict of interest to declare.