Remote sensing shows south-east Queensland koalas (Phascolarctos cinereus) prefer areas of higher tree canopy height within their home ranges

Keywords: ecology, habitat preference, habitat use, koala, LiDAR, radio telemetry, spatial clustering, spatial ecology, wildlife management.

An animal’s home range is broadly defined as ‘that area traversed by the individual in its normal activities of food gathering, mating, and caring for young’ (Burt 1943). Powell and Mitchell (2012) proposed that a better concept of the home range is ‘that part of an animal’s cognitive map of its environment that it chooses to keep updated’. For animal ecology researchers, understanding the cognitive map requires integration of spatial and temporal animal location data along with the spatial distribution of resource availability and complexity (Heit et al. 2021).

VHF radio telemetry is routinely used to track fauna species, with location recorded at daily to weekly intervals using a hand-held GPS. Location data are then used to define home ranges using one of several methods, e.g. grid cell counts, minimum convex polygon, kernel density estimators, or low convex hull (Gregory 2017). More recently, GPS collars have been used to automatically record location data at much higher temporal frequency, e.g. hourly, but no habitat information is captured. High-frequency locations introduce spatial and temporal autocorrelation (Noonan et al. 2019), and location errors are also more likely with GPS collar data (Fleming et al. 2020). Approaches to reduce the influence of GPS errors on home range estimation include: (1) removing data points with lower precision, defined by the estimated horizontal dilution of precision (HDOP), (Adams et al. 2013); (2) calibration of collars at known locations so that an error term can be used in home range estimates (Ellis et al. 2011; Ironside et al. 2017); and (3) removing fixes with unlikely movement speeds (Wilson-Aggarwal et al. 2021). Many studies do not address GPS error (Meyer et al. 2021), and Ironside et al. (2017) concluded that screening data by HDOP resulted in a major data reduction and had few benefits, whereas Frair et al. (2010) found that home range estimates were robust to GPS error. Some home range packages specifically incorporate GPS error (Fleming et al. 2020). Home range studies may have other shortcomings, including: (1) the small number of individuals tracked; (2) the number of location fixes used to define home ranges; (3) insufficient study length to capture seasonal shifts in home ranges; and (4) indiscriminate use of the term ‘home range’, e.g. 86% of studies lasting less than 3 months referred to estimated areas as home ranges (Goldingay 2015).

Within home ranges ‘core areas’ have been defined as areas of more intensive use (Samuel et al. 1985), but for some species the concept of core home ranges is contentious because it may have little relevance to an animal’s home range use, e.g. northern bettongs (Bettongia tropica) utilise their home ranges evenly (Vernes and Pope 2001), and Bengsen et al. (2012) found that feral cats had no distinct core areas. Other studies have shown that core areas contain critical resources, e.g., an Indonesian study found that core areas contained all Bornean southern gibbon (Hylobates albibarbis) sleeping trees (Cheyne et al. 2019); in the pacific northwest (Washington, USA), female fishers (Pekania pennanti) select core areas dominated by larger trees (Lewis et al. 2016), and a study on spider monkeys (Ateles geoffroyi) found that, because core areas had better habitat quality than non-core areas, they were key to a wider understanding of habitat preferences (Asensio et al. 2012).

Most international papers using kernel density estimators have employed a 50% isopleth to define home range cores (Goldingay 2015), but other studies have used slightly different isopleths; e.g. Stirrat (2003) in their study of agile wallabies (Macropus agilis) used utilisation distribution to define the core area as the 55% isopleth. Using the same methods, Telfer and Griffiths (2006) adopted the 50% isopleth in their short-eared rock wallaby (Petrogale brachyotis) study; and Moseby et al. (2009) adopted a 60% isopleth for their study on feral cats and foxes. For Koala (Phascolarctos cinereus, henceforth koala) studies, Whisson et al. (2020) employed an animal-determined method (Powell 2000) to use 60% and 70% isopleths to define core areas within urban areas; de Oliveira et al. (2014), Goldingay and Dobner (2014), Matthews et al. (2016), Phillips (2016), and Rus et al. (2021) all used 50% isopleths.

