Combining Citizen Science Data and Satellite Descriptors of Ecosystem Functioning to Monitor the Abundance of a Migratory Bird during the Non-Breeding Season

: Migratory birds are particularly exposed to habitat changes in their breeding and non-breeding grounds. Remote sensing technologies offer an excellent opportunity to monitor species’ habitats from space at unprecedented spatiotemporal scales. We analyzed if remotely sensed ecosystem functioning attributes (EFAs) adequately predict the spatiotemporal variation of the Woodcock’s ( Scolopax rusticola ) relative abundance in southwest Europe, during autumn migration and wintering periods. We used data gathered from Woodcock monitoring through citizen science (N = 355,654 hunting trips) between 2009 and 2018. We computed a comprehensive set of EFAs on a weekly basis from three MODIS satellite products: enhanced vegetation index (EVI), tasseled cap transformation (TCT), and land surface temperature (LST). We developed generalized linear mixed models to explore the predictive power of EFAs on Woodcock’s abundance during the non-breeding season. Results showed that Woodcock abundance is correlated with spatiotemporal dynamics in primary productivity (measured through the EVI), water cycle dynamics (wetness component of TCT), and surface energy balance (LST) in both periods. Our ﬁndings underline the potential of combining citizen science and remote sensing data to monitor migratory birds throughout their life cycles—an issue of critical importance to ensure adequate habitat management in the non-breeding areas.


