Thibaut Marin-Cudraz
email: thibaut.marincudraz@gmail.com

Bertrand Muffat-Joly
email: bertrand.muffat-joly@oncfs.gouv.fr

Claude Novoa
email: claude.novoa@oncfs.gouv.fr

Philippe Aubry
email: philippe.aubry@oncfs.gouv.fr

Jean-François Desmet
email: jfdesmet@orange.fr

Mathieu Mahamoud-Issa
email: mathieumahamoudissa@hotmail.fr

Florence Nicolè
email: florence.nicole@univ-st-etienne.fr

Mark H Van Niekerk
email: markh.vn@gmail.com

Nicolas Mathevon
email: mathevon@univ-st-etienne.fr

Frédéric Sèbe
email: frederic.sebe@univ-st-etienne.fr

Acoustic monitoring of rock ptarmigan: A multi-year comparison with point-count protocol

Keywords: acoustic monitoring, rock ptarmigan, point-count protocol, long-term study, individual information, methodology comparison

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Acoustic monitoring is becoming an effective means to assess wildlife diversity, resulting in minimal impact to the environment [START_REF] Towsey | Ecology and acoustics at a large scale[END_REF][START_REF] Sueur | Ecoacoustics: The Ecological Investigation and Interpretation of Environmental Sound[END_REF].

Importantly, it enables to focus on species used as bioindicators as well as of patrimonial or economic concern. The use of acoustic monitoring is usually motivated by the difficulty in observing the species because of its secretive behavior or the difficulty in accessing its habitat [START_REF] Hoodless | Vocal individuality in the roding calls of Woodcock Scolopax rusticola and their use to validate a survey method[END_REF][START_REF] Vögeli | Who are we sampling? Apparent survival differs between methods in a secretive species[END_REF][START_REF] Marques | Estimating cetacean population density using fixed passive acoustic sensors: An example with Blainville's beaked whales[END_REF][START_REF] Buxton | Measuring nocturnal seabird activity and status using acoustic recording devices: applications for island restoration: Acoustic Monitoring of Nocturnal Seabirds[END_REF]Dugan et al., 2013;[START_REF] Marques | Estimating animal population density using passive acoustics: Passive acoustic density estimation[END_REF][START_REF] Andreassen | Semi-automatic long-term acoustic surveying: A case study with bats[END_REF][START_REF] Ulloa | Screening large audio datasets to determine the time and space distribution of Screaming Piha birds in a tropical forest[END_REF]. Levels of investigation range from the simple assessment of species' presence/absence to more complex studies that aim at determining the number of individuals present on an area.

Depending on the species, these approaches may require the discrimination of individuals [START_REF] Terry | The role of vocal individuality in conservation[END_REF][START_REF] Pollard | Pre-screening acoustic and other natural signatures for use in noninvasive individual identification[END_REF], which is only possible when vocalizations contain individual vocal signatures based on morpho-physical, genetics and/or learning abilities [START_REF] Kroodsma | Production, perception, and design features of sounds, Acoustic communication in birds[END_REF][START_REF] Tibbetts | Individual recognition: it is good to be different[END_REF][START_REF] Catchpole | Bird song: biological themes and variations[END_REF][START_REF] Taylor | The contribution of source-filter theory to mammal vocal communication research: Advances in vocal communication research[END_REF][START_REF] Nowicki | The evolution of vocal learning[END_REF][START_REF] Tamura | Mating calls are a sensitive indicator of phylogenetic relationships in tropical tree squirrels (Callosciurus spp.)[END_REF].

Although acoustic monitoring is promising, it still raises several potential issues: high cost of monitoring material, design of monitoring protocols, analysis of long-duration recorded signals, weather conditions impacting the quality of the data, needs of bioacoustics experts for data analysis, and sensitivity to the density of the species [START_REF] Budka | Is it possible to acoustically identify individuals within a population?[END_REF]Linhart and Šȧlek, 2017). To the best of our knowledge, acoustic monitoring has not yet replaced other protocols. Point-counts protocols [START_REF] Lancia | Estimating the number of animals in wildlife populations[END_REF] are still largely used to provide estimates of the number of individuals. However, these conventional, human-operated methods are exposed to biases due to inter-individual differences between observers, unpredictability of field conditions, and biological parameters such as species abundance [START_REF] Tyre | Improving precision and reducing bias in biological surveys: estimating false-negative error rates[END_REF][START_REF] Bart | Density estimation in wildlife surveys[END_REF][START_REF] Lotz | Observer Bias in Anuran Call Surveys[END_REF][START_REF] Elphick | How you count counts: the importance of methods research in applied ecology[END_REF][START_REF] Fitzpatrick | Observer bias and the detection of low-density populations[END_REF].