The koala was listed as Endangered in New South Wales, Queensland and the Australian Capital Territory in February 2022 under the Environment Protection and Biodiversity Conservation Act 1999 (EPBC Act), with the major reason given as loss of habitat (DAWE 2022). Further research into how koalas select and use habitat is required, which can guide strategic planning and habitat management programs.

Knowledge of habitat use by koalas has so far relied on data collected at high resolution – with specific study aims guiding survey design, sampling intensity and scale. Data is typically collected from: (1) recorded koala locations (Gutteridge Haskins & Davey (GHD) 2009; Hanger et al. 2017); (2) plots (Phillips and Callaghan 2000); (3) transects (Hindell and Lee 1987; Dique et al. 2003); and (4) faecal pellet examination to determine food resources (Ellis et al. 1999). Current habitat management primarily relies on information derived from these data. For example, preferred tree species information obtained at plot scale is used to classify low-resolution, regional-scale mapped vegetation communities according to the proportion of preferred tree species within those communities (Lunney et al. 2000; Callaghan et al. 2011; Department of Environment and Science (DES) 2021a). Occasionally, other data, including landscape configuration and distance to roads, is incorporated into landscape or regional-scale studies (McAlpine et al. 2006; Januchowski et al. 2008; Crowther et al. 2014). Although these approaches are vital, regional-scale habitat maps may not capture habitat variability at more local scales, potentially leading to poor management decisions (Mitchell et al. 2021a).

Studying habitat preferences of individual koalas within their home ranges provides an additional approach to plot or landscape-scale studies. Koalas have high site fidelity (Thompson 2006), and home range sizes are often reported and used, e.g. for comparison within and between regions (Davies et al. 2013). Other data (e.g. tree species, tree height) are routinely collected with koala telemetry (Ellis et al. 2002; Janssen 2012; Davies et al. 2013; Goldingay and Dobner 2014; Matthews et al. 2016; Hanger et al. 2017). These data can also be used to assist with population-level inferences (Aarts et al. 2008) required for broader habitat management. For instance, Rhodes et al. (2005) showed that vegetation associations with a higher proportion of preferred tree species were used more often than other areas. Ellis et al. (2009) suggested that, although tree girth might guide koala foraging patterns, other factors including shelter tree availability and social networks also guide ranging behaviour. Foliar nitrogen availability is a strong driver for tree revisitation, including for isolated trees (Crowther et al. 2022). Crowther et al. (2014) also highlighted the importance of shelter trees for thermoregulation. Barth et al. (2020) found that isolated paddock trees provide ‘stepping stones’ between habitat patches and can be used for forage, particularly during the breeding season. In urban environments Goldingay and Dobner (2014) found that home ranges can require large areas of non-habitat to enable koalas to access patchy food resources. Davies et al. (2013) found that koala core areas in southwest Queensland were restricted to riparian habitats, and Goldingay and Dobner (2014) noted that core areas (50% fixed kernel) varied in proportion to home ranges (95% fixed kernel). Studies examining resource variability within home ranges are particularly rare, and we are aware of only one study (Gallahar et al. 2021) that assessed preferred tree species variability and tree size within home ranges. While providing useful information, few studies suggest how this knowledge can be specifically incorporated into habitat management (Goldingay 2015), exceptions being Crowther et al. (2014), who emphasised that shelter trees should be retained as well as food trees, and Gallahar et al. (2021), who recommended protection of remnant forest on higher-nutrient soils and retention of connecting corridors.

We acquired existing koala tracking datasets from three studies in south east Queensland covered by airborne LiDAR (Light Detection And Ranging) data: Coomera on the Gold Coast (14 koalas), and Redlands (10), both of which used Sirtrack koala-GPS collars (Ellis et al. 2016); and the Moreton Bay Rail Link (hereafter, MBR) research project, which tracked 471 koalas using a mixture of radiotelemetry locations recorded at least once every 3–4 days, and GPS collars (Hanger et al. 2017). After data examination (see below), we excluded some datasets and used locations from 135 koalas (Fig. 1) for our study.

Fig. 1.  Koala telemetry/GPS tracking datasets (solid black) in southeast Queensland used in this study (n = number of koala location datasets used in each area).

