Anaïs Charbonnel
email: anais.charbonnel@espaces-naturels.fr

Laëtitia Buisson
email: laetitia.buisson@univ-tlse3.fr

Marjorie Biffi
email: m.biffi@live.fr

Frank D ' Amico
email: frank.damico@univ-pau.fr

Aurélien Besnard
email: aurelien.besnard@cefe.cnrs.fr

Stéphane Aulagnier
email: stephane.aulagnier@toulouse.inra.fr

Frédéric Blanc
email: frederic.blanc@espaces-naturels.fr

François Gillet
email: francois.gillet@espaces-naturels.fr

Vincent Lacaze
email: vincent.l@ariegenature.fr

Johan R Michaux
email: michaux@supagro.inra.fr

Mélanie Némoz
email: melanie.nemoz@espaces-naturels.fr

Christian Pagé
email: christian.page@cerfacs.fr

José Miguel Sanchez-Perez
email: jose-miguel.sanchez-perez@univ-tlse3.fr

Sabine Sauvage
email: sabine.sauvage@univ-tlse3.fr

Pascal Laffaille
email: pascal.laffaille@ensat.fr

Integrating hydrological features and genetically validated occurrence data in occupancy modelling of an endemic and endangered semi-aquatic mammal, Galemys pyrenaicus, in a Pyrenean catchment

Keywords:

As freshwater habitats are among the most endangered, there is an urgent need to identify critical areas for conservation, especially those that are home to endangered species. The Pyrenean desman (Galemys pyrenaicus) is a semi-aquatic mammal whose basic ecological requirements are largely unknown, hindering adequate conservation planning even though it is considered as a threatened species. Species distribution modelling is challenging for freshwater species. Indeed, the complexity of aquatic ecosystems (e.g., linear and hierarchical ordering) must be taken into account as well as imperfect sampling. High-quality and relevant hydrological descriptors should also be used. To understand the influence of environmental covariates on the occupancy and detection of the Pyrenean desman, we combine both a robust sign-survey data set (i.e. with genetic validation ensuring true presence information) and a hydrological model to simulate the flow regime across a whole catchment. Markovian site-occupancy analysis, taking into account sign detection and based on spatially adjacent replicates, indicated a positive influence of heterogeneity of substrate and shelters, and a negative influence of flow variability on Pyrenean desman detection. This valuable information should help to improve monitoring programs for this endangered species. Our results also highlighted a spatially clustered distribution and a positive influence of stream flow and number of tributaries on occupancy. Hence, modifications of flow regime (e.g. hydropower production, irrigation, climate change) and habitat fragmentation appear to be major threats for this species, altering the connectivity between tributaries and the mainstream river as well as between adjacent sub-catchments.

Introduction

Freshwater habitats hold a notable biodiversity with for example, one third of vertebrate species being restricted to this ecosystem [START_REF] Dudgeon | Freshwater biodiversity: importance, threats, status and conservation challenges[END_REF]. However, freshwater habitats are one of the most endangered ecosystems in the world [START_REF] Vörösmarty | Global threats to human water security and river biodiversity[END_REF] and human-induced alterations of the natural river conditions strongly affect aquatic biodiversity [START_REF] Dudgeon | Freshwater biodiversity: importance, threats, status and conservation challenges[END_REF][START_REF] Vörösmarty | Global threats to human water security and river biodiversity[END_REF]. Extinction rates of freshwater fauna are currently extremely high [START_REF] Allan | Climate change and freshwater ecosystems[END_REF] with around 15,000 freshwater species worldwide already extinct or imperiled as a result of human activity [START_REF] Strayer | Freshwater biodiversity conservation: recent progress and future challenges[END_REF]. As a consequence, there is an urgent need to identify critical areas for conservation of freshwater biodiversity especially for rare, endemic and endangered species.

Among rare freshwater species, the Pyrenean desman (Galemys pyrenaicus) is one of the less well-known European mammals. The distribution of this small semi-aquatic species is restricted to the Pyrenees (Andorra, France and Spain), as well as parts of northern and central Spain and northern Portugal. In the French Pyrenees, the species lives in mountain brooks, cold and well oxygenated water courses from sea level to 2700 m [START_REF] Némoz | A French Conservation Action Plan for the Pyrenean Desman Galemys pyrenaicus[END_REF]. The Pyrenean desman is becoming increasingly threatened, triggering several conservation regulations [START_REF] Fernandes | Galemys pyrenaicus[END_REF]). Yet, even basic knowledge such as distribution range and habitat preferences that are essential for conservation planning are not complete for this species [START_REF] Aymerich | Els micromamıfers semiaquatics d'Andorra: distribucio i estat de conservacio[END_REF][START_REF] Barbosa | Transferability of environmental favourability models in geographic space: The case of the Iberian desman (Galemys pyrenaicus) in Portugal and Spain[END_REF][START_REF] Barbosa | Use of coarse-resolution models of species' distributions to guide local conservation inferences[END_REF][START_REF] Nores | Aproximación a la metodología y estudio del área de distribución, estatus de población y selección de hábitat del desmán (Galemys pyrenaicus) en la Península Ibérica[END_REF][START_REF] Nores | Informe sobre la situación del Desmán Ibérico (Galemys pyrenaicus) en España. Seminariode conservación de argaritifera margaritifera y Galemys pyrenaicus en la Península Ibérica[END_REF][START_REF] Palmeirim | Galemys pyrenaicus[END_REF][START_REF] Queiroz | Status and conservation of Desmaninae in Europe[END_REF].

The environmental factors influencing the spatial distribution of species can be identified by the use of Species Distribution Models (SDMs). They model the statistical relationships between species presence records and environmental variables, and may be used to predict habitat suitability for species in unsampled areas [START_REF] Guisan | Predictive habitat distribution models in ecology[END_REF]. A few studies have used SDMs for understanding the ecological requirements of the Pyrenean desman and have reported a positive influence of topographic variability, and a strong but contrasting influence of climatic variables, depending on the study area [START_REF] Barbosa | Transferability of environmental favourability models in geographic space: The case of the Iberian desman (Galemys pyrenaicus) in Portugal and Spain[END_REF][START_REF] Morueta-Holme | Climate change risks and conservation implications for a threatened small-range mammal species[END_REF][START_REF] Williams-Tripp | Modeling rare species distribution at the edge: the case for the vulnerable endemic Pyrenean desman in France[END_REF]. A negative effect of the density of urban areas was also identified [START_REF] Barbosa | Transferability of environmental favourability models in geographic space: The case of the Iberian desman (Galemys pyrenaicus) in Portugal and Spain[END_REF]. As these studies were applied at large scales (e.g., the whole Iberian Peninsula) and with coarse resolution (e.g., 10 km), they did not take into account the particular features of freshwater environments.

