Introduction

Climate projections indicate that drought frequency and severity is likely to increase over much of eastern Australia over the coming century (Hennessy et al. 2007; White et al. 2010). In evolutionary terms, species can respond in one of three ways to changing environmental conditions: extinction, adaptation or stasis. Evidence from the Pleistocene shows that species responses to climatic oscillations varied with topography, latitude and individual species characteristics (Hewitt 2004). However, the most typical response was evolutionary stasis in situ combined with changes in distribution and abundance (Stewart and Lister 2001; Byrne 2008; Magri 2008; Provan and Bennett 2008; Kearns et al. 2010). Locations where conditions are such that species can persist in situ while their populations are generally contracting in range or abundance have been termed cryptic refugia and micro-refuges. Given climate change projections, the ability to identify potential climate change micro-refuges is a research and conservation priority (Keppel et al. 2011; Sublette Mosblech et al. 2011).

Refugia are conceptualised as locations that provide protection from extreme or protracted climatic conditions and are potential source areas for population expansion if conditions outside the refuge become suitable again. There is ongoing debate about terminology (Rull 2009; Stewart et al. 2010; Keppel et al. 2011), with differences primarily based on spatial and temporal scale. Irrespective of the terminology used, refugia are necessarily specific to the type of disturbance (Berryman and Hawkins 2006) and the species of interest (Ashcroft 2010). The actual area required for populations to persist in situ will vary with the intensity and duration of the disturbance, and the life history attributes of the species.

Recent research has focussed on identifying topographically driven micro-climates as potential climate change refuges (Ashcroft et al. 2009, 2012; Ashcroft 2010; Dobrowski 2011). Landscape genetics and phylogeographic analyses are also being used to identify the locations of historical refugia (Hugall et al. 2002; Carnaval et al. 2009; Scoble and Lowe 2010). An alternative approach proposes that locations where mean plant productivity is relatively high and temporally stable compared to other locations of the same vegetation classification could potentially function as drought refuges (Mackey et al. 2012). The theoretical basis for this proposal rests (1) on the general relationship that exists between plant productivity, resources, population size and extinction risk, and (2) on the matching that occurs between species and the habitats that they occupy (Southwood 1988).

Population size is primarily determined by the interaction of the space–time distribution of resources, species’ physiological and life history attributes, and local environmental conditions (Gates 1980; Andrewartha and Birch 1984). The same principle applies to plants and animals because population size is mechanistically connected to resources through metabolism and allometric scaling laws (Enquist et al. 1998, 1999; Carbone and Gittleman 2002; West and Brown 2005). The specific mechanisms that link resources with population size vary depending on how and where the per capita effects on species’ survival and fecundity are the greatest (Huston 1994; Newton 2013). For example, the numbers of several species of migrant European songbirds fluctuate according to rainfall (and hence food supplies) in their African wintering grounds (Newton 2004); and the availability of nesting sites can limit the numbers of hollow nesting birds (Newton 1994). Irrespective of the specific mechanism, as populations become smaller they become more vulnerable to extinction through demographic and environmental stochasticity (Caughley 1994; Gaggiotti and Hanski 2004). Mackey et al. (2012) proposed that variability in the distribution and availability of habitat resources could be represented by space/time variability in vegetation productivity.

Gross primary productivity (GPP), the rate per unit area at which new biomass is produced by the vegetation cover, can be monitored remotely using time series of satellite images (Box 1989). NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS) sensor detects the energy reflected in distinct spectral bands from every part of the Earth’s surface every 1–2 days. Reflectance values are used to calculate the normalised difference vegetation index (NDVI) which has been shown to be sensitive to spatial and temporal variation in the amount of vegetation and its’ condition (Huete et al. 2002). The fraction of photosynthetically active radiation absorbed by a sunlit canopy (fPAR) which is a reliable proxy for GPP (Berry and Roderick 2004) can be derived from the NDVI.

