The construction of small‐scale, quasi‐mechanistic spatial models of insect energetics in habitat restoration: A case study of beetles in Western Australia

The management and restoration of ecological processes mediated by biotic interactions is now broadly advocated and may be achieved by targeting restoration towards key agents. Although theoretically examined, a practical approach to incorporating the physiology and energetics of insects into restoration planning is poorly articulated. I aimed to provide a case study using the thermal biology and energetics of beetles to identify the distribution of habitat suitability in a large restoration landscape.


| INTRODUC TI ON
Interactions between plants and other elements of the biome, including animals and fungi, are essential to many services that underpin plant community persistence and functional ecosystems (Perring et al., 2015). Habitat fragmentation, use of chemicals and proliferation of invasive species have led to a decline in native populations of many key taxa (Fahrig, 2003;Kearns, Inouye, & Waser, 1998) and a concurrent decline in the resilience and robustness of many ecological communities. Failure to consider fauna in habitat restoration planning and practice contributes to failures to achieve self-sustaining ecosystem dynamics due to the inability to reinstate fauna-mediated services (Cristescu, Frère, & Banks, 2012;Cross, Tomlinson, Craig, Dixon, & Bateman, 2019;Miller et al., 2017). In the southwest of Western Australia, as elsewhere in the world, insects play the dominant role in the provision of many services to ecological function (Keighery, 1980;Losey & Vaughan, 2006;Springett, 1978), but are amongst the most severely affected by landscape change (Ewers & Didham, 2006;Tylianakis, Didham, Bascompte, & Wardle, 2008).
Due to the obscurity and diversity of insects in Australia, the precise ecological relationships and even the insect taxa involved in these crucial interactions remain largely unknown (Mitchell, Williams, & Deesmond, 2002). Thus, the practicalities of restoring insect-mediated functions in disturbed habitats have been heavily debated (Dixon, 2009;Menz, Dixon, & Hobbs, 2013;Menz et al., 2011), given the paucity of information regarding the key threatening processes affecting insects, and lack of understanding of their resource requirements, microhabitat specificity and movement patterns in modified landscapes.
Recently, there has been strong advocacy for the value of incorporating thermal biology into ecological management, and specifically into remediation of landscape fragmentation (Tuff, Tuff, & Davies, 2016). One of the most advanced aspects of this advocacy argues that the basis of many constraints on animals in changing landscapes or changing climates could be understood by investigating their thermal physiology and energetics (Tomlinson et al., 2014). Tomlinson et al. (2014) provided a convincing theoretical construct, using a thermoenergetic model of key insect groups to empirically inform restoration planning. The practical approach to the construction of these models was, however, not well articulated. While the value of mechanistic models was advocated, it was not made clear which specific traits may be most informative, or most readily quantified through empirical data collection, nor how those traits could be transposed into a modelling environment. Several recent publications have described novel statistical techniques to estimate thermal performance of ectotherms using simple respirometry technology (Angilletta, 2006;Kovac, Käfer, Stabentheiner, & Costa, 2014;Kovac, Stabentheiner, Hetz, Petz, & Crailsheim, 2007;Marshall, Dong, McQuaid, & Williams, 2011;Tomlinson, Dixon, Didham, & Bradshaw, 2015;Tomlinson & Menz, 2015;Tomlinson & Phillips, 2015) and have made suggestions that these data could be extended to applied circumstances. Topoclimatic modelling has been explored as an accurate way of constructing small-scale, site-specific climate layers from which to construct model projections of management scenarios (Ashcroft & Gollan, 2012;Patsiou, Conti, Zimmermann, Theordoris, & Randin, 2014;Tomlinson, Webber, Bradshaw, Dixon, & Renton, 2017). Unfortunately, topoclimatic models are inherently data-referential, and substantial differences in local projections can result from small differences in the microclimatic modelling applied, or the extent of the environmental data subtending it. As yet, despite speculations by McCallum, McDougall, and Seymour (2013) and Tomlinson et al. (2014) that such an approach to understanding fauna-mediated service provision is possible with extant technology, a complete and standardized study of the spatial energetics of insect vectors has yet to be produced.
Here, I describe a complete ecological energetic approach to estimating the small-scale distributions of two species of scarab beetle in a highly fragmented, peri-urban landscape of Banksia woodland on the northern Swan Coastal Plain in Western Australia. I hypothesized that the species of scarab beetles would have different thermoenergetic constraints, implying different potential models of their smallscale distributions. I describe a process, starting with no specific data for either of these species, to construct mechanistic models of their physiological constraints. Further to this, I translate these physiological data to landscape-specific niche envelope models to estimate local distribution and potential barriers to movement in a restoration context. This is the first attempt of which I am aware to apply novel niche envelope modelling approaches (Kearney & Porter, 2004Kearney, Simpson, Raubenheimer, & Helmuth, 2010) to extremely high-resolution, site-specific projections. Finally, I use these models to estimate quantifiable energetic constraints of adult beetles and to provide empirical guidelines to accelerate their return and dispersal throughout the restoration trajectory.