Moreover, individuals of species emitting long-range signals can be counted by several observers simultaneously, leading to double counting and abundance overestimation.

Besides, human presence can disturb birds' behavior making them stop singing and leading to population underestimation.

While previous works have explored the technical feasibility of acoustic monitoring based on vocal individual signature [START_REF] Terry | Census and monitoring based on individually identifiable vocalizations: the role of neural networks[END_REF][START_REF] Hartwig | Individual acoustic identification as a non-invasive conservation tool: an approach to the conservation of the african wild dog Lycaon pictus (temminck, 1820)[END_REF][START_REF] Grava | Individual acoustic monitoring of the European Eagle Owl Bubo bubo[END_REF][START_REF] Policht | Potential for the use of vocal individuality as a conservation research tool in two threatened Philippine hornbill species, the Visayan Hornbill and the Rufous-headed Hornbill[END_REF][START_REF] Adi | Acoustic censusing using automatic vocalization classification and identity recognition[END_REF][START_REF] Feng | Individuality and Stability in Male Songs of Cao Vit Gibbons (Nomascus nasutus) with Potential to Monitor Population Dynamics[END_REF][START_REF] Budka | Is it possible to acoustically identify individuals within a population?[END_REF][START_REF] Budka | Vocal individuality in drumming in great spotted woodpecker-A biological perspective and implications for conservation[END_REF]Peri, 2018a) most field applications were based on sounds recorded from already known individuals [START_REF] O'farrell | A Comparison of Acoustic Versus Capture Techniques for the Inventory of Bats[END_REF][START_REF] Peake | Corncrake Crex crex census estimates: a conservation application of vocal individuality[END_REF][START_REF] Vögeli | Who are we sampling? Apparent survival differs between methods in a secretive species[END_REF][START_REF] Digby | A practical comparison of manual and autonomous methods for acoustic monitoring[END_REF][START_REF] Peri | Censusing a tawny owl (Strix aluco) population living at high density merging two consolidated techniques[END_REF]. To the best of our knowledge, there is no published study investigating the generalization and reliability of an acoustic monitoring approach based on individual vocal signatures aimed at estimating the number of individuals in real field conditions.

Focusing on the rock ptarmigan Lagopus muta in France, the present paper reports a multiyear study aimed at assessing the potential interests, as well as the caveats, of acoustic monitoring by comparing performances of this approach with a traditional point-count methodology.

The rock ptarmigan is a bird species that inhabits the northern parts of Eurasia and North America. In France, its range is restricted to sub-alpine and alpine habitats (altitude > 1800 m) of the Pyrenees and the Alps Mountain ranges [START_REF] Sale | Grouse Of The World[END_REF]. This species is secretive and difficult to access. It is highly adapted to its environment, has mimetic plumage and vocalizes at dawn and dusk [START_REF] Macdonald | The breeding behavior of the Rock Ptarmigan[END_REF]. Population abundances are decreasing in both the Alps [START_REF] Imperio | Climate Change and Human Disturbance Can Lead to Local Extinction of Alpine Rock Ptarmigan: New Insight from the Western Italian Alps[END_REF][START_REF] Furrer | Variable decline of Alpine Rock Ptarmigan (Lagopus muta helvetica) in Switzerland between regions and sites[END_REF][START_REF] Martinoli | Where is the pulse to have the finger on? A retrospective analysis of two decades of Alpine Galliforms (Aves: Galliformes) census and game bag data in Italy[END_REF] and the Pyrenees, where they are threatened by both climate change and habitat transformation [START_REF] Revermann | Habitat at the mountain tops: how long can Rock Ptarmigan (Lagopus muta helvetica) survive rapid climate change in the Swiss Alps? A multi-scale approach[END_REF][START_REF] Bech | Pyrenean ptarmigans decline under climatic and human influences through the Holocene[END_REF]. Due to these extreme environmental and selective pressures, rock ptarmigan is often considered as a bioindicator of the ecosystem health, a sentinel and umbrella species for biodiversity conservation of the alpine environments [START_REF] Sandercock | Demographic consequences of agestructure in extreme environments: population models for arctic and alpine ptarmigan[END_REF][START_REF] Hanser | Greater sage-grouse as an umbrella species for shrubland passerine birds: a multiscale assessment[END_REF][START_REF] Henden | Changed Arctic-alpine food web interactions under rapid climate warming: implication for ptarmigan research[END_REF].

