Elie Gaget
email: elie.gaget@gmail.com

Thomas Galewski

Jon E Brommer

Isabelle Le Viol

Frederic Jiguet

Nicola Baccetti

Tom Langendoen

Blas Molina

Filipe Moniz

Caroline Moussy

Marco Zenatello

Matthieu Guillemain

Habitat management favouring hunted waterbird species prevents distribution changes in response to climate warming

Keywords: climate change, duck, Community Temperature Index, Mediterranean, persistence, wintering Oxyura leucocephala Plant Clangula hyemalis Invertebrate Other hunted 9 -2.499599758 Somateria mollissima Invertebrate Other hunted 9 -7.528618341

Climate warming is driving species to shift their geographic distribution poleward to track suitable climatic conditions. Two strategies have been suggested to help species respond to climate warming: facilitating distribution change or improving persistence. We questioned whether habitat management in favour of duck hunting activities interacted with duck responses to climate warming. We studied non-breeding waterbird community changes (110 species) over 28 years at 851 sites in South-West Europe, where habitat management is a common practice to attract and hunt ducks. We hypothesized that duck species targeted by habitat management do not need to track temperature changes as much as non-hunted species, because of the availability of suitable habitats provided by hunting land managers. We used the Community Temperature Index (CTI) to assess the temporal responses of communities and species to climate warming. We compared the effect of hunting status to other functional traits on species responses, controlling for phylogenetic relatedness. The CTI trend has increased over the study period, indicating a community adjustment to climate warming. However, hunted ducks contributed to almost 40% of the negative contributions to this community adjustment, suggesting that hunted ducks do not shift their distribution as much as the other waterbirds do. Winter fidelity associated with the provision of attractive feeding grounds might explain why ducks did not seem to shift their distribution in response to climate warming. This study suggests the broad impact of human activities on wildlife, including on large-scale distribution processes, and questions the longterm consequences on duck populations.

Introduction

Global warming is a major threat to biodiversity, leading to massive species distribution changes [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF]. Global temperatures have increased by 1.1°C since the pre-industrial period (IPCC 2021), and are expected to reach +2°C to +5°C by 2100 [START_REF] Raftery | Less than 2 C warming by 2100 unlikely[END_REF]. A possible step backwards to pre-industrial temperatures in the next centuries is uncertain (IPCC 2021). Therefore, temperature changes are not only rapid, but also likely to last over time.

Species responses to climate warming are diverse but their intensity is largely insufficient [START_REF] Radchuk | Adaptive responses of animals to climate change are most likely insufficient[END_REF][START_REF] Lenoir | Species better track climate warming in the oceans than on land[END_REF]. Species may respond to climate warming through changes in migration distance, geographic shift, phenology, plasticity, or microevolution [START_REF] Jiguet | Climate envelope, life history traits and the resilience of birds facing global change[END_REF], Moller et al. 2010[START_REF] Dawson | Beyond predictions: biodiversity conservation in a changing climate[END_REF][START_REF] Elmberg | Interpreting seasonal range shifts in migratory birds: a critical assessment of 'short-stopping' and a suggested terminology[END_REF]. The most commonly observed and documented response is geographical distribution change, leading to changes in community composition in response to temperature changes at the local scale [START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF]. For instance, temperature increases might cause local population increases at the cold edge of species distributions, and decreases in population sizes at the warm distribution edge [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF][START_REF] Thomas | Range retractions and extinction in the face of climate warming[END_REF], Pavón-Jordan et al. 2015). Species can also persist locally, either by finding micro-refuges (e.g. at higher altitude, [START_REF] Massimino | Can microclimate offer refuge to an upland bird species under climate change?[END_REF] or by tolerating changes through phenotypic or genetic adaptation (plasticity, micro-evolution) [START_REF] Charmantier | Adaptive phenotypic plasticity in response to climate change in a wild bird population[END_REF][START_REF] Karell | Climate change drives microevolution in a wild bird[END_REF]. Functional traits play an important role to shape species responses to climate warming [START_REF] Brommer | Extent of recent polewards range margin shifts in Finnish birds depends on their body mass and feeding ecology[END_REF]. Species sharing the same characteristics, like climatic niche, body mass, migration strategies, might be more likely to respond in a similar way [START_REF] Jiguet | Climate envelope, life history traits and the resilience of birds facing global change[END_REF]) while facing the same vulnerabilities to temperature increase [START_REF] Aubin | Traits to stay, traits to move: a review of functional traits to assess sensitivity and adaptive capacity of temperate and boreal trees to climate change[END_REF]. Species responses to climate warming may also be modulated by local environmental factors, like habitat quality [START_REF] Dawson | Beyond predictions: biodiversity conservation in a changing climate[END_REF][START_REF] Lenoir | Species better track climate warming in the oceans than on land[END_REF]) and conservation efforts [START_REF] Gaget | Protected area characteristics that help waterbirds respond to climate warming[END_REF].