| MATERIAL S AND ME THODS
The methods described here follow a distribution modelling approach projected at very high resolution. The basis to such distribution modelling is the interaction of an organism model, in this case informed by the thermal performance of metabolic rates, with an environmental model. The resulting values are then used to estimate the suitability of each location through time. A workflow diagram of how I approached this process is included in Figure 1.

| Study location, animal source and maintenance
The Gnangara Mound defines a large, elevated area of sand north of Perth, Western Australia (Figure 2), subtended by an aquifer which is currently the single most important source of potable water for the city. Although the native vegetation is predominantly Banksia woodland (Mitchell et al., 2002), which is a nationally protected community, the area was extensively clear-felled in the late 1920s for the establishment of commercial pine plantations (Finn, Stock, & Valentine, 2009;Perry, 1948), the removal and restoration of which are planned by 2028 (Finn et al., 2009). Much of the rest of the region is semirural, with suburban residential estates developed in the early 1990s (Mitchell et al., 2002). Approximately 25% of the area of this study has also been exploited for commercial extraction of construction sands since the 1980s, which has been associated with extensive research into restoration back to native vegetation (Benigno, Dixon, & Stevens, 2013;Rokich & Dixon, 2007;Turner et al., 2006). Many of the remnant native fragments in the area are regionally significant (Anonymous, 2000;Mitchell et al., 2002), and their restoration and reconnection are considered a substantial challenge for Western Australian restoration biology (Grose, 2013;Rokich & Dixon, 2007).
Beetles are the most diverse group of insects globally (Hart & Bale, 1997;Hutchinson, 1959), and this diversity is congruent with their fulfilling many key niches and providing many key ecosystem services, including bioturbation (Garcia & Niell, 1991), nutrient cycling (Cobb et al., 2010;Nichols et al., 2008) and pollination (Proctor, Yeo, and Lack (1996) ;Whitehead, Giliomee, and Rebelo (1987)). Two species of beetle were remarkably common during surveys of the study area. Didham et al. (unpublished) captured in ultraviolet-reflective vane traps large numbers of an unidentified species given the manuscript name Colpochila "species 2". Colpochila is the second-largest genus in the tribe Liparetrini, and unique to Australia, with approximately half the known species unique to or occurring in Western Australia (Britton, 1986). Additionally, Phyllococerus purpurascens was relatively under-represented in the long-term trapping programme (Didham et al., unpublished), but commonly responded to light trapping. Very little information is available on these groups beyond their taxonomic descriptions, so their ecological roles and reproductive biology remain undocumented. Both species are crepuscular or nocturnal.
Twenty Colpochila "species 2" and thirteen Phyllococerus purpurascens were collected using light traps (Bulbert, Gollan, Donnelly, & Wilkie, 2007), consisting of a fluorescent tube and two blacklight tubes suspended in front of a sheet of white cotton, established at remnant native vegetation blocks across the target landscape and returned the UWA Crawley campus for measurement on the night of capture. Each individual was exposed to a randomly assorted repeated-measures design of seven ambient temperature treatments ranging between five and 40°C. The beetles were measured only once daily, and between-treatments were maintained in a Perspex tank (25 cm height × 15 cm width × 30 cm length) kept near windows in a laboratory room maintained at 23.5 ± 0.1°C. The beetles were kept on a bed of damp local soil and fed an ad libitum diet of commercially available plum jam (IXL, Australia) and water.