During the mating season (May-June), males display courtship rituals, which often includes simple, pulsatile vocalizations as well as a peculiar "singing in flight" behavior before dawn [START_REF] Macdonald | The breeding behavior of the Rock Ptarmigan[END_REF]. The flight is hyperbolic, and vocalization starts when the bird reaches the point of highest altitude. Point-count protocols rely on this acoustic behavior to evaluate the number of males after they have established their breeding territories [START_REF] Bossert | Bestandesaufnahme am Alpenschneehuhn Lagopus mutus im Aletschgebiet[END_REF].

Although vocalizations are loud and easily heard by an observer, low visibility, birds' mobility, frequent harsh weather conditions and other constraints associated with the alpine environment are likely to increase the possibility of double counting and overestimation of the monitored population. Conversely, this may also impair song perception by observers [START_REF] Andreev | Ecological characteristics of Lagopus mutus on the commander islands[END_REF]. The use of an acoustic monitoring technique could provide a feasible alternative as a response to these difficulties and potential biases. Despite these constraints, rock ptarmigan presents several advantages in terms of monitoring. Vocalizations are easily recognizable, population densities are generally low (around 5 males/km 2 ), and males are mostly located in stable territories which facilitates their localization and recording.

The present study proposes [START_REF] Fraley | Model-Based Clustering, Discriminant Analysis, and Density Estimation[END_REF] to assess the individual vocal signature embedded in rock ptarmigan males' calls, [START_REF] Bech | Pyrenean ptarmigans decline under climatic and human influences through the Holocene[END_REF] to test the feasibility and reliability of a bioacoustic monitoring approach over several years, and (3) to compare the results obtained using this approach to those obtained with a traditional point-count protocol and long-term field observations.

Material and Methods

We performed this study in 2015, 2016 and 2017, at the ski resort of Flaine (French Alps, Haute-Savoie, 45°59'32.8"N 6°43'44.2"E; altitude: 1600-2500 m).

Counting of singing males

We used three different methods to evaluate males' abundance: a point-count protocol, longterm observations and acoustic monitoring. The latter was followed by signal processing and statistical analysis.

Point-count protocol

The counting area was delimited empirically to cover roughly 100 hectares (Fig. 1). Previous field observations and literature reviews had suggested that male ptarmigans could be heard at distances of 500-1000 m [START_REF] Watson | A population study of ptarmigan (Lagopus mutus) in Scotland[END_REF][START_REF] Bossert | Bestandesaufnahme am Alpenschneehuhn Lagopus mutus im Aletschgebiet[END_REF][START_REF] Marty | Point-count method for estimating rock ptarmigan spring density in the Pyrenean chain[END_REF].

We selected three counting points, approximately 500 m apart from one another to maximize the listening area (Fig. 1). These counting points were located at precisely the same positions over the three-year period of study. The counting procedure followed that proposed by [START_REF] Léonard | Méthode de dénombrement des Lagopèdes alpins mâles au chant et présentation des résultats[END_REF], which was based on the work of [START_REF] Bossert | Bestandesaufnahme am Alpenschneehuhn Lagopus mutus im Aletschgebiet[END_REF]. Three experienced observers, with previous knowledge of the field (TMC, FS, BMJ; 1 per counting point), accompanied by volunteers, were placed at each point. As ptarmigan males mainly display their acoustic signal early in the morning, the observation periods started at 4.30 a.m. and ended at 5.30 a.m. The observers were positioned 15 min before the beginning of the observation period. Throughout the observation period, each observer noted on an observation form (with a drawn map of the area) the timing, the number and the approximate estimated locations of the vocalizing ptarmigans. At the end of the observation period, BMJ collected all observers' forms.

To estimate the total number of males in the area, we counted 1 male for each group of neighboring vocalizations indicated on the maps (each group had to be clearly separated from the others; see Fig. 1). Each male was confirmed by cross-checking the observers' data. To avoid double-counting by two different observers, vocalizations localized nearby and heard at approximately the same time (interval < 20 sec.), were considered identical. We considered the total number of males unambiguously localized as the minimum total number of males. To take into account ambiguous localizations (e.g. when one observer noted two birds on the same location while another observer noted a single one) we defined a maximum total number of birds by adding these ambiguous males to the minimum total number. We thus defined an estimated interval (min-max) of the total number of males. This procedure was repeated on several consecutive days namely: June 3rd and 4th in 2015 and 2016; June 6th, 7th, 9th in 2017.

Long-term observations and acoustic recordings

After the point-counting days, two bioacousticians (TMC and FS) remained in the field for a period of one month to observe and record all the birds present in the area (recording material: Sennheiser MKH70 shotgun microphones connected to Marantz PMD 660 recorders; sampling rate: 48000 Hz). This was an important long-term observation effort aimed at ensuring a thorough knowledge and identification of each bird present within the area, which may have been acoustically recorded during the point-count protocol. The resulting number of males spotted by this technique therefore represented the expected maximum value that could be deduced using the other methods.

Despite this effort however, we were only able to record birds during 7, 9, and 15 days in 2015, 2016 and 2017 respectively, primarily due to the harsh weather conditions. We used two different strategies in order to build up our bank of recordings:

1) Recording of non-identified birds before sunrise (4.30 a.m. -6 a.m.), i.e. during the time slot corresponding to the spontaneous singing activity peak. Due to the low ambient luminosity, the recorded males could not be visually observed and identified. The recording data sets were thus named: "unknown datasets". Each day, both bioacousticians recorded from different locations within the study area to sample a maximum number of males.