Helping species respond to climate warming is mainly considered under two complementary climate adaptation strategies: facilitation and persistence [START_REF] Bonebrake | Managing consequences of climate-driven species redistribution requires integration of ecology, conservation and social science[END_REF][START_REF] Prober | Shifting the conservation paradigm: a synthesis of options for renovating nature under climate change[END_REF][START_REF] Pearce-Higgins | A framework for climate change adaptation indicators for the natural environment[END_REF]. Facilitation focuses on helping species shifting their geographical distribution and settle in new areas in response to global warming (e.g. [START_REF] Maclean | Guidelines on the measures needed to help waterbirds adapt to climate change[END_REF][START_REF] Guillemain | Ready for climate change? Geographic trends in the protection status of critical sites for Western Palearctic ducks[END_REF]. It aims at adjusting the rate of species distribution change to the thermal isocline shift. For example, increasing connectivity between protected habitats facilitates colonization processes in response to climate warming [START_REF] Mcguire | Achieving climate connectivity in a fragmented landscape[END_REF]. Conversely, persistence aims at maintaining species in an environment despite temperature changes. It may indirectly be based on the intrinsic abilities of species to tolerate the changes (e.g. plasticity) or involve management measures to provide additional resources or create microrefugia benefiting a targeted species. For example, food supply may decrease the local extinction risk of species vulnerable to climate warming (Correira et al. 2015). Both facilitation and persistence strategies are likely to influence the response of species to climate warming [START_REF] Bonebrake | Managing consequences of climate-driven species redistribution requires integration of ecology, conservation and social science[END_REF][START_REF] Prober | Shifting the conservation paradigm: a synthesis of options for renovating nature under climate change[END_REF].

Empirical evidence of long-term species response to climate change has rarely been demonstrated under the persistence strategy (Correira et al. 2015[START_REF] Pearce-Higgins | Site-based adaptation reduces the negative effects of weather upon a southern range margin Welsh black grouse Tetrao tetrix population that is vulnerable to climate change[END_REF]. Indeed, as anthropogenic climate change is a relatively recent biodiversity pressure, conservation measures have not primarily aimed at improving long-term species persistence in response to temperature increase specifically. In contrast, efforts have historically aimed at improving species conservation, with consequences on their distribution regarding temperature change. For example, the creation of protected areas to preserve migratory waterbirds during cold spells in Western Europe [START_REF] Ogilvie | A migration study of the Teal (Anas crecca) in Europe using ring recoveries[END_REF][START_REF] Guillemain | The importance of protected areas as nocturnal feeding grounds for dabbling ducks wintering in Western France[END_REF] are also buffering the negative impact of temperature increase on some colddwelling waterbirds [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF]. Similarly, the provision of supplementary foods in British gardens changed the migration behaviour of the Eurasian blackcap (Sylvia atricapilla), of which many individuals are no longer migrating South [START_REF] Plummer | Is supplementary feeding in gardens a driver of evolutionary change in a migratory bird species?[END_REF]. Such anthropogenic influences can be examined a posteriori in the light of climate warming [START_REF] Pearce-Higgins | A framework for climate change adaptation indicators for the natural environment[END_REF].

Here, we question the effect of habitat management for hunting purposes on waterfowl responses to climate warming at international scale in South-West Europe. Duck hunting is a popular practice in this region yielding 2.4 million harvested ducks annually (EEA 2020). Indeed, habitat management for hunting purposes is a local, but widespread practice (FACE 2019), aiming at attracting targeted waterfowl species to hunting grounds during the non-breeding period. Pro-active habitat management actions for hunting purposes include the creation of man-made ponds, spreading of seeds as bait, artificial flooding, vegetation clearing including the creation of openings in reedbeds, and other actions (e.g. [START_REF] Mathevet | Que ferons-nous des canards sauvages ? Chasse, nature et gestion adaptative[END_REF][START_REF] Vallecillo | Expériences cynégétiques innovantes en Camargue : des pistes pour la chasse au gibier d'eau de demain ?[END_REF]. Economic incentive for such intensive management actions is that sites with high abundance of ducks can yield high income by renting out hunting rights (e.g. [START_REF] Niang | A cost-benefit analysis of rice field winter flooding for conservation purposes in Camargue, Southern France[END_REF]. Specific habitat management targeted to boost abundance of hunted species is often carried out in large private properties, e.g. in the Camargue, Southern France [START_REF] Tamisier | A review of habitat changes in the Camargue: an assessment of the effects of the loss of biological diversity on the wintering waterfowl community[END_REF] and provide prestigious social status and reputation of the land-owner, e.g. in Italy. However, if such management measures attract hunted duck species in large numbers, they may not be so attractive to other waterbird species, unless for species sharing common habitat or diet characteristics with hunted ducks, such as granivorous and herbivorous species (Tamisier andDehorter 1999, Mouronval et al. 2014). In particular, such habitat management are not expected to provide benefits to invertebrate or fish eating species (e.g. shorebirds, herons) [START_REF] Tourenq | Longterm population and colony patterns of four species of tree-nesting herons in the Camargue, South France[END_REF]. [START_REF] Gaget | Waterbird communities adjust to climate warming according to conservation policy and species protection status[END_REF] showed that in the south-western European Union Member States, waterbird responses to climate warming were correlated with their protection status. Hunted species had a lower response to climate warming through distribution change than the other species. Here, we further explore the reason for this interaction between hunting status and climate warming, by investigating the possible role of habitat management practices on hunted ducks.

We studied community reshufflings from 110 non-breeding waterbird species in south-western Europe over 28 years on 851 sites. Many of these waterbirds breed at northern latitudes and overwinter in wetlands in this region. Our main interest was to explore whether the distributional changes of duck species, which are actively targeted by hunting in this region, differ from those of other hunted and non-hunted waterbird species, while controlling for functional traits that differ between these species. To infer distribution changes we used the Community Temperature Index (CTI) framework, a community weighted mean index used to measure the average temperature affinity of all the species of a community, based on individual species temperature index (STI) [START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF]. We assessed temporal changes of CTI over years for the non-breeding waterbird community using abundance time series collected at site level through the International Waterbird Census [START_REF] Delany | Guidance on Waterbird Monitoring Methodology: Field Protocol for Waterbird Counting (Wetlands International)[END_REF] and STI calculated as the mean winter temperature values across the non-breeding distribution of the species [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF]. Under a climate warming scenario, we expected the avian community at a given site to change over time towards species with a higher thermal affinity. Thus, the CTI was expected to increase over time. Such community reshuffling has already been documented for nonbreeding waterbirds in the western Palearctic [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF] and in numerous other studies [START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF][START_REF] Auffret | Synergistic and antagonistic effects of land use and nonnative species on community responses to climate change[END_REF][START_REF] Fourcade | Climate and land-cover change alter bumblebee species richness and community composition in subalpine areas[END_REF]. We then assessed the species' contributions to this temporal increase in CTI. Our main interest was whether this contribution depends on the species' hunting status (hunted ducks, other hunted species and non-hunted species), while accounting for functional traits (species thermal affinity, body mass, generation length, diet and migration strategy) and phylogeny [START_REF] Jiguet | Climate envelope, life history traits and the resilience of birds facing global change[END_REF][START_REF] Brommer | Extent of recent polewards range margin shifts in Finnish birds depends on their body mass and feeding ecology[END_REF] to control for potential confounding effects. Because of the intensive management practices in the study region [START_REF] Tamisier | A review of habitat changes in the Camargue: an assessment of the effects of the loss of biological diversity on the wintering waterfowl community[END_REF], Pernollet et al. 2015a[START_REF] Zenatello | International Waterbird Census Report -Italy 2009-2018[END_REF], we expected that hunted ducks would be slower to adjust their distribution to climate warming compared to non-hunted waterbirds, because management of hunting areas may provoke a greater site fidelity. We thus expected that hunting status shapes a species' response to climate warming.

Methods

We studied non-breeding waterbird abundance time series collected from 1993 to 2020 at 851 sites in Portugal, Spain, France and Italy. The data were collected for the International Waterbird Census coordinated by Wetlands International [START_REF] Delany | Guidance on Waterbird Monitoring Methodology: Field Protocol for Waterbird Counting (Wetlands International)[END_REF], resulting in one survey per year performed by skilled ornithologists in mid-January on hundreds of wetlands. We filtered the original dataset (n = 3,288 sites) to ensure the quality of the data. To do so, we restricted the list of waterbird species to the list of the Agreement on the Conservation of African-Eurasian Migratory Waterbirds (www.unepaewa.org), we removed count events with less than two species to allow estimation of species contributions to CTI trend (see below) and we removed species with less than 100 individuals observed among all sites included in the study. Because of the particular focus on ducks, we removed sites without any hunted duck observation. We retained sites with at least five surveys and one survey per decade (1990s, 2000s and 2010s), resulting in a subset of 851 sites with 20.5±6.2 surveys on average (±SD) (see map Appendix 1). Considering taxonomic changes that occurred during the study period, Larus michahellis, L. cachinnans and L. argentatus were lumped into one 'species' L. michahellis. This lumping has a negligible impact on most of the study area where L. michahellis is largely dominant over the two other species, with the exception of the French Atlantic coast. The dataset used for the analyses resulted in 17,462 surveys of 851 sites over 28 years, with a cumulative record of 119 million birds from 110 species (Appendix 2).

We considered three hunting status categories between species: eleven duck species mainly targeted by habitat management, 25 other hunted waterbirds and 75 non-hunted species (Appendix 2). Hunted ducks are Northern Pintail (Anas acuta), Common Teal (A. crecca), Mallard (A. platyrhynchos), Eurasian Wigeon (Mareca penelope), Gadwall (M. strepera), Northern Shoveler (Spatula clypeata), Garganey (S. querquedula), Common Pochard (Aythya ferina), Tufted Duck (A. fuligula), Greater Scaup (A. marila) and Red-crested Pochard (Netta rufina). These species are all listed in Annex II of the Birds Directive, meaning that they can be hunted according to the regulation of the European Union or Member States (Red-crested Pochard is hunted in France and Spain only, and Greater Scaup is hunted in France only). The other hunted species are listed in Annex II of the Birds Directive and hunted at least in one Member State of the study area, and include species such as Greylag Goose (Anser anser), Coot (Fulica atra), Lapwing (Vanellus vanellus) and Eurasian Curlew (Numenius arquata). The non-hunted species are not listed in Annex II, or listed in Annex II but not hunted in the study area, including species such as Mute Swan (Cygnus olor), Ferruginous Duck (Aythya nyroca), Grey Heron (Ardea cinerea), Eurasian Spoonbill (Platalea leucorodia), Little Grebe (Tachybaptus ruficollis), Great Cormorant (Phalacrocorax carbo) and Mediterranean Gull (Larus melanocephalus).

We document species functional traits and phylogenetic proximities to control for potential confounding effects (Appendix 2). Functional traits were species thermal affinity, body mass, generation length, diet and migration strategy. These functional traits are expected to influence bird response to climate warming [START_REF] Jiguet | Climate envelope, life history traits and the resilience of birds facing global change[END_REF][START_REF] Brommer | Extent of recent polewards range margin shifts in Finnish birds depends on their body mass and feeding ecology[END_REF]. Species thermal affinity was defined as the Species Temperature Index (STI), calculated as the average temperature in January (1950-2000, www.worldclim.org, spatial resolution of 0.25°) on the non-breeding species distribution (BirdLife International & HBW 2018) from [START_REF] Gaget | Waterbird communities adjust to climate warming according to conservation policy and species protection status[END_REF]. Body mass (continuous, log transformation) and diet (Plant, Omnivore, Invertebrate, Vertebrate) were extracted from [START_REF] Wilman | EltonTraits 1.0: Species-level foraging attributes of the world's birds and mammals: Ecological Archives E095-178[END_REF] (where omnivore corresponds to a combination of less than 50% of all other categories), generation length (continuous) from IUCN (2021) and migration strategy (resident, short, long) from [START_REF] Storchová | Life-history characteristics of European birds[END_REF]. Migration strategies were represented in similar proportions between hunting status categories, with most of hunted ducks being short to long-distance migrants [START_REF] Storchová | Life-history characteristics of European birds[END_REF]; Appendix 2). To control for species phylogeny, we used phylogeny trees from BirdTree (https://birdtree.org).

We used temperature anomalies to assess changes in winter temperatures on the survey sites over the study period. We collected average monthly temperatures gridded at 0.5° resolution, from the Climatic Research Unit v4.05 [START_REF] Harris | Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset[END_REF]. Winter temperatures anomalies were computed each winter at site level as the average of the mean monthly temperatures for the 3-month period of November, December and January, i.e. the duration expected to impact waterbird distribution at winter time [START_REF] Gaget | Antagonistic effect of natural habitat conversion on community adjustment to climate warming in nonbreeding waterbirds[END_REF], minus the average monthly temperatures from November to January for the 1993-2020 period.

Statistical analyses

To assess and explain patterns of species responses to climate warming, we first calculated the longterm trend of CTI over years, then we assessed species contribution to the CTI trend, and finally we investigated the effect of hunting status and functional traits on species contribution, taking species phylogeny into account.

The CTI is a community weighted mean index used to measure the average temperature affinity of all the species included in the community based on individual STI [START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF]). The CTI is the average STI weighted by the loge(abundance + 1) of the species observed during a survey [START_REF] Gaget | Waterbird communities adjust to climate warming according to conservation policy and species protection status[END_REF]). We calculated one CTI value per count event, based on all the species observed at this time among the 110 waterbird species included in the study. Over time, a CTI can increase in response to climate warming by a relative increase of species with high STI, or a relative decrease of species with low STI.

We assessed CTI temporal trends with a linear mixed effect model (LMM) with the CTI as the response variable and year as a fixed effect. Sites, nested in countries, were used as a random effect to account for heterogeneity among sites and national monitoring schemes [START_REF] Gaget | Waterbird communities adjust to climate warming according to conservation policy and species protection status[END_REF]. Spatial autocorrelation was taken into account with an exponential structure on the geographic coordinates [START_REF] Gaget | Waterbird communities adjust to climate warming according to conservation policy and species protection status[END_REF]. We used the same model with temperature anomalies as response variable, to assess temporal changes in temperature anomalies. To visualize inter-annual changes in CTI and temperature anomalies, we also used the same models but with year as a categorical variable.

To assess the contribution of each species to CTI trends, we ran a jackknife analysis by removing each species in turn from the CTI before running LMMs on partial CTIs. Hence, species contribution to the CTI temporal trends was assessed as the difference between the original temporal trend (all species included) and the temporal trend of the model based on the partial CTI (all but one species). Species contributions were then reported as a percentage of influence, positive or negative, on the CTI trends [START_REF] Gaget | Waterbird communities adjust to climate warming according to conservation policy and species protection status[END_REF]. Note that species contributions are relative values, comparing species to the whole community. A negative contribution to a positive CTI trend means that a species is contributing relatively less than the other species to the community adjustment to climate warming.

Finally, we used phylogenetic generalized least squares (PGLS) models to assess differences between species contribution to CTI changes. We compared species contributions between hunting statuses and functional traits (species thermal affinity, body mass, generation length, diet and migration strategy), considering species phylogeny. We used a forward model selection approach to select the most parsimonious model with predictors explaining the most species contributions, i.e. selecting models with the lowest AIC among a delta AIC <2 [START_REF] Burnham | Model Selection and Multimodel Inference: a Practical Information-theoretic Approach[END_REF]. Phylogeny is a potential factor structuring adaptive bird response to climate change [START_REF] Hällfors | Shifts in timing and duration of breeding for 73 boreal bird species over four decades[END_REF]. Hunted ducks are from the same Anatidae family and often share the same genus. Their intrinsic sensitivity to climate warming and adaptive capacities are likely closer among each other than with species from other taxonomic clades. Consequently, contrasting species groups exposes the analytical framework to a potential confounding effect. To avoid this bias, we accounted for phylogenetic relatedness by adding the phylogeny tree as a Brownian correlation structure in the models. In addition, we assessed the multicollinearity between functional traits with the Variance Inflation Factor (VIF).

We performed two sensitivity analyses (Appendix 3) to ensure that the results were independent of 1) the site protection status, because protected areas can also interact with wintering waterbird responses to climate warming [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF], and of 2) the number of surveys per site, by selecting only sites with more than 25 surveys over the 28 years period.