Coomera and Redlands koala tracking data were acquired between 2008 and 2010, contemporaneous with LiDAR acquisition in 2009. For MBR, we used the Historical Imagery function in Google Earth Pro (Google 2021) to discard data for individual koalas where any land clearing had occurred, either after LiDAR acquisition or during koala tracking data acquisition (2013–2017), i.e. we did not use any MBR data obtained with GPS collars that replaced telemetry after 2015.

For the MBR telemetry data, we rejected koalas with fewer than 80 locations (i.e. koalas that had been tracked for fewer than 6 months), and removed outliers which we defined as single locations more than 200 m from other locations. For Coomera and Redlands GPS data, to assess the influence of GPS accuracy on home range size estimation (95% and 50% isopleths), we generated several home ranges (see below) for one koala (Redlands K1009, n locations = 1359) by partitioning location data. Tracking data are commonly reported with HDOP (estimated horizontal dilution of precision, a measure of GPS satellite configuration). Firstly, we removed data with HDOP >10 (Adams et al. 2013), then removed data in two further stages (HDOP ≤ 5, HDOP ≤ 2.5). We repeated the process using only data from a minimum of four satellites, then used one daily (midday) location, and finally, split the dataset into two 4-week periods.

This exercise demonstrated the differences between home range estimators incorporating spatial and temporal autocorrelation (AC) locations, and those estimators treating locations as independent and identically distributed (IID). For AC methods, the difference between isopleths was smaller compared with IID methods and suggests that AC models handle location error better than IID models (Fig. 2). So, because the GPS data, and potentially, the telemetry data, were autocorrelated, we retained all data to minimise information loss. Fig. 2 shows that AC-generated home ranges were approximately 20% larger than IID. Splitting the data into two 4-week periods had the greatest effect on home range size – for the AC 95% isopleth, the home range was twice the size of the second 4-week period, and for the IID 50% isopleth, approximately 60% larger.

Fig. 2.  GPS collar data thinning techniques using tracking data from one Redlands koala (K1009) for (a) entire home range, and (b) core home range. We investigated improving data quality by: (1) removing locations with poorest horizontal dilution of precision (HDOP) in increments (1–4, 10–13, see key); (2) removing data from less than 4 satellites (5, 14); (3) retaining one daily location only (6, 15); and (4) splitting the 16-week tracking period into two equal periods (7, 8, 16, 17). Removing data based on HDOP has a minor effect on the estimates, with no evident trend. For AC, daily fixes (6) provide comparable results to HDOP thinning, but, for IID, daily fixes (15) increased the estimate size. The use of shorter tracking periods (7, 8, 16, 17) can have a marked effect on the estimate and widens confidence limits (black bars) except for (8). Models incorporating spatial and temporal autocorrelation (AC) generate home ranges approximately 20% greater, and with wider confidence intervals than models treating tracking data as independent and identically distributed (IID).

To determine home range and core areas, we used ctmm (Continuous-Time Movement Modelling; R Core Team 2021) with R ver. Ri186 4.1.1 (R Core Team 2015), RStudio ver. 1.4.1717 (RStudio Team 2021) and ctmm ver. 0.6.1 (Fleming and Calabrese 2017). As we illustrated above (Fig. 2), ctmm specifically addresses spatial and temporal autocorrelation arising from the use of modern GPS collars which, while recording multiple locations per day, may also include gaps in data collection and short sampling periods (Meyer et al. 2021). ctmm recognises that animal movement is a stochastic process and employs geostatistical methods to: (1) improve kernel estimation; (2) identify candidate kernel models; (3) identify the best model using maximum likelihood fitting with Akaike information criterion (AIC) selection; and (4) generates kernel density estimates at any desired isopleth (Calabrese et al. 2016). The outlie function highlights unusually large animal movements, i.e. locations more likely due to GPS error than actual movement, and we removed several of these from different datasets. We also used a visual diagnostic, the variogram function, to reject approximately 20 koalas with no defined home range (Fig. 3).

Fig. 3.  To remove koalas with no defined home range from analysis, we used the variogram function in ctmm. The first koala (a) has reached the variogram sill after less than 1 month (arrow), after which, area semi-variance is asymptotic, i.e. the home range is defined. The second koala (b) also reaches a sill after 1 month (arrow). After this point, the monotonic increase in semi-variance indicates no home range is defined, even after 3 years (Calabrese et al. 2016). Approximately 20 potential koala location datasets were rejected for analysis on this basis.

