Thermal buffering effect of alpine boulder field microhabitats in Australia: Implications for habitat management and conservation

Biological Conservation 180 (2014) 278–287

Contents lists available at ScienceDirect

Biological Conservation journal homepage: www.elsevier.com/locate/biocon

Thermal buffering effect of alpine boulder field microhabitats in Australia: Implications for habitat management and conservation Haijing Shi a,⇑, David Paull a, Zhongming Wen b, Linda Broome c a

School of Physical, Environmental and Mathematical Sciences, University of New South Wales Canberra, P.O. Box 7916, Canberra BC, ACT 2610, Australia Institute of Soil and Water Conservation, Northwest A & F University, Yangling, Shaanxi Province 712100, China c Office of Environment and Heritage, P.O. Box 733, Queanbeyan, NSW 2620, Australia b

a r t i c l e

i n f o

Article history: Received 5 June 2014 Received in revised form 8 October 2014 Accepted 11 October 2014

Keywords: Microhabitat Boulder fields Extreme temperatures Refuges Climate change

a b s t r a c t Identifying and protecting microhabitats with conservation value is becoming a high priority for conservationists under projected climate change. Information on the buffering effects of microhabitats and its major drivers is a prerequisite to implementing effective conservation. We collected hourly temperature data from 70 alpine–subalpine sites nested within 9 boulder field clusters in an area measuring approximately 60 km  30 km in New South Wales, Australia, over a period of more than 2 years. We studied the thermal buffering effect of the boulder fields and investigated the factors driving the effect using multilevel (i.e., sampling points nested within boulder field clusters) regression modeling. We found remarkable seasonal and spatial variations in the thermal buffering effect of boulder fields. Boulder fields buffered the surface temperature maxima by 2.91 °C at a depth of 50 cm and 4.39 °C at a depth of 100 cm, while they buffered the surface temperature minima by 0.54 °C at the depth of 50 cm and 1.36 °C at the depth of 100 cm, with temperature range reduced by 3.45 °C and 6.48 °C in warmer period and by 1.23 °C and 2.05 °C in colder period at two depths. We defined thermal buffering in temperature ranges as the thermal buffering capacity (TBC) as an exemplar for multilevel modeling. We found that vegetation cover and elevation affected the TBC at the point level, whereas the aspect and inclination of slopes affected the TBC at the cluster level. These findings are instructive for the protection of habitats, as high priority should be given to those habitats that offer effective TBC. Furthermore, the environmental factors identified to be driving TBC provide important input to the management of species habitats. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction With accelerating global climate change, the conservation values of microhabitats are being increasingly recognized by ecologists (Scheffers et al., 2014; Williams et al., 2008). Microhabitats can partially decouple from regional climatic conditions, and species can persist in situ as regional climates become less suitable (Ashcroft, 2010; Bennie et al., 2008; Keppel and Wardell-Johnson, 2012). Microhabitats thus serve as thermal refuges and help species actively avoid climatic extremes (Ashcroft et al., 2012; Shoo et al., 2010). During the past decade, a number of studies have considered the thermal buffering effect of microhabitats of tree hollows (Isaac et al., 2008), roost cavities (Sedgeley, 2001), tropical boulder fields (Shoo et al., 2010) and various microhabitats within primary rainforests (Scheffers et al., 2014). In these studies, the microhabitats buffered temperature extremes by 1.6–10 °C. ⇑ Corresponding author. Tel.: +86 29 87012871; fax: +86 29 87012210. E-mail address: [email protected] (H. Shi). http://dx.doi.org/10.1016/j.biocon.2014.10.019 0006-3207/Ó 2014 Elsevier Ltd. All rights reserved.

Considering this microhabitat buffering effect is essential for robust assessments of the vulnerability of species to regional climate changes and specific conservation strategies (Suggitt et al., 2011; Williams et al., 2008). However, we facing considerable challenges when moving from research to specific conservation practice. One challenge comes from the fact that thermal buffering differs from one habitat to another (Shenbrot et al., 2002; Suggitt et al., 2011). Even within a few meters, the differences in microhabitat temperature is still great due to fine-scale changes in landforms, substrates and vegetation (Geiger et al., 2003; Suggitt et al., 2011). Capturing these fine-scale changes proves challenging because it requires accurate, fine-scale, modeled climate surfaces based on a wide variety of climate-forcing factors (Ashcroft, 2010; Suggitt et al., 2011). In most of cases, however, the effect of fine-scale microclimatic variation is overlooked in bioclimatic analyses (Suggitt et al., 2011). In recent years, the use of miniaturized temperature data loggers has provided opportunities to investigate the thermal buffering of micro-

H. Shi et al. / Biological Conservation 180 (2014) 278–287

habitats, but most studies conducted to date have been very site specific due to the cost of regional-scale sampling. Beside the cost of regional-scale sampling, spatial and temporal autocorrelations inherent in ecological data are an issue (Araújo et al., 2005; Dormann et al., 2007). Usually, locations in close proximity are more similar than those farther apart; thus, spatial autocorrelation occurs when samples are not independent (Dormann et al., 2007). In non-spatial comparisons, degrees of freedom and corresponding estimates of significance are often highly inflated due to an unstated assumption that residuals are independent; thus, type I errors are inflated (Anselin, 2002; Dormann et al., 2007). Moreover, the factors driving microclimatic variation may operate at various spatial scales (Keppel et al., 2012; Suggitt et al., 2011). Topographic variations may act to a greater extent on a slope scale, whereas vegetation cover and other driving factors operate on a finer scale, such as at a sampling-point. Scale variation thus makes it difficult to model the relationships between microhabitat temperatures and driving factors. The lack of long term monitoring data is another challenge for thermal buffering studies. The sampling duration of many studies has often been short, for example only one night (Paclík and Weidinger, 2007), several weeks (Ardia et al., 2006; Isaac et al., 2008), or a few months (Suggitt et al., 2011). In many instances, it is important to know whether the decoupling of microclimate is a consequence of divergence of temperature minima, maxima or both as this may have different implications for different taxa (i.e., in instances where the primary challenge is to avoid very low temperatures, namely ‘critical thermal minima’, or very high temperatures, namely ‘critical thermal maxima’). The lack of long term data may prevent such analyses, and seasonal variations are still not clear for most microhabitats. To address these knowledge gaps, we selected a site in the Australian alpine zone to investigate the capability of boulder fields inhabited by endangered and vulnerable small mammal species to buffer temperature extremes. We monitored hourly temperatures at various depths (from surface to subsurface) in multiple boulder fields across an extensive portion of the Australian alpine zone over a period of more than two years and thus obtained a long-term, high-temporal-resolution, and geographically extensive data set for analysis. To address the uncertainty and autocorrelation issues, and to use as small a number of data loggers as possible, we devised a multilevel/hierarchical data sampling and modeling approach. We used normal quantile–quantile (Q–Q) plots and Shapiro-Wilk tests to analyze the homogeneity of residuals and Moran’s I coefficients to test the residual spatial autocorrelation. This data set addressed two specific objectives: (1) to quantify the thermal buffering of alpine boulder fields, and (2) to determine the environmental factors that drive the thermal buffering of alpine boulder fields on various spatial and temporal scales. These two objectives are significant for the conservation of biodiversity under conditions of climate change.

