Highlights

  • The fitness of the Golden Eagle is influenced by ecosystem functioning.

  • Surface energy budget disruptions also affected the fitness of this apex predator.

  • Satellite-based ecosystem functions as indicators of the species fitness.

Introduction

Anthropogenic global change is severely affecting ecosystem functioning and biodiversity globally (Flombaum and others 2017; Rillig and others 2019; Díaz and others 2020). A growing body of work has shown how climate warming, land-use change and their joint effects have been—directly or indirectly—affecting biodiversity and ecosystem functioning, with difficult-to-predict side effects on their complex and dynamic relationships (Carnicer and others 2011; Dib and others 2020). A better understanding of the ecosystem functioning-biodiversity relationships at the different trophic levels is essential to help managers adequately respond to the threats that global change poses to biodiversity and ecosystem health (Thébault and Loreau 2006). It is also critical to develop cost-effective monitoring and early warning systems to proactively respond to these threats with effective management strategies (Cabello and others 2012; Anderson 2018; Aragón and others 2019).

Apex predators are at the top of food chain, playing a critical role in the population dynamics of lower trophic levels (Wallach and others 2015). The interactive effects of global climate warming, regional land-use changes and altered disturbance regimes can critically affect the population dynamics of these keystone species, with subsequent impacts on other trophic levels through predator–prey interactions (Guiden and others 2019). Until now, most ecological studies on apex predators are based upon modelling approaches supported by macroclimate data (obtained from interpolations of meteorological stations) and static habitat features (see for example, López-López and others 2007; Tapia and others 2007). Despite the increasing availability and quality of macroclimate datasets, these data cannot accurately represent the microclimatic conditions often driven by topographic and local vegetation patterns—failing to capture critical features of their ecological niche (Amiri and others 2020). These limitations highlight the need of using more proximal descriptors of the ‘real’ habitat conditions—in the broadest sense—that the species is experiencing at the ground level.

Ecosystem functional attributes offer several advantages to monitor species populations since they are integrative descriptors of the environmental change (see for example, Alcaraz-Segura and others 2017; Arenas-Castro and others 2019; Regos and others 2020), being closely related to the processes directly affecting food chains via trophic interactions (for example, energy budget and hydric balance or nutrient cycles) (Oksanen and Oksanen 2000; Letnic and others 2018). It is widely accepted that energy balance and net primary productivity affects local species richness because the more energy that is available, the more biomass per unit area that can be supported (Phillips and others 2008; Brasil and others 2019). More biomass enables more individuals to coexist in an area, potentially resulting in more energy available for reproduction, lowering extinction rates (Connell and Orians 1964). Thus, the non-inclusion of critical ecosystem processes such as primary productivity or surface energy balance can undermine our capacity to monitor these keystone species. For instance, the temporal variability in primary productivity was recently found to affect habitat-quality fluctuations and meta-population dynamics of consumer species by affecting the energy levels available as food resources (Fernández and others 2016; Letnic and others 2018). Sea ice phenology and primary productivity pulses shape breeding success in Artic seabirds (Ramírez and others 2017). Polar bear foraging behaviour was also found to vary temporally and spatially with observed variation in primary productivity and the health of prey and other species within their Arctic food web (Pagano and others 2018; Rode and others 2018). These are examples that illustrate the need of more functional perspective when assessing the ecological fitness and population dynamics of species at upper trophic levels (Fernández and others 2016; Aragón and others 2019).

The remote sensing discipline has largely evolved over the last decades, currently allowing monitoring of and assessing key components of ecosystem functioning (that is, the biogeochemical flow of energy and matter within ecosystems) (Lovett and others 2006; Jax 2010; Van den Broeke and others 2011; Lausch and others 2016). Despite the increasing availability of remotely sensed descriptors of ecosystem functions, their potential contribution to the monitoring of species at upper trophic levels remains poorly studied. Among others, land surface temperature and albedo, and seasonal productivity and phenology are good descriptors of ecosystem functions—being relatively easy to obtain from remote sensing images (Metz and others 2014; Zhao and others 2018b; Aragón and others 2019).

