Testing remotely-sensed predictors of meso-carnivore habitat use in Mediterranean ecosystems

The legacy of human use of Mediterranean ecosystems results in spatial and temporal heterogeneity of resources for wildlife. Understanding wildlife use of these ecosystems may be improved by including information on ecosystem type, structure, and function extracted from remote sensing data. To assess whether we can improve our understanding of wildlife-habitat use by including information on ecosystem type, structure and function. We tested whether remote sensing derived descriptors of ecosystem type, structure (tree cover and patch size) and function (productivity and stress) determine the habitat of stone martens (Martes foina), common genets (Genetta genetta), and European badgers (Meles meles) in southern Portugal. We linked radio-tracking data from five stone martens, five genets and eight badgers with aerial photography, and some spectra-selectivity to classify vegetation, its structure, productivity and drought stress. Statistically-derived generalized linear mixed regression models using combinations of remotely sensed descriptors of ecosystem type, structure and function, performed better than single ecosystem type descriptors. Inclusion of information on ecosystem functioning in predictive models of habitat use is more informative than ecosystem type alone, suggesting functional relationships between wildlife and their habitat. However, inclusion of both ecosystem type and function maybe limited to finer spatial resolutions. Our results illustrate the untapped potential of remote sensing to provide detailed descriptors of habitat at adequate spatial scales, now that they are freely available and are systematically collected over space and time. This information adds useful insights on wildlife-habitat relationships under changing patterns of land use and climate.


