Drivers of black grouse trends in the French Alps: The prevailing contribution of climate

Mountains host complex ecosystems whose wide range of ecological conditions over small geographical distances makes them biodiversity hotspots. To ensure their long‐term conservation, a better understanding of the interaction between climate change and modifications in land use is necessary. Most studies have focused on only one of these factors at a time, leading to incomplete predictions. In this study, we explored the relative contribution of both recreative activities and climate change on the population dynamics of the black grouse (Tetrao tetrix), an emblematic cold‐adapted species of the Alps.


| INTRODUC TI ON
Mountain regions cover around 25% of the Earth's surface (Kollmair et al., 2005 ). They are characterized by a high diversity of habitats at a relatively small scale, resulting from different topoclimates within narrow elevation gradients (Körner, 2004). As a consequence of such a variety of habitats, mountains are associated with high levels of species diversity, supporting one-quarter of terrestrial biodiversity and hosting nearly half of the world's biodiversity hotspots (Myers et al., 2000). At high altitude, species specialized for cold conditions inhabit small patches comparable to islands separated by unsuitable plains (Constanzi & Steifetten, 2019). As marked changes in environmental conditions occur over short distances in mountains, species tend to occupy a smaller range compared to lowland species. This characteristic restricted range, combined with isolation, makes them more prone to local extinction (Körner, 2007;Quintero & Jetz).
Compared to lowland ecosystems, mountains have been relatively spared from human activities (European Environment Agency, 2010). Nevertheless, the rapidity of current global changes, which are deeply modifying the structure and functioning of ecosystems, is expected to have major consequences on mountains in the years to come, notably in Europe (Beniston, 2003;Chapin et al., 2000). In the Alps, over the last half century, land use has dramatically changed. Agricultural activities, principally pastoral, have either strongly declined or increased, depending on a parcel's potential yield and accessibility (Fischer et al., 2008;Strebel & Bühler, 2015). In addition, recreative activities have exploded over the last three decades, with the introduction of new and ubiquitous summer and winter sports (Arlettaz et al., 2007;Jenni-Eiermann & Arlettaz, 2008). The intensity of these activities, however, varies from one massif or valley to another due to accessibility, conditions, substrate and climate (Braunisch et al., 2011).
In addition to these changing patterns of human land use and recreative activities, ongoing climate change is occurring faster in mountains than in plains and is expected to have a strong impact on mountain ecosystems (Beniston, 2003;Scridel et al., 2018). An increase in temperature is expected to lead to a shift of the upper edge of forests towards higher altitudes and a general increase in the dominance of woody deciduous shrubs at high elevations (Gehrig-Fasel et al., 2007;Walther et al., 2005). The major alterations in both the intensity and distribution of precipitation that is predicted will also strongly modify microclimates, notably through changes in snowfall and soil humidity (Walther et al., 2005). The consequences of precipitation changes on ecosystems are difficult to forecast and are likely to vary greatly locally. Nevertheless, significant changes in both vegetal and animal phenology are expected (Visser & Both, 2005).
Mountain biodiversity has already been deeply impacted by ongoing land use, recreative activities and climate change and the intensity of these anthropogenic factors are expected to increase in the years to come (Dirnböck et al., 2011;Hoorn et al., 2018). While an increasing number of studies are documenting the response of animal species to global change (Both et al., 2006;Chen et al., 2011;Stephens et al., 2016), these tend to focus on one factor at a time, leading to poor knowledge of the relative contribution of each potential factor and their interactions and thus to incomplete predictions (Harrington et al., 1999;Jetz et al., 2007). Indeed, for example, as climate change will reduce snow cover at the lowest elevation, we can expect that skiing will move to areas that remain well-preserved and might generate new disturbance (Arlettaz et al., 2017).
The aim of this study was to explore the relative contribution of global change on black grouse (Tetrao tetrix) population trends in the French Alps. The black grouse is an emblematic species of cold ecosystems, inhabiting a broad range of mountainous and boreal habitats (Cramp & Simmons, 1980;Storch, 2007). Despite its large global population, like all alpine Galliformes, the species is considered in need of conservation and is listed in Annex I and II of the EU Birds Directive (Storch, 2007). The black grouse is a particularly relevant case study for investigating the interactions and relative contribution of different global changes as the species is expected to be sensitive to pressures such as habitat fragmentation (Kurki et al., 2000;Kuki et Linden, 1995), forest maturation (Pearce-Higgins et al., 2006), either shrub encroachment of alpine pastures (Dettenmaier et al., 2017;Patthey et al., 2012) or intensification of grazing (Strebel & Bühler, 2015), climate change (Barnagaud et al., 2010;Moss, 2001;Novoa et al., 2008;Swenson et al., 1994;Wegge & Rolstad, 2017), development of recreational activities (Arlettaz et al., 2013;Pattey et al., 2008) and hunting (Zbinden et al., 2018).
Despite the value of disentangling the relative influence of these factors, no studies have simultaneously studied the impact of land use, recreative activities and climate on the population dynamics of this species, probably because like all high-altitude species, alpine populations of black grouse are often difficult to reach and thus to monitor (Imperio et al., 2013).
We based our study on data from a total of 47 monitoring sites (breeding grounds of displaying males) located throughout the French Alps. Over the monitoring period, sites were progressively added: about one-third of the sites were monitored from the 1980s onwards, the second third since the 1990s and the last third since at least 2005. We estimated site-specific population growth rates in the network of sites and used them to explore the spatial structure of trends. We then tested for the potential effect of selected climatic factors and recreative activities on both long-term population trends and inter-annual variation in population growth rates to investigate the relative contribution of these factors to the population dynamics of this species in the Alps.