Mackey et al. (2012) demonstrated the potential application of time series of fPAR to the identification of ecosystem greenspots, i.e., locations that maintain relatively high and stable levels of gross primary productivity (GPP) during drought. Integrating the spatial and temporal dimensions of productivity has also been applied to continental scale habitat analysis in relation to dispersive fauna in Australia (Berry et al. 2007), biodiversity monitoring in Canada (Coops et al. 2008), and seasonal dynamics of habitat quality of brown bears in Spain (Wiegand et al. 2008). Mackey et al. (2012) proposed that the ecosystem greenspot index could be used to identify locations that currently function as habitat refuges from fire and drought which in some bioregions are likely to become more persistent under future climatic conditions. If this is true, these locations may have the potential to function as climate change micro-refuges. A conceptual diagram of relationships is shown in Fig. 1. The ecosystem greenspots index, however, awaits validation with in situ data.

Fig. 1
figure 1

Conceptual diagram of the relationships between fPAR, the ecosystem greenspots index and potential micro-refuges. The fraction of photosynthetically active radiation absorbed by a sunlit canopy (fPAR), which is derived from NDVI, is a reliable proxy for gross primary productivity (GPP). We assume that for any given vegetation type (represented by vegetation maps), higher long term mean productivity and more temporally stable productivity (represented by fPAR) will be reflected in the biomass of specific framework species. In turn, the resources provided by these framework species may help sustain populations in situ during drought conditions

Full size image

Validation of the ecosystem greenspots index requires at least two lines of supporting evidence. First, vegetation data are needed to test whether the ecosystem greenspots index is related to vegetation specific gradients in the amount of habitat resources and not simply an artefact of classification error or weed infestation. Second, demographic and dispersal studies of potential beneficiary species are needed to test whether species persistence during drought is supported by the identified greenspot locations. In this study, we aimed to test whether the ecosystem greenspots index can accurately identify habitat specific gradients in the amount of vegetation based habitat resources. We also analysed how landscape and climatic variables were related to the spatial distribution of ecosystem greenspot classes. This study represents the first step towards validating ecosystem greenspots as a tool for identifying potential drought micro-refuges.

Data and methods

The study area

The Northern Midlands of Tasmania, Australia, was selected as the location for validating the ecosystem greenspot index as it is currently the focus of research into regional scale conservation planning tools (Fig. 2). Two critically endangered ecosystems and many vulnerable, endangered and critically endangered species occur in the Northern Midlands bioregion.

Fig. 2
figure 2

Location of study area. The shaded area indicates the Northern Midlands bioregion

Full size image

Calculating the ecosystem greenspot index

Vegetation specific greenspot indices were calculated by combining spatial vegetation data with a time series of fPAR. To calculate the fPAR time series we followed methods that were originally developed for a pre-cursor of the MODIS NDVI by Sellers et al. (1994) then followed by (Roderick et al. 1999; Berry and Roderick 2002, 2004; Mackey et al. 2012). Methods for calculating the ecosystem greenspot index are outlined in detail in Mackey et al. (2012). We used this method rather than the available fPAR MOD15A2 product because it provides higher geographic resolution, i.e., 250 m compared with the resolution of global data which is 500 m. Furthermore, the method of Mackey et al. (2012) does additional processing to remove cloud contamination that is present in the global data and has been corrected to derive a soil adjusted value for Australian conditions. The resulting index identifies potential greenspots within six percentiles for each specific vegetation type within a defined area of interest. The greenspot analysis was restricted to the Northern Midlands bioregion. Index thresholds corresponding to the 10, 25, 50, 75, 90 and 95th percentiles were calculated for each vegetation type.

The source for the NDVI data was NASA’s MODIS sensor, 16 day L3 Global 250 m MOD13Q1 for the period June 2000–July 2011. This period incorporated record low rainfall periods in the Northern Midlands including record low annual rainfall (2008), and record low monthly rainfall totals for February (2003), June (2007), October (2008) and November (2006) (Bureau of Meteorology 2013).

Vegetation data were sourced from TASVEG Version 2.0 (Tasmanian Vegetation Monitoring and Mapping Program 2009), a state-wide coverage at a cartographic scale of 1:25,000. Vector vegetation data were rasterized at 250 m resolution to match the resolution of the fPAR data and each grid cell was classified according to the vegetation type at the cell centre. The final output was a raster coverage in which each polygon had two attributes: (i) the vegetation type, and (ii) a greenspot index percentile calculated for the specific vegetation type within the bioregion.

