# A Demo-Genetic Analysis of a Small Reintroduced Carnivore Population: The Otter (Lutra lutra) in The Netherlands

Authors: Maelle Seignobosc; Lia Hemerik; Hans Peter Koelewijn
Journal: International Journal of Ecology
 (2011)

---

## Abstract

Assessing the viability of reintroduced animal populations is a complicated task. Reintroductions are usually carried out with a small number of individuals, thereby, limiting the possibilities for monitoring because of the possible negative effects of intensive monitoring on survival and reproduction. Moreover, reintroduction studies are part of a socioeconomic interplay of forces, thereby, also limiting monitoring possibilities. Also, knowledge of population demography and abundance can be incomplete or unattainable. Here, we illustrate how we combined traditional telemetry and novel non-invasive genetic methodology to construct a detailed life table of a small reintroduced otter population in The Netherlands. Combining an appropriate capture-mark-recapture framework with a matrix modelling approach provides, in general, useful insights for such populations. The data indicated that (i) male survival is lower than female survival, (ii) the reintroduced population is currently growing (estimatedλ=1.26: range [1.06, 1.42]) and seems viable, (iii) increasing adult survival is currently the critical stage at which efforts of field managers should concentrate, and (iv) the modelling framework allowed us to determine the boundary conditions for the vital rates under which the population would go extinct. The applied approach directs at measurements that help field managers to implement the right conservation strategy after reintroductions.


---

## Body

## 1. Introduction

Reintroductions of extinct or nearly extinct species in a country are becoming nowadays a more accepted tool for the restoration of biodiversity [1]. Reintroduction attempts are, however, not undisputed [2]. On the one hand, several stakeholders are involved, each with their own priorities (e.g., government, field managers, scientists, and the general and local public) making the decision process surrounding reintroductions sometimes cumbersome. On the other hand, the number of released individuals is usually limited, and a new small population is created, with all the accompanying risks of extinction again [3]. Extensive monitoring of the released individuals for a longer time period, in order to derive guidelines for success or failure of attempts, is, therefore, of utmost importance [4]. Unfortunately, many reintroduction studies suffer from a lack of detailed observations for a longer time span [4–6], causing difficulties in explaining the success or failure of reintroduction attempts and preventing adequate advice to field managers in charge of reintroductions. Currently, most inferences in reintroduction biology are made by induction, gained from post hoc interpretation of monitoring results or by exploratory comparative analyses [4]. The main information needed for a long-term perspective on reintroduced populations is on the vital statistics, for example, mortality and recruitment. This entails extensive monitoring, which is an obstacle for many projects because of its expense [7]. Yet from the many anecdotal and contradictory results [4], it is clear that postrelease monitoring is badly needed for evaluating the success of reintroduction efforts.Elusive and nocturnal animals are usually hard to trace, and it is, therefore, difficult to collect population information, especially when densities are low, as is often the case with reintroduced animals. These animals could be provided with a radio transmitter, but these usually only last for limited time. An alternative is the use of noninvasive genetic procedures, that is, extracting DNA from biological samples such as hair, faeces, or urine without handling, capturing, or even observing the animals, since they offer the opportunity to collect detailed information on the population for a longer period [8]. Genetic monitoring, that is, quantifying temporal changes in population metrics or other population data generated using molecular markers, is becoming increasingly important for monitoring the consequences of anthropogenic change on wild species [9–11]. Combining genetic monitoring and noninvasive sampling into noninvasive genetic monitoring offers an excellent supplementary technique for studying the structure of reintroduced populations, as this enables time-dependent processes such as mortality and recruitment to be monitored, especially if DNA is sampled from the onset of the reintroduction, without disturbing the released animals [9]. Vital statistics can be calculated from the demographic database that originates from the observed DNA profiles and be used in population viability modelling.Noninvasive genetic approaches have been applied for many purposes, with the largest contributions from studies that focus on (i) the identification of individuals and the estimation of population abundance, (ii) the estimation of population genetic parameters (gene flow, demographic history, kinship), (iii) identifying species diets, and (iv) assessment of mating systems and behavioural ecology (cf. Table  1 in [12]). The integration of DNA-based observations in a capture-mark-recapture (CMR) framework [13], although being recognized as promising for years, has only just started [14–16]. CMR models provide estimates of survival rates that, when combined with information on reproduction, could be used in matrix models [17] to predict the future fate of the population and to identify those factors that influence population growth the most.The current study aims at integrating the results from a six-year noninvasive genetic study that accompanied the reintroduction of the otter (Lutra lutra) into The Netherlands [18] with a CMR and population viability modelling approach to get reliable estimates of the vital rates of the reintroduced population. This enables us to evaluate the viability of the newly founded population on the long term, while identifying the most critical aspects influencing yearly population growth, thereby, aiding field managers to take appropriate measures to ensure the sustainability of the new Dutch otter population.


## 2. Material and Methods

### 2.1. Biology of the Otter

The Eurasian otter (Lutra lutra) is a semiaquatic and mainly nocturnal mammal from the mustelid family [19]. The majority of its diet consists of fish and sea food. Sometimes frogs or small mammals are eaten. The preferred habitat is a dense vegetation situated near water offering the tranquility and favorable environment for reproduction.The otter is solitary, and its territory, especially those of males, can be large. Searching for food makes animals cover large distances, at the same time increasing their exposure to road dangers. A dominant male can settle over up to 80 km of rivers, overlapping the territory of several females [19]. Males and females rarely interfere, except for mating. Litter size is generally two or three young. Cubs can be born year-round and stay with the mother until they reach adult age between one and two year old. While an otter in captivity can live up to 15 years and the oldest one found in the wild was 16, the otter’s life expectancy is rather low, between four and five years. Within a community individuals communicate and mark their home ranges by using scented marks, “spraints”, usually left in prominent sites, near water, on tree stumps, under bridges or at the junction of roads. It is these scent marks that, when freshly collected, provide DNA from which a genetic fingerprint can be obtained.


### 2.2. The Otter in The Netherlands