| Flow-through respirometry and thermal performance
A three-channel flow-through respirometry system was constructed following the guidelines of Withers (2001). Compressed air was passed through a mass flow controller (Aarlborg DFC-17, USA) at a F I G U R E 1 A schematic of the workflow used in this study. The habitat suitability model results from the interaction between an animal model characterized by the thermal performance of metabolic rates with an environment model resulting from quasi-mechanistic downscaling of a global climate model. Habitat suitability layers are generated by repeated iteration of these functions over each cell of a rasterized landscape rate of 50 ml/min (ATPD). The incurrent air was passed into a cylindrical, 5-ml glass chamber or a 20-ml plastic chamber, depending on the size of the species, containing a bed of soda lime (CaOH, Sigma chemicals) separated from the beetle by a plug of cotton wool. Following findings by White, Portugal, Martin, and Butler (2006) that Drierite (anhydrous calcium sulphate, W.A. Hammond Company Ltd) is an appropriate desiccant for steady-state measurements, the excurrent airstream was dried with Drierite and passed through Qubit S151 infra-red CO 2 analysers (Qubit systems Inc). The Qubits were calibrated and checked for linearity using a three-point calibration with N 2 gas (0 ppm), air (350 ppm) and a calibration gas mix (1,500 ppm; BOC Gases, Australia). Each respirometry system was maintained at constant nominal temperatures (T a ) of 5, 15, 25, 30, 35, 38 and 40°C, and actual temperatures were measured throughout the respirometry trials using a Vaisala HMI 33 probe (Vaisala Oyj). These temperatures were selected to span the range of environmental conditions that the beetles might experience naturally, with a specific focus around the temperature range where a pilot study suggested that the beetles show thermal inhibition of metabolism. A blocked random-stratified design was used to assign T a , where blocks of three beetles were randomly assigned to T a treatments until each beetle had been exposed to all treatments.
Natural attrition meant that some individuals were not measured at all treatments, and some trials were removed from the data set when the beetles did not rest. In total, 37 C. "species 2" (25%) and 49 P. purpurascens (54%) measurements were lost or omitted (Table 1).

F I G U R E 2 (a)
The location of the Gnangara Mound on the Swan Coastal Plain north of Perth in Western Australia. The red polygon indicates the project area. (b) The area is characterized by remnants of Banksia woodland (dark green) surrounded by silvicultural plantations (pale green) of introduced softwoods and low-density residential development (grey). Yellow and orange areas indicate agricultural land use TA B L E 1 Average mass and coefficients (±1 SEM in parentheses) for the repeated measures nonlinear models of thermal responses of metabolic rate (VCO 2 ) and evaporative water loss (EWL), where y 0 is the y-intercept, k is the scaling exponent, and T d the temperature of divergence from the base pattern of exponential increase. For CO 2 ,V̇CO 2 T MMR is estimated from the first derivative of the fitted function, and errors represent the precision of estimation for both. N indicates the number of individuals measured, and n indicates the total number of measurements analysed In these cases, new individuals were substituted into the experimental design. Analogue data signals from all equipment were interfaced to a computer via a DataQ 710 data acquisition board (DataQ Instruments) polled every 10 s using a custom-written Visual Basic version 6.0 data acquisition program. All beetles were weighed before and after respirometry trials, and exposed to each T a for an hour, with an additional minimum 30-min baseline before and after the respirometry trials.

