Study area

The study was conducted in Victoria, the southernmost state of mainland Australia (Fig. 1). We chose Victoria as our study system because in addition to a diverse reptile community it has a substantial gradient in both biophysical states and landscape structure; high levels of land clearing and urbanisation since European colonisation have resulted in marked landscape change. Victoria has a total land area of 227 496 km² and is Australia’s most densely populated state, with a total population of 6.3 million people (2017 figure) (Australian Bureau of Statistics 2018). Just over three-quarters (77%) of the population are concentrated in the state’s capital city of Melbourne (Australian Bureau of Statistics 2018). The total extent of native vegetation is estimated to have decreased by about half since European colonisation (i.e. pre 1750) (Eigenraam et al. 2013). Nevertheless, Victoria has a rich diversity of ecosystems including grasslands, wet and dry forest, coastal habitats, woodlands, and alpine and mallee communities (Cheal 2010). In most areas of the state the climate is temperate, with mild/warm-hot, dry summers (December to February) and cold, wet winters (Bureau of Meteorology 2016). The average minimum temperature for Victoria is 4.5°C in winter, with an average maximum of 27.4°C in summer. The annual average rainfall is 637 mm.

Fig. 1. 

Map showing the distribution of the five land cover types (cropping, grazing pasture, native vegetation, plantation forest and other) across Victoria, Australia.

Species occurrence data

We obtained reptile species occurrence records from the Victorian Biodiversity Atlas, a database that collates data of species observations across Victoria (State Government of Victoria 2018). We included records between 2005 and 2017 to increase the congruence between the sample data and the date of spatial layers (i.e. 2016/2017) used to represent landscape structure. Turtles and sea snakes were excluded from our dataset, as were any records that did not identify individuals by both the genus and species names. Records were also removed if they included descriptions suggesting that an individual may have been sighted outside its range (for example, records that were described as ‘beach washed’ or ‘found in predator scat’). To improve independence and spatial accuracy, data were thinned to one occurrence point per 100-m cell, and records with an accuracy of >100 m were excluded. Finally, we excluded occurrence points within 5 km of the state boundary where landscape structure could not be evaluated accurately (see below), resulting in a dataset of 5749 species occurrence records. Records with additional naming information were renamed to only include the genus and species (for example, ‘Liopholis whitii PLAIN BACK MORPH’ was renamed ‘Liopholis whitii’, and ‘Eulamprus tympanum marnieae’ and ‘Eulamprus tympanum tympanum’ were renamed ‘Eulamprus tympanum’).

Background data

To account for spatial sampling bias, we generated background locations using the target group sampling method. This involved using the species occurrence locations described above as background data, operating under the assumption that the same sampling biases will apply to taxonomically related species sampled using similar methods as those used to observe the focal species (Merow et al. 2013). Because the 5749 occurrence records that we extracted were fewer than the recommended number of background points for MaxEnt (Elith et al. 2011), an additional 5000 points were generated within the core distribution of the occurrence records. We achieved this by using the 5749 data points to calculate a 50% kernel density polygon and then selecting the additional points at random from within this space. This resulted in a background sample of approximately 10 000 locations that was used as background in all models. Kernel density estimation was conducted using the adehabitatHR package (Calenge 2006) in the R statistical environment (R Core Team 2020).

Predictor variables

Predictor variables were selected within two broad categories: landscape structure and biophysical variables.

Landscape structure variables

A map of land cover types for Victoria was obtained from the DataVic website (State of Victoria (Agriculture Victoria) 2018). The original layer included 15 land cover types, which we reclassified into five classes: (1) cropping; (2) grazing pasture; (3) native vegetation; (4) plantation forests; and (5) other (e.g. horticulture and urban areas) (Fig. 1). The layer was clipped to the study area boundary and converted to a raster layer with a pixel size of 250 × 250 m using ArcGIS Desktop 10.5. A moving window analysis of the land cover layer was performed using FRAGSTATS (ver. 4.2; McGarigal et al. 2012) to calculate five landscape structure metrics: (1) Contagion Index; (2) Euclidean nearest neighbour (mean); (3) Native Vegetation Cover; (4) Shape Index (arithmetic mean); and (5) Shannon’s Diversity Index (Table 1) (McGarigal 2015). We used a square window with a side length of 5 km, and accepted samples with a maximum 10% of border/no data for all metrics. This analysis calculates a specified metric within a 5-km window for each pixel defined by the base land cover layer, resulting in a spatial layer of the specified metric with the same pixel size (250 × 250 m). In preliminary analysis we calculated a range of metrics, but chose the five listed above because they were effective at characterising different aspects of landscape structure and were minimally correlated (Pearson’s correlation: −0.2 < r < 0.6).