The otter disappeared from The Netherlands in 1989. In 2002 a reintroduction programme was started with wild and captive otters, in order to start a new, self-sustainable, resident population (cf. [20]). From June 2002 until November 2008, a total of 31 animals, originating from Eastern and Northern Europe, have been released in a lowland peat marsh area in the north of the Netherlands (Appendix A), ca. 200 km2 in extent and geographically located at 52°38′–52°50′N, 5°53′–6°09′E. The area consists of a mosaic of peat grasslands, reed beds (Phragmites australis), and swamp woodland, mainly composed of Willow (Salix spp.), Black alder (Alnus glutinosa), and Birch species (Betula spp.) intersected by small ditches and pools. Most of the surrounding landscape is intensively farmed. The area is divided in three main sections: the Weerribben and the Rottige Meenthe solely consist of peat grassland, woodland and ditches, while in the Wieden also larger lakes are present (Figure 1).Figure 1
Geographical location of the study area. Within black borders the area where 31 otters (see AppendixA) were released between 2002 and 2008. The expanding population was monitored using telemetry and noninvasive genetic sampling.The released animals varied between 1 and 5 years old. Before release, tissue and blood samples were taken for DNA fingerprinting, all animals were tagged with a transponder, and a radio transmitter was implanted intraperitoneally (by a veterinarian at the Burgers Zoo (Arnhem, The Netherlands)). The study was conducted in accordance with Dutch legislation on the protection and welfare of vertebrate animals used for experimental and other scientific purposes.


### 2.3. Monitoring

Data collection in mark-recapture protocols requires individually marked (recognisable) animals. Here, identification was done according to two methods, without (re)capturing the individuals. Released otters were observed by either telemetry or noninvasive genetic sampling. Offspring (native borns) could only be observed by noninvasive genetic sampling [18]. It should be noted that the full population at the sample site consists only of released individuals and their offspring.Released otters were equipped with a radio transmitter, and their presence and movements could be recorded during the initial stages of the project. Each animal was searched at least once a week for almost one year, until transmitters failed. Based on the analysis of the blood samples taken before release, DNA fingerprinting could be performed by microsatellite analysis (AppendixB). Both released and native-born individuals from the release area were subsequently traced by applying noninvasive genetic methodology and using DNA as a tag [8, 9, 21].Fresh droppings consisting of scented marks (spraints) of the animals were collected at regular intervals during the winter periods (see for full details [18]). DNA was extracted from these spraints to build DNA profiles of the population (Appendix B). A profile could either reveal a new pattern or match a released otter’s pattern. In the latter case, it was considered a “resighting” in this individual’s encounter history. A new profile had to originate from a native-born otter, because no otters were present before the first one was released. For these native-born animals, the parents could be identified. For most animals several spraints were found during one year (range 1–27; [18]).In addition, dead animals were reported together with the cause of death (mainly traffic collision but also drowning, trapping in cages designed for muskrats and natural cause). These dead individuals were identified by their DNA fingerprint. Together, the three types of observations (telemetry, DNA profiles, dead individuals) allowed us to compile detailed encounter histories for the individuals. The multiple observations from both telemetry and DNA analysis are summarised in the event histories of the animals: every period of 6 months (from July 2002 to January 2009) is an observation interval in which they could either or not be seen alive (1/0) and either or not be recovered dead (1/0) (see Figure2). Per half-year interval, two places in the encounter history were filled as required for analysis with a joint model for live and dead data with the program MARK [22]. This implies that for the released animals both telemetry and genetic data were used, for the first two considered half-year periods per animal, and genetic data for the remainder of their lifespan; while native-born individuals only had encounter histories based on genetic data. In this way the rather small dataset of 81 animals (31 released and 50 native-born animals) and 13 encounter observations was constructed (each encounter history was a string of 26 zeros and ones, where at odd places information about resighting alive was summarized and at even place information on whether or not the individual was recovered dead).Figure 2
Fate diagram for the live-dead model. The convention used for the parameterization is the one developed by Seber [23], with S representing the probability of survival from one occasion to the next, p the probability of recapture in the following interval given the animal survived the interval, r the probability of being recovered dead, and F the site fidelity. The succession of events occurring during the period of observation is coded in an event history with information on alive or dead in the odd and even numbered places of the event history, respectively. For example, an individual marked and released on the first occasion, which survives the interval between occasion 1 and occasion 2 because it was recaptured at occasion 2, is subsequently recovered dead during the second interval. Its event history for these two intervals is 1011 (seen alive, not reported dead, seen alive, reported dead).


### 2.4. Statistical Analysis with RMark

