# Space-time species distribution modeling for opportunistic presence-only data: a case study of passerines in a protected area

Florian Lasgorceux, Julien J. Papaix, Yoann Bunz, Damien Combrisson, Thomas Opitz

► To cite this version:

Florian Lasgorceux, Julien J. Papaix, Yoann Bunz, Damien Combrisson, Thomas Opitz. Space-time species distribution modeling for opportunistic presence-only data: a case study of passerines in a protected area. 2024. hal-04616332

HAL Id: hal-04616332

https://hal.science/hal-04616332

Preprint submitted on 18 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Space-time species distribution modeling for opportunistic presence-only data: a case study of passerines in a protected area

Florian Lasgorceux, Julien Papaïx, Yoann Bunz, Damien Combrisson, Thomas Opitz

► To cite this version:

Florian Lasgorceux, Julien Papaïx, Yoann Bunz, Damien Combrisson, Thomas Opitz. Space-time species distribution modeling for opportunistic presence-only data: a case study of passerines in a protected area. 2024. hal-04616332v2

HAL Id: hal-04616332

https://hal.science/hal-04616332v2

Preprint submitted on 1 Jul 2024 (v2), last revised 5 Dec 2025 (v3)

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Space-time species distribution modeling for opportunistic presence-only data: a case study of passerines in a protected area

Florian Lasgorceux<fnref n="1" />, Julien Papaïx<fnref n="1" />, Yoann Bunz², Damien Combrisson², and Thomas Opitz<fnref n="1" />

https://doi.org/10.5802/fake.doi

## Abstract

In recent decades, Europe has experienced a significant decline in common bird species, particularly farmland species, due to anthropic pressures like agricultural intensification. Protected areas, such as the Écrins National Park in France, play a crucial role in mitigating these impacts. In this study, we utilized an opportunistic presence-only dataset of passerine occurrences to evaluate the impact of protected areas on passerine trends. The data, gathered by skilled rangers from Écrins National Park, were used for spatio-temporal modeling to assess the target-group relative abundance of 76 passerine species across a regular grid, with aggregated occurrences per spatio-temporal cell serving as a proxy for sampling effort. The model demonstrated good calibration for most species, with AUC values from cross-validation scenarios exceeding 0.8 in most cases. The model also effectively distinguished habitat preferences and migratory status. Additionally, we compared trends in target-group relative abundance over 2001–2019 in the Écrins National Park with those across France using the French Breeding Bird Survey (STOC program) for 59 common species. Our analysis revealed that forest specialist species benefited the most from protected areas, farmland species decreased more slowly, and high-elevation specialist species generally decreased. Though, some results should be interpreted with care since our assumption of a homogeneous reporting rate across target-group species may be too strong for some very common species. Our findings illustrate that the extensive spatio-temporal resolution provided by opportunistic presence-only data offers valuable insights into biological phenomena and their trends, while reducing the need for other data sources to determine sampling effort proxies. Similar datasets from various protected areas could serve as powerful tools for assessing the effectiveness of conservation policies.

**Keywords:** Écrins National Park, INLA, habitat preference, migratory status, Poisson process, opportunistic presence-only data, relative abundance, sampling effort, spatio-temporal model, STOC program, target-group

<footnote n="1">¹Biostatistics and Spatial Processes, INRAE, Avignon, France, ²Parc national des Écrins, Domaine de Charance, Gap, France</footnote>
# Space-time species distribution modeling with opportunistic presence-only data: a case study of passerines in a protected area

Florian Lasgorceux<fnref n="1" />, Julien Papaïx<fnref n="1" />, Yoann Bunz², Damien Combrisson², and Thomas Opitz<fnref n="1" />

https://doi.org/10.17180/RAM8-MV76

## Abstract

Over the recent decades, Europe has experienced a significant decline in common bird species, particularly farmland species, due to anthropic pressures like agricultural intensification. Protected areas, such as the Écrins National Park in France, play a crucial role in mitigating these impacts. In this study, we utilized an opportunistic presence-only dataset of passerine occurrences to evaluate the impact of protected areas on passerine trends. The data, gathered by trained rangers from Écrins National Park, were used for spatio-temporal modeling to assess the relative abundance of 76 passerine species on a regular grid, with aggregated occurrences per spatio-temporal cell serving as a proxy for sampling effort. The model demonstrated good calibration for most species, with AUC values from different cross-validation schemes exceeding 0.8 in most cases. The model also effectively distinguished habitat preferences and migratory status. Additionally, we compared trends in relative abundance over 2001–2019 in the Écrins National Park with the trends in abundance across France using the French Breeding Bird Survey (STOC program) for 59 common species. Our analysis revealed that forest specialist species benefited the most from protected areas, farmland species decreased more slowly, and high-elevation specialist species generally decreased. Though, some results should be interpreted with care since our assumption of a homogeneous reporting rate across target group species may be too strong for some very common species. Our findings illustrate that the high spatio-temporal resolution provided by opportunistic presence-only data offers valuable insights into biological phenomena and their trends while reducing the need for other data sources to determine sampling effort proxies. Similar datasets from other protected areas could serve as powerful tools for assessing the effectiveness of conservation policies.

**Keywords:** Écrins National Park, INLA, habitat preference, migratory status, Poisson process, opportunistic presence-only data, relative abundance, sampling effort, spatio-temporal model, STOC program, target group

<footnote n="1">¹Biostatistics and Spatial Processes, INRAE, Avignon, France, ²Parc national des Écrins, Domaine de Charance, Gap, France</footnote>
2

# Introduction

3

4

5

6

7

8

9

10

11

12

In recent decades, Europe has seen a significant decline in the populations of common bird species (Inger et al., 2015), particularly those inhabiting agricultural landscapes (Donald et al., 2001). This sharp decline in farmland bird species has prompted extensive research, leading to numerous country-specific studies (Heldbjerg et al., 2018; Kamp et al., 2021; Newton, 2004; Reif et al., 2008b; Sanderson et al., 2013; Traba and Morales, 2019; Wretenberg et al., 2006). A large-scale study in Europe by Rigal et al., 2023 recently quantified the anthropic pressures impacting common bird abundance. They identified agricultural intensification, particularly the use of pesticides and fertilizers, as the primary factor contributing to the general decline, especially among invertebrate feeders. In contrast, responses to changes in forest cover, urbanization, and temperature showed greater variation among different species.

13

14

15

16

17

18

19

20

21

22

23

24

25

The conclusions regarding the health status of forest species in Europe are less consistent. Gregory et al., 2007 estimated that, on average, common forest birds declined by 13% and common forest specialists by 18% from 1980 to 2003. In contrast, by considering a different indicator over the period 1980-2015, Gregory et al., 2019 found a relative stability in common European forest bird populations. A possible explanation for these differences could be the migratory status of the species. Schulze et al., 2019 demonstrated that non-migratory forest bird species in Central Europe have shown significant positive population trends due to beneficial changes in forest management and environmental conditions. In contrast, birds migrating across continents have declined in recent decades, suggesting significant, contrasting changes in bird populations in Europe. Additionally, habitat preference plays a critical role. Reif et al., 2008a concluded that species closely associated with lowland broad-leaved forests have experienced more positive population trends, whereas those associated with montane and coniferous forests have shown more negative trends.

26

27

28

29

30

31

Finally, mountain bird species are particularly vulnerable, experiencing an average decline of 10% in Europe (Lehikoinen et al., 2019). This decline is driven by climate and land use changes, such as grazing pressure and afforestation, which have led to distribution shifts towards mouttaintops. Studies by Flousek et al., 2015 and (Zamora and Barea-Azcón, 2015) support this, indicating that species breeding at higher altitudes are experiencing more negative population trends, with long-distance migrants being particularly affected.

32

33

34

35

36

37

38

39

40

41

42

Protected areas and conservation efforts play a crucial role in mitigating the decline of bird populations across Europe. According to the IPBES, 2019, conservation investments from 1996 to 2008 have resulted in an average reduction of 29% in the extinction risk for mammals and birds across 109 countries. Focusing on birds across Europe, the EU's Natura 2000 protected area network, while not preventing country-wide population declines in some threatened grassland bird species (Silva et al., 2018), has been effective for bird conservation, notably for threatened species (Duckworth and Altwegg, 2018). Studies by Barnes et al., 2023 also indicate positive associations between protected areas and bird abundance, particularly benefiting rare and declining species. Timmers et al., 2022 conducted a comprehensive meta-analysis and demonstrated a strong association between strict protection measures (International Union for Conservation of Nature [IUCN] categories I-IV) and higher bird occurrence in larger forest fragments.

43

44

45

In France, a National Park is formed of two zones with different regulatory statuses: the core (regulated area) - the national Park guarantees the protection of this area, as per its decrees, and manages any human activities carried out there, in line with its management objectives - and the
2

Introduction

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

Over the recent decades, Europe has experienced a significant decline in the populations of common bird species (Inger et al., 2015), particularly those inhabiting agricultural landscapes (Donald et al., 2001). This sharp decline in farmland bird species has led to extensive research across various countries (Heldbjerg et al., 2018; Kamp et al., 2021; Newton, 2004; Reif et al., 2008b; Sanderson et al., 2013; Traba and Morales, 2019; Wretenberg et al., 2006). The conclusions regarding the health status of forest species in Europe are however less consistent. Gregory et al., 2007 estimated a decline of 13% in common forest birds and 18% in specialists from 1980 to 2003. In contrast, by considering a different indicator over the period 1980–2015, Gregory et al., 2019 found a relative stability in common European forest bird populations. migratory status may explain these differences, with non-migratory birds benefiting from improved forest management and migratory species declining (Schulze et al., 2019). Habitat preferences also play a role; species in lowland broad-leaved forests fare better than those in montane and coniferous forests (Reif et al., 2008a). This has been confirmed by the study of Lehikoinen et al., 2019 affirming that mountain bird species in Europe have declined by 10% due to climate and land use changes, pushing distributions towards mountaintops (Flousek et al., 2015; Zamora and Barea-Azcón, 2015). A recent study by Rigal et al., 2023 identified agricultural intensification as the primary factor in the general decline of common bird abundance in Europe, especially among invertebrate feeders, with varied responses to forest cover, urbanization, and temperature changes.

22

23

24

25

26

27

28

29

30

31

32

33

Protected areas and conservation efforts therefore play a crucial role in mitigating the decline of bird populations across Europe. According to the IPBES, 2019, conservation investments from 1996 to 2008 have resulted in an average reduction of 29% in the extinction risk for mammals and birds across 109 countries. Focusing on birds across Europe, the EU's Natura 2000 protected area network, while not preventing country-wide population declines in some threatened grassland bird species (Silva et al., 2018), has been effective for bird conservation, notably for threatened species (Duckworth and Altwegg, 2018). Studies by Barnes et al., 2023 also indicate positive associations between protected areas and bird abundance, particularly benefiting rare and declining species. Timmers et al., 2022 conducted a comprehensive meta-analysis and demonstrated a strong association between strict protection measures (International Union for Conservation of Nature [IUCN] categories I-IV) and higher bird occurrence in larger forest fragments.

34

35

36

37

38

39

40

41

42

43

44

National Parks in France are one of those protected areas. A National Park is formed of two geographic zones with different regulatory statuses: the core (regulated area) – the National Park guarantees the protection of this area, as per its decrees, and manages all human activities carried out there, in line with its management objectives – and the surrounding area (projects) – a space designed to achieve ecological consistency and solidarity with the National Park core. In this area, the National Park plays an advisory role and can act as a partner in projects to preserve and promote the natural, historic, cultural and landscape heritage. The Écrins National Park is also characterized by significant elevation changes, including multiple peaks exceeding 4000 meters, extensive coniferous forests (particularly larch), and numerous alpine pastures. The park is responsible for multiple biodiversity monitoring programs within its territory (Bunz, 2022), which have led to extensive synthesis and research works (Noël et al., 2023).

45

In addition to data from monitoring programs, the Écrins National Park maintains a large
46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

surrounding area (projects) - a space designed to achieve ecological consistency and solidarity with the National Park core. In this area, the National Park plays an advisory role and can act as a partner in projects to preserve and promote the natural, historic, cultural and landscape heritage. The Écrins National Park is one of those protected areas. The park's geographic features include significant elevation changes, with several peaks over 4000 meters, extensive coniferous forests (especially larch), and numerous alpine pastures. The park is responsible for multiple biodiversity monitoring programs within its territory (Bunz, 2022), which have led to extensive synthesis and research works (Noël et al., 2023).

In addition to data from monitoring programs, the Écrins National Park possesses a large dataset of presence-only records collected opportunistically by rangers. The primary use of these data is to ensure ecological surveillance with broad spatio-temporal resolution and to serve as a basis for knowledge-sharing within the context of development projects. Using these data in statistical modeling presents a challenge due to sampling variability, which occurs in space and time, and between species and observers because of their opportunistic nature (Van Strien et al., 2013). This sampling variability translates into biases in species distribution models that do not account for it. In this study, we propose to harness this extensive dataset for modeling species habitat preferences, classifying species migratory status, and studying time variations in target-group relative abundance, meaning the relative abundance of a species against a pool of species constituting the so-called target-group of species.

Various estimation frameworks and algorithms exist for constructing species distribution models and calibrating their parameters from presence-only data. Widely used methods include Maxent, based on the maximum-entropy principle (Phillips et al., 2006), and various statistical regression modeling frameworks (Valavi et al., 2022). These methods are often based on point-process representations of data (Renner et al., 2015; Warton and Shepherd, 2010), where each observed occurrence of a species is viewed as a point in a cloud of points located in space and time. In this context, models such as Maxent and logistic regression are equivalent to a special case of the point-process setting (Fithian and Hastie, 2013; Renner and Warton, 2013). Due to the uncertainty in reported positions in the Écrins National Park dataset, we rather implement a degraded version of these point process models by aggregating occurrences into spatio-temporal cells, which leads to the use of a generalized additive regression model.

Random-effect modeling are often used to capture the complex spatio-temporal structures in data and ecological processes, accurately assessing various sources of uncertainty. This has led to the widespread adoption of Bayesian inference techniques such as Markov Chain Monte Carlo (MCMC, (Gilks et al., 1995; Link et al., 2002)) or Integrated Nested Laplace Approximations (INLA, Illian et al., 2013). We apply the INLA approach to perform complex Bayesian spatio-temporal inferences with relatively large datasets by combining generalized additive regression models with the Stochastic Partial Differential Equation (SPDE) approach, the latter used to represent spatial random effects (Lindgren et al., 2011). The construction and estimation of complex spatio-temporal ecological models are still recent developments (Belmont et al., 2024; Soriano-Redondo et al., 2019), and as far as we know, our study represents the first attempt to evaluate management policies of protected areas using solely opportunistic presence-only data for this type of model.
dataset of presence-only records collected opportunistically by rangers. The primary use of these data is to ensure ecological surveillance at high spatio-temporal resolution and to serve as a basis for knowledge-sharing within the context of development projects. Using these data in statistical modeling presents a challenge due to sampling variability, which concerns the spatial and temporal dimensions, and between species and observers because of their opportunistic nature (Van Strien et al., 2013). This sampling variability translates into biases in species distribution models that do not account for it. In this study, we propose to harness this extensive dataset for modeling species habitat preferences, classifying species migratory status, and studying time variations in relative abundance, by which we mean the ratio of the abundance of a species divided by the cumulated abundance of a pool of species constituting the so-called target group of species.

Various estimation frameworks and algorithms exist for constructing species distribution models and calibrating their parameters from presence-only data. Widely used methods include Maxent, based on the maximum-entropy principle (Phillips et al., 2006), and various statistical regression modeling frameworks (Valavi et al., 2022). These methods are often based on point-process representations of data (Renner et al., 2015; Warton and Shepherd, 2010), where each observed occurrence of a species is viewed as a point in a cloud of points located in space and time. In this context, models such as Maxent and logistic regression are equivalent to a special case of the point-process setting (Fithian and Hastie, 2013; Renner and Warton, 2013). Due to the uncertainty in reported positions in the Écrins National Park dataset, we rather implement a degraded version of these point process models by aggregating occurrences into spatio-temporal cells, which leads to the use of a generalized additive regression model.

Models including random effects are often used to capture the complex spatio-temporal structures in data and ecological processes, thus allowing for an accurate assessment of various sources of uncertainty. This has led to the widespread adoption of Bayesian inference techniques such as Markov Chain Monte Carlo (MCMC, (Gilks et al., 1995; Link et al., 2002)) or Integrated Nested Laplace Approximations (INLA, Illian et al., 2013). We apply the INLA approach to perform complex Bayesian spatio-temporal inferences with relatively large datasets by combining generalized additive regression models with the Stochastic Partial Differential Equation (SPDE) approach, the latter used to represent spatial random effects (Lindgren et al., 2011). The construction and estimation of complex spatio-temporal ecological models are still recent developments (Belmont et al., 2024; Soriano-Redondo et al., 2019), and as far as we know, our study represents the first attempt to evaluate management policies of protected areas using solely opportunistic presence-only data for this type of model.

# Materials and methods

## Data

*Passerines presence-only data.* The dataset comprises 102,513 opportunistically reported events of passerine detections between 1994 and 2021, recorded by rangers of the Écrins National Park with professional training in naturalistic expertise. We chose to focus our study on the *Passeriformes* because it is a taxon generally well known by all rangers, ensuring consistent reporting rates, unlike more specific taxa. Furthermore, we have excluded three species – Common Raven (*Corvus corax*), Red-billed Chough (*Pyrrhocorax pyrrhocorax*), and Alpine Chough (*Pyrrhocorax graculus*) – due to the high uncertainty associated with their reported positions compared to the rest
# Material and methods

88

89

## Data

90

91

92

93

94

95

96

97

98

Passerines presence-only data. The dataset of interest comprises 102,513 opportunistically reported events of passerine detections between 1994 and 2021, recorded by rangers of the Écrins National Park with professional skills in naturalistic expertise. We chose to focus our study on the Passeriformes because it is a taxon generally understood by all rangers, ensuring consistent reporting rates, unlike more specific taxa. Additionally, we have excluded three species – Common Raven (Corvus corax), Red-billed Chough (Pyrrhocorax pyrrhocorax), and Alpine Chough (Pyrrhocorax graculus) – due to the high uncertainty associated with their reported positions compared to the rest of the passerines. We narrowed the temporal window to the 1994–2021 period to ensure a minimum of 2000 reported events per year (see Appendix 1 for details).

99

100

101

102

103

104

105

106

107

The data were collected opportunistically, by sight and/or ear, without a predefined survey methodology, meaning that observation and reporting of observed individuals rely on the motivation of the rangers. Each row of the dataset contains a position, a date, the name of the reported species, the ID of the observer(s), and the number of individuals. However, we choose not to consider the number of individuals because it can vary greatly due to species behavior, particularly in gregarious species, and differences in observers' methods since some observers report species abundance, while others only note presence (see the eBird example in Horns et al., 2018). Thus, we consider one reported event as the presence of at least one reported individual, resulting in a total of 102,513 passerine occurrences.

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

Based on interviews with various rangers, we have gained insight that the reported position could be i) an approximation of the reported individual's position at detection time, ii) the observer's position at detection time, or iii) another, different position where the observer gathered several occurrences before adding them jointly to the data records. Though, when multiple occurrences are gathered at a station position, rangers ensure that the environment is consistent between the station position and the bird's individual position at detection time and that the distance between these two positions does not exceed 250m. This distance of 250m is obtained by considering at most 100m for the difference between the positions of the detected individual and the observer by sound recognition (Haupert et al., 2023), plus an additional 150m for the distance between the position at detection time and the station position. To account for this uncertainty surrounding the reported position, we aggregated the data by assigning each reported position to the center of a regular spatial grid with 500m resolution covering the Écrins National Park. We excluded boundary cells with an area smaller than 25m², which corresponds to the resolution of the environmental covariates we use for characterizing the spatial habitat preferences. This discretization results in dividing the Écrins National Park domain into 11,260 spatial cells.

124

125

126

127

128

129

Our objective is not only to infer species habitat preferences but also to analyze temporal trend variations in relative abundance. Therefore, we divided the 28 years of data into seven periods of 4 years each to capture long-term changes between such periods. Additionally, in order to take into account that some studied species are migratory, we further subdivided the data temporally by month to investigate intra-annual variations. This spatio-temporal discretization scheme results in a total of 7×12×11,260 = 945,840 spatio-temporal cells.
89

90

of the passerines. We narrowed the temporal window to the 1994–2021 period to ensure a minimum of 2000 reported events per year (Appendix A, Lasgorceux et al., 2024c).

91

92

93

94

95

96

97

98

The data were collected opportunistically, by sight and/or ear, without a predefined survey methodology. Each row of the dataset contains a position, a date, the name of the reported species, the ID of the observer(s), and the number of individuals. However, we choose not to consider the number of individuals because it can vary greatly due to species behavior, particularly in gregarious species, and differences in observers' methods since some observers report species abundance, while others only note presence (see the eBird example in Horns et al., 2018). Thus, we consider one reported event as the presence of at least one reported individual, resulting in a total of 102,513 passerine occurrences.

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

Based on interviews with various rangers, we have gained insight that the reported position could be i) an approximation of the reported individual's position at detection time, ii) the observer's position at detection time, or iii) another, different position where the observer gathered several occurrences before adding them jointly to the data records. Though, when multiple occurrences are gathered at a station position, rangers ensure that the environment is consistent between the station position and the bird's individual position at detection time and that the distance between these two positions does not exceed 250m. This distance of 250m is obtained by considering at most 100m for the difference between the positions of the detected individual and the observer by sound recognition (Haupert et al., 2023), plus an additional 150m for the distance between the position at detection time and the station position. To account for this uncertainty surrounding the reported position, we aggregated the data by assigning each reported position to the center of a regular spatial grid with 500m resolution covering the Écrins National Park. We excluded boundary cells with an area smaller than 25m², which corresponds to the resolution of the environmental covariates we use for characterizing the spatial habitat preferences. This discretization results in dividing the Écrins National Park domain into 11,260 spatial cells.

115

116

117

118

119

120

Our objective is not only to infer species habitat preferences but also to analyze temporal trends in relative abundance. Therefore, we divided the 28 years of data into seven periods of 4 years each to capture long-term changes between such periods. Additionally, in order to take into account that some studied species are migratory, we further subdivided the data temporally by month to investigate intra-annual variations. This spatio-temporal discretization scheme results in a total of 7×12×11,260 = 945,840 spatio-temporal cells.

121

122

123

124

125

126

127

128

129

*Passerine population trend synthesis in France.* To compare local and national trends in relative abundance variations, we also utilized the French Breeding Bird Survey synthesis, known as the STOC program, for the period 2001–2019 (Fontaine et al., 2020; Jiguet et al., 2012). These data present the long-term abundance variation percentages of common birds in France during this period. Comparing trends from opportunistic data and STOC data serves three main purposes: (i) assessing the consistency of modeling assumptions, (ii) comparing local indicators within a protected area to global indicators across the territory of the entire country, and (iii) providing insights into space-time population trends for species specialized in mountainous habitats, for which surveys such as STOC are quite limited in terms of spatio-temporal coverage.