Airborne LiDAR uses laser pulses to measure the distance between the aircraft sensor and objects on the ground and, after classification, delineates ground cover, mid-canopy, and upper canopy features. We obtained LiDAR data from ELVIS (Elevation Information System) at https://elevation.fsdf.org.au/. This site collects data from many Australian government agencies, and we downloaded the 2 km × 2 km tiles coincident with the telemetry and GPS datasets. Following the methods of Mitchell et al. (2021b), we derived an upper-canopy height model (CHM) using FUSION software (McGaughey 2020). Although several trees might be within each 10-m cell, we extracted the maximum height within cells to characterise the overall canopy height.

The ‘Biodiversity status of pre clearing regional ecosystems – Queensland series’ (ver. 12.1; Department of Environment and Science (DES) 2021b) maps preclearing Regional Ecosystem (RE) vegetation communities in south east Queensland at a nominal scale of 1:25 000 (Neldner et al. 2019). Intersection of this dataset with the CHM captures isolated trees, small clumps of trees, and disturbed vegetation not included in the companion Remnant Vegetation dataset (Neldner et al. 2019). Within the 135 home ranges we examined, there were 225 RE patches with 20 distinct REs. RE tree species composition is summarised in Table 1.

Table 1.  RE (Regional Ecosystem) tree species composition within home range patches. Dominant species are listed first (Department of Environment and Science (DES) 2021c).

Within RE nomenclature, the second digit refers to the Land Zone, e.g. RE 12.3.6 is in Land Zone 3. Land Zones describe geology, landforms and geomorphic processes (Wilson and Taylor 2012), and provide a convenient means to group REs. In our study, these include Land Zone 3 (alluvial river and creek flats), Land Zone 5 (old loamy and sandy plains), Land Zone 11 (hills and lowlands on metamorphic rocks), plus an additional RE (12.9-10.4), which occurs on a combination of Land Zones (fine, and coarse-grained, sedimentary rocks).

We clipped the canopy height model to each 95% home range kernel. Canopy height varies among individual REs, so mean heights are not directly comparable among RE patches. To incorporate this variability into our analysis, we treated each patch individually, and following the methods of Mitchell et al. (2021b), we used Geoda (Anselin et al. 2006) with the Redcap agglomerative clustering algorithm (Guo 2008) to classify canopy height within each patch. We generated four height classes, i.e. highest, medium-high, medium-low, and lowest. The lowest class included areas with few or no trees. This approach allowed us to compare height classes between each RE patch rather than actual heights. Fig. 4 illustrates our approach for processing multiple REs within individual home ranges.

Fig. 4.  For home ranges within multiple Regional Ecosystems (REs), the CHM for each RE patch is extracted and classified to generate four height classes in each patch. In this example, there are three patches (A, B, C), each with two core areas (a), (b), which for each patch, are combined for analysis. Most home ranges (66) had one patch, 52 had two patches, and the remainder had three or more patches.

We extracted the Geoda-classified canopy height cells within the core area, with remaining cells constituting non-core classified cells. We tabulated data from 225 patches within 135 home ranges, and calculated the proportion of each class as a proportion of the total number of cells in core and non-core areas (Fig. 5).

Fig. 5.  Method used to compare canopy height class distribution between core (50% kernel) and non-core (remainder of 95% kernel) home range area for one MBR koala, ‘Shorty’. Core and non-core class proportions are calculated as a proportion of all cells (n) in each kernel. In this example, Class 1 comprises a slightly higher proportion of non-core (0.90) area than core (0.82) area. Class proportions in core and non-core areas for all home ranges, or home range patches, are then compared using paired t-tests. A Google Earth image contemporaneous with LiDAR acquisition (2009) of Shorty’s home range is also shown (image © Maxar Technologies 2021).

