Modelling local habitat suitability
A total of 1083 presence records (all direct observations during visual surveys), across 10 study locations representing the range of habitats used by P. muralis (urban, suburban, rural), were used in producing relative habitat suitability maps and predictive models of local range expansion. These study sites encompassed heterogeneous land cover that helped in identifying variables affecting local habitat suitability and features acting as important corridors for range expansion. Data for six environmental variables at 2m resolution were used for the MaxEnt input and are summarised in Table 1. All variables were calculated and prepared in ArcGIS® (Esri 2017). We used the Phase One Habitat Survey Toolkit (Centre for Ecology Environment and Conservation, 2018) to create fine scale habitat type (categorical) data layers.
Modelling local range expansion (IBM )
Habitat suitability maps from our local scale MaxEnt models were prepared as habitat quality landscape layers by linear transformation of the MaxEnt logistic values (estimates between 0 and 1 of probability of presence) above the maximum test sensitivity plus specificity logistic threshold. This is the threshold at which the MaxEnt models maximize their discrimination of presences from background data (Jimenez-Valverde & Lobo, 2007; Glover-Kapfer, 2015). The resulting habitat quality landscape (scaled 0-100, and where cell values scale with cell carrying capacity in RangeShifter), provided the patch input for RangeShifter v1.1 (Bocedi, Palmer, et al., 2014), in addition to a cost layer to movement created by reclassifying (inverting) the habitat quality landscape layer. All inputs were resampled using bilinear interpolation to 15m x 15m cell size to reduce demands on computational memory whilst retaining biological relevance to wall lizard movement capabilities. A single cell in each landscape was identified as the initial species distribution (i.e., point of introduction for each population respectively) based on knowledge of the precise location of introduction when known, or by using the centre point of the current extent of sighting records for the population.
Parameterisation
Parameters of wall lizard demographics and behavioural attributes were based on empirical data in the published literature. Where published empirical data were not available, reasonable judgements and/or simplifying assumptions were made. The final parameter values used were biologically realistic and justifiably reflect the functional biology ofP. muralis (Table S2 in Appendix 1). Parameterisation was further refined through an iterative process, where simulations were repeated across all study sites with fine parameter adjustments within biologically meaningful limits until a single set of parameters was found that modelled as closely as possible the currently observed spatial extent of each study population (Fraser et al., 2015).
Initialisation
Simulations were initialised using known founder size where documented (Michaelides et al., 2015; Langham, 2019). Where founder size was unknown, we used a minimal founder size that resulted in reasonable simulation outputs as per the iterative process mentioned above. We assumed adult age class for all founders. Local extinction probability was set at a constant of 0.003 across populations. Simulations (50 replicates) of population range expansion for the 10 study populations were then run for the period of time since introduction (which varies among sites) up to the year 2040.
Analysis
We investigated how landscape characteristics might influence population size, rate of population growth and range expansion, by first obtaining standard population growth metrics: carrying capacity (K ), and intrinsic rate of increase (r ), from linear growth curves applied to mean yearly population size data taken across all simulation iterations in R Studio (R Core Team, 2017) using the package Growthcurver (Sprouffske, 2018). We then created binary habitat suitability layers from our MaxEnt outputs for a radius of 200m around introduction points which served as inputs for the programme FRAGSTATS v4 (McGarigal, Cushman, & Ene, 2002). We ran linear regression models with two FRAGSTAT metrics describing heterogeneity of suitable habitat patches within the landscape (Normalised Landscape Shape Index – a measure of patch aggregation; and Connectance – a measure of functional joinings of patches) and average habitat quality as explanatory variables, and the growth rate parameters (k, r ) and annual dispersal distance as response variables. We set the threshold distance within which patches are deemed ”connected” to an arbitrary 100m.