Statistical analyses were performed in R 4.1.2 (R Development Core Team 2021). LMM and PGLS were run with package 'nlme' [START_REF] Pinheiro | Package 'nlme'. Linear and nonlinear mixed effects models[END_REF], Nagelkerke pseudo R-squared were assessed with package 'rcompanion' [START_REF] Mangiafico | Summary and Analysis of Extension Program Evaluation in R, version 1[END_REF] and VIF with package 'performance' [START_REF] Lüdecke | Package 'performance[END_REF].

Results

The Community Temperature Index (CTI) showed a very clear linear increase over years (β = 0.045 ± 0.001 SE, p < 0.001), indicating that the avian community was, on average, changing according to species thermal affinity (Fig. 1). Winter temperatures increased by 0.6°C on average over the study period on the studied sites (β = 0.023 ± 0.001 SE, p < 0.001).

However, not all species contributed equally to this change in CTI. Species diet and hunting status were the main factors explaining species contribution (βPlant = -2.61 ± 1.19 (±SE) [-4.74;-0.48] (95%CI), βInvertebrate = -0.15 ± 1.19 [-0.39;0.09], βVertebrate = -1.27 ± 1.20 [-1.45;-1.09], βOmnivore = 1. 06 ± 1.19 [0.86;1.26]; βHunted ducks = -3.62 ± 1.24 [-5.09;-2.15], βOther hunted = 0.57 ± 1.19 [-2.20;1.57], βNon-hunted = 0.82 ± 1.19 [0.15;0.89]), whereas other predictors had little effect (Table 1, Appendix 4). Multicollinearity between species diet and hunting status was low (VIF = 1.05). Most of the variance explained by the model was related to species diet (82%), with plant eating species contributing predominantly negatively to the CTI change (Fig. 2). Indeed, 83% of the plant-eating species had negative contributions, compared to 48%, 44% and 64% for species eating invertebrates, vertebrates and being omnivores, respectively (omnivores potentially comprising some plant-eating diet). Hunting status also strongly affected species contributions to CTI change (Table 1, Fig. 2). Most (73%) of hunted duck species contributions were negative (summing to -39.8%). Only Mallard, Greater Scaup and Garganey, among 11 hunted duck species, made positive contributions to CTI change (summing to +15.5%), although the Mallard accounted for almost all positive contribution (+14.9%; Fig. 2). Among the other hunted and non-hunted species, 63% and 52% of the species contributed negatively, respectively.

The sensitivity analyses conducted 1) without sites included inside protected areas and 2) only on site with more than 25 surveys, returned fairly similar results, indicating the importance of the species hunting status to explain species contributions to CTI change (Appendix 3).

Discussion

Over the study period , the observed CTI increase demonstrates that non-breeding waterbird communities have changed, with warm-dwelling species becoming relatively more abundant than cold-dwelling species in the region. This community change was expected because of the temperature increase in the study area, suggesting that non-breeding waterbird communities are reshuffling in response to climate warming. This finding is fully consistent with other studies on nonbreeding waterbirds [START_REF] Godet | Waders in winter: long-term changes of migratory bird assemblages facing climate change[END_REF][START_REF] Pavón-Jordán | Effects of flyway-wide weather conditions and breeding habitat on the breeding abundance of migratory boreal waterbirds[END_REF][START_REF] Gaget | Waterbird communities adjust to climate warming according to conservation policy and species protection status[END_REF][START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF], breeding birds [START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF]) and other taxa [START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF][START_REF] Auffret | Synergistic and antagonistic effects of land use and nonnative species on community responses to climate change[END_REF][START_REF] Fourcade | Climate and land-cover change alter bumblebee species richness and community composition in subalpine areas[END_REF] in Europe. This reshuffling is a long-term process caused by species distribution changes in response to climate warming, with warm-dwelling species colonizing new sites at the leading distribution edge and cold-dwelling getting extirpated at the trailing distribution edge [START_REF] Thomas | Range retractions and extinction in the face of climate warming[END_REF][START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF]).

While the community on average has responded to climate change, individual species responses varied. Some species distribution strongly tracked climate warming, and such species had, in the CTI framework used in this study, positive contributions to the CTI trend. Examples include Common Sandpiper (Actitis hypoleucos), Goldeneye (Bucephala clangula) and White stork (Ciconia ciconia). However, other species barely adjusted their distribution to climate warming and these species had negative contributions to the CTI trend. We here showed that most ducks targeted by habitat management for hunting purposes (except Mallard, and to a lesser extent Garganey and Greater Scaup) had negative contributions to the CTI trend. This finding therefore suggested that habitat management may influence species response to climate warming. We believe there are two nonmutually exclusive explanations for this finding.

Firstly, managing sites to benefit hunted species may have lowered the propensity of these species to adjust their distribution to climate warming. Indeed, intensification of habitat management for hunting purposes has been suggested as one of the reasons for contemporary changes in duck diet, body condition and migration phenology (Mathevet & Tamisier 2002;[START_REF] Guillemain | Disentangling the drivers of changes in Common Teal migration phenology over 50 years: land use vs. climate change effects[END_REF]. Various habitat management actions are commonly used to attract wintering ducks (e.g. [START_REF] Mathevet | Que ferons-nous des canards sauvages ? Chasse, nature et gestion adaptative[END_REF][START_REF] Vallecillo | Expériences cynégétiques innovantes en Camargue : des pistes pour la chasse au gibier d'eau de demain ?[END_REF]. One potentially important management action is food supplementation, where seeds (wheat, corn or rice) are provided as bait, a common practice in South-West European countries [START_REF] Brochet | Diet and feeding habitats of Camargue dabbling ducks: what has changed since the 1960s?[END_REF]. For example, such bait may account for half of the Teal food intake in some areas [START_REF] Legagneux | Seeds consumed by waterfowl in winter: a review of methods and a new web-based photographic atlas for seed identification[END_REF], and feeding activity for this duck species has shifted these past decades from natural habitats to artificially managed freshwater marshes [START_REF] Brochet | Diet and feeding habitats of Camargue dabbling ducks: what has changed since the 1960s?[END_REF].

Over the same period, Teal body mass has increased by more than 10% at their major Camargue wintering site in southern France [START_REF] Guillemain | Wintering French Mallard and Teal are heavier and in better body condition than 30 years ago: effects of a changing environment?[END_REF], and arrival date at this wintering ground has advanced [START_REF] Guillemain | Disentangling the drivers of changes in Common Teal migration phenology over 50 years: land use vs. climate change effects[END_REF]. Intensive water level management is another management practice likely influencing duck distribution change or local persistence in response to climate warning. Wetlands are flooded in summer, including temporary ponds that otherwise would be dry. The aim of this summer flooding is to maximize plant biomass growth in these ponds to coincide with the arrival of ducks from their breeding grounds in autumn [START_REF] Tamisier | A review of habitat changes in the Camargue: an assessment of the effects of the loss of biological diversity on the wintering waterfowl community[END_REF][START_REF] Davis | Habitat and resource use by waterfowl in the northern hemisphere in autumn and winter[END_REF]). Supplementation of freshwater in brackish habitats (e.g., coastal lagoons), aims also to increase habitat diversity hence attracting freshwater duck species. All such management actions improve habitat quality and attractiveness to ducks, resulting in large alterations of duck distribution in the last decades, which now concentrate in areas of intensive habitat management for hunting purposes [START_REF] Zenatello | Risultati dei censimenti degli uccelli acquatici svernanti in Italia. Distribuzione, stima e trend delle popolazioni nel 2001-2010[END_REF][START_REF] Zenatello | International Waterbird Census Report -Italy 2009-2018[END_REF]. Habitat management for hunting purposes may hence have counterbalanced the effect of climate warming and "retain" hunting ducks to these wintering grounds, which they would otherwise abandon for more northern areas.

Secondly, hunted ducks may have been slower to adjust their distribution to climate change because agricultural practices were particularly beneficial for these species. In particular, rice fields are flooded and used as foraging grounds by ducks during the non-breeding season, with ducks feeding on nonharvested rice and weed seeds (Pernollet et al. 2015a). Other waterbird species use rice fields for feeding to a lesser extent, except for species than can specialize on rice, such as the Greater Flamingo (Phoenicopterus roseus) during the seeding period in spring (Pernollet et al. 2015a). In Portugal, Spain, France and Italy, some wetlands hosting a large number of ducks are close to large rice producing areas, e.g. Ebro delta, Albufera, Camargue, Pavia. In such regions, the number of wintering ducks was positively related to the area of flooded harvested ricefields (Pernollet et al. 2015b). Furthermore, rice farming practices and management of sites to promote hunting are not mutually exclusive, since harvested ricefields can also be flooded for hunting purposes [START_REF] Niang | A cost-benefit analysis of rice field winter flooding for conservation purposes in Camargue, Southern France[END_REF][START_REF] Pernollet | L'utilisation des rizières par les canards hivernants : vers une gestion des rizières en interculture favorable aux canards et aux agriculteurs[END_REF]. Indeed, in addition to agronomic benefits, this practice provides suitable wintering grounds for ducks and may therefore promote site fidelity.

Despite these explanations, our current analysis cannot definitely infer the reasons for the slower distribution adjustment of hunted ducks to climate warming compared to other waterbirds. In particular, we conducted a comparison between species and not between sites. To establish that hunting-oriented site management promotes persistence of hunted ducks, site-level information on management practices and intensity in South-West Europe is needed. Spatialized data on habitat management for hunting purposes are currently lacking at regional or national scale. Importantly, ducks are often mobile over large areas within their winter range [START_REF] Caizergues | Emigration rates and population turnover of Teal (Anas crecca) in two major wetlands of Western Europe[END_REF][START_REF] Gourlay-Larour | Movements of wintering diving ducks: new insights from nasal saddled individuals[END_REF]) and may hence be affected by large scale management driven by regional factors in addition to site-specific management, such as land-use policy or socio-cultural practices. The fact that ducks commute daily between a day-roost (generally protected sites where they can rest and where they are surveyed) and surrounding nocturnal feeding areas (functional units sensu [START_REF] Tamisier | L'existence d'unités fonctionnelles démontrée chez les sarcelles d'hiver en Camargue par la biotélémétrie[END_REF]) further inflates uncertainty in how these species are affected by habitat management in a given area [START_REF] Vallecillo | Modélisation de la répartition spatio-temporelle des oiseaux d'eau en Camargue -un outil pour la gestion[END_REF]). An experimental design contrasting the effect of site management on duck carrying capacity would be challenging because of the mismatch between diurnal bird monitoring and nocturnal feeding behaviour [START_REF] Brochet | Marsh management, reserve creation, hunting periods and carrying capacity for wintering ducks and coots[END_REF].