Introduction
Understanding the degree to which species can cope with the effect of land use and climate changes on its habitat is fundamental in ecology, conservation, and sustainability (Morrison et al. 1998;Scott et al. 2002;Boyce 2006).The last decades have seen a proliferation of modeling studies on species niche (Pulliam 2002;Guisan and Thuiller 2005), resource selection (Boyce 2006), habitat suitability (Hirzel et al. 2006;Zielinski et al. 2006), and how these relationships dictate species geographic distributions (Guisan and Zimmermann 2000;Rushton et al. 2004;Moisen et al. 2006).These inferences have shown that habitat can be described in many forms, from descriptors of ecosystem type (Scott et al. 2002), spatial and temporal structure (Zielinski et al. 2006;Santos 2010), and ecosystem functioning (Wulder et al. 2004), to inclusion of multiple (cumulative and dynamic) effects between ecosystems and climate (Tingley et al. 2012;Santos et al. 2014).This complicates both theoretical and empirical assessments as the number of variables to consider increases, as do the number of replicates needed for a comprehensive understanding of the underlying mechanisms that link species to the ecosystems they inhabit.Understanding habitat becomes therefore highly limited by the ability to measure enough parameters and the precision at which measurements of animal behavior and their environmental context can be obtained.
Currently, much of our knowledge on specieshabitat relationships comes from natural history assessments and observations, museum records, and, more recently, from non-invasive molecular approaches (Long et al. 2008), and animal locations obtained through remote assessments using technologies like VHF radio-telemetry, global positioning systems (GPS) and satellite tracking (Perras and Nebel 2012).For example, to understand large scale patterns of diversity and species distributions, museum and other type of observation records are highly valuable if they cover large spatial areas at sufficient detail (Tingley and Beissinger 2009).At this scale telemetry data becomes impractical because of inherent logistic and budget constraints of sampling large areas, except perhaps for migrating species covering large areas (for example see Aarts et al. 2008;Bischof et al. 2012).Inversely, the sparse distribution of museum data that exists at the scale of management units makes these data barely usable, except as a reference data set.At the spatial scale of management, higher precision data provided by telemetry is essential to understand the mechanisms through which species use the habitat.In fact, telemetry data can resolve the minutiae of movement, foraging and resting activities, and interactions with other species.However, these data require matching with adequate environmental context information.
Most commonly used measurements of environmental context include remote sensing descriptors of the Earth surface.These include climate variables, satellite imagery, derived land cover classifications and ecosystem function variables, and digital elevation models and their derivatives (aspect, elevation, and slope).Several authors suggested that these large scale metrics are sufficient for the goals of predicting species distributions at the national and global scales (Guisan and Zimmermann 2000;Guisan and Thuiller 2005;Meyer and Thuiller 2006).However, ultimately the management of species for which these global data sets are produced, occurs at local and regional scales.At these scales, the spatial resolution of such data sets ([1 km cell size) is often too broad to include detailed information.Since species responses to habitat are scaled to their body size (Gehring and Swihart 2003), there may be a mismatch in the spatial resolution at which a species relates to its habitat and that of some remote sensing products, especially for smaller species.A potential solution is to find imagery at appropriate spatial resolution (for example, Landsat at 30 m, or airborne aerial photography at\1-5 m).For example, Landsat data has been actively used to produce high resolution land cover maps, which have been applied to resolve wide ranging animals' habitat use (Schadt et al. 2002;Seoane et al. 2004;Gottschalk et al. 2005;Requena-Mullor et al. 2014).Landsat can provide semi-automated land cover classifications, and a wide variety of metrics to describe ecosystem functioning (Wulder et al. 2004), which have barely been used to understand wildlife-habitat relationships (Kerr and Ostrovsky 2003).Further, Landsat data is now available at no cost, with multi-temporal resolution (see Wulder and Masek 2012 and papers in the special issue).However, for some species finer resolution aerial photographs may better match the scale at which species relate to their habitat.
Land cover classification depicts only one of the components of a species habitat-its type.While ecosystem type can be a surrogate to the type of resources available to a species, it may not be sufficient to describe all the facets that make habitat for a given species.For example, species respond to heterogeneity (Cockburn and Lidicker 1983;Pickett and Cadenasso 1995), either within (canopy cover, structure, etc.) or across (edges, gaps) the ecosystem types they inhabit, which can be estimated using remote sensing products (Seixas 2000;Turner et al. 2003;Vega-Garcia and Chuvieco 2006;Goetz et al. 2007).Species also respond to habitat quality.For example, Saba et al. (2008) showed that the foraging and nesting frequency of the worldwide population of leatherback turtles (Dermochelys coriacea) is a function of resource quantity and persistence as estimated using satellite data.Habitat quality and persistence can be inferred from other satellite imagery products, including vegetation productivity (Oppelt and Mauser 2004), phenology (Di et al. 1994), and stress (Seghieri et al. 1995).Integration of these types of information (type, canopy cover, productivity and stress) can further aid understanding about why animals select given parts of the landscape for their movements and establishment of home ranges (Nielsen et al. 2005;Neumann et al. 2015).
Our goal was to understand whether ecosystem type, structure and function explained mesocarnivore habitat use.To do so we matched telemetry data to remote sensing derived descriptors of ecosystem type (land cover type), structure (canopy cover and patch size) and function (productivity and stress).More specifically we asked: (1) Which descriptors of ecosystem type, structure and function are the best to assess mesocarnivore habitat?And (2) Is there an effect of ecosystem type, structure and function on mesocarnivore use of their habitat?To answer these questions we selected an inherently heterogeneous landscape, the cork oak woodlands in the Mediterranean climates of southern Portugal.We selected as focal species three co-occurring mesocarnivores, the stone marten (Martes foina), the common genet (Genetta genetta), and the European badger (Meles meles), because these species have different patterns of habitat use (e.g.Santos and Santos-Reis 2009;Soto and Palomares 2015), at a scale that matches that of aerial photography and Landsat satellite data.The marten and the genet are arboreal and solitary, and the badger is cursorial, ground dwelling and social; all of them are nocturnal and omnivorous (Gittleman 1989).We predicted that incorporating different descriptors extracted from the Landsat imagery will be more informative than ecosystem type alone (from land cover classifications), since these species use multiple ecosystems for different resources/activities (Rosalino et al. 2004(Rosalino et al. , 2005b, c;, c;Santos-Reis et al. 2004;Santos and Beier 2008;Santos and Santos-Reis 2009), require cover for their movements and resting (Rosalino et al. 2004(Rosalino et al. , 2005c;;Loureiro et al. 2007), and use fruits as important food resources (Rosalino et al. 2005a;Santos et al. 2007;Rosalino andSantos-Reis 2002, 2008), all of which are features detectable by remote sensing imagery.

Study area
The study was conducted in a 20 km 2 area in Serra de Gra ˆndola (Alentejo, Portugal, Fig. 1), which is part of the LTSER Montado platform (http://www.ltsermontado.pt/).The area is dominated by cork oak woodlands (Quercus suber), with patches of holm oak woodland (Q.ilex), pastures, Tasmanian blue gum plantations (Eucalyptus globulus), riparian vegetation (dominated by alder Alnus glutinosa, elm Ulmus spp., and blackberry Rubus ulmifolius), orchards (mainly pear Pyrus bourgeana, fig Ficus carica, and loquat Eryobrotia japonica), olive yards (Olea europaea) and small urban areas and scattered farms.Topography is moderate, with gentle slopes and low altitude (159-238 m).Climate is Mediterranean with Atlantic influence, with mean annual precipitation levels of 500 mm.One temporary stream-Castelhanos-runs along the eastern border of the study area.Human activity in the study area is concentrated in one small village-Santa Margarida da Serra-and several isolated farmhouses, but their effects are extended to areas where cork extraction and livestock production occur, and to hunting areas and timber-producing stands.Legacy of human presence in these ecosystems has created a heterogeneous landscape in which humans and wildlife coexist (Santos and Thorne 2010).

Animal trapping and radio-tracking
From 1997 to 2001 we conducted two animal capture campaigns: 1997-1999 for genets and stone martens and 2000-2001 for European badgers.For stone martens and genets, we set a grid of live traps (Tomahawk Live Trap Co., Wisconsin, USA) baited with canned sardines.Traps were visited daily and rebaited when necessary.For badgers we used a similar approach complemented with soft-catch leg-hold traps (Victor #2) in the vicinity of their setts.Traps were visited multiple times a day and set at sufficient distance to assure that the captured animals could not enter a sett.Captured animals were tranquilized, measured (total length, weight), sexed and aged (tooth wear), and fitted with radio collars (Telonics for genets and martens (Telonics Inc., Arizona, USA), and Biotrack for badgers (BioTrack, Dorset, UK)), and later released into their capture location (Rosalino et al. 2005b).Captured animals were handled following all the recommendations of the Animal Welfare Protocol of the European Union and with the capture permission of the Portuguese Instituto da Conservac ¸a ˜o da Natureza e Biodiversidade (ICNB).
We radio-tracked five stone martens (2 females, 3 males), five genets (2 females, 3 males), and eight European badgers (4 females, 4 males) during the study period.Animals were monitored using two types of surveys, (1) focal samples-when the animal was located continuously for 24 h or while active; and (2) daily locations, either at day or night, to obtain locations of the animal in their resting sites or when active.The first type of survey was conducted using triangulation, whereas the second was through homing-in procedures (for details see Santos-Reis et al. 2004;Rosalino et al. 2005b;Santos and Santos-Reis 2009).Radio-tracking data was entered into Tracker (Camponotus AB and Radio Location Systems AB 1994) and bearings and distances were converted to easting and northing locations.Location data was then exported to ArcGIS 10.1.(Redlands, California, USA) to extract habitat variables (see Remote sensing data).Telemetry positional errors were assumed to be 1 m when animals were inactive, and errors during activity averaged 115 ± 30 m (range 10-269 m) (Santos-Reis et al. 2004;Rosalino et al. 2004Rosalino et al. , 2005b)).We took the pixel corresponding to the point at the center of the error polygon to match active locations positional errors with the minimum mapping unit of the habitat predictors (2 m for aerial photography and 30 m Landsat derived products).The alternative approach would be to create an error polygon and calculate average values for each polygon, but we opted to preserve the resolution in the remote sensing data as the analysis unit.
We used the radio-tracking data to create home ranges for each individual using the fixed kernel estimator method (Worton 1989).This method requires spatial and temporal independence of radiotracking locations and we therefore used Moran's I index (Moran 1950) to evaluate whether spatial autocorrelation was significant on species presence data.The index value is used to calculate the Moran's I statistic, which tests the null hypothesis that there is no spatial autocorrelation through comparing the I statistic to a normal distribution (Cheng and Stephens 1989).This autocorrelation coefficient measures the similarity in the spatial patterns of the variables (Fortin et al. 1989) and varies from -1 (perfect negative spatial autocorrelation) to 1 (perfect positive spatial autocorrelation), with values close to 0 representing no spatial autocorrelation.We calculated temporal independence by determining the time necessary for a species to cross its home range-time to independence-and then use this value as the time lag in between consecutive locations (Swihart andSlade 1985, 1986).We used the 95 % fixed kernel to estimate species home-ranges (to avoid potential outliers in species detections) and used the home range boundary to derive a set of random locations in the same number as the independent locations used to generate the home range.This would reflect our ''pseudo-absence'' data for the modeling approach (see data analysis section).We call these locations as pseudo-absences, because they do not correspond to true absences (surveyed locations where the animal was not detected).

Remote sensing data
We acquired data from two sources of remotely sensed information: low-elevation aerial photography and Landsat satellite imagery.Natural color aerial photography was acquired in 1999 at 1 m resolution.Two Landsat TM scenes (30 m ground resolution) over the study area were acquired for June of 1998 and 2000.We chose only June because we wanted to match (as much as possible) the remote sensing data to the time of acquisition of the radio-tracking data (1997)(1998)(1999)(2000)(2001), given the reduced availability of low cloud cover Landsat data in the archives for our study area in 1998.We used the 1998 Landsat data for genets and stone martens and the 2000 Landsat data for badgers.The Landsat scenes were preprocessed to convert from radiance to apparent reflectance using standard remote sensing tools available in ENVI v.4 (ITT, Boulder, Colorado USA).The different image dates were cocalibrated by selecting pseudo-invariant targets (very bright and very dark pixels in the image), determining the regression line between the two image dates and applying this invariant target regression to the image.

Ecosystem type: land cover type
To describe ecosystem type we used three different land cover classifications.First, aerial photography was photo-interpreted to ten main land cover classes: (1) dense ([50 % cover) cork oak woodland with and (2) without understory, (3) sparse (\50 % cover) cork oak woodland with and (4) without understory, (5) riparian vegetation, (6) pastureland, (7) orchards, (8) eucalyptus plantations, (9) reservoirs, and (10) urban areas and scattered farms.The photo-interpreted aerial photography allowed us to delineate a higher number of land cover classes, as for example the presence of riparian vegetation, grasslands and small orchards, which are often not identifiable by satellite remote sensing due to the small patch sizes of these cover types.Second, we downloaded the CORINE land cover data set (EEA 2002), which is a supervised classification of Landsat TM satellite imagery from imagery acquired in 1999 and 2000.This classification produced twenty three land cover types in Portugal, four of which were present in our study area: broad leaf forest (cork oak woodland), agroforestry, grasslands, and transition woodlandshrubland (http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2000-clc2000-seamless-vectordatabase).Third, the two Landsat scenes were used to calculate the subpixel composition of each of the 30 m pixels using Linear Spectral Unmixing (LSU; Ustin et al. 1986).This algorithm has been developed and applied successfully to numerous Landsat data sets to estimate the proportions in a mixture (Ustin et al. 1993).A set of endmembers (pure pixels) from each class is required to estimate the class proportions within the pixel, and the number of endmembers must be limited to a few spectrally distinct classes to avoid redundancy.We created four endmembers, forest, grassland, bare soil and urban, which represent the variability of structural land cover types in the study area.The inclusion of an additional endmember for shrubland did not improve the separation and therefore we decided not to include it in the final analysis.Further, the grassland endmember was most often mixed with bare soil (dry grassland and soil have very similar spectral signatures, at the Landsat spectral resolution) and thus shows negative values on the final LSU results.For simplification of the results we subtracted the grassland LSU value from the soil and show results for an aggregated grassland/soil class.

Ecosystem structure: canopy cover and patch size
We used tree canopy cover as a proxy for ecosystem structure, since these species require cover for their movements and resting (Rosalino et al. 2004(Rosalino et al. , 2005c;;Santos-Reis et al. 2004;Loureiro et al. 2007).We used the approach developed by Carreiras et al. (2006) to derive Tree Canopy Cover (TCC).These authors related TCC to raw reflectance, tasseled cap transform bands and vegetation indices (Eq.1).The best regression model for southern Portuguese oak woodlands used raw reflectance data from Landsat: where b i is the Landsat TM band and i is the band number (3 through 7).The second metric of ecosystem structure was patch size.We calculated patch area for patches delineated by each classification scheme as described in the section above.We calculated the area of the polygons from the aerial photo interpretation, and for the raster classification we counted the number of pixels per class type and multiplied it per pixel size.

Ecosystem function: productivity and stress
Plant productivity describes resource quantity and quality, which is important as these species often use fruits as food resources (Rosalino andSantos-Reis 2002, 2008;Rosalino et al. 2005a;Santos et al. 2007;Loureiro et al. 2009), and it can also be a surrogate for other food resources such as rodents and insects'abundance (Owen 1988).Previous work has shown a relation between vegetation productivity with the canopy reflectance measured by Landsat (Tucker 1979;Huete et al. 1997), especially in the red edge region (650-850 nm, Landsat bands 3 and 4; Vogelmann et al. 1993;Curran et al. 1995).This region of the electromagnetic spectrum is related to the fraction of intercepted photosynthetically active radiation (Gamon et al. 1995;Ludeke et al. 1996), and is related to the distribution of plant communities, vegetation biomass, land degradation and vegetation quality for herbivores and omnivores (Pettorelli et al. 2005;Wiegand et al. 2008).To measure productivity we calculated a suite of vegetation indices well established in the remote sensing literature (Table 1).We calculated the Normalized Difference Vegetation Index (NDVI), the Green NDVI (a modified version of the NDVI to account for the green reflectance peaks), the Simple Ratio Index (SRI), and the Enhanced Vegetation Index (EVI; Huete et al. 2002), which measure the vegetation fraction within a pixel (Tucker 1979).These indices measure the difference between the reflectance in the red and the near-infrared bands of Landsat (bands 3 and 4, respectively; and band 2 for Green NDVI).We also calculated the Atmospherically Resistant Vegetation Index (ARVI), which accounts for atmospheric interference in the reflectance values (Kaufman and Tanre 1996), the Soil Adjusted Vegetation Index (SAVI), which accounts for the soil influence in the reflectance signal (Huete et al. 1997), and the Greenness Tasseled Cap transformation (Greenness; Table 1) for the Landsat scene, which is the result of a Principal Component Analysis of the Landsat bands (Crist and Cicone 1984).We also measured plant stress as a proxy of the persistence of habitat quality over time.Plant stress  Enhanced Vegetation Index (EVI) (1 through 7) Crist and Cicone (1984) Landscape Ecol 123 corresponds to a reduction in the concentration of photosynthetic pigments, a decrease in water content and a relative increase in woody plant material, as a result of lack of water or nutrients, adverse climatic conditions, plant diseases, and insect damage (Jackson 1986;Grime 1993).If aggravated, plant stress will result in senescence and death, which could indicate unsuitable feeding areas and potentially less suitable resting sites for the mesocarnivores (Santos-Reis et al. 2004;Loureiro et al. 2007).From a remote sensing perspective, plant stress results in changes in leaf color measured in the visible range of the electromagnetic spectrum (Pen ˜uelas et al. 1995a), changes in water content measured in the water absorption features in the near-infrared (NIR) (Pen ˜uelas et al. 1995b) and shortwave-infrared (SWIR), as well as changes in cellulose and lignin measured in the NIR (Jackson 1986).To measure changes in the pigment contents we calculated two indices of the ratio of carotenoids to chlorophyll, the Structure Insensitive Pigment Index (SIPI; Pen ˜uelas et al. 1995a) and the Plant Senescence Reflectance Index (PSRI; Merzylak et al. 1999).To measure a relative indicator of water stress we calculated two indices of canopy water content (Table 1), the Moisture Stress Index (MSI) (Hunt and Rock 1989;Ceccato et al. 2001) and the Normalized Difference Water Index (NDWI) (Jackson et al. 2004).We also calculated the Wetness (water content) Tasseled Cap transformation (Table 1; Crist and Cicone 1984).
To extract the information from each of the remote sensing derived products corresponding to the animal locations we used STARSPAN (Center for Spatial Technologies and Remote Sensing; http://starspan.casil.ucdavis.edu/doku/doku.php),which allows extracting raster data overlaid with vector data.

Data analysis
The first step was to normalize all the variables to avoid model convergence problems and so that model coefficients are comparable.Normalization was done by dividing the difference between the index value at a given location and the mean index value, by the index standard deviation.Secondly, we tested for correlations between the continuous variables selected among the different sets for ecosystem type, patch size, canopy cover, productivity and stress using the Spearman's correlation coefficient value (all r [ 0.70 were considered highly correlated; Hosmer and Lemeshow 2000).Third, we selected which of the metrics of ecosystem type, structure and function best explain animal locations obtained by radio-tracking.We developed a series of univariate generalized linear mixed models (GLMM; Zuur et al. 2009) with presence and pseudo-absence data (binomial distribution and logit link function) for each of the variables reflecting ecosystem type and function (i.e., ecosystem type: aerial photo, CORINE, LSU; productivity: NDVI, NDVIg, SRI, EVI, ARVI, SAVI, Greenness; and stress: MCI, NDWI, Wetness; Table 1), and using individuals as a random factor.For ecosystem type and function (both productivity and stress), we used a multi-model selection procedure creating a model for each of the predictors, and used the Akaike's Information Criteria (AIC; Eq. 2) to determine the best performing model (and therefore the best predictor of ecosystem type, productivity and stress) by calculating the change in AIC from model to model (DAIC; Eq. 3).
where L(h) is the maximized likelihood value and k is the number of parameters.Models with DAIC values less than 2 have equal empirical support and were considered as best models (Burnham and Anderson 2002).We also calculated the Akaike's weight (w i ; Eq. 4) to determine if we could advocate a single top model or not.
where i is the model run, and D is the difference in AIC for every pair of models within the M set of models.If more than one top model was selected, we used w i to calculate model averaged coefficients.This allowed us to select the variables that best resolved animal presence; calculations were made for each species.Since all productivity and stress variables previously selected with the univariate generalized linear mixed models (see above) were highly correlated (r [ 0.7) and correlated with total canopy cover, we opted to create linear combinations of the variables using Principal Component Analysis (PCA; Zuur et al. 2007).The PCA was run using the best variables to describe ecosystem type and function as selected by the univariate GLMMs.We chose to do these two steps instead of just do the PCA's because while we expected some correlation between remote sensing indices within and among categories (for example, correlations between indices for productivity, and correlation between productivity and stress indices), we did not know how to select them.Using first the univariate approach allowed us to select among indices for each category (indices for productivity and indices for stress), and then using the PCA allowed us to include correlated indices representing different categories.Therefore, the PCA reduces the effect of multi-collinearity among indices of ecosystem productivity and stress, and ensures the orthogonality of the predictors.We used the PCA components that accounted for more than 90 % of data variability as predictor variables in the models.We created two sets of PCA-transformed variables to indicate ecosystem function (which we equate to productivity and stress): (1) ecosystem function including Landsat-specific indices and (2) ecosystem function including only ''universal indices''.For the set of Landsat-specific indices we allowed the inclusion of indices customized for Landsat data as the tasseled cap transformations as well as TCC.The tasseled cap transformations and the total canopy cover are multivariate combinations of Landsat bands, and their coefficients are sensor specific.While tasseled cap transformations have been developed for other sensors (Zhang et al. 2002;Scho ¨nert et al. 2014), and so have canopy cover algorithms (Hansen et al. 2002), here we used parameters that were Landsat-specific.For the set of ''universal'' indices we did not include the indices deemed Landsat specific, that is, it included all the indices in Table 1, except the tasseled cap transformations and the total canopy cover.These ''universal'' indices are those whose formulation is not dependent on coefficients specific for each sensor, although their calculation is still dependent on sensor specifications, such as band width.The ''universal'' indices values should be related across sensors, and therefore test whether this approach could be expanded to other types of data.
We used the main PCA axes, land cover type and patch size as variables to a series of GLMM models corresponding to all possible combinations of those variables.For each case, we used the multi-model selection procedure described above creating all possible combinations of predictors (ecosystem type, structure, and function).The influence of each predictor variable on the dependent variable was assessed by the significance of Wald z-statistic test, i.e., variables which tests resulted in P \ 0.05 were considered to have a significant influence on the dependent variable.The influence of the predictors was also assessed by the 95 % Confidence Interval (CI) around the coefficients for each predictor variable, i.e., the influence of a predictor variable was deemed reliable when the confidence interval around the mean estimate of its coefficient did not cross zero; through the analysis of the CI we could determine if a given predictor variable had a positive (positive CI) or negative (negative CI) effect on the dependent variable.Model's performance was evaluated by calculating the Area Under the Curve (AUC), derived from the Receiver Operating Characteristic (ROC) curve (van Erkel and Pattynama 1998).AUC values *0.7 to 0.9 indicate useful applications of the model results (Manel et al. 2001).Analysis were performed in R 3.1.2(R Core Team 2014) using the packages ''ape'' (Paradis et al. 2004),''lme4'' (Bates et al. 2014), ''AICcmodavg'' (Mazerolle 2015), ''MuMIn'' (Barton 2015) and ''pROC'' (Robin et al. 2011).

Animal locations, ecosystem type and resolution
We obtained a total of 3296 temporally independent locations from the radio-tracked stone martens (n = 441), genets (n = 401), and European badgers (n = 2454) during the period of the study.Spatial autocorrelation was significant for badgers (I = -0.037,P value = 0.001) and stone martens (I = -0.019,P value ( 0.001), but not for genets (I = 0.004, P value = 0.06).
The three land cover data sets showed consistent higher numbers of locations on forest/woodland for all the three species (Table 2).The second most used ecosystem type was dependent on the land cover classification, with the aerial photography identifying the high use of pastureland by all species, which is probably related to the CORINE grassland class.The spectral unmixing did not add much to the aerial photo and Landsat land cover classifications, because species presence was divided between forests (72-76 %) and urban areas (14-22 %; Table 2).

Animal location data and predictor variable selection
The continuous indices of productivity and stress were highly correlated.From the original set of variables and along with the land cover classifications, the univariate models selected NDVI, NDVI g , SAVI, EVI, Greenness, SIPI, and PSRI for stone marten, all variables for genet, and NDVI, SAVI, Greenness, and MSI for badger (Online Appendix 1).In Online Appendix 1 we show a list of all the variables and their DAIC for the univariate GLMM for each variable set (i.e.land cover, productivity and stress) and for each carnivore species.Variables were selected when DAIC \ 2. The selected variables were combined using a PCA for Landsat only variables and for ''universal variables'' (see methods for the description of these two sets of candidate variables).Because PCA1 explained most of the variation in the dataset we used it for all the models, with the exception of the genet universal model where we also included PCA2 (Stone marten: Landsat PCA1 = 95 %, Universal PCA1 = 97.1 %; badger: Landsat PCA1 = 95.6 %, Universal PCA1 = 96.7 %; genet: Landsat PCA1 = 91.2%, Universal PCA1 ?PCA2 = 95.4 %; Table 3).In addition, we added patch size and TCC as metrics of ecosystem structure.These were the variables used in the competing models below.

Competing models
The competing models analysis showed that the best model for stone martens and badgers included ecosystem type, function and structure, while the model for the genet included only ecosystem type (Table 4); the effect of each of the parameters was species-specific.The presence of stone marten was significant and positively affected by dense oak woodlands with understory, grasslands, and riparian vegetation, but negatively affected by sparse oak woodland without understory and a combination of productivity and stress, and TCC-PCA1 (Table 4).PCA component loadings show a positive influence of productivity (Greenness) and canopy cover (TCC), and negative of influence of stress (PSRI) on stone marten's presence according to the ''Landsat'' model.The ''Universal'' model PCA loadings indicate a positive association between productivity (SAVI, NDVI, NDVIg and EVI) and stone martens and a negative association with stress (PSRI and SIPI; Table 3).
The presence of badgers was significantly and positively affected by areas of dense oak woodland with understory, sparse oak without understory and orchards, and a combination of productivity (NDVI and SAVI) and stress (MSI; Table 4, Online Appendix 2).NDVI and SAVI productivity indices had a positive influence on PCA1 and stress index MSI had a negative one (Online Appendix 2).On the other hand, badger presence was negatively related to eucalyptus plantations and grassland, patch area and a combination of productivity, stress and canopy cover metrics, which indicated a negative influence of stress index MSI and a positive of productivity index Greenness and TCC on badger presence (Online Appendix 2).
The presence of genets was poorly predicted, and we only found a positive significant effect of dense oak with understorey and eucalyptus plantations.All the models were significant and had a predictive power greater than 60 % (AUC values in Table 4).Genet models have low accuracy (AUC \ 0.7) and therefore should be carefully interpreted.

Discussion
Anthropogenic activities result in habitat loss and fragmentation, and changes in climate, which can produce irreversible modifications in the environmental conditions that allow the persistence of biological communities.Assessments of habitat and its use that can move beyond traditional approaches to measure habitat properties that species respond to, and that are themselves responsible for ecosystem functioning, should be of increasing value.In this study we assessed the use of remote-sensing derived information about ecosystem type, structure, and function to describe habitat use for three species of mesocarnivores.For all the studied species there was a consistent higher use of forest/woodland land cover types.Dense cork oak woodlands with understorey was the ecosystem type that all species use.In addition, stone martens and badgers can also be found in orchards, grasslands and riparian vegetation.Our results suggest that finer spatial resolution (lower pixel size) and higher categorical resolution (more land cover classes) provided better explanation of used ecosystem types.For the European badger and the stone marten, including information on ecosystem structure (patch size and canopy cover), and function (productivity and stress) substantially improved predictions of habitat use based on ecosystem type alone.No conclusive results were found for the genet.
Previous studies in the Iberian Peninsula, showed that cork oak woodlands with understory are important for badgers (Revilla et al. 2001;Virgo ´s 2002;Rosalino et al. 2004Rosalino et al. , 2005cRosalino et al. , 2007;;Santos and Beier 2008) and stone martens (Virgo ´s and Casanovas 1998; Virgo ´s and Garcia 2002; Santos-Reis et al. 2004;Santos and Santos-Reis 2009).The shrub layer in these oak woodlands is removed every 3-5 years to enhance cork production and quality.Apparently, this management activity does not cause avoidance of the altered areas (Santos-Reis et al. 2004).Cork oak woodlands are a source of food for badgers, as they provide acorns which badgers consume in great quantities; in these areas badgers also find beetles (Coleoptera), one of their other main prey items (Rosalino et al. 2005a;Loureiro et al. 2009).The stone marten, as an arboreal species, utilizes the resting sites that cork oak trees provide (Santos-Reis et al. 2004).Cork oak trees live up to 250 years, and as trees senesce cavities are created in the wood (Aronson et al. 2009;Carvalho et al. 2014).
We found a high likelihood of finding stone martens in riparian vegetation and badgers in orchards, and opposing directions of influence of grasslands.There are few studies that show the importance of riparian vegetation for carnivores (Matos et al. 2008;Virgo ´s 2001;Santos et al. 2011).In our study area, riparian areas provide water and food resources (blackberries, small mammals, cray fish etc.; Matos et al. 2008;Pereira and Rodriguez 2010), and are also attractive resting locations in the summer because of their cool microclimates (Santos-Reis et al. 2004;Loureiro et al. 2007).One of badgers main diet component are fruits which are found in orchards (Rosalino et al. 2003;Santos et al. 2007;Rosalino and Santos-Reis 2008).Many orchards are abandoned, representing a past legacy of when the area was highly productive and farms, today abandoned, were inhabited by montado workers and their families (Antrop 1993;Pinto-Correia and Fonseca 2009).Badgers seem to use these quite often and the use of these small orchard areas can probably explain the negative effect of patch size, i.e. a negative relation between the probability of badgers being presence and the size of the patch.Thus, in this landscape smaller patches such as orchards and riparian areas, have higher probability of being used by badgers than larger ones (e.g.pastureland).We found that badgers were negatively affected by grasslands while stone martens' presence was positively related to this land cover type.This is surprising as in other areas of badger distribution grasslands are a main source of earthworms; however, in Mediterranean ecosystems earthworm consumption by European badgers is rare (Rosalino et al. 2003).
No conclusive results were found for habitat use by genets, for land cover type, structure and function.Genet model performance had relatively low AUC values (Manel et al. 2001).Genets more often found in dense cork oak woodland with understory; however, the model runs were statistically inconclusive in spite of the number of locations used for model runs being similar to those used for the stone marten and the radio-tracking occurred simultaneously for both species, eliminating potential yearly variations.Perhaps the genet is the most habitat generalist of the three species, which is potentially a result of its being a non-native species (Dobson 1998).These types of models tend to perform worse for generalist species because it becomes difficult to detect specific habitat features, as they match closely the availability represented in by the random locations (Elith et al. 2006).The genet is a naturalized non-native species and perhaps is not selecting any habitat specifically in our study area.
Our results suggest that adding information about structure (canopy cover) and function (plant The PCAs used in the badger and stone marten models were mostly dictated by plant stress (with the exception of the universal model for the badger which used mostly productivity), which explains the negative and the positive relationships that our models revealed.
The negative relation to plant stress is probably related to avoidance of dry (and hot) areas, especially during the summer.This could be because these areas do not confer protection against the high summer temperatures.Further, and because these results used the Universal model PCA they suggest that the use of these indices can be transferable to other sensors beyond Landsat.Despite the diluted effect of canopy cover in our models, during the summer the outer layer of cork plays an important role in the thermal isolation of resting sites during the day (Santos-Reis et al. 2004), and cover is also an important feature as species prefer to move in sheltered rather than open areas (Rosalino et al. 2004(Rosalino et al. , 2005b, c), c).Our analysis showed that there are some ecosystem function predictors that improved model performance more than others.NDVI (Tucker 1979) is probably one of the most widely used vegetation index in remote sensing and other applications to measure plant greenness (Ludeke et al. 1996;Pettorelli et al. 2011) but is sensitive to differences in background reflectance.To compensate for some of these errors, the soil adjusted vegetation index (SAVI) includes a parameter L (see Table 2 for formulation), which accounts for the effects of soil in the overall pixel reflectance  (Huete et al. 1997).Other productivity indices are also alternative metrics to productivity to respond to the inherent limitations of NDVI.Our results showed that most models used both SAVI and NDVI (Online Appendix 1), suggesting a complementary effect of both indices.This is probably because SAVI gives better estimates of vegetation cover when soil reflectance varies considerably, complementing the information from NDVI.Other authors have also used NDVI as measures of seasonal productivity to which wildlife species may or not respond.For example, Wiegand et al. (2008)  The inclusion of ecosystem function metrics improved predictions of habitat use based on ecosystem type alone but at the expense of a reduction on spatial scale and categorical resolution of the land cover map.There is a wide range of spatial resolutions that have been used in habitat-species studies.For example, some studies looking at the influence of environmental covariates on European badgers' presence have been based on coarser spatial resolutions than those used in our study (e.g.Revilla et al. 2000;50 m resolution;Hammond et al. 2001-100 m resolution;Requena-Mullor et al. 2014).However, others have also used Landsat Thematic Mapper (TM) satellite images or aereal photos/on-site data for assessing habitat use or predict badgers distribution with similar resolutions to ours (e.g.Virgo ´s and Casanovas 1999-on-site; Wright et al. 2000-30 m;Rosalino et al. 2004-2 m;Newton-Cross et al. 2007-25 m).This could lead to a mismatch between the habitat composition within the error polygon and the habitat category associated with the central point of that polygon, which may be a source of uncertainty.However, this uncertainty is inversely related with landscape heterogeneity, i.e. in more homogeneous areas the uncertainty may be low as there is a higher probability that the central point of the error polygon will represent the dominant habitat in the polygon.In highly heterogeneous areas this bias may affect results and should be considered in the analysis.Our study area is relatively homogeneous so we expect this bias not to considerably affect the results.
Our results corroborate our predictions that adding additional information on ecosystem structure, and function would be more informative than ecosystem type alone (even at finer spatial and categorical resolution), providing a unique insight into animal habitat use.The advantages of using such an approach are the inclusion of the habitat predictors that represent different ecosystem functionality and go beyond the traditional wildlife-habitat studies that reflect species use of land cover types.This insight can bring about new ways in which to describe the ecological niche of species approximating the underlying mechanisms and processes, which can be measured using indices that are applicable across a range of sensors.However, such applicability may be restricted as our second main result is that the spatial resolution of the land cover classification is very important for predicting habitat use.Therefore, these results suggest a trade-off between the inclusion of metrics of ecosystem functioning and the scale of the analysis-a spatial resolution that probably matches that of the perception abilities of mesocarnivores.The Landsat archive and analytical software are currently freely available, allowing inferences to current and past dates, and predictions into the future potential responses of species to effects of climate changes in their habitat (Santos et al. 2014).Climate change predictions for southern European countries indicate decrease in rainfall and in number of rainfall events, and increased summer temperatures (Santos et al. 2002;Reid 2006).This will likely result in changes in plant phenology and productivity, and plants will show higher stress levels, decrease in canopy cover, which are expected to affect future mesocarnivore habitat use.

Fig. 1
Fig. 1 Study area location in southern Portugal.Plus signs indicate animal locations obtained from radio-tracking vegetation cover.Values range from cover modified to account for green reflectance.Values range from -

Table 1
Vegetation indices related with vegetation physiology and water content

Table 2
Land cover type in locations of stone martens, genets and European badgers

Table 3
Principal Component loadings for each species and model set

Table 4
Averaged model parameters of habitat use of stone martens, genets and European badgers in cork oak woodlands of southern Portugal Hunt and Rock 1989;Ceccato et al. 2001 productivity were linked to brown bear (Ursus arctos) population decreases in northern Spain.Despite all the benefits of NDVI, including its wide use, simplicity, easiness to understand, and reliability, it also has the disadvantage of being influenced by the underlying soil reflectance, and at high values of canopy closure NDVI reaches saturation.For example, Requena-Mullor et al. (2014) used the Enhanced Vegetation Index (EVI) as an alternative to NDVI to assess the effects of ecosystem function on European badger habitat use.Our results were mixed in the selection of the stress metric, while stress was an important descriptor of habitat use.PSRI and SIPI were selected for stone marten and genet models, and MSI for badger models.Plant stress can be manifested in many physiological alterations, including changes in pigment concentration, with a decrease in chlorophyll and an increase in carotenoids, a decrease in water content, and increase in woody material and leaf litter.PSRI, Pigment Sensitive Vegetation Index(Merzylak et al. 1999), and SIPI, Structure Insensitive  Pigment Index (Pen ˜uelas et al. 1995a) as the names indicate, measure changes in pigment concentrations whereas MSI (Moisture Stress Index;Hunt and Rock 1989;Ceccato et al. 2001) measures changes in moisture content.This indicates that the species may be responding to different sources of plant stress, with badger responding to water content, and stone martens and genets responding to leaf color changes and senescence.The specific mechanisms by which each species is associated with indices of plant productivity and stress is unknown, but worthy of future study.