Paired t-tests were used to examine whether differences between core and non-core proportions were statistically significant for: (1) all patches; (2) patches grouped by Land Zone; and (3) patches within individual REs. Anticipating that some differences would not meet the standard t-test significance criterion (P = 0.05), we also used Cohen’s d, a secondary statistical measure that standardises effect size across datasets with different degrees of freedom (Eqn 1), where values between 0.3 and 0.5 are regarded as having medium effect size (Cohen 1988).

where M₁ = core class mean, M₂ = non-core class mean, , s.d. = standard deviation of the mean.

The entire data processing workflow is illustrated in Fig. 6.

Fig. 6.  Data processing workflow for one koala. The procedure is repeated for all koalas and tabulated into a final dataset for analysis using paired t-tests to compare core and non-core areas.

Collection of Coomera/Redlands koala location data (Ellis et al. 2016) was approved by the University of Queensland animal ethics committee (CMLR/937/08/ARC/RIOTINTO, ZOO/ENT/115/04/RT) and San Diego Zoo Global IACUC (no. 09-006), and conducted under permits from Queensland Department of Wildlife and Heritage (WITK05609808), WISP05609708, WISP00491303, WISP05609708). Collection of Moreton Bay Rail Link location data (Hanger et al. 2017) was approved by the Sunshine Coast (USC) Animal Ethics Committee (Animal ethics number AN/A/13/80) and by the Queensland Government (Scientific Purposes Permit, WISP11532912).

Mean home range size for our south east Queensland koala dataset was 13.9 ha (n = 135, median = 7.9 ha, range 0.8–113.5 ha). Home ranges on Land Zone 3 (alluvial soils) and Land Zone 5 (plains) appeared to be substantially smaller than those on Land Zone 11 (hillier metamorphic substrates); home range sizes in Land Zone 11 were also highly variable (Table 2). We did not test for significant differences among Land Zones because of unequal sample sizes, unequal variances, non-normal data distribution, and outliers.

Table 2.  Home range size (ha) varies by Land Zone.

Mean core home range area was 3.3 ha (n = 135, median 2.0 ha, range 0.19–22.5 ha), an average 23.8% of all home ranges. There were slight differences among Land Zones (range 23.4–26.4%), but none were significant (P_ANOVA = 0.48).

Class 1 (highest canopy) comprised 33.6% of entire home ranges, followed by Class 2 (27.1%) and Class 3 (21.0%), with Class 4 (lowest canopy height class) last at 18.5% (Fig. 7). Class 1 comprised 42.3% of core areas compared with 30.7% of non-core areas, Classes 2 and 3 were evenly distributed, and Class 4 occupied a greater proportion of non-core areas (20.3%) compared with core areas (11.0%).

Fig. 7.  Distribution of canopy height classes within 225 entire home range patches (grey), with core (50% kernel density estimate; blue) and non-core (remainder of 95% kernel density estimate; orange) components of the same patches. The mean proportion of Class 1 within entire patches is 0.34, but in the core patch component, Class 1 occupies a higher proportion of the patch area (0.42) compared with the non-core component (0.31). Classes 2 and 3 are fairly evenly distributed, and for Class 4 (lowest canopy height), the proportion is higher in non-core areas (0.20) than core areas (0.11).

Differences between core and non-core area proportions were significant for Class 1. Class 2 proportions were not significantly different, and Class 3 proportions were similar but still significantly different. Class 4 (which, as well as areas of lower canopy height, might include disturbed areas) differences were the most significant (Table 3). Variance (not shown) was similar in all classes, and Classes 1 and 4 had largest Cohen’s d.

Table 3.  For the combined dataset, canopy height Class 1 (highest canopy) is a significantly higher component of core areas compared with non-core areas (t-test for paired means).

The same trends were evident for RE patches grouped by Land Zones. Patches in Land Zone 11 (with lower-nutrient soils) have the highest proportion of Class 1 in core areas (50.2%) compared with the mean of other Land Zones, which have higher-nutrient soils (40.5%) (Fig. 8).

Fig. 8.  Distribution of height classes within Land Zones (geomorphic units). Height classes are shown within entire patches (grey), within core areas (blue), and non-core areas (orange). Land Zone 11 patches have the highest proportion of Class 1 (highest canopy height), both within entire patches (24% higher) and 25% higher than both the patch mean and core mean of other Land Zones combined.