Introduction
In the face of rapid global environmental changes, it is increasingly important to understand the main factors affecting species' habitats throughout their annual life cycle [1]. The status of habitats in wintering and migratory seasons becomes critical for species that spend part of the year in different regions. These species are particularly exposed to changes in conditions and resources. Several aspects of ecosystem functioning relevant for species have been severely affected by climate and land-use change over recent decades [2]. As the annual life cycle of migratory species involves movements between different locations, it becomes more complex to evaluate the influence of interannual environmental changes in these species [3,4]. Managers and decision makers need to be updated with the best information available on the status of both breeding and non-breeding populations to protect them effectively.
Advances in geographic information systems (GIS) and remote sensing for Earth observation have significatively improved our ability to understand, monitor, and forecast species distribution changes over recent decades [5,6]. Combining both technologies provides a broad spectrum of environmental data at unprecedented spatiotemporal scales offering increasing amounts of information about the entire planet [7][8][9]. In addition, ecological modeling has also received greater attention in recent years due to its wide range of applications for the study and conservation of biodiversity [5,10]. Species distribution models (SDMs), based on the assessment of species' ecological niches, link georeferenced observations of a biotic response variable such as the occurrence or abundance of species with several environmental predictors through a wide range of statistical or machine learning algorithms. These techniques are powerful tools for conservation biology, allowing predictions of habitat suitability and the probability of species occurrence with relatively high accuracy [11] even when complete information of their entire distributional range is not available [6,12]. Therefore, SDMs are very useful to fill knowledge gaps about the geographic distribution of species. SDMs have been used to identify priority areas for conservation, assess environmental impacts, and predict future environmental changes, making it possible to manage critical resources or conditions that affect species habitat [8,13]. Additionally, species abundance models (SAMs) can provide critical information for species monitoring and conservation management. However, these models have received less attention from the scientific community, possibly due to difficulties in obtaining abundance data related to sampling protocols' costs and specificity [14]. In general, species occurrence data are more accessible than population abundance data because they only require recording the presence or evidence of at least one individual rather than estimates of the absolute or relative number of individuals [15]. However, SAMs have the advantage of being much more informative and capable of providing relevant information on the distribution of species and the size of populations, which can reflect the role of essential demographic and environmental factors [15][16][17].
A major challenge in developing SDMs and SAMs relates to the choice of environmental variables. Predictor variables should match the spatial resolution of the species data to capture habitat characteristics at finer scales, thus including a wider range of information that is ecologically relevant for the species [12,18,19]. Environmental predictor variables derived from remote sensing data are commonly used to model species distributions. However, satellite images have been used mainly in conservation biology to classify, describe, and map vegetation and habitats' structure [20]. A recent review points out that SDMs and SAMs have suffered from the lack of spatially explicit predictor variables capturing the species' habitats dynamics throughout their annual life cycle [21]. Nevertheless, emerging remote sensing technologies face these challenges and contribute to a new generation of distribution and abundance models [13,17]. New satellite products can improve SDMs' and SAMs' performance by providing essential information to predict species ranges.
An example of such satellite products is the ecosystem functioning attributes (EFAs), which are biophysical descriptors of ecosystem functioning that describe exchanges of matter and energy between the biota and the environment [12,22]. EFAs are calculated from satellite recorded time series and offer a more integrated and faster assessment of ecosystem responses to environmental factors and changes than macroclimatic databases or structural attributes (e.g., vegetation height and density, landscape composition, or spatial configuration) [23]. Additionally, since EFAs are remotely detected in a standardized and synoptical fashion, species habitats' spatial and temporal (seasonal and interannual) variability can be easily included in SAM or SDM workflows [22]. Despite these advances, they are still relatively unexplored for these purposes [13,24]. The possibility of including predictors of the seasonal and interannual habitat dynamics in SAMs can be an opportunity for studying migratory species whose annual life cycle develops in different areas of the globe with specific phenological events marking their departure, arrival, or permanence periods.
The Woodcock (Scolopax rusticola) is a migrant wader (Charadriiformes), distinguished from other scolopacids by its association with forest areas [25,26]. Environmental conditions influence the Woodcock's behavior, especially in the non-reproductive period (i.e., wintering period), when the distribution and abundance of its populations seem to be particularly affected by temperature and rainfall [27][28][29]. The Woodcock is sensitive to the thermal regime [27,29,30] and seasonal variations in soil moisture [31,32], as well as habitat features related to vegetation cover [33][34][35]. Landscape and forest habitats' heterogeneity can significantly affect Woodcock abundance since habitat requirements appear to vary with different stages of its annual life cycle [35]. The Woodcock's distribution is also strongly conditioned by food availability [31,36]. It has a diet specialized in soil microfauna: arthropods, annelids (mainly worms), and slugs [36]. Low temperatures decrease its capacity to regulate body temperature and affect food accessibility since frozen water makes it harder to penetrate deeper into the soil to access the worms [37]. Under these conditions, the Woodcock shows a more irregular distribution and variable density, with stochastic events shaping its population dynamics [27]. When rainfall is strong, and the temperature is mild, food is more abundant, and the Woodcock distributes more evenly in relatively low densities.
The Woodcock is a species with great game interest in Europe, where it is estimated that 3 to 4 million individuals are hunted annually [38]. There, Woodcock hunter associations from different European countries seek to involve their affiliates in managing and conserving the species, encouraging them to collect data from their activity. Thus, a large volume of long-term data has been collected during the hunting season across Europe. This citizen science data source offers an excellent opportunity to study the factors affecting the Woodcock's distribution and abundance during the winter season across a wide geographic area. As referred by Runge et al. [39], management and conservation actions for migratory birds need to be coordinated across different regions, habitat types, seasons, and jurisdictions. Additionally, there is an increasing consensus on the need for sustainable hunting, which depends on a deeper understanding of the factors driving the game species' interannual population fluctuations [40].
The present study aims to understand better how ecosystem functioning influences the population abundance of migratory birds throughout the non-breeding season. To do so, we focused on the Woodcock population across southwestern Europe (Portugal, Spain, and France) during the autumn and winter seasons (migrating and wintering birds), using hunting data.
We hypothesized that remotely sensed variables depicting several dimensions of ecosystem functioning (i.e., EFAs expressing vegetation and edaphoclimatic conditions) and their seasonal dynamics could explain the spatiotemporal variation of the Woodcock's abundance. In particular, we took advantage of a long-term time series obtained from citizen science (N = 355,654 hunting trips) gathered between 2009 and 2018. We computed a comprehensive set of zonal statistical parameters for three EFAs, for 8-day periods, from MODIS satellite products: enhanced vegetation index (EVI), tasseled cap transformation (TCT), and land surface temperature (LST). Finally, we analyzed to what extent these remotely sensed indicators of primary productivity, water, and energy balance predict  [44]; study area (Portugal, Spain, and France) delimited by the black line. In each country, the respective regions of the Nomenclature of Territorial Units for Statistics of level 3 of 2016, NUTS-III (data adapted from [45] are delimited) and spatial variation (by NUTS-III) of the number of Woodcock hunting trips analyzed, carried out by hunters (during the respective hunting seasons, between 2009 and 2018), are represented in five classes by the number quantiles of journeys, i.e., each class represents 20% of the total number of observations. NUTS-III white color-no data.
During the breeding season, the Woodcock selects habitat mosaics with specific characteristics such as deciduous forest habitats [46,47], mixed forest, or conifers [25,48,49]. During the wintering period, it shows less specific habitat requirements, selecting different habitats such as hedged forests during the day, for refuge, and meadows at night, for feeding [25,31,50], where it can find high invertebrate biomass, particularly annelids [26,51].
The current Woodcock breeding population in the Western Palearctic is estimated to be between 10 and 26 million individuals, most of which spend the winter in western and  [44]; study area (Portugal, Spain, and France) delimited by the black line. In each country, the respective regions of the Nomenclature of Territorial Units for Statistics of level 3 of 2016, NUTS-III (data adapted from [45] are delimited) and spatial variation (by NUTS-III) of the number of Woodcock hunting trips analyzed, carried out by hunters (during the respective hunting seasons, between 2009 and 2018), are represented in five classes by the number quantiles of journeys, i.e., each class represents 20% of the total number of observations. NUTS-III white color-no data.
During the breeding season, the Woodcock selects habitat mosaics with specific characteristics such as deciduous forest habitats [46,47], mixed forest, or conifers [25,48,49]. During the wintering period, it shows less specific habitat requirements, selecting different habitats such as hedged forests during the day, for refuge, and meadows at night, for feeding [25,31,50], where it can find high invertebrate biomass, particularly annelids [26,51].
The current Woodcock breeding population in the Western Palearctic is estimated to be between 10 and 26 million individuals, most of which spend the winter in western and southern Europe and northern Africa [52,53]. However, these estimates are mainly supported by specialists' opinions and not on objective data collected in the field [54]. The Woodcock is a solitary, elusive, and cryptic species, limiting classic survey techniques usually applied to assess bird abundance [34,35,55,56]. Thus, citizen science data obtained from hunting activity offer an important opportunity to study this species.
The Woodcock appears in Part A of Annex II and Part B of Annex III of the "Birds Directive", which means it can be hunted in the geographical area where the directive is applicable. However, the effects of hunting on their populations are still poorly assessed, with evidence of an additive effect on mortality and causing the use of more extensive areas during the day [33,57]. In addition to hunting, land-use changes may also cause some disturbance in the Woodcock's annual life cycle. Although there are some studies about its ecology, these are mainly concentrated in the breeding period and/or in restricted areas of its distribution, such as those carried out using ringing data [28], stable isotopes [43], and telemetry [49,58].
Its current global conservation status is evaluated as Least Concern; the only available regional trend assessment concerns Europe, where the population seems to be stable [59].