Table 1. 

Description of landscape structure variables used in reptile species distribution modelling (Source: McGarigal 2015).

Species distribution modelling

We used MaxEnt to develop species distribution models in the R statistical environment (R Core Team 2020) using the dismo package (Hijmans et al. 2020). MaxEnt uses presence–background data to model species distributions, and performs well compared to other algorithms (Elith et al. 2006) and on datasets with a range of sample sizes (Wisz et al. 2008). We used MaxEnt to develop species distribution models for 40 reptile species with ≥30 occurrence records. We chose a minimum of 30 records as the cut off for analysis, because MaxEnt performs well for samples of this size or greater (Wisz et al. 2008).

Our aim was to compare the relative influence of different predictor variables on species distributions (rather than to identify the parsimonious model), so distribution models were generated for all 40 species using the same set of predictors, features and regularisation parameters. We used hinge features only to simplify the model while enabling the representation of non-linear responses, and chose the default regularisation value as a practical option for modelling multiple species (Merow et al. 2013). We selected the logistic output to allow for comparison of species responses to the different predictor variables (Merow et al. 2013). Area under the curve (AUC) values were used to assess the predictive ability of the models. AUC measures the area under the Receiver Operating Characteristic curve, and in the context of presence–background data represent the ability of a model to discriminate between species presences and background sample points. We followed the thresholds set out in Thuiller et al. (2003), where AUC values greater than 0.9 represent excellent predictive ability, 0.8–0.9 represent good predictive ability, 0.7–0.8 represent fair predictive ability, 0.6–0.7 represent poor predictive ability, and less than 0.6 represents models that perform little better than random. We used five-fold cross-validation to estimate the uncertainty associated with predictive ability, presenting the mean, minimum and maximum AUC values from the five runs. Maps of the relative likelihood of occurrence predicted by MaxEnt were produced for all 40 species in R using the raster (Hijmans 2020), tidyverse (Wickham et al. 2019) and cowplot (Wilke 2020) packages. Graphs representing species responses to individual predictors were produced in R using tidyverse. Relative likelihood of occurrence was represented on scale from 0 to 1 to allow for comparison between predictor variables and species.

Permutation importance values were used to assess the relative influence of the environmental variables to the species distribution models. Permutation importance is calculated by randomly permuting the values of each variable between the presence and background points, and measuring the resulting decrease in AUC. The bigger the decrease, the greater importance of that variable to the model. Permutation importance is considered a strong measure of variable influence because it does not depend on the pathway used to obtain the final model (Phillips 2017). We considered predictors influential if permutation importance was ≥5%.

Influence of landscape structure on species distributions

Our results indicate that landscape structure plays an important role in driving the distributions of some reptile species in Victoria. The extent of native vegetation was the most influential landscape structure variable, and was important for species within all seven reptile families included in the study. Extent of native vegetation cover was an important driver of the distributions of three-quarters of the species in this study. Most of these species had a positive response to native vegetation cover, such as Diplodactylus vittatus, which appears to occur more commonly in regions with warm, dry climates and high vegetation cover. Likewise, Ctenophorus fordi was found to have a strong positive relationship with extent of native vegetation, which may be due to their association with hummock grasses (Triodia sp.) (Verdon et al. 2020). In semiarid parts of Australia, many reptiles, including C. fordi, are associated with the extent cover of hummock grasses, which are a dominant ground cover and an important resource for thermoregulation and protection from predators (Verdon et al. 2020). Habitat availability is critical to the persistence of most species, and habitat loss and modification are a clear driver of reptile declines in eastern Australia (Driscoll 2004; Brown et al. 2008) and worldwide (Doherty et al. 2020). Most reptiles, and particularly skinks (which represent most species in this study), are ground dwelling and therefore sensitive to loss and disturbance of ground and understorey layers including leaf litter (Brown 2001).

Interestingly, several species had negative responses to the extent of vegetation cover. For example, D. impar, a small, threatened legless lizard (~300 mm long), preferred areas with less native vegetation cover, potentially reflecting this species’ reliance on threatened native grassland habitat, rather than a preference for low levels of vegetation cover per se (Kutt 1993). Since European colonisation, grasslands have suffered extensive loss and degradation, and now remain primarily in isolated remnants (Scroggie et al. 2019). Disturbances such as fire and grazing may promote the persistence of D. impar when applied in isolation by promoting structural and floristic diversity of grasslands and increasing connectivity (Scroggie et al. 2019). However, when applied simultaneously, fire and grazing were associated with a high risk of local extinction, and any further loss of native vegetation cover is likely to exacerbate the decline of D. impar (Scroggie et al. 2019).