Responses to climate warming differ between species for many reasons, including, among others, species movement ability or adaptive capacities, [START_REF] Dawson | Beyond predictions: biodiversity conservation in a changing climate[END_REF]). Our best model only explained 5.8% of the differences between species contributions, suggesting that an important part of the variation was not captured by the tested functional traits. This might have been expected because waterbirds share more ecological characteristics than we could get by analyzing a more phylogenetically diverse bird community (e.g. [START_REF] Jiguet | Climate envelope, life history traits and the resilience of birds facing global change[END_REF][START_REF] Brommer | Extent of recent polewards range margin shifts in Finnish birds depends on their body mass and feeding ecology[END_REF]. Consequently, accurate traits might be neglected, e.g. at population-level across species flyways, but needed to reflect local adaptations to environmental parameters, including climate warming (DeMarche et al.

2019).

Waterbird diet was an important predictor of species contribution to CTI trend, but other functional traits were not. Plant eating species mainly contributed negatively to the community changes. One likely reason might be that species sharing similar feeding ecology to hunted ducks, mostly plant eaters, were also affected by habitat management, like coots and swans [START_REF] Brochet | Marsh management, reserve creation, hunting periods and carrying capacity for wintering ducks and coots[END_REF]). Thus, habitat management for hunting purposes may also benefit plant eaters in general, to some extent. Generation length actually explained a fair amount of variation, but was not selected according to parsimony. It may therefore be an additional factor explaining why some species have a delayed distribution change [START_REF] Maclean | Guidelines on the measures needed to help waterbirds adapt to climate change[END_REF]. Indeed, juvenile settlement can influence over-winter location when mature, meaning that long-lived species might need more time to shift their distribution than short-lived species [START_REF] Godet | Waders in winter: long-term changes of migratory bird assemblages facing climate change[END_REF].

Mallards did not follow the same pattern as the other hunted ducks. We found that their distribution was tracking climate warming. Indeed, this species has started to shorten its migration distance decades ago along its European flyway, likely in response to temperature increases [START_REF] Sauter | Evidence of climate change effects on withinwinter movements of European Mallards Anas platyrhynchos[END_REF][START_REF] Gunnarsson | Direct and indirect effects of winter harshness on the survival of Mallards Anas platyrhynchos in northwest Europe[END_REF], Guillemain et al. 2015). Previous work highlighted that the Mallard clearly responds to winter temperatures, both in North America [START_REF] Schummer | Weather-related indices of autumn-winter dabbling duck abundance in middle North America[END_REF]) and in Europe [START_REF] Dalby | Temperature does not dictate the wintering distributions of E uropean dabbling duck species[END_REF]). Its metabolism and opportunistic feeding strategy may explain these findings [START_REF] Schummer | Weather-related indices of autumn-winter dabbling duck abundance in middle North America[END_REF]): its larger size and weight allow Mallard to endure negative temperatures more easily than smaller species, disadvantaged by their greater surface area to volume ratio. It is also a relative habitat generalist species exploiting numerous wetland habitats, small, even degraded or urban and can switch from water to terrestrial feeding while the other duck species are more wetland dependent. Accordingly, Mallard are likely to winter closer to their temperature limits than other duck species, resulting in a different trade-off between saving energy from migration and thermal regulation, and likely a quicker response to climate warming by a distribution change. It cannot be excluded, however, that the fact that millions of mallards are released annually in Europe [START_REF] Champagnon | Assessing the genetic impact of massive restocking on wild mallard[END_REF]) also blurs the relationship between environmental conditions and the geographic distribution of this species [START_REF] Söderquist | Longevity and migration distance differ between wild and hand-reared mallards Anas platyrhynchos in Northern Europe[END_REF]. Garganey and Greater Scaup were the two other duck species with a positive contribution to the CTI trend. This result should however be considered with great caution, given both species are rare winter visitors in the study area, the former usually wintering in Africa and the latter at more northern latitudes. Indeed, Garganey and Greater Scaup numbers had a cumulative abundance of 983 and 37,634 counted individuals in our dataset, respectively, compared to 15.3 million Mallards.

Perspectives

The economic incentive provided by hunting is expected to drive landowners to manage habitat to locally boost hunted species abundance. This kind of habitat improvements could locally release hunted species from some of the adverse effects of climate change, but also impede current or future distributional changes. Consistent with this prediction, and based on community changes in 110 waterbird species over a large region (south-west Europe), we demonstrated here that hunted ducks are indeed slower in their response to climate warming compared to other waterbirds. Although our analysis did not directly infer that site management for hunting purposes drives this outcome, we suggest the pattern we found may be general. Because hunting is a global phenomenon, more studies exploring the interactions between site management and climate change as drivers of altered distributional changes are needed to understand whether hunted ducks are more persistent against climate change because of habitat management, or trapped there despite the ongoing climate warming. While hunted ducks may in the short-term benefit from such artificial land management, this may already have negative consequences for other elements of biodiversity [START_REF] Tamisier | A review of habitat changes in the Camargue: an assessment of the effects of the loss of biological diversity on the wintering waterfowl community[END_REF], and may be detrimental for the ducks themselves in the longer run.

In addition to the potential effect of habitat management on waterbird responses to climate warming, harvesting can also have a direct impact on population dynamics through removal of individuals [START_REF] Cooch | The effects of harvest on waterfowl populations[END_REF]). Among studied populations of non-breeding hunted ducks, several were stable or increasing over both short and long terms, such as Northern Pintail, Common Teal, Gadwall and Redcrested Pochard, while for others the trends were declining at least in one studied country, such as in Mallard, Eurasian Wigeon, Common Pochard and Tufted Duck (EEA 2020). Sustainable management of duck populations is required to maintain these species at a favorable conservation status. Management planning needs international collaborations at species distribution level, using hunting data and demographic parameters [START_REF] Holopainen | Sustainable management of migratory European ducks: finding model species[END_REF]. Whether the distribution of harvested populations is changing or not in response to climate warming should be carefully considered, to anticipate ecologic, cultural and economic challenges [START_REF] Guillemain | Effects of climate change on European ducks: what do we know and what do we need to know?[END_REF]. 2) Sites with more than 25 surveys over 28 years We performed a sensitivity analysis to ensure that the results from the main text were not affected by the number of surveys per site. We removed all sites with less than 25 surveys from this sensitivity analysis and performed the same analyses as in the main text. We investigated temporal change in Community Temperature Index (CTI), species contributions and assessed the differences of contribution between species, depending on hunting status and functional traits. We identified 288 monitored sites outside a protected area, totalling 7,703 surveys. The CTI temporal trend was positive and significant (β = 0.053 ± 0.002 SE, p < 0.001). Species contributions to the CTI trends were fairly similar to the full analysis (Fig. A3.2). Species diet, hunting status and generation length were the main factors explaining species contribution, whereas other predictors had little effect (Table A3.2). Multicollinearity was low between species diet and hunting status (VIF = 1.05). Consequently, this sensitivity analysis confirmed that species diet and hunting status are important factors explaining species contribution to CTI trend independently from site survey effort.

Bucephala clangula

Figure 1 :

 1 Figure 1: Temporal trends of (a) Community Temperature Index (CTI) ±95% CI and (b) temperature anomalies ±95% CI. Note that although temperature and CTI share the same unit, the scales are not proportional, i.e. a temperature increase of +1°C does not lead to the same CTI increase (Gaget et al. 2021).

Figure 2 :

 2 Figure2: Species contributions (%) to the Community Temperature Index (CTI) temporal trend, depending on species diet (plant, invertebrate, vertebrate, omnivore) and hunting status (hunted ducks, other hunted and non-hunted), with species ordered by their contribution value (n = 110 species). The horizontal dotted line represents a null contribution. Not all species names are shown for esthetic reason.

Figure

  Figure Appendix 3.2 : Violin plot of species contributions (%) to the Community Temperature Index (CTI) temporal trend, depending on species diet. Color scale indicates hunting status; hunted ducks (blue), other hunted species (green), non-hunted species (black). Not all species names are shown for esthetic reason.