Study design

To test whether the productivity gradient that underlies the ecosystem greenspot index is related to vegetation specific gradients of habitat resources we quantified vegetation structure within two vegetation types: (i) Eucalyptus viminalis grassy forest/woodland (henceforth woodland); and (ii) lowland Themeda triandra grassland (henceforth grassland). These vegetation types represent two broad vegetation categories that dominated the pre-European and extant native vegetation in the Northern Midlands of Tasmania (Fig. 3). We selected sites to sample across the ecosystem greenspot index within each vegetation type subject to site accessibility (Table 1).

Fig. 3
figure 3

Variation within the two study ecosystems. Photos represent the (a) top, (b) middle, and (c) low greenspot groups for (1) Eucalyptus viminalis grassy forest/woodland, and (2) lowland Themeda triandra grassland respectively

Full size image
Table 1 Sample size of vegetation types and greenspot percentiles
Full size table

We used the biomass of relevant framework species as a surrogate for vegetation specific gradients in the amount habitat resources. Framework species, also referred to as foundation species (Ellison et al. 2005) and structural species (Huston 1994), are the species that dominate the structure and function of an ecosystem. Because of their structural dominance, the identity of framework species is often used as the basis for vegetation classification. Framework species influence the micro-climatic conditions and provide the bulk of vegetation based habitat resources used by interstitial species. The relevant framework species were E. viminalis in the woodland and T. triandra in the grassland.

Within each site, we sampled vegetation within square plots of 200 m2 located to avoid pixel boundaries. For each plot the positions of fifteen random sampling points were generated and uploaded into a GPS for location in the field. At each sampling point, we collected data for three a priori vegetation layers: (i) tussock grasses/sedges; (ii) a middle layer of woody vegetation; and (iii) a canopy layer of woody vegetation. The middle and canopy vegetation layers were defined differently for grasslands and woodlands because of differences in vegetation structure. In grasslands, the middle layer was defined as woody vegetation <2 m in height whereas in forest/woodland sites the middle layer was defined as woody vegetation <5 m in height. Data were collected using the point-quarter method and a radius truncated to 15 m. The point-quarter method is a plotless sampling method that is widely used in vegetation sampling. The method is described and illustrated in Krebs (1997).

In each quarter we collected the following data for each vegetation layer: distance to individual; species identity; height; canopy length; and canopy width. In addition, we measured diameter at the base for tussock grasses/sedges and diameter at breast height (DBH) for the canopy layer of woody vegetation. This method potentially provides 60 distance measures per vegetation layer per site. For the middle and canopy layers, distance to the closest individual was measured. For tussock grasses/sedges, distance to the second-closest individual was measured. Distance to the second-closest individual is preferred as it provides a smaller sampling variance than distance to the closest individual (Pollard 1971; Engeman et al. 1994), however, densities of shrubs and trees were too low to use the second-closest individual method. All vegetation data were collected during February 2013.

Data analysis

Site means were calculated for each vegetation layer and where necessary, data were transformed to meet assumptions of normality. We calculated stem densities using Morisita’s point-based estimator which is more robust to spatial non-randomness than other commonly used procedures (Mitchell 2007; Bouldin 2008). The fPAR variables used to represent the amount of productivity and the temporal stability of productivity were long term mean fPAR and coefficient of variation of monthly mean fPAR. We used linear regression to test whether site vegetation variables were related to site fPAR variables. To test for differences between groups of woodland greenspots we used ANOVA. Woodland greenspots were grouped as follows for ANOVA, top = 10 and 25th percentiles (n = 3), middle = 50 and 75th percentile (n = 4), and low = 90 and 95th percentile (n = 4). All data analysis was conducted in R (R Development Core Team 2005).