Variation in body mass of individual beetles between species and
T a treatment categories was tested using a linear mixed effects model (LME) with Gaussian error distribution, where individual beetle identities were held as a random factor to account for repeated measures. For this analysis, T a was treated as a categorical variable to identify any inadvertent differences in the body mass amongst treatment categories. Where no patterns in the residual distribution against the fitted values were found, indicating normality and heteroscedasticity, no transformations were applied to normalize the mass data prior to analysis. A model averaging approach was undertaken based on the second-order AIC (AICc) to find the most parsimonious model of the univariate and multivariate interactions of factors influencing body mass patterns after Burnham and Anderson (2002). Where mass was found to differ significantly between species, all metabolic measurements were allometrically corrected using mass 0.75 after Chown et al. (2007), and evaporative water loss (EWL) measurements were allometrically corrected using mass 0.67 after Edney (1977); Chown, Pistorius, and Scholtz (1998), and all interspecific comparisons made using these allometrically corrected values.
Respirometry data were analysed using a custom-written Visual Basic version 6.0 program, following the calculations articulated by Withers (2001). In order to estimate the resting metabolic rates (RMR) for each temperature, the lowest, stable, concurrent 10-min period of the respirometry measurements was used to calculate the RMR of the beetles. The different thermal performance curves of each species' metabolic rate (V̇CO 2 ) were modelled by nonlinear least squares regression, fitting a three-parameter biexponential function after Tomlinson and Menz (2015); Ayton, Tomlinson, Phillips, Dixon, and Withers (2016) and .
The theoretical basis to these types of relationships is broadly accepted (Angilletta, 2009;Angilletta, Niewiarowski, & Navas, 2002), and has been soundly articulated in relation to these functions (Angilletta, 2006;Tomlinson, 2019). The three-parameter biexponential curve articulates several physiologically significant effects of temperature, including the intercept of the curve at 0°C, y 0 , the metabolic scaling exponent k, describing the magnitude of the temperature effect on metabolic rate and T d , the temperature where performance begins to diverge from the base exponential pattern of increase. The function also allows the calculation of the lethal upper temperature, CT max , by solving for the point where the function is zero, and a sublethal point of maximal tolerance, T MMR , using the first order derivative following Fornberg and Sloan (1994) using the numDeriv package (Gilbert & Varadhan, 2012), and again solve for the upper point in the function where the instantaneous rate of change is zero using the "uniroot" function after Brent (1973) where precision is set to the precision of the temperature regulator (0.1°C).
The effects of species and T a on EWL were modelled using an exponential model (Ayton et al., 2016;Tomlinson & Menz, 2015). Unique thermal performance curves were fitted for each species to account for the repeated measures following Ritz and Streibig (2008), Peek, Russek-Cohen, Wait, and Forseth (2002) and Pinheiro and Bates (2000), which precludes statistical comparisons between the resulting models. Statistical analyses were performed in R version 3.0.3 using the nlme (Pinheiro, Bates, DebRoy, Sarkar, & R Core Team, 2014), numDeriv (Gilbert & Varadhan, 2012), AICcmodavg (Mazerolle, 2013) and generic packages. Data are presented as mean ± 1 standard error of measurement (SEM).

| Thermo-energetic landscape context
The study region has been subject to earlier microclimatic characterization , using a topoclimatic modelling approach following Ashcroft and Gollan (2012) which integrated local measurements from temperature loggers (El-USB1; Lascar Electronics) made at 1,024 randomly stratified locations across the study area. In an effort to standardize the statistical underpinnings of the climatic modelling, and to facilitate the production of distribution models subsequent to the microclimatic modelling, in this study I used the micro_global algorithm in the "NicheMapR" package (Kearney & Porter, 2017) to achieve this and tested the accuracy of the model using the data collected by the temperature loggers. The The resulting microclimatic model was used as the basis for the ectotherm algorithm to project species distributions and habitat suitability for each species, again in the "NicheMapR" package (Kearney & Porter, 2017). The ectotherm algorithm is essentially a biophysical model that solves thermal mass-balance equations to estimate heat, water and energy budgets on the basis of specified morphology (body shape), physiology (thermal tolerance thresholds) and behaviour (diurnality or nocturnality), and can be extended to the calculation of dynamic energy budgets across the whole life cycle (e.g. Kearney (2012); Kearney, Simpson, Raubenheimer, and Kooijman (2012)). The thermal performance curves were used to estimate key parameters for the ectotherm model, including the point where the performance function intersects the x-axis (CT max ), the point of peak metabolic activity (T MMR ) was used to parameterize the voluntary thermal maximum (VTMAX), and T d of each species was used as an of estimate T pref .
Basking temperature (T bask ) was arbitrarily set at 15°C for the two species, as was the temperature associated with minimum metabolic rates (VTMIN). Both species of beetle were estimated as ellipses (lometry = 2), and the species' mass was informed by the average mass of each species measured during metabolic trials. I used the resulting models to estimate hourly body temperatures (T b ), metabolic rates according to the biexponential thermal performance functions described above, which I standardized to comparable energetic equivalents (cal/hr) assuming a respiratory exchange ratio of one from a pure carbohydrate diet (Withers, 1992), and the activity of each beetle (0 = inactive, 1 = basking and 2 = active). Due to the substantial amount of data that resulted (approximately 6.22 × 10 10 bits) and the storage and data handling limits that this implied, the resulting data were av-