Study Area and Hunting Data
The study area covers central and southwestern Europe (France, Spain, and Portugal), which represents one of the main wintering areas for Woodcock [25,43] (see Figure 1). Due to the long hunting traditions, population monitoring programs were developed to support a more sustainable hunting activity in these countries. In Europe, Woodcock hunter associations from different countries are gathered in the Federation of Western Palearctic Woodcock Hunter Association (Féderation des Associations Nationales des Bécassiers du Palearctique Occidental (FANBPO)), which facilitates their cooperation and communication.
We analyzed data collected between 2009 and 2018 by Woodcock hunters, members of the Club National des Bécassiers, the Club de Cazadores de Becada, and the Associação Nacional de Caçadores de Galinholas, during Woodcock hunting trips in mainland France, Spain, and Portugal, respectively ( Figure 1). Woodcock hunting seasons occur from September to February of the following year, with differences in the beginning date between these countries (France: 14 September; Spain: 8 October; Portugal: 1 November), and in the number of hunting days and bag limits. Only data collected from hunting with pointing dogs were considered in the analyses. Hunters registered the date, location, duration, and number of Woodcocks seen for each hunting trip.
A total of 355,654 hunting trip records were analyzed: 326,519 (~92%) from France, 25,746 (~7%) from Spain, and 3389 (~1%) from Portugal ( Figure 1), totaling 536,873 contacts with Woodcock during 1,154,222 h of hunting. We estimated a relative abundance index corresponding to the number of different Woodcock observed per hour in each hunting trip (ICA or ICA1 [60][61][62]). This abundance index is strongly related to another relative abundance index obtained from night ringing sessions (IAN or NIA) and is considered a good indicator to evaluate winter abundance variation [60][61][62].
Since the location of the hunting trips was provided with variable detail among the three countries, the relative abundance index values were aggregated at level 3 of NUTS classification (NUTS-III), which corresponds to departments in France, provinces in Spain, and intermunicipal entities in Portugal. The average area (±standard deviation) of these units is 6320.4 ± 8054.7 km 2 , 8576 ± 5698.6 km 2 , and 3676 ± 2020.7 km 2 , for France (n = 88), Spain (n = 41), and Portugal (n = 23), respectively.