We used recursive binary partitioning to analyse the spatial distribution of woodland greenspots in relation to environmental variables. Recursive partitioning embeds tree structured regression models into conditional inference procedures and thereby reduces over-fitting and biased variable selection (Hothorn et al. 2006). The resulting tree and leaf nodes mean that groups can be identified based on interactions between explanatory variables. Recursive partitioning procedures improve on linear regression by allowing for interactions and non-linearities when there are multiple explanatory variables and are useful for spatial mapping (Prasad et al. 2006). The analysis included all 250 m pixels classified as E. viminalis grassy forest/woodland within the Northern Midlands bioregion (n = 5,428). Woodland greenspots were grouped as follows for recursive partitioning: top = 10 and 25th percentiles (n = 149), middle = 50 and 75th percentiles (n = 1,852), low = 90 and 95th percentiles (n = 2,592) and bottom = all remaining sites lower than the 95th percentile (n = 835). Potentially explanatory variables included in the analysis were: annual mean rainfall, topographic slope angle, topographic aspect, topographic position and plant available water capacity (Table 2). The procedure was performed with the rpartv4.1–8 package in R using the default values.

Table 2 Potentially explanatory variables used to analyse the spatial distribution of greenspots
Full size table

Results

Woodlands greenspot index

Field observations confirmed that the mapped vegetation data were accurate in that sites that had been mapped as woodlands were found to be woodlands. Within woodland sites, however, there was variation in vegetation structure and species turnover between sites. In particular, the composition of the shrub layer varied between woodland sites and included exotic species, especially Gorse (Ulex europaeus). The woodland greenspot index correctly identified a gradient in woodland vegetation structure. E. viminalis formed a canopy at sites in the high woodland greenspot group but occurred as isolated emergent trees at sites in the low woodland greenspot group (Fig. 3). Site values for the basal area of the framework species, E. viminalis, and height of the tree layer were significantly positively related to long term mean fPAR but not to the coefficient of variation of monthly mean fPAR (Table 3). Differences between groups of woodland greenspots were significant for the basal area of E. viminalis (F = 9.59, p < 0.007) and for the height of the canopy layer (F = 4.65, p < 0.05) (Fig. 4).

Table 3 Results from linear regressions on fPAR
Full size table
Fig. 4
figure 4

Comparison of woodland vegetation structure. Box and whisker plots (showing the range, the interquartile range and the median) comparing canopy height and basal area of Eucalyptus viminalis between woodland greenspot groups

Full size image

Grasslands greenspot index

There was a high error rate of vegetation classification for grasslands as many sites that were mapped as grasslands were found to have a substantial component of shrubland and woodland vegetation. Sites in the highest percentile of the grassland greenspot index attained their high values because of the presence of shrublands and woodlands in sites that had been mapped as grasslands. After removing sites that included patches of shrubland and woodland we only had data for 11 grassland sites (Table 1). For these sites, values for the basal area of the grassland framework species, T. triandra and basal area of tussock grasses were positively related to the coefficient of variation of monthly mean fPAR but not long term mean fPAR (Table 3). It is noteworthy that even in the subset of “pure” grassland sites there was still considerable spatial variation in the composition and structure of the grassland greenspots with the inclusion of small patches of saline herbland, wetland and lowland grassland complex within sites mapped as lowland T. triandra grasslands. Given the vegetation classification error and the spatial heterogeneity within the remaining grassland sites our grassland vegetation data probably do not accurately represent a productivity gradient of Lowland T. triandra grasslands.

Covariates of the woodland greenspot index

The spatial distribution of woodland greenspot groups in the Northern Midlands was significantly correlated with climatic and landscape variables. Rainfall explained most of the variation followed by the index of plant available water capacity which contains information about soil texture and soil depth among other attributes (Table 2). The scaled variable importance values were rainfall (89), plant available water capacity (7), slope (3) and topographic position (1). There was interaction between explanatory variables so that the distribution of woodland greenspots varied as a function of rainfall, plant available water capacity, slope and topographic position. There were four distinct sets of conditions in which there was >25 % chance that woodland greenspots would be in the highest group. For example, in locations where annual mean rainfall exceeded 786 mm, high woodland greenspots were associated with low, i.e., <0.2 plant available water capacity. In contrast, where annual mean rainfall was less than 578 mm, high woodland greenspots occurred in locations where topographic aspect was southeast and 0.36 < plant available water capacity <0.4. Between 578 and 786 mm annual mean rainfall, high woodland greenspots were more likely to occur in locations where topographic slope >14° and topographic aspect was south–southeast.