Statistical analysis has been conducted with RMark [24], an R package using the program MARK [22] under the R interface.Data included both live resightings and recoveries of dead animals. These were combined into encounter histories as input for the Burnham model for joint data [25], allowing for increased precision in the estimation of the parameters compared to models treating live or dead data separately.For the parameters we followed the convention as developed by Seber [23], with S representing the probability of survival from one occasion to the next, p the probability of recapture in the following interval given the animal survived the interval, r the probability of being recovered dead, and F the site fidelity. Our focus is on survival as it is needed for the construction of a matrix model. The other estimations are not explicitly used, but allow accounting for all possible fates of the animals under the assumptions of the model and avoiding confusion between them. For example, if one does not know the probability that an animal stays in the sample area (F), the events of dying and leaving the area are both treated as mortality.For fitting the models to the data, we assumed for the parameters (a) the site fidelity is 1 since we are able to identify the animals anywhere in the area, even when they leave the release region, because there is only one (reintroduced) otter population in The Netherlands, (b) all parameters are time independent, because of the low number of animals and years involved, and (c)S and p can either be dependent on both sex and stage (defined by the age when released and time elapsed since release or birth) or sex or stage or be constant among sexes and stages.We built (st)age structure into the models, because in advance it is not obvious whether survival of differently aged young otters is different from adult survival. For that reason we started by considering three stages of the animals: juveniles (from birth to one year old), subadults (between one and two years of age), and adults (older than two years). This was possible because the age of the animal when released was known, and for the native borns we could estimate the most likely birth date based on the first encounter and the inferred ancestry of the parents. Native borns are at least ca. 6 months old before they start producing spraints. We also considered only two age classes (making a distinction for individuals less than one year old and older than one year) to allow for one (common) survival probability for subadults and adults. Fitting a model for every parameter combination resulted in fitting 32 different models (four possibilities of dependency forS(sex, stage, sex * stage, or constant), four for p(same as for S), and two for r(sex dependent or constant)).The first step in the modelling process is to find the most general and biologically relevant model (hereS(~sex * stage) p(~sex * stage) r(~sex) see Table 1) and to test for its quality of adjustment, stating that, if this adequately fits the data, then the nested models also do. Estimations for the Burnham model are computed in MARK by likelihood maximization assuming that data can be appropriately represented by a multinomial distribution. If variability is greater than what can be explained by model structure and the multinomial assumption, they are called overdispersed. This occurs when there is interdependence in the data or when the assumption of parameter homogeneity among individuals of the same group is violated. As a consequence variances of the parameters are underestimated, and model selection favours complicated models. These effects can be accounted for with a reliable estimate of overdispersion. We used the variance inflation factor c-hat(ĉ) [26]. If ĉ  is greater than 3, the data show a relatively important overdispersion [27]; lower values are not problematic.Table 1
Model selection (two-stage models) for survival (S), resighting probability (p), and recovery probability (r) with the site fidelity (F) of the otters set to 1. For each model the corrected quasi likelihood Akaike’s information criterion (QAICc), the difference in QAICc between the current model and the best fitting model (ΔQAICc), the number of estimable parameters (np), their weight (QIACc weight), and deviance (Qdev) are given. Fitted models are sorted by increasing QAICc, the first being the one that best represents the data. We present here only those models for which QAICc weight was greater than 0.01. These results are obtained after adjusting for overdispersion with ĉ=1.338.
FormulanpQAICcΔQAICcQAICc weightQDevS(~sex * stage) p(~stage) r(~1)7524.36700.478386.088S(~sex * stage) p(~stage) r(~sex)8525.6071.2400.257385.191S(~sex * stage) p(~sex * stage) r(~1)9527.0742.7060.123384.504**S(~sex * stage) p(~sex * stage) r(~sex)10528.2363.8690.069383.495S(~sex) p(~stage) r(~1)5531.7097.3420.012397.648S(~stage) p(~stage) r(~1)5532.0567.6890.010397.995**denotes the most general model and bold characters the two models used for estimating the parameter values (see Table2).Testing for goodness of fit is usually done with a parametric bootstrapping procedure [26]. The values of the parameters found for the model are used to generate a new dataset (simulated encounter history) according to the assumptions of the model, so that the simulated data present no overdispersion. For every simulation the deviance and other adjustment criteria are computed, sorted by increasing order, and compared to those of the tested model. The percentage of simulated deviances larger than the observed deviance gives us the probability that the tested model presents lack of fit. In addition, this method gives us a measure of the magnitude of the overdispersion of data, the “variance information criterion” or adjusted ĉ. It is obtained by dividing the observed ĉ and the mean ĉ from bootstraps. This adjusted value is subsequently used to have a better model fit [28]Usually, model selection of the fitted models is based on Akaike’s Information Criterion corrected for small sample size (AICc): the model with minimum AICc is the best [29]. When, however, there is evidence of overdispersion, the model selection is based on the Quasi likelihood Akaike Information Criterion (QAICc), proposed by Barker and White [28]. Let L be the likelihood of the model,ĉ the inflation criterion, n the number of observations, and k the number of parameters; we have(1)QAICc=-2log⁡⁡(L)ĉ+2k+2k(k+1)n-k-1.We also used QAICc weights which are the relative percentage of “support” from the data to a given model in the whole set of considered models. These weights sum to one. Models withwi<0.01 have a bad support from the data [26].


### 2.5. Matrix Model Development

The next step is to determine how the population evolves in the long term given the demographic parameters, reproduction and survival. We build a stage-structured Leslie matrix model for females only [17], considering that they drive the population dynamics as there are enough males to fertilize females [18]. Based on the results of the survival analysis, we develop a matrix model for either three or two stages (see Appendix C, (C.1)). We assume a 1 : 1 sex ratio at birth, density independence of the population growth rate as the carrying capacity of the environment is most probably not yet reached [18], and a birth-flow population (breeding season is continuous in the population, and female otters can have cubs all year round; [19]). We used a one-year time step, because we easily could obtain yearly estimates of survival and fertility parameters from our dataset.The fertility estimate, used in the matrix model, was based on the field observations during each year [18]. From the combined genetic and demographic database, it was possible to estimate (i) the number of identified offspring per successful mating, (ii) the total number of breeding pairs, and (iii) the number of adult females that did not reproduce. Fertility was calculated as the total number of new offspring divided by the total number of adult females present, both reproducing and nonreproducing. For the dominant eigenvalue we performed a bootstrapping procedure to assess the range in which the mean value of this yearly population growth factor falls. The 1,000 simulations are based upon intervals that span from the (mean – 2 * SE of the mean) to the (mean + 2 * SE of the mean) for the fertility (Ra=2fa), the juvenile (sj), and adult survival (sa).


### 2.6. Elasticity Analysis

Elasticity analysis is used to identify which matrix element influences population development most. Therewith, it is indicated which element should require closer attention in the perspective of improving conservation of the species. Elasticity of the yearly growth rateλ (the dominant eigenvalue of the matrix) to each of the demographic rates is the proportional change in λ when the parameter undergoes a very small variation, while other elements being held equal. Because matrix elements can be formulas of different survival and/or reproduction parameters, we also have to analyse the elasticity of these lower-level parameters (see Appendix C).


## 2.1. Biology of the Otter

The Eurasian otter (Lutra lutra) is a semiaquatic and mainly nocturnal mammal from the mustelid family [19]. The majority of its diet consists of fish and sea food. Sometimes frogs or small mammals are eaten. The preferred habitat is a dense vegetation situated near water offering the tranquility and favorable environment for reproduction.The otter is solitary, and its territory, especially those of males, can be large. Searching for food makes animals cover large distances, at the same time increasing their exposure to road dangers. A dominant male can settle over up to 80 km of rivers, overlapping the territory of several females [19]. Males and females rarely interfere, except for mating. Litter size is generally two or three young. Cubs can be born year-round and stay with the mother until they reach adult age between one and two year old. While an otter in captivity can live up to 15 years and the oldest one found in the wild was 16, the otter’s life expectancy is rather low, between four and five years. Within a community individuals communicate and mark their home ranges by using scented marks, “spraints”, usually left in prominent sites, near water, on tree stumps, under bridges or at the junction of roads. It is these scent marks that, when freshly collected, provide DNA from which a genetic fingerprint can be obtained.


## 2.2. The Otter in The Netherlands