Remote Sensed Ecosystem Functioning Variables
To evaluate how ecosystem functioning influences the Woodcock's relative abundance, we computed a comprehensive set of zonal statistical parameters for three ecosystem functioning attributes from MODIS satellite data products: enhanced vegetation index (EVI), tasseled cap transformation (TCT), and land surface temperature (LST) [63,64] (Table 1). These variables are suited to characterizing ecosystem functioning in a repeatable, continuous, and standardized way, being related to primary productivity, water, and energy balance [12,17,22]. Table 1. List of variables estimated from the Terra/MODIS satellite imaging products. For the list of acronyms of the remote sensed ecosystem functioning variables, see Table A3.

Components of Ecosystem Functioning
Ecosystem Functioning Attributes-EFAs

MODIS Product(s)/Pixel Size and Time Frequency Zonal Statistical Parameters (NUTS-III)
Carbon cycle (primary productivity) Water cycle (water in soil/vegetation)
EVI value can vary between −1 and 1, with the lowest values being associated with artificial cover types and the highest values, closer to one, are related to higher levels of biomass, vegetation cover, leaf area index, or photosynthetic activity. The EVI was obtained from the MODIS MOD13Q1 product with 250 m spatial resolution and 8-day temporal resolution.
TCT allows, similarly to principal component analysis, performance of a linear transformation to satellite spectral data to highlight specific aspects of the land surface. These components are related to "brightness" (which approximately translates to the albedo coefficient), "greenness" (as a proxy of vegetation biomass), or the "wetness" component (as a proxy of water content in the soil and vegetation). The coefficients used in the transformation to the wetness component were derived from Zhang et al. [66] (Table A1). The TCTwet index (i.e., TCT transformation for the wetness component) value can vary between −1 and 1 and indicates, in a relative way, the amount of water in the soil or vegetation, being able to, under continuous observation, show the availability and water balance in a given region. TCTwet was calculated from the MODIS MOD09A1 product with 500 m spatial resolution and 8-day temporal resolution.
To calculate LST from the MODIS product MOD11A2 (with 1 km spatial resolution and 8-day temporal resolution), a multiplicative conversion factor of 0.02 (units in Kelvin) was applied.
Using metadata of MODIS images, the pixels marked as clouds were removed to avoid spurious values. The data collection was carried out through original images for each variable (Table 1) and by applying zonal statistics (i.e., aggregating all the pixels included in a given area) for each NUTS-III region (see the available list in Table A2 in the Appendix A). Given the large size of the spatial units (NUTS-III), in addition to the average value for each unit, we also computed the minimum and the maximum value (Table 1) to characterize the extreme conditions at the NUTS level. We also calculated each variable's range to characterize the spatial heterogeneity inside each NUTS-III region (the larger the range, the greater the heterogeneity within the NUTS-III region). Thus, four zonal statistic measures were used to characterize the data distribution at the NUTS level: minimum, average, maximum, and range. These values were estimated for 8-day periods from daily images extracted between September and February of the following year, covering the annual hunting seasons from 2009 (September) to 2018 (February). The EVI values at the NUTS-III level changes with the latitude and 8-day period. All image reprocessing and calculation processes were performed in the Google Earth Engine cloud-based platform [67]. We used "dplyr" and "tidyr" R packages [68,69] for data management and preprocessing.