Land Surface Temperature (LST) is a key descriptor of the radiative energy budget of the Earth’s surface (Hulley and others 2019). The intra- and interannual variation of this descriptor is often used as proxy for sensible heat dynamics (Li and others 2013; Metz and others 2014). Satellite-derived LST is typically retrieved by estimating the surface-emitted radiance that is obtained by atmospherically correcting the at-sensor radiance. LST measures the emission of thermal radiance from the land surface where the incoming solar energy interacts with and heats the ground, or the surface of the canopy in vegetated areas (Oyler and others 2019). LST is therefore a good indicator of energy partitioning at the land surface-atmosphere boundary and especially sensitive to variations in land-cover type, vegetation density or soil moisture (Hulley and others 2019). In fact, LST is largely shaped by climatic gradients and vegetation productivity patterns, being also locally influenced by topographic and microclimatic conditions (Oyler and others 2019). In addition, both natural and anthropogenic disturbances (for example, wildfires) can be another important source of spatiotemporal variation due to abrupt changes in LST patterns.

Albedo, defined as the reflected solar radiation to the incoming solar radiation, is another key component of ecosystem functioning related to surface energy balance. It can be measured from space by integrating information retrieved from specific multi-band sensors based on the anisotropy properties defined by bidirectional reflectance distribution functions. Albedo varies over large areas due to variations in land surface properties such as land cover type, green vegetation cover, surface roughness or soil moisture. It can be therefore affected by human-induced land-use changes and climate change feedback loops (Hu and others 2019). In addition, fire disturbance can abruptly change land surface with long-standing effects on albedo (Dintwe and others 2017).

Both LST and Albedo are therefore critical descriptors of species habitat conditions (of both hunting and breeding areas), being partially influenced by vegetation productivity—another key component of ecosystem functioning. Vegetation productivity is often measured from spectral vegetation indices such as the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI), among others. Spatiotemporal variations in phenology and vegetation productivity are driven by climate (including extreme events; Ma and others 2015), land use (Paruelo and others 2001a) and their combined effects (Zhang and others 2014). Forest disturbance such as wildfires, insect outbreaks or defoliators has the potential to greatly impact vegetation productivity but differently depending on the frequency and severity of the disturbance event, and the ecoregion (Sparks and others 2018). In addition, primary consumers play a critical role in the spatiotemporal dynamics of vegetation productivity (via consumption), indirectly affected by the top-down effects of predators on herbivore populations (Letnic and others 2018). Simultaneously, temporal fluctuations in primary productivity were found to affect biomass of plants and consumers (bottom-up effects), with side-effects on predator populations via trophic cascades (Letnic and Ripple 2017; Letnic and others 2018)—relationship that depends on the spatial and temporal gradients in primary productivity (see the exploitation ecosystems hypothesis; Oksanen and Oksanen 2000).

This work aims to improve our understanding on how ecosystem functioning influence keystone species at upper trophic levels by leveraging the great potential of remote sensing to track matter and energy fluxes from space. We expect that remotely sensed ecosystem functional attributes (derived from remote sensing products such as LST, Albedo and vegetation productivity) can—at least—partially explain the ecological fitness of an apex predator—the Golden Eagle. To do so, we took advantage of a long-term time series database of the reproductive success of the Golden Eagle over a 17-year period across a bioclimatic gradient. We computed a comprehensive database of ecosystem functional attributes from three MODIS (Moderate Resolution Imaging Spectroradiometer) satellite-products related to the carbon cycle, heat dynamics and radiative balance. We also assessed possible time-lag in the response of the Golden Eagle to fire, a critical disruptor of the surface energy budget in our region.

Material and Methods

Study Region and Target Species

The study area is Galicia, region located in the Northwestern Spain (c. 29,575 km2) (Figure 1). The climate is transitional between Atlantic and Mediterranean (Rodríguez Guitián and Ramil-Rego 2007; Rodríguez-Lado and others 2016). Average annual rainfall ranges between 700 and 2000 mm, and annual mean temperature between 8.5 and 15 °C (Rodríguez-Lado and others Rodríguez-Lado). Altitudes are maximal (up to 2071 m.a.s.l.) in the eastern mountains and minimal in the coast and fluvial valleys. Forest areas cover over 70% of the study region, being currently one of the most densely forested areas of Spain. In mountain areas, extensive traditional agropastoral systems have been preserved until today, although they have gradually been abandoned since the mid twentieth century, being recolonized by natural forests and fast-growing tree plantations (Regos and others 2015; Salaverri and others 2018). The study area is relatively sparsely populated, with some large areas having less than 10 inhabitants/km2 in the eastern mountainous areas. Human-driven fire regime (including negligence and arson causes) exposes the study area to an unusually high frequency of wildfires (Chas-Amil and others 2010). Between 1968 and 2016, more than 256,000 fires took place in Galicia, burning more than 1,931,000 ha (c.a. 63% of the territory) (Regos 2018).

