Choosing what is left: the spatial structure of a wild herbivore population within a livestock-dominated landscape

Study area

The study was carried out in Península Valdés (PV, Fig. 1), a provincial protected area and also a UNESCO World Heritage Site located in the Argentine Patagonia. The area presents a temperate semi-arid climate with a mean annual temperature of 13.6 °C (Barros, 1983), and an average annual rainfall of about 230 mm with a high interannual variation (100–300 mm; Coronato, Pessacg & Alvarez, 2017). The vegetation is represented by the southern Monte and the northern Patagonian Phytogeographic Provinces (León et al., 1998). Shrubs and grass-shrubs steppes dominate northern and central PV with a vegetation cover that varies between 40% and 60%, while grass steppes predominate in the southern part of the area with an average cover of 70% (Fig. 1; Bertiller et al., 2017). The most abundant shrub species are Chuquiraga avellanedae and Chuquiraga erinacea, while the predominant perennial grasses are species with high forage value as Nassella tenuis, Piptochaetium napostaence, Poa ligularis and Poa Lanuginosa (Bertiller et al., 2017).

Extensive sheep ranching for wool production predominates the surface of PV (SHEEP areas; Fig. 1), which is divided by fences into more than 60 properties. Each ranch is subsequently fenced into a series of irregular paddocks of variable size (100–2,500 ha) where the sheep graze on the native vegetation. There is usually one building per ranch permanently occupied by a rural worker. Although there are no permanent residents in areas where sheep ranching has ceased (NOSHEEP areas; Fig. 1), fences and other infrastructure such as water tanks and windmills remain in place.

Field surveys

We conducted ground, line transect surveys (Buckland et al., 1993; Laake et al., 1993) of guanacos at the beginning of the autumn in 2017, totaling 383.4 km surveyed along secondary dirt-roads and tracks (average transect length: 7 km). Surveys were conducted as previously described in Antún et al. (2018).

Selection of predictor variables

To test our hypotheses, we identified natural and anthropic variables as potential predictors of L. guanicoe abundance (Table 1). As a correlate of primary productivity, we calculated the mean Normalized Difference Vegetation Index (NDVI; Travaini et al., 2007) for the spring-summer season of 2016–2017 (from September 21st to March 21st). As described in Antún & Baldi (2019) we calculated the CV of the NDVI to account for variation in vegetation physiognomy, and found that it was larger across shrub steppes than in mixed and grass steppes where grasses with high forage value predominate (see Supplemental Information 1). The variables derived from the NDVI were based on MODIS MOD13Q1 satellite images of 250 m spatial resolution available at https://lpdaac.usgs.gov. Altitude values were obtained from the Digital Elevation Model for South America (resolution of about 220 m) at https://lta.cr.usgs.gov/SRTM1Arc. Subsequently, we calculated the CV of the altitude as an estimate of terrain irregularity, relevant in terms of early detection of predators by guanacos (Bank et al., 2002; Pedrana et al., 2010). The current sheep stocking per paddock was obtained by consulting landowners and rural workers between 2015 and 2017. Data on the location of ranches, permanent water sources, and wire fences delimiting the paddocks were available from the database at our institute but checked and updated in the field between 2013 and 2017 (Antún, 2018). Additionally, we included latitude and longitude as geographic variables. Although the geographical predictors do not have ecological significance, they could account for the effect of other variables not available to be included (Travaini et al., 2007) such as the spatial variation in precipitation (Coronato, Pessacg & Alvarez, 2017; Antún et al., 2018; Table 1). We obtained the values for each variable using the QGIS Open Source Geographic Information System (QGIS Development Team, 2016) and packages reshape2 version 1.4.2 (Wickham, 2007), raster version 2.5.8 (Hijmans et al., 2016) and ggplot2 version 2.2.1 (Wickham & Chang, 2016; R software, version 3.2.1, R Development Core Team, 2015). The range of values of each variable across the study area was included as far as possible in the surveyed tracks. Multicollinearity in predictor variables could make it difficult to separate the effects on the response variable and to compare alternative models (Lennon, 1999), thus we evaluated the collinearity between pairs of covariates taking the values measured at each segment (see below). We considered two predictors not to be collinear when Pearson’s correlation coefficients were <0.7. The variables CV of NDVI and geographic latitude showed collinearity (|r| > 0.7). Thus, we kept the former as we considered its ecological significance and possible effect on the spatial structure of the guanaco.

Estimating the detection function

Using standard distance sampling methodology (Buckland et al., 1993), we fitted detection functions to account for the probability of detecting guanacos for the whole of PV, and for SHEEP and NOSHEEP areas. We evaluated the half-normal, hazard-rate and uniform functions as candidate detection functions (Thomas et al., 2010). The final functions were selected following the procedure detailed in Antún et al. (2018). All analyses were performed using the “Distance” package version 0.9.6 (Miller, 2017) for R.