Variable Selection and Ranking
The data on the relative abundance of Woodcock were grouped by week to match as close as possible with the time scale of the available remote sensing data ( Figure 2). Both data types were divided into two periods, corresponding to two phases of the species' annual life cycle, following [54,70]: Remote Sens. 2022, 14, x FOR PEER REVIEW 7 of 24 latitude and 8-day period. All image reprocessing and calculation processes were performed in the Google Earth Engine cloud-based platform [67]. We used "dplyr" and "tidyr" R packages [68,69] for data management and preprocessing.

Variable Selection and Ranking
The data on the relative abundance of Woodcock were grouped by week to match as close as possible with the time scale of the available remote sensing data ( Figure 2). Both data types were divided into two periods, corresponding to two phases of the species' annual life cycle, following [54,70]: 1. Autumnal migration, from early September (first week of the month, week 35) to mid-December (second week of the month, week 50). 2. Wintering season, from mid-December (week 51) to late February (week 9). All remote sensing variables were standardized (2), considering the respective mean and standard deviation values using the "scale" function available in the "arm" R package [71].

1.
Autumnal migration, from early September (first week of the month, week 35) to mid-December (second week of the month, week 50).

2.
Wintering season, from mid-December (week 51) to late February (week 9). All remote sensing variables were standardized (2), considering the respective mean and standard deviation values using the "scale" function available in the "arm" R package [71].
To avoid multicollinearity problems, we calculated Pearson's correlation coefficient between each pair of variables (see the list and respective acronyms in Table A3). The correlations were represented graphically ( Figure A1). The variance inflation factor (VIF) was also calculated with the "usdm" R package [72] to estimate how much a regression coefficient's variance is inflated by multicollinearity. The variables with Pearson's correlation coefficient (with absolute value) greater than 0.7 and a VIF value greater than three were excluded from further analyses [73,74].
To explore the effects of remotely sensed ecosystem functioning variables (fixed factors) on the spatiotemporal variation of the Woodcock's relative abundance (dependent variable), mixed generalized linear models (GLMMs) were developed separately for the migration and wintering period. To avoid outliers that could affect model inference, we removed abundance values greater than the interquartile distance multiplied by 1.5 [75,76]. The "country" and the "week" nested with the annual "hunting season" were considered as random effects to account for repeated measurements within and across units of time ("week") and space ("countries"). Each country has a hunting period due to climatic and hunting legal specificities, so a "country" level effect must be accounted for. The "week" was also included as a random effect nested with the annual "hunting season" because abundance data have been collected every week of the migratory and wintering season over the years. Since the response variable (i.e., the relative abundance index) is a continuous variable, our models were fitted assuming a normal distribution and using the "link" identity function. Previous analyses with Poisson and negative binomial error distributions and zero-inflated models were finally ruled out due to their lower explanatory power (conditional R 2 value ranged between 0.07 and 0.12, and marginal R 2 between 0.04 and 0.07). The analyses were performed with the "glmmTMB" R package [77]. Model performance was assessed using the "check_model" function, available in the "performance" R package [78]. This complementary analysis confirms various model assumptions: normality of residuals, normality of random effects, heteroscedasticity, homogeneity of variance, and multicollinearity (see Figures A2 and A3 in the Appendix A).
In addition, a multi-model inference approach was performed for both migratory and wintering seasons, running models from all possible combinations of environmental variables (predictors) to verify the importance of each one of them through the "dredge" function of the "MuMIn" R package [79]. Finally, only those models with a delta Akaike information criterion (AIC) value of less than four were considered relevant [80]. The variables were considered significant for p-values < 0.05. Šidák's correction [81] was computed to adjust p-values for multiple comparisons.