| Patterns in body mass
On average, C. "species 2" weighed 1,039.9 ± 15.22 mg and P. purpurascens weighed 145.2 ± 2.05 mg. Model averaging resolved only a single model of the variation in body mass, indicating significant differences in mean body mass between species (F 1,124 = 1,283.8, p = 2.20 × 10 -16 ), but no other significant effectors on variation in body mass of the beetles.

| Thermal performance of metabolic rate and evaporative water loss
The thermal performance curves of CO 2 production (V̇CO 2 ) to temperature (T a ) for each species were characteristically hump-shaped ( Figure 3a), with V̇CO 2 increasing exponentially up to a peak, and then declining rapidly. Colpochila "species 2" had the lower k and higher T d (Table 1, Figure 3a). The estimated T MMR (homologous with M T R of Tomlinson and Phillips (2012)) values reflected the same patterns between species, with C. "species 2" being the most thermally tolerant (Table 1, Figure 3a). There were exponential increases in evaporative water loss (EWL) with increasing T a for both species (Figure 3b). The interspecific patterns were congruent with the performance of metabolic rate, in that C. "species 2" had the higher y 0 (  Figure 3b), and was more robust to the effects of T a .

| Thermo-energetic extrapolations
There was a highly significant difference between the hourly esti- in comparison with the high requirements between 3.6 and 6.7 calories per hour required by C. "species 2." Much of the landscape, however, fell between approximately 2 and 20 calories per hour for the two species (Figure 4c,d). Notably, energetic requirements were heavily influenced by different types of vegetation cover, with pine plantations requiring the lowest energy expenditure (Table 2; Supplementary Material).
Despite the highly variable energetic requirement, the models predicted that both species were ubiquitously able to forage for approximately 50% of the time across the entire study area in summer (Figure 4e,

| D ISCUSS I ON
My hypotheses concerning the thermal performance of both beetles were supported to the extent that they displayed the same pattern of exponential increase in metabolic rate to a critical point. For P. purpurescens, this was followed by a rapid decline in metabolic rate, whereas for C. "species 2" this decline was not strongly evident, despite a levelling off of metabolic rate, and a departure from constant exponential increase with increasing T a . The metabolic rates for both species, however, supported a thermal performance model characterized by critical tolerance points. Furthermore, the critical points were different for each species, with the most naturally abundant species also being the most thermally tolerant. Indeed, C. "species 2" is more thermally tolerant than I expected during the design of my experiments, exceeding even the thermal tolerance of my equipment. As a result, the T MMR value resulting from the nonlinear fits here might be an artefact, and not accurately representative of real thermal tolerances of C. "species 2." If so, it is most likely that I underestimated the thermal tolerance of the species.
Underestimates of thermal tolerance are less damaging in a management and restoration context than overestimates, because they are likely to imply more rigorous management intervention than a species actually requires. What is also evident, however, is the strikingly high-thermal tolerance of both species, relative to their nocturnal habits, and the climatic conditions of the local area (see below

| The management utility of quasi-mechanistic projections
My second aim for this study was to provide a complete guide for one approach to constructing project-specific ecological F I G U R E 4 Thermo-energetic landscapes generated using the ectotherm algorithm within NicheMapR (Kearney & Porter, 2017) for the two scarab beetles studied during the austral summer. (a,b) Body temperatures were calculated as a daily average, incorporating both daylight and nocturnal estimates to account for the tendency of both species to retreat underground during the day. These do not exceed the tolerance limits of either of these species as defined by T MMR (Figure 3). (c,d) average standard metabolic rates projected by solving the thermal performance functions of each species for each projected T b , converted from V̇CO 2 to large-unit calories by assuming a respiratory exchange ratio of 1.
(e,f) Nocturnal habitat suitability for both species as a function of the amount of time that each species was projected to be "active" according to the conventions of NicheMapR at each location TA B L E 2 Seasonal average estimates of body temperature (T b ), metabolic rate (SMR) and proportion of time spent foraging in each habitat type by each beetle species in the project area.
Values are mean ± 1 SEM ×10 0.00 ± 0.00 × 10 +00 0.00 ± 0.00 × 10 +00 TA B L E 2 (Continued) energetic trajectories for environmental management (Tomlinson et al., 2014). Overall, the niche envelope models presented here provide useful and relatively accurate first-pass guidelines to the thermal tolerance limits and energetic requirements of a guild of insects that have not previously been studied in Australian systems. My approach also shows the relative simplicity with which such data can be assembled once an appropriately landscape-contextualized design has been constructed. The small-scale distributional models that I projected from metabolic thermal performance data are constructed on the basis of freely available, high-resolution spatial data. Since the microclimatic model reported here was constructed on the basis of a global model trained on biophysical properties without reference to a specific, local data set, the microclimatic projections developed here are less subject to the errors incorporated into the topoclimatic model previously produced for the study region . While the ambient temperatures projected for the region by the micro_global algorithm (Kearney & Porter, 2017) were significantly different from the measurements taken in the region , this significance is probably highly inflated by the large number of measurements and estimates available. The actual magnitude of these differences is relatively small, and unlikely to be highly biologically relevant for the beetles in question, especially at temperatures so far below their thermal tolerance thresholds. The intrinsic value of mechanistically based models, however, is their flexibility to the incorporation of new or more data Tomlinson et al., 2014). If the inaccuracy in the estimated environmental temperatures should propagate crucial thermo-energetic errors, the models presented here and by  these can be refined and redirected in the future as more accurate climatic data become available.
There is an emerging literature on the value of incorporating mechanistic data into the study of ecological management (Fedriani et al., 2018;Seebacher & Franklin, 2012 With increased scale of disturbance, and increasingly longsighted environmental management aspirations, there is increasing recognition that environmental variation and change make "success" of management programs unlikely, no matter how clearly goals are defined and how realistic they are (Weins & Hobbs, 2015).
Nevertheless, there has long been a recognition of the essential value of good hypotheses as the basis of adaptive management (Fontaine, 2011;McLain & Lee, 1996), and models such as the niche envelope models that I present here are essentially highly detailed, flexible hypotheses which are highly suitable to this purpose.
Furthermore, spatially explicit hypotheses in the form of niche envelope models incorporate mechanisms of change into empirically guided environmental management, which is also strongly advocated across a number of fora (Buckley et al., 2010;Helmuth, Kingsolver, & Carrington, 2005;Lavorel et al., 2015). In a rapidly changing environment, using historical conditions to set targets for restoration ensures that the targets will be missed (Weins & Hobbs, 2015), whereas small-scale, mechanistic or quasi-mechanistic modelling approaches can better incorporate temporal and spatial dynamics into restoration and conservation planning (congruent with the dynamic reference concept; Hiers et al. (2012) & Beyers, 2018). To be genuinely incorporated into an adaptive management programme, however, model hypotheses of any level of complexity remain in need of real and tangible testing. Thermoenergetic models of the type articulated here imply testing by measuring the energetics of free-ranging animals, which have previously indicated that unexpected patterns can emerge between the model and readlity (Tomlinson, Dixon, Didham, & Bradshaw, 2017).

| CON CLUS IONS
I conclude by reiterating arguments that the incorporation of wellfounded mechanistic models into project-level environmental management practice has powerful utility (Ashcroft & Gollan, 2012;Tomlinson et al., 2014), both in developing flexible, adaptive, empirical guidelines and in identifying novel vectors for overcoming management challenges. At project-specific size and resolution, the development of niche envelope models and other high-resolution projections provides substantial insight into interactions between ecology and environmental and landscape contexts (Ashcroft & Gollan, 2012;Tuff et al., 2016). These models can be rapidly expanded to provide flexible and empirical guidance to genuine adaptive management with appropriate and relatively simple, rapid and low-cost mechanistic characters  scientifically obscure, and with substantial insights as reward.
Furthermore, similar models have been applied to other ectotherms to direct management for public health benefits (Kearney, Porter, Williams, Ritchie, & Hoffmann, 2009) and conservation management (Kearney et al., 2008;Mitchell et al., 2013), and the thermal biology of many processes in plants is constrained by identical expectations (Rajapakshe, Turner, Cross, & Tomlinson, 2020), suggesting that the potential utility of modelling habitat suitability at high resolution is substantial.

ACK N OWLED G EM ENTS
The data underpinning this were the result of Australian Research