Table 1 :

 1 Selection of the best models based on AIC and parsimony (Delta AIC <2 and lowest DF), identifying predictors the most related to species contributions. The F-value corresponds to the sum of the effects and DF is degrees of freedom. The best model to explain species contribution considers both species hunting status and diet (in bold). Models are fitted with phylogenetic generalized least squares (PGLS). Diet, hunting status and migration are factors and STI (Species Temperature Index), body mass and generation length are quantitative variables (See Appendix 4 for full model selection).

	Formula	AIC	Delta AIC Pseudo R-square	F-value DF
	Diet + Hunting status	697.2	1.4	0.058	2.329	6
	Hunting status + Diet + Generation length	695.8	0.0	0.097	6.835	7
	Hunting status + Diet + Body mass	696.3	0.5	0.077	4.410	7
	Hunting status + Diet + Migration	696.5	0.7	0.060	2.410	8

  Species contributions (%) to the Community Temperature Index (CTI) temporal trend, depending on species diet (plant, invertebrate, vertebrate, omnivore) and hunting status (hunted ducks, other hunted and non-hunted), with species ordered by their contribution value (n = 110 species). The horizontal dotted line represents a null contribution. Not all species names are shown for esthetic reason.

	Mareca strepera Bubulcus ibis Scolopax rusticola Pluvialis apricaria Figure Appendix 3.1:	Plant Invertebrate Non-hunted Hunted duck Invertebrate Other hunted Invertebrate Non-hunted	7.1 8.1 6.3 5.2	11.73496464 22.96031792 4.242742 5.356544857	Res Res Short Short	6.82 5.91 5.73 5.37
	Anas acuta Casmerodius albus Gallinago gallinago Fulica cristata	Plant Vertebrate Non-hunted Hunted duck Invertebrate Other hunted Plant Non-hunted	6.8 9.1 4.8 7	16.89718621 21.85591657 18.58383523 22.58950745	Short Short Short Res	6.85 6.77 4.74 6.72
	Anas crecca Egretta garzetta Numenius arquata Porzana porzana	Plant Invertebrate Non-hunted Hunted duck Omnivore Other hunted Omnivore Non-hunted	6.3 6.6 5 2.7	12.72483539 21.26252234 18.88781824 23.546673	Short Long Short Long	5.84 5.75 6.69 4.48
	Anas platyrhynchos Ixobrychus minutus Numenius phaeopus Gallinula chloropus	Omnivore Hunted duck Invertebrate Non-hunted Invertebrate Other hunted Omnivore Other hunted	6.6 4.1 9.1 5.9	-0.020462632 24.01168239 22.05833801 16.30709877	Short Long Long Res	6.74 4.61 5.90 5.83
	Spatula querquedula Botaurus stellaris Larus minutus Fulica atra	Omnivore Hunted duck Vertebrate Non-hunted Invertebrate Non-hunted Plant Other hunted	6.5 5.5 10.5 7	23.06512822 18.39468171 3.7 5.864663198	Long Res Short Res	5.79 7.19 4.78 6.58
	Spatula clypeata Phalacrocorax aristotelis Vertebrate Non-hunted Invertebrate Hunted duck Larus genei Omnivore Non-hunted Porphyrio porphyrio Plant Non-hunted	6.5 8.8 10.5 5.9	14.91201976 2.298892 11.89835354 22.11989646	Short Res Short Res	6.42 7.48 5.64 6.65
	Branta bernicla Phalacrocorax pygmeus Larus ridibundus Rallus aquaticus	Plant Vertebrate Non-hunted Non-hunted Invertebrate Other hunted Omnivore Other hunted	10.9 8.1 9.6 4.6	2.863382582 2.744268717 6.581840595 4.523513575	Short Short Short Res	7.15 6.52 5.65 4.72
	Branta leucopsis Phalacrocorax carbo Larus canus Grus grus	Plant Vertebrate Non-hunted Non-hunted Omnivore Non-hunted Plant Non-hunted	10.5 11.3 9.8 13.4	1.685984156 18.40826409 1.617542777 14.71592338	Short Short Short Short	7.43 7.84 6.02 8.61
	Anser albifrons Lymnocryptes minimus Larus marinus Phoenicopterus roseus	Plant Invertebrate Other hunted Non-hunted Omnivore Non-hunted Omnivore Non-hunted	11.3 5.4 12 16.3	2.533119452 20.02804147 -2.491862384 20.58416977	Short Short Short Short	7.83 3.93 7.41 8.02
	Anser fabalis Limosa limosa Larus michahellis Tachybaptus ruficollis	Plant Invertebrate Other hunted Other hunted Vertebrate Non-hunted Invertebrate Non-hunted	11.4 8.6 13 5.4	-2.514040367 21.01453888 4.57847 18.5962031	Short Long Short Short	7.92 5.67 7.00 5.14
	Anser brachyrhynchus Limosa lapponica Larus fuscus Podiceps auritus	Plant Invertebrate Non-hunted Non-hunted Vertebrate Non-hunted Vertebrate Non-hunted	11.4 8.9 13.9 7.1	2.021202386 19.22660394 18.5800917 3.393542644	Short Short Short Short	7.88 5.68 6.64 6.12
	Anser anser Arenaria interpres Larus audouinii Podiceps grisegena	Plant Invertebrate Non-hunted Other hunted Vertebrate Non-hunted Omnivore Non-hunted	11.6 7.3 13 7.1	4.465428485 17.86138751 11.45294748 2.273349328	Short Short Short Short	8.10 4.92 6.28 6.93
	Cygnus columbianus Calidris alpina Larus melanocephalus Podiceps cristatus	Plant Invertebrate Non-hunted Non-hunted Vertebrate Non-hunted Vertebrate Non-hunted	12.7 8.1 10.1 7.1	2.436899298 11.83058422 9.305989705 10.08300404	Short Short Short Short	8.76 3.97 5.55 6.60
	Cygnus cygnus Calidris maritima Rissa tridactyla Podiceps nigricollis	Plant Invertebrate Non-hunted Non-hunted Omnivore Non-hunted Invertebrate Non-hunted	12.3 7.6 12.9 7.1	-1.490706698 -2.482943324 4.798985422 17.99795243	Short Short Res Short	9.14 4.19 6.03 6.02
	Cygnus olor Calidris minuta Sterna caspia	Invertebrate Other hunted Plant Non-hunted Invertebrate Non-hunted Vertebrate Non-hunted	8 12.2 6.8 12.2	-1.10518807 1.269147723 22.77386732 21.40258537	Short Res Long Long	6.82 9.28 3.10 6.49
	Melanitta nigra Platalea leucorodia Calidris alba Sterna hirundo	Invertebrate Other hunted Invertebrate Non-hunted Invertebrate Non-hunted Vertebrate Non-hunted	7.5 7.2 8.1 11.5	2.279566759 17.95014398 18.86150911 25.24947703	Short Short Long Long	6.96 7.53 3.96 4.87
	Melanitta fusca Plegadis falcinellus Calidris temminckii Sterna sandvicensis	Invertebrate Other hunted Invertebrate Non-hunted Invertebrate Non-hunted Vertebrate Non-hunted	7.5 6.7 7.1 10.8	0.669883387 23.13696657 23.31781746 12.38123263	Short Long Long Short	7.50 6.44 3.18 5.30
	Tadorna tadorna Ciconia nigra Calidris ferruginea Chlidonias hybrida	Invertebrate Non-hunted Vertebrate Non-hunted Invertebrate Non-hunted Vertebrate Non-hunted	10.9 15.9 7.6 9.9	4.30007485 22.98346583 23.43021111 21.65184602	Short Long Long Long	7.05 7.98 4.08 4.44
	Tadorna ferruginea Ciconia ciconia Calidris canutus Charadrius alexandrinus Invertebrate Non-hunted Plant Non-hunted Vertebrate Non-hunted Invertebrate Other hunted	10.9 16.5 6.8 5	8.879664277 22.94906346 19.07810312 19.1981442	Res Long Long Short	7.12 8.15 4.96 3.77
	Netta rufina Gavia stellata Tringa ochropus Vanellus vanellus	Plant Vertebrate Non-hunted Hunted duck Invertebrate Non-hunted Invertebrate Other hunted	7 9.8 5.6 9	5.352488457 3.054206132 21.09896305 4.516681624	Short Short Short Short	7.02 7.30 4.28 5.39
	Aythya fuligula Gavia arctica Tringa erythropus Charadrius dubius	Omnivore Hunted duck Vertebrate Non-hunted Invertebrate Other hunted Invertebrate Non-hunted	7.3 9.8 5.6 5	10.45232776 0.763783716 21.03462626 23.79976584	Short Short Long Long	6.55 7.72 5.07 3.68
	Aythya marila Gavia immer Tringa nebularia Charadrius hiaticula	Invertebrate Hunted duck Vertebrate Non-hunted Invertebrate Other hunted Invertebrate Non-hunted	8.2 9.8 6.3 5.1	0.428197756 -0.043080125 23.15678175 22.10771097	Short Short Long Short	6.91 8.51 5.24 4.18
	Aythya ferina Nycticorax nycticorax Tringa glareola Recurvirostra avosetta	Plant Vertebrate Non-hunted Hunted duck Invertebrate Non-hunted Invertebrate Non-hunted	7.6 8.8 5.2 8.7	11.33614757 23.47202582 23.5308322 22.04492438	Short Long Long Short	6.71 6.70 4.14 5.72
	Aythya nyroca Ardeola ralloides Tringa totanus Himantopus himantopus Invertebrate Non-hunted Plant Non-hunted Vertebrate Non-hunted Invertebrate Other hunted	7.6 6.7 6.2 7.3	9.773350072 23.78935548 15.48429434 22.12582625	Short Long Short Long	6.35 5.66 4.87 5.18
	Marmaronetta angustirostris Ardea cinerea Tringa stagnatilis Haematopus ostralegus Invertebrate Other hunted Plant Non-hunted Vertebrate Non-hunted Invertebrate Non-hunted	7 10.3 5.6 13.7	5.953747563 18.73373742 23.22543672 15.34138378	Short Short Long Short	6.17 7.28 4.36 6.27
	Mareca penelope Ardea purpurea Actitis hypoleucos Pluvialis squatarola	Plant Vertebrate Non-hunted Hunted duck Omnivore Non-hunted Invertebrate Other hunted	6.4 10.5 6.8 6	16.53006205 23.86808726 23.15678175 18.83630341	Short Long Long Short	6.65 6.97 3.89 5.53

Table Appendix 3

 Appendix .1. Model selection identifying predictors the most related to species contributions based on AIC (delta AIC<2) and parsimony. The F-value corresponds to the sum of the effects and DF is degrees of freedom. The best model to explain species contribution considers both species hunting status, generation length and diet (in bold). Models are fitted with phylogenetic generalized least squares (PGLS). Diet, hunting status and migration are factors and STI (Species Temperature Index), body mass and generation length are quantitative variables.

				Pseudo R-
	Formula	AIC	Delta AIC	squared	F-value	DF
	Hunting status + Diet + Generation
	length	699.7	0.0	0.105	7.574	7
	Hunting status + Diet + Body mass	700.5	0.8	0.083	4.901	7
	Hunting status + Diet + Migration	701.1	1.4	0.062	2.512	8
	Diet + Hunting status	701.8	2.1	0.062	2.514	6
	Diet + Generation length	704.9	5.2	0.085	6.095	5
	Diet + Migration	705.5	5.8	0.050	1.830	6
	Diet + Body mass	705.7	6.0	0.063	3.389	5
	Diet	706.2	6.5	0.048	1.794	4
	Hunting status + Diet + STI	707.3	7.6	0.062	2.516	7
	Hunting status	709.5	9.8	0.010	0.521	3
	Diet + STI	711.7	12.0	0.049	1.836	5
	Migration	712.4	12.7	0.005	0.294	3
	Null model	713.5	13.8	0.000	0.000	1
	Body mass	714.2	14.5	0.005	0.574	2
	Generation length	714.4	14.7	0.019	2.041	2
	STI	719.1	19.4	0.001	0.020	2