Patterns of Spatiotemporal Variation in Woodcock Abundance
Our results confirm the expected variation in the Woodcock's relative abundance over the hunting season ( Figure 3). The relative abundance index initially increases over the migration period until reaching a maximum peak, followed by a relatively stable phase that coincides with the wintering season. There is a slight time lag in this pattern between countries; initially, the abundance increases in France, Spain, and Portugal. Regarding the spatial variation of Woodcock abundance, the species is more abundant in the north of France, the southwest of Spain, and in the center and the south of Portugal, for both wintering and migration periods (Figure 4). Some sampling units have no abundance data available, particularly in southern Spain and northeastern France. Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 24       [85]. NUTS-III white color-no data.

Autumnal Migration
The GLMM fitted for the autumn migration period, including all environmental variables, presented a moderate explanatory ability, with a conditional R 2 value of 0.33 and a marginal R 2 value of 0.09. Six of the seven variables considered for analysis had a significant effect (Table 2). Those with the most significant effect were the medium LST (Z = −13.766; p < 0.001) and the maximum EVI (Z = −10.512; p < 0.001), both with a negative effect on species abundance ( Table 2). The EVI amplitude and medium value (Z = 9.541; p < 0.001; Z = 8.982; p < 0.001) also stand out, having a positive effect on abundance. The medium value of the TCTwet component also had a significant and negative effect on the abundance of the species (Z = −8.220; p < 0.001) and the maximum LST had a positive effect (Z = 4.659; p < 0.001). All these significant variables are present in all models considered relevant in AIC multimodel inference (i.e., delta AIC < 4; Table 2). LST amplitude was the only descriptor of ecosystem functioning that was not significant (p = 0.225; Table 2). Table 2. Generalized linear mixed model (GLMM) results for the autumn migration and wintering periods. Estimate of regression coefficients results from an averaged model obtained from those with delta AIC < 4. The importance of each ecosystem functioning variable (last column) is represented by the number of times that each variable appears contained in the explanatory models (delta AIC < 4). EFAs: EVI-enhanced vegetation index; TCTwet-tasseled cap transformation-water component; LST-land surface temperature; respective maximum (max), minimum (min), average (avg), and range (range) values. Šidák's correction [81] was computed to adjust p-values for multiple comparisons. Significant results, at p < 0.05, are shown in bold and underlined if significant after correction for seven comparisons.

Wintering
The GLMM built for the wintering period, including all environmental variables, also presented a moderate explanatory capacity, with a conditional R 2 value of 0.16 and a marginal R 2 value of 0.09. Five of the seven variables considered for analysis had a significant effect (Table 2). Those with the most significant effect were the average EVI (Z = 11.692; p < 0.001) and the maximum EVI (Z = −10.774; p < 0.001), the first with a positive effect and the second with a negative effect on the species abundance (Table 2). EVI amplitude (Z = 7.723; p < 0.001) showed a positive effect, while the average and maximum value of the LST (Z = −5.277; p < 0.001; Z = 5.130; p < 0.001, respectively) had for the first a negative effect and the second a positive effect. The maximum LST component positively affected the species relative abundance (Z = 5.130; p < 0.001). These variables are shown in all the models considered, supported by AIC-based multimodel inference (delta AIC < 4; Table 2). The average value of the TCTwet component (Z = −2.129; p = 0.333) and the amplitude of LST (Z = 1.350; p = 0.177; Table 2) are the only descriptors of ecosystem functioning that were not significant after adjustment for multiple comparisons.