Shape index was the most influential configuration metric considered in this study, and reflects the perimeter-to-area ratio and therefore the amount of edge environment. Four species responded positively to shape index, including P. textilis which was more likely to occur in irregularly shaped habitat patches and where there was an intermediate level of native vegetation cover. More irregular, elongated patches have a higher perimeter-to-area and are therefore more exposed to surrounding land cover types, which may advantage species that engage in landscape complementation and require resources in multiple land cover types (Nimmo et al. 2019). Pseudonaja textilis are a widely distributed species known to exploit agricultural areas where prey species (specifically Mus musculus) are abundant (Whitaker and Shine 2003). Similarly, for two species (N. scutatus and L. bougainvillii), relative likelihood of occurrence was higher when the contagion index was low (representing low clumping and high intermixing of patch types), yet had hump-shaped responses to vegetation cover extent, such that occurrence was more likely when vegetation cover was intermediate. This suggests that although vegetation cover is important, these species may be able to take advantage of resources associated with increased intermixing of different land use types.

In contrast, several species, including L. punctatovittata, were more likely to occur in more regular-shaped habitat. Irregular and elongated patches may also have undesirable implications for species that inhabit them, including increased exposure to introduced plants and animals (Brown et al. 2006) and disease (Jellinek et al. 2014; Uriostegui-Velarde et al. 2018). Negative relationships have been found between exotic plant cover and both the abundance and richness of lizard species (Jellinek et al. 2004). Because weed cover is expected to be higher in irregular shaped habitat, such areas may represent less suitable habitat for reptiles. Additionally, more elongated patches have a higher edge-to-area ratio; edges are subject to more to extreme climatic conditions and more extreme variations in temperature, which may impact the ability of reptiles to thermoregulate (Jellinek et al. 2004). Shape index has also been found to negatively influence plant species richness, possibly due to the greater availability of resources in more regular-shaped patches (Li et al. 2020), which may in turn afford more resources for animals.

Shannon’s Diversity Index was important to the distributions of three reptile species, but the directions of relationships were variable. Eulamprus tympanum responded negatively to increasing diversity, and positively to native vegetation cover. This is probably at least partly due to this species’ preference for shelter sites with logs (Langkilde et al. 2003), which are likely to be more common in native vegetation than other land cover types. In general, the presence of other land cover types is thought to negatively impact reptile species persistence, and a global meta-analysis found most types of habitat modification (including mining, grazing, agriculture, plantations and patch size reduction) have an overall negative impact on reptile species abundance (Doherty et al. 2020). In south-eastern Australia, agriculture is thought to negatively affect reptiles (Driscoll 2004), and pine plantations have been shown to act as a barrier to connectivity in some species and increase connectivity for others (Mortelliti et al. 2015). In our study, T. rugosa, a habitat generalist, was more likely to occur in areas with moderate levels of diversity and native vegetation cover. This species is known to favour open habitat and commonly found on roadsides, and may utilise both natural and artificial shelter (Norval and Gardner 2020), characteristics that may provide it with the flexibility to exploit multiple land cover types.

Although our study shows that native vegetation cover and configuration influences the distributions of some species at a statewide scale, it is important to note other reptiles may respond to aspects of landscape structure at finer spatial scales. Habitat cover and configuration, as well as other aspects of landscape structure such as the edge-to-area ratio, influence microclimate and the availability of shade, which is required for thermoregulation. South-eastern Australia is also subject to low-frequency, high-intensity fires (Murphy et al. 2013) that create mosaics of vegetation age classes in the landscape. Where species select different age classes, fire may represent an important additional source of habitat fragmentation, which is likely to influence distributions and threaten species persistence (Sitters and Di Stefano 2020). Furthermore, it is anticipated following future changes to the climate, many species worldwide will alter their present distributions (Bezeng et al. 2018). However, such range expansions will not be possible if the landscape structure does not facilitate this movement. Incorporating landscape structure, as well as dynamic predictor variables such as fire (Senior et al. 2021; Swan et al. 2021), into models may improve our ability to predict future species distributions and enable more informed management strategies.