Figure 1
figure 1

Location of the study area in Europe. Hue-Saturation-Value (HSV colour space) image computed from the combination of MODIS products depicting phenology (date of annual minimum Albedo (ALB), in the hue circular component), seasonality or intra-annual variation (standard deviation of Land Surface Temperature (LST), saturation component), and primary productivity (annual average of Enhanced Vegetation Index (EVI), value component). Polygons represent the mountain areas wherein historical and current breeding territories are located.

Full size image

The Golden Eagle is an apex predator, placed at the upper level of the trophic cascade. This raptor can be described as both generalist and opportunistic predator, feeding on a wide variety of species across its worldwide distribution range (Watson 2010; Bedrosian and others 2017). Typically, this species feeds on medium-sized mammals ranging from 0.5 to 4.0 kg, although it expands its diet when preferred prey species are scarce (Steenhof and Kochert 1988; Watson 2010). An apex predator such us the Golden Eagle have high food and energy requirements so its performance and ecological fitness, measured by the reproductive success, strongly depends on interannual fluctuations of prey availability—mostly primary consumers—which in turn depend on surface energy balance dynamics and ecosystem functionality (Fernández and others 2016; Preston and others 2017).

The study area constitutes the Northwestern border of the Golden Eagle’s distribution range in the Iberian Peninsula, where species is considered currently ‘Near Threatened’ (Madroño and others 2005). In Galicia, the Golden Eagle has a population of 14 confirmed breeding pairs, with other 9 pairs shared with neighbouring regions of northern Portugal and Spain (Gil-Carrera and others 2016). Suitable areas for nesting in the study site are characterized by being located in mountain ranges between 255 and 1443 m.a.s.l., river canyons of steep slopes, Mediterranean climatic influence and low human pressure areas (Lado and Tapia 2012). The main conservation problems are associated with habitat limitations such as inadequate food-supply, poor availability of nesting sites, human disturbance in breeding areas such as electrocution and land-use changes (Madroño and others 2005; Tapia and others 2007). Moreover, this small population inhabits a transitional bioclimatic area relatively isolated from other nuclei of the Iberian Peninsula (Tapia and others 2007; Lado and Tapia 2012), and an increase in mortality can have a major impact on its demographic parameters (Whitfield and others 2004).

Field Surveys

Annual field survey campaigns (foot and road searching from a 4 × 4 vehicle) have been conducted between 2001 and 2017 on areas suitable for nesting in order to monitor the known breeding territories and detect new ones (see (Tapia and others 2007, 2009; Gil and Tapia 2009) for details on the field survey and inventory of the former breeding territories). The distribution of the Golden Eagle in the study area is well known, being the rate of pseudo-absences low (N = 397 observations). The final dataset includes 59 historical and currently occupied breeding territories in Galicia and the neighbouring regions (up to 10 km from the Galician border). These 59 breeding territories are located across different mountain systems with particular geomorphological and bioclimatic conditions (Figure 1). We characterized the home range of each breeding pair by considering radii of 5 km from the nest (that is, including the nest site, post-fledging and foraging areas; Tapia and Zuberogoitia 2018). Over the 17-year period, 168 breeding pairs were recorded in the 59 territories, although only 81 were confirmed. Some breeding pairs have changed their nests during the study period within the same breeding territory. In these cases, we selected the central point of nest location in each territory for the analysis.

Together with immigration and emigration, natality and mortality are the most important demographic parameters that determine the year-to-year trends in raptor local populations (Steenhof and Newton 2007). In the present study, we focused on the reproductive success for each territory and year of monitoring. Reproductive success is defined as the proportion of nesting and laying pairs that raise young to the age of fledging. We consider pairs to be successful when well-grown young are observed in the nest at some point prior to fledging. Studies that consider nests with young of any age to be successful tends to overestimate nest success because they fail to consider mortality that may occur late in the brood-rearing period (Steenhof and Newton 2007).

Remotely Sensed Descriptors of Energy Balance