Discussion

The woodland greenspot index was a good predictor of variation in the basal area, and hence biomass, of the framework species, E. viminalis. The long term mean fPAR was a better predictor of variation in vegetation structure than the coefficient of variation of monthly mean fPAR. We conclude that the ecosystem greenspot method is informative about the biomass of specific vegetation based habitat resources in woody vegetation. The ability to provide high resolution information about gradients of specific vegetation based habitat resources has immediate applications for conservation planning. Conservation planning is informed, among other things, by the scientific criteria of comprehensiveness, adequacy and representativeness (Commonwealth of Australia 1999). The criterion of representativeness includes consideration of the need to protect variability in the quality of habitats within ecosystems. The information provided by the ecosystem greenspots index has potential application in assessing how representative the protected area estate is with respect to habitat quality. A similar approach has been used to assess spatial bias in protection of and threats to biodiversity in Canada (Andrew et al. 2011). The ability to identify the most productive locations of specific vegetation based habitat resources may also be a useful tool for species level conservation planning.

The ecosystem greenspot method as applied here was not successful for grasslands. To begin, there was a high error rate in grassland classification. Furthermore, there were no gradients in the biomass of the framework species, T. triandra. Most of the identification error was due to the use of the vegetation data at the spatial scale of 250 m. Use of the data at this resolution meant that there were errors in the location of vegetation boundaries and in vegetation classification. Some classification error was also caused by the method we used to classify vegetation when transforming it from vector data into raster data. These combined sources of error resulted in the inclusion of non-grassland vegetation which distorted the ecosystem greenspots index for grasslands. More work is needed to address sources of error and field validation methods for grassy ecosystems in the Midlands of Tasmania.

The finding that the coefficient of variation of mean monthly fPAR was positively correlated with the basal area of T. triandra and the basal area of tussock grasses is counter to our expectation. We predicted that the basal area of framework species would increase as temporal stability increased i.e., as the coefficient of variation decreased. There are a few factors that may contribute to this finding. First, the vegetation data in grasslands is confounded by land management with variation in grazing regimes between sites. Second, specific components of the total vegetation cover may account for different fractions of the fPAR signal. Different leaf functional types exhibit different temporal dynamics which account for different components of the total fPAR signal (Berry and Roderick 2002). Given the amount of variation in grassland vegetation a larger sample size is needed to clarify the relationships. Furthermore, any future field work will need to consider testing different fractions of the fPAR signal.

Conversion of land to agriculture in the Northern Midlands has primarily occurred on valley fill, valley slopes and colluvial fill. This land use history results in a spatial bias in our analysis of factors driving the spatial distribution of greenspots. Nevertheless, our analysis of the environmental correlates of greenspots indicates that soil texture, soil depth, topographic slope and topographic position interact with regional climatic conditions. These interactions modify light, temperature and soil moisture regimes at a topographic scale in ways that enhance plant soil water availability, a key constraint on rates of photosynthesis. Thus the spatial distribution of ecosystem greenspots is to some extent spatially fixed by intrinsic variables. Given that precipitation deficit and surplus conditions, relative to the 1961–1990 baseline period, are projected to occur more frequently over much of Tasmania over the twenty-first century (A2 scenario) (White et al. 2010), intrinsic variables that can buffer sites from moisture stress will be potentially more important in the future.

The concepts of spatial insurance, conservation capacity and vulnerability to climate change that have been discussed in relation to conservation planning (Gillson et al. 2013), rely on a largely spatial approach to conservation. An important feature of the ecosystem greenspots method is that it provides information about temporal dynamics in vegetation productivity. The ability to identify locations that exhibit low temporal variation in vegetation productivity during drought conditions is potentially a powerful conservation planning tool in the context of climate change. Given the projections of increasing variability in climatic conditions, locations that remain productive during drought conditions may function as drought micro-refuges. Supporting evidence of the importance of temporal variation in productivity to habitat quality comes from faunal demographic studies (Gunnarsson et al. 2005; Wiegand et al. 2008). We suggest that conservation planning should explicitly consider temporal variation in habitat quality, particularly when exploring options for mitigating the impacts on biodiversity of projected future climate.