2. Research design 2.1. Alpine boulder field habitat We selected the Australia alpine zone to investigate the thermal buffering effect of boulder field microhabitat (Figs. A1–A2). This habitat has been identified as an important refuge (Green and Osborne, 2012) because it provides den and nest sites for a range of endemic small mammals. In particular, boulder fields modify the local thermal microclimate and shelter the faunal inhabitants from wind, precipitation and predators (Green and Osborne, 2012). Noteworthy species that rely on Australian Alpine boulder fields include the nationally endangered mountain pygmy-possum

279

Burramys parvus and the broad-toothed rat Mastacomys fuscus, which is listed as vulnerable in New South Wales. When snow lies on the ground, an extensive subnivean space is created by boulders through which these small mammals can move freely (Sanecki et al., 2006). Considerable variation has been found in boulder field characteristics, such as vegetation cover, rock size, cavity dimensions and rock structure (Broome et al., 2012). We focused on alpine–subalpine regions in New South Wales, Australia. The boundary of the study area was based on the mapped distribution of boulder fields at elevations of >1100 m, resulting in a study area measuring approximately 60 km (north– south) by 30 km (east–west). The 12 primary boulder field clusters within this area included those of Mt. Kosciuszko, Summit Road, Charlotte Pass, Mt. Townsend, Paralyser, Blue Cow, Whites River, Gungartan, Gungartan Pass, Gungarlin, Toolong Range and Rough Creek (Fig. 1). 2.2. Data collection 2.2.1. Sampling design The field observation program began by identifying an initial set of 645 boulder field units originally mapped by the NSW Office of Environment and Heritage, which were then spatially corrected and ground-truthed by the first author (Fig. 1). According to the nested structure of boulder field clusters and multilevel modeling theory (Gelman, 2007), a stratified sampling design was used given that the elevation (m above sea level), aspect (bearing to the sun, in degrees) and inclination (angle of the surface to the horizontal, in degrees) of hill slopes are three fundamental factors that determine a site’s microclimate, particularly in mountainous areas (Busing et al., 1993; Geiger et al., 2003; Jin et al., 2008; Stage and Salas, 2007). Because simple aspect and inclination values do not account for the effects of topographic shading and other shadows, the amount of incident solar radiation was used for the stratification in place of aspect and inclination. A digital elevation model (DEM) of the study area with 25 m pixel resolution was used to compute potential direct solar radiation (under clear sky conditions) using a commercially available GIS software product (Arc/ Info) with AML script. The details of this method are described by Kumar et al., (1997). Based on the mean values of elevation and solar radiation, the study area was divided into 4 sampling strata (Fig. 2). The boulder field clusters within each stratum were then identified using a digitized boulder field map. Within each stratum, boulder field clusters were randomly selected within the constraints of budget and accessibility. The number of sampling points at each selected boulder field cluster was determined according to the size of the cluster. In total, we selected 91 sampling points in 38 boulder field units in 12 boulder field clusters distributed throughout the study area (Fig. 2). 2.2.2. Temperature sampling We collected temperature data at all 91 sampling points from November 2009 to February 2012 using 250 i-buttons (DS1922L Thermochron) and 30 Tinytag temperature sensors (TG-4100). The calibrations of all of these sensors were tested in the laboratory by immersion in an ice-water slurry prior to field deployment. The sensors were programmed to start at the same time, and each synchronously logged the hourly temperature for the duration of the sampling period. At each sampling point, three loggers were tied to a steel wire to record temperatures at the ground surface and at subsurface depths of 50 cm and 100 cm. The data loggers were covered by a white polyvinyl chloride (PVC) tube to prevent direct solar radiation effects. The data were downloaded in situ every 6 months (November and March, i.e., when the boulder fields were free of snow), to ensure that the sensors contained sufficient reserve memory to continue recording data. During the tempera-

280

H. Shi et al. / Biological Conservation 180 (2014) 278–287

Fig. 1. Contour map of the study area and the distribution of boulder field habitats. Boulder field units and clusters are delineated by digitized polygons. Original data sourced from the NSW Office of Environment and Heritage. Boulder field clusters consisted of those at Mt. Kosciuszko (2162 m), Summit Road (1894 m), Charlotte Pass (1774 m), Mt. Townsend (1955 m), Paralyser (1813 m), Blue Cow (1838 m), Whites River (1685 m), Gungartan (1920 m), Gungartan Pass (1933 m), Gungarlin (1214 m), Toolong Range (1644 m) and Rough Creek (1569 m).

ture sampling, some sensors were dismissed because the PVC tubes used to suspend the temperature loggers became loose and surface sensors were heated partially by the direct sun (i.e., surface sensors identified by abnormally high readings at Rough Creek, Toolong Range and Gungartan). Finally, we got the quality data from 70 alpine–subalpine sites nested within 9 boulder field clusters. 2.2.3. Environmental data sampling At each temperature sampling point, we characterized topographic variables, vegetation cover and rock structure. Topographic variables derived from a 25 m resolution digital elevation model (DEM) included elevation, aspect and inclination. Vegetation cover was measured by visual estimation (Klimeš, 2003) and supported by a library of photos obtained at each sampling point from November 21–29, 2009. The photographs were categorized in increments of 10% (0%, 10%, 20% ... 100% vegetation cover) with the assistance of five people who were professional ecologists or PhD candidates with at least three years of experience in vegetation and habitat studies. The rock structure of the boulder fields was characterized based on three environmental variables: rock layers, rock sizes and cavity sizes. The number of rock layers was measured in accordance with the method of Broome et al. (2005). Rock layers were counted during fieldwork based on the following categories: 0 = rock protruding from soil, 1 = one layer, rock sitting on soil; 2 = two layers, substantial rocks sitting on other rocks; 3 = three layers or more