The otter disappeared from The Netherlands in 1989. In 2002 a reintroduction programme was started with wild and captive otters, in order to start a new, self-sustainable, resident population (cf. [20]). From June 2002 until November 2008, a total of 31 animals, originating from Eastern and Northern Europe, have been released in a lowland peat marsh area in the north of the Netherlands (Appendix A), ca. 200 km2 in extent and geographically located at 52°38′–52°50′N, 5°53′–6°09′E. The area consists of a mosaic of peat grasslands, reed beds (Phragmites australis), and swamp woodland, mainly composed of Willow (Salix spp.), Black alder (Alnus glutinosa), and Birch species (Betula spp.) intersected by small ditches and pools. Most of the surrounding landscape is intensively farmed. The area is divided in three main sections: the Weerribben and the Rottige Meenthe solely consist of peat grassland, woodland and ditches, while in the Wieden also larger lakes are present (Figure 1).Figure 1
Geographical location of the study area. Within black borders the area where 31 otters (see AppendixA) were released between 2002 and 2008. The expanding population was monitored using telemetry and noninvasive genetic sampling.The released animals varied between 1 and 5 years old. Before release, tissue and blood samples were taken for DNA fingerprinting, all animals were tagged with a transponder, and a radio transmitter was implanted intraperitoneally (by a veterinarian at the Burgers Zoo (Arnhem, The Netherlands)). The study was conducted in accordance with Dutch legislation on the protection and welfare of vertebrate animals used for experimental and other scientific purposes.


## 2.3. Monitoring

Data collection in mark-recapture protocols requires individually marked (recognisable) animals. Here, identification was done according to two methods, without (re)capturing the individuals. Released otters were observed by either telemetry or noninvasive genetic sampling. Offspring (native borns) could only be observed by noninvasive genetic sampling [18]. It should be noted that the full population at the sample site consists only of released individuals and their offspring.Released otters were equipped with a radio transmitter, and their presence and movements could be recorded during the initial stages of the project. Each animal was searched at least once a week for almost one year, until transmitters failed. Based on the analysis of the blood samples taken before release, DNA fingerprinting could be performed by microsatellite analysis (AppendixB). Both released and native-born individuals from the release area were subsequently traced by applying noninvasive genetic methodology and using DNA as a tag [8, 9, 21].Fresh droppings consisting of scented marks (spraints) of the animals were collected at regular intervals during the winter periods (see for full details [18]). DNA was extracted from these spraints to build DNA profiles of the population (Appendix B). A profile could either reveal a new pattern or match a released otter’s pattern. In the latter case, it was considered a “resighting” in this individual’s encounter history. A new profile had to originate from a native-born otter, because no otters were present before the first one was released. For these native-born animals, the parents could be identified. For most animals several spraints were found during one year (range 1–27; [18]).In addition, dead animals were reported together with the cause of death (mainly traffic collision but also drowning, trapping in cages designed for muskrats and natural cause). These dead individuals were identified by their DNA fingerprint. Together, the three types of observations (telemetry, DNA profiles, dead individuals) allowed us to compile detailed encounter histories for the individuals. The multiple observations from both telemetry and DNA analysis are summarised in the event histories of the animals: every period of 6 months (from July 2002 to January 2009) is an observation interval in which they could either or not be seen alive (1/0) and either or not be recovered dead (1/0) (see Figure2). Per half-year interval, two places in the encounter history were filled as required for analysis with a joint model for live and dead data with the program MARK [22]. This implies that for the released animals both telemetry and genetic data were used, for the first two considered half-year periods per animal, and genetic data for the remainder of their lifespan; while native-born individuals only had encounter histories based on genetic data. In this way the rather small dataset of 81 animals (31 released and 50 native-born animals) and 13 encounter observations was constructed (each encounter history was a string of 26 zeros and ones, where at odd places information about resighting alive was summarized and at even place information on whether or not the individual was recovered dead).Figure 2
Fate diagram for the live-dead model. The convention used for the parameterization is the one developed by Seber [23], with S representing the probability of survival from one occasion to the next, p the probability of recapture in the following interval given the animal survived the interval, r the probability of being recovered dead, and F the site fidelity. The succession of events occurring during the period of observation is coded in an event history with information on alive or dead in the odd and even numbered places of the event history, respectively. For example, an individual marked and released on the first occasion, which survives the interval between occasion 1 and occasion 2 because it was recaptured at occasion 2, is subsequently recovered dead during the second interval. Its event history for these two intervals is 1011 (seen alive, not reported dead, seen alive, reported dead).


## 2.4. Statistical Analysis with RMark