130

131

*Environmental Data.* The Écrins National Park is a high-mountain park where the terrain structure plays a predominant role in the ecosystems. We utilized the Digital Elevation Model data
Passerine population trends synthesis in France. To compare the local and national trends in relative abundance variations, we also utilized the French Breeding Bird Survey synthesis, known as the STOC program, for the period 2001–2019 (Fontaine et al., 2020; Jiguet et al., 2012). These data present the long-term abundance variation percentages of common birds in France during this period. Comparing trends from opportunistic data and STOC data serves three main purposes: (i) assessing the consistency of modeling assumptions, (ii) comparing local indicators within a protected area to global indicators across the territory of the entire country, and (iii) providing insights into space-time population trends for species specialized in mountainous habitats, for which surveys such as STOC are quite limited in terms of spatio-temporal coverage.

Environmental Data. The Écrins National Park is a high-mountain park where the terrain structure plays a predominant role in the ecosystems. We utilized the Digital Elevation Model data from RGE (IGN, 2018) to extract the elevation (McVicar and Körner, 2013) on a 5m regular grid. Subsequently, we calculated the slope at the same resolution using the ‘terrain’ function from the R package terra (Hijmans, 2024). Ultimately, we derived two spatial covariates: the mean elevation and the mean slope values per spatial cell at 500m resolution.

Regarding the land cover and land use data, we make use of the OSO dataset (Inglada et al., 2018). Originally at a resolution of 10/20m, we applied the nearest neighbor strategy to downscale the values to a 5m resolution, for alignment with the RGE data resolution. To avoid statistical issues due to uncommon land cover and land use types in the Écrins National Park, we aggregated data towards 10 representative categories: Urban, Crops, Meadows, Deciduous forest, Coniferous forest, Grasslands, Woody Heaths, Mineral Surfaces, Glaciers and snow, and Water (see Appendix 2 for further details). We then calculated the percentage of coverage of each category in each spatial cell at 500m resolution.