Density surface model (DSM)

Each transect line was divided into smaller segments of 1.8 km in length and 2 km wide (Miller et al., 2013; Antún & Baldi, 2019), totaling 213 segments in all the study area of PV, 164 located in SHEEP and 49 in NOSHEEP. Subsequently, each observation was assigned to its corresponding segment according to its location. The size of the segment was defined according to the information available for the population of guanacos of PV (a sedentary population with an average home range of 5.5 km²; Burgi, 2005), the detection functions and the length of the transects (Schroeder et al., 2014). When the detection function did not include other covariates than distance, the probability of detection was constant for all segments, and therefore we modeled the species abundance per segment using generalized additive models (GAMs; Wood, 2006) with the “count method” (Hedley & Buckland, 2004). Contrary, when the detection function included the group size as a covariate, the abundance was estimated by GAMs but with the Horvitz–Thompson-like estimator (Hedley & Buckland, 2004). The abundance in each segment was estimated including the probability of detection and described as the sum of smooth functions of uncorrelated predictor variables measured at the segment. Restricted Maximum Likelihood (REML) was used for smoothness selection of all the GAMs performed (Reiss & Ogden, 2009; Wood, 2011). The concurvity (the non-linear multicollinearity) could lead GAMs to produce unstable or even wrong estimates of the covariates’ functional effects (Peng, Dominici & Louis, 2006; Gu, Kenney & Zhu, 2010). Consequently the concurvity degree of the smoothing terms (Wood, 2006) was evaluated after fitting the models (Miller et al., 2018) and so we guarantee that any smoothing term was approximated by one or more of the other smoothing terms in the model. Following Miller et al. (2013) for each GAM we explored three response distributions including: Tweedie, negative binomial and quasi-Poisson. For each distribution we built a “base model” and considered all the covariates as univariate smooths. We performed the covariate selection in each base model by removing the non-significant covariates (with approximate P-values > 0.001; Marra & Wood, 2011). We obtained three models as final candidates for each area. Finally we selected the best fitting DSM for the whole PV, SHEEP and NOSHEEP areas based on the inspection of residual plots (Miller et al., 2013). The concurvity measures were very small in all the models evaluated, suggesting negligible concurvity (Wood, 2006). Residual autocorrelation was checked as previously described in Antún & Baldi (2019). Models were fitted using the ‘dsm’ package version 2.2.16 for R (Miller et al., 2018; modeling procedure can be checked at http://github.com/DistanceDevelopment/dsm).

Abundance estimation

Our study area was split into 4 km² cells. Subsequently, we excluded zones with values beyond the surveyed range of the significant covariates and areas adjacent to the coastal limits of the study site or inside the salt pans as they represent marginal habitat of the study area that have not been surveyed. Then, we obtained a prediction surface of 3,196 km² for all the PV area, 2,616 km² for SHEEP and 580 km² for NOSHEEP areas. Finally, we predicted the number of animals for each cell according to the selected DSMs and subsequently obtained an overall estimate of abundance in each area.

Competitive interactions and conflict with humans

As intermediate feeders (sensu Jarman, 1974) able to forage on grasses and dicotyledonous plants, guanacos and sheep largely overlap their diets increasing the potential for interspecific competition (Baldi et al., 2004). Spatial and temporal segregation between species have been previously recorded at different sites and guanacos were less abundant where the proportion of preferred plant species was higher, supporting the hypothesis of competition for forage (Baldi, Albon & Elston, 2001; Burgi et al., 2011). Changes in the abundance and displacement of wild ungulates in presence of livestock have been documented across a range of socio ecological systems such as African savannas (Fritz, Garine-Wichatitsky & Letessier, 1996; Hibert et al., 2010), North American grasslands (Loft, Menke & Kie, 1991) and Trans Himalayan rangelands (Mishra et al., 2004) among others.

Although the observed patterns in this study—such as contrasting densities of guanacos in sympatric and allopatric conditions in relation to sheep—are consistent with the hypothesis of interspecific competition, the influence of human presence is significant. Other studies reported that the probability of finding guanacos decreased as human presence increased (Pedrana et al., 2010; Schroeder et al., 2014), reflecting the threats imposed by chasing and hunting and therefore the additional limitations to select habitats. Although Península Valdés is a protected area, guanacos are still harassed by landowners and rural workers due to the conflict with sheep ranching for forage resources and water (Antún, 2018). To what extent is the observed pattern of variation in guanaco abundance the result of competitive interactions with sheep or human harassment? In this study, the predictive variables “sheep stocking rate” and “distance to the nearest ranch building” allowed us to assess their effects independently (see collinearity and concurvity in subsections Selection of predictor variables and Density surface model). Additionally, our results showed that the effect of sheep stocking rate is stronger than that of the distance to inhabited buildings (see Fig. 3), but we cannot rule out other human-related variables potentially associated with sheep stocking rate (e.g., hunting pressure). However, as some landowners are currently implementing protocols of coexistence among sheep ranching and wildlife populations, there is an opportunity to conduct studies under controlled conditions of no hunting or harassment to understand the mechanisms underlying the competitive interactions between guanacos and sheep.

What’s left to choose?

Our results suggest that guanacos would still be able to select open habitats—in terms of vegetation physiognomy—within a modified landscape where human activities predominate. Across sheep-free areas where the average density of guanacos almost triples compared to the rest of PV, we found that they were more abundant at grassland and mixed steppe habitats. Guanacos have been described as “intermediate feeders” with perennial grasses and evergreen shrubs accounting for up to 60% of their diet (Bonino & Pelliza-Sbriller, 1991; Puig, 1995; Puig et al., 2001; Somlo, 1997; Baldi et al., 2004). Also, open habitats favored the formation of larger groups of guanacos compared to tall shrub dominated habitats. Marino & Baldi (2008) found that group size was negatively related to individual time spent vigilant, while individual time invested in foraging increased. Therefore, the availability of mono and dicotyledonous plants and the benefits derived from grouping behavior can be influencing selection of mixed steppe habitats by guanacos across NOSHEEP areas.

The average density of guanacos in areas where sheep were removed was almost three times higher than in areas where the activity continued. Although guanacos selected habitat types as shown above, we found a significant influence of the remaining fences on their spatial variation in abundance across NOSHEEP areas. In Península Valdés, fences dividing properties and paddocks add up 2,135 km forming 284 polygons inside our prediction area of 3,196 km² (Nabte et al., 2013; Antún, 2018). The occurrence of guanaco seasonal migration has been associated to the absence of fences and other anthropogenic barriers in the Payunia Reserve in NW Patagonia (Schroeder et al., 2014). In contrast, movements of radio-collared guanacos were found to be limited by fences, which added up 175 km within a 400 km² area comprised by private properties in NW Patagonia (Rey, Novaro & Guichón, 2012; Carmanchahi et al., 2015). Additionally, entanglement was reported to account for up to 6.7% of total annual mortality estimated for the same population, affecting mainly juveniles (Rey, Novaro & Guichón, 2012). Whether or not fences affect the ability of guanacos to move around NOSHEEP areas in Península Valdés, hence imposing further restrictions to habitat selection, is a matter requiring further investigation. Detailed studies on guanaco individual movement can throw light on the processes and mechanisms underlying habitat selection and its relationship with productivity gradients and barriers imposed by human activities.

Methodological aspects

The use of density surface models allows to evaluate the spatial variation and estimate the abundance of populations for either a whole area or any sub-region. The combination of distance sampling methods with spatial modeling techniques provide unbiased estimates of abundance independent of the sampling design (Hedley & Buckland, 2004). However, the distribution of the available roads and tracks might limit the range of values of a given variable surveyed and therefore the extent of the prediction area. In our study area the transects surveyed covered most of the range of each significant variable, thus the extent of the prediction is nearly the whole study area (Fig. 2).

Future prospects

It has been argued that high-density populations of guanacos in Patagonia occur on the few protected areas or on land that is abandoned or where sheep ranching was terminated (Novaro, Funes & Walker, 2000; Baldi et al., 2010). Despite facing threats, the number of guanacos in Península Valdés increased markedly during the last 25 years. The average density reported in this study (11.71 ± 0.99 guanacos.km⁻²) is markedly higher than previous estimates for Península Valdés, that turned out to be as low as 1 guanaco.km² (Baldi, Campagna & Saba, 1997), and similar to the density reported for another large protected area, the Payunia Reserve in NW Patagonia (12.28 ± 3.69 guanacos.km⁻²; Schroeder et al., 2014). Although 80% of the land in Península Valdés is still devoted to sheep ranching, the overall stock is probably the lowest in decades (Evolución Existencia de Ganado Ovino 2005–2014, 2016). Both the consolidation of Península Valdés as a protected area, declared as World Heritage Site in 1999, and the decline of sheep ranching could have contributed to the recovery of the population of guanacos. However human—carnivore conflict persists and pumas (Puma concolor) are still chased and killed by rural workers, the occasional presence of the native predator of guanacos has been reported back in PV during the last few years (D’Agostino, 2018). If management actions oriented to the coexistence of wildlife and human activities are implemented, we believe there is an opportunity to restore functional populations of native herbivore and carnivore species in Península Valdés.