Statistical analysis has been conducted with RMark [24], an R package using the program MARK [22] under the R interface.Data included both live resightings and recoveries of dead animals. These were combined into encounter histories as input for the Burnham model for joint data [25], allowing for increased precision in the estimation of the parameters compared to models treating live or dead data separately.For the parameters we followed the convention as developed by Seber [23], with S representing the probability of survival from one occasion to the next, p the probability of recapture in the following interval given the animal survived the interval, r the probability of being recovered dead, and F the site fidelity. Our focus is on survival as it is needed for the construction of a matrix model. The other estimations are not explicitly used, but allow accounting for all possible fates of the animals under the assumptions of the model and avoiding confusion between them. For example, if one does not know the probability that an animal stays in the sample area (F), the events of dying and leaving the area are both treated as mortality.For fitting the models to the data, we assumed for the parameters (a) the site fidelity is 1 since we are able to identify the animals anywhere in the area, even when they leave the release region, because there is only one (reintroduced) otter population in The Netherlands, (b) all parameters are time independent, because of the low number of animals and years involved, and (c)S and p can either be dependent on both sex and stage (defined by the age when released and time elapsed since release or birth) or sex or stage or be constant among sexes and stages.We built (st)age structure into the models, because in advance it is not obvious whether survival of differently aged young otters is different from adult survival. For that reason we started by considering three stages of the animals: juveniles (from birth to one year old), subadults (between one and two years of age), and adults (older than two years). This was possible because the age of the animal when released was known, and for the native borns we could estimate the most likely birth date based on the first encounter and the inferred ancestry of the parents. Native borns are at least ca. 6 months old before they start producing spraints. We also considered only two age classes (making a distinction for individuals less than one year old and older than one year) to allow for one (common) survival probability for subadults and adults. Fitting a model for every parameter combination resulted in fitting 32 different models (four possibilities of dependency forS(sex, stage, sex * stage, or constant), four for p(same as for S), and two for r(sex dependent or constant)).The first step in the modelling process is to find the most general and biologically relevant model (hereS(~sex * stage) p(~sex * stage) r(~sex) see Table 1) and to test for its quality of adjustment, stating that, if this adequately fits the data, then the nested models also do. Estimations for the Burnham model are computed in MARK by likelihood maximization assuming that data can be appropriately represented by a multinomial distribution. If variability is greater than what can be explained by model structure and the multinomial assumption, they are called overdispersed. This occurs when there is interdependence in the data or when the assumption of parameter homogeneity among individuals of the same group is violated. As a consequence variances of the parameters are underestimated, and model selection favours complicated models. These effects can be accounted for with a reliable estimate of overdispersion. We used the variance inflation factor c-hat(ĉ) [26]. If ĉ  is greater than 3, the data show a relatively important overdispersion [27]; lower values are not problematic.Table 1
Model selection (two-stage models) for survival (S), resighting probability (p), and recovery probability (r) with the site fidelity (F) of the otters set to 1. For each model the corrected quasi likelihood Akaike’s information criterion (QAICc), the difference in QAICc between the current model and the best fitting model (ΔQAICc), the number of estimable parameters (np), their weight (QIACc weight), and deviance (Qdev) are given. Fitted models are sorted by increasing QAICc, the first being the one that best represents the data. We present here only those models for which QAICc weight was greater than 0.01. These results are obtained after adjusting for overdispersion with ĉ=1.338.
FormulanpQAICcΔQAICcQAICc weightQDevS(~sex * stage) p(~stage) r(~1)7524.36700.478386.088S(~sex * stage) p(~stage) r(~sex)8525.6071.2400.257385.191S(~sex * stage) p(~sex * stage) r(~1)9527.0742.7060.123384.504**S(~sex * stage) p(~sex * stage) r(~sex)10528.2363.8690.069383.495S(~sex) p(~stage) r(~1)5531.7097.3420.012397.648S(~stage) p(~stage) r(~1)5532.0567.6890.010397.995**denotes the most general model and bold characters the two models used for estimating the parameter values (see Table2).Testing for goodness of fit is usually done with a parametric bootstrapping procedure [26]. The values of the parameters found for the model are used to generate a new dataset (simulated encounter history) according to the assumptions of the model, so that the simulated data present no overdispersion. For every simulation the deviance and other adjustment criteria are computed, sorted by increasing order, and compared to those of the tested model. The percentage of simulated deviances larger than the observed deviance gives us the probability that the tested model presents lack of fit. In addition, this method gives us a measure of the magnitude of the overdispersion of data, the “variance information criterion” or adjusted ĉ. It is obtained by dividing the observed ĉ and the mean ĉ from bootstraps. This adjusted value is subsequently used to have a better model fit [28]Usually, model selection of the fitted models is based on Akaike’s Information Criterion corrected for small sample size (AICc): the model with minimum AICc is the best [29]. When, however, there is evidence of overdispersion, the model selection is based on the Quasi likelihood Akaike Information Criterion (QAICc), proposed by Barker and White [28]. Let L be the likelihood of the model,ĉ the inflation criterion, n the number of observations, and k the number of parameters; we have(1)QAICc=-2log⁡⁡(L)ĉ+2k+2k(k+1)n-k-1.We also used QAICc weights which are the relative percentage of “support” from the data to a given model in the whole set of considered models. These weights sum to one. Models withwi<0.01 have a bad support from the data [26].


## 2.5. Matrix Model Development

The next step is to determine how the population evolves in the long term given the demographic parameters, reproduction and survival. We build a stage-structured Leslie matrix model for females only [17], considering that they drive the population dynamics as there are enough males to fertilize females [18]. Based on the results of the survival analysis, we develop a matrix model for either three or two stages (see Appendix C, (C.1)). We assume a 1 : 1 sex ratio at birth, density independence of the population growth rate as the carrying capacity of the environment is most probably not yet reached [18], and a birth-flow population (breeding season is continuous in the population, and female otters can have cubs all year round; [19]). We used a one-year time step, because we easily could obtain yearly estimates of survival and fertility parameters from our dataset.The fertility estimate, used in the matrix model, was based on the field observations during each year [18]. From the combined genetic and demographic database, it was possible to estimate (i) the number of identified offspring per successful mating, (ii) the total number of breeding pairs, and (iii) the number of adult females that did not reproduce. Fertility was calculated as the total number of new offspring divided by the total number of adult females present, both reproducing and nonreproducing. For the dominant eigenvalue we performed a bootstrapping procedure to assess the range in which the mean value of this yearly population growth factor falls. The 1,000 simulations are based upon intervals that span from the (mean – 2 * SE of the mean) to the (mean + 2 * SE of the mean) for the fertility (Ra=2fa), the juvenile (sj), and adult survival (sa).


## 2.6. Elasticity Analysis

Elasticity analysis is used to identify which matrix element influences population development most. Therewith, it is indicated which element should require closer attention in the perspective of improving conservation of the species. Elasticity of the yearly growth rateλ (the dominant eigenvalue of the matrix) to each of the demographic rates is the proportional change in λ when the parameter undergoes a very small variation, while other elements being held equal. Because matrix elements can be formulas of different survival and/or reproduction parameters, we also have to analyse the elasticity of these lower-level parameters (see Appendix C).


## 3. Results

### 3.1. Statistical Analysis with RMark