Ecosystem Functioning Attributes and Woodcock Abundance
Our results confirm that remotely sensed habitat descriptors (namely EFAs and their seasonal dynamics) have a significant explicative ability for the spatiotemporal patterns of abundance of migratory birds. In the autumn migration period, the maximum EVI and average LST (descriptors of primary productivity and energy balance, respectively) correlated significantly and negatively with Woodcock abundance. The maximum EVI stands out in the wintering period and the average EVI has a positive effect. Additionally, significantly correlated, in both periods, are the EVI amplitude values.
The highest maximum EVI values signal those areas with the highest primary productivity [65], such as forest areas with dense vegetation. These characteristics seem to be associated with the Woodcock's lower abundance in the two studied periods ( Table 2). This result is consistent with previous knowledge about the species' habitat requirements indicating a preference for mixed areas of forests and hedges during the day and fields and meadows at night, rather than areas of homogeneous and dense forest [31,33,49]. On the other hand, higher EVI amplitude values indicate a greater heterogeneity in the vegetation cover. This situation will benefit the Woodcock, which, in these periods of its annual life cycle, needs open areas for night feeding (i.e., pastures, natural or artificial meadows, agricultural fields) and forest and bush areas for shelter during the day [31,33,49]. Our results are aligned with recent studies that support a relationship between primary productivity and the migration of long-distance migratory birds, pointing out that their migratory movement is influenced by seasonal changes in the landscape's resource availability [86].
The maximum value of land surface temperature (LST) showed a significant and positive correlation with abundance but was slightly higher in the wintering period ( Table 2). The average value of LST is negatively related to Woodcock abundance. This effect may indicate that the maximum temperature at the surface (reached during the day) is important. These low temperatures imply a greater difficulty for birds regarding body temperature regulation with subsequent increases in energy requirements and, therefore, food intake [30]. Péron et al. [37] demonstrated that Woodcock make great movements in winter conditioned by the air temperature. On the other hand, it is also known that the low temperatures reached during the wintering period hamper accessibility to their primary food source, earthworms, obtained by probing the soil [37]. Thus, the maximum value of the temperature at the surface, reached during the daytime, can allow favorable conditions for accessing food, at least in a certain period of the day. Finally, LST amplitude does not seem to influence Woodcock abundance in both periods significantly. Our results reinforce the interest in satellite products of surface temperature, more related to the conditions experienced by individuals on the ground in critical periods of their life cycle [87].
The average value of the TCT wetness component (TCwet; a descriptor of the water cycle in the soil/vegetation) made a significant contribution, negatively affecting the species' abundance in the migration period (Table 2). This result may indicate that this species avoids areas with high levels of waterlogging (flooding), detrimental to obtaining food. The Woodcock, unlike other waders highly dependent on water dynamics such as the Common Snipe (Gallinago gallinago) [88], is not frequently seen in habitats with these characteristics, such as peat bogs and swamps [25].