Influence of biophysical variables on species distributions

We found climate variables, especially annual precipitation and the maximum temperature of the warmest month, were the strongest drivers of species distributions for almost all reptile species. Multiple studies have shown climate is a key factor determining both the current and future distributions of reptiles (Ihlow et al. 2012; Cabrelli et al. 2014; Hosseinzadeh et al. 2020) and their nesting habitat (Pike 2013). Environmental temperature plays a particularly important role in regulating the metabolic rates and activity patterns of ectotherms (Buckley et al. 2012), and regulates the reproductive success and sex determination in many species (Witt et al. 2010; Ihlow et al. 2012; Castelli et al. 2021).

Models of both reptile and amphibian species distributions have predicted that, on average, species are likely to experience increases in their climatic range due to climate change (Araújo et al. 2006; Hosseinzadeh et al. 2020). However, both groups are generally considered to be poor dispersers, and this factor, combined with extensive fragmentation, means that we should not assume that species have the ability and pathways to disperse and colonise new habitat (Araújo et al. 2006). Indeed, modelling of European amphibians and reptiles has shown that when models assume no dispersal, virtually all species are projected to experience range contractions (Araújo et al. 2006). Therefore, it is important that other non-climatic factors like dispersal (Sahlean et al. 2014; Rodrigues and Lima-Ribeiro 2018) and changes to landscape structure (Uriostegui-Velarde et al. 2018) are incorporated into models if we are to properly understand the factors driving current and future species distributions.

Soil texture was an important driver of the distributions of a fifth of the species in this study, including common species such as T. rugosa and L. guichenoti, which both showed an aversion to clay loams. Clay soils experience seasonal changes, becoming waterlogged during the wet season and cracked during the dry season (Woinarski et al. 1999). Although clay soils may provide an important temporary water source for frogs (Ferreira et al. 2018), this seasonality may be unsuitable for many reptiles (Woinarski et al. 1999). Notable differences have also been reported among reptile communities on sandy and clay soils that support different plant communities in arid and semiarid parts of Australia (Driscoll et al. 2012). Amphisbaenians, which likely share traits with other burrowing reptiles with absent or reduced limbs, prefer sandy, loose soils and avoid compact, heavy soils – presumably due to the reduced energetic costs of moving underground in sandy soils (Martín et al. 2013). Interestingly, however, we found no relationships between sandy textured soils and either Lerista species, which are semifossorial and known for their ability to ‘swim’ through sand (Morinaga and Bergmann 2020).

The importance of scale

Scale is an important factor in ecological studies because species respond to different drivers at different spatial scales (Cushman and McGarigal 2004). In the context of presence–background modelling, the area within which background locations are sampled (landscape extent) is an important consideration because it can influence model outputs (Elith et al. 2011; Merow et al. 2013). We chose the state of Victoria in south-eastern Australia as the landscape extent for all modelled species because it has a substantial gradient in both biophysical characteristics and landscape structure – the two groups of predictor variables of primary interest. Although using a common landscape extent to model the distributions of multiple species is a common approach (Guerin and Lowe 2013; Liu et al. 2013; Hageer et al. 2017; Swan et al. 2021), the importance of particular predictor variables is likely to be influenced by the range size of each species. For example, in our study, the relative influence of climate compared with landscape structure may be reduced for species whose range is larger than the landscape extent. Results should be interpreted with this in mind because the geographic range size of the species we modelled is variable.

Management implications

Although native vegetation is generally considered important for reptiles, most studies that model species distributions at large spatial scales focus on species responses to climate variables. Nevertheless, landscape structure is expected to be important to reptiles because the composition and configuration of habitat influence access to resources including shelter, food and mates (Fischer and Lindenmayer 2007). Our findings show that the extent of native vegetation had an important influence on the distribution of many terrestrial reptile species in our data set, including endangered species such as D. impar, E. tympanum and V. varius. Therefore, management that focuses on minimising the removal of native vegetation and increasing revegetation is fundamental to maintaining species persistence. Our study also showed that the configuration of native vegetation was important for some reptiles, including common species such as L. guichenoti and P. textilis. This suggests there is potential to use landscape structure metrics such as configuration to help identify patches of native vegetation where habitat restoration and revegetation are likely to have the greatest conservation outcomes (Jellinek 2017). Strategic revegetation may be especially pertinent in regions such as Victoria, where increases in human population and urban development over the next three decades (Commissioner for Environmental Sustainability Victoria 2018a) are likely to have further negative consequences for biodiversity.