The data were first analysed with a model allowing for stage and sex dependence onS and p and sex dependence on r (S(~sex * stage) p(~sex * stage) r(~sex)). A bootstrap procedure with 5,000 simulations showed that the model did not describe the data well (P<0.0002), which might be due to capture heterogeneity. However, ĉ was only slightly larger than 1 (ĉ=1.338) suggesting that variance inflation was not too important. Therefore, we followed the procedure as suggested by Barker and White [28] and adjusted the likelihood parameters over this value, giving better model support. Because of this adjustment, the model selection procedure was based on the quasi likelihood Akaike’s Information Criterion (QAICc; see (1) above).From the general model with three stages (juvenile, subadult, and adult), we fitted several reduced models (see Section2.4), for instance, a two-stage model (juvenile and adult). A likelihood ratio test comparing the two general models with two and three stages confirmed that the most parsimonious model (with two stages) was fitting best (χ2=3.05, df=4, P=0.55). Therefore, we report in Table 1 only the results of the two-stage models with a distinction between juveniles (less than one year old) and adults (older than one year). Moreover, only models with relative high QAICc weights are presented in Table 1.Of the two-stage models, the two models which were best supported by the data (minimum QAICc) both had a stage- and sex-dependent survival parameterS, a stage-dependent recapture probability p which was the same for males and females, and either a constant or a sex-dependent recovery parameter r. These two models together have a weight of 73% (for all considered models, column weight sums to one; Table 1). As these two models differ in QAICc by less than 2, they were considered equivalent [30], and parameter estimates (Table 2) were obtained by model averaging (cf. [26]). Females have a greater survival rate than males, given any age group (Table 2). On the contrary, older animals were recaptured more than younger ones irrespective of sex (recapture probability p). Males (r=0.620) are more easily recovered dead than females (r=0.564).Table 2
parameter value estimates after averaging of the two best fitting models from Table1.
ParameterEstimateStandard errorSurvival rateS males juveniles0.8440.079S males adults0.4600.082S females juveniles (sj)0.9460.061S females adults (sa)0.7860.054Recapture probabilityp juveniles0.3300.074p adults0.6320.039Recovery probabilityr males0.6200.090r females0.5640.110


### 3.2. Population Growth Rate and Elasticity Analysis

Since native-born otters stay within the den for the first half year and do not produce spraints during this time, we can only estimate the net reproduction per female per year of cubs that survived the first half year (Ra). This mean number is calculated by dividing the total number of otters born in each year by the cumulated sum of adult females that were present in the year before. Our estimated value was 1.38 (±SD 2.08, for 40 otters) offspring per adult female per year with a minimum of no offspring per female and a maximum of four. In the model (2), the parameter fa is the number of female offspring per adult female per year and, because only half of the offspring is female [18], half the value of Ra. The model represented in Figure 3 and developed in Appendix C resulted in the matrix in (2). The estimated annual population growth rate λ of the matrix model was 1.26. This value of λ implies that the number of animals doubles after three years. The bootstrap simulation with the intervals (as stated in Section 2.5) resulted in a range for λ of [1.06, 1.42] implying that the population is sustainable in a nonchanging environment:(2)(0fasjsasa)=(00.6900.9460.7860.786).Figure 3
Two-stage life cycle graph for females. Animals spend one year in the juvenile stage before reaching maturity at the end of their second year.The formulas resulting from the elasticity analysis are given in (C.4) in the appendix. The numerical calculation indicates that the parameter affecting λ the most is adult female survival (e(sa)=0.59, a change of 1% in adult female survival rate will increase λ by 0.59%). The reproductive rate has an influence on λ of e(fa)=0.273. Juvenile survival has the smallest effect on  λ(e(sj)=0.137), probably because its values are already high (Table 2). Figure 4 illustrates how annual growth rate varies with adult female survival sa and adult fertility fa when juvenile survival is fixed. We give contour curves for sj equal 0.12 and 0.946 (covering the range of possible values for a growing population). The line where λ equals one separates the area where the population declines and the area where it shows exponential growth. With the current juvenile survival rate of 0.946 (Figure 4(b)), adult female survival could decrease as low as 0.52 with fertility being held equal, and yet the population could maintain itself. Reversely, with constant adult female survival, fertility can be 0.25 and λ still equals one.Sensitivity ofλ to varying adult demographic rates (range [0,1]), when juvenile survival (sj) is set to (a) 0.12 and (b) 0.964 (observed juvenile mortality, see Table 2). Bold line (λ=1) represents the separation between set of parameters resulting in decline of the population and those resulting in growth. Dotted lines indicate the current values for sa (adult female survival) and fa (net reproduction of female native-born cubs per adult female per year).
(a)
Sjuv=0.12(b)
Sjuv=0.946


## 3.1. Statistical Analysis with RMark

The data were first analysed with a model allowing for stage and sex dependence onS and p and sex dependence on r (S(~sex * stage) p(~sex * stage) r(~sex)). A bootstrap procedure with 5,000 simulations showed that the model did not describe the data well (P<0.0002), which might be due to capture heterogeneity. However, ĉ was only slightly larger than 1 (ĉ=1.338) suggesting that variance inflation was not too important. Therefore, we followed the procedure as suggested by Barker and White [28] and adjusted the likelihood parameters over this value, giving better model support. Because of this adjustment, the model selection procedure was based on the quasi likelihood Akaike’s Information Criterion (QAICc; see (1) above).From the general model with three stages (juvenile, subadult, and adult), we fitted several reduced models (see Section2.4), for instance, a two-stage model (juvenile and adult). A likelihood ratio test comparing the two general models with two and three stages confirmed that the most parsimonious model (with two stages) was fitting best (χ2=3.05, df=4, P=0.55). Therefore, we report in Table 1 only the results of the two-stage models with a distinction between juveniles (less than one year old) and adults (older than one year). Moreover, only models with relative high QAICc weights are presented in Table 1.Of the two-stage models, the two models which were best supported by the data (minimum QAICc) both had a stage- and sex-dependent survival parameterS, a stage-dependent recapture probability p which was the same for males and females, and either a constant or a sex-dependent recovery parameter r. These two models together have a weight of 73% (for all considered models, column weight sums to one; Table 1). As these two models differ in QAICc by less than 2, they were considered equivalent [30], and parameter estimates (Table 2) were obtained by model averaging (cf. [26]). Females have a greater survival rate than males, given any age group (Table 2). On the contrary, older animals were recaptured more than younger ones irrespective of sex (recapture probability p). Males (r=0.620) are more easily recovered dead than females (r=0.564).Table 2
parameter value estimates after averaging of the two best fitting models from Table1.
ParameterEstimateStandard errorSurvival rateS males juveniles0.8440.079S males adults0.4600.082S females juveniles (sj)0.9460.061S females adults (sa)0.7860.054Recapture probabilityp juveniles0.3300.074p adults0.6320.039Recovery probabilityr males0.6200.090r females0.5640.110


## 3.2. Population Growth Rate and Elasticity Analysis