| Study species and area
In the French Alps, the black grouse mainly lives in semi-open habitats close to the treeline in the subalpine zone between 1,400 and 2,300 m, especially in transition zones between forests and alpine meadows. In the northern part of the Alps, the species mainly occupies dwarf-shrub moors dominated by alpenrose (Rhododendron ferrugineum) and subalpine meadows with a few scattered conifers, while in the southern Alps it occupies more diverse habitats, from open habitats to forests dominated by larch (Larix decidua), arolla pine (Pinus cembra) or mountain pine (P. uncinata) (Caizergues, 1997). Over its lifecycle, the black grouse requires a mosaic of different habitats: it selects areas that have open, sparsely vegetated land for courtship displays, good shelter for roosting, and shrubs or trees for feeding above the snow in winter (Schweiger et al., 2012;Zeitler, 2003). Populations' dynamics are also strongly dependent on the presence of an herbaceous or ericaceous stratum for brood rearing. Various studies showed the importance of structural heterogeneity and composition of breeding habitats, with the association of dense herbaceous meadows or undergrowth provided with mosaic grasses and ericaceous shrubs for the cover and feeding of the broods, and of scattered trees and shrubs for the cover and feeding of adults (Bernard-Laurent, 1994;Braunisch et al., 2016).

| Male counts
Black grouse have a very unique and well-recorded courtship ritual that involves the males strutting around at dawn in a breeding arena, or "lek" (Glutz von Blotzheim et al., 1994). They have a highly distinctive mating call that allows relatively easy counts. In the subalpine zone, birds form leks from March to June in open microhabitats such as ridges, clearings and meadows with few or no trees (Bernard-Laurent et al., 1994;Magnani, 1988). While the Black Grouse is known to be primarily a lekking species, as described in northern Fennoscandia (Borecha et al., 2017), studies showed that solitary display in southern populations can be very common and may represent up to 50% of singing males (Chamberlain et al., 2012;Nelli et al., 2016). In addition to males that can be mobile in space-especially sub-adults-previous studies showed that lek attendance, as well as the intensity of displaying, can vary along the year and in relation to weather (Baines et al., 1996;Gregersen and Gregersen, 2014). Count of individuals has been based as much as possible on visual observations of birds.
In France, monitoring of the black grouse has been coordinated by the Observatory of Mountain Galliformes (Observatoire des Galliformes de Montagne: OGM) since 1998. They have established 53 counting sites of displaying males in accessible areas known to host black grouse leks. Site selection was based on subjective criteria, some sites have been selected because easily accessible and occupied with large leks in very good habitats, others only because they were close to observers that wanted to count leks despite being in a marginal zone, and finally some because hunters feel that a local population was decreasing so they wanted to document the trend. It seems then difficult to predict what effect this site selection could have on population trends at regional or national scale, as we suspect opposing patterns from the different situations listed above. Nevertheless, the frame of reference is not limited to the sites that are already known to be best for black grouse in the core of their range as it includes sites at the southern edges. Starting year of monitoring varied between 1979 and 2018 depending on the site. Since then, counts have been conducted annually or every two years, following a standardized protocol established for black grouse monitoring in France (described in Leonard, 1989 andMontadert, 2016). As many observers as required to cover the studied area patrol the sector walking or staying in a predetermined fixed position (depending on the site structure) and count, from a certain distance, displaying males on or around leks (Ellison & Magnani, 1985). Counts are performed from early morning, before sunrise, until no later than two hours after sunrise between the last week of April and the end of May, as close as possible to the activity peak (Montadert, 2016).
Intra-annual replicates (2 or 3) have also been performed on most sites since 1995. Counting is always performed when weather conditions are good, that is avoiding heavy rain, strong winds or dense fog (Baines, 1996;Leonard, 1989;Montadert, 2016).
For this study, we focused on the count data of a subset of 47 sites (see Figure 1. From the complete network of 53 sites, we excluded six sites as these had been monitored for too few years (<7 years). Note that 4 of these sites stopped being counted a few years after population reached zero, despite possible impact of discarding such sites when numerous, we chose not to take them into consideration in our analysis as this would need a different approach (for example occupancy models). The size of the 47 selected sites varied from 650 to 7,300 ha (mean = 1862 ha). We defined four regions based on bio-geographic characteristics: northern pre-Alps (13 sites), northern internal Alps (16 sites), southern pre-Alps (6 sites) and southern internal Alps (12 sites). We analysed the time series of the 47 sites from 1996 to 2018 (hereafter, "22-year analysis") as well as a subset of the 11 sites for which counts started the earliest, from 1979 to 2018 ("39-year analysis"). While the first analysis provided good coverage to test the effect of spatial covariates on trends, the second analysis, relying on a longer time series, provided a more precise view of the long-term trends.

| Anthropogenic and environmental variables
We assembled variables related to environmental conditions and human recreative activities that we postulated might affect black grouse population trends. All variables Table 1 & Appendix 8 were scaled (x = x-mean/sd), either on the full matrix or site by site when data were not comparable between sites. We tested for an effect of each covariate on median trends or on inter-annual variation in population growth rates or both.
The surface area and mean elevation of each site were calculated based on the contour and centroids of the monitoring site.
Site sizes were relatively homogenous, except for one site that was much larger than the others (7,000 ha). We tested for an effect of the surface area on initial population size, observation error and median trend (see "Analysis" and Table 1. Regarding elevation (median = 1780 m, range = 1065-2098 m), range at the site scale (elevation max-min) were comparable between sites, so we chose to use the mean.
Previous studies have shown that the timing of snowmelt influences the date of initiation of plant growth and consequently the body condition of grouse hens and their breeding success (Moss & Watson, 1984;Watson, 1965). We extracted snowmelt dates from the moderate resolution imaging spectroradiometer (MODIS) on board the Terra platform's Earth Observing System. MODIS data was available from 2000 onwards. We defined the date of snowmelt as the end of the period of continuous snow cover. Date_100 corresponds to total snowmelt at the site (100% bare ground), and date_50 to 50% bare ground (see (Novoa et al., 2016) for the use of such covariates on another Galliformes). In parallel, we used MODIS normalized difference vegetation index (NDVI) data with the "greenbrown" package in R (Forkel & Wutzler et al., 2015;R Core Team, 2018) to estimate two metrics that summarize vegetation phenology: start of growing season (SOS) and length of growing season (LOS). Both the start and end of the growing season were defined from the derivative of the seasonal curve, as the mid-points of spring green up and autumn senescence, respectively . We considered that MODIS data was of comparable quality across sites,hence, we tested for a linear and a quadratic effect of date_50, date_100, SOS and LOS on both median trends and inter-annual variation in population growth rates.
We based our weather indices on meteorological parameters estimated by SAFRAN (Système d'analyse fournissant des renseignements atmosphériques à la neige), a mesoscale atmospheric analysis system for surface variables using ground data observations from Météo France (Durand et al., 1993). While the temporal resolution of SAFRAN is very good, with data available from 1958 onwards, its spatial resolution of 8 km was not well suited to our small mountainous sites. Thus, we chose to define only two synthetic variables, the first representing the summer temperature (mean temperature in June, July and August) and the second summarizing precipitation, calculated as the sum of monthly precipitation in June, July and August divided by the number of days with available data. While weather covariate values might be biased due to cell resolution, bias was considered constant over time,hence, we tested for a linear and a quadratic effect of temperature and precipitation only on the interannual variation in growth rates.
We extracted yearly ski-related cable density as a proxy of winter sports activity (public data from OGM: https://obser vatoi re-galli forme s-monta gne.com). We used a 2-km buffer zone around monitoring sites, defined to correspond to the black grouse's altitudinal range (1400-2300 m). In France, hunting of adult cocks is practised according to regional and local regulations, which vary by locality and have changed over time. Hunting bag sizes have been available since 1998, coinciding with the establishment of hunting guidelines for black grouse in France (Com. pers. OGM). The number of cocks killed per year has progressively decreased ever since, from 1,000 birds in the late 1990s to about 400 birds in recent years. The annual harvest rate is now thought to represent about 5% of the estimated French black grouse population, with rates that can rise to 10% in some locations (Com. pers. OGM). We gathered hunting bag data from departmental hunting organizations. Data collection of hunting bags started on different dates depending on the site (between 1994 and 2006). For 17 sites, hunting bag data matched monitoring sites exactly (in terms of geographical area). For these sites we created a variable for harvest density (number of birds killed per ha) and tested its effect on both median trends and the inter-annual variation in growth rates. For the remaining 19 sites for which we had hunting data, the surface area on which hunting bags were collected could differ to a certain extent from that of the monitoring site, but remained constant over the study period. We then built a second variable for the 36 sites (for which we had data on the number of harvested birds) and tested its effect only on the inter-annual variation in growth rates.

F I G U R E 1
Map showing the 47 black grouse monitoring sites in the French Alps used in the 22-year trend analysis (1996-2018) (the sites are shown with their OGM reference number). Symbols represent the four regions (triangle = northern pre-Alps, circle = northern internal Alps, star = southern pre-Alps, square = southern internal Alps). The white symbols correspond to the subset of 11 sites used for the 39-year time series . Elevation is shown in grey-scale, and the grey lines delineate the French departments

| Statistical analysis
We analysed population trends using state-space models (Buckland et al., 2004;de Valpine et al., 2002;Dennis et al., 2006) within the Bayesian framework (Kéry & Schaub, 2012), based on the work of Furrer et al. (2016) on the rock ptarmigan. This type of model allowed us to separate the dynamics of the state over time (the abundance) from the observation process (the counts). As our populations are expected to be mainly declining and to occur at low densities (BirdLife International, 2016), we assumed stochastic exponential population growth and we chose not to consider density dependence. Population growth rate is the logarithm of the ratio of two successive population sizes, hence the size of the monitored site (which was constant over the years) had no impact on the population growth rate.
We built a general model on which we applied several decompositions to address different questions. In all cases, the state process of N i,t+1 the population size of site i in year t + 1 was defined as follows: with r i,t the annual population growth rate of site i from year t to t + 1 and i,t , reflecting the variability in the population growth rates around the median trend r i defined by a stochastic process based on a Normal distribution with 2 proc the residual variance. Then, the observation process was defined as a log-normal abundance distribution, with y i,t,k the k counts on site i at year t and 2 obs the residual variance (observation error): This observation model was based on the assumption that the counts were correct on average and adjusted for random variation in counting errors that can come either from double counting or missed individuals.
We first assessed site-specific population growth rates, with i,t reflecting temporal variability in the site-specific growth rate around the mean annual growth rate r i . We computed the probability that r i was larger than zero, that is the proportion of iterations from the posterior distribution of r i over zero.
We then applied different decompositions to the general model presented above to test the effect of explanatory variables on population growth rates, depending on the structure of the variables and the question Table 1. First, to test the effect of the surface area (S) of obs TA B L E 1 Explanatory variables used to investigate the potential effects of biotic and abiotic factors on median trends and/or on interannual variation in population growth rates of black grouse in the French Alps between 1979 and 2018 the site on initial population size and observation errors, we applied to the general model the following linear equation, with N i,fc (i) the initial population size at site i, i the intercept, and the slope of the effect: To estimate the impact of explanatory variables on site median trends, we used the following two equations, with r the common trend across sites, X i and X i,t the explanatory variable X at site i, in the first case time independent and in the second case recorded for year t: Finally, we tested the impact of explanatory variables on the detrended inter-annual variation in growth rates, that is the variations in mean annual growth rate around the median trend. When covariate values were comparable between sites, we applied the following model: Otherwise, we treated each site separately: We implemented models in the Bayesian context using the R packages "rjags" (Plummer et al., 2019) and "R2jags" (Su et al., 2015).
We had no prior knowledge about the growth of the populations, thus we specified vague priors for all parameters. We ran 3 chains of 100,000 iterations with 50,000 burn-in iterations and a thinning of 9. We checked that convergence was satisfactory by visualizing the mixing of chains and using the Brooks-Rubin-Gelman diagnostic (Brook & Gelman, 1998). To evaluate whether the explanatory variables had an effect on population growth, we computed the probability that the regression coefficients i were below zero as well as examining the size of the effect (Kéry & Schaub, 2012).

| 22-year trends of 47 sites (1996-2018)
The findings showed that the annual population growth rates obtained for the 47 sites were contrasting Figure  CI width mean = 0.151). Inter-annual variation in mean annual growth rates around the median trend was, on average, one-quarter lower than variation associated with observation error (sd proc median = 0.

| 39-year trends of 11 sites (1979-2018)
Conclusions from the 11 longer time series differed slightly from those obtained with the 47 sites monitored over 22 years Figure 3 & Appendix 2. The credibility intervals of mean annual growth rate estimations were, as expected, narrower with these longer time series, although they remained relatively wide (80% CI width mean = 0.070,95% CI width mean = 0.111). Once again, inter-annual variation in mean annual growth rates around the median trend were, on average, one-quarter lower than variation associated with observation error (sd proc median = 0.15, 95% CI 0.13-0.17,sd obs median = 0.19, 95% CI 0.18-0.31).
Considering an 80% CI threshold, no site showed any significant population increase, while three sites presented a significant decline (sites 2, 10 and 9, with median r ranging from −0.02 to −0.05, corresponding to a 54% and 86% decline over 39 years respectively).
Definitive conclusions could not be drawn for the remaining eight sites. Of these eight sites, three showed a median r between −0.02 and −0.01 (corresponding to a 54% and 32% decline over 39 years respectively), four sites between 0.009 and 0.016 (corresponding to a 40% and 83% increase) and one site close to 0. All three sites in the northern pre-Alps showed a median growth rate below −0.016 ( Figure 3 & Appendix 2).

| Explanatory variable effects on median trends (Figures 4 & 6), see all result graphs in Appendix 4
While the effect was not significant, our results suggest that median trends were negatively related to site surface area (

| Explanatory variables effects on the detrended inter-annual variation in growth rates (Figures 6, 7 & Appendix 5)
We found a strong positive linear relationship between inter-annual variation in growth rates and summer temperatures Figure 6d:

| D ISCUSS I ON
Our results showed that black grouse population trends were strongly heterogeneous, not just at a national scale, but also at a relatively small scale, with no significant general trends observed or even spatial autocorrelation of site-specific trends (for more details on spatial autocorrelation analysis refer to Appendix 3). Nevertheless, a general decline was found in the pre-Alps. While elevation and the start date of the vegetation season had a strong impact on site-specific median trends, no significant effect of snowmelt dates and length of vegetation season was found. Climatic covariates also strongly influenced inter-annual variation in growth rates. We found no significant effect of harvesting or cable density on site-specific median trends,however, there was a surprisingly positive correlation of inter-annual variation in growth rates with cable density. Finally, the impact of the number of harvested birds on inter-annual variation in growth rates differed between sites, with only a few sites showing strong negative or positive effects. Thus, during our study period, climatic conditions seemed to have had a much stronger influence than recreational activities on black grouse population growth rates in the French Alps.
F I G U R E 6 Impact of covariates on inter-annual variation in growth rates (mean ± 95%CI) estimated by the statespace models applied on the 47 time series of black grouse males counted at leks scattered throughout the French Alps. Panels show results from models that looked at different covariates: 100% snowmelt date (Date_100, from  The large CI width appears to be the consequence of both high inter-annual variation in population growth rates and high measurement errors (Auger-Méthé, 2016). This measurement error is probably linked to the behaviour of both observers and birds, the latter strongly influenced by meteorological conditions (Baumgardt, 2017).
Counting errors might be reduced either by more rigorous protocol standardization or by including variables associated to the observation process in the models. Unfortunately, we lacked any systematic data to be used as covariates that could help improve inference of the detection bias and model variance more accurately (Kéry et al., 2009). Moreover, as monitoring is partly conducted by volunteers, imposing more constraints on protocols might lead to a lower number of participants (Couvet et al., 2008;Dinckinson et al., 2010).
Our results showed that local population trends were strongly heterogeneous in space, even within the same region, with a balanced number of sites declining or increasing. While some sites lost up to 60% of their initial abundance over the 22 years study period, others remained stable and some even doubled their population size.
This result is surprising given the overall decline of black grouse in Europe (Storch, 2007). However, monitored sites we studied were not selected based on a probabilistic survey design but rather on subjective criteria, hence the overall trend is likely to be biased (Fournier et al., 2019;Yoccoz et al., 2001). Indeed, the setup of new sites occurred at various scales, based upon local needs. Historically, actors tended to select sites where access was convenient and black grouse populations large (or locally large, for example in the southern pre-Alps), often with good quality habitats. These sites may show different trends than the overall trend. Unfortunately, we have no precise quantitative data about setup choice,hence, no extrapolation is feasible with these time series to get national nor regional trends.
Besides, because lek counts are subject to numerous potential errors (Baines et al., 1996;Baumgardt et al., 2017), results have to be interpreted with caution. Nevertheless, all sites located at the southern edge of the distribution of the species also showed a clear decline in population size. The elevated temperature and shorter persistence of snow in these southern sites, probably partly explain this general decline in the black grouse, which is a glacial relict species adapted to cold temperatures (Bionda et al., 2017).
Other studies on different species of grouse have observed similar strong spatial heterogeneity of population growth rates even among close populations (Cattadori & Hudson, 1999;Martinoli, 2016). For example, a recent study on rock ptarmigan in the Swiss Alps found similar spatially heterogonous trends with no synchrony between populations (Furrer et al., 2016). Such strong spatial heterogeneity in population growth rates might be explained by the significant differences in local conditions even over short distances in mountains, which might overcome synchronizing drivers expected from large-scale climate patterns (Bjørnstad, 1999;Brambilla et al., 2019;Paradis et al., 2000).
Such a major role of local conditions seems to be confirmed by our study, which shows that, despite high heterogeneity in local trends, abiotic variables had a strong impact on both local trends and interannual variation in local population growth rates. On a given site, years with warm temperatures and early snowmelt and vegetation growth were associated with higher black grouse population growth rates. Numerous studies on grouse have demonstrated that vegetation and snowmelt phenology influence population growth rates (Clark & Johnson, 1992;Liebezeit et al., 2014;Novoa et al., 2016;Pulliainen & Tunkkari, 1991). Warm temperatures probably favour an early start to snowmelt and the vegetative growth season (Visser et al., 2009). An early vegetation season allows hens to benefit from protein-rich vegetation during egg production, along with higher resource availability during the chick-rearing period, factors that might improve chick body condition and in turn increase winter survival (Brittas, 1988;García-González et al., 2016). Regarding weather covariates, we also found a negative correlation between summer precipitation levels and mean annual population growth rate, which is consistent with previous studies of several species of grouse showing that high levels of precipitation and adverse climatic conditions during the hatching period and the two first weeks of a chick's life can negatively affect survival (Erikstad & Spidso, 1982), (Novoa et al., 2008;SAEther et al., 2004).
When looking at median population trends, we also found a significant quadratic effect of an early vegetation start. Consistent with the results on inter-annual variation described above, sites with late vegetation phenology were associated with lower population growth rates. However, sites that exhibited a very early phenology were also more likely to show decline. This discrepancy between different analyses is probably related to the fact that the start of the vegetation season is strongly negatively correlated with elevation (Appendix 7, rho = 0.79). Thus, sites with very early phenology were mostly those located in the southern pre-Alps, at the southern F I G U R E 7 Regression coefficients (mean ± 95%CI) of covariates' effects on long-term annual population growth rates, estimated by the state-space models applied on the 47 time series of black grouse males counted at leks scattered throughout the French Alps. Both α and β coefficients were represented for covariates for which the best model from selection was a quadratic relation (see full results in appendix 6) margin of the French distribution, where despite a long vegetative season, summer conditions might be too harsh, as suggested by the general decline observed at these sites (see above).
While our time series were probably too short to detect any strong signal of climatic covariates on long-term trends, the strong weather signals on inter-annual variation of population size highlighted in our study suggests that ongoing and future climate change should have a strong impact on population trends. Studies carried out in the French Alps using data from 1958 to 2005 showed that, temperature drops in late summer and early winter, together with slightly rising precipitation, resulted in earlier and stronger snowfalls, particularly at high elevation. Nevertheless, strongly rising temperatures in late winter and early summer together with decreasing precipitation caused less snowfall and thus a reduced snow depth and earlier melt, particularly at low and medium elevations. Besides, they showed that the number of days with snow on the ground has generally declined at low elevation, but showed little change at medium and high elevations (except for some central and southern massifs) from 1958 to (Durant et al., 2019. This study illustrates very well the complexity and high spatial and temporal climatic variability in mountainous areas. Nevertheless, large-scale climatic data have already shown long-term trends of temperature and precipitation, particularly marked in mountainous areas (Benniston, 2003;Scridel et al., 2018). Hence, a signal of climate impact on black grouse should progressively emerge as time series continue to lengthen, especially since climate scenario models predict an intensification of these changes (Dirnböck et al., 2011). The distribution of black grouse is expected to contract towards higher elevations due to ambient temperature warming (Braunisch et al., 2016;Chamberlain & Pearce-Higgins, 2013), our results on the impact of the start of the vegetation season as well as the observation of a general decrease in southern populations at lower latitude and altitude seem to confirm such a pattern. Upward displacements of leks have been reported in the Swiss Alps since 1970 following the altitudinal shift of the treeline induced by warming temperatures (Bionda et al., 2017;Zurell et al., 2012). Yet the strong impact of precipitation on population trends we found leads us to suspect that alterations in the distribution of black grouse populations will be difficult to predict as shifts in rainfall regimes are especially difficult to forecast in mountains (Tingley et al., 2012). Moreover, climate change will result in profound changes in both vegetation and peaks of insect communities, along with an increase in the frequency of extreme events (Benniston, 2003), making it difficult to predict accurately how the black grouse will respond to this interaction of direct and indirect factors.
Studies conducted on boreal black grouse showed that climate change can have complex consequences on the lifecycle of the species, with warmer springs resulting in an advanced onset of breeding. Yet summers were also colder and wetter, leading to a phenological mismatch between chick growth and the peak of food availability, and thus a negative relationship between an advanced spring and the population growth rate (Ludwig et al., 2006). But another study on Norwegian populations showed the opposite pattern, with population growth rates positively correlated with an advanced spring (Wegge & Rolstad, 2017). Finally, as for potential forest succession impact on populations, in opposition to what can be speculated in some other areas, such as Scotland or boreal forests, vegetation change in the Alps at the alpine limit is very slow, and because black grouse can occur in both meadows and young forests, the expected pattern is not so clear (Patthey et al., 2012;Ramazin et al., 2000). On most sites, it seems that habitat encroachment is not sufficient to be unfavourable to black grouse yet and could actually be beneficial on some sites (Braunisch et al., 2016;Scridel et al., 2017).
Concerning anthropogenic activities, contrary to our expectation, cable density, used as a proxy of recreational activity near the studied sites, was not correlated to population trends, but was positively correlated to variation in local growth rates (i.e. at specific sites local growth rate increased with cable density). These results are in contradiction with previous research carried out in Switzerland that showed a negative impact of ski lifts and related outdoor winter sports on black grouse abundance, after accounting for the effect of habitat type (Patthey et al., 2008). Our counterintuitive result may be explained by the fact that ski resort development in the French Alps started before World War II and continued into the 1980s (Jenni-Eiermann & Arlettaz, 2008), while our earliest dataset started in 1979 and most data is from the 1980s and 1990s. So the strong negative impact of ski resort development may have occurred before monitoring in our study area. Nonetheless, it is more complicated to explain the increase in local population growth rate with the increase in cable density on a specific site. Our study may suggest a certain tolerance, or even a local benefit, of human recreational activities for black grouse. A study carried out in the Bavarian Alps showed that some black grouse populations actually seemed to use quiet winter refuges near ski areas when favourable habitat is available and disturbance is predictable (Zeitler, 2007), such a pattern may also hold true in our area. Alternative hypothesis to explain the surprising positive correlation with cable density could be either an opening up of the forest structure that might benefit black grouse since they like mires and open patches in forests as display sites, or cascading effects of habitat change on predators' occurrence (Lyly et al., 2014).
Note that in our study, only ski-related cables were taken into consideration, it would be interesting to also explore a potential impact of powerlines on black grouse populations, as they are known to impact many bird species (Hovick et al., 2014).
In terms of hunting, our analysis did not find any significant longterm effects of harvest levels on population trends since 1996. The impact of site-specific harvest on inter-annual variation in growth rates during our study period differed between sites, but seemed to be equally distributed around zero. Current harvest rates are considered to be very low compared to what they were before the 1990s.
Today approximately 5% of the male population is estimated to be harvested each year; the number of harvested birds was about 10 times higher around 1986-88 (Magnani, 1995) and was 20 times higher in 1965 (Revue nationale de la chasse, 1966). Previous studies investigating the impact of such high hunting rates on local population dynamics showed a highly distorted adult sex ratio,however, this did not seem to lead to negative effects on breeding success, probably thanks to immigration from nearby populations (Ellison et al., 1988).
Today, hunting quotas are adjusted every year depending on the reproductive success of the previous year. This protocol may limit the negative impact of harvesting on population growth rates, and may even explain the positive relationship observed at some sites.
Nevertheless, our results do not demonstrate that the current level of hunting is sustainable at the national scale, as local harvesting may be compensated by immigration from unharvested areas (Little et al., 1993;Novaro et al., 2005). Further information regarding source-sink dynamics of black grouse populations in the study area would be needed to answer this question and determine is mortality is compensatory or additive (Lundberg & Jonzén, 1999).

| CON CLUS ION
Our study highlighted strong spatial heterogeneity, even at a relatively small scale, of black grouse population trends in the French Alps, while finding that populations at the margin of the distribution showed a strong decline. These results suggest that mountain bird species may be drastically affected by ongoing climate change, perhaps even more than by increased recreational activities. The black grouse, like other mountain birds, is expected to experience even more marked changes in abundance in response to profound changes in trophic networks through cascading effects linked to increasing climate change (Larson et al., 2019;Lyly et al., 2016;Ockendon et al., 2014;Pearce-Higgins et al., 2019). This conclusion, however, should be tempered by the fact that major components of bird population dynamics are habitat quality and its dynamics in space and time-information that was unfortunately not available in this study (Angelstam, 2004;Ludwig et al., 2009;Pearce-Higgins, 2006). While human recreative activities did not seem to be the major factor directly affecting the population dynamics of the black grouse in our study areas, its strong indirect effects on habitats, through forest dynamics and agricultural practices, compared to the slower predicted landscape changes due to climate change, were not taken into consideration in our study and merit further investigation for a better understanding (Ockendon et al., 2014).

ACK N OWLED G EM ENTS
We would like to thank all the participants involved in the monitoring of black grouse in the French Alps, which included volunteers and professional observers who are members of the Observatory of Mountain Galliformes (Observatoire des Galliformes de Montagne: OGM). We are grateful to the protected areas in which the monitoring took place (ASTERS Nature Reserve of Contamines Montjoie, Vanoise National Park, Ecrins National Park, Mercantour National Park, Bauges Regional Nature Reserve, Chartreuse Regional Nature Reserve, Vercors Regional Nature Reserve), to the French National Hunting and Wildlife Agency and the French National Forest Agency, and to the hunting organizations (FDC04, FDC05, FDC06, FDC26, FDC38, FDC73, FDC74) that provided hunting bag data. We would also like to thank A. Bernard-Laurent, D. Maillard, E. Menoni and C. Novoa, whose encouragement, useful advice and expertise in developing the study and throughout the project were greatly appreciated. M. Garel and C. Bernard provided technical assistance with MODIS and SAFRAN data, C. Calenge advised on the FPCA analysis, and N. Yoccoz, J.Y.
Barnagaud and J. Chiffard offered valuable help and advice.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available in DRYAD at https://doi.org/10.5061/dryad.w9ghx 3fnp. Note that because the studied species is considered susceptible in France, data have been randomized; hence, site numbers do not reflect original OGM site numbers.