We also used the Historical Monthly Weather Data from WorldClim (CRU-TS 4.06 (Harris et al., 2020) downscaled with WorldClim 2.1 (Fick and Hijmans, 2017)) at 2.5 minutes spatial resolution (roughly 5km) between the years 1994 and 2021. The available variables were average minimum temperature (°C), average maximum temperature (°C), and total precipitation (mm). We extracted the corresponding values in each spatial cell, for each month, and each year. Since the inter-annual resolution of our spatio-temporal cells consists of periods of 4 years, we computed the mean of the extracted values of the 4 corresponding years for a given spatial cell and a given month. Finally, we removed the average minimum temperature covariate from our study due to its high correlation with the average maximum temperature (Pearson's R > 0.99).

Preprocessing of model covariates. Elevation shows strong variability in the Écrins National Park with certain climate and land cover configurations arising predominantly within specific elevation ranges. As a result, both land cover and climatic variables are correlated.

Hence, we began by decomposing the climatic variables into their spatial means and spatio-temporal anomalies with respect to their means. We then performed a Principal Component Analysis (PCA, Janžekovič and Novak, 2012) using the R package ade4 (Dray and Dufour, 2007). The input variables were: the mean elevation, the mean slope, the percentages of coverage computed for each of the ten representative categories derived from the OSO classification, as well as the spatial means of average maximum temperature and total precipitation. The first 6 axes of the PCA, accounting for 75% of the information, allow for useful interpretation and have
from RGE (IGN, 2018) to extract the elevation (McVicar and Körner, 2013) on a 5m regular grid. Subsequently, we calculated the slope at the same resolution using the 'terrain' function from the R package *terra* (Hijmans, 2024). Ultimately, we derived two spatial covariates: the mean elevation and the mean slope values per spatial cell at 500m resolution.

Regarding land cover and land use data, we make use of the OSO dataset (Inglada et al., 2018). Originally at a resolution of 10/20m, we applied the nearest neighbor strategy to downscale the values to a 5m resolution, for alignment with the RGE data resolution. To avoid statistical issues due to uncommon land cover and land use types in the Écrins National Park, we aggregated data towards 10 representative categories: Urban, Crops, Meadows, Deciduous forest, Coniferous forest, Grasslands, Woody Heaths, Mineral Surfaces, Glaciers and snow, and Water (see Appendix B in Lasgorceux et al., 2024c for further details). We then calculated the percentage of coverage of each category in each spatial cell at 500m resolution.

We also used the Historical Monthly Weather Data from WorldClim (CRU-TS 4.06 (Harris et al., 2020) downscaled with WorldClim 2.1 (Fick and Hijmans, 2017)) at 2.5 minutes spatial resolution (roughly 5km) between the years 1994 and 2021. The available variables were average maximum temperature (°C) and total precipitation (mm). We extracted the corresponding values in each spatial cell, for each month, and each year. Since the inter-annual resolution of our spatio-temporal cells consists of periods of 4 years, we computed the mean of the extracted values of the 4 corresponding years for a given spatial cell and a given month.

*Preprocessing of model covariates.* Elevation shows strong variability in the Écrins National Park with certain climate and land cover configurations arising predominantly within specific elevation ranges. As a result, both land cover and climatic variables are correlated.

Hence, we began by decomposing the climatic variables into their spatial means and spatio-temporal anomalies with respect to those means. We then performed a Principal Component Analysis (PCA, Janžekovič and Novak, 2012) using the R package *ade4* (Dray and Dufour, 2007). The input variables were: the mean elevation, the mean slope, the percentages of coverage computed for each of the ten representative categories derived from the OSO classification, as well as the spatial means of average maximum temperature and total precipitation. The first 6 axes of the PCA, accounting for 75% of the information, allow for useful interpretation and have therefore been incorporated into the species distribution models (see Appendix C in Lasgorceux et al., 2024c for details). The first principal component axis explains a significant 31.1% of the variation in habitat types, primarily driven by high-elevation environments. The second axis is characterized by vegetation at medium to high altitudes, particularly in coniferous forests, especially in steep areas. The third axis is influenced by grasslands, and the fourth axis is primarily driven by wetlands and glaciers. However, since we observe a significant predominance of bird occurrences in "water" habitats within urbanized valleys (Figure 14 in Appendix C, Lasgorceux et al., 2024c) compared to "glaciers and snow" habitats (751 occurrences versus 80), we interpret the fourth axis as "urban wetlands". The interpretation of the fifth axis is also associated with high-elevation vegetation, as it is primarily driven by dense, low-growing shrubs and small trees; henceforth, we will refer to these environments as "woody heaths". The sixth axis is driven by low-elevation deciduous forests.

This approach results in the implementation of eight covariates in our statistical model: six spatial covariates from the PCA analysis, and two spatio-temporal covariates defined as the anomalies of average maximum temperature and anomalies of total precipitation.
therefore been incorporated into the species distribution models (Appendix 3). The first principal component axis explains a significant 31.1% of the variation in habitat types, primarily driven by high-elevation environments, especially in steep areas. The third axis is influenced by grasslands, and the fourth axis is primarily driven by wetlands and glaciers. However, since we observe a significant predominance of bird occurrences in “Water” habitats within urbanized valleys (Figure 13 in Appendix 3) compared to “Glaciers and snow” habitats (751 occurrences versus 80), we interpret the fourth axis as ‘urban wetlands’. The interpretation of the fifth axis is also associated with high-altitude vegetation, as it is primarily driven by dense, low-growing shrubs and small trees; henceforth, we will refer to these environments as “woody heaths”. The sixth axis is driven by low-elevation deciduous forests.

This results in the implementation of eight covariates in our statistical model: six spatial covariates from the PCA analysis replicated for each month and each period, and two spatio-temporal covariates defined as the anomalies of average maximum temperature and anomalies of total precipitation.

## Statistical modeling

We applied the same Bayesian hierarchical model (Wikle, 2003) to each species with more than 100 occurrences over the period from 1994 to 2021, excluding *Bombycilla garrulus* where 90% of its 146 occurrences were reported in 2005, due to a brief invasion event of this species within the Écrins National Park that year. This resulted in a total of 76 species being analyzed.

**Data model.** We denote by $s$ a spatial cell, by $m$ a month, and by $p$ a temporal period of 4 years, with $(s, m, p)$ representing the corresponding space-time cell. The response variable, denoted as $Y_i(s, m, p)$, is the number of reported occurrences of species $i$ from the opportunistic presence-only dataset. We model $Y_i(s, m, p)$ with a hierarchical generative additive mixed model (Knape, 2016; Pedersen et al., 2019) using a Poisson distribution:

$$ (1) \qquad Y_i(s, m, p) | \mu_i(s, m, p) \sim \text{Poisson}(\mu_i(s, m, p)), $$

where $\mu_i(s, m, p)$ is the average number of reported occurrences of species $i$ in the cell $(s, m, p)$.

We assume that $\mu_i(s, m, p)$ is obtained through the combination of two factors: (i) the abundance of species $i$ in the cell $(s, m, p)$, denoted by $\Lambda_i(s, m, p)$, and (ii) the sampling effort in cell $(s, m, p)$ for species $i$, denoted by $E_i(s, m, p)$, which we assume to both act in a multiplicative way on $\mu_i(s, m, p)$ (Giraud et al., 2016). These relationships lead us to formulate the following equation:

$$ (2) \qquad \mu_i(s, m, p) = \Lambda_i(s, m, p) \times E_i(s, m, p) $$

**Sampling effort.** A major challenge is now to quantify or eliminate the term representing the sampling effort. When working with passerine bird survey data, parameters such as the number of listening points, the duration of time spent at each point, distances, study areas, etc., are pre-defined and meticulously recorded. This ensures that data sampling remains consistent across different zones or years. With opportunistic data, this is not the case. Worse, sampling methods differ among observers, causing variations in the frequency of species reporting. This occurs because observers have the freedom to report sightings as they wish, resulting in significant differences in detection rates between observers and species. If not taken into account in the statistical modeling, this heterogeneous sampling will result in bias and wrong conclusions.
176 Statistical modeling

We applied the same Bayesian hierarchical model (Wikle, 2003) to each species with more than 100 occurrences over the period from 1994 to 2021, excluding *Bombycilla garrulus* for which 90% of its 146 occurrences were reported in 2005, due to a brief invasion event of this species within the Écrins National Park that year. This resulted in a total of 76 species being analyzed.

181 *Data model.* We denote by *s* a spatial cell, by *m* a month, and by *p* a temporal period of 4 years, 182 with the triplet (*s*, *m*, *p*) representing the corresponding space-time cell. The response variable, 183 denoted as $Y_i(s, m, p)$, is the number of occurrences of species *i* reported for that cell in the op- 184 portunistic presence-only dataset. We model $Y_i(s, m, p)$ with a hierarchical generalized additive 185 mixed model (Knape, 2016; Pedersen et al., 2019) using a Poisson response distribution:

$$ (1) \qquad Y_i(s, m, p) | \mu_i(s, m, p) \sim \text{Poisson}(\mu_i(s, m, p)), $$

where $\mu_i(s, m, p)$ is the average number of reported occurrences of species *i* in the cell ($s, m, p$).

We assume that $\mu_i(s, m, p)$ is obtained through the combination of two factors: (i) the abundance of species *i* in the cell ($s, m, p$), denoted by $\Lambda_i(s, m, p)$, and (ii) the sampling effort in cell ($s, m, p$) for species *i*, denoted by $E_i(s, m, p)$, which we assume to both act in a multiplicative way on $\mu_i(s, m, p)$ (Giraud et al., 2016). These relationships lead us to formulate the following equation:

$$ (2) \qquad \mu_i(s, m, p) = \Lambda_i(s, m, p) \times E_i(s, m, p) $$

192 *Sampling effort.* A major challenge is now to quantify or eliminate the term representing the 193 sampling effort. Within our spatio-temporal discretized framework, the first step is to obtain in- 194 formation about where and when no observations occurred. This helps distinguish between cells 195 with zero occurrences due to the absence of observer presence and those with zero occurrences 196 due to the absence of the species. The strategy involves utilizing the extensive opportunistic 197 presence-only dataset, using reported occurrences of certain species as indicators of observer 198 presence, as proposed by Botella et al., 2020. However, it is necessary to assume that these 199 species are uniformly sampled by all the observers, corresponding to homogeneity in reporting 200 rates across the species pool, as well as in the probability of detection. This explains one of our 201 reasons to primarily restrict our study to passerines (see the subsection *Passerine data*). Following 202 existing literature, we refer to this set of chosen species with homogeneous sampling conditions 203 as the target group (TG) of species. Therefore, we denote the sampling effort term as *E* instead 204 of $E_i$, and we write $\Lambda_{TG}(s, m, p)$ for the cumulative abundance of all species in the target group 205 in cell ($s, m, p$). We can then reformulate the average number of occurrences in terms of the 206 ratio of the abundance of species *i* to the cumulative abundance of the target group, and the 207 product of the sampling effort with the cumulative abundance of the target group. (Equation 4):

$$ (3) \qquad \mu_i(s, m, p) = \Lambda_i(s, m, p) \times E_i(s, m, p) $$

$$ (4) \qquad = \frac{\Lambda_i(s, m, p)}{\Lambda_{TG}(s, m, p)} \times \Lambda_{TG}(s, m, p)E(s, m, p). $$

To statistically identify the first term that indicates the abundance of species *i* relative to 208 the cumulative abundance of the target group of species -hereafter referred to as relative 209 abundance- we need a proxy for the second term $\Lambda_{TG}(s, m, p)E(s, m, p)$. The most natural ap- 210 proach is to use the number of target-group reported occurrences, denoted as $y_{TG}(s, m, p)$. We 211
Within our spatio-temporal discretized framework, the first step is to obtain information about where and when no observations occurred. This helps distinguish between cells with zero occurrences due to the absence of observer presence and those with zero occurrences due to the absence of the species. The strategy involves utilizing the extensive opportunistic presence-only dataset, using reported occurrences of certain species as indicators of observer presence, as proposed by Botella et al., 2020. However, it is necessary to assume that these species are uniformly sampled by all the observers, corresponding to homogeneity in reporting rates across the species pool, as well as in the probability of detection. This explains one of our reasons to primarily restrict our study to passerines (see the Subsection *Passerines data*). Following existing literature, we refer to this set of chosen species with homogeneous sampling conditions as the target-group (TG) of species. Therefore, we denote the sampling effort term as *E* instead of *E*<sub>i</sub>, and we write Λ<sub>TG</sub>(s, m, p) for the cumulative abundance of all species in the target-group in cell (s, m, p). We can then reformulate the average number of occurrences in terms of the abundance of species *i* compared to the cumulated abundance of the target-group, and the product of the sampling effort by the cumulated abundance of the target-group (Equation 4).

$$ (3) \qquad \mu_i(s, m, p) = \Lambda_i(s, m, p) \times E(s, m, p) $$

$$ (4) \qquad = \frac{\Lambda_i(s, m, p)}{\Lambda_{TG}(s, m, p)} \times \Lambda_{TG}(s, m, p)E(s, m, p). $$

To be able to statistically identify the first term that informs us about target-group relative abundance of species *i*, we now need a proxy of the second term Λ<sub>TG</sub>(s, m, p)E(s, m, p). The most natural approach is to use the number of reported occurrences of the target-group occurrences, denoted as y<sub>TG</sub>(s, m, p). We chose to use the combined occurrences of all passerine species, which was consistent with algorithms similar to that of Botella et al., 2020 since they were all selected except for 4 species observed only once [à modifier : (DOI, Data, scripts, code, supplementary materials.)] We kept only the cells with at least one occurrence of the target group, reducing the number of studied cells to 32,486.

*Target-group relative abundance.* Employing a spatio-temporal distribution model enables the implementation of numerous spatio-temporally structured latent fixed or random effects within the linear structure that we will assume for the predictor log(μ<sub>i</sub>(s, m, p)). We opted to use only one generic model for each of the 76 species to automate model fitting and facilitate inter-species comparisons. In Poisson models, it is common practice to distribute the Poisson intensity, denoted here as μ<sub>i</sub>(s, m, p), using a log-link function to ensure that μ<sub>i</sub>(s, m, p) remains positive. Focusing on the parameter of interest, the target-group relative abundance (Λ<sub>i</sub>(s,m,p)/Λ<sub>TG</sub>(s,m,p)), we constructed our linear regressor by incorporating the proxy of the sampling effort, y<sub>TG</sub>(s, m, p), as an offset (Equation 5).

By studying the abundance of a passerine compared to all others, we compare species with different biological characteristics, such as habitats and migration status, which need to be taken into account in the linear regression. We denoted α<sub>i</sub> as the intercept. We examined a linear effect of the spatial covariates X<sup>(PCA)</sup>(s) (PCA components) and the spatio-temporal covariates X<sup>(anomalies)</sup>(s, m, p) (anomalies of climatic covariates) to estimate habitat preferences of species *i*, relative to other passerines. The term W<sub>i</sub><sup>(m)</sup>(s), m = 1, ..., 12, represents a monthly latent spatio-temporal Gaussian field. Concretely, this field aims to model intra-annual spatial variations not
chose to use the combined occurrences of all passerine species, which was consistent with algorithms similar to that of Botella et al., 2020 since they were all selected except for 4 species observed only once. We kept only the cells with at least one target-group occurrence, reducing the number of studied cells to 32,486.

**Relative abundance.** Employing a spatio-temporal distribution model enables the implementation of various spatio-temporally structured latent fixed or random effects within the linear structure that we will assume for the predictor $log(\mu_i(s, m, p))$. We opted to use only one generic model for each of the 76 species to automate model fitting and facilitate inter-species comparisons. In Poisson models, it is common practice to represent the Poisson intensity, denoted here as $\mu_i(s, m, p)$, using a log-link function to ensure that $\mu_i(s, m, p)$ remains positive. Focusing on the parameter of interest, the relative abundance $\frac{\Lambda_i(s, m, p)}{\Lambda_{TG}(s, m, p)}$, we constructed our linear regressor by incorporating the proxy of the sampling effort, $y_{TG}(s, m, p)$, as a fixed offset (Equation 5).

By studying the abundance of a passerine species compared to all others, we compare species with different biological characteristics, such as habitats and migratory status, which need to be taken into account in the linear regression. We therefore allowed for a linear effect of the spatial covariates $X^{(PCA)}(s)$ (PCA components) and the spatio-temporal covariates $X^{(anomalies)}(s, m, p)$ (anomalies of climatic covariates) to estimate habitat preferences of species $i$, relative to other passerines. We denoted $\alpha_i$ as the intercept. The term $W_i^{(m)}(s)$, $m = 1, ..., 12$, represents a monthly latent spatio-temporal Gaussian field. Concretely, this field aims to model intra-annual spatial variations not captured by the covariates, typically variations in relative abundance due to the migration of certain passerines within the territory of the Écrins National Park. Lastly, we denote by $f_i(p)$ a non-spatialized inter-annual effect, aiming to capture the variations in relative abundance across each period $p$ of 4 years.

$$ (5) \quad \log(\mu_i(s, m, p)) = \log(y_{TG}(s, m, p)) + \log\left(\frac{\Lambda_i(s, m, p)}{\Lambda_{TG}(s, m, p)}\right) $$

$$ (6) \quad = \log(y_{TG}(s, m, p)) + \alpha_i + X^{(PCA)}(s)\beta_i + X^{(anomalies)}(s, m, p)\gamma_i + W_i^{(m)}(s) + f_i(p). $$

**Statistical inference.** We employed the R-INLA package for the statistical inference of our models (Rue et al., 2009). INLA is a fast but accurate estimation method based on deterministic Laplace approximations –by contrast with approaches based on stochastic simulation, such as MCMC– and requires that the latent layer ($log(\mu_i(s, m, p))$ in our model) follows a multivariate Gaussian distribution, possibly in very high dimension, as with our space-time fields. In this setting, INLA enables the implementation of a wide range of Gaussian random effects, including spatial fields implemented through the SPDE approach, which provides a numerically convenient representation of the Matérn covariance function used routinely in spatial statistics (Lindgren et al., 2011), and of priors for hyperparameters such as variances or correlation ranges. In Appendix D, we precisely describe the prior distributions of all parts of $log(\mu_i(s, m, p))$, i.e., the structure of the model components before updating them with information from observation data to obtain their posterior distributions (Lasgorceux et al., 2024c).

## Validation and post-processing

The model presented previously is quite comprehensive, allowing for a thorough examination of information on a species-by-species basis. In the following sections, we delineate the post-processing phase undertaken to synthesize this information. We address questions pertaining
captured by the covariates, typically variations in target-group relative abundance due to the migration of certain passerines within the territory of the Écrins National Park. Lastly, we denote by $f_i(p)$ a non-spatialized inter-annual effect, aiming to capture the variations in target-group relative abundance across each period $p$ of 4 years.

$$ (5) \quad \log(\mu_i(s, m, p)) = \log(y_{TG}(s, m, p)) + \log\left(\frac{\Lambda_i(s, m, p)}{\Lambda_{TG}(s, m, p)}\right) $$

$$ (6) \quad = \log(y_{TG}(s, m, p)) + \alpha_i + \mathbf{X}^{(PCA)}(s)\beta_i + \mathbf{X}^{(anomalies)}(s, m, p)\gamma_i + W_i^{(m)}(s) + f_i(p). $$

**Statistical inference.** We employed the INLA package for the statistical inference of our models (Rue et al., 2009). INLA is a fast but accurate estimation method based on deterministic Laplace approximations –by contrast with approaches based on stochastic simulation, such as MCMC– and requires that the latent layer ($\log(\mu_i(s, m, p))$ in our model) follows a multivariate Gaussian distribution, possibly in very high dimension. In this setting, INLA enables the implementation of a wide range of Gaussian random effects, including spatial fields implemented through the SPDE approach, which provides a numerically convenient representation of the Matérn covariance function used routinely in spatial statistics (Lindgren et al., 2011), and of priors for hyperparameters such as variances or correlation ranges. In Appendix 4, we precisely describe the prior distributions of all parts of $\log(\mu_i(s, m, p))$, i.e., the structure of the model components before updating them with information from observation data to obtain their posterior distributions.

## Validation and post-processing

The model presented previously is quite comprehensive, allowing for a thorough examination of information on a species-by-species basis. In the following sections, we delineate the post-processing phase undertaken to synthesize this information. We address questions pertaining to the model's validation and also delve into the three components of the linear predictor (Equation 6): those concerning the fixed effects – specifically, habitat preferences compared to the target-group denoted by $\beta_i$ and $\gamma_i$; the spatial intra-annual component, characterized by the month-dependent Gaussian field $W_i^m$, which is associated with migratory status; and the inter-annual component, represented by the random i.i.d. Gaussian effect $f_i(p)$, which captures long-term trends in target-group relative abundance.

**Computational cost.** For each species, the model comprises 32,486 space-time cells with observations, and a total of 7,676 latent variables to estimate across the prediction cells. We conducted analyses using the R programming language (R Core Team, 2024) and utilized version 24.02.09 of INLA (Rue et al., 2009), leveraging the PARDISO 8.2 library to accelerate computations (Eftekhari et al., 2021; Gaedke-Merzhäuser et al., 2023; Pasadakis et al., 2023). This configuration typically yields computation times averaging around one minute for each model fit on an 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz processor.

**Validation.** We evaluated whether our model effectively extracted meaningful insights from the opportunistic data by comparing its predictive performance against a baseline that considered only sampling effort. We focus on the prediction of the binary event of presence or absence of a species in a given cell and month and compute the Area-Under-the-Curve (AUC, Fawcett, 2006; Jin Huang and Ling, 2005) (where the curve is the so-called Receiver Operating Characteristic curve) as a simple scalar prediction score for each species using two predictors: (i) the sampling effort proxy offset and (ii) the full linear predictor (Equation 6). Next, we assessed the model's
to the model's validation and also delve into the three components of the linear predictor (Equation 6): those concerning the fixed effects – specifically, habitat preferences compared to the target group denoted by $\beta_i$ and $\gamma_i$; the spatial intra-annual component, characterized by the month-dependent Gaussian field $W_i^m$, which is associated with migratory status; and the inter-annual component, represented by the random i.i.d. Gaussian effect $f_i(p)$, which captures long-term trends in relative abundance.

**Computational cost.** For each species, the model comprises 32,486 space-time cells with observations and a total of 7,676 latent variables to estimate across the prediction cells. We conducted analyses using the R programming language (R Core Team, 2024) and utilized version 24.02.09 of R-INLA (Rue et al., 2009), leveraging the PARDISO 8.2 library to accelerate computations (Eftekhari et al., 2021; Gaedke-Merzhäuser et al., 2023; Pasadakis et al., 2023). This configuration typically yields computation times averaging around one minute for each model fit on a standard workstation with 2.80GHz processor.

**Validation.** We evaluated whether our model effectively extracted meaningful insights from the opportunistic presence-only data by comparing its predictive performance against a baseline that considered only the sampling effort proxy, i.e. the target-group occurrences as an offset. We focus on the prediction of the binary event of the reported presence or non-reported presence of a species in a given cell and month and compute the Area-Under-the-Curve (AUC, Fawcett, 2006; Jin Huang and Ling, 2005) (where the curve is the so-called Receiver Operating Characteristic curve) to obtain a simple scalar prediction score for each species, where we compare results for two predictors: (i) the sampling effort proxy offset and (ii) the full linear predictor (Equation 6).

Next, we assessed the model's generalization skill for predicting new data using k-fold, spatial, and temporal cross-validation scenarios. In k-fold cross-validation, the dataset is divided into k equal parts, with one part used for validation while the rest is used for training (Jung, 2018). This approach yields $k + 1$ AUC scores per species: $k$ scores from predicting each part separately, and one from predicting the entire dataset. We selected $k = 10$ to ensure comparability with temporal and spatial scenarios, where the training data represented between 3.9% and 13.7% of the total volume. This choice aimed to maintain consistency in the quantity of test data across scenarios, ensuring equitable evaluation conditions. In the spatial scenario, we assessed model robustness by alternately removing all spatio-temporal cells within a given municipality of the Écrins National Park. Seven municipalities were selected based on geographic zones and data availability, resulting in seven AUC scores per species. For the temporal scenario, we alternately randomly removed half of the data for a given period, resulting in seven AUC scores per species. Further details on the size of the training datasets are given in Appendix E in Lasgorceux et al., 2024c.

**Fixed effects.** We estimated the *a posteriori* distribution of coefficients $\beta_i$ and $\gamma_i$ separately for each species *i* (Equation 6). The posterior mean of each parameter and its standard quantity the importance of the covariates in explaining occurrences of species *i* relative to the occurrences of all the passerines taken together. The credibility interval of each parameter helps determine how strongly the corresponding covariate influences the species; if the interval does not contain zero, we can consider the covariate as being “significant” (by analogy with the frequentist notion of significance). To synthesize the estimated mean coefficients and assess parameter similarities across species occupying comparable ecological habitats, we conducted Principal Component
generalization ability for predicting new data using k-fold, spatial, and temporal cross-validation scenarios. In k-fold cross-validation, the dataset is divided into k equal parts, with one part used for validation while the rest is used for training (Jung, 2018). This approach yields $k + 1$ AUC scores per species: k scores from predicting each part separately, and one from predicting the entire dataset. We selected $k = 10$ to ensure comparability with temporal and spatial scenarios, where the training data represented between 3.9% and 13.7% of the total volume. This choice aimed to maintain consistency in the quantity of test data across scenarios, ensuring equitable evaluation conditions. In the spatial scenario, we assessed model robustness by alternately removing all spatio-temporal cells within a given municipality of the Écrins National Park. Seven municipalities were selected based on geographic zones and data availability, resulting in seven AUC scores per species. For the temporal scenario, we alternately randomly removed half of the data for a given period, resulting in seven AUC scores per species. Further details on the size of the training datasets are given in Appendix 5.

**Fixed effects.** We estimated the *a posteriori* distribution of coefficients $\beta_i$ and $\gamma_i$ separately for each species *i* (Equation 6). The posterior mean of each parameter and its standard quantile quantify the importance of the covariates in explaining occurrences of species *i* relative to the occurrences of all the passerines taken together. The credibility interval of each parameter helps determine how strongly the corresponding covariate influences the species; if the interval does not contain zero, we can consider the covariate as being “significant” (by analogy with the frequentist notion of significance). To synthesize the estimated mean coefficients and assess parameter similarities across species occupying comparable ecological habitats, we conducted Principal Component Analysis (PCA) on estimated parameters followed by species clustering based on the estimated mean effects. This clustering was performed using a hierarchical agglomerative clustering (HAC) algorithm implemented through the *hclust* function, employing the *ward.D2* method (Murtagh and Legendre, 2014).

**Month-based spatio-temporal Gaussian field.** The month-based spatio-temporal Gaussian field $W_i^{(m)}$ adds flexibility to the model but is challenging to summarize graphically since we have 12 maps per species. Therefore, we opted to reduce dimensionality by spatially averaging each field. For each species, this yields a curve of 12 points depicting the intra-annual variations in occurrences relative to the target-group. Instead of plotting all 76 curves on a single graph, we conducted a Functional Principal Component Analysis (FPCA, Ramsay and Silverman, 2005). The mechanism of the approach is analogous to PCA but adapted to functions instead of vectors. In our case, FPCA decomposes the intra-annual mean effect of a given species into (i) the mean intra-annual effect across all species plus (ii) the remaining signal represented by scores associated with empirically identified harmonics. In a functional principal component analysis (FPCA), harmonics represent the most significant modes of variation in the observed functions. These harmonics are akin to principal components in traditional PCA but are applied to functions rather than conventional variables, where the difference is that functions are represented as vectors with a very large number of components (larger than the sample size of functions), and typically show correlation for components at near positions within the vector. We utilized the *fda* package to conduct these analyses (Ramsay et al., 2009).

**Inter-annual effects.** Our motivation for introducing the inter-annual effect in the linear predictor (Equation 6) lies in comparing the target-group relative abundance trends of passerines in the
Analysis (PCA) on estimated parameters followed by species clustering based on the estimated mean effects. This clustering was performed using a hierarchical agglomerative clustering (HAC) algorithm implemented through the `hclust` function, employing the `ward.D2` method (Murtagh and Legendre, 2014).

Month-based spatio-temporal Gaussian field. The month-based spatio-temporal Gaussian field $W_i^{(m)}$ adds flexibility to the model but is challenging to summarize graphically since we have 12 maps per species. Therefore, we opted to reduce dimensionality by spatially averaging each field. For each species, this yields a curve of 12 points depicting the intra-annual variations in occurrences relative to the target group. Instead of plotting all 76 curves on a single graph, we conducted a Functional Principal Component Analysis (FPCA, Ramsay and Silverman, 2005). In our case, FPCA decomposes the intra-annual mean effect of a given species into (i) the mean intra-annual effect across all species plus (ii) the remaining signal represented by scores associated with empirically identified harmonics. In a functional principal component analysis, harmonics represent the most significant modes of variation in the observed functions. These harmonics are akin to principal components in traditional PCA but are applied to functions rather than conventional variables, where the difference is that functions are represented as vectors with a very large number of components (larger than the sample size of functions), and typically show correlation for components at near positions within the vector. We utilized the `fda` package to conduct these analyses (Ramsay et al., 2009).

Inter-annual effects. Our motivation for introducing the inter-annual effect in the linear predictor (Equation 6) lies in comparing the relative abundance trends of passerines in the Écrins National Park with the abundance trends of common birds in the STOC. If we denote by $\lambda_i^{(STOC)}(y)$ the estimated abundance of species i in year y from the STOC, then the percentage variation in abundance computed by the STOC (Fontaine et al., 2020) can be expressed as

$$ (\hat{\lambda}_i^{(STOC)}(2019) - \hat{\lambda}_i^{(STOC)}(2001)) \times \frac{100}{\hat{\lambda}_i^{(STOC)}(2001)} $$

In the model presented in Equation 6, the term capturing the effect of the period p on the relative abundance of species i is $\exp(f_i(p))$. The INLA method allows us to estimate a posterior distribution of each $\exp(f_i(p))$. To mimic the STOC period of study, we assume the relative abundance to be constant in each period p. This allows us to attribute an estimated relative effect $f_i(y)$ to each year y between 2001 and 2019. We then calculated a regression line using a linear model with the $f_i(y)$ as the response variable and y the explanatory variable. If we denote by $\hat{f}_i^{(ENP)}(y)$ the fitted value of this linear regression, we could calculate the percentage variations in relative abundance of species i as

$$ (\hat{f}_i^{(ENP)}(2019) - \hat{f}_i^{(ENP)}(2001)) \times \frac{100}{\hat{f}_i^{(ENP)}(2001)} $$

To enhance robustness and take into account statistical uncertainty, we generated 1000 posterior distributions for each relative abundance trend parameter and repeated this procedure. Then, we used the median percentage variation in relative abundance as a reference for each species and compared it to the STOC.
Écrins National Park with the relative abundance trends of common birds in the STOC. If we denote $\lambda_i^{(STOC)}(y)$ as the estimated abundance of species i in year y by the STOC, then the percentage variation in abundance computed by the STOC (Fontaine et al., 2020) can be expressed as:

$$ (\hat{\lambda}_i^{(STOC)}(2019) - \hat{\lambda}_i^{(STOC)}(2001)) \times \frac{100}{\hat{\lambda}_i^{(STOC)}(2001)} $$

In the model presented in Equation 6, the term capturing the effect of the period p on the target-group relative abundance of species i is $exp(f_i(p))$. The INLA method allows us to estimate a posterior distribution of each $exp(f_i(p))$. To mimic the STOC period of study, we assume the target-group relative abundance to be constant in each period p. This allows us to attribute an estimated relative effect $f_i(y)$ to each year y between 2001 and 2019. We then calculated a regression line using a linear model with the $f_i(y)$ as the response variable and y the explanatory variable. If we denote $\hat{f}_i^{(ENP)}(y)$ the fitted value of this linear regression, we could calculate the percentage variations in target-group relative abundance of species i as

$$ (\hat{f}_i^{(ENP)}(2019) - \hat{f}_i^{(ENP)}(2001)) \times \frac{100}{\hat{f}_i^{(ENP)}(2001)} $$

To enhance robustness and take into account statistical uncertainty, we generated 1000 posterior distributions for each target-group relative abundance trend parameter and repeated this procedure. Then, we used the median percentage variation in target-group relative abundance as a reference for each species and compared it to the STOC.

348

# Results

349

## Validation

Figure 1 displays the boxplots of AUC values computed species by species as described in section Validation and post-processing. We notice a clear improvement in AUC when using the full linear predictor compared to using only the sampling-effort proxy as offset. This confirms that opportunistic data are not solely constrained by sampling methods and that valuable information can be extracted from them.

Regarding the boxplots of AUC values obtained using cross-validation (with one value for each hold-out dataset), we observe a few outliers among the species, some of which occur in all scenarios, as it is the case for *Cinclus cinclus* and *Motacilla cinerea*. These two species with poor AUC values relative to others are stream specialist species, for which covariates do not capture the habitat well. It is worth noting that a species appearing multiple times as an outlier within the same scenario is listed only once for readability.

The AUC values of models with cross-validation are typically slightly lower than those without cross-validation, as expected due to the reduction in the amount of training data and the more challenging prediction setting with new data not used during training. However, they remain generally stable, suggesting that the model does not overfit. Furthermore, despite this relative decrease, the AUC values generally remain high, with approximately 75% of AUC values above 0.8 across all scenarios. Consequently, we conclude that the model demonstrates a very good fit to the data, except for two species: *Cinclus cinclus* and *Motacilla cinerea*, which we have excluded from further analysis.

361
362
363
364
365
366
367
368
# Results

330

331

## Validation

332

333

334

335

336

Figure 1 displays the boxplots of AUC values computed species-by-species as described in section *Validation and post-processing*. There is a clear improvement in AUC when using the full linear predictor compared to using only the sampling-effort proxy as offset. This confirms that sampling methods do not solely constrain opportunistic presence-only data and that valuable information can be extracted from them.

337

338

339

340

341

342

343

344

345

The AUC values of models with cross-validation (with one value for each hold-out dataset) are typically slightly lower than those without cross-validation. This is expected due to the reduction in training data and the more challenging prediction setting with new data not used during training. Despite this slight decrease, the AUC values generally stay high, with approximately 75% of values above 0.8 across all scenarios. Consequently, the model demonstrates a very good fit to the data, except for two species: *Cinclus cinclus* and *Motacilla cinerea*, which are consistent outliers in all scenarios. These two species with poor AUC values relative to others are stream specialist species, for which covariates do not capture the habitat well so we have excluded them from further analysis.

346

## Fixed Effects

347

348

349

350

351

352

353

354

355

356

357

358

359

The mean estimated effects of each covariate and their significance are detailed in Tables 5-7 in Appendix F (Lasgorceux et al., 2024c). To summarize the information on habitat spatial covariates, we provide a classification into three clusters based on Principal Component Analysis conducted on the estimated mean coefficients for each habitat spatial covariate (see subsection 'Materials and Methods - Validation and Post-processing - Fixed Effects'). We did not include any species with less than two significant covariates since we consider that the model does not provide enough information on them. These species are generally those with very few occurrences, as 10 out of these 11 species have fewer than 551 occurrences over the 28 years of data, where the value 551 is below the 25th-percentile of species occurrences. The list of removed species is given in Table 8 in Appendix F (Lasgorceux et al., 2024c). The PCA reveals two principal axes explaining 62.2% of the classification: the first axis distinguishes vegetated habitat (coniferous and deciduous forests VS urban wetlands), while the second axis highlights an altitudinal gradient.

360

361

362

363

364

365

The three clusters are presented in Figure 2. We interpret them as follows. Cluster 1: Species mostly found in closed habitats such as coniferous and deciduous forests (green); Cluster 2: Species mostly found in high-elevation environments (orange); Cluster 3: Species mostly found in open habitats and valleys (violet). These results are coherent when comparing a map of elevation with a map of clusters with the highest occurrences on a 500m×500m regular grid in the Écrins National Park (see Figure 15 in Appendix F in Lasgorceux et al., 2024c).

366

367

368

369

370

371

372

The ellipse associated with forest species represents an altitudinal gradient, with species found in low-elevation forests (e.g. *Pyrrhula pyrrhula*), to medium (e.g. *Anthus trivialis*) and high-elevation forests (e.g. *Turdus torquatus*). We also observe this altitudinal gradient in the cluster associated with species from open habitats and valleys, with rocky habitat species such as *Tichodroma muraria* and *Phoenicurus ochruros*, at the boundary of the high-elevation species cluster. The cluster associated with species specialized to high-elevation habitats is smaller but well-defined, with typical species such as *Prunella collaris* and *Montifringilla nivalis*.
<figure><img src="image_2.png" /></figure>

**Figure 1** – AUC values calculated based on various linear-predictor and cross-validation scenarios. Boxplots summarize the set of AUC values for all species, and in the case of cross-validation for each combination of species and hold-out dataset. Notations are as follows: “Offset” – model where linear predictor contains only the offset; “Linearpredictor” – model with full linear predictor; “K-folds” – model fitted by randomly removing 10% of the cells, ten times; “Concatenation of K-folds” – AUC values computed by concatenating the predictions of the K-folds scenarios; “Temporal” – model fitted by randomly removing half of the data in a given time period, for each time period; “Spatial” – model fitted by removing all the data which occurred in a given city, for seven cities. The species with multiple outlier AUC values in “K-folds”, “Temporal” and “Spatial” scenarios are only noted once. Details on the four cross-validation scenarios (where models always have the full linear predictor) can be found in the “Post-processing and Validation - Validation” section.

369

Fixed Effects

370

371

372

373

374

375

376

The mean estimated effects of each covariate and their significance are detailed in Tables 5-7 in Appendix 6. To summarize the information on habitat spatial covariates, we provide a classification into three clusters based on Principal Component Analysis conducted on the estimated mean coefficients for each habitat spatial covariate (see subsections “Data - Covariates” and “Validation and Post-processing - Fixed Effects”). We did not include any species with less than two significant covariates since we consider that the model does not provide enough information on them. These species are generally those with very few occurrences, as 10 out of these
<figure><img src="image_2.png" /></figure>

**Figure 1** – AUC values calculated for various linear-predictor and cross-validation scenarios. Boxplots summarize the set of AUC values for all species, and in the case of cross-validation for each combination of species and hold-out dataset. Notations are as follows: “Offset” – model where linear predictor contains only the offset; “Linear predictor” – model with full linear predictor; “K-folds” – model fitted by randomly removing 10% of the cells, ten times; “Concatenation of K-folds” – AUC values computed by concatenating the predictions of the K-folds scenarios; “Temporal” – model fitted by randomly removing half of the data in a given time period, for each time period; “Spatial” – model fitted by removing all the data which occurred in a given city, for seven cities. The species with multiple outlier AUC values in “K-folds”, “Temporal” and “Spatial” scenarios are only noted once. Details on the four cross-validation scenarios (where models always have the full linear predictor) can be found in the ‘Materials and Methods – Post-processing and Validation – Validation’ section.

To compare the contribution of spatio-temporal covariates with spatial ones in the linear pre-
dictor, we computed the percentage of explained variance for each covariate with respect to
the overall variance of the linear part $(\alpha_i + \mathbf{X}^{(PCA)}(s)\beta_i + \mathbf{X}^{(anomalies)}(s, m, p)\gamma_i)$ for each species.
The results are presented in Figure 3. Our analysis revealed that, for most species, PCA Axis 1,
interpreted as being related to high-elevation environments, explains the majority of the vari-
ance (Median=54%). Axis 2 and 3, interpreted as environments with medium to high altitude
vegetation and grasslands, respectively, explain a substantial portion of the variance for about
twenty species (Q3=27% and 15%).

Regarding spatio-temporal covariates, precipitation anomalies had the smallest percentage
of explained variance, with a median below 10⁻⁵%. This is further supported by only 12 out of 76
11 species have fewer than 551 occurrences over the 28 years of data, where the value 551 is below the 25th-percentile of species occurrences. The list of removed species is given in Table 8 in Appendix 6. The PCA reveals two principal axes explaining 62.2% of the classification: the first axis distinguishes vegetated habitat (coniferous and deciduous forests VS urban wetlands), while the second axis highlights an altitudinal gradient.

The three clusters are presented in Figure 2. We interpret them as follows. Cluster 1: Species mostly found in closed habitats such as coniferous and deciduous forests (green); Cluster 2: Species mostly found in high-elevation environments (orange); Cluster 3: Species mostly found in open habitats and valleys (violet). These results are coherent when comparing a map of elevation with a map of clusters with the highest occurrences on a 500m×500m regular grid in the Écrins National Park (see Figure 14 in Appendix 6).

The ellipse associated with forest species represents an altitudinal gradient, with species found in low-elevation forests (*Pyrrhula pyrrhula*), to medium (*Anthus trivialis*) and high-elevation forests (*Turdus torquatus*). We also observe this altitudinal gradient in the cluster associated with species from open habitats and valleys, with rocky habitat species such as *Tichodroma muraria* and *Phoenicurus ochruros*, at the boundary of the high-altitude species cluster. The cluster associated with species specialized to high-elevation habitats is smaller but well-defined, with typical species such as *Prunella collaris* and *Montifringilla nivalis*.

To address the spatio-temporal covariates, we began by examining their impact on the linear predictor across all species. We computed the standard deviation of each covariate $X_j \in \mathbf{X}^{(PCA)}$ (and $\mathbf{X}^{(anomalies)}$), multiplied by the mean estimated associated parameter $\beta_{i,j}$ (or $\gamma_{i,j}$) for each species, and then averaged these values. Our analysis revealed that the precipitation anomalies covariate had the smallest impact on the linear predictor, with a Mean(SD) of 0.09. This finding was further supported by only 12 species exhibiting a significant parameter associated with precipitation anomalies. In contrast, the temperature anomalies covariate had a substantial effect on the linear predictor, with a Mean(SD) of 0.52, which is significant when compared to the elevation PCA axis (Mean(SD) = 0.79). Notably, 45 species demonstrated a significant parameter associated with maximum temperature anomalies. Among these, we identified species with positive parameters—indicating they tend to be more prevalent when deviations in maximum temperatures from the mean are positive—such as *Emberiza hortulana* and *Lanius collurio*. Conversely, we also identified species with negative parameters—indicating they are more frequently observed when deviations in maximum temperatures from the mean are negative—such as *Coccothraustes coccothraustes* and *Fringilla montifringilla*. These results align with the ecological characteristics of these species, confirming that thermophilic species thrive in warmer conditions, while cold-resistant species are more abundant in cooler conditions.

## Intra-annual effect

The first two extracted harmonics encompass a significant portion of the intra-annual variability, explaining 83% of the variance. While the primary harmonic exhibits a greater positive effect during winter months compared to summer months, the second harmonic shows a positive effect in summer months and a negative effect in winter months (see Figure 15 in Appendix 7). Consequently, species with positive scores on the first harmonic tend to be overrepresented in winter, whereas species with positive scores on the second harmonic are more likely to be overrepresented in spring and summer, which hints at a potential link with the migratory status of species.
<figure><img src="image_3.png" /></figure>

**Figure 2** – Classification of 63 species in three clusters based on the estimated mean effects of spatial covariates. A Principal Component Analysis was conducted on the species with at least two credibility intervals that do not include 0. The two main axes explain 62.2% of the variance. Cluster 1: Species mostly found in closed habitats such as coniferous and deciduous forests (green); Cluster 2: Species mostly found in high-elevation environments (orange); Cluster 3: Species mostly found in open habitats and valleys (valley).

species showing a significant parameter associated with precipitation anomalies. Similarly, the explained variance for maximum temperature anomalies is low for most species (Median < 10⁻³%); however, 45 species showed a significant parameter. It notably explained a significant amount of information for five species: *Coccothraustes coccothraustes*, *Fringilla montifringilla*, *Lanius colur* (*lurio*, *Spinus spinus*, and *Turdus pilaris*. These results align with the ecological characteristics of these species, confirming that the model effectively correlates temperature anomalies with the presence of thermophilic and cold-resistant species.

## Intra-annual effect

The first two extracted harmonics capture a significant portion of the intra-annual variability, explaining 83% of the variance. The primary harmonic shows a stronger positive effect in winter months compared to summer months, while the second harmonic has a positive effect in spring and summer and a negative effect in winter (see Figure 16 in Appendix G in Lasgorceux et al., 2024c). Thus, species with positive scores on the first harmonic, referred to as the winter harmonic, are predominantly overrepresented in winter. Conversely, species with positive scores on the second harmonic, now termed the warm months harmonic, are more likely to be
<figure><img src="image_3.png" /></figure>

**Figure 2** – Classification of 63 species in three clusters based on the estimated mean effects of spatial covariates. A Principal Component Analysis was conducted on the species with at least two credibility intervals that do not include 0. The two main axes explain 62.2% of the variance. Cluster 1: Species mostly found in closed habitats such as coniferous and deciduous forests (green); Cluster 2: Species mostly found in high-elevation environments (orange); Cluster 3: Species mostly found in open habitats and valleys (valley).

Figure 3 illustrates the scores associated with the two first harmonic for each species, categorized by migratory status within the territory of the Écrins National Park, as provided by experts from the park. Colored ellipses are provided for readability. This analysis confirms that species overrepresented in winter and underrepresented in spring and summer primarily consist of sedentary species, while those underrepresented in winter and overrepresented in spring and summer correspond to migratory species. The details of the scores for each species are given in Table 9 in Appendix 7. These harmonics effectively distinguish between sedentary and migratory species within the Écrins National Park, affirming that the month-based spatio-temporal Gaussian field $W_i^{(m)}$ adequately accounts for spatial and intra-annual species variations not explained by the physical covariates across the territory.
<figure><img src="image_4.png" /><figcaption>**Figure 3** – Percentage of explained variance of each covariate for 74 species. The co-variates include six spatial factors derived from PCA analysis and two spatio-temporal climatic factors. The interpretation of the PCA axes is detailed in the 'Materials and Methods - Data - Preprocessing of model covariates'. Covariate boxplots are ordered by their median values.</figcaption></figure>

overrepresented in spring and summer, suggesting a potential link with migratory status.

Figure 4 displays the harmonic scores for each species, categorized by migratory status within the territory of the Écrins National Park, as provided by experts from the park. This analysis confirms that species overrepresented in winter are mainly sedentary, while those overrepresented in warm months are migratory. Detailed scores for each species are provided in Table 9 in Appendix G (Lasgorceux et al., 2024c). These harmonics effectively distinguish between sedentary and migratory species in the Écrins National Park, affirming that the month-based spatio-temporal Gaussian field $W_i^{(m)}$ accurately accounts for spatial and intra-annual species variations not explained by physical covariates across the territory.

Inter-annual variability and STOC comparison

The estimated percentage variations in relative abundance in the Écrins National Park between 2001 and 2019, based on 1000 *a posteriori* samples, are illustrated in Figure 5. Species associated with the high-elevation habitats cluster (Figure 2) exhibit a mean decline of –16.5% in relative abundance during 2001–2019. This is mainly due to the strong decrease of three species:
<figure><img src="image_4.png" /></figure>

**Figure 3** – Point plot of scores associated with the first two harmonics of the FPCA computed on the month-based spatio-temporal Gaussian field $W_i^{(m)}$ estimated for each species. The two first harmonics explain 83% of the variance. The species are categorized by their main migratory status within the territory of the Écrins National Park (ENP), as provided by experts from the park: breeding migratory (red); resident (green); wintering (blue).

431 **Inter-annual variability and STOC comparison**

432 The estimated percentage variations in target-group relative abundance in the Écrins National Park between 2001 and 2019, based on 1000 *a posteriori* samples, are illustrated in Figure 4. 433 Species associated with high-elevation habitats exhibit a mean decline of -16.5% in target- 434 group relative abundance during 2001–2019. This is mainly due to the strong decrease of three 435 species: *Monticola saxatilis* (-53.5%), *Prunella collaris* (-26.4%) and *Montifringilla nivalis* (-23.4%). 436 Conversely, species affiliated with the forest cluster experienced a growth of 12.4% over the 437 same period. 438

439 Figure 4 shows that, out of the 19 forest species, 6 exhibited a positive trend while only 2 440 showed a negative trend based on the 95% credibility interval. For comparison, it is noteworthy 441 that the STOC reported a decrease of 9.7% in abundance for specialist forest species in France. 442 When focusing specifically on passerines considered forest specialists by the STOC, the differ- 443 ence becomes even more pronounced, with a growth of 23.2% in the Écrins National Park, while 444 the STOC observed a growth of 2% for these species during the same period in France. It is im- 445 portant to note that this percentage is nevertheless skewed by *Regulus ignicapilla*, which shows 446 a growth of 168% in relative abundance over the period. 448
<figure><img src="image_5.png" /><figcaption>**Figure 4** – Point plot of scores associated with the first two harmonics of the FPCA computed on the month-based spatio-temporal Gaussian field $W_i^{(m)}$ estimated for each species. The two first harmonics explain 83% of the variance. The species are categorized by their main migratory status within the territory of the Écrins National Park (ENP), as provided by experts from the park: breeding migratory (red); resident (green); wintering (blue).</figcaption></figure>

Monticola saxatilis (–53.5%), Prunella collaris (–26.4%) and Montifringilla nivalis (–23.4%). Conversely, species affiliated with the forest cluster experienced a growth of 12.4% over the same period. Figure 5 shows that, out of the 19 forest species, 6 exhibited a positive trend while only 2 showed a negative trend based on the 95% credibility interval. For comparison, it is noteworthy that the STOC reported a decrease of 9.7% in abundance for specialist forest species in France. When focusing specifically on the passerines considered as forest specialists by the STOC, the difference becomes even more pronounced, with a growth of 23.2% in the Écrins National Park, while the STOC observed a growth of 2% for these species during the same period in France. It is important to note that this percentage is nevertheless skewed by Regulus ignicapilla, which shows a growth of 168% in relative abundance over the period.

Species linked to the valleys cluster experienced a mean decrease of -3.2% in relative abundance. However, comparing this figure with STOC trends is challenging due to the mix of generalist species with specialists of built and agricultural environments within this cluster. To address this, we compare below the trends of passerines in the Écrins National Park, classified by the STOC as generalist species or species specialized in built and agricultural environments, with our estimated trends in relative abundance. Built-environment passerines demonstrate stability
Species linked to the valleys cluster experienced a mean decrease of -3.2% in target-group relative abundance. However, comparing this figure with STOC trends is challenging due to the mix of generalist species with specialists of built and agricultural environments within this cluster. To address this, we compare below the trends of passerines in the Écrins National Park, classified by the STOC as generalist species or species specialized in built and agricultural environments, with our estimated trends in target-group relative abundance. Built-environment passerines demonstrate stability in the Écrins National Park (0.26%) compared to a slight decrease in France (-5.68%). Generalist passerines show an increase of 7.4% in target-group relative abundance in the park, aligning with the 7% growth reported by the STOC for these species (and 19.4% for generalist common birds, not only passerines). For farmland species, both show a decrease, but at a slower rate of -16.5% compared to -22.6% (and 29.5% for farmland common birds, not only passerines).

<figure><img src="image_5.png" /><figcaption>**Figure 4** – Illustration of the estimated median percentage variations in target-group relative abundance and credibility intervals for 74 retained species in the Écrins National Park spanning the period 2001–2019. These trends are derived from 1000 samples of the *a posteriori* distribution of the inter-annual effect $exp(f_i(p))$. Details can be found in the “Post-processing and Validation - Inter-annual effects” section. Credibility intervals shown in red (respectively blue) correspond to intervals below (respectively above) zero, while yellow intervals include zero. Species names are color-coded according to their cluster affiliation in Figure 2. Species with fewer than two significant covariates were not classified and are shown in black.</figcaption></figure>
in the Écrins National Park (0.26%) compared to a slight decrease in France (-5.68%). Generalist passerines show an increase of 7.4% in relative abundance in the park, aligning with the 7% growth reported by the STOC for these species (and 19.4% for generalist common birds, not only passerines). For farmland species, both show a decrease, but at a slower rate of -16.5% compared to -22.6% (and -29.5% for farmland common birds, not only passerines).

<figure><img src="image_6.png" /></figure>

**Figure 5** – Illustration of the estimated median percentage variations in relative abundance and credibility intervals for 74 retained species in the Écrins National Park spanning the period 2001–2019. These trends are derived from 1000 samples of the *a posteri*
*ori* distribution of the inter-annual effect $exp(f_i(p))$. Details can be found in the ‘Materials and Methods – Post-processing and Validation – Inter-annual effects’ section. Credibility intervals shown in red (respectively blue) correspond to intervals below (respectively above) zero, while yellow intervals include zero. Species names are color-coded according to their cluster affiliation in Figure 2. Species with fewer than two significant covariates were not classified and are shown in black.

Figure 6 displays more precisely the percentage variations of relative abundance in the Écrins National Park against the percentage variations of abundance in France computed by the STOC, between 2001 and 2019. More details on the trends are given in Table 10, Appendix H (Lasgorceux et al., 2024c). Each species lies in one of the four quadrants. Quadrants (ii) and (iii) correspond to species whose trends locally within the Écrins National Park align with those in France: increasing for quadrant (ii) and decreasing for quadrant (iii). Quadrant (i) corresponds to species trending upward nationally but downward in the Écrins National Park, whereas quadrant (iv) represents species increasing in the Écrins National Park but decreasing nationally. The number of species in each quadrant is as follows: 8 for (i), 14 for (ii), 24 for (iii) and 13 for (iv).
Figure 5 displays more precisely the percentage variations of target-group relative abundance in the Écrins National Park against the percentage variations of abundance in France computed by the STOC, between 2001 and 2019. More details on the trends are given in Table 10, Appendix 8. Each species lies in one of the four quadrants. Quadrants (ii) and (iii) correspond to species whose trends locally within the Écrins National Park align with those in France: increasing for quadrant (ii) and decreasing for quadrant (iii). Quadrant (i) corresponds to species trending upward nationally but downward in the Écrins, whereas quadrant (iv) represents species increasing in the Écrins but decreasing nationally. The number of species in each quadrant is as follows: 8 for (i), 14 for (ii), 24 for (iii) and 14 for (iv). Quadrant (iv) is primarily composed of forest species, while quadrant (i) consists mainly of generalist species, particularly common ones.

Furthermore, while the national park status does not prevent the decline of certain species, it is noteworthy that for 38 species out of 60, trends are more favorable within the national park compared to national trends for all of France as represented by the y = x slope in Figure 5. However, some species for which this is not the case are very common (quadrant (i)) and may be subject to reporting bias (see Discussion).

<figure><img src="image_6.png" /><figcaption>**Figure 5** – Point plot of the estimated median percentage variations of target-group relative abundance in the Écrins National Park between 2001 and 2019, against the percentage variations of relative abundance in France between 2001 and 2019 computed by the STOC for 60 species. The names of species with less than 0.1% variation in target-group relative abundance are not plotted for readability.</figcaption></figure>
443
444

Quadrant (iv) is primarily composed of forest species, while quadrant (i) consists mainly of generalist species, particularly common ones.

445
446
447
448
449

Furthermore, while the national park status does not prevent the decline of certain species, it is noteworthy that for 38 species out of 59, trends are more favorable within the national park compared to national trends for all of France. The names of these species are represented under the diagonal in Figure 6. However, some species for which this is not the case are very common (quadrant (i)) and may be subject to reporting bias (see Discussion).

<figure><img src="image_7.png" /></figure>

**Figure 6** – Point plot of the estimated median percentage variations of relative abundance in the Écrins National Park between 2001 and 2019, against the percentage variations of abundance in France between 2001 and 2019 computed by the STOC for 59 species. The names of species with less than 0.1% variation in relative abundance are not plotted for readability.

450

## Discussion

451
452
453
454
455
456
457
458

In our study, we demonstrated that under the assumption of sampling homogeneity across species (but where sampling effort is still allowed to vary in space and time), opportunistic presence-only data could be leveraged to estimate habitat preferences, migratory status, and relative abundance trends of passerine species. This finding is not novel *per se*, as numerous studies have utilized such data to model species distributions or abundance trends (Botella et al., 2020; Bradter et al., 2018; Kéry et al., 2010; Phillips et al., 2009; Valavi et al., 2022; Van Strien et al., 2013). However, it is one of only a few approaches achieving this within a unique spatio-temporal model. To our knowledge, the most similar work has been done by
Discussion

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

In our study, we demonstrated that under the assumption of sampling homogeneity across species (but where sampling effort is still allowed to vary in space and time), opportunistic presence-only data could be leveraged to estimate habitat preferences, migratory status, and target-group relative abundance trends of passerine species. This finding is not novel *per se*, as numerous studies have utilized such data to model species distributions or abundance trends (Botella et al., 2020; Bradter et al., 2018; Kéry et al., 2010; Phillips et al., 2009; Valavi et al., 2022; Van Strien et al., 2013). However, it is one of only a few approaches achieving this within a unique spatio-temporal model. To our knowledge, the most similar work has been done by Giraud et al., 2016 and Coron et al., 2018 modeling spatial relative abundance (abundance divided by the abundance of an arbitrarily chosen cell) of passerines using a Poisson distribution. Their work provides theoretical results regarding parameter identifiability, as well as scenarios for data simulation and parameter estimation. However, the use of their models requires access to appropriate survey data based on a standardized sampling protocol, whereas our study relied solely on a single set of opportunistic data. We also found the inclusion of the month-based spatio-temporal latent Gaussian field $W_i^{(m)}$ and the use of Functional Principal Component Analysis to be pertinent for assessing phenological effects when modeling target-group relative abundances of species with different migratory statuses. Another potential application of such Gaussian fields and FPCA could involve implementing a nonlinear additive predictor component for altitude to study species' phenological changes as responses to climate: altitude ascents, migrations, etc.

The STOC has demonstrated the effectiveness of public policies on protected areas. Outside of reserves, the average decline of 56 species between 2004 and 2018 was -6.6%, whereas within reserves, an increase of 12.5% was observed (Fontaine et al., 2020). Additionally, our study offers new insights into the role of protected areas in mitigating the decline of passerine populations, though our trends should be interpreted with caution. Indeed, several very common species show a negative estimated target-group relative abundance in the park, while the STOC indicates positive trends for these same species. Thus, we believe that our assumption of homogeneous reporting is overly optimistic and partially inaccurate. Despite this issue affecting 7 out of 60 species, we can still draw meaningful conclusions.

Notably, we demonstrate that forest passerines benefit the most in the Écrins National Park. While not reversing the trend, the decline of farmland passerines in the park is slower than in the rest of France, consistent with findings from previous studies (Palacín and Alonso, 2018; Silva et al., 2018). Furthermore, our study provides novel insights into high mountain species not previously investigated by the STOC. Our indicator, the target-group relative abundance, shows a global decline for species specialist to high-elevation environments. This trend aligns with the global decrease of mountain specialist species in Europe (-10% during 2002–2014) as reported by Lehikoinen et al., 2019. These findings will be valuable for comparison with upcoming analyses of the STOM, a survey analogous to the STOC for mountain species, which has only been operational since 2012, with its first results yet to be analyzed.

We identify three areas for future improvements in our modeling approach using solely opportunistic data:

(i) The assumption of uniform sampling, which could be questionable due to differences in detectability, species commonness/rarity (Snäll et al., 2011), or agent identification skills and reporting prevalence (Bradter et al., 2018), albeit reasonable in practice and lacking better options.
459

460

461

462

463

464

465

466

467

468

469

470

Giraud et al., 2016 and Coron et al., 2018 modeling spatial relative abundance (abundance divided by the abundance of an arbitrarily chosen site) of passerines using a Poisson distribution. Their work provides theoretical results regarding parameter identifiability, as well as scenarios for data simulation and parameter estimation. However, the use of their models requires access to appropriate survey data based on a standardized sampling protocol, whereas our study relied solely on a single set of opportunistic data. We also found the inclusion of the month-based spatio-temporal latent Gaussian field $W_i^{(m)}$ and the use of Functional Principal Component Analysis to be pertinent for assessing phenological effects when modeling relative abundances of species with different migratory statuses. Another potential application of such Gaussian fields and FPCA could involve implementing a nonlinear additive predictor component for altitude to study species' phenological changes as responses to climate: altitude ascents, migrations, etc.

471

472

473

474

475

476

477

478

479

The STOC has demonstrated the effectiveness of public policies on protected areas. Outside of reserves, the average decline of 56 species between 2004 and 2018 was −6.6%, whereas within reserves, an increase of 12.5% was observed (Fontaine et al., 2020). Additionally, our study offers new insights into the role of protected areas in mitigating the decline of passerine populations, though our trends should be interpreted with caution. Indeed, several very common species show a negative estimated relative abundance in the park, while the STOC indicates positive trends for these same species. Thus, we believe that our assumption of homogeneous reporting is overly optimistic and partially inaccurate. Despite this issue affecting 8 out of 59 species, we can still draw meaningful conclusions.

480

481

482

483

484

485

486

487

488

489

Notably, we demonstrate that forest passerines benefit the most in the Écrins National Park. While not reversing the trend, the decline of farmland passerines in the park is slower than in the rest of France, consistent with findings from previous studies (Palacín and Alonso, 2018; Silva et al., 2018). Furthermore, our study provides novel insights into high mountain species not previously investigated by the STOC. Our indicator, the relative abundance, shows a global decline for species specialist to high-elevation environments. This trend aligns with the global decrease of mountain specialist species in Europe (−10% during 2002–2014) as reported by Lehikoinen et al., 2019. These findings will be valuable for comparison with upcoming analyses of the STOM, a survey analogous to the STOC for mountain species, which has only been operational since 2012, with its first results yet to be analyzed.

490

491

We identify three areas for future improvements in our modeling approach using solely opportunistic data:

492

493

494

495

496

497

498

499

500

501

502

(i) The assumption of homogeneous sampling of different species, which could be questionable due to differences in detectability, species commonness/rarity (Snäll et al., 2011), or agent identification skills and reporting prevalence (Bradter et al., 2018), albeit reasonable in practice and lacking better options. Heterogeneous detectability remains an issue without data based on standardized sampling protocols and seems difficult to correct except for trying to include prior expert knowledge into models. Different reporting behavior due to species commonness/rarity can introduce bias in estimated intercept terms, such as $\alpha_i$ in Equation 6, but would not have any impact on the other terms varying in space and/or time, as long as this bias in reporting behavior is homogeneous across space, time and physical predictors. Effects related to agent identification skills and reporting prevalence could potentially be identified by including random effects for the rangers in the model. As a general recommendation to reduce biases and improve
Heterogeneous detectability remains an issue without data based on standardized sampling protocols and seems difficult to correct except for trying to include prior expert knowledge into models. Different reporting behavior due to species commonness/rarity can introduce bias in estimated intercept terms, such as $\alpha_i$ in Equation 6, but would not have any impact on the other terms varying in space and/or time, as long as this bias in reporting behavior is homogeneous across space, time and physical predictors. Effects related to agent identification skills and reporting prevalence could potentially be identified by including random effects for the rangers in the model. As a general recommendation to reduce biases and improve the exploitability of opportunistic presence-only data, we suggest that when observers report a species, all species of the same taxon should also be reported for better use of opportunistic presence-only data in species distribution modeling.

(ii) The acquired information is relative to a target-group of species since we used the target-group occurrences as a proxy for the sampling effort (Equation 4). This implies that the potential for interpretations of results in terms of absolute abundance and its variation across covariates, space, and time remains relatively limited.

(iii) There are multiple ways to define this target-group, which impacts the interpretation of abundances relative to it. The composition of this target-group primarily relies on available data sources, ranging from studies covering a single species (Farr et al., 2021) to those encompassing thousands (Botella et al., 2021). Some studies focus exclusively on species within the same taxonomic group (Van Strien et al., 2013), while others consider species across different taxa (Escamilla Molgora et al., 2022). Botella et al., 2020 recommends choosing a set of species that are consistently abundant across a broad range of environmental sub-regions, that we did in this study.

We have recognized a pressing need for more comprehensive descriptions of what we refer to as “sampling effort” when dealing with opportunistic presence-only data, and for methodologies to effectively incorporate it into modeling endeavors. Various strategies have been proposed to address sampling efforts depending on the studied taxa, the available datasets and the type of model. These include the utilization of random or target-group pseudo-absences with Maxent (Phillips et al., 2006, 2009), target-group background occurrences (Botella et al., 2020), integration of information from detection/nondetection data (Dorazio, 2014; Fithian et al., 2015; Giraud et al., 2016), or questionnaires (Bradter et al., 2018). Additionally, consideration of proxy variables such as road density or population density (Kervellec et al., 2023), and covariates representing the accessibility of a site, such as distance to roads, cities, coastlines, and protected areas (El-Gabbas and Dormann, 2018; Henckel et al., 2020; Moreira et al., 2024; Warton et al., 2013) or refinements between these methods and the target-group approach (Chauvier et al., 2021). Another commonly used proxy, especially in studies focusing on cetaceans and fish, is the logbooks of vessels (Alglave et al., 2022; Ver Hoef et al., 2021).

Our findings suggest that, despite the statistical challenges arising from taking sampling biases into account, opportunistic presence-only data provide unparalleled spatio-temporal coverage, allowing for the depiction of consistent biological indicators that would be challenging to obtain using traditional methods. This includes monitoring trends in target-group relative abundance over extended periods and developing species distribution models across vast geographical areas. Such information can serve as valuable quantitative indicators for assessing the effectiveness of conservation policies in protected areas.
the exploitability of opportunistic presence-only data, we suggest that when observers report a species, all species of the same taxon should also be reported for better use of opportunistic presence-only data in species distribution modeling.

(ii) The acquired information is relative to a target group of species since we used the target-group occurrences as a proxy for the sampling effort (Equation 4). This implies that the potential for interpretations of results in terms of absolute abundance and its variation across covariates, space, and time remains relatively limited, as long as we do not have additional external information on the absolute abundance of the target group as a whole.

(iii) There are multiple ways to define this target group, which impacts the interpretation of abundances relative to it. The composition of this target group primarily relies on available data sources, ranging from studies covering a single species (Farr et al., 2021) to those encompassing thousands (Botella et al., 2021). Some studies focus exclusively on species within the same taxonomic group (Van Strien et al., 2013), while others consider species across different taxa (Escamilla Molgora et al., 2022). Botella et al., 2020 recommends choosing a set of species that are consistently abundant across a broad range of environmental sub-regions, which is what we did in this study.

We have recognized a pressing need for more comprehensive descriptions of what we refer to as “sampling effort” and “observer bias” when dealing with opportunistic presence-only data, as well as for methodologies to effectively incorporate it into modeling endeavors. Various strategies have been proposed to address sampling efforts and observer bias depending on the studied taxa, the available datasets, and model type. These strategies include the use of random or target-group pseudo-absences with Maxent (Phillips et al., 2006, 2009), proxy variables representing site accessibility (Henckel et al., 2020; Moreira et al., 2024; Warton et al., 2013) and refinements (Chauvier et al., 2021), integration of information from detection/nondetection data (Dorazio, 2014; Fithian et al., 2015; Giraud et al., 2016), and questionnaires (Bradter et al., 2018). In our case study, very large occurrence numbers were available for the target group, such that uncertainties about the cumulative abundance of the target group can be expected to be relatively small. For modeling more precisely the cumulative abundance of the target group, for example in cases with smaller occurrence numbers or when large areas of the study domain are only very weakly sampled, we could develop an alternative approach wherein $Y_{TG}(s, m, p)$ could be initially modeled with a Poisson distribution using a mean parameter $\mu_{TG}(s, m, p)$, fitted to the data to estimate $\mu_{TG}(s, m, p)$. Subsequently, we would use the estimated Poisson mean as a proxy for $\Lambda_{TG}(s, m, p)E(s, m, p)$, rather than the raw count. This modification could potentially enhance robustness and improve the assessment of uncertainties inherent to our modeling approach.

Our findings suggest that, despite the statistical challenges arising from taking sampling biases into account, opportunistic presence-only data provide unparalleled spatio-temporal coverage, allowing for the depiction of consistent biological indicators that would be challenging to obtain using traditional methods. This includes monitoring trends in relative abundance over extended periods and developing species distribution models across vast geographical areas. Such information can serve as valuable quantitative indicators for assessing the effectiveness of conservation policies in protected areas.

545
In our case study, very large occurrence numbers were available for the target-group, such that uncertainties about the target-group abundance can be expected to be relatively small. For modeling more precisely the target-group abundance, for example in cases with smaller occurrence numbers or when large areas of the study domain are only very weakly sampled, we could develop an alternative approach wherein $Y_{TG}(s, m, p)$ could be initially modeled with a Poisson distribution using a mean parameter $\mu_{TG}(s, m, p)$, fitted to the data to estimate $\mu_{TG}(s, m, p)$. Subsequently, we would use the estimated Poisson mean as a proxy for $\Lambda_{TG}(s, m, p)E(s, m, p)$, rather than the raw count. This modification could potentially enhance robustness and improve the assessment of uncertainties inherent to our modeling approach.

## Acknowledgements

The authors would like to thank all the staff of the Écrins National Park and its scientific team for contributing to the dataset. They also appreciate advice of Denis Allard on using factorial principal component analysis, and Loïc Houde for technical support with the high-performance computing cluster used to implement the models of this case study.

## Fundings

The authors declare that they have received no specific funding for this study.

## Conflict of interest disclosure

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.

## Data, script, code, and supplementary information availability

Data are available online (https://doi.org/10.24072/fake1[Replace by the DOI of the webpage hosting the data]; cite[Replace by the citation of the data])
Script and codes are available online (https://doi.org/10.24072/fake2[Replace by the DOI of the webpage hosting the script and code]; cite[Replace by the citation of the script and code])
Supplementary information is available online (https://doi.org/10.24072/fake3[Replace by the DOI of the webpage hosting the Supplementary information]; cite[Replace by the citation of the Supplementary information])

The DOI hyperlinks should be active. They should also be present in the reference list and cited in the text. For the reference section below, do not forget to add a doi for each reference (if available). Do not forget to add the reference of the recommendation, the reference of the data, scripts, code and supplementary material to your bib file, if appropriate.

## References

Alglave B, Rivot E, Etienne MP, Woillez M, Thorson JT, Vermard Y (2022). Combining scientific survey and commercial catch data to map fish distribution. *ICES Journal of Marine Science* **79**, 1133–1149. https://doi.org/10.1093/icesjms/fsac032.

Barnes AE, Davies JG, Martay B, Boersch-Supan PH, Harris SJ, Noble DG, Pearce-Higgins JW, Robinson RA (2023). *Rare and declining bird species benefit most from designating protected areas for conservation in the UK*. *Nature Ecology & Evolution* **7**, 92–101. https://doi.org/10.1038/s41559-022-01927-4.
546

## Acknowledgements

547

548

549

550

The authors would like to thank all the staff of the Écrins National Park and its scientific team for contributing to the dataset. They also appreciate advice of Denis Allard on using functional principal component analysis, and Loïc Houde for technical support with the high-performance computing cluster used to implement the models of this case study.

551

## Fundings

552

The authors declare that they have received no specific funding for this study.

553

## Conflict of interest disclosure

554

555

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.

556

## Data, script, code, and supplementary information availability

557

558

Data, scripts, and code can be accessed online (https://doi.org/10.5281/zenodo.12545052; cite Lasgorceux et al., 2024a).

559

560

Outputs from fitted models under various cross-validation scenarios are also available online (https://doi.org/10.5281/zenodo.12544212; cite Lasgorceux et al., 2024b).

561

562

Supplementary information is accessible online (https://doi.org/10.5281/zenodo.12541412; cite Lasgorceux et al., 2024c).

563

564

## References

565

566

567

568

Barnes AE, Davies JG, Martay B, Boersch-Supan PH, Harris SJ, Noble DG, Pearce-Higgins JW, Robinson RA (2023). *Rare and declining bird species benefit most from designating protected areas for conservation in the UK*. *Nature Ecology & Evolution* **7**, 92–101. https://doi.org/10.1038/s41559-022-01927-4.

569

570

Belmont J, Martino S, Illian J, Rue H (2024). *Spatio-temporal Occupancy Models with INLA*. https://doi.org/10.48550/arXiv.2403.10680.

571

572

573

Botella C, Joly A, Bonnet P, Munoz F, Monestiez P (2021). *Jointly estimating spatial sampling effort and habitat suitability for multiple species from opportunistic presence-only data*. *Methods in Ecology and Evolution* **12**, 933–945. https://doi.org/10.1111/2041-210X.13565.

574

575

576

Botella C, Joly A, Monestiez P, Bonnet P, Munoz F (2020). *Bias in presence-only niche models related to sampling effort and species niches: Lessons for background point selection*. *PLOS ONE* **15**(5), e0232078. https://doi.org/10.1371/journal.pone.0232078.

577

578

579

Bradter U, Mair L, Jönsson M, Knape J, Singer A, Snäll T (2018). *Can opportunistically collected Citizen Science data fill a data gap for habitat suitability models of less common species? Methods in Ecology and Evolution* **9**, 1667–1678. https://doi.org/10.1111/2041-210X.13012.

580

581

582

583

Bunz Y (2022). *Les programmes scientifiques de la faune vertébrée*. Parc national des Écrins, Domaine de Charance, Gap, France, p. 423. URL: https://www/documentation.eauetbiodiversite.fr/notice/les-programmes-scientifiques-de-la-faune-vertebree-edition-20220.
601

602

Belmont J, Martino S, Illian J, Rue H (2024). Spatio-temporal Occupancy Models with INLA. https://doi.org/10.48550/arXiv.2403.10680.

603

604

605

Botella C, Joly A, Bonnet P, Munoz F, Monestiez P (2021). Jointly estimating spatial sampling effort and habitat suitability for multiple species from opportunistic presence-only data. Methods in Ecology and Evolution **12**, 933-945. https://doi.org/10.1111/2041-210X.13565.

606

607

608

Botella C, Joly A, Monestiez P, Bonnet P, Munoz F (2020). Bias in presence-only niche models related to sampling effort and species niches: Lessons for background point selection. PLOS ONE **15**(5), e0232078. https://doi.org/10.1371/journal.pone.0232078.

609

610

611

Bradter U, Mair L, Jönsson M, Knape J, Singer A, Snäll T (2018). Can opportunistically collected Citizen Science data fill a data gap for habitat suitability models of less common species? Methods in Ecology and Evolution **9**, 1667-1678. https://doi.org/10.1111/2041-210X.13012.

612

613

614

615

Bunz Y (2022). Les programmes scientifiques de la faune vertébrée. Parc national des Écrins, Domaine de Charance, Gap, France, p. 423. URL: https://www documentation.eauetbiodiversite.fr/notice/les-programmes-scientifiques-de-la-faune-vertebree-edition-20220.

616

617

618

Chauvier Y, Zimmermann NE, Poggiato G, Bystrova D, Brun P, Thuiller W (2021). Novel methods to correct for observer and sampling bias in presence-only species distribution models. Global Ecology and Biogeography **30**, 2312-2325. https://doi.org/10.1111/geb.13383.

619

620

621

Coron C, Calenge C, Giraud C, Julliard R (2018). Bayesian estimation of species relative abundances and habitat preferences using opportunistic data. Environmental and Ecological Statistics **25**, 71-93. https://doi.org/10.1007/s10651-018-0398-2.

622

623

624

Donald PF, Green RE, Heath MF (2001). Agricultural intensification and the collapse of Europe's farmland bird populations. Proceedings of the Royal Society of London. Series B: Biological Sciences **268**, 25-29. https://doi.org/10.1098/rspb.2000.1325.

625

626

627

Dorazio RM (2014). Accounting for imperfect detection and survey bias in statistical analysis of presence-only data. Global Ecology and Biogeography **23**, 1472-1484. https://doi.org/10.1111/geb.12216.

628

629

Dray S, Dufour AB (2007). The ade4 Package: Implementing the Duality Diagram for Ecologists. Journal of Statistical Software **22**, 1-20. https://doi.org/10.18637/jss.v022.i04.

630

631

Duckworth GD, Altwegg R (2018). Effectiveness of protected areas for bird conservation depends on guild. Diversity and Distributions **24**, 1083-1091. https://doi.org/10.1111/ddi.12756.

632

633

634

Eftekhari A, Pasadakis D, Bollhöfer M, Scheidegger S, Schenk O (2021). Block-enhanced precision matrix estimation for large-scale datasets. Journal of Computational Science **53**. https://doi.org/10.1016/j.jocs.2021.101389.

635

636

637

Escamilla Molgora JM, Sedda L, Diggle P, Atkinson PM (2022). A joint distribution framework to improve presence-only species distribution models by exploiting opportunistic surveys. Journal of Biogeography **49**, 1176-1192. https://doi.org/10.1111/jbi.14365.

638

639

640

Farr MT, Green DS, Holekamp KE, Zipkin EF (2021). Integrating distance sampling and presence-only data to estimate species abundance. Ecology **102**, e03204. https://doi.org/10.1002/ecy.3204.

641

642

Fawcett T (2006). An introduction to ROC analysis. Pattern Recognition Letters **27**, 861-874. https://doi.org/10.1016/j.patrec.2005.10.010.
584

Chauvier Y, Zimmermann NE, Poggiato G, Bystrova D, Brun P, Thuiller W (2021). Novel methods to correct for observer and sampling bias in presence-only species distribution models. *Global Ecology and Biogeography* **30**, 2312-2325. https://doi.org/10.1111/geb.13383.

585

586

587

Coron C, Calenge C, Giraud C, Julliard R (2018). Bayesian estimation of species relative abundances and habitat preferences using opportunistic data. *Environmental and Ecological Statistics* **25**, 71-93. https://doi.org/10.1007/s10651-018-0398-2.

588

589

590

Donald PF, Green RE, Heath MF (2001). Agricultural intensification and the collapse of Europe's farmland bird populations. *Proceedings of the Royal Society of London. Series B: Biological Sciences* **268**, 25-29. https://doi.org/10.1098/rspb.2000.1325.

591

592

593

Dorazio RM (2014). Accounting for imperfect detection and survey bias in statistical analysis of presence-only data. *Global Ecology and Biogeography* **23**, 1472-1484. https://doi.org/10.1111/geb.12216.

594

595

596

Dray S, Dufour AB (2007). The ade4 Package: Implementing the Duality Diagram for Ecologists. *Journal of Statistical Software* **22**, 1-20. https://doi.org/10.18637/jss.v022.i04.

597

598

Duckworth GD, Altwegg R (2018). Effectiveness of protected areas for bird conservation depends on guild. *Diversity and Distributions* **24**, 1083-1091. https://doi.org/10.1111/ddi.12756.

599

600

Eftekhari A, Pasadakis D, Bollhöfer M, Scheidegger S, Schenk O (2021). Block-enhanced precision matrix estimation for large-scale datasets. *Journal of Computational Science* **53**. https://doi.org/10.1016/j.jocs.2021.101389.

601

602

603

Escamilla Molgora JM, Sedda L, Diggle P, Atkinson PM (2022). A joint distribution framework to improve presence-only species distribution models by exploiting opportunistic surveys. *Journal of Biogeography* **49**, 1176-1192. https://doi.org/10.1111/jbi.14365.

604

605

606

Farr MT, Green DS, Holekamp KE, Zipkin EF (2021). Integrating distance sampling and presence-only data to estimate species abundance. *Ecology* **102**, e03204. https://doi.org/10.1002/ecy.3204.

607

608

609

Fawcett T (2006). An introduction to ROC analysis. *Pattern Recognition Letters* **27**, 861-874. https://doi.org/10.1016/j.patrec.2005.10.010.

610

611

Fick SE, Hijmans RJ (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. *International Journal of Climatology* **37**, 4302-4315. https://doi.org/10.1002/joc.5086.

612

613

614

Fithian W, Elith J, Hastie T, Keith DA (2015). Bias correction in species distribution models: pooling survey and collection data for multiple species. *Methods in Ecology and Evolution* **6**, 424-438. https://doi.org/10.1111/2041-210X.12242.

615

616

617

Fithian W, Hastie T (2013). Finite-sample equivalence in statistical models for presence-only data. *The Annals of Applied Statistics* **7**, 1917-1939. https://doi.org/10.1214/13-A0AS667.

618

619

Flousek J, Telenský T, Hanzelka J, Reif J (2015). Population Trends of Central European Montane Birds Provide Evidence for Adverse Impacts of Climate Change on High-Altitude Species. *PLOS ONE* **10**, e0139465. https://doi.org/10.1371/journal.pone.0139465.

620

621

622

Fontaine B, Moussy C, Chiffard Carricaburu J, Dupuis J, Schmaltz L, Lorrillière R, Loïs G, Gaudard C, Couzi L (2020). Suivi des oiseaux communs en France 1989-2019 : 30 ans de suivis participatifs. MNHN- Centre d'Ecologie et des Sciences de la Conservation, LPO BirdLife France - Service Connaissance, Ministère de la Transition écologique et solidaire, p. 46. URL: https://www.vigienature.fr/sites/vigienature/files/atoms/files/syntheseoiseauxcommuns2020_final.pdf.

623

624

625

626

627
643

644

645

Fick SE, Hijmans RJ (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37, 4302-4315. https://doi.org/10.1002/joc.5086.

646

647

648

Fithian W, Elith J, Hastie T, Keith DA (2015). Bias correction in species distribution models: pooling survey and collection data for multiple species. Methods in Ecology and Evolution 6, 424-438. https://doi.org/10.1111/2041-210X.12242.

649

650

Fithian W, Hastie T (2013). Finite-sample equivalence in statistical models for presence-only data. The Annals of Applied Statistics 7, 1917-1939. https://doi.org/10.1214/13-A0AS667.

651

652

653

Flousek J, Telenský T, Hanzelka J, Reif J (2015). Population Trends of Central European Montane Birds Provide Evidence for Adverse Impacts of Climate Change on High-Altitude Species. PLOS ONE 10, e0139465. https://doi.org/10.1371/journal.pone.0139465.

654

655

656

657

658

659

Fontaine B, Moussy C, Chiffard Carricaburu J, Dupuis J, Schmaltz L, Lorrillière R, Loïs G, Gaudard C, Couzi L (2020). Suivi des oiseaux communs en France 1989-2019 : 30 ans de suivis participatifs. MNHN- Centre d'Ecologie et des Sciences de la Conservation, LPO BirdLife France - Service Connaissance, Ministère de la Transition écologique et solidaire., p. 46. URL: https://www.vigienature.fr/sites/vigienature/files/atoms/files/syntheseoiseauxcommuns2020_final.pdf.

660

661

662

El-Gabbas A, Dormann CF (2018). Improved species-occurrence predictions in data-poor regions: using large-scale data and bias correction with down-weighted Poisson regression and Maxent. Ecography 41, 1161-1172. https://doi.org/10.1111/ecog.03149.

663

664

665

Gaedke-Merzhäuser L, Niekerk J, Schenk O, Rue H (2023). Parallelized integrated nested Laplace approximations for fast Bayesian inference. Statistics and Computing 33, 25. https://doi.org/10.1007/s11222-022-10192-1.

666

667

Gilks WR, Richardson S, Spiegelhalter D (1995). Markov Chain Monte Carlo in Practice. New York, NY: Chapman and Hall/CRC. 512 pp. https://doi.org/10.1201/b14835.

668

669

670

Giraud C, Calenge C, Coron C, Julliard R (2016). Capitalizing on opportunistic data for monitoring relative abundances of species. Biometrics 72, 649-658. https://doi.org/10.1111/biom.12431.

671

672

673

Gregory RD, Skorpilova J, Vorisek P, Butler S (2019). An analysis of trends, uncertainty and species selection shows contrasting trends of widespread forest and farmland birds in Europe. Ecological Indicators 103, 676-687. https://doi.org/10.1016/j.ecolind.2019.04.064.

674

675

676

Gregory RD, Vorisek P, Van Strien A, Gmelig Meyling AW, Jiguet F, Fornasari L, Reif J, Chylarecki P, Burfield IJ (2007). Population trends of widespread woodland birds in Europe. Ibis 149 (s2), 78-97. https://doi.org/10.1111/j.1474-919X.2007.00698.x.

677

678

679

Harris I, Osborn TJ, Jones P, Lister D (2020). Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. Scientific Data 7, 109. https://doi.org/10.1038/s41597-020-0453-3.

680

681

682

683

Haupert S, Sèbe F, Sueur J (2023). Physics-based model to predict the acoustic detection distance of terrestrial autonomous recording units over the diel cycle and across seasons: Insights from an Alpine and a Neotropical forest. Methods in Ecology and Evolution 14, 614-630. https://doi.org/10.1111/2041-210X.14020.

684

685

686

Heldbjerg H, Sunde P, Fox AD (2018). Continuous population declines for specialist farmland birds 1987-2014 in Denmark indicates no halt in biodiversity loss in agricultural habitats. Bird Conservation International 28, 278-292. https://doi.org/10.1017/S0959270916000654.
628

629

630

Gaedke-Merzhäuser L, Niekerk J, Schenk O, Rue H (2023). *Parallelized integrated nested Laplace approximations for fast Bayesian inference*. *Statistics and Computing* **33**, 25. https://doi.org/10.1007/s11222-022-10192-1.

631

632

Gilks WR, Richardson S, Spiegelhalter D (1995). *Markov Chain Monte Carlo in Practice*. New York, NY: Chapman and Hall/CRC. 512 pp. https://doi.org/10.1201/b14835.

633

634

635

Giraud C, Calenge C, Coron C, Julliard R (2016). *Capitalizing on opportunistic data for monitoring relative abundances of species*. *Biometrics* **72**, 649–658. https://doi.org/10.1111/biom.12431.

636

637

638

Gregory RD, Skorpilova J, Vorisek P, Butler S (2019). *An analysis of trends, uncertainty and species selection shows contrasting trends of widespread forest and farmland birds in Europe*. *Ecological Indicators* **103**, 676–687. https://doi.org/10.1016/j.ecolind.2019.04.064.

639

640

641

Gregory RD, Vorisek P, Van Strien A, Gmelig Meyling AW, Jiguet F, Fornasari L, Reif J, Chylarecki P, Burfield IJ (2007). *Population trends of widespread woodland birds in Europe*. *Ibis* **149** (s2), 78–97. https://doi.org/10.1111/j.1474-919X.2007.00698.x.

642

643

644

Harris I, Osborn TJ, Jones P, Lister D (2020). *Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset*. *Scientific Data* **7**, 109. https://doi.org/10.1038/s41597-020-0453-3.

645

646

647

648

Haupert S, Sèbe F, Sueur J (2023). *Physics-based model to predict the acoustic detection distance of terrestrial autonomous recording units over the diel cycle and across seasons: Insights from an Alpine and a Neotropical forest*. *Methods in Ecology and Evolution* **14**, 614–630. https://doi.org/10.1111/2041-210X.14020.

649

650

651

Heldbjerg H, Sunde P, Fox AD (2018). *Continuous population declines for specialist farmland birds 1987-2014 in Denmark indicates no halt in biodiversity loss in agricultural habitats*. *Bird Conservation International* **28**, 278–292. https://doi.org/10.1017/S0959270916000654.

652

653

654

655

Henckel L, Bradter U, Jönsson M, Isaac NJB, Snäll T (2020). *Assessing the usefulness of citizen science data for habitat suitability modelling: Opportunistic reporting versus sampling based on a systematic protocol*. *Diversity and Distributions* **26**, 1276–1290. https://doi.org/10.1111/ddi.13128.

656

657

Hijmans RJ (2024). *terra: Spatial Data Analysis*. R package version 1.7-71. URL: https://CRAN.R-project.org/package=terra.

658

659

660

Horns JJ, Adler FR, Şekercioğlu ÇH (2018). *Using opportunistic citizen science data to estimate avian population trends*. *Biological Conservation* **221**, 151–159. https://doi.org/10.1016/j.biocon.2018.02.027.

661

IGN (2018). *RGE ALTI® [Data set]*. URL: https://geoservices.ign.fr/rgealti.

662

663

664

665

Illian JB, Martino S, Sørbye SH, Gallego-Fernández JB, Zunzunegui M, Esquivias MP, Travis JMJ (2013). *Fitting complex ecological point process models with integrated nested Laplace approximation*. *Methods in Ecology and Evolution* **4**, 305–315. https://doi.org/10.1111/2041-210x.12017.

666

667

668

Inger R, Gregory R, Duffy JP, Stott I, Voříšek P, Gaston KJ (2015). *Common European birds are declining rapidly while less abundant species' numbers are rising*. *Ecology Letters* **18**, 28–36. https://doi.org/10.1111/ele.12387.

669

670

Inglada J, Vincent A, Thierion V (2018). *Theia OSO Land Cover Map 2018 [Data set]*. https://doi.org/10.5281/zenodo.3613415.
687

688

689

690

Henckel L, Bradter U, Jönsson M, Isaac NJB, Snäll T (2020). Assessing the usefulness of citizen science data for habitat suitability modelling: Opportunistic reporting versus sampling based on a systematic protocol. Diversity and Distributions **26**, 1276-1290. https://doi.org/10.1111/ddi.13128.

691

692

Hijmans RJ (2024). terra: Spatial Data Analysis. R package version 1.7-71. URL: https://CRAN.R-project.org/package=terra.

693

694

695

Horns JJ, Adler FR, Şekercioğlu ÇH (2018). Using opportunistic citizen science data to estimate avian population trends. Biological Conservation **221**, 151-159. https://doi.org/10.1016/j.biocon.2018.02.027.

696

IGN (2018). RGE ALTI® [Data set]. URL: https://geoservices.ign.fr/rgealti.

697

698

699

700

Illian JB, Martino S, Sørbye SH, Gallego-Fernández JB, Zunzunegui M, Esquivias MP, Travis JMJ (2013). Fitting complex ecological point process models with integrated nested Laplace approximation. Methods in Ecology and Evolution **4**, 305-315. https://doi.org/10.1111/2041-210x.12017.

701

702

703

Inger R, Gregory R, Duffy JP, Stott I, Voříšek P, Gaston KJ (2015). Common European birds are declining rapidly while less abundant species' numbers are rising. Ecology Letters **18**, 28-36. https://doi.org/10.1111/ele.12387.

704

705

Inglada J, Vincent A, Thierion V (2018). Theia OSO Land Cover Map 2018 [Data set]. https://doi.org/10.5281/zenodo.3613415.

706

707

708

IPBES (2019). Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. Bonn, Germany. 56 pages. URL: https://doi.org/10.5281/zenodo.3553579.

709

710

Janžekovič F, Novak T (2012). PCA – A Powerful Method for Analyze Ecological Niches. In: Principal Component Analysis. 8. Rijeka: IntechOpen. https://doi.org/10.5772/38538.

711

712

713

Jiguet F, Devictor V, Julliard R, Couvet D (2012). French citizens monitoring ordinary birds provide tools for conservation and ecological sciences. Acta Oecologica **44**, 58-66. https://doi.org/10.1016/j.actao.2011.05.003.

714

715

716

Jin Huang, Ling C (2005). Using AUC and accuracy in evaluating learning algorithms. IEEE Transactions on Knowledge and Data Engineering **17**, 299-310. https://doi.org/10.1109/TKDE.2005.50.

717

718

Jung Y (2018). Multiple predicting K -fold cross-validation for model selection. Journal of Nonparametric Statistics **30**, 197-215. https://doi.org/10.1080/10485252.2017.1404598.

719

720

721

722

Kamp J, Frank C, Trautmann S, Busch M, Dröschmeister R, Flade M, Gerlach B, Karthäuser J, Kunz F, Mitschke A, Schwarz J, Sudfeldt C (2021). Population trends of common breeding birds in Germany 1990-2018. Journal of Ornithology **162**, 1-15. https://doi.org/10.1007/s10336-020-01830-4.

723

724

725

726

Kervellec M, Milleret C, Vanpé C, Quenette PY, Sentilles J, Palazón S, Jordana IA, Jato R, Elósegui Irurtia MM, Gimenez O (2023). Integrating opportunistic and structured non-invasive surveys with spatial capture-recapture models to map connectivity of the Pyrenean brown bear population. Biological Conservation **278**, 109875. https://doi.org/10.1016/j.biocon.2022.109875.

727

728

729

730

Kéry M, Royle JA, Schmid H, Schaub M, Volet B, Häfliger G, Zbinden N (2010). Site-occupancy distribution modeling to correct population-trend estimates derived from opportunistic observations. Conservation Biology **24**, 1388-1397. https://doi.org/10.1111/j.1523-1739.2010.01479.x.
671

672

673

IPBES (2019). Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. Bonn, Germany. 56 pages. URL: https://doi.org/10.5281/zenodo.3553579.

674

675

Janžekovič F, Novak T (2012). PCA – A Powerful Method for Analyze Ecological Niches. In: Principal Component Analysis. 8. Rijeka: IntechOpen. https://doi.org/10.5772/38538.

676

677

678

Jiguet F, Devictor V, Julliard R, Couvet D (2012). French citizens monitoring ordinary birds provide tools for conservation and ecological sciences. Acta Oecologica 44, 58–66. https://doi.org/10.1016/j.actao.2011.05.003.

679

680

681

Jin Huang, Ling C (2005). Using AUC and accuracy in evaluating learning algorithms. IEEE Transactions on Knowledge and Data Engineering 17, 299–310. https://doi.org/10.1109/TKDE.2005.50.

682

683

Jung Y (2018). Multiple predicting K -fold cross-validation for model selection. Journal of Nonparametric Statistics 30, 197–215. https://doi.org/10.1080/10485252.2017.1404598.

684

685

686

687

Kamp J, Frank C, Trautmann S, Busch M, Dröschmeister R, Flade M, Gerlach B, Karthäuser J, Kunz F, Mitschke A, Schwarz J, Sudfeldt C (2021). Population trends of common breeding birds in Germany 1990–2018. Journal of Ornithology 162, 1–15. https://doi.org/10.1007/s10336-020-01830-4.

688

689

690

691

Kéry M, Royle JA, Schmid H, Schaub M, Volet B, Häfliger G, Zbinden N (2010). Site-occupancy distribution modeling to correct population-trend estimates derived from opportunistic observations. Conservation Biology 24, 1388–1397. https://doi.org/10.1111/j.1523-1739.2010.01479.x.

692

693

694

Knape J (2016). Decomposing trends in Swedish bird populations using generalized additive mixed models. Journal of Applied Ecology 53, 1852–1861. https://doi.org/10.1111/j.1365-2664.12720.

695

696

697

Lasgorceux F, Papaïx J, Bunz Y, Combrisson D, Opitz T (2024a). Data, scripts and code for 'Space-time species distribution modeling with opportunistic presence-only data: a case study of passerines in a protected area'. https://doi.org/10.5281/zenodo.12545052.

698

699

700

701

Lasgorceux F, Papaïx J, Bunz Y, Combrisson D, Opitz T (2024b). Outputs from fitted models across the cross-validation scenarios for 'Space-time species distribution modeling with opportunistic presence-only data: a case study of passerines in a protected area'. https://doi.org/10.5281/zenodo.12544212.

702

703

704

Lasgorceux F, Papaïx J, Bunz Y, Combrisson D, Opitz T (2024c). Supplementary information for 'Space-time species distribution modeling with opportunistic presence- only data: a case study of passerines in a protected area'. https://doi.org/10.5281/zenodo.12541412.

705

706

707

708

Lehikoinen A, Brotons L, Calladine J, Campedelli T, Escandell V, Flousek J, Grueneberg C, Haas F, Harris S, Herrando S, Husby M, Jiguet F, Kålås JA, Lindström Å, Lorrillière R, Molina B, Pladevall C, Calvi G, Sattler T, Schmid H, et al. (2019). Declining population trends of European mountain birds. Global Change Biology 25, 577–588. https://doi.org/10.1111/gcb.14522.

709

710

711

712

Lindgren F, Rue H, Lindström J (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society Series B: Statistical Methodology 73, 423–498. https://doi.org/10.1111/j.1467-9868.2011.00777.x.
731

732

733

Knape J (2016). Decomposing trends in Swedish bird populations using generalized additive mixed models. Journal of Applied Ecology 53, 1852-1861. https://doi.org/10.1111/j.1365-2664.12720.

734

735

736

737

Lehikoinen A, Brotons L, Calladine J, Campedelli T, Escandell V, Flousek J, Grueneberg C, Haas F, Harris S, Herrando S, Husby M, Jiguet F, Kålås JA, Lindström Å, Lorrillière R, Molina B, Pladevall C, Calvi G, Sattler T, Schmid H, et al. (2019). Declining population trends of European mountain birds. Global Change Biology 25, 577-588. https://doi.org/10.1111/gcb.14522.

738

739

740

741

Lindgren F, Rue H, Lindström J (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society Series B: Statistical Methodology 73, 423-498. https://doi.org/10.1111/j.1467-9868.2011.00777.x.

742

743

744

Link WA, Cam E, Nichols JD, Cooch EG (2002). Of Bugs and Birds: Markov Chain Monte Carlo for Hierarchical Modeling in Wildlife Research. The Journal of Wildlife Management 66, 277. https://doi.org/10.2307/3803160.

745

746

747

McVicar TR, Körner C (2013). On the use of elevation, altitude, and height in the ecological and climatological literature. Oecologia 171, 335-337. https://doi.org/10.1007/s00442-012-2416-7.

748

749

750

Moreira GA, Menezes R, Wise L (2024). Presence-Only for Marked Point Process Under Preferential Sampling. Journal of Agricultural, Biological and Environmental Statistics 29, 92-109. https://doi.org/10.1007/s13253-023-00558-x.

751

752

753

Murtagh F, Legendre P (2014). Ward's Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward's Criterion? Journal of Classification 31, 274-295. https://doi.org/10.1007/s00357-014-9161-z.

754

755

756

Newton I (2004). The recent declines of farmland bird populations in Britain: an appraisal of causal factors and conservation actions. Ibis 146, 579-600. https://doi.org/10.1111/j.1474-919X.2004.00375.x.

757

758

759

Noël F, Combrisson D, Geoffroy JJ, Nicolas MG (2023). Synthèse des connaissances sur les Isopodes et Diplopodes terrestres du Parc national des Écrins (SE France). Naturae, 151-170. https://doi.org/10.5852/naturae2023a9.

760

761

762

Palacín C, Alonso JC (2018). Failure of EU Biodiversity Strategy in Mediterranean farmland protected areas. Journal for Nature Conservation 42, 62-66. https://doi.org/10.1016/j.jnc.2018.02.008.

763

764

765

Pasadakis D, Bollhöfer M, Schenk O (2023). Sparse Quadratic Approximation for Graph Learning. IEEE Transactions on Pattern Analysis and Machine Intelligence 45, 11256-11269. https://doi.org/10.1109/TPAMI.2023.3263969.

766

767

Pedersen EJ, Miller DL, Simpson GL, Ross N (2019). Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7, e6876. https://doi.org/10.7717/peerj.6876.

768

769

770

Phillips SJ, Anderson RP, Schapire RE (2006). Maximum entropy modeling of species geographic distributions. Ecological Modelling 190, 231-259. https://doi.org/10.1016/j.ecolmodel.2005.03.026.

771

772

773

Phillips SJ, Dudík M, Elith J, Graham CH, Lehmann A, Leathwick J, Ferrier S (2009). Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecological Applications 19, 181-197. https://doi.org/10.1890/07-2153.1.
713

714

715

Link WA, Cam E, Nichols JD, Cooch EG (2002). Of Bugs and Birds: Markov Chain Monte Carlo for Hierarchical Modeling in Wildlife Research. The Journal of Wildlife Management **66**, 277. https://doi.org/10.2307/3803160.

716

717

718

McVicar TR, Körner C (2013). On the use of elevation, altitude, and height in the ecological and climatological literature. Oecologia **171**, 335-337. https://doi.org/10.1007/s00442-012-2416-7.

719

720

721

Moreira GA, Menezes R, Wise L (2024). Presence-Only for Marked Point Process Under Preferential Sampling. Journal of Agricultural, Biological and Environmental Statistics **29**, 92-109. https://doi.org/10.1007/s13253-023-00558-x.

722

723

724

Murtagh F, Legendre P (2014). Ward's Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward's Criterion? Journal of Classification **31**, 274-295. https://doi.org/10.1007/s00357-014-9161-z.

725

726

727

Newton I (2004). The recent declines of farmland bird populations in Britain: an appraisal of causal factors and conservation actions. *Ibis* **146**, 579-600. https://doi.org/10.1111/j.1474-919X.2004.00375.x.

728

729

730

Noël F, Combrisson D, Geoffroy JJ, Nicolas MG (2023). Synthèse des connaissances sur les Isopodes et Diplopodes terrestres du Parc national des Écrins (SE France). *Naturae*, 151-170. https://doi.org/10.5852/naturae2023a9.

731

732

733

Palacín C, Alonso JC (2018). Failure of EU Biodiversity Strategy in Mediterranean farmland protected areas. Journal for Nature Conservation **42**, 62-66. https://doi.org/10.1016/j.jnc.2018.02.008.

734

735

736

Pasadakis D, Bollhöfer M, Schenk O (2023). Sparse Quadratic Approximation for Graph Learning. IEEE Transactions on Pattern Analysis and Machine Intelligence **45**, 11256-11269. https://doi.org/10.1109/TPAMI.2023.3263969.

737

738

Pedersen EJ, Miller DL, Simpson GL, Ross N (2019). Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ **7**, e6876. https://doi.org/10.7717/peerj.6876.

739

740

741

Phillips SJ, Anderson RP, Schapire RE (2006). Maximum entropy modeling of species geographic distributions. Ecological Modelling **190**, 231-259. https://doi.org/10.1016/j.ecolmodel.2005.03.026.

742

743

744

Phillips SJ, Dudík M, Elith J, Graham CH, Lehmann A, Leathwick J, Ferrier S (2009). Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecological Applications **19**, 181-197. https://doi.org/10.1890/07-2153.1.

745

746

R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria. URL: https://www.r-project.org.

747

748

Ramsay JO, Silverman BW (2005). Functional Data Analysis. Springer Series in Statistics. New York, NY: Springer. https://doi.org/10.1007/b98888.

749

750

Ramsay J, Hooker G, Graves S (2009). Functional Data Analysis with R and MATLAB. New York, NY: Springer. https://doi.org/10.1007/978-0-387-98185-7.

751

752

753

Reif J, Storch D, Voříšek P, Šťastný K, Bejček V (2008a). Bird-habitat associations predict population trends in central European forest and farmland birds. Biodiversity and Conservation **17**, 3307-3319. https://doi.org/10.1007/s10531-008-9430-4.

754

755

756

Reif J, Voříšek P, Šťastný K, Bejček V, Petr J (2008b). Agricultural intensification and farmland birds: new insights from a central European country. *Ibis* **150**, 596-605. https://doi.org/10.1111/j.1474-919X.2008.00829.x.
774
775

R Core Team (2024). *R: A Language and Environment for Statistical Computing*. R Foundation for
Statistical Computing. Vienna, Austria. URL: https://www.r-project.org.

776
777

Ramsay JO, Silverman BW (2005). Functional Data Analysis. Springer Series in Statistics. New York, NY: Springer. https://doi.org/10.1007/b98888.

778
779

Ramsay J, Hooker G, Graves S (2009). Functional Data Analysis with R and MATLAB. New York, NY: Springer. https://doi.org/10.1007/978-0-387-98185-7.

780
781
782

Reif J, Storch D, Voříšek P, Šťastný K, Bejček V (2008a). Bird-habitat associations predict population trends in central European forest and farmland birds. Biodiversity and Conservation **17**, 3307-3319. https://doi.org/10.1007/s10531-008-9430-4.

783
784
785

Reif J, Voříšek P, Šťastný K, Bejček V, Petr J (2008b). Agricultural intensification and farmland birds: new insights from a central European country. Ibis **150**, 596-605. https://doi.org/10.1111/j.1474-919X.2008.00829.x.

786
787
788

Renner IW, Elith J, Baddeley A, Fithian W, Hastie T, Phillips SJ, Popovic G, Warton DI (2015). Point process models for presence-only analysis. Methods in Ecology and Evolution **6**, 366-379. https://doi.org/10.1111/2041-210X.12352.

789
790
791

Renner IW, Warton DI (2013). Equivalence of MAXENT and Poisson Point Process Models for Species Distribution Modeling in Ecology. Biometrics **69**, 274-281. https://doi.org/10.1111/j.1541-0420.2012.01824.x.

792
793
794
795
796

Rigal S, Dakos V, Alonso H, Auniņš A, Benkő Z, Brotons L, Chodkiewicz T, Chylarecki P, Carli E, Moral JC, Domşa C, Escandell V, Fontaine B, Foppen R, Gregory R, Harris S, Herrando S, Husby M, Ieronymidou C, Jiguet F, et al. (2023). Farmland practices are driving bird population decline across Europe. Proceedings of the National Academy of Sciences **120**, e2216573120. https://doi.org/10.1073/pnas.2216573120.

797
798
799
800

Rue H, Martino S, Chopin N (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society Series B: Statistical Methodology **71**, 319-392. https://doi.org/10.1111/j.1467-9868.2008.00700.x.

801
802
803

Sanderson FJ, Kucharz M, Jobda M, Donald PF (2013). Impacts of agricultural intensification and abandonment on farmland birds in Poland following EU accession. Agriculture, Ecosystems & Environment **168**, 16-24. https://doi.org/10.1016/j.agee.2013.01.015.

804
805
806
807

Schulze ED, Craven D, Durso AM, Reif J, Guderle M, Kroiher F, Hennig P, Weiserbs A, Schall P, Ammer C, Eisenhauer N (2019). Positive association between forest management, environmental change, and forest bird abundance. Forest Ecosystems **6**, 3. https://doi.org/10.1186/s40663-019-0160-8.

808
809
810

Silva JP, Correia R, Alonso H, Martins RC, D'Amico M, Delgado A, Sampaio H, Godinho C, Mora eira F (2018). EU protected area network did not prevent a country wide population decline in a threatened grassland bird. PeerJ **6**, e4284. https://doi.org/10.7717/peerj.4284.

811
812
813

Snäll T, Kindvall O, Nilsson J, Pärt T (2011). Evaluating citizen-based presence data for bird monitoring. Biological Conservation **144**, 804-810. https://doi.org/10.1016/j.biocon.2010.11.010.

814
815
816
817

Soriano-Redondo A, Jones-Todd CM, Bearhop S, Hilton GM, Lock L, Stanbury A, Votier SC, Il-lian JB (2019). Understanding species distribution in dynamic populations: a new approach using spatio-temporal point process models. Ecography **42**, 1092-1102. https://doi.org/10.1111/ecog.03771.
757 Renner IW, Elith J, Baddeley A, Fithian W, Hastie T, Phillips SJ, Popovic G, Warton DI (2015). Point process models for presence-only analysis. Methods in Ecology and Evolution 6, 366-379. https://doi.org/10.1111/2041-210X.12352.

760 Renner IW, Warton DI (2013). Equivalence of MAXENT and Poisson Point Process Models for Species Distribution Modeling in Ecology. Biometrics 69, 274-281. https://doi.org/10.1111/j.1541-0420.2012.01824.x.

763 Rigal S, Dakos V, Alonso H, Auniņš A, Benkő Z, Brotons L, Chodkiewicz T, Chylarecki P, Carli E, Moral JC, Domşa C, Escandell V, Fontaine B, Foppen R, Gregory R, Harris S, Herrando S, Husby M, Ieronymidou C, Jiguet F, et al. (2023). Farmland practices are driving bird population decline across Europe. Proceedings of the National Academy of Sciences 120, e2216573120. https://doi.org/10.1073/pnas.2216573120.

768 Rue H, Martino S, Chopin N (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society Series B: Statistical Methodology 71, 319-392. https://doi.org/10.1111/j.1467-9868.2008.00700.x.

772 Sanderson FJ, Kucharz M, Jobda M, Donald PF (2013). Impacts of agricultural intensification and abandonment on farmland birds in Poland following EU accession. Agriculture, Ecosystems & Environment 168, 16-24. https://doi.org/10.1016/j.agee.2013.01.015.

775 Schulze ED, Craven D, Durso AM, Reif J, Guderle M, Kroiher F, Hennig P, Weiserbs A, Schall P, Ammer C, Eisenhauer N (2019). Positive association between forest management, environmental change, and forest bird abundance. Forest Ecosystems 6, 3. https://doi.org/10.1186/s40663-019-0160-8.

779 Silva JP, Correia R, Alonso H, Martins RC, D'Amico M, Delgado A, Sampaio H, Godinho C, Mora eira F (2018). EU protected area network did not prevent a country wide population decline in a threatened grassland bird. PeerJ 6, e4284. https://doi.org/10.7717/peerj.4284.

782 Snäll T, Kindvall O, Nilsson J, Pärt T (2011). Evaluating citizen-based presence data for bird monitoring. Biological Conservation 144, 804-810. https://doi.org/10.1016/j.biocon.2010.11.010.

785 Soriano-Redondo A, Jones-Todd CM, Bearhop S, Hilton GM, Lock L, Stanbury A, Votier SC, Il-lian JB (2019). Understanding species distribution in dynamic populations: a new approach using spatio-temporal point process models. Ecography 42, 1092-1102. https://doi.org/10.1111/ecog.03771.

789 Timmers R, Kuijk M, Verweij PA, Ghazoul J, Hautier Y, Laurance WF, Arriaga-Weiss SL, Askins RA, Battisti C, Berg Å, Daily GC, Estades CF, Frank B, Kurosawa R, Pojar RA, Woinarski JC, Soons MB (2022). Conservation of birds in fragmented landscapes requires protected areas. Frontiers in Ecology and the Environment 20, 361-369. https://doi.org/10.1002/fee.2485.

793 Traba J, Morales MB (2019). The decline of farmland birds in Spain is strongly associated to the loss of fallowland. Scientific Reports 9, 9473. https://doi.org/10.1038/s41598-019-45854-0.

795 Valavi R, Guillera-Arroita G, Lahoz-Monfort JJ, Elith J (2022). Predictive performance of presence-only species distribution models: a benchmark study with reproducible code. Ecological Monographs 92, e01486. https://doi.org/10.1002/ecm.1486.

798 Van Strien AJ, Van Swaay CA, Termaat T (2013). Opportunistic citizen science data of animal species produce reliable estimates of distribution trends if analysed with occupancy models. Journal of Applied Ecology 50, 1450-1458. https://doi.org/10.1111/1365-2664.12158.
818

819

820

821

822

823

824

825

826

827

828

829

830

831

832

833

834

835

836

837

838

839

840

841

842

843

844

845

846

847

848

Timmers R, Kuijk M, Verweij PA, Ghazoul J, Hautier Y, Laurance WF, Arriaga-Weiss SL, Askins RA, Battisti C, Berg Å, Daily GC, Estades CF, Frank B, Kurosawa R, Pojar RA, Woinarski JC, Soons MB (2022). *Conservation of birds in fragmented landscapes requires protected areas. Frontiers in Ecology and the Environment* **20**, 361–369. https://doi.org/10.1002/fee.2485.

Traba J, Morales MB (2019). *The decline of farmland birds in Spain is strongly associated to the loss of fallowland. Scientific Reports* **9**, 9473. https://doi.org/10.1038/s41598-019-45854-0.

Valavi R, Guillera-Aroita G, Lahoz-Monfort JJ, Elith J (2022). *Predictive performance of presence-only species distribution models: a benchmark study with reproducible code. Ecological Monographs* **92**, e01486. https://doi.org/10.1002/ecm.1486.

Van Strien AJ, Van Swaay CA, Termaat T (2013). *Opportunistic citizen science data of animal species produce reliable estimates of distribution trends if analysed with occupancy models. Journal of Applied Ecology* **50**, 1450–1458. https://doi.org/10.1111/1365-2664.12158.

Ver Hoef JM, Johnson D, Angliss R, Higham M (2021). *Species density models from opportunistic citizen science data. Methods in Ecology and Evolution* **12**, 1911–1925. https://doi.org/10.1111/2041-210X.13679.

Warton DI, Renner IW, Ramp D (2013). *Model-Based Control of Observer Bias for the Analysis of Presence-Only Data in Ecology. PLoS ONE* **8**, e79168. https://doi.org/10.1371/journal.pone.0079168.

Warton DI, Shepherd LC (2010). *Poisson point process models solve the “pseudo-absence problem” for presence-only data in ecology. The Annals of Applied Statistics* **4**, https://doi.org/10.1214/10-A0AS331.

Wikle CK (2003). *Hierarchical Bayesian Models for Predicting the Spread of Ecological Processes. Ecology* **84**, 1382–1394. https://doi.org/10.1890/0012-9658(2003)084[1382:HBMFPT]2.0.CO;2.

Wretenberg J, Lindström Å, Svensson S, Thierfelder T, Pärt T (2006). *Population trends of farmland birds in Sweden and England: similar trends but different patterns of agricultural intensification. Journal of Applied Ecology* **43**, 1110–1120. https://doi.org/10.1111/j.1365-2664.2006.01216.x.

Zamora R, Barea-Azcón JM (2015). *Long-Term Changes in Mountain Passerine Bird Communities in the Sierra Nevada (Southern Spain): A 30-Year Case Study. Ardeola* **62**, 3–18. https://doi.org/10.13157/arla.62.1.2015.3.
801

802

803

Warton DI, Renner IW, Ramp D (2013). *Model-Based Control of Observer Bias for the Analysis of Presence-Only Data in Ecology*. PLoS ONE 8, e79168. https://doi.org/10.1371/journal.pone.0079168.

804

805

806

Warton DI, Shepherd LC (2010). *Poisson point process models solve the “pseudo-absence problem” for presence-only data in ecology*. The Annals of Applied Statistics 4. https://doi.org/10.1214/10-AOAS331.

807

808

809

Wikle CK (2003). *Hierarchical Bayesian Models for Predicting the Spread of Ecological Processes*. Ecology 84, 1382-1394. https://doi.org/10.1890/0012-9658(2003)084[1382:HBMFPT].2.0.CO;2.

810

811

812

813

Wretenberg J, Lindström Å, Svensson S, Thierfelder T, Pärt T (2006). *Population trends of farmland birds in Sweden and England: similar trends but different patterns of agricultural intensification*. Journal of Applied Ecology 43, 1110-1120. https://doi.org/10.1111/j.1365-2664.2006.01216.x.

814

815

816

Zamora R, Barea-Azcón JM (2015). *Long-Term Changes in Mountain Passerine Bird Communities in the Sierra Nevada (Southern Spain): A 30-Year Case Study*. Ardeola 62, 3-18. https://doi.org/10.13157/arla.62.1.2015.3.
1. Passerines data

<figure><img src="image_7.png" /></figure>

Source: Écrins National Park

**Figure 6** – Number of reported opportunistic passerine occurrences in the Écrins National Park between 1974 and 2021. Three passerines were excluded from the barplot: Common Raven (Corvus corax), Red-billed Chough (Pyrrhocorax pyrrhocorax), Alpine Chough (Pyrrhocorax graculus) for a total of 102,513 occurrences.
## 2. OSO covariates

The OSO data are raster data produced at 10m and 20m spatial resolution. These two outputs are provided because the Sentinel-2 spectral bands used in the production process are acquired at 10m (visible and near-infrared) and 20m (red-edge, near and medium infrared) (Inglada, 2018). The nomenclature includes 24 classes. We then apply the nearest neighbor strategy to adjust the values to a 5m resolution, aligning with the elevation data resolution used in this study (IGN, 2018). We present in Table 1 the percentage of coverage for each OSO class in the Écrins National Park. The percentage of coverage for a chosen class is computed as the number of spatial cells of this class over the total number of spatial cells, based on a 5m regular grid of the Écrins National Park domain.

<table><caption>**Table 1** – OSO data nomenclature including percentage coverage at a 5m resolution Écrins National Park. Classes with a percentage of coverage greater than 1% are in bold for easier readability.</caption><thead><tr><td>Number</td><td>Nomenclature</td><td>Percentage</td></tr></thead><tbody><tr><td>1</td><td>Dense buildings</td><td>&lt;0.1%</td></tr><tr><td>2</td><td>Sparse buildings</td><td>0.6%</td></tr><tr><td>3</td><td>Industrial and commercial areas</td><td>0.3%</td></tr><tr><td>4</td><td>Road surfaces</td><td>&lt;0.1%</td></tr><tr><td>5</td><td>Rapeseed</td><td>&lt;0.1%</td></tr><tr><td>6</td><td>Cereal with straw</td><td>0.1%</td></tr><tr><td>7</td><td>Protein crops</td><td>&lt;0.1%</td></tr><tr><td>8</td><td>Soybean</td><td>&lt;0.1%</td></tr><tr><td>9</td><td>Sunflower</td><td>&lt;0.1%</td></tr><tr><td>10</td><td>Corn</td><td>0.1%</td></tr><tr><td>11</td><td>Rice</td><td>&lt;0.1%</td></tr><tr><td>12</td><td>Tubers/roots</td><td>&lt;0.1%</td></tr><tr><td>13</td><td>Meadows</td><td>2.1%</td></tr><tr><td>14</td><td>Orchards</td><td>&lt;0.1%</td></tr><tr><td>15</td><td>Vineyards</td><td>&lt;0.1%</td></tr><tr><td>16</td><td>Deciduous forests</td><td>3.7%</td></tr><tr><td>17</td><td>Coniferous forests</td><td>16.0%</td></tr><tr><td>18</td><td>Grasslands</td><td>40.1%</td></tr><tr><td>19</td><td>Woody heaths</td><td>5.3%</td></tr><tr><td>20</td><td>Mineral surfaces</td><td>27.7%</td></tr><tr><td>21</td><td>Beaches and dunes</td><td>0.1%</td></tr><tr><td>22</td><td>Glaciers or snow</td><td>2.8%</td></tr><tr><td>23</td><td>Water</td><td>1.1%</td></tr><tr><td>0</td><td>Other</td><td>&lt;0.1%</td></tr></tbody></table>

Some classes are nearly absent in the Écrins National Park. Hence, the first four classes: Dense buildings, Sparse buildings, Industrial and commercial areas and Road surfaces are gathered into the Urban terminology. The following classes: Rapeseed, Cereal with straw, Legumes, Soybean, Sunflower, Corn, Rice, Tubers/roots, Orchards, Vineyards are gathered into the Crops terminology. Due to their very low percentage of surface area, the classes Beaches and dunes and Other are excluded from the analysis.

Hence, the land cover land used classification used in the article is composed of 10 classes and is summarized in Table 2.
<table><caption>**Table 2** – OSO data nomenclature adjusted for the Écrins National Park, including percentage coverage at a 5m resolution in the Écrins National Park.</caption><thead><tr><td>Nomenclature</td><td>Percentage</td></tr></thead><tbody><tr><td>Urban</td><td>0.9%</td></tr><tr><td>Crops</td><td>0.3%</td></tr><tr><td>Meadows</td><td>2.1%</td></tr><tr><td>Deciduous forests</td><td>3.7%</td></tr><tr><td>Coniferous forests</td><td>16.0%</td></tr><tr><td>Grasslands</td><td>40.1%</td></tr><tr><td>Woody heaths</td><td>5.3%</td></tr><tr><td>Mineral surfaces</td><td>27.7%</td></tr><tr><td>Glaciers or snow</td><td>2.8%</td></tr><tr><td>Water</td><td>1.1%</td></tr></tbody></table>

Using a spatial grid with a resolution of 500m, we calculated the coverage percentage of these 10 classes in each spatial cell to construct the OSO covariates, which formed the basis of the Principal Component Analysis.

## References

871 Inglada, J, Vincent, A, & Thierion, V (2018). Theia OSO Land Cover Map 2018 [Data set]. Zenodo.

872 https://doi.org/10.5281/zenodo.3613415

874 IGN (2018). RGE ALTI®. Version 2.0.
3. PCA figures

<figure><img src="image_8.png" /><figcaption>Figure 7 – Correlation circles of axis 1 and 2 of the PCA.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8 – Contribution of each spatial covariate to the axis 1 (a) and 2 (b) of the PCA.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9 - Correlation circles of axis 3 and 4 of the PCA.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10 - Contribution of each spatial covariate to the axis 3 (a) and 4 (b) of the PCA.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11 - Correlation circles of axis 5 and 6 of the PCA.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12 - Contribution of each spatial covariate to the axis 5 (a) and 6 (b) of the PCA.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13** – Map of passerine occurrences in wetlands (red circles) on OSO land cover land used classification in the Écrins National Park. Each number correspond to a class described in Table 1. The observations mainly occur in urban wetlands in the valleys.</figcaption></figure>
# 4. Parameter models – Priors

We will now describe each prior used, i.e., the structure of the model components before updating them with information from observation data to obtain their posterior distributions. Each of the parameters $\beta_i$ or $\gamma_i$ associated with the fixed covariate effects has a very weakly informative Gaussian prior distribution centered at zero and with a precision hyperparameter $\tau_\beta = 10^{-3}$. For each month $m$, the spatial field $W_i^{(m)}$ is a priori represented using a 2-dimensional Gaussian random field based on the SPDE approach, where the Matérn covariance function has a fixed shape parameter $\nu^{(m)} = 1$, while the effective range, $\rho_i^{(m)}$, and the standard deviation, $\sigma_i^{(m)}$, are endowed with so-called penalized complexity priors (PC priors, Simpson et al., 2017): a priori, the probability of having range less than $\rho_0^{(m)} = 1$ is set to $\alpha_\rho^{(m)} = 0.5$, and the probability of having standard deviation larger than $\sigma_0^{(m)} = 1$ is set to $\alpha_\sigma^{(m)} = 0.5$. Since the spatial field $W_i^{(m)}$ is replicated month by month, a cyclic AR(1) correlation structure (Box et al., 2015) is applied such that the value of $W_i^{(m)}$ is informed by the values of $W_i^{(m-1)}$ and $W_i^{(m+1)}$, where $m - 1 = 12$ for $m = 1$ and $m + 1 = 12$ for $m = 12$. Regarding each period $p$, the inter-annual random effects $f_i(p)$ have independent and identically distributed Gaussian prior distributions with a precision hyperparameter $\tau_f = 10^{-3}$.

# References

Box G.E.P., Gwilym M.J., Reinsel G.C., Ljung G.M. (2015). Time Series Analysis : Forecasting and Control. 5th edition. Hoboken, New Jersey: John Wiley and Sons Inc. 712 pp.

Simpson D., Rue H., Riebler A., Martins T.G., Sørbye S.H. (2017). Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors. Statistical Science 32, 1–28.
https://doi.org/10.1214/16-STS576.
## 5. Temporal and spatial cross-validation scenarios

**Table 3** – Table of temporal cross-validation scenarios. The total number of space-time cells with at least one occurrence in the target-group is 32,486. For each scenario, half of these cells of the corresponding period of four years and their associated occurrences are removed.

<table><thead><tr><td>Scenario</td><td>Period</td><td>Removed cells</td><td>Percentage</td></tr></thead><tbody><tr><td>13</td><td>1994-1997</td><td>1782</td><td>5.5%</td></tr><tr><td>14</td><td>1998-2001</td><td>1274</td><td>3.9%</td></tr><tr><td>15</td><td>2002-2005</td><td>1512</td><td>4.7%</td></tr><tr><td>16</td><td>2006-2009</td><td>1525</td><td>4.7%</td></tr><tr><td>17</td><td>2010-2013</td><td>2670</td><td>8.2%</td></tr><tr><td>18</td><td>2014-2017</td><td>4477</td><td>13.8%</td></tr><tr><td>19</td><td>2018-2021</td><td>3001</td><td>9.2%</td></tr></tbody></table>

**Table 4** – Table of spatial cross-validation scenarios. The total number of space-time cells with at least one occurrence in the target-group is 32,486. For each scenario, the cells of a given city and their associated occurrences are removed.

<table><thead><tr><td>Scenario</td><td>Period</td><td>Removed cells</td><td>Percentage</td></tr></thead><tbody><tr><td>20</td><td>Saint-Christophe-en-Oisans</td><td>4407</td><td>13.6%</td></tr><tr><td>21</td><td>Valjouffrey</td><td>2473</td><td>7.6%</td></tr><tr><td>22</td><td>La Chapelle-en-Valgaudémar</td><td>2700</td><td>8.3%</td></tr><tr><td>23</td><td>Orcières</td><td>2167</td><td>6.7%</td></tr><tr><td>24</td><td>Châteauroux-les-Alpes</td><td>1839</td><td>5.7%</td></tr><tr><td>25</td><td>L'Argentiére-la-Bessée</td><td>1345</td><td>4.1%</td></tr><tr><td>26</td><td>La Grave</td><td>2656</td><td>8.2%</td></tr></tbody></table>
## 6. Habitat preferences

<table><caption>**Table 5** – Mean estimated parameters on spatial covariates : PCA Axis 1 to 3. Significant parameters based on the 95% interval of credibility are in bold.</caption><thead><tr><th>Species (Occurrences)</th><th>High Elevation</th><th>Coniferous forests</th><th>Grasslands</th></tr></thead><tbody><tr><td>Acanthis flammea (464)</td><td>0.03</td><td>0.1</td><td>0.39</td></tr><tr><td>Acrocephalus palustris (138)</td><td>-0.39</td><td>0.25</td><td>0.45</td></tr><tr><td>Aegithalos caudatus (1599)</td><td>-0.29</td><td>0.15</td><td>-0.29</td></tr><tr><td>Alauda arvensis (1068)</td><td>0.98</td><td>-0.54</td><td>0.73</td></tr><tr><td>Anthus spinoletta (1934)</td><td>1</td><td>-0.57</td><td>0.65</td></tr><tr><td>Anthus trivialis (1292)</td><td>-0.13</td><td>0.53</td><td>0.47</td></tr><tr><td>Bombycilla garrulus (146)</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>Carduelis carduelis (1549)</td><td>-0.21</td><td>-0.12</td><td>0.24</td></tr><tr><td>Carduelis citrinella (730)</td><td>0.49</td><td>0.04</td><td>0.33</td></tr><tr><td>Certhia brachydactyla (449)</td><td>-0.62</td><td>0.07</td><td>-0.35</td></tr><tr><td>Certhia familiaris (984)</td><td>-0.09</td><td>0.71</td><td>-0.27</td></tr><tr><td>Chloris chloris (668)</td><td>-0.36</td><td>-0.4</td><td>-0.11</td></tr><tr><td>Cinclus cinclus (1887)</td><td>-0.67</td><td>0.11</td><td>-0.05</td></tr><tr><td>Coccothraustes coccothraustes (486)</td><td>-0.18</td><td>-0.18</td><td>-0.11</td></tr><tr><td>Corvus corone (1745)</td><td>-0.21</td><td>-0.04</td><td>0.24</td></tr><tr><td>Corvus corone cornix (187)</td><td>-0.66</td><td>0.08</td><td>1.02</td></tr><tr><td>Corvus monedula (152)</td><td>-0.43</td><td>-0.05</td><td>-0.71</td></tr><tr><td>Cyanistes caeruleus (1892)</td><td>-0.27</td><td>0.01</td><td>-0.17</td></tr><tr><td>Delichon urbicum (734)</td><td>0.21</td><td>-0.13</td><td>-0.04</td></tr><tr><td>Emberiza cia (1548)</td><td>0.1</td><td>0.01</td><td>0.05</td></tr><tr><td>Emberiza cirlus (202)</td><td>-0.69</td><td>0.16</td><td>0.11</td></tr><tr><td>Emberiza citrinella (1734)</td><td>-0.21</td><td>-0.07</td><td>0.39</td></tr><tr><td>Emberiza hortulana (134)</td><td>0.02</td><td>0.24</td><td>0.91</td></tr><tr><td>Erithacus rubecula (3354)</td><td>-0.35</td><td>0.26</td><td>-0.18</td></tr><tr><td>Ficedula hypoleuca (273)</td><td>-0.55</td><td>0.13</td><td>0.05</td></tr><tr><td>Fringilla coelebs (5447)</td><td>-0.26</td><td>0.3</td><td>0.03</td></tr><tr><td>Fringilla montifringilla (551)</td><td>-0.22</td><td>-0.01</td><td>0.04</td></tr><tr><td>Garrulus glandarius (2201)</td><td>-0.37</td><td>0.31</td><td>-0.13</td></tr><tr><td>Hirundo rustica (642)</td><td>-0.38</td><td>-0.11</td><td>0.38</td></tr><tr><td>Lanius collurio (834)</td><td>-0.24</td><td>-0.14</td><td>0.35</td></tr><tr><td>Linaria cannabina (1369)</td><td>0.47</td><td>-0.31</td><td>0.55</td></tr><tr><td>Lophophanes cristatus (2160)</td><td>-0.16</td><td>0.63</td><td>-0.17</td></tr><tr><td>Loxia curvirostra (1890)</td><td>0.22</td><td>0.83</td><td>0.34</td></tr><tr><td>Lullula arborea (321)</td><td>-0.15</td><td>0.48</td><td>0.94</td></tr><tr><td>Luscinia megarhynchos (376)</td><td>-0.46</td><td>0.12</td><td>0.29</td></tr><tr><td>Monticola saxatilis (566)</td><td>1.16</td><td>-0.35</td><td>0.47</td></tr><tr><td>Montifringilla nivalis (999)</td><td>1.48</td><td>-1.13</td><td>0.44</td></tr><tr><td>Motacilla alba (1088)</td><td>-0.55</td><td>-0.22</td><td>0.22</td></tr><tr><td>Motacilla cinerea (697)</td><td>-0.37</td><td>-0.05</td><td>-0.08</td></tr></tbody></table>
<table><tr><td>Motacilla flava (144)</td><td>-0.24</td><td><b>-0.59</b></td><td><b>0.37</b></td></tr><tr><td>Nucifraga caryocatactes (2275)</td><td><b>0.39</b></td><td><b>0.49</b></td><td>0.02</td></tr><tr><td>Oenanthe oenanthe (2263)</td><td><b>0.94</b></td><td><b>-0.68</b></td><td><b>0.66</b></td></tr><tr><td>Parus major (3165)</td><td><b>-0.31</b></td><td>0</td><td><b>-0.16</b></td></tr><tr><td>Passer domesticus (947)</td><td><b>-0.45</b></td><td><b>-0.44</b></td><td>-0.09</td></tr><tr><td>Passer montanus (209)</td><td><b>-0.9</b></td><td>0.22</td><td><b>0.8</b></td></tr><tr><td>Periparus ater (3595)</td><td><b>-0.18</b></td><td><b>0.45</b></td><td><b>-0.21</b></td></tr><tr><td>Petronia petronia (197)</td><td><b>-0.48</b></td><td><b>-0.42</b></td><td>0.24</td></tr><tr><td>Phoenicurus ochruros (4183)</td><td><b>0.48</b></td><td><b>-0.43</b></td><td><b>0.11</b></td></tr><tr><td>Phoenicurus phoenicurus (668)</td><td><b>-0.27</b></td><td><b>-0.24</b></td><td>0.07</td></tr><tr><td>Phylloscopus bonelli (1668)</td><td><b>-0.17</b></td><td><b>0.17</b></td><td><b>-0.28</b></td></tr><tr><td>Phylloscopus collybita (2944)</td><td><b>-0.27</b></td><td><b>0.3</b></td><td>-0.03</td></tr><tr><td>Pica pica (1061)</td><td><b>-0.32</b></td><td><b>-0.18</b></td><td><b>0.24</b></td></tr><tr><td>Poecile montanus (2433)</td><td>0.07</td><td><b>0.7</b></td><td>-0.01</td></tr><tr><td>Poecile palustris (774)</td><td><b>-0.35</b></td><td>0.09</td><td><b>-0.21</b></td></tr><tr><td>Prunella collaris (1775)</td><td><b>1.2</b></td><td><b>-0.78</b></td><td>0.05</td></tr><tr><td>Prunella modularis (1386)</td><td><b>0.2</b></td><td><b>0.29</b></td><td><b>0.24</b></td></tr><tr><td>Ptyonoprogne rupestris (1695)</td><td><b>0.26</b></td><td><b>-0.11</b></td><td><b>-0.12</b></td></tr><tr><td>Pyrrhula pyrrhula (1816)</td><td><b>-0.32</b></td><td><b>0.44</b></td><td><b>-0.29</b></td></tr><tr><td>Regulus ignicapilla (686)</td><td><b>-0.37</b></td><td><b>0.74</b></td><td><b>-0.33</b></td></tr><tr><td>Regulus regulus (1030)</td><td><b>-0.27</b></td><td><b>0.79</b></td><td><b>-0.22</b></td></tr><tr><td>Saxicola rubetra (944)</td><td><b>-0.2</b></td><td><b>-0.15</b></td><td><b>0.74</b></td></tr><tr><td>Saxicola rubicola (313)</td><td>-0.14</td><td><b>-0.13</b></td><td><b>0.75</b></td></tr><tr><td>Serinus serinus (736)</td><td><b>-0.46</b></td><td><b>-0.09</b></td><td>-0.01</td></tr><tr><td>Sitta europaea (1841)</td><td><b>-0.47</b></td><td>0.03</td><td><b>-0.51</b></td></tr><tr><td>Spinus spinus (848)</td><td><b>-0.07</b></td><td><b>0.18</b></td><td><b>0.17</b></td></tr><tr><td>Sturnus vulgaris (353)</td><td><b>-0.56</b></td><td><b>-0.29</b></td><td><b>0.46</b></td></tr><tr><td>Sylvia atricapilla (1904)</td><td><b>-0.46</b></td><td><b>0.13</b></td><td><b>-0.24</b></td></tr><tr><td>Sylvia borin (576)</td><td><b>-0.49</b></td><td><b>0.14</b></td><td><b>0.23</b></td></tr><tr><td>Sylvia communis (143)</td><td><b>-0.39</b></td><td>0.36</td><td><b>0.71</b></td></tr><tr><td>Sylvia curruca (753)</td><td><b>0.26</b></td><td><b>0.35</b></td><td><b>0.28</b></td></tr><tr><td>Tichodroma muraria (1318)</td><td><b>0.86</b></td><td><b>-0.27</b></td><td><b>-0.46</b></td></tr><tr><td>Troglodytes troglodytes (2903)</td><td><b>-0.27</b></td><td><b>0.47</b></td><td>0.01</td></tr><tr><td>Turdus merula (2938)</td><td><b>-0.39</b></td><td><b>0.18</b></td><td>-0.04</td></tr><tr><td>Turdus philomelos (1209)</td><td><b>-0.37</b></td><td><b>0.56</b></td><td>-0.06</td></tr><tr><td>Turdus pilaris (1270)</td><td>0.04</td><td><b>0.2</b></td><td><b>0.29</b></td></tr><tr><td>Turdus torquatus (1665)</td><td><b>0.51</b></td><td><b>0.64</b></td><td><b>0.36</b></td></tr><tr><td>Turdus viscivorus (2585)</td><td><b>-0.11</b></td><td><b>0.43</b></td><td><b>0.21</b></td></tr></table>
<table><caption>**Table 6** – Mean estimated parameters on spatial covariates : PCA Axis 4 to 6. Significant parameters based on the 95% interval of credibility are in bold.</caption><thead><tr><th>Species (Occurrences)</th><th>Urban wetlands</th><th>Woody Heaths</th><th>Deciduous forests</th></tr></thead><tbody><tr><td>Acanthis flammea (464)</td><td>-0.26</td><td>-0.09</td><td>0.13</td></tr><tr><td>Acrocephalus palustris (138)</td><td>-0.24</td><td>-0.04</td><td>-0.25</td></tr><tr><td>Aegithalos caudatus (1599)</td><td>0</td><td><strong>0.07</strong></td><td><strong>-0.2</strong></td></tr><tr><td>Alauda arvensis (1068)</td><td>0.11</td><td><strong>-0.34</strong></td><td>-0.01</td></tr><tr><td>Anthus spinoletta (1934)</td><td>0.02</td><td><strong>0.1</strong></td><td>0.01</td></tr><tr><td>Anthus trivialis (1292)</td><td><strong>-0.43</strong></td><td><strong>-0.37</strong></td><td><strong>0.14</strong></td></tr><tr><td>Bombycilla garrulus (146)</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>Carduelis carduelis (1549)</td><td><strong>0.15</strong></td><td>-0.01</td><td><strong>-0.08</strong></td></tr><tr><td>Carduelis citrinella (730)</td><td>-0.02</td><td>0.04</td><td>0.03</td></tr><tr><td>Certhia brachydactyla (449)</td><td>0.05</td><td>0.03</td><td><strong>-0.22</strong></td></tr><tr><td>Certhia familiaris (984)</td><td>-0.27</td><td><strong>-0.56</strong></td><td>0.08</td></tr><tr><td>Chloris chloris (668)</td><td><strong>0.23</strong></td><td>-0.01</td><td><strong>-0.17</strong></td></tr><tr><td>Cinclus cinclus (1887)</td><td><strong>0.14</strong></td><td>0.04</td><td><strong>-0.16</strong></td></tr><tr><td>Coccothraustes coccothraustes (486)</td><td><strong>0.15</strong></td><td>0</td><td><strong>-0.1</strong></td></tr><tr><td>Corvus corone (1745)</td><td>-0.02</td><td>-0.02</td><td><strong>-0.11</strong></td></tr><tr><td>Corvus corone cornix (187)</td><td>-0.09</td><td>0.08</td><td><strong>-0.41</strong></td></tr><tr><td>Corvus monedula (152)</td><td>-0.1</td><td><strong>0.67</strong></td><td>-0.23</td></tr><tr><td>Cyanistes caeruleus (1892)</td><td>-0.06</td><td>0.05</td><td><strong>-0.12</strong></td></tr><tr><td>Delichon urbicum (734)</td><td>-0.02</td><td><strong>0.26</strong></td><td>0.04</td></tr><tr><td>Emberiza cia (1548)</td><td><strong>-0.2</strong></td><td><strong>0.16</strong></td><td>-0.05</td></tr><tr><td>Emberiza cirlus (202)</td><td>-0.19</td><td>0</td><td>-0.06</td></tr><tr><td>Emberiza citrinella (1734)</td><td>-0.11</td><td><strong>-0.18</strong></td><td><strong>-0.18</strong></td></tr><tr><td>Emberiza hortulana (134)</td><td>-0.29</td><td>0.09</td><td>-0.1</td></tr><tr><td>Erithacus rubecula (3354)</td><td><strong>-0.1</strong></td><td>0.01</td><td>-0.01</td></tr><tr><td>Ficedula hypoleuca (273)</td><td>0.09</td><td>0.09</td><td><strong>-0.13</strong></td></tr><tr><td>Fringilla coelebs (5447)</td><td><strong>-0.12</strong></td><td><strong>-0.1</strong></td><td><strong>0.07</strong></td></tr><tr><td>Fringilla montifringilla (551)</td><td>-0.09</td><td>-0.02</td><td>-0.08</td></tr><tr><td>Garrulus glandarius (2201)</td><td><strong>-0.16</strong></td><td>-0.03</td><td><strong>-0.06</strong></td></tr><tr><td>Hirundo rustica (642)</td><td><strong>0.27</strong></td><td>-0.07</td><td><strong>-0.14</strong></td></tr><tr><td>Lanius collurio (834)</td><td>-0.14</td><td><strong>-0.2</strong></td><td><strong>-0.14</strong></td></tr><tr><td>Linaria cannabina (1369)</td><td>0.08</td><td>0.04</td><td><strong>0.12</strong></td></tr><tr><td>Lophophanes cristatus (2160)</td><td>-0.11</td><td><strong>-0.15</strong></td><td><strong>0.24</strong></td></tr><tr><td>Loxia curvirostra (1890)</td><td>-0.23</td><td><strong>-0.46</strong></td><td><strong>0.51</strong></td></tr><tr><td>Lullula arborea (321)</td><td><strong>-0.21</strong></td><td>-0.05</td><td>0.01</td></tr><tr><td>Luscinia megarhynchos (376)</td><td>-0.14</td><td>-0.05</td><td>-0.11</td></tr><tr><td>Monticola saxatilis (566)</td><td><strong>-0.61</strong></td><td>0.09</td><td><strong>0.35</strong></td></tr><tr><td>Montifringilla nivalis (999)</td><td><strong>-0.18</strong></td><td>-0.09</td><td>0.01</td></tr><tr><td>Motacilla alba (1088)</td><td><strong>0.34</strong></td><td>-0.02</td><td><strong>-0.16</strong></td></tr><tr><td>Motacilla cinerea (697)</td><td><strong>0.25</strong></td><td>0.07</td><td><strong>-0.2</strong></td></tr><tr><td>Motacilla flava (144)</td><td><strong>0.83</strong></td><td><strong>-0.72</strong></td><td>-0.02</td></tr><tr><td>Nucifraga caryocatactes (2275)</td><td><strong>0.15</strong></td><td>-0.01</td><td><strong>0.13</strong></td></tr></tbody></table>
<table><tr><td>Oenanthe oenanthe (2263)</td><td><strong>-0.15</strong></td><td>0.06</td><td>0.03</td></tr><tr><td>Parus major (3165)</td><td>-0.03</td><td><strong>0.11</strong></td><td><strong>-0.17</strong></td></tr><tr><td>Passer domesticus (947)</td><td><strong>0.38</strong></td><td>-0.3</td><td>-0.05</td></tr><tr><td>Passer montanus (209)</td><td>0.07</td><td>-0.27</td><td>0</td></tr><tr><td>Periparus ater (3595)</td><td><strong>-0.23</strong></td><td><strong>-0.24</strong></td><td><strong>0.26</strong></td></tr><tr><td>Petronia petronia (197)</td><td>0.1</td><td>-0.32</td><td><strong>-0.33</strong></td></tr><tr><td>Phoenicurus ochruros (4183)</td><td><strong>-0.15</strong></td><td><strong>0.11</strong></td><td>-0.04</td></tr><tr><td>Phoenicurus phoenicurus (668)</td><td><strong>0.25</strong></td><td>0</td><td><strong>-0.16</strong></td></tr><tr><td>Phylloscopus bonelli (1668)</td><td>-0.06</td><td><strong>0.14</strong></td><td><strong>-0.15</strong></td></tr><tr><td>Phylloscopus collybita (2944)</td><td>0</td><td><strong>0.1</strong></td><td><strong>-0.14</strong></td></tr><tr><td>Pica pica (1061)</td><td>-0.02</td><td><strong>-0.12</strong></td><td><strong>-0.12</strong></td></tr><tr><td>Poecile montanus (2433)</td><td>-0.17</td><td><strong>-0.37</strong></td><td><strong>0.19</strong></td></tr><tr><td>Poecile palustris (774)</td><td>-0.05</td><td>0.05</td><td><strong>-0.11</strong></td></tr><tr><td>Prunella collaris (1775)</td><td>-0.08</td><td><strong>0.18</strong></td><td>0.09</td></tr><tr><td>Prunella modularis (1386)</td><td><strong>-0.17</strong></td><td>0.02</td><td><strong>0.22</strong></td></tr><tr><td>Ptyonoprogne rupestris (1695)</td><td>-0.07</td><td><strong>0.32</strong></td><td>-0.03</td></tr><tr><td>Pyrrhula pyrrhula (1816)</td><td>-0.1</td><td><strong>-0.24</strong></td><td><strong>0.14</strong></td></tr><tr><td>Regulus ignicapilla (686)</td><td><strong>-0.56</strong></td><td><strong>-0.42</strong></td><td><strong>0.35</strong></td></tr><tr><td>Regulus regulus (1030)</td><td>-0.22</td><td><strong>-0.24</strong></td><td><strong>0.3</strong></td></tr><tr><td>Saxicola rubetra (944)</td><td>-0.07</td><td><strong>-0.17</strong></td><td><strong>-0.22</strong></td></tr><tr><td>Saxicola rubicola (313)</td><td><strong>-0.21</strong></td><td>-0.1</td><td><strong>-0.26</strong></td></tr><tr><td>Serinus serinus (736)</td><td><strong>0.14</strong></td><td><strong>0.11</strong></td><td>-0.08</td></tr><tr><td>Sitta europaea (1841)</td><td>0</td><td><strong>-0.08</strong></td><td><strong>-0.16</strong></td></tr><tr><td>Spinus spinus (848)</td><td>0.09</td><td>0.04</td><td><strong>0.13</strong></td></tr><tr><td>Sturnus vulgaris (353)</td><td><strong>0.32</strong></td><td><strong>-0.67</strong></td><td>0.03</td></tr><tr><td>Sylvia atricapilla (1904)</td><td>0.01</td><td><strong>0.16</strong></td><td><strong>-0.09</strong></td></tr><tr><td>Sylvia borin (576)</td><td>0</td><td>0.06</td><td><strong>-0.21</strong></td></tr><tr><td>Sylvia communis (143)</td><td><strong>-0.3</strong></td><td>0.09</td><td>-0.1</td></tr><tr><td>Sylvia curruca (753)</td><td><strong>-0.49</strong></td><td>-0.03</td><td><strong>0.2</strong></td></tr><tr><td>Tichodroma muraria (1318)</td><td><strong>-0.16</strong></td><td><strong>0.4</strong></td><td><strong>0.13</strong></td></tr><tr><td>Troglodytes troglodytes (2903)</td><td><strong>-0.23</strong></td><td><strong>-0.13</strong></td><td>0.05</td></tr><tr><td>Turdus merula (2938)</td><td><strong>-0.07</strong></td><td><strong>0.05</strong></td><td><strong>-0.05</strong></td></tr><tr><td>Turdus philomelos (1209)</td><td>0.03</td><td><strong>-0.26</strong></td><td>0.05</td></tr><tr><td>Turdus pilaris (1270)</td><td><strong>-0.23</strong></td><td><strong>-0.15</strong></td><td>0.06</td></tr><tr><td>Turdus torquatus (1665)</td><td><strong>-0.6</strong></td><td>-0.3</td><td><strong>0.41</strong></td></tr><tr><td>Turdus viscivorus (2585)</td><td><strong>-0.22</strong></td><td><strong>-0.19</strong></td><td><strong>0.08</strong></td></tr></table>
<table><caption>**Table 7** – Mean estimated parameters on spatio-temporal covariates : anomalies of maximum temperature and anomalies of precipitation. Significant parameters based on the 95% interval of credibility are in bold.</caption><thead><tr><th>Species (Occurrences)</th><th>Anomalies of Tmax</th><th>Anomalies of Precipitation</th></tr></thead><tbody><tr><td>Acanthis flammea (464)</td><td>-0.002</td><td>-0.0215</td></tr><tr><td>Acrocephalus palustris (138)</td><td>-0.0044</td><td><strong>0.2126</strong></td></tr><tr><td>Aegithalos caudatus (1599)</td><td>0.001</td><td><strong>-0.0633</strong></td></tr><tr><td>Alauda arvensis (1068)</td><td>0.003</td><td><strong>0.1266</strong></td></tr><tr><td>Anthus spinoletta (1934)</td><td>0.0026</td><td>0.0376</td></tr><tr><td>Anthus trivialis (1292)</td><td>0.0019</td><td>0.0237</td></tr><tr><td>Bombycilla garrulus (146)</td><td>NA</td><td>NA</td></tr><tr><td>Carduelis carduelis (1549)</td><td>6e-04</td><td>-0.0013</td></tr><tr><td>Carduelis citrinella (730)</td><td>-0.0019</td><td><strong>-0.0407</strong></td></tr><tr><td>Certhia brachydactyla (449)</td><td>-0.001</td><td>-0.0068</td></tr><tr><td>Certhia familiaris (984)</td><td>0.0012</td><td><strong>-0.0321</strong></td></tr><tr><td>Chloris chloris (668)</td><td><strong>0.0071</strong></td><td><strong>0.1904</strong></td></tr><tr><td>Cinclus cinclus (1887)</td><td>-5e-04</td><td><strong>-0.0223</strong></td></tr><tr><td>Coccothraustes coccothraustes (486)</td><td>-0.0031</td><td><strong>-0.272</strong></td></tr><tr><td>Corvus corone (1745)</td><td>-0.0016</td><td><strong>-0.0549</strong></td></tr><tr><td>Corvus corone cornix (187)</td><td>-0.0052</td><td><strong>-0.0928</strong></td></tr><tr><td>Corvus monedula (152)</td><td>-0.0061</td><td>-5e-04</td></tr><tr><td>Cyanistes caeruleus (1892)</td><td>1e-04</td><td><strong>-0.0823</strong></td></tr><tr><td>Delichon urbicum (734)</td><td>5e-04</td><td><strong>0.1054</strong></td></tr><tr><td>Emberiza cia (1548)</td><td>-5e-04</td><td>0.0086</td></tr><tr><td>Emberiza cirlus (202)</td><td>0.0024</td><td><strong>0.0605</strong></td></tr><tr><td>Emberiza citrinella (1734)</td><td>-0.0012</td><td>-0.0128</td></tr><tr><td>Emberiza hortulana (134)</td><td>0.0097</td><td><strong>0.2859</strong></td></tr><tr><td>Erithacus rubecula (3354)</td><td>-0.0013</td><td>-0.0195</td></tr><tr><td>Ficedula hypoleuca (273)</td><td>-5e-04</td><td>0.0073</td></tr><tr><td>Fringilla coelebs (5447)</td><td>-0.0011</td><td>-0.0236</td></tr><tr><td>Fringilla montifringilla (551)</td><td>0.0031</td><td><strong>-0.2194</strong></td></tr><tr><td>Garrulus glandarius (2201)</td><td>2e-04</td><td>0.0014</td></tr><tr><td>Hirundo rustica (642)</td><td><strong>0.0073</strong></td><td>0.0208</td></tr><tr><td>Lanius collurio (834)</td><td><strong>0.0065</strong></td><td><strong>0.2762</strong></td></tr><tr><td>Linaria cannabina (1369)</td><td>-0.0021</td><td><strong>0.0514</strong></td></tr><tr><td>Lophophanes cristatus (2160)</td><td>-2e-04</td><td>-0.0279</td></tr><tr><td>Loxia curvirostra (1890)</td><td>0.0013</td><td><strong>-0.0489</strong></td></tr><tr><td>Lullula arborea (321)</td><td>3e-04</td><td>0.0426</td></tr><tr><td>Luscinia megarhynchos (376)</td><td><strong>-0.0094</strong></td><td>0.1224</td></tr><tr><td>Monticola saxatilis (566)</td><td>0.006</td><td><strong>0.1575</strong></td></tr><tr><td>Montifringilla nivalis (999)</td><td>0.0014</td><td><strong>-0.1305</strong></td></tr><tr><td>Motacilla alba (1088)</td><td>-0.0024</td><td>0.0163</td></tr><tr><td>Motacilla cinerea (697)</td><td>-0.0032</td><td><strong>0.0632</strong></td></tr><tr><td>Motacilla flava (144)</td><td>-0.0051</td><td>0.1618</td></tr></tbody></table>
<table><tr><td>Nucifraga caryocatactes (2275)</td><td>0</td><td>4e-04</td></tr><tr><td>Oenanthe oenanthe (2263)</td><td>-6e-04</td><td><strong>0.129</strong></td></tr><tr><td>Parus major (3165)</td><td>2e-04</td><td><strong>-0.0695</strong></td></tr><tr><td>Passer domesticus (947)</td><td>-0.001</td><td>-0.0128</td></tr><tr><td>Passer montanus (209)</td><td>-0.0014</td><td>-0.0342</td></tr><tr><td>Periparus ater (3595)</td><td>4e-04</td><td><strong>-0.0424</strong></td></tr><tr><td>Petronia petronia (197)</td><td>-0.0036</td><td><strong>0.0778</strong></td></tr><tr><td>Phoenicurus ochruros (4183)</td><td>-1e-04</td><td>0.0478</td></tr><tr><td>Phoenicurus phoenicurus (668)</td><td>-0.0015</td><td><strong>-0.0212</strong></td></tr><tr><td>Phylloscopus bonelli (1668)</td><td>-4e-04</td><td>0.0714</td></tr><tr><td>Phylloscopus collybita (2944)</td><td><strong>-0.0031</strong></td><td>0.0105</td></tr><tr><td>Pica pica (1061)</td><td><strong>0.0026</strong></td><td><strong>-0.0426</strong></td></tr><tr><td>Poecile montanus (2433)</td><td>0.0013</td><td><strong>-0.0437</strong></td></tr><tr><td>Poecile palustris (774)</td><td><strong>0.0036</strong></td><td><strong>-0.1009</strong></td></tr><tr><td>Prunella collaris (1775)</td><td>0.0014</td><td><strong>-0.0971</strong></td></tr><tr><td>Prunella modularis (1386)</td><td><strong>-0.0047</strong></td><td>-0.0339</td></tr><tr><td>Ptyonoprogne rupestris (1695)</td><td>0.001</td><td><strong>0.0746</strong></td></tr><tr><td>Pyrrhula pyrrhula (1816)</td><td>9e-04</td><td><strong>-0.0381</strong></td></tr><tr><td>Regulus ignicapilla (686)</td><td>-0.004</td><td>-0.0081</td></tr><tr><td>Regulus regulus (1030)</td><td>0.002</td><td><strong>-0.0646</strong></td></tr><tr><td>Saxicola rubetra (944)</td><td>0.0027</td><td><strong>0.2046</strong></td></tr><tr><td>Saxicola rubicola (313)</td><td>0.0021</td><td><strong>0.1124</strong></td></tr><tr><td>Serinus serinus (736)</td><td>-7e-04</td><td><strong>0.1231</strong></td></tr><tr><td>Sitta europaea (1841)</td><td><strong>-0.0027</strong></td><td><strong>-0.0289</strong></td></tr><tr><td>Spinus spinus (848)</td><td>1e-04</td><td><strong>-0.1358</strong></td></tr><tr><td>Sturnus vulgaris (353)</td><td>0.005</td><td>0.0386</td></tr><tr><td>Sylvia atricapilla (1904)</td><td><strong>-0.0034</strong></td><td>0.0182</td></tr><tr><td>Sylvia borin (576)</td><td>-0.0041</td><td><strong>0.1405</strong></td></tr><tr><td>Sylvia communis (143)</td><td><strong>0.0193</strong></td><td><strong>0.2733</strong></td></tr><tr><td>Sylvia curruca (753)</td><td>0.0049</td><td><strong>0.1386</strong></td></tr><tr><td>Tichodroma muraria (1318)</td><td>-0.0019</td><td><strong>-0.0576</strong></td></tr><tr><td>Troglodytes troglodytes (2903)</td><td>-1e-04</td><td>0.0162</td></tr><tr><td>Turdus merula (2938)</td><td>6e-04</td><td><strong>-0.0262</strong></td></tr><tr><td>Turdus philomelos (1209)</td><td>-5e-04</td><td>-0.0198</td></tr><tr><td>Turdus pilaris (1270)</td><td>-0.0012</td><td><strong>-0.1565</strong></td></tr><tr><td>Turdus torquatus (1665)</td><td>-8e-04</td><td>-0.0144</td></tr><tr><td>Turdus viscivorus (2585)</td><td><strong>0.0014</strong></td><td><strong>-0.0779</strong></td></tr></table>
<table><caption>**Table 8** – List of species removed from the habitat preference clustering based on the estimated mean parameters associated with the 6 spatial covariates corresponding to the first 6 axes of the principal component analysis detailed in Appendix 3.</caption><thead><tr><td>Species (Occurrences)</td><td>Reason</td></tr></thead><tbody><tr><td>Acanthis flammea (464)</td><td>1 significant covariate</td></tr><tr><td>Acrocephalus palustris (138)</td><td>2 significant covariates</td></tr><tr><td>Bombycilla garrulus (146)</td><td>131 occurrences in 2005</td></tr><tr><td>Carduelis citrinella (730)</td><td>2 significant covariates</td></tr><tr><td>Cinclus cinclus (1887)</td><td>Poor AUCs</td></tr><tr><td>Corvus monedula (152)</td><td>1 significant covariate</td></tr><tr><td>Emberiza cirlus (202)</td><td>1 significant covariate</td></tr><tr><td>Emberiza hortulana (134)</td><td>1 significant covariate</td></tr><tr><td>Ficedula hypoleuca (273)</td><td>2 significant covariates</td></tr><tr><td>Fringilla montifringilla (551)</td><td>1 significant covariate</td></tr><tr><td>Luscinia megarhynchos (376)</td><td>1 significant covariate</td></tr><tr><td>Motacilla cinerea (697)</td><td>Poor AUCs</td></tr><tr><td>Passer montanus (209)</td><td>2 significant covariates</td></tr><tr><td>Sylvia communis (143)</td><td>2 significant covariates</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure 14** – a) Map of elevation in the Écrins National Park on a 500m×500m regular grid. b) Map of clusters with the highest occurrences on a 500m×500m regular grid. The highest number of occurrences is computed from the 63 classified species: the 77th ones with more than 100 occurrences, excluding the species listed in Table 8. Cluster 1 is associated with forest species, Cluster 2 with high-elevation species, and Cluster 3 with species in the valleys.</figcaption></figure>