Since native-born otters stay within the den for the first half year and do not produce spraints during this time, we can only estimate the net reproduction per female per year of cubs that survived the first half year (Ra). This mean number is calculated by dividing the total number of otters born in each year by the cumulated sum of adult females that were present in the year before. Our estimated value was 1.38 (±SD 2.08, for 40 otters) offspring per adult female per year with a minimum of no offspring per female and a maximum of four. In the model (2), the parameter fa is the number of female offspring per adult female per year and, because only half of the offspring is female [18], half the value of Ra. The model represented in Figure 3 and developed in Appendix C resulted in the matrix in (2). The estimated annual population growth rate λ of the matrix model was 1.26. This value of λ implies that the number of animals doubles after three years. The bootstrap simulation with the intervals (as stated in Section 2.5) resulted in a range for λ of [1.06, 1.42] implying that the population is sustainable in a nonchanging environment:(2)(0fasjsasa)=(00.6900.9460.7860.786).Figure 3
Two-stage life cycle graph for females. Animals spend one year in the juvenile stage before reaching maturity at the end of their second year.The formulas resulting from the elasticity analysis are given in (C.4) in the appendix. The numerical calculation indicates that the parameter affecting λ the most is adult female survival (e(sa)=0.59, a change of 1% in adult female survival rate will increase λ by 0.59%). The reproductive rate has an influence on λ of e(fa)=0.273. Juvenile survival has the smallest effect on  λ(e(sj)=0.137), probably because its values are already high (Table 2). Figure 4 illustrates how annual growth rate varies with adult female survival sa and adult fertility fa when juvenile survival is fixed. We give contour curves for sj equal 0.12 and 0.946 (covering the range of possible values for a growing population). The line where λ equals one separates the area where the population declines and the area where it shows exponential growth. With the current juvenile survival rate of 0.946 (Figure 4(b)), adult female survival could decrease as low as 0.52 with fertility being held equal, and yet the population could maintain itself. Reversely, with constant adult female survival, fertility can be 0.25 and λ still equals one.Sensitivity ofλ to varying adult demographic rates (range [0,1]), when juvenile survival (sj) is set to (a) 0.12 and (b) 0.964 (observed juvenile mortality, see Table 2). Bold line (λ=1) represents the separation between set of parameters resulting in decline of the population and those resulting in growth. Dotted lines indicate the current values for sa (adult female survival) and fa (net reproduction of female native-born cubs per adult female per year).
(a)
Sjuv=0.12(b)
Sjuv=0.946


## 4. Discussion

### 4.1. Noninvasive Genetic Monitoring and Vital Rates

Noninvasive genetic methods have proven to be a valuable tool in wildlife research [9, 12, 31]. One area of wildlife research, in particular, that has benefited from molecular genetics is the estimation of animal abundance. Historically, estimation of demographic parameters on species that are rare, elusive, difficult, or expensive to capture has been limited due to small sample sizes. Advances in molecular genetics allow individual identification from the collection of noninvasive samples (mainly hairs and faeces) and the use of DNA as a “molecular tag.” This has been applied to study the population biology of, amongst others, bears [15], cougars [32], coyotes [33], otters [18], wolverines [34], and wolves [35].In this study we take the genetic approach one step further by applying noninvasive genetic monitoring, that is, adding a time component, in estimating the viability of reintroduced populations (cf. [9]). By using tools from demography and genetics, we obtained estimates of both survival and reproduction for a reintroduced otter population. Using these estimates in a matrix model framework allowed us to assess the population viability since the start of the reintroduction and to determine the boundary conditions for the vital rates under which the population would not be any more sustainable. This provided us with the necessary information for the field managers on how to manage the population; for example, how large mortality of females or males can we still afford, or should infrastructural changes prioritise on juveniles or adults?The statistical model used for analysing the field observations revealed overdispersion of the data under the assumptions of the Burnham model. This model, like all mark-recapture models, assumes that every marked animal in the population at timet has the same probability of recapture. Because the number of observed spraints per individual varies between 1 and 27 [18], this assumption is most likely violated. We, therefore, decided to use only one resighting per individual per half year interval. Even so overdispersion was present, but to the extent that according to Barker and White [28] it was possible to account for, thereby preventing a bias in favour of overadjusted models.We modelled the population dynamics using a simple Leslie matrix. We observed that the annual growth rate of the population (λ) of 1.26 is more sensitive to adult survival than to fertility. The latter findings are in accordance with the results of Björklund and Arrendal [36] on a Swedish otter community. Thus, litter size and the ability of cubs to survive are less important for population growth than adult female survival once a female is able to mate several times in her life.The yearly growth rate, as determined with the matrix model, is larger than 1. If we take the variation in the estimates of survival and fertility into account with a bootstrap model, the range for the yearly growth rate becomes [1.06,1.42], implying a growing population. This is confirmed by the realised population that shows growth like that in our “mean field” modelled otter population [18].Our approach should be equally applicable to other elusive carnivore species. Recently, several other researchers have more or less applied the same approach [14, 15] to estimate the growth rate of small bear populations. De Barba et al. [15] used their CMR estimates as input for the widely used VORTEX model [37]. We deliberately decided to build our own matrix model, since this allows us to incorporate mating structure and inbreeding depression in more advanced models (Seignobosc, in prep.). Marucco et al. [16] summarised the necessities for a full noninvasive genetic CMR analysis and emphasized the need for (i) the design of an appropriate sampling scheme, (ii) established genetic protocols, and (iii) checking the assumptions of the CMR models. If either one of these requirements fails, the results of the CMR analysis will produce less reliable estimates. The first requirement is hard to achieve in reintroduction studies, since in most cases no replications are performed and sampling is constrained by the field conditions. Concerning the third requirement, CMR models assume equal detection probability for all individuals and are sensitive to individual detection heterogeneity (IDH) which easily leads to an underestimation of the population abundance [16, 38].


### 4.2. The Dutch Otter Reintroduction Project