Table Appendix 3

 Appendix .2. Model selection identifying predictors the most related to species contributions based on AIC (delta AIC<2) and parsimony. The F-value corresponds to the sum of the effects and DF is degrees of freedom. The best model to explain species contribution considers both species hunting status, generation length and diet (in bold). Models are fitted with phylogenetic generalized least squares (PGLS). Diet, hunting status and migration are factors and STI (Species Temperature Index), body mass and generation length are quantitative variables.

	Hunting status + Diet + STI	696.2	5.5	0.116	5.483
	Diet + Body mass	696.4	5.7	0.100	4.541
	Diet + Generation length	697.2	6.5	0.108	5.558
	Diet + STI	701.1	10.4	0.098	4.326
	Hunting status	704.4	13.7	0.010	0.543
	Migration	707.9	17.2	0.001	0.022
	Null model	708.3	17.6	0.000	0.000
	Body mass	709.5	18.8	0.001	0.073
	Generation length	711.0	20.3	0.003	0.365
	STI	712.9	22.2	0.001	1.061
	782
	783
				Pseudo R-
	Formula	AIC	Delta AIC	squared	F-value	DF
	Hunting status + Diet + Migration	690.7	0.0	0.112	4.669	8
	Hunting status + Diet + Body mass	691.0	0.3	0.122	6.173	7
	Diet + Hunting status	691.5	0.8	0.108	4.528	6
	Hunting status + Diet + Generation
	length	692.1	1.4	0.129	7.001	7
	Diet + Migration	695.5	4.8	0.094	3.654	6
	Diet	696.2	5.5	0.092	3.567	4

Acknowledgments

We are grateful to the thousands of volunteers and professionals involved in the International Waterbird Census, making this research possible. We thank M. Sorrenti and J. I. Dies for sharing their knowledge on hunting management in Italy and Spain. We thank the Turku Collegium for Science, Medicine and Technology for funding and 2 anonymous reviewers for comments that improved earlier drafts of the manuscript.

Author Contributions Statement

EG, TG and MG conceived the ideas and designed methodology; NB, TL, BL, FM, CM and MZ collected the data; EG analysed the data; EG and MG led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Appendix 2. Species functional traits and hunting status. We considered three hunting status categories between species: eleven duck species listed in Annex II of the Birds Directive, the other species listed in Annex II of the Birds Directive hunted in at least one Member State, and the non-hunted species not listed in Annex II, or listed in Annex II but not hunted in the study area. The STI was calculated as the average temperature in January (1950-2000, www.worldclim.org, spatial resolution of 0.25°) on the non-breeding species distribution (BirdLife International & HBW 2018) from [START_REF] Gaget | Waterbird communities adjust to climate warming according to conservation policy and species protection status[END_REF]. Body mass (continuous, log transformation) and diet (Plant, Omnivore, Invertebrate, Vertebrate) were extracted from [START_REF] Wilman | EltonTraits 1.0: Species-level foraging attributes of the world's birds and mammals: Ecological Archives E095-178[END_REF], generation length (continuous, in years) from IUCN (2021) and migration strategy (resident, short, long) from [START_REF] Storchová | Life-history characteristics of European birds[END_REF] Appendix 3. Sensitivity analysis.

1) Protected areas

We performed a sensitivity analysis to ensure that the results from the main text were not affected by site protection status. Indeed, protected areas can affect wintering waterbird response to climate warming, by increasing species persistence [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF]). The risk is then to observe the effect of site protection on wintering waterbird community change rather than the potential effect of hunting management. It is important however to note that hunting activities can occur inside protected areas. This means that we cannot dissociate the effect of the protected area and the potential effect of hunting management inside protected areas. To ensure that our results were not affected by protected areas, we removed all protected sites from this sensitivity analysis. Protected sites were identified following the methodology described in [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF]. We performed the same analyses, without protected sites, as in the main text. We investigated temporal change in Community Temperature Index (CTI), species contributions and assessed the differences of contribution between species, depending on hunting status and functional traits.

We identified 411 monitored sites outside a protected area, totalling 8,229 surveys. The CTI temporal trend was positive and significant (β = 0.042 ± 0.002 SE, p < 0.001). Species contributions to the CTI trends were fairly similar to the full analysis (Fig. A3.1). Species diet, hunting status and generation length were the main factors explaining species contribution, whereas other predictors had little effect (Table A3.1). Multicollinearity was low between species diet and hunting status (VIF = 1.09) and between species diet and generation length (VIF = 1.14). Consequently, this sensitivity analysis confirmed that species diet and hunting status are important factors explaining species contribution to CTI trend independently from site protection. In addition, generation length significantly contributed to explain species contributions in this subset of sites. Generation length can indeed be an important factor together with site fidelity, because juvenile settlement can influence over-winter location at adult stage, meaning that long-lived species might need more time to shift their distribution than short-lived species (see general discussion). on AIC (delta AIC<2) and parsimony. The F-value corresponds to the sum of the effects and DF is degree 785 of freedom.


Hugo Jaquard

Michaël Fanuel

Pierre-Olivier Amblard

Rémi Bardenet

Simon Barthelmé

Nicolas Tremblay

SMOOTHING COMPLEX-VALUED SIGNALS ON GRAPHS WITH MONTE-CARLO

Keywords: magnetic laplacian, spanning forests, determinantal point processes, graph smoothing, angular synchronization, ranking

We introduce new smoothing estimators for complex signals on graphs, based on a recently studied Determinantal Point Process (DPP). These estimators are built from subsets of edges and nodes drawn according to this DPP, making up trees and unicycles, i.e., connected components containing exactly one cycle. We provide a Julia implementation of these estimators and study their performance when applied to a ranking problem.

INTRODUCTION