Three MODIS (Moderate Resolution Imaging Spectroradiometer) satellite-products were selected to describe three dimensions of ecosystem functioning: (1) Enhanced Vegetation Index as a surrogate for the carbon cycle dynamics; (2) Land Surface Temperature as a surrogate of sensible heat dynamics; and (3) Albedo as a surrogate for the radiative balance (Paruelo and others 2001b; Cabello and others 2012). In addition, we included (4) Burnt Area, as a descriptor of fire disturbance—which can dramatically affect these three biophysical descriptors of ecosystem functions (Bowman and others 2009; Archibald and others 2018). The four MODIS selected products were:

  1. (1)

    Enhanced Vegetation Index (EVI; MOD13Q1.006; 16-Day L3 Global 250 m). EVI is an descriptor of carbon gains since it is known to be more reliable in both low and high vegetation cover situations than Normalized Difference Vegetation Index (NDVI), and resistant to both soil influences, canopy background signals, and atmospheric effects on vegetation index values (Potter and others 2007). EVI values ranged from -1 to 1, with healthy vegetation generally holding values between 0.20 and 0.80.

  2. (2)

    Land Surface Temperature (LST; MOD11A2.006; 8-Day L3 Global 1 km). LST is a good descriptor of the energy balance at the Earth’s surface, and one of the key parameters in the physics of land-surface processes from regional to global scales. In addition, LST is directly linked to the primary environmental regimes and to habitat quality attributes (for example, productivity, vegetation structure, land cover type; Cord and Rödder 2011). Temperatures (LST) ranged from -25 °C to 45 °C.

  3. (3)

    Albedo (ALB; MCD43A3.006; Daily L3 Global 500 m). ALB is surrogate for surface properties such as the extent and nature of the vegetation cover, and it is affected by the change of the land surface bio-physical factors such as vegetation, land surface temperature and soil moisture (Zhao and others 2018a). ALB values ranged from 0 to 1 (fresh snow and bare soil usually fall around 0.9).

  4. (4)

    Burnt Area (BA; Collection 6 MODIS Burned Area product). BA is a good descriptor of the surface energy budget disruptions since wildfires represent a critical exchange of energy and matter between the surface and the atmosphere via combustion (Bowman and others 2009; Archibald and others 2018). BA values ranged from 0 to 1 (values up to 1 for those areas totally affected by fire).