Despite the limited number of observations due to the small population size, our genetic CMR-modelling approach gave valuable insights on how to manage the otter reintroduction programme. The results revealed female-biased survival. Subadult males had a higher incidence of being recovered dead, mostly as victims of traffic incidents, than subadult females [18]. The mating structure of the population is such that a few dominant males patronise the release area [18] and subdominant and subadult males are pushed to the borders of the area or leave. Subadult females, on the other hand, stay in the area close to their mothers’ territory. These observations caused a change in the original release plan: a new release of otters should in first instance concentrate on releasing females in neighbouring areas since the males will come in time.We used a simple model (Figure3) to calculate the growth rate resulting from the projection of the population matrix. However, this does not take into account the counteracting effects that can be at stake in the population but are not visible yet. In particular, inbreeding, due to increased kinship among the population members [18], probably becomes problematic in the next years. Thus, it should be included in the modelling process. Furthermore, detailed biological insight on this population could be gained by improving the matrix model with social structure. For example, the mating behaviour within the population suggests polygyny, with few males mating with several females whereas some never mate [18]. Integration of inbreeding and polygyny in the model permits us to study their effect in the long term on population growth and repartition among the different sex and age groups.Despite its simplicity our model yielded applicable results that can help orientating the conservation decisions that should be taken in the reintroduction programme. Some results were already known qualitatively (protection from road kills is necessary), but we provide the quantitative background support for implementing such measures and give an upper limit for yearly mortality under the assumption that we want to maintain the population. Adult females play an important role in the population dynamics, and, even if their mortality is low, they are nonetheless heavily exposed to traffic incidents. Thus, the need to prevent them from approaching dangerous roads is real.


## 4.1. Noninvasive Genetic Monitoring and Vital Rates

Noninvasive genetic methods have proven to be a valuable tool in wildlife research [9, 12, 31]. One area of wildlife research, in particular, that has benefited from molecular genetics is the estimation of animal abundance. Historically, estimation of demographic parameters on species that are rare, elusive, difficult, or expensive to capture has been limited due to small sample sizes. Advances in molecular genetics allow individual identification from the collection of noninvasive samples (mainly hairs and faeces) and the use of DNA as a “molecular tag.” This has been applied to study the population biology of, amongst others, bears [15], cougars [32], coyotes [33], otters [18], wolverines [34], and wolves [35].In this study we take the genetic approach one step further by applying noninvasive genetic monitoring, that is, adding a time component, in estimating the viability of reintroduced populations (cf. [9]). By using tools from demography and genetics, we obtained estimates of both survival and reproduction for a reintroduced otter population. Using these estimates in a matrix model framework allowed us to assess the population viability since the start of the reintroduction and to determine the boundary conditions for the vital rates under which the population would not be any more sustainable. This provided us with the necessary information for the field managers on how to manage the population; for example, how large mortality of females or males can we still afford, or should infrastructural changes prioritise on juveniles or adults?The statistical model used for analysing the field observations revealed overdispersion of the data under the assumptions of the Burnham model. This model, like all mark-recapture models, assumes that every marked animal in the population at timet has the same probability of recapture. Because the number of observed spraints per individual varies between 1 and 27 [18], this assumption is most likely violated. We, therefore, decided to use only one resighting per individual per half year interval. Even so overdispersion was present, but to the extent that according to Barker and White [28] it was possible to account for, thereby preventing a bias in favour of overadjusted models.We modelled the population dynamics using a simple Leslie matrix. We observed that the annual growth rate of the population (λ) of 1.26 is more sensitive to adult survival than to fertility. The latter findings are in accordance with the results of Björklund and Arrendal [36] on a Swedish otter community. Thus, litter size and the ability of cubs to survive are less important for population growth than adult female survival once a female is able to mate several times in her life.The yearly growth rate, as determined with the matrix model, is larger than 1. If we take the variation in the estimates of survival and fertility into account with a bootstrap model, the range for the yearly growth rate becomes [1.06,1.42], implying a growing population. This is confirmed by the realised population that shows growth like that in our “mean field” modelled otter population [18].Our approach should be equally applicable to other elusive carnivore species. Recently, several other researchers have more or less applied the same approach [14, 15] to estimate the growth rate of small bear populations. De Barba et al. [15] used their CMR estimates as input for the widely used VORTEX model [37]. We deliberately decided to build our own matrix model, since this allows us to incorporate mating structure and inbreeding depression in more advanced models (Seignobosc, in prep.). Marucco et al. [16] summarised the necessities for a full noninvasive genetic CMR analysis and emphasized the need for (i) the design of an appropriate sampling scheme, (ii) established genetic protocols, and (iii) checking the assumptions of the CMR models. If either one of these requirements fails, the results of the CMR analysis will produce less reliable estimates. The first requirement is hard to achieve in reintroduction studies, since in most cases no replications are performed and sampling is constrained by the field conditions. Concerning the third requirement, CMR models assume equal detection probability for all individuals and are sensitive to individual detection heterogeneity (IDH) which easily leads to an underestimation of the population abundance [16, 38].


## 4.2. The Dutch Otter Reintroduction Project

Despite the limited number of observations due to the small population size, our genetic CMR-modelling approach gave valuable insights on how to manage the otter reintroduction programme. The results revealed female-biased survival. Subadult males had a higher incidence of being recovered dead, mostly as victims of traffic incidents, than subadult females [18]. The mating structure of the population is such that a few dominant males patronise the release area [18] and subdominant and subadult males are pushed to the borders of the area or leave. Subadult females, on the other hand, stay in the area close to their mothers’ territory. These observations caused a change in the original release plan: a new release of otters should in first instance concentrate on releasing females in neighbouring areas since the males will come in time.We used a simple model (Figure3) to calculate the growth rate resulting from the projection of the population matrix. However, this does not take into account the counteracting effects that can be at stake in the population but are not visible yet. In particular, inbreeding, due to increased kinship among the population members [18], probably becomes problematic in the next years. Thus, it should be included in the modelling process. Furthermore, detailed biological insight on this population could be gained by improving the matrix model with social structure. For example, the mating behaviour within the population suggests polygyny, with few males mating with several females whereas some never mate [18]. Integration of inbreeding and polygyny in the model permits us to study their effect in the long term on population growth and repartition among the different sex and age groups.Despite its simplicity our model yielded applicable results that can help orientating the conservation decisions that should be taken in the reintroduction programme. Some results were already known qualitatively (protection from road kills is necessary), but we provide the quantitative background support for implementing such measures and give an upper limit for yearly mortality under the assumption that we want to maintain the population. Adult females play an important role in the population dynamics, and, even if their mortality is low, they are nonetheless heavily exposed to traffic incidents. Thus, the need to prevent them from approaching dangerous roads is real.



---

*Source: 870853-2011-11-29.xml*