Advantages, Limitations, and Future Perspectives
Currently, most studies indicate that habitat loss and degradation (such as the degradation and fragmentation of forest mosaics or the intensification of agricultural practices) and climate change are significant threats to species and ecosystems [89,90]. For birds, climate change influences their phenology, population dynamics, abundance, and distribution [91][92][93].
However, although these changes affect bird communities, it is unclear which factors are most prevalent in each case and which pathways. Different factors are challenging to analyze independently. Ecosystem functioning attributes (EFAs) are multipurpose descriptors that offer a more integrative and rapid assessment of responses to environmental changes. Our results allow us to verify that the abundance of Woodcock in autumn migration and winter periods is affected by several dimensions linked to primary productivity dynamics, the water cycle (water in the soil and vegetation), and ecosystems' energy balance. These results are consistent with those obtained in previous studies, highlighting the validity of these remotely sensed indicators for species distribution and abundance assessments [12,17,22,23,88,94,95]. Several studies have already demonstrated the advantages of using models calibrated with functional variables over those based on models calibrated with exclusively climatic variables since the former allowed for capture of the joint effect of changes in climate and land cover/use on the availability of habitat for a broad group of threatened species [22,95]. Other studies about the variation in threatened plant species abundance and various bird species (especially for migrant, forest specialist species) also confirmed a greater predictive power of models calibrated with functional variables related to productivity and energy balance compared to models based on macroclimatic data or landscape composition variables [17,95].
Migratory bird species, especially those considered for hunting, such as the Woodcock, represent a greater challenge for management and conservation policies. It is vital to consider territorial continuity (beyond administrative countries' boundaries) to analyze which factors affect species and an adequate definition of population management goals. In this sense, our study considered a broader and more representative area of the distribution in autumn migration and wintering of Woodcock than previous studies (e.g., [31,96]).
The use of data collected from hunting activity to study game species' population dynamics is often hampered by the temporal (between hunting seasons) and the spatial variation in their quantity and quality [96][97][98][99]. The data collected often do not include metadata that allow quantifying the hunting effort (e.g., time dedicated to hunting on each journey or the area covered). Our relative abundance index considers an effort variable (sampling effort, i.e., the time spent hunting). In addition, it reports the number of birds seen and not the number of birds hunted, which reduces the influence of certain factors such as legal limitations due to the catch limits in each country. Thus, previous works have recognized this index as a useful tool [60][61][62].
On the other hand, hunting data present some limitations that must be considered. NUTS-III regions were chosen due to the location of the hunting trips that were provided with variable detail among the three countries. The large size of the NUTS-III regions leads to high environmental heterogeneity within and among them. Being aware of this limitation, we think that this dataset still provides vital information since it consists of more than 350,000 records, collected over a long time interval (nine hunting seasons) and a large geographical extent.
Still, the spatial resolution used for analyses prevents a more refined characterization of the species habitat and its ecological requirements, partly explaining the low values of marginal R 2 . An increase in the spatial resolution at which the abundance data are collected will improve future studies since alternative satellites such as Landsat or Sentinel (e.g., [100]) would allow calculating EFAs at a finer spatial resolution (up to 30 m and 10 m, respectively) although with less temporal resolution than MODIS. In addition, the application of predictive modeling techniques based on artificial intelligence such as machine learning algorithms (e.g., random forest or artificial neural networks) can also improve our ability to predict the species' abundance from space.

Final Considerations
The European Union (EU) seeks to curb biodiversity loss, a task that needs updated knowledge of the general conservation status of species (including wild birds) and habitats of community interest under the 1992 "Habitats Directive" (92/43/EEC) and the 1979 "Birds Directive" (2009/147/EC). However, the information on birds is still considered insufficient, as is the case in Iberian territory [10].
The present work used a new explanatory analysis for assessing the variation of Woodcock abundance in autumn and winter at large spatial and temporal scales, based on ecosystem functioning attributes obtained by Earth observation satellites. Additionally, it demonstrates the importance of citizen science programs in collecting data essential for monitoring biodiversity.
The approach developed can be applied to other taxa since it establishes strong relationships between the species abundance and essential ecosystem functions and vital ecological processes closely related to the water and carbon cycles and the energy balance. Considering the increasing availability of remote sensing products, these data provide a promising opportunity for integrating ecosystem dynamics into habitat monitoring to support management and species conservation policies. This is especially relevant given the current context of global environmental changes and the new environmental, agricultural, and nature restoration policies to be implemented in Europe in the coming decades in areas critically important for Woodcock.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
We would like to thank all the members and collaborators of Club National des Bécassiers (France), Club de Cazadores de Becada (Spain), and Associação Nacional de Caçadores de Galinholas (Portugal), that made this study possible through a citizen science network.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.        Model performance for the autumn migration period, obtained with the R package "performance" [78]. Figure A2. Model performance for the autumn migration period, obtained with the R package "performance" [78].