Differences in Class 1 distribution were statistically significant for all Land Zones except 9–10, with a small sample size (d.f. = 9). Cohen’s d indicates a small-to-moderate effect for this Land Zone (Table 4). Generally, there were no significant differences for Classes 2 and 3, and for Class 4 there were significant differences across all Land Zones (Table 4).

Table 4.  Canopy height class distribution in core and non-core areas by Land Zone.

The same trends (Class 1 proportions higher in core areas, Classes 2 and 3 distributed equally, Class 4 more prominent in non-core areas) were graphically evident across individual REs, and these differences were significant for most REs, e.g. for RE 12.3.11/12.3.6/12.3.5 (Fig. 9, Table 5; see Supplementary material for all REs).

Fig. 9.  Height class distribution for RE 12.3.11/12.3.6/12.3.5. Class 1 occupies a greater proportion of core areas (49.6%) than non-core areas (37.9%), Class 2 is evenly distributed, and Classes 3 and 4 occupy a higher proportion of non-core areas.

Table 5.  For RE 12.3.11/12.3.6/12.3.5, there were significant differences in the distribution of Classes 1 and 4 in core and non-core areas, but not for Class 2, and Class 3 approaches significance.

The only exception in these trends was for RE 12.3.5, where there were no apparent (Fig. 10) or significant (Table 6) differences in class distribution between core and non-core areas.

Fig. 10.  Distribution of canopy height classes RE 12.3.5. Class distribution is similar for core and non-core areas. This is the only RE showing this trend.

Table 6.  For RE 12.3.5, class area proportions were very similar between core and non-core areas, with no significant differences among any height classes.

Our aim was to investigate whether koalas prefer forest areas with higher canopy height within their home ranges. We took a straightforward approach by firstly using spatial clustering to generate contiguous height classes from LiDAR within individual home ranges, and then examined differences in height class allocation between core and non-core home range areas.

Within our dataset of 225 Regional Ecosystem (REs – mapped forest communities) patches within 135 home ranges, we found that areas of highest canopy height (Class 1) comprised a significantly higher proportion of core areas (42.3%) compared with non-core areas (30.7%). For Classes 2 and 3 there was no significant difference between the two areas, but the lowest height class (Class 4) occupied 11.0% of core areas and 20.3% of non-core areas. Our results were consistent across both REs grouped by Land Zones (Wilson and Taylor 2012) and individual REs.

Our study examines use versus availability (Johnson 1980) of a potential resource within estimated koala home ranges. We achieved this by comparing forest canopy height distribution within non-core and core home ranges, where koalas can be expected to be found 50% of the time. Core areas, occupying 25% of the total home range, have a higher proportion of the highest-canopy class compared with non-core areas. We conclude that areas of higher canopy height are therefore likely to be an important habitat resource for koalas. Our results concur with previous findings that show koalas generally prefer larger or taller trees (Hindell and Lee 1987; Phillips and Callaghan 2000; Moore and Foley 2005; Matthews et al. 2007; Ellis et al. 2009; Smith et al. 2013; Office of Environment and Heritage (OEH) 2018).

Our results were similarly consistent across Land Zones, but we note that for Land Zone 11, the proportion of Class 1 (highest canopy) in core areas was 24% higher than the average of remaining Land Zones. This Land Zone has lower-fertility soils, therefore we suggest that these localised areas of higher canopy are more likely to be on soils with slightly higher nutrient availability, an important factor in habitat use (Januchowski et al. 2008; Law et al. 2017). Another explanation might be that core areas have a higher density of preferred food trees, e.g. E. tereticornis on alluvial soils (White 1999; Phillips et al. 2000; Office of Environment and Heritage (OEH) 2018). This species is generally taller than co-occurring tree species such as Corymbia spp. (Mitchell et al. 2021b) and in undisturbed areas, any patches with higher density of E. tereticornis might be captured by the highest canopy height class. We also note that Class 4 (lowest canopy height class) showed the greatest difference in use between core and non-core areas. Why koalas use non-core areas less is a complex issue – in some home ranges it might be because of habitat disturbance, in others because of lower nutrient status, local variations in tree species composition, or social factors. Nevertheless, these areas form a necessary component of home ranges. For individual REs, core area Class 1 proportions were also greater than non-core proportions. The only exception was RE 12.3.5, where we did not distinguish any significant difference between core and non-core areas for any height class. RE 12.3.5 is composed of Melaleuca quinquenervia open forest with Eucalyptus robusta and E. tereticornis as scattered individuals or clumps (Department of Environment and Science (DES) 2021c). These species are generally higher than the M. quinquenervia canopy (D.M., pers. obs), but do not have the density required to generate a distinct height class within home ranges.