Indeed, SDMs for aquatic species often fail to account fully for links between organism occurrence and environmental constraints imposed by river networks [START_REF] Jähnig | Modelling of riverine ecosystems by integrating models: conceptual approach, a case study and research agenda[END_REF] as they do not take into account the linear configuration of the river network (SDMs built on grid cells with both aquatic and terrestrial ecosystems considered together; e.g. [START_REF] Blank | Using ecological niche modeling to predict the distributions of two endangered amphibian species in aquatic breeding sites[END_REF][START_REF] Domisch | Modelling distribution in European stream macroinvertebrates under future climates[END_REF]. Not accounting for these particular features seems to be appropriate for broad-scale studies as the use of climatic or topographic variables might be more meaningful [START_REF] Bucklin | Comparing species distribution models constructed with different subsets of environmental predictors[END_REF]. As aquatic species' movements are constrained by the spatial orientation of the watercourses and by the connectivity between streams and sub-catchments, such hydrological parameters should however be considered [START_REF] Ottaviani | Modelling semi-aquatic vertebrates' distribution at the drainage basin scale: the case of the otter Lutra lutra in Italy[END_REF], especially for small-scale predictions (e.g. across a river catchment). The importance of hydrological variables on freshwater species ecology and distribution is well known (e.g. [START_REF] Kuemmerlen | Integrating catchment properties in small scale species distribution models of stream macroinvertebrates[END_REF] even for river birds [START_REF] Royan | River birds' response to hydrological extremes: new vulnerability index and conservation implications[END_REF] or semi-aquatic mammals [START_REF] Pedroso | The response of otters to environmental changes imposed by the construction of large dams[END_REF][START_REF] Toner | Muskrat abundance responses to water level regulation within freshwater coastal wetlands[END_REF]. Despite this, hydrological variables (e.g. stream flow) are often ignored in SDMs due to the lack of fine scale spatial data available for studies conducted in large areas. A solution to counterbalance this lack of data may be to simulate flow variables using a hydrological model. One of the most commonly applied, the Soil and Water Assessment Tool (SWAT), is a catchment-scale, physically based model [START_REF] Arnold | Large area hydrologic modeling and assessment Part I: model development[END_REF], running on a daily time step and capable of continuous simulation over a long time period at different spatial scales [START_REF] Gassman | The soil and water assessment tool: historical development, applications, and future research directions[END_REF]. By using spatial information (i.e. topography, climate, soil and land-use), SWAT simulates the hydrological cycle both in space and time (see [START_REF] Neitsch | Soil and Water Assessment Tool Theoretical Documentation, Version 2013[END_REF] for more details). To our knowledge, very few studies have coupled this tool with SDMs to understand the influence of hydrological parameters on the presence of aquatic species, and to predict habitat suitability (but see [START_REF] Jähnig | Modelling of riverine ecosystems by integrating models: conceptual approach, a case study and research agenda[END_REF][START_REF] Kuemmerlen | Integrating catchment properties in small scale species distribution models of stream macroinvertebrates[END_REF].

Another important challenge in SDMs is the quality of species presence-absence data. When surveys are based on the recording of indirect signs, such as faeces, ambiguous signs could lead to the risk of species misidentification, leading to false presences (i.e. wrongly attributed to the species of interest; [START_REF] Miller | Improving occupancy estimation when two types of observational error occur: non-detection and species misidentification[END_REF] or false absences (i.e. wrongly attributed to species other than the species of interest). To overcome this issue, techniques that identify species using faecal DNA analysis are increasingly used [START_REF] Waits | Noninvasive genetic sampling tools for wildlife biologists: a review of applications and recommendations for accurate data collection[END_REF]. Species detection is another major issue. It is well known that the absence of records in the field is a combination of undetected presences (i.e. false absences) and true absences [START_REF] Gu | Absent or undetected? Effects of non-detection of species occurrence on wildlife-habitat models[END_REF]. Semi-aquatic mammals, such as the Pyrenean desman, are particularly sensitive to this detection issue as monitoring is usually based on faeces searches in heterogeneous environments [START_REF] González-Esteban | A comparison of methodologies used in the detection of the Pyrenean desman Galemys pyrenaicus (E. Geoffroy, 1811)[END_REF]. Site occupancy models have been developed to deal with species detection issues at a large scale. They model the probability that a species occupies some sites even though the species has not been detected with any certainty when the sites were visited [START_REF] Mackenzie | Estimating site occupancy rates when detection probabilities are less than one[END_REF]. This class of models requires replication of detection-non detection data at sampling sites. Usually based on temporal replication, recent developments of site occupancy models now allow the use of spatial instead of temporal replicates. Among them, the Markovian occupancy model can be applied when spatially adjacent replicates are available at sites, to test for spatial correlation of occupancy probabilities between replicates [START_REF] Charbonnel | Spatial replicates as an alternative to temporal replicates for occupancy modelling when surveys are based on linear features of the landscape[END_REF][START_REF] Hines | Tigers on trails: occupancy modeling for cluster sampling[END_REF]. To date, this recent model has rarely been applied to investigate the influence of covariates on species distribution (see however [START_REF] Barber-Meyer | Influence of prey depletion and human disturbance on tiger occupancy in Nepal[END_REF][START_REF] Karanth | Monitoring carnivore populations at the landscape scale: occupancy modelling of tigers from sign surveys[END_REF][START_REF] Thorn | Brown hyaenas on roads: Estimating carnivore occupancy and abundance using spatially autocorrelated sign survey replicates[END_REF].

In this study, we applied the Markovian occupancy model to a genetically validated dataset for the Pyrenean desman, and used a hydrological model to simulate flow in the river network of a single catchment in the French Pyrenees. Our aim was to highlight the environmental factors, including hydrological variables, influencing detection and occupancy of this threatened semi-aquatic species for better implementation of conservation plans.

Methods

Study area

The upper river Salat (42-43°N, 0-1°W; drainage area = 1156 km 2 ) is a tributary of the river Garonne, located in the French Pyrenean Mountains (Fig. 1). This catchment has a stream length of 1388 km (CARTHAGE Ó DB, French database on thematic mapping of the Water Agencies, 2011 version http:// www.sandre.eaufrance.fr). Elevation varies between 350 and 2870 m (mean elevation = 1200 m). Mean (±sd) annual rainfall, air temperature and stream flow are 1360.83 mm ± 216.72 (range = 1013.72-1998.49), 9.51 °C ± 1.76 (range = 4.21-12.11;[START_REF] Pagé | Nouvelles projections climatiques à échelle fine sur la France pour le 21ème siècle : les scénarii SCRATCH[END_REF] and 0.78 m 3 /s ± 2.76 (range = 0.00-31.96; SWAT simulations, see 2.4), respectively. Land cover (Corine Land Cover Ó DB, map of the European environmental landscape, version 2006) is dominated by forests (50%) and herbaceous and shrubby vegetation (25%). We focused our study on this catchment because of (i) its representativeness of other French Pyrenean catchments, (ii) its conservation status (this catchment is part of a Natura 2000 site) and (iii) the known presence of Pyrenean desman [START_REF] Bertrand | Répartition géographique et écologie alimentaire du desman des Pyrénées Galemys pyrenaicus (Geoffroy, 1811) dans les Pyrénées françaises[END_REF]. The stream network (CARTHAGE Ó DB) was divided into 1388 1-km-long sections (hereafter simply called sections) for the computation of environmental covariates.

Field sampling

One hundred thirty-one sites (i.e. river transects) were surveyed for this study (Fig. 1). Sites were selected using a spatially balanced, Generalized Random Tessellation Stratified sampling which is known to produce survey designs suitable for aquatic systems [START_REF] Stevens | Spatially balanced sampling of natural resources[END_REF]. Searches for Pyrenean desman faeces were conducted along these river transects which were waded by pairs of skilled observers, meticulously inspecting each emergent rock, tree root or branch in the riverbed (i.e. river banks were not inspected). The number of observers was limited as much as possible across transects to minimize observer bias for sign detection. All faeces detected and suspected of being left by a Pyrenean desman based on colour, size, and position were harvested for genetic analysis. Surveys were conducted between 2011 and 2013 during the summer, when faeces appeared to persist longest, to maximize detection [START_REF] Bertrand | Répartition géographique et écologie alimentaire du desman des Pyrénées Galemys pyrenaicus (Geoffroy, 1811) dans les Pyrénées françaises[END_REF]. We tried not to conduct surveys during or after a period of fluctuating water levels or heavy rainfall to minimize variations in sign detection probabilities (e.g., removal of faeces by rising water levels). Each site was a riverbed transect 500 m-long, which approximately matches the mean home range of the species [START_REF] Melero | New insights into social and space use behaviour of the endangered Pyrenean desman (Galemys pyrenaicus)[END_REF]. Each site was divided into five sub-units (i.e. segments) of equal length (i.e. 100 m) that constituted adjacent spatial replicates, as this segment length appears appropriate for the Pyrenean desman when analyzed with the Markovian occupancy model [START_REF] Charbonnel | Spatial replicates as an alternative to temporal replicates for occupancy modelling when surveys are based on linear features of the landscape[END_REF]. For each segment, information of detection or non-detection of faeces was thus available.

Genetic validation of faeces identification

Genomic DNA from faeces samples was extracted using the Stool Mini Kit (Qiagen Inc., Hilden, Germany). DNA extraction was conducted in a separate room with a UV-sterilised platform where no Pyrenean desman tissue samples had previously been treated. In order to identify the owner species of the faeces sample, we amplified a small cytochrome b fragment of approximately 400 bp using specific primers designed for this study (GPYRF1: 5 0 -TTGTAGAATGGAKCTGAGG-3 0 , GPYRF2: 5 0 -TTCCTTCACGAAA-CAGGATC-3 0 and GPYRR1: 5 0 -GTCGGCTGCTAAAAGTCAGAATA-3 0 ).

PCRs were carried out in a volume of 9 ll containing 0.17 ll of forward primer GPYRF1 and 0.17 ll of reverse primer GPYRR1 (10 lM), 2.89 ll of sterile water, 0.58 ll of dNTPs (10 lM), 1.70 ll of MgCl 2 (25 mM), 3.40 ll of 5X GoTaq Ò buffer reaction (Promega Inc., Madison, USA), 0.09 ll of GoTaq Ò DNA polymerase (Promega Inc., Madison, USA) and 8 ll of DNA. Amplifications were performed in a thermal cycler VWR Unocycler using one activation step at 94 °C for 5 min followed by 40 cycles (denaturation at 94 °C for 50 s, annealing at 52 °C for 45 s, extension at 72 °C for 45 s) and final extension step at 72 °C for 10 min. Three microlitres of the PCR product were amplified in a nested PCR with 14 ll of the PCR mixture described above with additions of 5 ll of sterile water and 0.17 ll of forward primer GPYRF2 (10 lM) in place of GPYRF1.

PCR products were sequenced on an Applied Biosystems Ò 3730 DNA analyzer and verified using CHROMASPRO v 1.5 (http://technelysium.com.au). Sequences were then submitted to the BLAST Ò functionality available on the NCBI website (http://blast.ncbi.nlm. nih.gov). For a detailed description of the protocol, see Gillet et al. (in press).

Simulation of stream flow using SWAT model

The combined use of hydraulic and distribution models involves hydrological information available at the same spatial resolution (here, 1-km river sections). SWAT requires several input datasets using the ArcSWAT interface in ArcGIS 10.0 [START_REF] Winchell | ArcSWAT interface for SWAT user's guide[END_REF]. SWAT uses a topography map to delineate the watershed that was divided into 1165 sub-basins (mean surface area = 100 ha ± 82; mean reach length = 873 m ± 704) with a discretization scale of 50 ha (Fig. 2). In this study, we used (i) a 25 m resolution Digital Elevation Model (ALTI Ó DB -IGN), (ii) a 1 km resolution Digital Soil Map of the World (FAO, 2007), and (iii) a 250 m resolution land cover map (BDOS Ó, Regional Natural Parks of Midi-Pyrénées). The climatic variables used to calibrate the SWAT models included daily rainfall, maximum and minimum air temperature, solar radiation, wind speed and relative humidity. Climatic variables were derived from the SAFRAN Ó DB (1992-2011) which has a spatial resolution of 8 km and accounts for the influence of the topography on the atmospheric variables (Habets et al., 2008). As our study area was located in a mountainous region, we modified snow parameters to calibrate SWAT and improved flow simulations (Appendix A, Supplementary Materials). Observed monthly output flows of three hydropower reservoirs (Fig. 2; Electricity of France) were also included in simulations to increase their accuracy. After running SWAT, simulated stream flows were available in the 1165 sub-basins with a monthly time step, between 1992 and 2011. For each of the 1388 1-km sections, the flow value assigned was that of the sub-basin it was included in. A simulation period of 20 years was chosen to reduce the influence of years with extreme hydrological events (e.g. flood, low water). Average monthly stream flow data recorded from 1992 to 2011 were available at five gauging stations (Fig. 2) and used to calibrate and evaluate the performance of the SWAT simulations using three different metrics: the Spearman correlation coefficient (rho), the coefficient of determination (R 2 ) and the Nash-Sutcliffe Efficiency (NSE) calculated between measured and simulated stream flow (e.g. [START_REF] Kiesel | Incorporating landscape depressions and tile drainages of a northern German lowland catchment into a semi-distributed model[END_REF][START_REF] Moriasi | Model evaluation guidelines for systematic quantification of accuracy in watershed simulations[END_REF]. As is typically done [START_REF] Qi | GIS-based hydrologic modeling in the Sandusky watershed using SWAT[END_REF], the model was calibrated using the gauging station located at the outlet (''Saint-Lizier'') and validated at the four other gauging stations (Fig. 2).

Covariates influencing detection and occupancy

Covariates related to detection probability

Bias in sign detection may arise due to a variety of factors including weather, habitat structure and observer. First, species detection is known to be influenced by the experience of the observer [START_REF] Mackenzie | Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence[END_REF]. Hence, we used the pairs of observers who inspected each transect as a first covariate (OBS; 3 categories). For the Pyrenean desman, we then hypothesized that substrate heterogeneity influences faeces detection as emergent items are supports for faeces deposits. We could expect that the greater the heterogeneity, the higher the detection. Hence, during sampling, observers visually assessed the percentage of heterogeneity of substrate and shelters (e.g. rocks, tree roots or branches) in the riverbed (SUBSTRATE) that was classified in four categories: 0-25%, 26-50%, 51-75% and 76-100%. The 0-25% category indicated a riverbed with homogeneous substrate and few or no shelters whereas the 76-100% category indicated a riverbed with a very high diversity of substrate and shelters (e.g. rocks, tree roots or branches). In addition, rainfall may negatively influence the faeces detection by washing out emergent items, as already shown for the European otter (Lutra lutra; [START_REF] Reid | Detecting detectability: identifying and correcting bias in binary wildlife surveys demonstrates their potential impact on conservation assessments[END_REF]. It could thus be more difficult to find faeces in areas with high annual rainfall as items are regularly washed out. Flow variability might also influence detection with higher variation in stream flow resulting in faeces regularly removed by the fluctuating water level, and thus a lowered detection. Consequently, we used the mean annual rainfall (RAIN, mean of the annual rainfall from 1992 to 2011, mean across the 131 river transects = 1252.80 mm ± 154.75) and the inter-monthly flow variability (FLOW VAR, variance of the twelve monthly flows simulated from SWAT, averaged from 1992 to 2011, mean = 1.81 ± 7.40). All these covariates were calculated for each 1-km section.

Covariates related to occupancy probability

We used four covariates that were not highly correlated (|r| 6 0.60) in the Markovian occupancy model. First, we used the mean monthly flow (FLOW, mean of the twelve monthly flows, averaged from 1992 to 2011, simulation of SWAT models, mean = 1.48 m 3 /s ± 2.76). Second, we calculated eight climatic covariates over the period 1992-2011 [START_REF] Pagé | Nouvelles projections climatiques à échelle fine sur la France pour le 21ème siècle : les scénarii SCRATCH[END_REF]: the mean (±sd) annual temperature (10.4 °C ± 1.34), the mean temperature of the coldest month (2.38 °C ± 1.23), the mean temperature of the warmest month (18.74 °C ± 1.32), the inter-annual temperature variability (32.39 ± 6.17), the mean annual rainfall (1253 mm ± 1.55), the mean rainfall of the driest month (32.33 mm ± 4.57), the mean rainfall of the wettest month (214.4 mm ± 2.41) and the inter-annual rainfall variability (359.50 ± 9.76). As they were quite strongly correlated (|r| P 0.72), we used the first axis of a principal component analysis (PCA) which explained 93.80% of the variation of the climatic covariates as a synthetic covariate (CLIMATE). Its values increased as mean annual rainfall increased while mean annual temperature decreased. For both FLOW and CLIMATE, linear and quadratic terms were included in occupancy models. Third, we calculated the number of tributaries (TRIBUTARIES, derived from CARTHAGE Ó DB, mean = 3.11 ± 1.99) for each 1-km section and its proximal upstream and downstream sections. Finally, the influence of the three main hydrographic sub-sectors was also tested (SUB-SEC-TOR; 3 categories; CARTHAGE Ó DB; Fig. 1). To improve the convergence of occupancy models, all non-categorical covariates were log-transformed and normalized.

Data analysis

We applied to our dataset the Markovian occupancy model recently developed by [START_REF] Hines | Tigers on trails: occupancy modeling for cluster sampling[END_REF] which estimated four parameters: p, the probability of detecting the species conditional on the presence of the species on the site (i.e. probability of detection); w, the probability that a site is occupied or used by a species (i.e. probability of site occupancy); h 0 , the probability that a species is present on a segment given that the site is occupied but the species was absent on the previous adjacent segment; and h 1 , the probability that a species is present on the segment given that the site is occupied and that the species was present on the previous adjacent segment (see [START_REF] Charbonnel | Spatial replicates as an alternative to temporal replicates for occupancy modelling when surveys are based on linear features of the landscape[END_REF][START_REF] Hines | Tigers on trails: occupancy modeling for cluster sampling[END_REF] for more details about the parameters estimated).

To evaluate the predictive accuracy of our final model predictions, we used a jackknife iterative procedure, best suited for small data sets (n = 131) than the traditional partitioning in 70% for calibration and 30% for validation. Thus, the model selection process described below was repeated 131 times. This model selection process initially focused on determining a suitable covariates model structure for detection (p) and subsequently used this model structure to test combinations of covariates for occupancy (w). We first defined a full occupancy model (i.e. including all four occupancy covariates) based on the recommendations of [START_REF] Burnham | Model selection and multi-model inference: a practical information-theoretic approach[END_REF] and [START_REF] Mackenzie | Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence[END_REF]. Then, we formulated covariates for detection, either without any covariates, or individually or in additive combination, restricting models to a maximum of two covariates to reduce convergence problems [START_REF] Burnham | Model selection and multi-model inference: a practical information-theoretic approach[END_REF]. This resulted in eleven different models. All model comparisons were based on Akaike Information Criterion (AIC) values [START_REF] Burnham | Model selection and multi-model inference: a practical information-theoretic approach[END_REF]. The Akaike weights (xi) were also calculated. To assess the relative importance of each detection covariate, the sum of xi of models (Rxi) that included each covariate was calculated. As this process was repeated 131 times (see above), we then summed the Akaike weights obtained for the 131 iterations. Based on this global Akaike weight, the detection covariates were ranked and the ones with the highest rank were selected to fix the model structure for p. Thereafter, we kept this model structure component unchanged and ran further occupancy analyses to compare models involving either no covariates or all combinations of covariates for w, resulting in sixteen models. No covariates were included for local occupancy parameters h 0 and h 1 to reduce the number of parameters to estimate. Like for detection covariates, the global Akaike weight of each occupancy covariate was calculated. For each iteration, model averaging was used to determine the effect size (b regression coefficient) of each covariate across the top set of models (DAIC 6 2; [START_REF] Burnham | Model selection and multi-model inference: a practical information-theoretic approach[END_REF]. Finally, a global average coefficient was computed for each covariate across the 131 iterations, and used to build covariate response curves for occupancy and detec-tion. Ninety-five percent confidence intervals were calculated through the jackknife procedure. Estimates obtained through model averaging for each iteration were then used to predict occupancy probabilities over the whole Salat catchment. A final prediction map was thus produced through the computation of average probabilities across the 131 iterations. The predictive accuracy was evaluated using the area under the ROC curve (AUC) which is an index of classification accuracy independent of species prevalence and arbitrary threshold effects [START_REF] Manel | Evaluating presence absence models in ecology: the need to account for prevalence[END_REF]. We fitted all the models using the freeware PRESENCE v. 6.2 [START_REF] Hines | PRESENCE 2: software to estimate patch occupancy and related parameters[END_REF] and R.2.14.1.

Results

Desman detection

A total of 579 faeces were collected from 94 out of the 131 sampled sites. Sixty-nine percent of these faeces were genetically confirmed to be Pyrenean desman, sampled in 54 sites (i.e. 41% of sampling sites with at least one detection event during the survey). DNA included in 27% of the collected faeces was too degraded to accurately identify the author of the faeces. Four percent of the 579 collected faeces were allocated to another species (e.g. 2% to Neomys spp.). Among the 54 sites where the presence of desman was genetically confirmed, 13 had only one 100-m segment with detection, 11 had two segments with detection, 9 had three segments with detection, 12 had four segments with detection and 9 had all the five segments with detection. Eighty-seven percent of sites with detection were located in the Salat sub-sector, 13% in the Lez while no presence was recorded in the Baup sub-sector (Fig. 1).

SWAT flow simulation

SWAT simulations of stream flow were accurate, as indicated by the model evaluation statistics computed between measured and simulated monthly stream flow at the gauging station used for calibration (rho = 0.89; R 2 = 0.78; NSE = 0.73; ''Saint Lizier''; Figs. 2 and3a). According to [START_REF] Moriasi | Model evaluation guidelines for systematic quantification of accuracy in watershed simulations[END_REF], the evaluation statistics were also considered to be high at the gauging stations used for the validation step (0.85 6 rho 6 0.91; 0.62 6 R 2 6 0.85; 0.60 6 NSE 6 0.8; Fig. 3c-e), except the Balaguères station (rho = 0.81; R 2 = 0.50; NSE = 0.33; Fig. 3b). Simulated mean monthly flow over the 1992-2011 period ranges spatially from 0.01 to 31.96 m 3 /s, with a mean of 0.78 m 3 /s (Fig. 2).

Influence of covariates on detection and occupancy

The covariate that influenced desman detection most was FLOW VAR (global Akaike weight = 130.56) followed by SUB-STRATE (global Akaike weight = 70.26). OBS had a moderate influence (global Akaike weight = 44.83) while RAIN (global Akaike weight = 0.60) did not explain desman detection at all. To avoid statistical convergence issues when building occupancy models, only the first two covariates, FLOW VAR and SUBSTRATE, were retained. The probability of detecting Pyrenean desman faeces decreased sharply with increasing flow variability (Fig. 4a) and was higher in streams dominated by heterogeneous substrates and shelters in spite of a large variability in detection probabilities for stream reaches with low heterogeneity (Fig. 4b).

For desman occupancy, the covariate SUB-SECTOR exerted the strongest influence (global Akaike weight = 131.00). FLOW was also important (global Akaike weight = 88.99) followed by TRIBU-TARIES (global Akaike weight = 76.20) and then CLIMATE (global Akaike weight = 57.07). The Pyrenean desman occupancy was very different between sub-sectors with the highest occupancy probabilities for stream sections located in the Salat sub-sector (Figs. 1 and4c). As no detection event has been reported in the Baup sub-sector (Fig. 1), occupancy probability in this area was estimated to be null by the model. Both covariates FLOW and TRIBU-TARIES had a positive influence on occupancy probability (Fig. 4d ande) suggesting that the Pyrenean desman has a higher occupancy probability in stream sections with high mean monthly flow and several tributaries. For example, in the Salat sub-sector, occupancy probability raises rapidly over 0.9 from 5 m 3 /s flow values. Finally, it appears that the Pyrenean desman occupancy was higher in areas with more abundant annual rainfall and colder annual temperature (i.e. high values of CLIMATE; Fig. 4f), although this covariate was the least influencing.

Expected local spatial dependence was highlighted by the average model estimates that showed that the probability of Pyrenean desman sign presence on a segment, given absence on the previous segment (h 0 = 0.48, 95% CI = 0.46-0.49) was lower than the probability of Pyrenean desman sign presence on a segment given presence on the previous segment (h 1 = 0.72, 95% CI = 0.71-0.73).

Predictive occupancy map

The average predicted occupancy probability for the Pyrenean desman across the study area had a moderate accuracy compared to observations given that AUC value was 0.74. Occupancy probability estimates ranged from 0 to 0.97 suggesting that some streams are very suitable for the Pyrenean desman while others are not suitable at all in the upper Salat catchment. There was a strong contrast among occupancy probabilities predicted on sections of the three sub-sectors, with a mean occupancy of 0.63 ± 0.18 in the Salat, 0.15 ± 0.12 in the Lez and 0.00 ± 0.00 in the Baup (Figs. 1, 4c and5a). Higher occupancy probabilities were predicted for major rivers of the Lez and Salat sub-sectors while lower occupancy probabilities were predicted for small tributaries (Fig. 5a), underlining the positive relationship with stream flow.

These predictions indicate a potential linear distribution of 176 km (12.5%) with predicted occupancy probability P0.80, and 462 km (32.7%) with predicted probabilities P0.60 in the whole upper Salat river catchment. Last, the area with the highest variability in predicted occupancy across the 131 iterations was located on the Lez river, upstream of the Lez and the Salat confluence, and also in the headwaters of the Lez sub-sector (Fig. 5b).

Discussion

Detection probability of the Pyrenean desman -implications for monitoring

We have emphasized that the probability of detecting desman faeces decreases in areas with high flow variability, which may regularly submerge emergent items where the desman usually leaves its faeces and thus limit the accumulation of signs. In agreement with [START_REF] Aymerich | La prospección de excrementos como metodología para el estudio de la distribución de los musgaños[END_REF], we thus suggest that periods of heavy rain or high water flow as well as the period shortly following flood events should be avoided when survey involves faeces detection. We also showed that the detection probability rises in streams with a larger diversity of substrate and shelters. This result is consistent with [START_REF] Nores | Aproximación a la metodología y estudio del área de distribución, estatus de población y selección de hábitat del desmán (Galemys pyrenaicus) en la Península Ibérica[END_REF] who suggested that the availability of emergent items may influence the faeces detectability of the Pyrenean desman. When streams with a homogeneous substrate and shelters or high flow variability (natural or artificial) have to be monitored, more efforts in sampling should therefore be directed towards compensating the lowest detection, or other survey methods with higher detection efficiency (e.g. live trapping) should be applied [START_REF] González-Esteban | A comparison of methodologies used in the detection of the Pyrenean desman Galemys pyrenaicus (E. Geoffroy, 1811)[END_REF].

Variables influencing occupancy of the Pyrenean desman

We found that stream flow has a strong and positive influence on the occupancy of the Pyrenean desman. Some authors also 5 1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 1992 1994 1996 1998 2000 2002 2004 2006 suggested a preference of this species for fast flowing waters [START_REF] Nores | Aproximación a la metodología y estudio del área de distribución, estatus de población y selección de hábitat del desmán (Galemys pyrenaicus) en la Península Ibérica[END_REF][START_REF] Queiroz | The small hydro plants: predicted impacts on the Pyrenean desman populations (Galemys pyrenaicus, Geoffroy)[END_REF][START_REF] Queiroz | Status and conservation of Desmaninae in Europe[END_REF][START_REF] Ramalhinho | Habitat of the Pyrenean Desman: assessment of running water quality. Monitoring pollution[END_REF]. One hypothesis to explain this preference for high water flow in mountainous regions would be that abundance and richness of invertebrates, which are the main prey of the Pyrenean desman [START_REF] Bertrand | Répartition géographique et écologie alimentaire du desman des Pyrénées Galemys pyrenaicus (Geoffroy, 1811) dans les Pyrénées françaises[END_REF], are known to be positively related to flow [START_REF] Aymerich | Els micromamıfers semiaquatics d'Andorra: distribucio i estat de conservacio[END_REF][START_REF] Dewson | A review of the consequences of decreased flow for instream habitat and macroinvertebrates[END_REF]. This finding is also in accordance with [START_REF] Nores | Aproximación a la metodología y estudio del área de distribución, estatus de población y selección de hábitat del desmán (Galemys pyrenaicus) en la Península Ibérica[END_REF] who suggested that the factors that influence the Pyrenean desman the most are those related to hydrology characteristics. They proposed a rather mechanical explanation for this preference arguing that high water velocity may help the desman float in the water for optimal use of its resources. [START_REF] Nores | Aproximación a la metodología y estudio del área de distribución, estatus de población y selección de hábitat del desmán (Galemys pyrenaicus) en la Península Ibérica[END_REF] found that the presence of the Pyrenean desman was favoured for a water velocity greater than 0.2 m/ s. In our study, highest occupancy was obtained for river transects with flow above 5 m 3 /s and 15 m 3 /s in the Salat and Lez subsectors, respectively. But, we can observe that for these two subsectors, the most rapid increase of occupancy is obtained up to 5 m 3 /s, and grows beyond this value much more slowly, meaning that this species can tolerate quite high flow. On a larger scale (e.g. Iberian Peninsula), [START_REF] Barbosa | Otter (Lutra lutra) distribution modeling at two resolution scales suited to conservation planning in the Iberian Peninsula[END_REF] and [START_REF] Morueta-Holme | Climate change risks and conservation implications for a threatened small-range mammal species[END_REF] highlighted the importance of water abundance, owing to the amphibious lifestyle of the Pyrenean desman, as more water offers better habitat conditions and greater food availability. As suggested by [START_REF] Aymerich | Els micromamıfers semiaquatics d'Andorra: distribucio i estat de conservacio[END_REF], artificial reduction of flow regime (e.g. due to hydropower production or irrigation) likely has a negative impact on Pyrenean desman occurrence, especially in mountain rivers which are more sensitive to flow reduction. Our results also revealed the positive influence of the number of tributaries on the presence of the Pyrenean desman. A high number of tributaries may act as refuge areas in case of natural or artificial flooding [START_REF] Lake | Disturbance, patchiness, and diversity in streams[END_REF], which can facilitate the recolonization of disturbed sites by the species. The strong influence of the geographic area that we found may indicate that Pyrenean desman occurrences show spatial clustering, with an influence of terrestrial barriers, and that dispersal is mainly constrained by aquatic environments [START_REF] Filipe | Biogeography of Iberian freshwater fishes revisited: the roles of historical versus contemporary constraints[END_REF]. Spatial structuring of distribution for other semi-aquatic species less restricted to aquatic habitats, such as the European otter (Lutra lutra), has already been emphasized and is thought to result from population dynamics [START_REF] Barbosa | Otter (Lutra lutra) distribution modeling at two resolution scales suited to conservation planning in the Iberian Peninsula[END_REF]. We never found Pyrenean desman faeces in the Baup catchment, which we assume to be caused by extremely low sign detection or possible extirpation of desman from this area due to unsuitable environmental conditions within this catchment (i.e. agricultural land, drier and warmer climate). The spatial clustering of Pyrenean desman occurrences and the high probability of suitable habitats in the Salat sub-sector both suggest that this area may represent one large population connected by movement and gene flow. Stream connectivity should thus be favoured as habitat fragmentation could potentially decrease demographic support for tributaries and lead to isolation of populations, as mentioned by [START_REF] Nores | Informe sobre la situación del Desmán Ibérico (Galemys pyrenaicus) en España. Seminariode conservación de argaritifera margaritifera y Galemys pyrenaicus en la Península Ibérica[END_REF] and [START_REF] Queiroz | Status and conservation of Desmaninae in Europe[END_REF]. Until more data on movement and dispersal of the Pyrenean desman are available, the most cautious conservation strategy is thus to protect habitat quality and connectivity throughout the entire sub-sector. The predicted occupancy map would therefore be useful as a decision-making instrument for future developments within this catchment. Lastly, studies previously conducted on a larger scale revealed the importance of climate to explain the Pyrenean desman distribution. Indeed, it has been shown that mean annual precipitation has a positive influence and annual air temperature a negative influence in Spain and Portugal [START_REF] Barbosa | Transferability of environmental favourability models in geographic space: The case of the Iberian desman (Galemys pyrenaicus) in Portugal and Spain[END_REF][START_REF] Barbosa | Use of coarse-resolution models of species' distributions to guide local conservation inferences[END_REF][START_REF] Morueta-Holme | Climate change risks and conservation implications for a threatened small-range mammal species[END_REF]. By contrast, a negative relationship between precipitation and the presence of Pyrenean desman as well as a positive relationship with air temperature were found in France [START_REF] Williams-Tripp | Modeling rare species distribution at the edge: the case for the vulnerable endemic Pyrenean desman in France[END_REF]. At a finer scale (i.e. catchment scale), we found a pattern similar to that reported in Spain and Portugal, but our results suggest that climatic variables appear to be much less influential than the three other covariates considered. This contrasted result is consistent with the fact that hydrology is known to be more important in structuring the distribution of aquatic species within the species range than climate which drives spatial distribution at larger scales [START_REF] Bucklin | Comparing species distribution models constructed with different subsets of environmental predictors[END_REF][START_REF] Morueta-Holme | Climate change risks and conservation implications for a threatened small-range mammal species[END_REF].

Better accounting for the hydrology when modelling the distribution of aquatic species

In aquatic ecosystems, abiotic factors that structure freshwater species distribution are different from those that influence species in terrestrial environments [START_REF] Jähnig | Modelling of riverine ecosystems by integrating models: conceptual approach, a case study and research agenda[END_REF]. However, hydrological variables are rarely taken into account into SDMs due to a lack of detailed data for the whole study area. Applying SDMs at fine scales (e.g. catchment scale) may result in an incomplete description of species ecological niche [START_REF] Ottaviani | Modelling semi-aquatic vertebrates' distribution at the drainage basin scale: the case of the otter Lutra lutra in Italy[END_REF] as these factors are known to strongly influence the distribution of freshwater species. Hence, as mentioned by [START_REF] Jähnig | Modelling of riverine ecosystems by integrating models: conceptual approach, a case study and research agenda[END_REF], integrated modelling approaches are needed to provide habitat suitability predictions of aquatic organisms using adequate predictors. Some hydrological models such as the SWAT model make it possible to simulate stream flow across an entire stream network and then to use flow as an input variable in SDMs. For example, by using hydrological and hydraulic models, [START_REF] Jähnig | Modelling of riverine ecosystems by integrating models: conceptual approach, a case study and research agenda[END_REF] simulated water levels, flow velocities and sediment processes in a river section, to serve as inputs into an SDM applied to a freshwater mollusc. All these hydrological and hydraulic parameters appeared to explain significantly its occurrence. Moreover, [START_REF] Kuemmerlen | Integrating catchment properties in small scale species distribution models of stream macroinvertebrates[END_REF] used SWAT to simulate hydrological factors (e.g. flow seasonality, flow range) and applied SDMs to invertebrate species at the catchment scale. They showed more accurate predictions of habitat suitability when hydrological factors were included among the predictors. Coupling between a hydrological model and SDMs thus appears promising to get more accurate predictions of species distribution from SDMs in freshwater ecosystems at fine resolution [START_REF] Kuemmerlen | Integrating catchment properties in small scale species distribution models of stream macroinvertebrates[END_REF]. This should be applied to a wider range of aquatic and semi-aquatic species (e.g. mammals, birds, reptiles, amphibians, fish, stream invertebrates), especially for studies at local scales to improve ecological knowledge about these species.

Accounting for the imperfect detection of species

A major issue in monitoring wild animal populations is the difficulty to detect species. This problem can be taken into account using occupancy models to avoid false absence data [START_REF] Mackenzie | Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence[END_REF]. The value of these statistical methods has been largely acknowledged during the last few years with many applications for a wide diversity of species such as insects [START_REF] Kéry | Predicting species distributions from checklist data using site-occupancy models[END_REF], fish [START_REF] Comte | Species distribution modelling and imperfect detection: comparing occupancy versus consensus methods[END_REF], reptiles [START_REF] Durso | Needles in haystacks: estimating detection probability and occupancy of rare and cryptic snakes[END_REF] and birds [START_REF] Kéry | Analysing and mapping species range dynamics using occupancy models[END_REF]. In our study, in spite of quite high detection probabilities found for the Pyrenean desman, the use of a site-occupancy model has limited the risk of underestimating the Pyrenean desman occupancy. More generally, occupancy models may avoid introducing errors into distribution patterns [START_REF] Comte | Species distribution modelling and imperfect detection: comparing occupancy versus consensus methods[END_REF] or introducing bias into wildlifehabitat relationships [START_REF] Gu | Absent or undetected? Effects of non-detection of species occurrence on wildlife-habitat models[END_REF][START_REF] Kéry | Predicting species distributions from checklist data using site-occupancy models[END_REF], two issues that are crucial for species with low detectability. Indeed, for rare, secretive or nocturnal species, it may be difficult to detect individuals visually and logistics constraints may appear for surveys in rough terrain [START_REF] Aing | A Bayesian hierarchical occupancy model for track surveys conducted in a series of linear, spatially correlated, sites[END_REF]. Consequently, survey methods based on recording indirect signs have become standard practice for many species [START_REF] Heinemeyer | Natural sign: tracks and scat[END_REF]. Depending on the species ecology, if the data are based on sign surveys and sampled following a design with spatially adjacent replicates, they could be analyzed with the Markovian occupancy model. This may concern species that use trails to move as do many large mammals [START_REF] Barber-Meyer | Influence of prey depletion and human disturbance on tiger occupancy in Nepal[END_REF][START_REF] Karanth | Monitoring carnivore populations at the landscape scale: occupancy modelling of tigers from sign surveys[END_REF][START_REF] Thorn | Brown hyaenas on roads: Estimating carnivore occupancy and abundance using spatially autocorrelated sign survey replicates[END_REF], using a linear home range (e.g. species living along rivers or ecotones) or species for which surveys follow a linear configuration (e.g. roads, trails) such as those carried out in remote systems (e.g. mountains, forests). Moreover, spatially adjacent replicates have the advantage of being more effective and straightforward to implement than temporal replication (see [START_REF] Charbonnel | Spatial replicates as an alternative to temporal replicates for occupancy modelling when surveys are based on linear features of the landscape[END_REF] for more details).

Whereas the importance of accounting for false negative errors is frequently recognized, much less attention has however been given to false positives [START_REF] Miller | Improving occupancy estimation when two types of observational error occur: non-detection and species misidentification[END_REF]. A specific occupancy model developed by [START_REF] Royle | Generalized site occupancy models allowing for false positive and false negative errors[END_REF] can account for false positives when genetic analyses are too expensive and not appropriate for some survey techniques (e.g. listening points). In our case, DNA analyses were appropriate to avoid false positives that would have occurred by misidentifying the faeces of other species (e.g. Neomys spp., Glis spp., Myotis spp, Turdus spp., Pordarcis spp.) as being Pyrenean desman. For species with shrinking ranges, false-positive observations may thus result in an underestimation of the population decline. However, the use of modern DNA techniques makes sure that a faeces sample belongs to the species of interest, thus resulting in a reduced risk of overestimating occupancy probabilities. Given the serious consequences of inaccurate estimates of the status of rare species for conservation and management decisions, accounting for false-positive errors should be an important component when designing and analyzing monitoring programs for rare species [START_REF] Miller | Improving occupancy estimation when two types of observational error occur: non-detection and species misidentification[END_REF] or species for which surveys are based on the recording of indirect cues (e.g. scats) like many carnivore species [START_REF] Heinemeyer | Natural sign: tracks and scat[END_REF].

Fig. 1 .

 1 Fig. 1. Map of the upper Salat river catchment (study area) with hydrographic sub-sectors (light grey: Lez, medium grey: Salat and dark grey: Baup) and sampling site locations (dots, n = 131). The size of dots indicates the number of segments with detection for each site.

Fig. 2 .

 2 Fig. 2. Map of SWAT sub-basin distributions, hydropower reservoirs (black squares) and gauging stations (red triangles) used to calibrate and validate SWAT modelling. The average simulated stream flow (1992-2011) is also shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .

 3 Fig.3. Measured (grey) and simulated (black) Soil and Water Assessment Tool (SWAT) stream flow at the downstream gauging station used for the calibration step (a) and at gauging stations used for the validation step (b-e; see Fig.1for the location). Model simulations were evaluated with the Spearman coefficient correlations (rho), coefficient of determination (R 2 ), and Nash-Sutcliffe efficiency (NSE) between measured and simulated stream flow.

Fig. 4 .

 4 Fig. 4. Relationship between detection (a, and b) and occupancy (c-f) probability, and covariates for the Pyrenean desman in the upper Salat river catchment. All the figures show the average probabilities predicted along the range of each covariate (solid black lines) and their 95% confidence intervals (dotted lines), except Fig.4bwhich plots the estimated detection probabilities in the 131 sites according to substrate heterogeneity categories, and Fig.4cwhich plots the predicted occupancy probability for the whole study area (1388 1-km long sections) according to the sub-sector. Predictions were computed with the other covariates at their mean values.

Fig. 5 .

 5 Fig. 5. Maps of (a) the average predicted occupancy probability of the Pyrenean desman in the upper Salat river catchment across the 131 iterations and (b) the prediction variability (i.e. coefficient of variation) across the iterative jackknife procedure. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Acknowledgements

We thank all the field investigators who helped for the field work: R. Lassus, F. Julien, F. Gilbert, S. Perré and A. Denis. We are also grateful to SWAT modellers for their precious help: Y. Grusson, A. Uhart, X. Sun, J. Payoux, as well as the R modeller A. Maire, and J. Brown and M. Williams-Tripp for the implementation of the GRTS sampling. Thank to S. Danflous who corrected the English text, and the Parc Naturel Régional des Pyrénées ariégeoises for the availability of the land use database (BDOS Ó 2013). This study was funded by ANRT (Cifre n°2011/1018), EDF (Electricité de France) and European Union (FEDER) and is part of the French conservation Action Plan for the Pyrenean desman (2010-2015) supervised by DREAL MP (Direction Régionale pour l'Environnement, l'Aménagement et le Logement de Midi-Pyrénées) and coordinated by the CEN MP (Conservatoire d'Espaces Naturels Midi-Pyrénées). Appendix A. Supplementary material

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.biocon.2015.01. 019.


Ludovic Dutreve
email: ludovic.dutreve@liris.cnrs.fr

Alexandre Meyer
email: alexandre.meyer@liris.cnrs.fr

Saïda Bouakaz
email: saida.bouakaz@liris.cnrs.fr

Feature Points Based Facial Animation Retargeting

Keywords: CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Animation, I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling, Facial Animation, Retargeting, Performance-Driven Facial Animation, Skeleton-Subspace Deformation

Figure 1: The bottom row shows the virtual face animated by retargeting expressions from the source face (top row).

Introduction 1.Motivation

Facial animation is an important aspect in 3D environment featuring human characters and is highly demanded by many kinds of applications, such as 3D games, interactive human/computer softwares and movies. However, a good facial animation requires a lot of time for a skilled artist: the complexity of the human faces may provide subtlety expressions, and each spectator is able to know instinctively if an animation is correct or not. Since pioneering work of [START_REF] Parke | Computer generated animation of faces[END_REF]], many efforts were provided to increase realism and satisfy the critical eyes of humans. Nevertheless, two aspects of facial animation are still stimulating active research: the process of rigging a face and the transfer of an expression from one person to another.

The process of rigging is tedious and requires many hours of manual work to an artist. Indeed, even simple method like shapes blending needs an artist to create the key shapes. The common skinning approach needs to manually define the joints (also called bones in the case of body) of the skeleton such as left of the right eyebrow, top of the lip, etc.on the mesh. And, then the artist needs to specify which parts of the surface are attached to which joint(s). More sophisticated models based on simulation of the facial tissues and muscles need also human interventions to attach tissue to skeleton and to tune complex muscles configuration. Furthermore, the animation produced for a specific face cannot be reused directly to another.

The second active aspect of facial animation domain is the interesting idea to produce realistic animation by transferring an existing one to a new character, instead of generating it from scratch. Performance-Driven Facial Animation systems introduced the idea of capturing animation of a real actor, often from 2D video(s) and transferring it to a virtual character. The transfer can also be applied from an existing 3D animation created by an artist. Another advantage of the transfer is to easily manage the amount of memory to store on sensitive applications like video games by storing only one expression for many face models. Despite all these advantages, the transfer suffers from an inherent difficult: the target model should still look like the original person, but smoothly reflect the source expression. The techniques of transfer have to adapt the source deformation to the target morphology.

Related Work

Many approaches exist to animate a face in real-time. Most common are linear weighting with blendshape models [START_REF] Deng | Animating blendshape faces by cross-mapping motion capture data[END_REF], skeleton-based with linear blend skinning [START_REF] Magnenat-Thalmann | Joint-dependent local deformations for hand animation and object grasping[END_REF][START_REF] Lewis | Pose space deformation: a unified approach to shape interpolation and skeleton-driven deformation[END_REF], physics-based muscles models [Kähler et al. 2001;[START_REF] Zacur | Transferring a labeled generic rig to animate face models[END_REF]], etc. Blendshape consists on an artist to create the key shapes and then linearly blend between those shapes to obtain a fluid motion. However, obtaining an exhaustive set of key shapes is time and memory consuming, as the human expression is very rich. Linear blend skinning (LBS), also known as skeleton subspace is the process of binding the face mesh to a skeleton. Each vertex is assigned a set of joint influences and their associated weights. To deform the mesh, each vertex is rigidly transformed by all of its influences and the blending weights are used to compute the weighted average of these positions. Despite more recent and sophisticated approaches such as muscle-based approaches, skinning deformation remains the most popular method used in practice for real-time animations of characters. Mostly, because of its simplicity and efficiency (with simple GPU implementation). The process of obtaining influences can be initiated by an automatic procedure based on the vertex proximity to the bone, as the one proposed by Maya or Blender. Because this automatic approach lead to artifacts, accurate influences have to be tuned manually. Notice that for body, the Pinocchio system [Baran and Popović 2007] proposed last year an automatic rigging of 3D characters based on heat simulation. And dedicated to professional animators, [START_REF] Zacur | Transferring a labeled generic rig to animate face models[END_REF]] proposed a transfer of their complex rigging system based on muscle from a face to another. For all the advantages of linear blend skinning mentioned before, we chose to use it to deform our face and we improved it by a rapid and near automatic influences computation.

After Gleicher published in [START_REF] Gleicher | Retargetting motion to new characters[END_REF]] the idea of reusing an animation in the case of body, and since animating a face from scratch needs hard manual efforts, many authors proposed retargeting techniques [START_REF] Deng | Animating blendshape faces by cross-mapping motion capture data[END_REF][START_REF] Vlasic | Face transfer with multilinear models[END_REF][START_REF] Na | Hierarchical retargetting of fine facial motions[END_REF][START_REF] Pyun | An example-based approach for facial expression cloning[END_REF]]. [START_REF] Chai | Vision-based control of 3d facial animation[END_REF]] presented an interesting technique consisting on using low-quality vision based 2D motion data to animate a face with high-quality 3D precomputed motion capture data. The Radial Basis Functions (RBF) [START_REF] Powell | Radial basis functions for multivariable interpolation: a review[END_REF]] were used to adapt motion vectors of a mesh to an other [START_REF] Noh | Expression cloning[END_REF]. Motion vectors are the displacements of vertices along time. By putting some landmarks on the target and the source faces, a data interpolation provided by RBF allows to compute a dense matching between target vertices and source triangles. Vectors are modified in scale and orientation to adapt with the target morphology. But some limitations occur: the target animation is precomputed and cannot be used in a runtime application, vertices motion vectors for each frame require a high memory space. Furthermore, the whole expression is retargeted, without any possibility of tuning. And, the need of source motion vectors prevents the use of data from motion capture. Our approach may be seen as an extension to this work in order to accept any kind of source animation, to make the setup more attractive to novice user by providing an almost automatic rigging and by using the skinning technique to deform the target mesh which is more common on 3D engine.

A more complete state of the art of the facial animation retargeting problem was proposed by [START_REF] Pighin | Facial motion retargeting[END_REF]]. We entrust [START_REF] Deng | Computer facial animation: A survey. Data-Driven 3D Facial Animations[END_REF] to present a more general survey about the facial animation.

Overview

We present a method to transform an animation, obtained from an animated 3D face or a recorded 2D video, to another 3D face. An overview of our technique is illustrated in Figure 2. After initializing the source and the target faces by manually select between 15 or 30 feature points, an automatic procedural rigging is applied on the target mesh. Then, any deformation/animation applied to the source is transferred to the target mesh using our RBF approach applied on feature points. Finally, deformation of the target mesh is deduced from the transferred features points by linear blend skinning on the GPU. Thus, all the process of transfer and mesh animation is real-time. Notice that the manual work needed by our system is tractable to a novice user. Indeed, a simple mesh can be loaded, quickly rigged by positioning the feature points and then be animated in real-time according to the source animation, which can be a 2D video motion capture. To summarize, the main contributions of our system are:

• Quasi-automatic rigging of the facial mesh in neutral position;

• Real-time retargeting;

• Partial retargeting: user could choose to partially transfer the facial animation by selecting feature points;

• Compatibility with various animation sources: 2D video or 3D motions of feature points. The remainder of this paper is organized as follows. In Section 2, we present our general usage of radial basis function in the context of space transformation. In Section 3, we describe our mesh deformation model and our automatic rigging. In Section 4, we present the retargeting method based on a RBF network and how we use it to adapt a motion from a face to an other. In Section 5, we present some results and conclude this paper.

RBF Space Transformation

Radial Basis Functions (RBF) are widely used in computer graphics for surface approximation and interpolation, but we may use them as a space transformation. In order to do this, we define two spaces with two set of feature points. Let S0 (resp. T0) be the set of source feature points (resp. target feature points) and N the size of these sets (see Figure 2). Each control points si ∈ S0 has its dual ti ∈ T0.

After the training of the RBF with these two sets of feature points, the RBF can transform a position from the source space to the target space.

The RBF equation is:

F ( sj) = N i=1 wi.h (|| sj -si||) (1)
with h (|| sj -si||) = (|| sj -si||) 2 + sc 2 j (multi-quadric function) and scj = min j =i (|| sj -si||). Training the network consists on solving the 3 linear systems of N equations (in the 3 dimensional case) such as tj = F ( sj)

(2)

Let H the matrix such as Hij = h (|| sj -si||) and Tx = ( t x 1 t x 2 ...t x N ) t , t x j is the x coordinate of tj. Using equations 1 and 2, the system can be defined by

Tx = H.Wx (3)
with the weights Wx = ( w x 1 w x 2 ...w x N ) t . Then, to solve the system, we have to compute Wx = H -1 Tx. Once the RBF network is trained for each axe, the position in the target space t for each point s of the source space is obtained by applying the transformation F ( s) provided by Equation 1.

Target Mesh Deformation

Animating a mesh involves using a mesh deformation technique. In our context, this deformation has to verify two constraints. First, it has to be computed easily in real-time. Second, it has to be guided by the deformations of a small number of feature points because our method aims to retarget animations providing by any kind of sources: 2D/3D tracking, 3D animation. To meet these constraints, we chose to use a skeleton-based deformation: linear blend skinning widely used by 3D engine for animating characters.

Facial Skinning

Linear blend skinning deforms a mesh according local transformations of the skeletal bones. In order to adapt skinning deformation technique with feature points displacements, we consider feature points as joints (bones). We define weights for each vertex such as wi = 1 with the i th bone's weight wi ∈ [0, 1]. A weight of 0 means that the feature point does not influence the current vertex, while a value's weight of 1 means that the current vertex is influenced only by this feature point. The position of the vertex v at frame f is defined by

v f = v0 + N i=1 wi. d if N
is the number of feature points, and d if is the displacement of the i th feature point at frame f . We can notice that for body animation, transformations are defined by only joint's rotations mainly because of the limb length constant constraint. However, faces have specific small deformations which can be described by joints translations.

Quasi Automatic Weighting

To tackle the tedious process of influences assignment in skinning, we propose a quasi automatic approach. The only manual work required is positioning the feature points on the target face. Our solution is procedural. The main idea is to compute for each vertex an influence value (also called a weight) to a feature point according to their distance. This weight function has to respect some criteria like smooth borders and smooth gradient from 1 when the vertex is on the joint to 0 when the vertex is far from the joint. We base our weight function on the ellipsoid function which allows to define independent parameters for each axe since feature points influence area is not regular. Let vi a vertex of the mesh, tj a feature point position and d = vi -tj = (dx, dy, dz) t . For each feature points, rx, ry and rz are the radius of the ellipsoid for the 3 axes which are defined for each part or the face using feature points relationship. For example, the top borderline of the eyes is large compared to its height, so its influence is mainly spread over the width. For this case, rx depends on distance between eye's corners, and ry depends on the distance with the eyebrow. vi is then influenced by the feature point tj as follow:

wij = max 0, 1 -d 2 x /r 2 x + d 2 y /r 2 y + d 2 z /r 2 z .
For each vertex, we compute a similar weight for each feature points. Then, we consider only the maximum weights with significant contribution and normalize them in order to have a sum of weights of 1. This simple method provides an automatic rigging. However, we have to admit that this technique can lead to lack of precision on some specific kind of face like non human faces. Thus, we offer to the user the possibility to modify and adapt the rigging if necessary by a painting process.

Animation Retargeting

In this section, we describe our animation transfer technique.The morphologies between two human faces can greatly vary, and much more when the 3D characters are cartoons, monsters or animals. Therefore, motions of some feature points on a source face along an animation can not be used directly. They have to be adapted in scale and in orientation. The RBF space transformation provides a good solution to this problem, feature points of the source face define the source space while feature points of the target face define the target space.

The original animation may come from motion capture or landmarks (labelled vertices) on an animated 3D face. In the case of video tracking, only displacements on the x and y directions are known. We assume that it is sufficient to describe a facial animation (human eye does not need 3D informations to recognize an actor's expression). Thus, we set the z coordinates of feature points to 0 when source animation is in 2D. Each feature points on the target model is tracked or followed in the source data. The retargeting method we describe now is illustrated in Figure 2. As transfer preprocess (usually at first frame), a RBF is trained between the 3D source feature points and the target feature points, both need to be in neutral expression. Let define FRBF = RBF (S0, T0), with S0 (resp. T0) the source (resp. target) feature points set (with a neutral expression). Then, at each frame f of the animation, coordinates of tracked feature points are adapted to the target virtual face by T f = FRBF (S f ). Thus, target feature points move according to the expression of the source face and the morphology of the target face. The target mesh is then deformed on the GPU using the linear blend skinning as described in Section 3.

Our retargeting provides the capacity of transferring only a part of an animation as illustrated by the Figure 3. Indeed, for some reasons, animators may want to transfer only eyes or mouth motions etc.Thus, he can disable feature points not owned by these areas and then the retargeting works only with active parts of the face. Notice that this feature may be difficult to realize with a global approach. Figure 1 and4 show transfer between virtual heads: 21 feature points are used, 3 on each eyebrow, 2 on each eye, 3 on the nose and 8 around the mouth. The transfer and target animation is applied at screen frame rate on a current laptop (2,2GHz dual core, 2GB of ram and a nVidia GeForce 8600M GT). Indeed, the time of the RBF-training preprocessing is negligible since the main computation time is due to the inversion of the N × N matrix where N is the number of feature points, usually about 25. The mesh dimension has almost no impact on the speed of transfer, since vertices displacements are computed by skinning on GPU.

To test our method in an interactive situation, we developed a tracking application based on the Pyramidal Lucas Kanade Feature Tracker [START_REF] Bouguet | Pyramidal implementation of the lucas kanade feature tracker: Description of the algorithm[END_REF]]. 24 white markers on the actor's face. To make the tracking more robust by increasing the contrast of markers, we draw black circles on them. We use a common commercial webcam with a resolution of 640x480 pixels. Our system tracks the 2D facial movements of the markers and sends, at each frame, their positions to the retargeting viewer. Our system runs at 30 frames per second which is the camera frame rate, on a single PC or on two PC sending tracked position through the network.

Figure 1 shows the result of a retargeting between motion capture data and a virtual face.

Though our technique requires few manual intervention, we plan to work on a full-automatic detection of feature points on 3D virtual faces and to improve the quality of our automatic rigging [START_REF] Zacur | Transferring a labeled generic rig to animate face models[END_REF]]. We also would like to remove the need of marker in the video tracker in order to provide a complete helpful system for animators in their animations prototyping tasks or for novice users in playful applications. An other improvement of our approach would be to capture and transfer fine details such as wrinkles and small deformations of skin.

Figure 2 :

 2 Figure 2: Overview of our retargeting method.

Figure 3 :

 3 Figure 3: Example of partial transfer. The left image shows the source expression. The second presents a full retargeting. The two last faces show partial transfer: chin and mouth points, eyes and eyebrows points.

Figure 4 :

 4 Figure 4: The top row shows the source face's expressions. The left column shows source and target faces with a neutral expression. Expressions of the source face are retarget to the two other faces.

http://www.playall.fr

Acknowledgements

This work has been supported by the PlayAll 1 consortium, a Cap Digital project.