2) Recording of identified birds (6 a.m. -around 10 a.m.). In 2017, we equipped 5 males present on the area with GPS solar tags (e-obs GmbH, Grünwald, Germany), and used these tags to pinpoint the males with their individual UHF (Ultra High Frequency) radiofrequencies. One male had already been equipped with a VHF (Very High Frequency) radio-emitter collar since 2015. Two additional males were identified using visual cues only.

Both had mated with females and remained within stable and well-defined territories. The combination of direct observations and GPS localizations greatly minimized the risk of misidentifying or confounding males during recordings. A ninth male had its territory on a neighboring summit (3 km away). As this male was alone on this site and easy to identify, it was thus added to the recording database (total = 9 birds in 2017).

The recording strategy was as follows: after sunrise, when the peak of males' vocal activity ended, we played back calls from an individual recorded in another area to elicit the focus male's territorial response. This allowed us to record each focus male while double-checking for its individual identity.

In the preceding years (2015 and 2016), we used the same approach although the results of the field effort were weaker:

-June 2016: 5 males recorded. Two of the males were equipped with GPS tags and 3 other males were identified using visual observations only.

-June 2015: 7 males recorded. One male was equipped with a VHF radiotransmitter necklace; 6 males were identified using visual observations only.

The 2015, 2016 and 2017 recording data sets obtained with this method were named "known datasets".

Acoustic analysis of recorded signals

Data bank of calls

Rock ptarmigan vocalizations are sequences of pulse trains, with energy spread over a frequency spectrum ranging from 900 to 3700 Hz. There are two major types of calls, namely short and long calls. These differ by the number of successive pulse trains, namely 3 and 4 respectively [START_REF] Macdonald | The breeding behavior of the Rock Ptarmigan[END_REF][START_REF] Watson | A population study of ptarmigan (Lagopus mutus) in Scotland[END_REF]. For the present study, we focused on the short calls, which are the most frequently recorded (Fig. 2). Our annual data bank of calls consisted of the following:

-2015: 183 short calls, of which 100 were of sufficient quality (in terms of signal-tonoise ratio) to be analyzed ("unknown" dataset: 75 calls; "known" dataset: 25 calls, 3.6 ± 2.4 calls/male, min = 1 call, max = 8 calls).

-2016: 249 short calls, with 98 of sufficient quality ("unknown" dataset: 66; "known" dataset: 32 calls, 6.4 ± 2.3 calls/male, min = 4 calls, max = 10 calls).

-2017: 180 short calls, with 133 of sufficient quality ("unknown" dataset: 52 calls; "known" dataset: 81 calls, 7.1 ± 3.9 calls/male, min = 3 calls, max = 24 calls).

Automatic detection of group of pulses