By using spatial clustering instead of a continuous canopy height model, we have avoided two potential complications. Adjoining cells are more likely to belong to the same class, thereby reducing the influence of under-canopy GPS errors, which can be in the order of 30 m (Frair et al. 2010). Secondly, all spatial layers we used were classified (core, non-core, and four height classes), enabling us to use simple paired t-tests for our analysis.

Our research relied on existing koala tracking data, so we could not design a study capturing a more representative sample of koala habitats in southeast Queensland, or other regions with koalas, nor could we observe fine-scale aspects of forest structure and tree use by tracked koalas. We only examined the potential influence of canopy height class stratified by low-resolution pre-clearing Regional Ecosystems (REs) maps with a nominal 1:25 000 accuracy (Neldner et al. 2019); for smaller home ranges encompassing more than one RE, this limitation could introduce error into our analysis. We recognise that many other variables might influence core home range establishment, e.g. the proportion of preferred tree species within REs (McAlpine et al. 2006) or the availability of shelter trees (Crowther et al. 2014). Other evidence suggests that social interaction influences koala movement (Ellis et al. 2009; Jiang et al. 2022) and therefore core area extent, but we did not attempt to account for this variable because so many core home ranges overlapped. Core extents are likely influenced by both seasonal use of shelter trees (Ellis et al. 2009) and use of different tree species during breeding and non-breeding seasons (Martin 1985; Goldingay 2015), but few koalas were tracked for 12 months or longer so we could not capture seasonal home range shifts. For Coomera and Redlands, telemetry data was contemporaneous with LiDAR, but for the Moreton Bay Rail data this was not possible. We addressed this issue by using Google Earth historical imagery to reject data with habitat disturbance between LiDAR and telemetry acquisition.

Where koalas are present, the major factor in habitat occupancy is the density of preferred food and shelter tree species (McAlpine et al. 2006; Ellis et al. 2009; Davies et al. 2013), but low-resolution maps depicting these resources do not capture spatial variability at finer scales where koala habitat management decisions are made (Mitchell et al. 2021a). To our knowledge only one study has mapped all trees, and tree species, within an area occupied by koalas (Moore et al. 2010), which would be the ideal basis for habitat management, but across large areas this approach is impractical. Our study exploits another known koala habitat preference, the use of taller trees, and this information is available over large areas within the koala’s geographic range using LiDAR.

Higher soil fertility (Phillips et al. 2000; Callaghan et al. 2011), gentler slopes and lower elevation (Gutteridge Haskins & Davey (GHD) 2009), and landscape position (Callaghan et al. 2011) are other abiotic factors influencing koala habitat quality; one study in Redlands links all these (Biolink 2019). To a large degree these factors also determine the growth rate and size of trees (Bell and Williams 1997), as well as local variation in tree species diversity and forest structure (Neldner et al. 2019), which can confound the utility of current broad-scale habitat maps at the higher spatial resolution required for habitat management (Mitchell et al. 2021a). We surmise that areas of higher canopy, which can be mapped at high resolution, encapsulate some of these unmapped and highly variable abiotic factors that contribute to habitat quality, and hence habitat use by koalas.