Graph signal processing (GSP, [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending highdimensional data analysis to networks and other irregular domains[END_REF]) usually considers real data defined over the nodes of a graph G = (V, E), and classically relies on the graph Laplacian. For instance, a typical task in GSP consists in smoothing (denoising) a signal g ∈ R V by solving the penalized (Tikhonov) problem

argmin f ∈R V q f -g 2 + f Lf (1)
where L ∈ R V ×V is the graph Laplacian. We will consider throughout weighted and undirected graphs. In Equation [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending highdimensional data analysis to networks and other irregular domains[END_REF], f Lf penalises the squared-norm of the discrete derivative on the graph, i.e.:

f Lf = e=(v,v ) w e (f (v ) -f (v)) 2 ,
where the sum runs over all edges of the graph and w e > 0 is the edge weight. Thus, the quadratic form f Lf can be This work was supported by the ANR grant GRANOLA (ANR-21-CE48-0009), the ERC grant BLACKJACK (ERC-2019-STG-851866), the ANR AI chair BACCARAT (ANR-20-CHIA-0002), as well as the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01), the Grenoble Data Institute (ANR-15-IDEX-02), MIAI@Grenoble Alpes (ANR-19-P3IA-0003), the LIA CNRS/Melbourne Univ Geodesic, and the IRS (Initiatives de Recherche Stratégiques) of the IDEX Université Grenoble Alpes. thought of as i/ computing the squared difference between signal values along neighbouring edges, ii/ computing a weighted sum.

For multivariate signals, defining a discrete derivative becomes less obvious, as additional geometry enters the picture. For instance, taking f (v ) -f (v) along an edge e = (v, v ) assumes that f (v ) and f (v) use the same coordinate system, and introducing a local change of basis along the edge (v, v ) is a way to precise the relation between f (v) and f (v ) [START_REF] Singer | Vector diffusion maps and the connection laplacian[END_REF]. We focus on complex-valued signals f ∈ C V , for which a multiplication by e iθe , where θ e ∈ [0, 2π], can for instance represent a known phase offset between measurements f (v ) and f (v). In this setting, we can define a magnetic Laplacian L ∈ C V ×V [START_REF] Kenyon | Spanning forests and the vector bundle Laplacian[END_REF] (see Appendix A) which acts as follows (f * denotes the conjugate transpose of f ):

f * Lf = e=(v,v ) w e |f (v ) -e iθe f (v)| 2 .
(

) 2
Note that this equation supposes an orientation of each edge: each θ e is thus given with an orientation of e. This choice of orientation is however arbitrary: for a given edge e, one orientation associated to θ e is equivalent to the other orientation associated to -θ e as |f (v )-e iθe f (v

)| 2 = |e -iθe f (v )-f (v)| 2 .
The set {e iθe } e describes a unitary connection between the nodes of G [START_REF] Kenyon | Spanning forests and the vector bundle Laplacian[END_REF]. Such Laplacians have applications to synchronization [START_REF] Singer | Angular synchronization by eigenvectors and semidefinite programming[END_REF] and ranking problems [START_REF] Yu | Angular embedding: from jarring intensity differences to perceived luminance[END_REF][START_REF] Yu | Angular embedding: A robust quadratic criterion[END_REF][START_REF] Cucuringu | Sync-rank: Robust ranking, constrained ranking and rank aggregation via eigenvector and sdp synchronization[END_REF], as we explain in Section 4.

GSP algorithms that use the graph Laplacian often have O(|V | 3 ) scaling when implemented exactly, due to the matrix inversions or factorisations that are used (e.g., the solution of Equation (1) reads f o = q(L + qI) -1 g). In large graphs approximate methods are necessary, and [START_REF] Yusuf | Graph Tikhonov regularization and interpolation via random spanning forests[END_REF] introduced a Monte-Carlo estimator for the Tikhonov smoothing problem of Equation (1), with favourable asymptotic runtime. In this work we generalise the estimator of [START_REF] Yusuf | Graph Tikhonov regularization and interpolation via random spanning forests[END_REF] to complex signals, using a process recently introduced in [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF].

In Section 2, we introduce a slight variation of the random process of [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF], defined over graphs with a unitary connection. Our main theoretical results are given in Section 3, where we derive estimators for the solution of the Tikhonov smoothing problem. We describe a practical application to ranking in Section 4. Some technical definitions and proofs are deferred to Appendices A, B and C.

A PROCESS OVER MULTI-TYPE SPANNING FORESTS

We are interested in a distribution generalizing both the uniform distribution over spanning trees (UST) of a graph, and the random spanning forests distribution [START_REF] Avena | On some random forests with determinantal roots[END_REF]. A spanning tree is a subset of edges φ ⊆ E such that the graph with nodes V and edges φ is both connected and without cycles (cycle-free). A rooted spanning forest (RSF) φ ⊆ E ∪ V is the combination of a spanning forest and a set of distinguished nodes called the roots (one root per tree). A variation of these distributions can be defined over the set of spanning forests of unicycles (FU) [START_REF] Kenyon | Spanning forests and the vector bundle Laplacian[END_REF][START_REF] Forman | Determinants of Laplacians on graphs[END_REF], subsets of edges containing exactly one cycle per connected component, and spanning all the nodes in V . The generalization we consider instead draws its samples from the set of rooted multi-type spanning forests (MTSF). A rooted MTSF is a spanning subset of edges and nodes φ ⊆ E ∪ V whose connected components are made up of rooted trees and unicycles. These different structures are illustrated in Fig. 1. We will use the following distribution, introduced in [9], over rooted MTSFs φ = φ • ∪ ρ(φ) of G:

P(φ) ∝ r∈ρ(φ) q r e∈φ• w e C (2 -2 cos(θ C )) , (3)
with ρ(φ) the set of roots of φ, φ • its set of edges and q v ∈ R * + positive parameters associated to each node v. The third product is over the cycles of the unicycles in φ, where θ C = e∈C θ e . One may wonder why we do not consider the uniform distribution over MTSFs and instead favor sampling unicycles for which 2 -2 cos(θ C ) is large, meaning that these cycles are inconsistent. Part of the answer lies in the fact that Eq. 3 is known to describe a Determinantal Point Process (DPP) over E∪V , a useful property for proving the results in Section 3. More details are available in Appendix B.1.

Aforementioned distributions such as USTs and RSFs are conveniently sampled using (a variant of) Wilson's algorithm [START_REF] Bruce | Generating random spanning trees more quickly than the cover time[END_REF], based on random walks on the graph. These sampling procedures have been generalized to FUs [START_REF] Kassel | Random curves on surfaces induced from the Laplacian determinant[END_REF] and, recently, to MTSFs [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF], under the sampling condition:

cos(θ γ ) ≥ 0 for all cycles γ. (4)
Note that this condition applies to all cycles γ in G, and not only in some MTSF φ. The sampling algorithm is recalled in Appendix B.2. Bounds on the running time of the sampling algorithm are discussed in [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF][START_REF] Kassel | Random curves on surfaces induced from the Laplacian determinant[END_REF]. Here, we only mention that the expected running time is linear in the number of edges, and that it decreases as min v q v increases.

ESTIMATORS FOR SMOOTHING

Monte-Carlo estimators for smoothing on graphs have been developped using RSFs [START_REF] Yusuf | Graph Tikhonov regularization and interpolation via random spanning forests[END_REF]. The main idea is to sample a rooted RSF φ before propagating the value of the function at the roots to the other nodes in their associated trees. In the following, we show how this can be generalized to graphs with a unitary connections using the distribution in Equation (3). Proofs are given in Appendix C. Specifically, given g ∈ C V we derive estimators of f o = (L + Q) -1 Qg. Here Q is the diagonal matrix of the q v 's and L = D -A θ is the Hermitian magnetic Laplacian matrix, with D the diagonal degree matrix and

(A θ ) v ,v = e iθe if e = (v, v ) ∈ E (0 if e / ∈ E). When q v = q > 0 for all v ∈ V , f o
is the optimal solution of the Tikhonov problem:

argmin f ∈C V q f -g 2 + f * Lf.
The first estimator f is built by propagating values through the transport maps z → e iθe z on a rooted multitype spanning forest φ sampled according to Equation (3):

f (v) =    ψ r φ (v)→v g(r φ (v)) if the connected component of v is a rooted tree, 0 otherwise,
where r φ : V → V maps nodes to the root of the tree containing them, a → b denotes the unique path from a to b in φ and ψ a→b = e∈a→b e iθe .

Proposition 1. f is an unbiased estimator of f o :

E φ ( f ) = f o .
As a consequence of the Central Limit Theorem, Monte-Carlo estimators converge at a O( σ √ n ) rate, with σ the standard deviation of the estimator and n the number of samples. Instead of increasing the number of samples one may instead focus on decreasing the variance, and in our setting this can be done at little additional cost. A first approach is to use a Rao-Blackwell version of the estimator f (see [START_REF] Blackwell | Conditional expectation and unbiased sequential estimation[END_REF][START_REF] Rao | Information and the accuracy attainable in the estimation of statistical parameters[END_REF]) by conditioning on the set of unrooted connected components π ⊆ E of the MTSF. Consider the estimator:

f (v) =    ψ r φ (v)→v h(r φ (v)) if the connected component of v is a rooted tree, 0 otherwise,
where h(u) = w∈Tu qwψw→ug(w)

w∈Tu qw

, with T v the tree containing v in φ. Variance reduction is achieved solely from computing a mean over the nodes of the rooted trees.

Proposition 2. We have ∀v ∈ V, f (v) = E φ ( f (v)|φ • = π) and, as a consequence: E π (f ) = E φ ( f ) = f o . Also, by the law of total variance, f (v) has a lower variance than f (v): Var(f (v)) = Var( f (v))-E(Var( f (v)|φ • = π)) ≤ Var( f (v))
The method of control variates is another classical variancereduction technique for Monte-Carlo estimators, adding a term with zero mean to obtain a modified estimator with the same expectation but lower variance [START_REF] Botev | Variance reduction[END_REF]. The following is proved in [START_REF] Pilavcı | Variance reduction in stochastic methods for large-scale regularised leastsquares problems[END_REF], when q v = q for all v ∈ V (i.e. Q = qI.)

Proposition 3. Set α = 2q
q+2dm , where d m is the maximum degree in the graph. Then, the estimator

f = f -α(q -1 (L + qI)f -g) is an unbiased estimator of f o and verifies ∀v ∈ V Var( f (v)) ≤ Var(f (v))
For a graph with heterogeneous degree distribution, with e.g. a maximum degree much larger than the mean degree, f is only a marginal inprovement over f . However, on graphs with a nearly-homogeneous degree distribution, which is the case for the graphs considered in Section 4, we obtain substantial improvements over f .

RANKING FROM CORRUPTED MEASUREMENTS

In order to study the actual performance of our estimators, we describe an application to the ranking problem. We focus on the ordinal ranking problem, which asks to linearly order a set X of n elements according to an incomplete, possibly incoherent, set of pairwise ordinal comparisons C i,j ∈ {-1, 1}, for j > i. From this data, build a graph G with n nodes, and a directed edge from i to j if C i,j = 1 (resp. from j to i if C i,j = -1). Ranking according to [START_REF] Yu | Angular embedding: A robust quadratic criterion[END_REF][START_REF] Cucuringu | Sync-rank: Robust ranking, constrained ranking and rank aggregation via eigenvector and sdp synchronization[END_REF] then suggests to define the unitary connection θ i,j = πδCi,j n with δ ∈ (0, 1), and to perform angular synchronization [START_REF] Singer | Angular synchronization by eigenvectors and semidefinite programming[END_REF][START_REF] Afonso | Ten lectures and forty-two open problems in the mathematics of data science[END_REF] by solving:

argmin (ωv) v∈V e=(v,v ) |e iω v -e iθe e iωv | 2 .
(5)

The optimal arguments ω i then describe an embedding of the n points onto the unit circle, from which we can extract a ranking (see [START_REF] Cucuringu | Sync-rank: Robust ranking, constrained ranking and rank aggregation via eigenvector and sdp synchronization[END_REF] for more details).

In practice, solving such a non-convex optimization problem is difficult, and a spectral relaxation is considered instead:

min f ∈C n , f 2 =n f * Lf,
for L = D -1/2 LD -1/2 the normalized graph Laplacian, with D the diagonal degree matrix. The solution of this classical problem is the eigenvector associated to the smallest eigenvalue of L, which can for instance be computed by iterating the map x → M x/ M x for M = q( L + qI) -1 , this is the power method [START_REF] Saad | Numerical methods for large eigenvalue problems[END_REF]. Computing x → M x can either be performed directly by solving a linear system, or using the estimators f , f and f .

If we set δ = 1 4 , the sampling condition in Equation ( 4) is satisfied and fast sampling can be achieved using the variation of Wilson's algorithm recalled in Appendix B.2.

Experimental results

We illustrate in Fig. 2 some results regarding the performance of this approach obtained with our Julia implementation 1 . We work with comparisons C i,j randomly obtained from a ground-truth ranking r GT according to the Erdös-Rényi Outliers model [START_REF] Cucuringu | Sync-rank: Robust ranking, constrained ranking and rank aggregation via eigenvector and sdp synchronization[END_REF]: a comparison is observed with probability s; if a comparison is observed, it follows the ranking r GT with probability p, or is chosen uniformly in {-1, 1} otherwise. We focus on the computationally challenging dense case and set s = 0.8. Unless otherwise specified, we use q = 0.1. . The performance of the estimator f in recovering the underlying ranking, using m = 5 MTSFs for each of the k = 10 iterations of the power method, is illustrated in Figs. 2a and2b for two graphs of size n = 300 and p ∈ {0.6, 0.9}. The eigenvector of L computed from the power method without using the estimator is also plotted. Here, the initialization vector y 0 used for the power method is a random embedding of the n points in the unit circle, spaced out with angle π 2n . Runtime benchmarks are available in Figs. 2c and2d for the smoothing problem (computation of M y 0 ) solved either using f or using a Cholesky decomposition, and for the computation of the eigenvector of M for the power method iteration using f or a Cholesky decomposition, and for the Lanczos method. We display the mean running time over 100 measurements on fixed graphs of size n ∈ {10, 10 2 , 10 3 , 10 4 }. We plot Kendall's τ coefficients [START_REF] Kendall | A new measure of rank correlation[END_REF] (the larger the better) for the rankings recovered from the power method in Figs. 2e, 2f and 2g for n = 3000 and varying q, averaging over 20 realisations.

Reconstruction errors M y -e(f o ) across varying m, averaged over 20 runs for each of the estimators e(f o ) of f o , are in Fig. 2h with n = 3000.

Discussion

The results in Figs. 2a and2b compare the performance of the reconstructions obtained by the power method with and without using f , on only one realization of the graph. For these noise regimes, the ground truth ranking τ GT can no longer be recovered (which would result in a straight diagonal line). In both cases, the recovered ranking is slightly more spread out when using the estimator instead of the exact power method.

We obtain faster runtime in the smoothing problem than a standard Cholesky decomposition, in graphs with more than n ∼ 10 3 nodes (Fig. 2c). For the eigenvector computation in Fig. 2d, cross-over occurs at n ∼ 10 4 against the power method, while we do not outperform a direct Lanczos iteration. Note that performance for the eigenvector computation may be increased by sampling only one set of m = 5 MTSFs, used for all iterations of the power method. Moreover, our current implementation is far from optimal.

As q decreases, the τ coefficient increases until it reaches a fixed value in Figs. 2e, 2f and 2g. Comparing of Figs. 2e and2f shows that this behavior actually reflects the convergence of the power method, as the spectral gap of M is lower for larger q, requiring more iterations to converge. This suggests that there is a trade-off for the best choice of q: it should be as large as possible in order to reduce sampling time, and not too large so that it allows fast convergence of the power method. Note also that Kendall's τ seems to reach the value of p when convergence occurs.

Regarding the convergence of the estimators, all three versions of the estimators in Fig. 2h have linear decay in log-logscale, which is characteristic of the O(1/ √ n) convergence rate of Monte-Carlo estimators. In our simulations, the estimator f improves on the regular estimator by a factor of 10.

Throughout these experiments, we sampled m = 5 MTSFs to compute the estimators. Here, this choice was arbitrary, but we notice nonetheless that this is of the order of the O(log |V |) uniform spanning trees necessary to obtain spectral sparsifiers for connection-free graphs [START_REF] Kaufman | Scalar and matrix chernoff bounds from l ∞independence[END_REF], recently adapted to MTSFs in [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF].

CONCLUSION

We define new estimators built by propagating values along edges sampled according to a recently introduced DPP, for the smoothing problem on graphs endowed with a unitary connection, thus generalizing previous existing approaches on graphs. The evaluation of these estimators on the ranking problem, using our Julia implementation, show that they can be advantageous for smoothing starting from moderately sized graphs, as compared to a Cholesky decomposition, while they do not improve on a Lanczos iteration for eigenvector computation. Nonetheless, the proposed estimators exhibit a potentially useful computational property that is uncommon among deterministic algorithms: their computation can be carried out without any prior global knowledge of the graph, from local neighbor queries only, which in some sce-narii may be the only practical interaction. The choice of the parameters q and m is application-dependent, and requires further investigation. Possible extensions of this work include designing new applications of the proposed estimators, as well as generalizing this approach to higher-dimensional signals, where current arguments do not carry over.

A. TECHNICAL DEFINITIONS

We provide some technical definitions used in Appendices B and C.

Recall the definition of the graph Laplacian L = ∇ * ∇ [START_REF] Rk | Spectral graph theory[END_REF], with ∇ ∈ R E×V the (weighted) |E| × |V | edge-vertex incidence matrix. The graph Laplacian reveals many properties of functions f ∈ R V defined on G. For instance, if the graph is connected, the one-dimensional kernel of L is generated by the constant function.

When considering complex signals f ∈ C V on G, one first defines a unitary connection given by a family ψ v,e = e iθv,e of unitary complex numbers, extended by ψ v,e = ψ * e,v and ψ v,v = ψ e,v ψ v,e = e iθv,v . This connection gives rise to a magnetic Laplacian L = ∇ * ∇ [START_REF] Kenyon | Spanning forests and the vector bundle Laplacian[END_REF], with ∇ a twisted discrete differential:

∇ e,v =      - √ w e ψ v,e if v = a, √ w e ψ v,e if v = b, 0 otherwise,
for e = (a, b). Unlike the usual graph Laplacian, the magnetic Laplacian has a trivial kernel (ker L = 0), unless the connection is trivial [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF], which is a key difference in the analysis of the DPP described in Appendix B.1. this matrix is the marginal kernel of the process. We say that X ⊆ X is distributed according to this DPP if P(A ⊆ X) = det K A,A for every A ⊆ X , where K S,T denotes the submatrix matrix of K whose rows (resp. columns) have been restricted to S (resp. T ). If K is a projection matrix (i.e. its eigenvalues are either 0 or 1), additional structural properties are ensured. For instance, a DPP based of K will only give non-zero probability to samples of size rk(K), the rank of K. DPPs were introduced as repulsive processes [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF], and have been sought after in machine learning [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF] and randomized linear algebra [START_REF] Derezinski | Determinantal point processes in randomized numerical linear algebra[END_REF].

B. DPP AND SMOOTHING ESTIMATORS

K m = ∇ √ Q (L + Q) -1 ∇ * √ Q , with Q the |V | × |V | diagonal matrix of the q v 's.
Proof. See [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF]. This can also be proved by a straightforward extension of the computations in [START_REF] Kenyon | Spanning forests and the vector bundle Laplacian[END_REF].

This property further relates our study to that of USTs, RSFs and FUs, which have also been showed to be DPPs [START_REF] Kenyon | Spanning forests and the vector bundle Laplacian[END_REF][START_REF] Avena | On some random forests with determinantal roots[END_REF].

When working with DPPs, the Cauchy-Binet formula is an important technical tool, and states that for an m × m matrix M = AB, with A m × n, B n × m and where n > m we have:

det M = T det A :,T det B T,: ,
where the sum is over T ⊆ {1, ..., n} of size |T | = m.

B.2. Sampling

We briefly recall the sampling algorithm from [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF]. Under the condition in Equation ( 4), a MTSF φ distributed according to Equation is built by iterating the following until all nodes belong to φ.

a) Start a random walk from any node v not in φ, while keeping track of the path p built from this random walk. The transitions are as follows: with probability proportional to q, set the current node as a root. Otherwise, choose a successor u in the neighborhood of the current node w with probability proportional to w (w,u) .