All native, and derived, MODIS predictors were computed for the 2001–2017 time period at 1-km resolution and then aggregated (mean values) at the 5-km buffer level by using the cloud-based Google Earth Engine computational platform (https://earthengine.google.com) (Gorelick and others 2017). For each of these 3 dimensions of ecosystem functioning (EVI, LST and ALB), we derived the following 7 summary metrics of their seasonal dynamics: annual mean (surrogate of annual total amount); annual maximum and minimum (descriptors of the annual extremes); seasonal standard deviation and coefficient of variation (descriptors of seasonality); and dates of maximum and minimum (descriptors of phenology) (see Table 1). The complete dataset included 21 descriptors of ecosystem functioning (7 metrics × 3 dimensions) as competing candidate predictors. Each of these 21 remotely sensed ecosystem functional attributes were then computed for each year to measure the interannual variability. In addition, we created variables describing a time-lag in the response of the Golden Eagle to fire. To do so, we estimated the burnt area of each year in each 5-km buffer around the nest. Then, lagged versions of the time series were created by shifting the year in which the fire event took place (from 1 to 5 years after a fire event).

Table 1 Remotely Sensed Metrics Computed From Each MODIS Product to Describe Ecosystem Functioning
Full size table

Model Fitting

To evaluate the effect of ecosystem functioning on reproductive success of the Golden Eagle, we performed logistic-exposure models of nest survival (Brown and others 2013) by using generalized linear mixed-effects models (GLMMs; Bolker and others 2009) with R packages ‘lmer4’ (Bates and others 2014).

These logistic-exposure nest survival models are similar to logistic regression models in that the response variable is binomial (breeding success or failure), but the link function is modified from the logit link to consider nest exposure days (Shaffer 2004)—101 days for the Golden eagle (see Brown and others 2013 and references therein). GLMMs were used to account for the hierarchical structure of the species dataset, whereby each underlying dataset covers a different breeding season (that is, year) in each mountain system (Bolker and others 2009; Harrison and others 2018). We included the remotely sensed descriptors of ecosystem functions as fixed factors and ‘year’ and ‘mountain system’ as nested random effects. Based on this random effect structure, the fixed effect component was then modified using the ‘dredge’ function available from the R package ‘MuMIn’ to run automated model selection with subsets of all (valid) combinations of explanatory variables (Barton 2016). We calculated each model delta AIC (ΔAIC) to verify the strength of evidence in favour of each model. Delta AIC values below 4 suggest substantial evidence for the model (Burnham and Anderson 2002; Arnold 2010). An average final model was computed from those with ΔAIC below 4 by using the ‘model.avg’ function available from the R package ‘MuMIn’. The importance of each predictor was obtained by adding the Akaike weights (Wi) to the models in which that variable is present (Burnham and Anderson 2002; Arnold 2010).

To avoid including highly correlated variables in model fitting, we removed remotely sensed variables with Pearson’s correlation coefficients above 0.7 or below − 0.7 and Variable Inflation Factor greater than 3. Thus, from the initial dataset of 21 remotely sensed ecosystem functional attributes, we only retained 12 non-correlated variables related with carbon cycle (4 variables from EVI time series), sensitive heat dynamics (4 variables from LST time series) and radiance balance (4 variables from Albedo time series) (see Figure S1 for pairwise correlations between explanatory variables).

We start by fitting GLMMs with the 4 non-correlated variables considered for each of the 3 dimensions of ecosystem functioning (LST, ALB and EVI; GLMMs 1, 2 and 3 in Table 2). From this first round of GLMMs, we selected a subset of the 4 most important remotely sensed ecosystem functional attributes according to the Akaike weights (ΣWi). We then fit a fourth GLMM with these 4 variables to assess the performance of the nest survival model fitted with the three components of ecosystem functioning (hereafter ‘integrative model’; GLMM 4 in Table 2). Although all models were fitted with the full dataset, their performance was evaluated by using a split-sample approach to produce predictions independent of the training data. Each model was fitted using 70% of the data and evaluated using the area under the curve (AUC) of a receiver operating characteristic (ROC) (Fielding and Bell 1997). The AUC is a threshold independent metric and commonly varies between 0.5 (for chance performance) and 1 (perfect fit).

Table 2 Results of the Generalized Linear Mixed Models (GLMMs) Performed to Explore the Influence of Ecosystem Functioning on the Golden Eagle’s Reproductive Success
Full size table

Finally, we added the “time-lag response” (‘TLagFi’) as fixed factor to explore the effect of fire disturbance on the species’ reproductive success. We fit the ‘integrative model’ plus the ‘TLagFi’—iteratively from 1 to 5 years after fire—to identify possible time-lagged responses of the species to surface energy budget disruptions.

Results

We analysed nesting success from 168 Golden Eagle nesting attempts representing 59 territories over a 17-year period. Models fitted with remotely sensed descriptors of ecosystem functions related to carbon cycle (that is, vegetation productivity) yielded the highest performance (AUC values of 0.78; Table 2). Models calibrated with ecosystem functional variables related to surface energy balance (that is, seasonal heat dynamics and radiance balance) also showed good predictive accuracy (AUC of 0.76 and 0.71, respectively; Table 2). The ‘integrative model’ (that is, the model fitted with the most important descriptors of the ecosystem functions) showed good performance (AUC value of 0.76; Table 2). Four remotely sensed ecosystem functional descriptors were correlated with the breeding success of the Golden Eagle in our study region (Table 2 and Figure 2). The standard deviation and the annual mean of LST—descriptors of the seasonal heat dynamics—were found to be important variables for the ecological fitness of the Golden Eagle (βSD-LST = 0.30, ΣWi SD-LST = 1; βMEAN-LST = -0.15, ΣWi MEAN-LST = 0.76; Table 2). The annual mean of EVI time series (surrogate of annual vegetation productivity) and the date of minimum Albedo (descriptor of radiance balance) also correlated with the reproductive success of the species (βMEAN-EVI = -0.18; ΣWi MEAN-EVI = 0.87; βDMi-ALB = -0.13; ΣWi DMi-ALB = 0.66; respectively; Table 2).

Figure 2
figure 2

Inter-annual variability between 2001 and 2017 for the most important explanatory variables of the reproductive performance of the Golden Eagle: annual mean EVI (EVI_Mean), standard deviation and annual mean of LST (LST_SD and LST_Mean) and the date of the minimum Albedo in the growing season (ALB_DMi). For all plots, colored lines indicate mean values across the mountain systems while the transparent colored areas indicate the error limits defined by the median range values. All variables were standardized before modelling (z-score values).

Full size image

The performance of the ‘integrative model’ increased up to AUC values of 0.80 after including the time-lagged response to surface energy budget disruptions. The burnt area—surrogate of the surface energy budget disruptions—was found to affect the ecological fitness of the species, but only 3 years after the fire events took place (ΣWi BA-T3 = 0.62).

Discussion

Our study provides evidence for the influence of the ecosystem functioning on the fitness of an apex predator (that is, a predator at the upper trophic level). From the initial dataset of 12 remotely sensed ecosystem functional descriptors, 4 of them (the standard deviation and annual mean of LST, the annual mean of EVI and the date of the minimum Albedo) were found to affect the reproductive performance of this residual population of Golden Eagle in NW Iberia over the last 17 years (see Table 2 and Figures 1 and 2). Fire disturbance also affected ecological fitness of this apex predator, but with a limited positive effect at 3 years after fire (that is, a time-lagged response to surface energy budget disruptions). Despite the moderate predictive capacity of our models (AUC values from 0.71 to 0.8), we must be aware of the well-known relevance of other abiotic and biotic factors—which are not included here—that determine the reproductive success of this apex predator (such as human pressure or outbreaks of rabbit disease).

The standard deviation of the intra-annual LST dynamics—a descriptor of the seasonal heat dynamics—positively affected the reproductive performance of the species (see Figure 2 and Table 2). According to our results, the higher the interannual variation in LST, the higher the probability of reproductive success (see Table 2). This pattern might be explained by the marked bioclimatic gradient of our study area, in the transition between the Eurosiberian and Mediterranean regions (Rodríguez Guitián and Ramil-Rego 2007). The breeding territories have been traditionally more productive in the eastern part of our region (Gil and Tapia 2009; Lado and Tapia 2012), characterized by stronger intra-annual climatic variation (that is, dry summers and mild, wet winters) than the Eurosiberian Atlantic zone (characterized by a humid climate moderated by the influence of the ocean, with cold winters and the lack of a distinct dry season). This correlation between the seasonal variation in LST and the reproductive performance can be also explained by altitudinal and topographic gradients, since the Golden Eagle tends to occupy and breed in mountainous areas on cliff ledges with Mediterranean climatic influence (Tapia and others 2007; Lado and Tapia 2012), characterized by strong thermal gradients. In addition, a large interannual variation in LST was observed over the last 17 years in our region (see Figures 2 and 3), due likely to the joint effect of interannual climate fluctuations and fire disturbance (Goetz and others 2007; Ueyama and others 2014) (Figure 3). Thus, positive feedbacks between climate warming and increased fire disturbance could favour reproductive performance of the species over the upcoming decades in our region, as previously forecasted for other eagles in Spain (Muñoz and others 2013).

Figure 3
figure 3

Inter-annual variability between 2001 and 2017 for the most important explanatory variables of the reproductive performance of the Golden: annual mean EVI (EVI_Mean), standard deviation and annual mean of LST (LST_SD and LST_Mean), the date of the minimum Albedo in the growing season (ALB_DMi) and the burnt area across the 9 mountain systems. For all plots, colored lines indicate mean values across breeding territories while the transparent colored areas indicate the error limits defined by the median range values. All variables were standardized before modelling (z-score values).

Full size image

The annual mean of LST was another relevant descriptor of the seasonal heat dynamics to explain the breeding success of the Golden Eagle. The higher the annual mean of LST, the lower the probability of reproductive success (see Table 2). Despite the preference of this population for areas with Mediterranean climatic influence (that is, with higher summer temperatures than more Atlantic areas), the Golden Eagle brood survival in exposed nests was found to decline as the number of days with maximum temperature at least 32.2 °C increased (Steenhof and others 1997; Kochert and others 2019). Our results confirm that surface temperature can be a limiting factor for the reproductive performance of the Golden Eagle.

The annual mean of EVI time series—a surrogate of vegetation productivity—was found to negatively influence the reproductive performance of the Golden Eagle. Thereby, the lower the vegetation productivity, the higher the probability of breeding success (Table 2 and Figure S2). Although it is widely accepted that higher primary productivity rates enables more energy available for reproduction (Connell and Orians 1964; Brasil and others 2019), areas with low primary productivity in the breeding territories of the Golden Eagle are usually associated with open habitats with sparsely vegetated areas—which might indirectly inform about hunting areas close to the breeding sites. Open shrubland areas are positively selected by the Golden Eagle as hunting habitats because the vegetation structure favours prey detection and hunting success (Tapia and others 2007; Lado and Tapia 2012). Further carrion feeding by the Golden Eagle is frequent when natural prey is scarce (Watson 2010), and this raptor lives in free-range mountainous areas of livestock farming in the study region (Tapia and others 2007). In fact, open areas dominated by low scrublands and pastures also represent the typical habitat selected by the more common prey species in the Iberian Peninsula: European Wild Rabbit (Oryctolagus cuniculus), Iberian Hare (Lepus granatensis) and Red-legged Partridge (Alectoris rufa) (Tapia and Domínguez 2007; Tapia and others 2010, 2014; Watson 2010). Previous research have already showed that the interannual dynamics in primary productivity affect habitat colonization patterns and occupancy dynamics of both primary and secondary consumers (Requena-Mullor and others 2014; Fernández and others 2016). Interannual changes in primary productivity driven by, for instance, climatic fluctuations or natural succession processes as well as abrupt changes caused by wildfires might therefore significatively impact the species’ performance (see Figures 2 and 3).

The date of the minimum albedo in the growing season—a descriptor of the net radiation on land surfaces—was associated with lower probability of breeding success of the Golden Eagle (Table 2 and Figure 2). As mentioned above, cliff ledges are often selected by the eagles for nesting (López-López and others 2007; Tapia and others 2007). Nesting sites are usually located in rugged terrains at high altitudes on bare soils with low vegetation covers. These areas are exposed to high surface temperatures and low primary productivity rates, which translates into high albedo values. Thus, lower probability of breeding success was found to be associated with the date when albedo values are minimum in the growing season (that is, highest values of primary productivity). Land-use change is recognized as one of the major factors driving surface albedo dynamics (see for example, Zhai and others 2015; Hu and others 2019). However, fire disturbance—that involves a critical exchange of energy and matter between land and atmosphere—can also affects surface energy balance by initially reducing surface albedo (Gatebe and others 2014) before recovering (at different rates depending on time since fire and post-fire vegetation recovery (Lyons and others 2008; Gatebe and others 2014; Dintwe and others 2017; Torres and others 2018) to their pre-fire albedo conditions. In fact, our study points out to a 3-year time-lagged response of Golden Eagle to surface energy budget disruptions caused by wildfires over the 17-year time period. Previous studies in Mediterranean Spain showed that European wild rabbit and other potential prey species progressively colonized burnt areas (Moreno and Villafuerte 1995; Tapia and others 2010, 2014), where their abundance increased for at least 5 years after the fire (Rollan and Real 2011). The effect of wildfires on Golden Eagle’s fitness can be explained by their joint impact on the spatiotemporal variation of LST, vegetation productivity and albedo (Figure 3).

In conclusion, the reproductive performance of the Golden Eagle—an apex predator at the upper trophic level—is influenced by ecosystem functioning (Table 2 and Figures 2 and 3). Remotely sensed ecosystem functional descriptors such as the standard deviation of LST (as surrogate of the seasonal heat dynamics), annual mean EVI (surrogate of primary productivity) and the date of minimum albedo in the growing season (surrogate of net radiation on land surfaces) were found to be relevant factors affecting its ecological fitness in a wide bioclimatic transition area in southern Europe (Table 2). Fire disturbance also affected ecological fitness of this apex predator by potentially affecting habitat availability of species at lower trophic levels—with a limited positive effect at 3 years after fire (that is, a time-lagged response to surface energy budget disruptions). Our study confirms the potential usefulness in considering the matter and energy fluxes between land surface and atmosphere to explain the reproductive success of species at the upper trophic levels. Although locally measured factors such as prey abundance and geographic nest location are known to be good predictors of breeding success, annual in-field measurements of prey abundance are very costly and time consuming. In this sense, our study confirms that remotely sensed ecosystem functional variables are good proxies of the reproductive success of the Golden Eagle allowing semi-automatic, standardized, and systematic measurements of its ecological fitness over long time periods. Therefore, it illustrates how remote sensing can improve our knowledge on the key role that ecosystem functioning play on the species fitness and actual habitat carrying capacity, boosting our ability to assess ecosystem functioning effects on apex predators from space.