Due to harsh weather conditions (wind and rain), recorded signals were frequently corrupted by noise. Before performing the automatic detection of pulses, we first filtered the signals with a 950-2700 Hz bandpass filter, and then processed a wavelet continuous transform (WaveleComp R package, [START_REF] Roesch | WaveletComp: Computational Wavelet Analysis[END_REF] to optimize the signal-to-noise ratio (see Supplementary Material for details).

After denoising, amplitude pulses were detected using a customized script (Seewave R package, [START_REF] Sueur | Seewave: a free modular tool for sound analysis and synthesis[END_REF]. The absolute amplitude of the signal was first smoothed using a Daniell kernel (length = 100). The time position of the pulses was then identified using an amplitude threshold fixed at 3 % of the maximum amplitude of the considered signal.

Pulses were gathered in "groups" (G1, G2, G3, see Fig. 2) by automatically measuring intervals between pulses (Fig. 2c) and computing the ratios between two successive intervals (Fig. 2d). Ratios superior to 1.7 characterized boundaries between groups of pulses (Sil1 and Sil2; see Fig. 2b).

Measurement of acoustic parameters

We measured 12 acoustic parameters from groups of pulses G1 and G2 (see list in Table 1).

We chose to ignore the G3 group, as the signal-to-noise ratio of this part of the call was usually very weak (this was mostly due to males flying away from the recorder while singing).

The mean acceleration (Acc.G1) was calculated as follows:

= 1 -2 1 - 1 -
Where n is number of pulses in G1; Pr is the pulse rate and t is the time of occurrence of the pulse's maximum amplitude.

The normalized Pairwise Variability Index (nPVI) is an index commonly used in phonetic studies [START_REF] Grabe | Durational variability in speech and the rhythm class hypothesis[END_REF] to illustrate the variability between consecutive pairs of intervals:

= 100 × ∑ | ( ! ! ) # | .
We used continuous wavelet transformation to calculate the peak frequency parameters (Fq1.G1 and Fq2.G1). Since pulse locations in the signal had already been calculated, each pulse was isolated from the original sound, filtered with an 800 Hz -3000 Hz passband filter.

For each pulse, wavelet power spectrum was then calculated. The two scales with the maximum power were then selected. The scales were further multiplied by the Fourier factor 6 / (2π) to obtain the classical Fourier periods (Aguiar-Conraria and Soares, 2011) with the corresponding frequencies. The median of each peak frequency was then calculated across the entire pulse train. Overall, medians were preferred to means since they are more conservative and less influenced by possible outliers resulting from analytical errors (e.g. due to rain drops occurring within a pulse group).

Acoustic space and individual vocal signatures

In order to build a functional tool allowing the acoustic discrimination between rock ptarmigan males, we proceeded as follows: 1) we built a 4-D acoustic space with a dedicated discriminant analysis which optimized separation between males using the 2017 "known" dataset (reference dataset); 2) we used this acoustic space to perform an unsupervised clustering analysis using the reference dataset for tuning the clustering hyperparameters; 3) we applied the workflow (centering, projection in the acoustic space and then tuned clustering) on the 2015, 2016 and 2017 complete data sets to further estimate the number of males present each year.

Acoustic space definition