b) If a node r is designated as a root, add it to the forest φ together with the path p and stop the random walk. Note that the θ C a priori depends on the orientation of the loop, but cos(θ C ) is orientation-agnostic. Such random walk based sampling algorithms are a noteworthy instance of efficient DPP sampling, which usually require a diagonalization of the kernel K and would result here in O((|E| + |V |) 3 ) sampling time.

C. PROOFS FOR SECTION 3

We provide straightforward proofs of the Propositions 1 and 2 in Section 3. Unless otherwise specified, we work with rooted MTSFs sampled according to Equation (3).

C.1. Proof of Proposition 1

The main part of the proof consists in rewriting (L + Q) -1 as an expectation (see Equation ( 7)). We first use Cramer's rule to express the coefficients of (L + Q) -1 :

(L + Q) -1 i,j = (-1) i+j det(L + Q) -j,-i det(L + Q) ,
where -k = V \ {k}.

For the denominator, we obtain the following sum over MTSFs (see e.g. [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF])

det(L + Q) = φ r∈ρ(φ) q r e∈φ• w e C (2 -2 cos(θ C )),
which is the normalizing constant in Equation ( 3). The Cauchy-Binet formula further allows to express det(L + Q) -j,-i as:

φ⊆E∪V,|φ|=|V |-1 det ∇ * √ Q -j,φ det ∇ √ Q φ,-i ,
which can be re-arranged to obtain:

φ⊆E∪V,|φ|=|V |-1 (-1) i+j det ∇ * √ Q :,φ δ j det   ∇ √ Q φ,: δ * i   , (6)
with the canonical basis vector δ i ∈ R V .

The product of determinants in Equation ( 6) can be split along the connected components φ k ⊆ E ∪ V of φ. Using similar computations as [START_REF] Kenyon | Spanning forests and the vector bundle Laplacian[END_REF], we find this product is non-zero if one of these connected components φ k is a tree containing i and j with no root for which the product of determinants is ψ j→i e∈φ k • w e and, if this component does not span the entire graph, the remaining subgraph contains disjoint connected components which have to be taken among the following subgraphs:

• a unicycle with cycle C for which the product is (2 -2 cos(θ C )) e∈φ k

• w e ,

• a rooted tree with root r for which the product is q r e∈φ k • w e . In the end, when computing the sum, this translates to:

(L + Q) -1
i,j = 1 q j φ 1 i∼j ψ j→i P(φ)

= 1 q j E(1 i∼j ψ j→i ), (7)
where 1 i∼j is the indicator that i belongs to a tree whose root is j.

Writing out the formula for the optimal solution f o finally yields:

f o (i) = δ * i (L + Q) -1 Qg = j q j (L + Q) -1
i,j g(j)

= j E(1 i∼j ψ j→i g(j))

= E( f (i)), (8)
where the second to last equation is obtained from Equation [START_REF] Cucuringu | Sync-rank: Robust ranking, constrained ranking and rank aggregation via eigenvector and sdp synchronization[END_REF], and the last line follows from the law of total expectation.

C.2. Proof of Proposition 2

Suppose that φ • = π, denote by r k the root of the connected component φ k • , and by V k the set of vertices in φ k . Then, if v belongs to π k a tree:

E φ k ( f (v) | φ k • = π k ) = r∈V k P(r k = r) f (v) = r∈V k q r
w∈V k q w ψ r→v g(r)

= f (v).

If π k is a unicycle, then E( f (v)) = 0.

The two consequences follow from the laws of total expectation and variance.

Fig. 1 :

 1 Fig. 1: Different subsets of edges and nodes, in red.

5 Fig. 2 :

 52 Fig. 2: Experimental results.

B. 1 .

 1 Background on DPP and the MTSF process A (discrete) Determinantal Point Process (DPP) is a distribution over subsets of a finite set X parameterized by a symmetric matrix K ∈ M |X | (C) whose eigenvalues all lie in [0, 1],

Proposition 4 .

 4 The probabilities in Equation (3) define a (projective) DPP over E ∪ V , with marginal kernel

c)

  If the path p self-intersects at w ∈ V before reaching a node in φ, and forms a loop C: keep this loop in p with probability 1 -cos(θ C ), add p to the forest φ and stop the random walk. Otherwise, erase this loop and continue the random walk from w. d) If p reaches a node in φ, stop the random walk and add the path p to the forest φ.

https://gricad-gitlab.univ-grenoble-alpes.fr/ tremblan/mtsf_for_graph_smoothing