than 2 layers but ground still visible; and 4 = more than 3 layers, ground not visible. The size of rock cavities was assigned to one of 5 categories based on the cavity diameter and the interconnection of tunnels among the rocks. Rock cavities typically had tunnels connected to one another, indicating that no rock cavity was completely closed. Large cavities, typically formed by large rocks, had long, wide tunnels, whereas small rock cavities had short, narrow tunnels. To quantify the rock cavity structure, the cavity diameters were classified as <30 cm, 30–60 cm, 60–120 cm and >120 m. If the cavity diameter was <30 cm and the cavity was closed at the bottom, it was assigned a value of 5. If the cavity diameter was greater than 120 cm and the cavity was open at the bottom, it was assigned a value of 1. The rock sizes were calculated using cubature formulae based on the profiles of the records. At each temperature sampling point, all visible rocks that surrounded a cavity were measured in three dimensions so that the rock size could be calculated. Snow cover duration during the cold months was measured at one point on the surface of each boulder field unit using a Hobo light sensor (UA-002-64-HOBO pendant, temperature/light data logger; http://www.onetemp.com.au/) based on the assumption that this point would be representative of the entire boulder field unit, provided that inclination, aspect and incident solar radiation receipt were consistent across the boulder field unit. The methodology of inferring the presence of snow cover from light attenuation measurements was validated by placing additional light sensors at sites monitored manually for snow depth by Snowy Hydro Limited.

H. Shi et al. / Biological Conservation 180 (2014) 278–287

281

Fig. 2. Sampling strata and selected sampling points. Yellow dots denote sampling points where temperatures were measured; red stars denote sites where snow cover was measured. Four sampling strata were defined based on incident solar radiation and elevation (Stratum 1: radiation <700  107 J/m2, elevation <1600 m; Stratum 2: radiation <700  107 J/m2, elevation P1600 m; Stratum 3: radiation P700  107 J/m2, elevation <1600 m; Stratum 4: radiation P700  107 J/m2, elevation P1600 m). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2.3. Data analysis 2.3.1. Seasonal patterns of temperature Prior to data analysis, we divided each year of monitoring into four periods based on gross patterns of temperature variation and snow cover (Fig. 3). During the summer (Period A), the boulder fields were free of snow and temperatures within the cavities were dynamic, that is, air flowed through them as part of the general air circulation in contact with the outside atmosphere (Geiger et al., 2003; Green, 2010). When light snow covered the boulder fields (Period B), temperatures among the boulders were not greatly buffered by the partially insulating snow and, thus, varied as a function of the ambient air temperature (Green and Osborne, 1998). Extremely low temperatures (below 0 °C) occurred during this period. During mid-winter (Period C), deep, persistent snow blanketed the boulder fields and decoupled the air space in the boulder interstices from the ambient air (Green and Osborne, 1998); the subsurface air then reached a constant temperature of approximately 0 °C. When snow started to melt during Period D, without sufficient insulating snow, widely fluctuating and low temperatures below 0 °C were again measured. Because the temperatures within all of the boulder fields at all of the measured depths were stable at approximately 0 °C during Period C, we only focused on calculating and modeling temperature during Periods A, B and D.

2.3.2. Quantification of thermal buffering effect According to previous studies, the thermal buffering effect of microhabitats can be measured by the differences in temperature maxima, minima or temperature ranges. In most cases, only the differences in temperature maxima or minima have been used because temperature extremes have been the focus of most studies. However, a stable microclimate is also important because dramatic temperature changes may be more harmful for species. In this study, we quantified thermal buffering by the differences in temperature maxima (for Period A) and minima (for Period B and D) as well as temperature ranges different seasons. For habitats like boulder fields in the Australian alpine region, surface temperatures can be as extreme as those of a hard, artificially paved surface because rocks absorb and emit radiation better than do soils or vegetation and they do not absorb and transpire moisture in the same way as soils and vegetation. Surface temperature changes may, therefore, be non-uniform with temperature changes within the boulder fields due to a time lag, which is an importance aspect that need to be addressed in data analysis (Scheffers et al., 2014). We tested this non-uniformity of temperature using the method of Scheffers et al. (2014). It is the temperature of a substrate rather than the general overhead air temperature that has a direct impact on the activities of ground-dwelling and subsurface species (Geiger et al., 2003). When quantify the thermal buffering effect of different microhab-

282

H. Shi et al. / Biological Conservation 180 (2014) 278–287

Fig. 3. Temporal variations in the temperature of alpine boulder fields at three depths (surface, 50 cm, and 100 cm) and annual snow cover variations during the period of 2009–2012. During Period A, the boulder fields were free of snow; during Period B, light snow covered the boulder fields; during Period C, persistent snow covered the boulder fields; during Period D, snow melting occurred. Thermal buffering capacities (TBCs) were calculated for Periods A, B, and D.

itats, previous studies have chosen various top layers as a baseline. For example, Shoo et al. (2010) measured near-surface (at a height of 0.7 m) temperature as a baseline for a boulder field on a tropical mountain in north-eastern Australia. Sedgeley (2001) measured ambient temperature below the canopy at 17 m above the ground to investigate the roost cavity microclimate in Fiordland National Park in New Zealand while Scheffers et al. (2014) measured air temperature using data loggers suspended above the rainforest canopy at Mt. Banahaw in southern Luzon in Philippines. In the Australia alpine zone, the surface of boulder fields is either bare rock or covered with prostrate vegetation. Therefore, we used the surface temperature of the rocks as a baseline to investigate how much the boulder field microhabitat can buffer at depth against surface temperature extremes. For the multilevel modeling and spatial comparison, we used thermal buffering in temperature ranges as an example due to the complexity of the modeling process. To provide a quantitative response variable, we defined the thermal buffering in temperature ranges as thermal buffering capacity (TBC) and quantified it by the following equation:

TBCijk ¼ TRijs  TRijk ;

ð1Þ

where TBCijk is the TBC at point i during period j at depth k, TRijs is the s surface temperature range at point i during period j, and TRijk is the temperature range at point i during period j at the depth k. 2.3.3. Multilevel modeling and data preparation Multilevel modeling is a statistically sound method for analyzing hierarchically structured data that allows for the introduction of explanatory (independent) variables at various scales and explicitly takes variability at various levels into account (Gelman, 2007; Overmars and Verburg, 2006). Multilevel models are regarded as extensions of regressions in which the data are structured in groups and coefficients can vary by group (Gelman, 2007). These models can be used to simultaneously estimate withingroup regression coefficients of predictor variables and the global means and variances of these variables among the groups (Jeffress et al., 2013), thus providing an efficient and flexible way to analyze hierarchically structured data. All of the analyses were performed in the R software environment (R Development Core Team, 2012). The software packages used included the multilevel package for data processing and multilevel modeling (Bliese, 2013), the nlme package for multilevel modeling (Pinheiro et al., 2011), the mgcv package for generalized additive modeling

(Wood, 2010), and the ape package for determining the Moran’s I coefficients (Paradis et al., 2004). All of the other analyses were conducted using the base package (R Development Core Team, 2012). The data used in multilevel modeling were arranged at two levels, i.e., at the point and cluster levels. The boulder field points in each cluster resembled one another in terms of their aspect to the sun and downslope inclination (Fig. 1). In the multilevel modeling, each level of the model was assigned its own matrix of potential predictors. At the sampling point level, the predictors of temperature primarily included vegetation cover, cavity size, rock layers and rock size (Table A1). Elevation, aspect and inclination of slopes were also associated with local effects at points within each cluster due to minor variations in the topography. Therefore, it was reasonable to incorporate these topographic variables at the point level. Aspect and inclination were combined using a potential solar insolation variable (Table A1) calculated as cos(aspect in radians)sin(inclination in radians) (Jeffress et al., 2013; Stage, 1976). This solar insolation index ranged from 1 to 1, with steep southfacing slopes being represented by large positive values (‘‘southness’’) and steep north-facing slopes represented by large negative values (‘‘northness’’). At the cluster level, the predictors included the ground surface elevation, aspect and inclination, which were aggregated from values determined at the sampling points within each cluster. 2.3.4. Multilevel model specification Before constructing the multilevel models, we fitted the variance component models of the TBC variables without substituting explanatory variables but, instead, by resorting to random intercepts (‘‘null’’ or ‘‘empty’’ variables), as in the following equations:

iTBCij ¼ c00 þ l0j þ cij ;

ð2Þ

VarðTBCij Þ ¼ r2l þ r2c ;

ð3Þ

where iTBCij is the TBC at boulder field point i in boulder field cluster j, c00 represents an overall average of the TBCs in all of the clusters (the ‘‘grand mean’’ c00). l0j is the difference between the average TBC of boulder field points in the jth cluster from that overall average. The value cij represents the difference between the TBC of a boulder field point and its cluster mean. In this model, rl2 denotes the TBC variation among the boulder field clusters and rc2 denotes the TBC variation between the sampling points.

283

H. Shi et al. / Biological Conservation 180 (2014) 278–287

Using Eqs. (2) and (3), it was a straightforward process to estimate the proportion of the overall variation in the TBC attributable to the clusters and to identify correlations between observations within the same cluster by calculating the intra-class correlation r2 (ICC), which was defined as r2 þlr2 (Gelman, 2007). Based on the varl c iance component models, it was a simple matter to introduce selected explanatory variables at various levels to construct multilevel models of the TBC. The form of the full multilevel model of the TBC, in the R software environment, is given by the equation

form < formulaðTBC ij  1 þ Cavityij  RockLayerij  RockSizeij þ Vegetationij  Insolarij  Elevationij þ G:elevation0j  G:Aspect 0j  G:Inclination0j Þ;

ð4Þ

We used the Akaike information criterion (AIC), Bayesian information criterion (BIC) and likelihood ratio test as criteria to select the best models (Bozdogan, 1987; Burham and Anderson, 2002). For the likelihood test, the model was refitted with maximum likelihood (ML) values instead of restricted maximum likelihood (REML) values because the nested fixed effects were expected to have the same random structure. The variables that had no significant effects were dropped from the final models. After the significant predictor variables were selected, the model was refitted with the REML values to obtain the optimal model corresponding to various temperatures. After fitting the optimal model with the various temperature variables, we used a normal Q–Q plot and Shapiro-Wilk normality test to inspect the homogeneity of the residuals cij of the optimal model. To test for residual spatial autocorrelation, we also estimated Moran’s I coefficients using the ape library in R software (Paradis et al., 2004) with a weights matrix of inverse distances. In the fitted models in the present study, no heterogeneity of the residuals cij or spatial autocorrelation were observed. A possible reason is that the correlation between sampling points within a cluster was accounted for in the multilevel models and the least distance among clusters was 1100 m, noting that distances >1000 m were given 0 weight, thus no spatial autocorrelation was observed. 3. Results 3.1. Thermal buffering of alpine boulder fields across the study area The temperature at the depths of 50 cm and 100 cm was on average cooler than surface temperature during Period A, and the temperature range was significantly smaller at the depths of 50 cm and 100 cm. Boulder fields buffered the surface temperature maxima by 2.91 °C at a depth of 50 cm and 4.39 °C at a depth of

100 cm, and reduced the temperature range by 3.45 °C at a depth of 50 cm and 6.48 °C at a depth of 100 cm (Table 1). In Period B, the temperature at the depths of 50 and 100 cm was on average warmer than the surface temperature, and the temperature range was significantly smaller at the depths of 50 cm and 100 cm. The boulder fields buffered the surface temperature minima by 0.54 °C at the depth of 50 cm and 1.36 °C at the depth of 100 cm while they reduced the temperature range by 1.23 °C and 2.05 °C at the two depths. We found a similar trend in Period D. The temperature maxima and minima at the depths of 50 cm and 100 cm increased at a slower rate than surface temperature maxima and minima, indicating that the temperature changes between surface and different depths were nonuniform (Table 1). With a 1 °C increase in temperature maxima, the temperature at the depth of 50 cm increased by 0.68 °C (R2 = 0.81, P < 0.001, 0.68x + 2.25), and at 100 cm depth by 0.67 °C (R2 = 0.67, P < 0.001, 0.54x + 2.97). In Periods B and D, with a 1 °C increase in temperature, the temperature minima at the depth of 50 cm increased by 0.71 °C and 0.94 °C (R2 = 0.61, P < 0.001, 0.71x + 0.01; R2 = 0.99, P < 0.001, 0.94x + 0.16), respectively, and at 100 cm depth by 0.68 °C and 0.87 °C (R2 = 0.94, P < 0.001, 0.68x + 0.67; R2 = 0.98, P < 0.001, 0.87x + 0.57), respectively. The quantified variable TBC showed that the buffering effect differed among boulder fields whereby certain fields displayed greater buffering than others (Fig. 4). The boulder fields at Summit Road, Paralyser, Mt. Kosciuszko and Charlotte Pass appeared to be more effective than other boulder fields at buffering against surface temperatures during Period A, as the TBCs were found to be high in those locations. The boulder fields at Summit Road, Gungartan and Paralyser exhibited higher TBCs than others during Periods B and D. 3.2. Estimated variance components of thermal buffering capacity (TBC) The established null models yielded the variance components of TBCs at the point and cluster levels (Table 2). The variance at the point level (rc2) was much greater than that at the cluster level (rl2), meaning that the TBC exhibited higher local variations. However, when two depths were compared, the intra-cluster correlations (ICC) at 100 cm depth were much higher than the ICC at 50 cm depth. During Period A, for example, the ICC was 12.6% at 50-cm depth, whereas it was 28.5% at 100 cm depth (Table 2). Similar trends were also observed during Periods B and D. This higher ICC meant that the TBCs at a depth of 100 cm were less variable than those at a 50 cm depth within a cluster. In addition to changes with depth, the ICCs displayed strong seasonality. During Period A, the ICC at 100 cm depth was 28.50%, whereas it was 40.9% and

Table 1 The differences in temperature maxima, minima and ranges at different depths with one-way ANOVA; the temperature maxima and minima are represented by mean daily maxima and minima in different time periods (refer to Fig. 3). The temperature range is the range of average daily temperatures at different depths. The rate of change refers to the rate at which the maxima and minima temperatures change at depths of 50 cm and 100 cm with a 1 °C increase in ambient temperature. Periods and temperatures

Surface (°C)

0.5 m (°C)

1 m (°C)

F

p

Period A Maximum Range Rate of change (°C) (surface:depths)

16.11 ± 2.28 8.07 ± 2.30 –

13.20 ± 1.53 4.62 ± 1.72 1:0.68

11.72 ± 1.19 1.59 ± 0.91 1:0.54

117.1 243.8 –

<0.001⁄⁄⁄ <0.001⁄⁄⁄ –

Period B Minimum Range Rate of change (°C) (surface:depths)

2.13 ± 0.64 3.31 ± 1.52 –

1.59 ± 0.65 2.08 ± 0.82 1:0.71

0.77 ± 0.78 1.26 ± 0.69 1:68

67.69 64.93 –

<0.001⁄⁄⁄ <0.001⁄⁄⁄ –

Period D Minimum Range Rate of change (°C) (surface:depths)

1.77 ± 0.81 3.90 ± 2.08 –

1.99 ± 0.82 2.32 ± 1.28 1:0.94

2.01 ± 0.77 1.59 ± 0.91 1:0.87

3.22 42.79 –

<0.04⁄⁄⁄ <0.001⁄⁄⁄ –

284

H. Shi et al. / Biological Conservation 180 (2014) 278–287 Table 2 Estimated variance components of the thermal buffering capacity (TBC) of boulder fields at two depths during three periods. The data in parentheses are the 95% CRIs (credible intervals) of the standard deviance. ICC is intra-class correlation (i.e., the correlation between sampling points within boulder field clusters). For the likelihood ratio test, p < 0.05 means that there is significant intercept variation between various clusters. Variables

Points rc2

Clusters rl2

ICC (%)

LRatio

p

TBCA50 TBCA100 TBCB50 TBCB100 TBCD50 TBCD100

2.69 (1.37, 1.96) 3.42 (1.55, 2.21) 1.07 (0.87, 1.23) 1.2 (0.92, 1.31) 1.02 (0.85, 1.20) 1.26 (0.94, 1.34)

0.39 (0.23, 1.67) 1.36 (0.59, 2.30) 0.34 (0.28,1.22) 0.83 (0.50, 1.65) 0.38 (0.32,1.19) 0.9 (0.53, 1.71)

12.60 28.50 24.30 40.90 27.20 41.70

2.41 8.48 6.03 16.03 10.5 20.18

0.12 <0.05 <0.05 <0.001 <0.05 <0.001

Note: In the TBC variables, A, B, and D represent Periods A, B, and D, respectively; 50 denotes the 50 cm depth, and 100 denotes the 100 cm depth within the boulder fields. For example, TBCA50 is the TBC at a depth of 50 cm in the boulder fields during Period A, and TBCA100 is the TBC at a depth of 100 cm depth in the boulder fields during Period A.

habitats at higher elevations exhibited higher thermal buffering capacities and thus constituted more stable microclimates. Multilevel models were fitted for other TBC variables (Table 3). The results indicated that elevation and vegetation cover were two factors determining the TBC at the point level, whereas aspect and inclination were two driving factors at the cluster level (Table 3). At the point level, elevation had a negative primary effect on the TBCB50, TBCB100 and TBCD50, whereas vegetation cover had a positive primary effect on the TBCA100, TBCB100 and TBCD100. At the cluster level, aspect positively affected the TBC during Period B (TBCB50 and TBCB100). With variations in aspect from south to north, the TBCB50 and TBCB100 increased. Inclination, like aspect, was a significant variable at the cluster level. Inclination evidently negatively affected the TBCD50 and TBCD100. With increasing inclination, the TBC decreased. The boulder field habitats on gentle slopes were found to be better at buffering against surface temperature variations.

4. Discussion 4.1. The thermal buffering effect of boulder fields in the Australian alpine zone

Fig. 4. Mean thermal buffering capacities (±standard errors) of boulder field clusters during three periods (Periods A, B and D). B C represents Blue Cow, C P represents Charlotte Pass, Gun represents Gungartan, M K represents Mt. Kosciuszko, Par represents Paralyser, S R represents Summit Road, G P represents Gungartan Pass, M T represents Mt. Townsend, and W R represents Whites River.

41.7%, respectively, during Periods B and D. This trend indicated that the differences in the TBCs among the clusters were greater during Periods B and D than during Period A. 3.3. Driving factors of thermal buffering capacities at various levels There was no significant random intercept effect in the TBCA50 (TBC at a depth of 50 cm during Period A). The TBCA50 was thus fitted with the linear equation

TBCA50 ¼ 11:67  0:0043  Elevationi þ ei ;

ð5Þ

where i = 1, 2, 3 . . . 70, ei  N (0, 1.622). Eq. (5) indicated that the TBC at 50 cm depth during Period A primarily related to variations in elevation. With a 1000 m increase in elevation, temperature range was reduced by 4.3 °C at the depth of 50 cm. The boulder field

Whether a microhabitat serves as a refuge from extreme climates depends on the extent to which it is decoupled from regional temperatures. Our study indicates that boulder field habitats in the alpine zone of Australia can be well buffered from their corresponding rock surface temperatures. For example, at a depth of 100 cm, the boulder field microhabitat was on average warmer by 1.36 °C during Period B and cooler by 4.39 °C during Period A than surface temperatures, and the temperatures became more stable (temperature ranges reduced by 2.05 °C in Period B and 6.48 °C in Period A). Although the depths of some boulder field cavities may exceed 500 cm (personal observation), and no extremely deep temperatures were measured due to the difficulty of deploying sensors amid the complex rock cavity morphology, our data suggested a greater thermal buffering effect would occur in deep cavities. These buffering effects must be taken into consideration when evaluating the vulnerability of species to climate change (Williams et al., 2008) and overlooking these fine-scale effects may overestimate the vulnerability of species to climate change (Luoto and Heikkinen, 2008; Randin et al., 2009; Scheffers et al., 2013). For example, Brereton et al. (1995) concluded that the bioclimatic range of B. parvus in Victoria will disappear with only a 1 °C rise in mean air temperature. However it may not be true because the temperature changes between surface and different

285

H. Shi et al. / Biological Conservation 180 (2014) 278–287

Table 3 Multilevel models of thermal buffering capacities (TBCs) at various depths and periods using the restricted maximum likelihood (REML) method. Standard errors are in parentheses.

* ** ***

Parameter

TBCA100

TBCB50

TBCB100

TBCD50

TBCD100

Fixed effect Intercept (c00) Vegetation cover Elevation G. aspect G. inclination

5.83 (0.48)*** 1.79 (0.60)** – – –

5.49 (2.37)* – 0.0030 (0.0014)* 0.35 (0.13)* –

6.35 (2.65)* 0.87 (0.35)* 0.0036 (0.0015)* 0.55 (0.15)** –

9.83 (1.67)*** – 0.0026 (0.0009)** – 0.23 (0.08)*

8.62 (1.75)*** 0.79 (0.37)* – – 0.44 (0.12)**

Variance components rc2 (point level) r2l (cluster level) % Explained deviance (point level) % Explained deviance (cluster level)

3.02 1.35 8.4 3.3

1.05 0.15 14.5 40.1

1.04 0.27 35.42 58.96

1.02 4.89e09 26.67 74.20

1.22 0.17 35.65 69.12

Metrics of model support AIC BIC 2LL

293.89 302.76 285.89

231.25 242.27 221.248

233.55 246.69 221.55

225.39 236.41 215.39

229.92 240.94 219.92

p < 0.05. p < 0.01. p < 0.001.

depths were nonuniform, and the temperature maxima and minima at the depths of 50 cm and 100 cm increased at a slower rate than surface temperature maxima and minima. The thermal buffering effect of the boulder fields displayed strong seasonal changes. The boulder field habitats, for example, had a higher buffering effect during summer, when extremely hot temperatures occurred (during Period A, the TBC was 6.48 °C) compared to other seasons. Species inhabiting various thermal niches may be threatened by extreme temperatures during different seasons. Therefore, these key seasons crucial to the survival of particular species must be identified. Erroneous predictions of extinction risk may occur if species distribution models are based on the wrong seasonal temperature values (Ashcroft et al., 2009). The buffering effect against temperature variations not only displayed strong seasonality but also differed among boulder field clusters, with some being higher than others. For example, the boulder fields at Summit Road, Paralyser, Mt. Kosciuszko and Charlotte Pass appeared to be more effective than others at buffering against surface temperatures during Period A (Fig. 4). The boulder fields at Summit Road, Gungartan and Paralyser exhibited higher TBCs than others during Periods B and D. As the buffering effect differed among boulder clusters it would be interesting to test whether the thermal buffering effect corresponds to density patterns in small mammals in the future studies. Moreover, with the predicted decline in snow cover in the Australian alpine region (Hennessy et al., 2003, 2008), the duration of persistent snow cover (Period C) may become progressively shorter with time. In this event, the thermal buffering effect of boulder fields during the times of year equivalent to Periods B and D will become particularly important. 4.2. Driving factors of the thermal buffering effect in boulder fields of the Australian alpine zone The factors that play important roles during certain seasons and at certain scales are closely related to the practical evaluations of refuges of alpine species. In Australian alpine boulder fields, vegetation cover, elevation, aspect and inclination were identified as the primary factors driving the thermal buffering effect, and these factors were found to act differently during various seasons and at various scales. Vegetation, which grows on and across the substrate (i.e., soil and rocks) and varies from one point to another, acts at the point level (Table 3). By reducing the surface temperature of rocks, the increase in vegetation cover may ameliorate temperature changes

by restricting air flow between the air mass above the vegetation cover and that below the rock surface (Geiger et al., 2003), thus reducing the temperature variation in the boulder field and creating a positive effect in terms of the thermal buffering of boulder fields. Our results indicated that elevation affects the thermal buffering effect at the point level, whereas aspect and inclination affect it at the cluster level (Table 3). Elevation exerted a negative influence on the thermal buffering effect at 50 cm depth during Periods B and D. This effect is likely due to elevation patterns that affect air temperatures in mountainous regions (Beniston, 2006). In general, because of the adiabatic lapse rate, the air temperature decreases with altitude (Laughlin, 1982; Rolland, 2003). Temperatures within boulder fields therefore tend to track the outside air temperature when they are free of snow (Green, 2010) and elevation thus has a negative effect on thermal buffering. Aspect and inclination act together at the cluster level because, within a cluster, they are essentially constant, as indicated by field observations; one cluster seldom encompasses different aspects or displays notable changes in inclination. In our study, aspect positively affected the thermal buffering effect of boulder fields at depths of 50 cm and 100 cm during Period B (Table 3). This effect is likely related to solar radiation in the Australian alpine region; in the southern hemisphere, north-facing slopes receive more solar radiation than south-facing slopes, resulting in more dramatic temperature changes on north-facing than south-facing slopes (Kirkpatrick et al., 1988). Inclination had a negative effect on thermal buffering at the depths of 50 cm and 100 cm during Period D. The boulder field habitats on gentle slopes are, therefore, better at buffering against surface air temperature changes. The specific reason for this negative relationship is unknown, but is possibly due to the effect that inclination has on air mass movements and wind velocity (Geiger et al., 2003), which merits future study. Rock structure variables, such as rock size, rock layer and rock cavity size, were expected to influence the thermal buffering effect, but they were not selected in the multilevel models. Boulder structure has recently been reported to affect the population densities of small mammals in the Australian alpine region (Broome et al., 2013). Therefore, the effect that rock structure exerts on microclimates in boulder field habitats also merits further study. 4.3. Conservation implications Alpine boulder fields generally create complex landscapes, and microclimate is driven by various factors at different scales. Iden-

286

H. Shi et al. / Biological Conservation 180 (2014) 278–287

tifying which factors drive the thermal buffering of microhabitats at particular scales and times provides valuable information for the assessment of vulnerability of species, and can contribute substantially to conservation planning. However, spatial autocorrelations inherent in these ecological data constitute a considerable issue (Araújo et al., 2005; Dormann et al., 2007). The multilevel modeling approach presented in this study may provide an option for similar studies in other regions, especially in mountainous areas. The findings of this study are also instructive for the protection of boulder field habitats in the Australia alpine zone under the influence of continuing climate change, as high priority should be given to the protection of habitats that exhibit the best thermal buffering effect. Understanding the factors that drive thermal buffering will provide important input to the management of habitats relied upon by threatened species. Acknowledgements We thank the University of New South Wales for funding travel expenses incurred during the fieldwork; the School of Physical, Environmental and Mathematical Sciences (UNSW, Canberra) for supplying maps, computing and library resources; and Julie Kesby (research officer) for editorial advice and support. We sincerely thank Hayley Bates for advice, encouragement and support during the fieldwork and for sharing her knowledge of boulder field habitats. We also thank the NSW Office of Environment and Heritage for additional logistical support and providing digital map layers of the study area, and we thank Snowy Hydro Limited for providing snow depth data. Three anonymous reviewers provided helpful comments on the manuscript, which were greatly appreciated. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.biocon.2014. 10.019. References Anselin, L., 2002. Under the hood issues in the specification and interpretation of spatial regression models. Agric. Econ. 27, 247–267. Araújo, M.B., Pearson, R.G., Thuiller, W., Erhard, M., 2005. Validation of species– climate impact models under climate change. Glob. Change Biol. 11, 1504– 1513. Ardia, D., Pérez, J., Clotfelter, E., 2006. Nest box orientation affects internal temperature and nest site selection by tree swallows. J. Field Ornithol. 77, 339–344. Ashcroft, M.B., 2010. Identifying refugia from climate change. J. Biogeogr. 37, 1407– 1413. Ashcroft, M.B., Chisholm, L.A., French, K.O., 2009. Climate change at the landscape scale: predicting fine-grained spatial heterogeneity in warming and potential refugia for vegetation. Glob. Change Biol. 15, 656–667. Ashcroft, M.B., Gollan, J.R., Warton, D.I., Ramp, D., 2012. A novel approach to quantify and locate potential microrefugia using topoclimate, climate stability, and isolation from the matrix. Glob. Change Biol. 18, 1866–1879. Beniston, M., 2006. Mountain weather and climate: a general overview and a focus on climatic change in the Alps. Hydrobiologia 562, 3–16. Bennie, J., Huntley, B., Wiltshire, A., Hill, M.O., Baxter, R., 2008. Slope, aspect and climate: spatially explicit and implicit models of topographic microclimate in chalk grassland. Ecol. Model. 216, 47–59. Bliese, P., 2013. Multilevel Modeling in R (2.5) – A Brief Introduction to R, The Multilevel Package and the NLME package. (viewed 10.04.13). Bozdogan, H., 1987. Model selection and Akaike’s information criterion (AIC): the general theory and its analytical extensions. Psychometrika 52, 345–370. Brereton, R., Bennett, S., Mansergh, I., 1995. Enhanced greenhouse climate change and its potential effect on selected fauna of south-eastern Australia: a trend analysis. Biol. Conserv. 72, 339–354. Broome, C., Dawson, M., Ford, F., Green, K., Little, D., McElhinney, N., 2005. Reassessment of Mountain Pygmy-possum Burramys parvus population size and distribution of habitat in Kosciuszko National Park. Biodiversity Conservation Section, Department of Environment and Conservation, Queanbeyan (unpublished report).

Broome, L., Archer, M., Bates, H., Shi, H., Geiser, F., McAllan, B., Heinze, D., Hand, S., Evans, T., Jackson, S., 2012. A brief review of the life history of, and threats to, Burramys parvus with a prehistory-based proposal for ensuring that it has a future. In: Lunney, D., Hutchings, P. (Eds.), Wildlife & Climate Change: Towards Robust Conservation Strategies for Australian Fauna. Royal Zoological Society of New South Wales, Mosman, NSW, Australia, pp. 114–126. Broome, L., Ford, F., Dawson, M., Green, K., Little, D., McElhinney, N., 2013. Reassessment of Mountain Pygmy-possum Burramys parvus population size and distribution of habitat in Kosciuszko National Park. Aust. Zool. 36, 381–403. Burham, K.P., Anderson, D.R., 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Springer-Verlag, New York. Busing, R., White, P., MacKenzie, M., 1993. Gradient analysis of old spruce-fir forests of the Great Smoky Mountains circa 1935. Can. J. Bot. 71, 951–958. Dormann, C.F., McPherson, J.M., Araújo, M.B., Bivand, R., Bolliger, J., Carl, G., Davies, R.G., Hirzel, A., Jetz, W., Daniel Kissling, W., Kühn, I., Ohlemüller, R., Peres-Neto, P.R., Reineking, B., Schröder, B., Schurr, F.M., Wilson, R., 2007. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. Ecography 30, 609–628. Geiger, R., Aron, R., Todhunter, P., 2003. The Climate Near the Ground, 6 edn. Rowman & Littlefield Publishers, New York. Gelman, A., 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press. Green, K., 2010. The aestivation sites of Bogong Moths, Agrotis Infusa (Boisduval) (Lepidoptera: Noctuidae), in the Snowy Mountains and the projected effects of climate change. Aust. Entomol. 37, 95–104. Green, K., Osborne, W.S., 1998. A winter niche: the subnivean space. In: Green, K. (Ed.), Snow: A Natural History, An Uncertain Future. Australian Alps Liaison Committee, Canberra, pp. 125–140. Green, K., Osborne, W., 2012. A Field Guide to Wildlife of the Australian Snow Country. New Holland Publishers, Chatswood. Hennessy, K., Whetton, P., Smith, I., Bathols, J., Hutchinson, M., Sharples, J., 2003. The Impact of Climate Change on Snow Conditions in Mainland Australia. CSIRO Atmospheric Research, Aspendale, Victoria, CSIRO Atmospheric Research, Aspendale, Victoria. Hennessy, K., Whetton, P., Walsh, K., Smith, I., Bathols, J., Hutchinson, M., Sharples, J., 2008. Climate change effects on snow conditions in mainland Australia and adaptation at ski resorts through snowmaking. Clim. Res. 35, 255–270. Isaac, J., De Gabriel, J., Goodman, B., 2008. Microclimate of daytime den sites in a tropical possum: Implications for the conservation of tropical arboreal marsupials. Anim. Conserv. 11, 281–287. Jeffress, M.R., Rodhouse, T.J., Ray, C., Wolff, S., Epps, C.W., 2013. The idiosyncrasies of place: geographic variation in the climate-distribution relationships of the American pika. Ecol. Appl. 23, 864–878. Jin, X., Zhang, Y., Schaepman, M., Clevers, J., Su, Z., 2008. Impact of elevation and aspect on the spatial distribution of vegetation in the Qilian mountain area with remote sending data. Int. Arch. Photogramm., Remote Sens. Spatial Inform. Sci. 37, 1385–1390. Keppel, G., Wardell-Johnson, G.W., 2012. Refugia: keys to climate change management. Glob. Change Biol. 18, 2389–2391. Keppel, G., Van Niel, K.P., Wardell-Johnson, G.W., Yates, C.J., Byrne, M., Mucina, L., Schut, A.G., Hopper, S.D., Franklin, S.E., 2012. Refugia: identifying and understanding safe havens for biodiversity under climate change. Glob. Ecol. Biogeogr. 21, 393–404. Kirkpatrick, J.B., Fensham, R.J., Nunez, M., Bowman, D.M.J.S., 1988. Vegetation– radiation relationships in the wet-dry tropics: granite hills in northern Australia. Vegetatio 76, 103–112. Klimeš, L., 2003. Scale-dependent variation in visual estimates of grassland plant cover. J. Veg. Sci. 14, 815–821. Kumar, L., Skidmore, A., Knowles, E., 1997. Modelling topographic variation in solar radiation in a GIS environment. Int. J. Geogr. Inform. Sci. 11, 475–497. Laughlin, G., 1982. Minimum temperature and lapse rate in complex terrain: influencing factors and prediction. Arch. Meteorol. Geophys. Bioclimatol. Ser. BTheoret. Appl. Climatol. 30, 141–152. Luoto, M., Heikkinen, R., 2008. Disregarding topographical heterogeneity biases species turnover assessments based on bioclimatic models. Glob. Change Biol. 14, 483–494. Overmars, K.P., Verburg, P.H., 2006. Multilevel modelling of land use from field to village level in the Philippines. Agric. Syst. 89, 435–456. Paclík, M., Weidinger, K., 2007. Microclimate of tree cavities during winter nightsimplications for roost site selection in birds. Int. J. Biometeorol. 51, 287–293. Paradis, E., Claude, J., Strimmer, K., 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics 20, 289–290. Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., 2011. Nlme: Linear and Nonlinear Mixed Effects Models. R Package Version 3.1-98, R Foundation for Statistical Computing, Vienna. . R Development Core Team, 2012. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, ISBN 3900051-07-0, http://www.R-project.org/. Randin, C.F., Engler, R., Normand, S., Zappa, M., Zimmermann, N.E., Pearman, P.B., Vittoz, P., Thuiller, W., Guisan, A., 2009. Climate change and plant distribution: local models predict high-elevation persistence. Glob. Change Biol. 15, 1557– 1569. Rolland, C., 2003. Spatial and seasonal variations of air temperature lapse rates in alpine regions. J. Clim. 16, 1032–1046.

H. Shi et al. / Biological Conservation 180 (2014) 278–287 Sanecki, G.M., Green, K., Wood, H., Lindenmayer, D., 2006. The implications of snowbased recreation for small mammals in the subnivean space in south-east Australia. Biol. Conserv. 129, 511–518. Scheffers, B.R., Brunner, R.M., Ramirez, S.D., Shoo, L.P., Diesmos, A., Williams, S.E., 2013. Thermal buffering of microhabitats is a critical factor mediating warming vulnerability of frogs in the Philippine biodiversity hotspot. Biotropica 45, 628– 635. Scheffers, B.R., Edwards, D.P., Diesmos, A., Williams, S.E., Evans, T.A., 2014. Microhabitats reduce animal’s exposure to climate extremes. Glob. Change Biol. 20, 495–503. Sedgeley, J., 2001. Quality of cavity microclimate as a factor influencing selection of maternity roosts by a tree-dwelling bat, Chalinolobus tuberculatus, in New Zealand. J. Appl. Ecol. 38, 425–438. Shenbrot, G., Krasnov, B., Khokhlova, I., Demidova, T., Fielden, L., 2002. Habitatdependent differences in architecture and microclimate of the burrows of Sundevall’s jird (Meriones crassus) (Rodentia: Gerbillinae) in the Negev Desert, Israel. J. Arid Environ. 51, 265–279.

287

Shoo, L., Storlie, C., Williams, Y., Williams, S., 2010. Potential for mountaintop boulder fields to buffer species against extreme heat stress under climate change. Int. J. Biometeorol. 54, 475–478. Stage, A.R., 1976. Notes: an expression for the effect of aspect, slope, and habitat type on tree growth. For. Sci. 22, 457–460. Stage, A.R., Salas, C., 2007. Interactions of elevation, aspect, and slope in models of forest species composition and productivity. For. Sci. 53, 486–492. Suggitt, A.J., Gillingham, P.K., Hill, J.K., Huntley, B., Kunin, W.E., Roy, D.B., Thomas, C.D., 2011. Habitat microclimates drive fine-scale variation in extreme temperatures. Oikos 120, 1–8. Williams, S., Shoo, L., Isaac, J., Hoffmann, A., Langham, G., 2008. Towards an integrated framework for assessing the vulnerability of species to climate change. PLoS Biol. 6, 2621–2626. Wood, S., 2010. Mgcv: GAMs with GCV/AIC/REML Smoothness Estimation and GAMMs by PQL. R Package Version: 1.6-2., .