We analyzed the differences between calls from the individuals of the 2017 "known" dataset (9 identified males) using powered partial least squares discriminant analysis (PPLS-DA, [START_REF] Liland | Powered partial least squares discriminant analysis[END_REF]; "pls" R package, [START_REF] Mevik | Partial Least Squares and Principal Component Regression[END_REF]. PPLS-DA enables more accurate analysis of a small sample size with a high number of acoustic parameters, better than the linear discriminant analysis signatures [START_REF] Hervé | Multivariate Analysis of Multiple Datasets: a Practical Guide for Chemical Ecology[END_REF], commonly used in studies on animal vocalizations. PPLS-DA calculates new variables as combinations of all centered acoustic variables, leading to a new acoustic space optimizing the discrimination between individuals. The number of dimensions was chosen by model cross validation [START_REF] Szymańska | Double-check: validation of diagnostic statistics for PLS-DA models in metabolomics studies[END_REF]. The mean classification error rate was established after 100 model cross validations for each number of dimensions (varying between 2 to 11, Fig. 3). We followed an analog method of the elbow method [START_REF] Cattell | The Scree Test For The Number Of Factors[END_REF] to assess the optimal number of dimensions and we considered the first four PPLS-DA variables, which define a 4-D acoustic space and explained 36.52 % of the variability (see Results and Fig. 4).

We tested the statistical significance of our PPLS-DA model with a procedure implemented in the RVAideMemoire package [START_REF] Westerhuis | Assessment of PLSDA cross validation[END_REF][START_REF] Hervé | RVAideMemoire: Testing and Plotting Procedures for Biostatistics[END_REF]. The PPLS-DA significance validation is composed of two steps. Firstly, a set of discriminant functions is obtained from a training data set and secondly, these functions are used to test the classification on a validation set. The measure of standard error is obtained by analyzing the correct assignment percentage of 999 random selections of the original data set, which have been divided into a fitting and testing set.

Clustering analysis and bootstrap reliability testing

We used an unsupervised classification method (high dimensional data clustering, HDDC, [START_REF] Bouveyron | High-dimensional data clustering[END_REF] to estimate the number of individual males present within the datasets. HDDC has already been applied on acoustic data with some success [START_REF] Ulloa | Estimating animal acoustic diversity in tropical environments using unsupervised multiresolution analysis[END_REF].

HDDA is known to be consistent and reliable with unbalanced datasets because it is based on gaussian mixture models [START_REF] Fraley | Model-Based Clustering, Discriminant Analysis, and Density Estimation[END_REF]. It is more parsimonious and flexible than gaussian mixture modelling by adding a noise term within the model covariance parametrization. The mixture model aims at identifying the meaningful variables for each cluster and is fitted with the E-M algorithm. The number of mixture components of the model maximizing the Bayesian information Criterion (BIC, Shwarz, 1978) is set as the number of clusters. The E-M algorithm is sensible to the selected random points during its initialization. Thus, we ran the clustering algorithm several times in order to obtain a reliable value for the number of clusters.

The 2017 'known' reference dataset was used for tuning the HDDC hyperparameters (K = 9 clusters in 2017; the covariance model M; the threshold t used to parametrize the dimension of each cluster; see R package HDclassif, Bergé et al., 2012 for details). Each call was represented by its 4 acoustic dimensions previously calculated through PPLS-DA. We tested 10 values of t namely: 0.000001, 0.00001, 0.0001, 0.01, 0.03, 0.05, 0.07, 0.1, 0.15, 0.2 (adapted from [START_REF] Ulloa | Estimating animal acoustic diversity in tropical environments using unsupervised multiresolution analysis[END_REF] and the 14 possible models of covariance parametrization. Each association of t and M value were tested.

The clustering algorithm was run 100 times for each association. For each run, we measured the similarity between the clustering output and the clustering membership with the adjusted Rand Index (ARI, [START_REF] Hubert | Comparing partitions[END_REF]package mclust, Scrucca et al., 2016). The ARI ranges from -1 to 1 and is an indicator of the concordance of two classifications for the same dataset: when ARI = -1, the classifications are totally opposed, or different. When ARI = 0, the classifications are considered random; when ARI = 1, they are identical. The mean ARI was then calculated for the 100 values and the tuning parameters associated with the highest mean were selected. The maximum mean adjusted Rand Index (ARI = 0.91) was found for the simplest covariance model ("abqd") and a threshold value of 0.1 was assigned.

We thus used these tuning parameters.

Once M and t fitted with our data, the reliability of the clustering process was further tested using sub-sets of the 2017 reference 'known' dataset. Subsets were built by randomly sampling several various males (random sampling of 1 to 9 males; 900 subsets; 100 trials/subset) or a various total number of calls (random sampling of 20 to 81 calls; 6200 subsets; 100 trials/subset). We tested models with K values ranging from 1 to 20.

Comparison between counting methods

To assess the number of males through the acoustic analysis of calls, we performed the clustering analysis on each year separately, using the entire datasets obtained by pooling "known" and "unknown" calls (100 trials/year). The male of 2017 that was geographically isolated was however excluded because it was located outside the point-counting area. We calculated the 4 dimensions of each call using the PPLS-DA functions previously defined with the 2017 "known" dataset for each year separately. Each acoustic dataset was centered before its projection in the 2017 acoustic space by subtracting the means of each acoustic variable calculated on the 2017 "known" dataset.

The number of males (i.e. the number of acoustic clusters) estimated for each year through the clustering analysis was then compared with the number of males estimated through the two other counting methods, i.e. 1) the point-count protocol and 2) the long-term observation.

All the acoustic and statistical analysis was performed on R version 3.5.0 (R core team, 2018).

Results

Acoustic discrimination between males and definition of the acoustic space

The PPLS-DA identified significant acoustic differences between ptarmigan males, with 4 significant functions that allowed maximizing individual separation (Fig. 4, mean classification rate = 79 %; min-max classification rate per individual = 0 -100 %; p = 0.001). Table 2 shows the correlation between each of the 4 components, the acoustic variables and shows that all parameter types (pulse number, pulse rate, durations, frequency parameters) contribute towards separating the males.

Reliability of High Dimensional Data Clustering

To test the reliability of HDDC, we compared the median number of males obtained through resampled HDDC with the actual number of males of each sub-data set. As displayed in Fig. 5, HDDC gives a reliable estimate of the number of recorded males if this number does not exceed 5 individuals. HDDC underestimates the number of recorded males when 6 or more individuals were included in the sub-dataset. It also underestimates the number of recorded males, when the number of vocalizations in the sub-datasets are sampled, and consistently predicts 7 clusters (i.e. 7 individuals) for sub-datasets composed of a minimum of 33 vocalizations (i.e. 41 % of the total number of calls) (Fig. 5b).

Comparison of counting methods' reliability

In 2017 and 2016, the most congruent results were given by the acoustic monitoring and long-term observation. In both years, the point-count protocol resulted in a lower estimation than the two other counting methods. Still, the long-term results were reached by the pointcount intervals for at least one day per year. Estimation through the point-count protocol appears to be highly dependent on the day of observation (this is particularly obvious in 2017, with an estimate of 5-8 males on the first day versus 4-5 males on the second day).

The 2015 results differed significantly from those of 2016 and 2017, with an apparent underestimation of the number of males through the acoustic monitoring method compared to long-term observations. However, the distribution is widespread and looks bimodal, with the second mode (6 clusters) being close to the actual number of males (7 individuals). This can be clearly seen in Fig. 6 which displays the number of males estimated by each counting method (point-count protocol, long-term observation, acoustic monitoring).

Discussion

Does the acoustic space built from recordings encompass the vocal variability of rock ptarmigan males?

The relative inconsistency of individual males' vocal signature might limit the bioacoustics approach. Although the mean PPLS-DA classification rate of recorded calls was around 80%, individual rates differed greatly among males (from 0% to 100%). Moreover, the discriminant functions used to build the acoustic space explained only 36% of the total acoustic variability of calls. A significant proportion of the variability remains out of reach, suggesting that individual identity is not the only factor driving the calls' structure. Rock ptarmigan' vocalizations are sequences of stereotyped pulses with few frequency modulations. The acoustic variation between individuals may thus be reduced when compared to other bird species with more complex signals. Ptarmigan are non-oscine birds [START_REF] Kroodsma | Acoustic Communication in Birds: Song learning and its consequences[END_REF][START_REF] Slater | Bird song learning: causes and consequences[END_REF], and their vocalizations thus lack then the individual variability that could have been induced by song learning. Inter-individual differences in ptarmigan acoustic signals mostly result from differences in their genetic background and their physiological conditions. It is known that ptarmigan males are highly philopatric and closely related genetically in the Alps at large scale [START_REF] Caizergues | Population genetic structure of rock ptarmigan Lagopus mutus in Northern and Western Europe[END_REF]. In addition to this, a genetic study of a closely related species, the red grouse (Lagopus lagopus scoticus), showed that males were highly related at local scales [START_REF] Piertney | Local genetic structure in red grouse (Lagopus lagopus scoticus): evidence from microsatellite DNA markers[END_REF]. The genetic variability between rock ptarmigan males is thus rather low. Moreover, the highly variable alpine weather conditions should promote great annual variations in food availability, especially due to snow cover and the timing of snow melt [START_REF] Körner | Alpine Plant Life: Functional Plant Ecology of High Mountain Ecosystems[END_REF][START_REF] Edwards | Changes in the seasonal snow cover of alpine regions and its effect on soil processes: A review[END_REF]Jonas, 2008). Thus, males' physiological state might be different both between individuals (e.g. depending on the individual food intake in each territory) and from year to year within individuals (depending on the availability of resources). The variability from year to year is of special concern as vocalizations from the same male could be very different each year, thus impairing recapturing males over consecutive years using acoustics only. For this reason, we considered each year independently within the present study.

Our analysis demonstrates that, despite these potential limitations and thanks to their low densities in France (few males present on the same area), it is possible to rely on individual vocal signatures to identify local rock ptarmigan males (Linhart and Šȧlek, 2017) and, ultimately, to count them. The first requirement is to include only high-quality recordings in the analysis step (recordings are regularly corrupted by background noise, mostly induced by wind). Besides, analysis should be mainly performed on the temporal acoustics parameters, since these are less influenced by noise than the frequency cues. These conservative choices and the fact that, by design, PPLS-DA optimizes the separation between males and not the explained variability, can partly explain the low percentage found. Nevertheless, we were still able to separate the males well enough for our purpose.

High dimensional data clustering

A second potential limitation of the bioacoustics method may arise in cases where some males are represented by only a few recordings, resulting in unbalanced recording datasets.

HDDC is a model-based clustering, fitted by maximizing log-likelihood estimation (MLE), itself based on the probabilities of clusters membership. The likelihood will tend to favour clusters forming large and homogenous groups of points (Fraley and Raftery, 1998). When the recording dataset is strongly unbalanced between males, individuals with few vocalizations can be confounded (i.e. included in the same cluster) with other males that show close vocalization characteristics. One of the males was overrepresented in the 2017 "known" dataset and its calls formed a cluster that incorporated vocalizations from some other males. This resulted in an underestimated total number of individuals. After removing this male from the dataset (for comparing between methods), the number of estimated clusters was higher and more reliable. Unbalanced situations are more likely to arise when the recording effort is mitigated, as was the case in 2015. This caused HDDC to underperform, resulting in an underestimation of the number of males present in the observation area. Such a difference can be explained when looking at the call's clusters obtained through the HDDC method for each year, using both 'known' and 'unknown' data sets (Fig. 7). In 2015, the vocalizations were mostly from unknown emitters. As a result, clusters strongly overlapped each other. Conversely, 2016 and 2017 vocalizations are more clearly separated into homogenous groups. The recording sampling effort must then be homogeneous over the whole study area and cautiously planned to optimize the representativity of recording banks [START_REF] Heupel | Automated acoustic tracking of aquatic animals: scales, design and deployment of listening station arrays[END_REF].

Sampling effort and balance

An adequate acoustic dataset needs a major field effort, due to frequent harsh weather conditions and difficulties associated with approaching males' territories. Despite these constraints, comparable results were found between bioacoustics and long-term methods.

Long-term observations seem the most reliable approach, although not feasible on a regular schedule given the required workforce. However, most of the vocalizations were obtained using playbacks and males were recorded directly within their territories ("known" datasets).

Marginal males were included (Fig. 1) as we assumed that they could be heard and could fly inside the study area. The area of interest was therefore slightly larger than the area covered by the point-count protocol. This increased the probability of male detection during the acoustic monitoring compared to the point-count protocol. Conversely, point-count monitoring appears to be less accurate, with a greater variability of males' abundance estimations between counting days.

The reduced reliability of the point-count census is not surprising, given that observers focus generally on only one observation day. We extended the counting period for a few days to show the variability of this method in this study. Weather conditions (e.g. wind speed) and variation in males' motivation to vocalize may impair the detectability of males. Moreover, the number of males present in the area fluctuates through the season and between observation days. In contrary to northern latitudes [START_REF] Unander | Behaviour and Social Structure in Svalbard Rock Ptarmigan Lagopus mutus hyperboreus[END_REF][START_REF] Cotter | The Reproductive Biology of Rock Ptarmigan (Lagopus mutus) in the Central Canadian Arctic[END_REF], ptarmigan males differ in their arrival dates depending on their mating status. Mated males return to their territory in late summer / autumn, while single males arrive only in spring of the subsequent year. Snow cover plays an important role in the availability of territories and reproductive success of the species [START_REF] Novoa | Effect of weather on the reproductive rate of Rock Ptarmigan Lagopus muta in the eastern Pyrenees: Weather and rock ptarmigan reproductive rate[END_REF] by delaying the males' arrival (one of the 2015 males equipped with VHF was not present in 2016 -a year during which snow covered its territory-but arrived in late June in 2017, when the snow cover of its territory started to dissipate).

Point-count census is a "one-shot" process: it estimates the number of males at a given day whereas bioacoustics and long-term observation estimations are obtained over larger time periods. The point-count protocol is therefore not able to capture changes throughout the mating season. In practice, only long-term monitoring using direct observations or acoustic recordings could overcome this variability.

In addition, our study was able to test the reliability of the census method based on the number of individuals and the number of signals taken into consideration. We showed that a minimum number of sampled vocalizations were necessary (approximately 33 vocalizations for 9 individuals, Fig. 5B) to ensure a consistent assessment of males' number. Although the efficiency of the clustering method is density dependent, this is on a smaller scale than those of classical counting methods [START_REF] Budka | Is it possible to acoustically identify individuals within a population?[END_REF]Linhart and Šȧlek, 2017). To apply our analysis methodology across populations in different locations, we would probably need to consider a range of higher signal diversities. The generalization of our study would thus need to train our statistical model (PPLS-DA) with recordings from identified males, from other rock ptarmigan populations, to encompass greater signal variability and to avoid staying at local scale variability.

Is bioacoustics monitoring a good solution for rock ptarmigan population monitoring?

The choice of a monitoring method is the result of a balance between the scale of the study and the expected results. At a fixed cost, the same number of automatic recorders allow to gather precise information regarding males' abundance and locations within a restricted area or could provide simple presence/absence survey of a wider region. Besides, the

Acknowledgements

This research was supported by the University of Lyon / Saint-Etienne (PhD stipend to TMC, research support), the ONCFS (Office National de la Chasse et de la Faune Sauvage), CeLyA (Centre Lyonnais d'Acoustique) and the Institut universitaire de France (NM). The Domaine Skiable de Flaine logistically supported the study. We would like to warmly thank Fabrice Antoine, Marc Arvin-Berod, Thomas Betton, Loïc Berger, Sébastien Bernard, Michel Bouchard, Félicien Bros, Léna De Framond Benard, Etienne Marlé, Joël Prince, Aymeric Richard, Fanny Ryback, and Théophile Yeme who willingly volunteered to participate into the point-counts. We thank two anonymous reviewers for their helpful comments. The authors declare no conflict of interest.

bioacoustics approach could enable the censusing of more areas without requiring an increase of the number of observers and/or days of observation. This may be particularly advantageous in large and remote mountain massifs where rock ptarmigan habitats can be situated in remote areas.

Author contribution statement

TMC, FS, BMJ, JFD, MMI, FN, PA, CN and NM conceived the ideas and designed the methodology; BMJ captured the males, equipped the GPS collar and organized the pointcount protocols. TMC and FS collected the acoustic data; TMC, MMI, FN analyzed the data; TMC, FS and NM wrote the manuscript. All authors contributed critically to the drafts and gave final approval for publication.