Our calculated home ranges may be larger than comparable home ranges shown in the Moreton Bay Rail Link (MBR) study (Hanger et al. 2017) or other studies, for three reasons. Firstly, we used telemetry datasets with a minimum of 80 locations (i.e. a minimum of 6 months) compared with a minimum of 20 locations used in the MBR study (Hanger et al. 2017). Short-term studies may only capture part of a larger home range (Goldingay 2015); in our exploratory analyses, we detected temporary shifts (not shown) within larger home ranges for some koalas over a period of a few weeks using variogram plots within ctmm. Secondly, using GPS locations for one koala and all suitable kernel density estimators, we confirmed that estimators incorporating autocorrelation (AC) increase calculated home range size by approximately 20% compared with estimators assuming location independence (IID), similar to results obtained by Fleming et al. (2015). Lastly, many telemetry studies assume that, when data is collected daily or weekly, or by thinning GPS data (Matthews et al. 2016), then temporal autocorrelation need not be considered. This may be the case, but spatial autocorrelation may still occur. For example, we used daily-to-weekly telemetry data from 114 koalas tracked during the MBR study, and AC estimators for 45 koalas (39%) were ranked higher than IID estimators. Therefore, reported home range estimates relying on IID likely underestimate home range size.

A study in Lismore (New South Wales) by Goldingay and Dobner (2014) found that core areas were proportionate to home range areas, occupying an average 4.7% of home ranges. These home ranges included cleared and developed land, i.e. areas with no food or shelter resources. In our study, home ranges occupied relatively undisturbed habitats, albeit within peri-urban landscapes, and core areas occupied an average of 23.8% of all home ranges. This contrast highlights that consideration of core areas may be a useful addition to future koala ecology studies.

Our study has shown that higher canopy height is a koala habitat resource, but how important is it? We looked at canopy height in isolation, in a future study we will insert this resource into landscape-scale multivariate mixed-effect models (McAlpine et al. 2006; Januchowski et al. 2008; Crowther et al. 2014) to establish the relative importance of canopy height compared with known predictor variables.

Previous remote sensing research at home range scales has concentrated on image classification (for koalas, see Ream 2013; Wu et al. 2019), which does not produce contiguous map classes to present readily interpretable research findings, nor encourage subsequent use in habitat management maps. For the first time, we have used spatial clustering techniques with remotely sensed data to identify and quantify the variability of a resource within home ranges. Our methods demonstrate the utility of this approach and provide a template for future research. Ellis et al. (2009) established the importance of lower-height shelter trees to koalas, and Eyre (2006) identified forest structure as a resource for greater gliders (Petauroides volans); both resources could be mapped using the classified standard deviation of LiDAR-derived canopy heights. Crowther et al. (2014) identified sheltered gullies, taller trees, lower elevations, and canopy shelter as important koala habitat variables – these can all be mapped at any desired resolution using LiDAR. The clustering algorithm detected no differences for RE 12.3.5 (Melaleuca quinquenervia open forest with Eucalyptus robusta and E. tereticornis), a vegetation community important to koalas along the New South Wales north coast (Phillips et al. 2000), as well as in south east Queensland. These species are generally higher than the M. quinquenervia canopy (D.M., pers. obs), and can be considered as emergent trees requiring a different approach, e.g. using LiDAR and non-spatial clustering to identify areas with higher trees. Ultimately, combined remotely sensed data types could be integrated with existing habitat maps to form a resource model (Di Stefano et al. 2011) and, when extended to population level, provide the quantitative framework required to better inform fauna habitat management.

Core and non-core home ranges provide a convenient means to examine resource availability and resource use within home ranges. We used this approach to confirm that areas of higher canopy are indeed important to koalas. Relatively, core areas have 30% more of the highest canopy class compared with non-core areas. This is a useful addition to our knowledge about how koalas use their home ranges, and has wider implications for habitat management, e.g. Callaghan et al. (2011) suggested a height overlay map added to existing habitat maps would be beneficial for koala management. Our approach, using core and non-core areas, might be useful to examine other factors determining habitat quality.

Supplementary material is available online.

We are not authorised to redistribute koala location data, but our derived data (core and non-core home ranges, canopy height model etc.) can be shared upon reasonable request to the corresponding author. On request, we can provide the R scripts required to facilitate the use of ctmm for home range estimation.

The authors declare no conflicts of interest. Dave L. Mitchell is employed by the Australian Koala Foundation as a Landscape Ecologist.

We are grateful to the Australian Koala Foundation (AKF), who funded our research. The funder had no role in study design, data analysis, decision to publish, or preparation of the manuscript.

W. (Bill) Ellis, Research Fellow, School of Agriculture and Food Sciences, Faculty of Science, University of Queensland.