Climate and habitat configuration limit range expansion and patterns of dispersal in a non‐native lizard

Abstract Invasive species are one of the main causes of biodiversity loss worldwide. As introduced, populations increase in abundance and geographical range, so does the potential for negative impacts on native communities. As such, there is a need to better understand the processes driving range expansion as species become established in recipient landscapes. Through an investigation into capacity for population growth and range expansion of introduced populations of a non‐native lizard (Podarcis muralis), we aimed to demonstrate how multi‐scale factors influence spatial spread, population growth, and invasion potential in introduced species. We collated location records of P. muralis presence in England, UK through data collected from field surveys and a citizen science campaign. We used these data as input for presence‐background models to predict areas of climate suitability at a national‐scale (5 km resolution), and fine‐scale habitat suitability at the local scale (2 m resolution). We then integrated local models into an individual‐based modeling platform to simulate population dynamics and forecast range expansion for 10 populations in heterogeneous landscapes. National‐scale models indicated climate suitability has restricted the species to the southern parts of the UK, primarily by a latitudinal cline in overwintering conditions. Patterns of population growth and range expansion were related to differences in local landscape configuration and heterogeneity. Growth curves suggest populations could be in the early stages of exponential growth. However, annual rates of range expansion are predicted to be low (5–16 m). We conclude that extensive nationwide range expansion through secondary introduction is likely to be restricted by currently unsuitable climate beyond southern regions of the UK. However, exponential growth of local populations in habitats providing transport pathways is likely to increase opportunities for regional expansion. The broad habitat niche of P. muralis, coupled with configuration of habitat patches in the landscape, allows populations to increase locally with minimal dispersal.


| INTRODUC TI ON
The global rise in the number of species introduced to regions beyond their native range via human-mediated translocation shows no sign of reaching saturation point (Seebens et al. 2017). While many species fail to establish or have little negative effect following introduction, a subset of these do spread and can have significant impact on economies, human health, native biodiversity, and ecosystem services (Keller et al. 2011;Kolar & Lodge, 2001;Vila et al. 2010).
The severity of potential negative impacts (e.g., local or even global extinction of native species) is such that invasive non-native species (INNS) are justifiably regarded as one of the most significant threats to biodiversity worldwide (Genovesi, 2009;Simberloff et al. 2013).
For non-native species to become widespread and potentially damaging following introduction to new regions, introduced populations must negotiate the three stages of an introductionestablishment-invasion continuum (Blackburn et al. 2011).
Evaluation of the likelihood of a species to be transported, establish, and to spread, as well as the potential for having ecological, economical, and health impacts, forms the basis of "invasive" risk assessment for alien species (Bacher et al. 2018;Roy et al. 2019). Although it has been argued that the term "invasive" does not always necessarily equate with a species' negative impact (Ricciardi & Cohen, 2007), the potential for damaging effects inherently increases as introduced species increase in population size and spread across novel landscapes, thus affecting broader areas and more ecological communities (Crooks, 2005). As such, there is great interest in understanding patterns and rates of expansion of introduced species, and the environmental factors which limit their distributions (Gallien et al. 2010;Roy et al. 2019).
Following introduction and successful establishment beyond native ranges, species can further expand their range through local dispersal processes and/or by jump dispersal events that may be human-mediated (i.e., deliberate or accidental movement of individuals between habitats). Invading species typically exhibit several phases in the rate of spread. Firstly, there is an initial establishment phase where rate of spread is slow, secondly, an expansion phase typified by increasing rates of spread, and finally, a saturation phase when available space is occupied and expansion rates reach a plateau (Arim et al. 2006).
A suite of factors influences patterns and rates of range expansion during these phases: propagule size, dispersal mode, matching of physiological and ecological traits of invading species with environmental conditions at the receptor site, vital rates (births and deaths), species interactions, evolutionary processes, spatial heterogeneity, and temporal variability (Enders et al. 2020). Furthermore, our ability to assess and predict the temporal dynamics of invasions is often complicated by the phenomenon of lag phases; wherein, an introduced species remains at low population levels in the early stages of establishment for a protracted period of time before the sudden onset of rapid range expansion (see Crooks (2005) for review of the causes of temporal lags at all stages in the invasion process). Introduced populations of the northern Raccoon (Procyon lotor), for example, remained small for a number of years following introduction to Europe before a population explosion in the mid 1990s (Salgado, 2018). Similarly, landscape complexity can result in temporal and spatial patterns of invasion dynamics that deviate from classic theory of symmetrical, radial expansion from a central point (diffusion theory) (Kinezaki et al. 2010;Shigesada et al. 1995;Skellam, 1951). The effects of landscape heterogeneity on patterns and rates of expansion have been demonstrated in the quick colonization of areas of suitable habitat in the early stages of the American mink (Neovison vison) invasion, compared with uptake in areas of low habitat suitability in Scotland (Fraser et al. 2015), and the fluctuating rates in range expansion of Cane toad (Rhinella marina) in response to changing environmental conditions in newly invaded areas of Australia (Urban et al. 2008). Consideration of dispersal processes across heterogeneous landscapes is therefore central to predicting potential for range expansion during the invasion process Grayson & Johnson, 2018;Travis et al. 2011). The development of platforms for spatially explicit individual-based modeling (Bocedi, Zurell, et al., 2014;Samson et al., 2017) has enabled the nested interactions between dispersal, landscape properties, and population dynamics to be considered in predicting species distributions, increasing the ecological realism of range expansion models (Andrew & Ustin, 2010;Ferrari et al. 2014;Hunter-Ayad & Hassall, 2020;Mang et al. 2018).
In this study, we determine the potential for range expansion of the non-native common wall lizard (Podarcis muralis) in the UK.
Podarcis muralis has a long history of introductions beyond its native range which covers most of Western and Southern Europe (Gassert et al. 2013). Many of these introductions continue to extend the species' range throughout continental Europe (Oskyrko et al. 2020;Šandera, 2017;Santos et al. 2019;Wirga & Majtyka, 2015), but the species also has several populations established in the New World, both in the United States (Brown et al. 1995) and Canada (Allan et al. 2006). Introduced to Vancouver Island, British Columbia, in 1970, the species persisted in isolated populations until 2006, but has since spread with alarming speed due to jump dispersal (humanmediated) and natural radial dispersal of 40-70 m a year in urban areas (Engelstoft et al. 2020).
To date, there is no empirical evidence of negative ecological impacts of P. muralis introductions in the UK, and there is mixed social perception and opinion toward the species' presence (Williams et al. 2019). However, suspected declines in native lizards through interference and/or exploitation competition have been reported following introductions of P. muralis, to both Germany (Kühnis & Schmocker, 2008;Münch, 2001;Schulte, 2009) and the UK (Mole, 2010).
There have been multiple introduction events of P. muralis to the UK both as deliberate releases of captive animals and as cargo stowaways, with some extant populations having been established on the UK mainland as early as the 1970s (Michaelides et al. 2013). More recent introductions (1980s onward) have mostly arisen from movement of individuals from already established populations (secondary introduction) or captive-bred animals, rather than directly sourced from the native range . The UK populations represent the species at the northern extent of its range, with sites having markedly different climatic conditions compared to source regions of the native range. For example, air temperatures experienced by populations during the main activity season at sites in England are 5-10°C lower than their source regions in Tuscany and western France .
We investigated the potential for range expansion of P. muralis in the UK with models highlighting different (but complementary) parameters likely to influence spread at two spatial scales. Firstly, since long-distance jump dispersal via translocation is important in facilitating spread of this species, we aimed to predict the national extent of the area potentially available for further colonization by running a species distribution model (SDM) based on current climatic suitability at these northern extremes. As has been speculated elsewhere, the ability to survive cold winters is likely limiting to the spread of introduced Podarcis populations ( Burke et al. 2002). We, therefore, hypothesized that latitudinal clines in climate would restrict the area available for northward expansion of P. muralis via long-distance human-assisted translocation in the UK. Climate change scenarios may significantly affect the future extent of habitat suitability for P. muralis. However, in this study, we are primarily interested in trying to understand what is limiting the current distribution and not necessarily trying to predict where the species' will be in the future.
Second, to make predictions of population growth and dispersal patterns, as well as identify environmental features important to range expansion at a local level, we took a hybrid model approach combining SDMs, informed by variables characterizing 10 local landscapes (i.e., microclimate, proximity to geographic features, and habitat type), with a high resolution (15 × 15 m) spatiotemporal individual-based model (IBM) simulating local population and dispersal dynamics. We expected that landscape characteristics (i.e., configuration and connectivity of suitable habitat patches) would result in asymmetrical patterns of predicted dispersal within populations, which in turn would result in spatial and temporal variance in patterns of population growth and range expansion between populations. These analyses allow us to investigate the proximate and ultimate barriers to spread, as well as simulating the potential for invasion lag in each population.

| ME THODS
The locations of known established P. muralis populations were obtained from data collated on the Wall Lizard Project hosted on the Surrey Amphibian and Reptile Group website (Langham, 2019).
We determined the current geographic extent of as many of these populations as logistically possible (10 sites) using a combination of visual surveys, canvassing of the local public at sites of interest, and press releases in local and regional media encouraging members of the public to report their wall lizard sightings (see Appendix S1). Of the 30 extant populations recorded on the UK mainland, we visited 21 between three field seasons (April-September) of 2016, 2017, and 2018. We did not visit the remaining nine locations because the lizard populations were either known to be very small, had no recent confirmed sightings and no accurate location data, access was restricted, and/or site locations were otherwise logistically challenging (i.e., distance from other populations). Specific attention was given to assessing the extent of P. muralis presence along railway habitat at West Worthing, Sussex (50.818°N, 0.390°W) during a five-week period in June-July 2018-the railway acting as a linear transect, along which we could assess the utility of railway habitat as a corridor for dispersal ( Figure S1 Appendix S1). We also collected wall lizard location data from populations at two additional sites, Eastbourne (50.768°N, 0.291°E) and Kingswear (50.349°N, 3.568°W). The former was confirmed to be an established colony through a site visit, and the latter was reported to RW by a member of the public responding to a citizen science campaign in regional media. Sightings were confirmed for a further eight new locations as a result of the citizen science campaign ( Figure 1) (see Appendix S1 for detailed methods). Overall, we collated a total of 1,331 lizard sightings (76 from online portal, 52 from postcard returns, and 1,203 from visual surveys) across 25 sites ( Figure 1).

| Modeling climate suitability across the UK
We used our UK P. muralis presence data and records from the native range (GBIF.org, 2020) to develop relative habitat suitability maps at the UK national extent using MaxEnt v3.3.3k software . We minimized spatial autocorrelation (due to possible clustered sampling effort) between data points by thinning all data to a minimal distance of 10 km between points using the R package SpThin (Aiello-Lammens et al. 2015), resulting in 25 UK records and 1,542 from the native range to be used in the models. Podarcis muralis has demonstrated rapid adaptive responses following introduction to cool climates, with ability to evolve broader thermal tolerance at dispersal fronts (Litmer & Murray, 2019), and prolonged embryo retention and faster embryonic growth at low temperaturescompared with ancestral states . To best reflect the current climatic tolerance of the species, we therefore compared models that used P. muralis presence records from the introduced UK range only (acknowledging the possibility for regional overfitting with such data), with models constructed using only native occurrence data, and the UK and native occurrence data combined.

| Environmental variables
We selected 10 climatic variables that have most relevance to wall lizard biology and therefore likely to influence distribution (Wirga & Majtyka, 2015) (Table 1), obtained from WorldClim-Global Climate Data (Fick & Hijmans, 2017) and E-OBS datasets from the EU-FP6 project ENSEMBLES (Cornes et al. 2018). Due to the different resolution in data from these sources, we up-scaled E-OBS climatic variables by using the bilinear interpolation to a spatial resolution of 0.0083°. These 10 variables were refined from an initial input of 13 climate variables through an iterative process of removal/retention to limit covariate correlation (Spearman's rank correlation; correlated if r s ≥ 0.6) and maximize model performance (Glover-Kapfer, 2015).
We kept parameter settings in MaxEnt the same for modeling at the national and local levels (see Appendix S1 for detailed method).

| Modeling local habitat suitability
A total of 1,083 presence records (all direct observations during visual surveys), across 10 study locations representing the range of habitats used by P. muralis (urban, suburban, and 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 2 m resolution were used for the MaxEnt input and are summarized in Sampling bias is often inherent in occurrence data (Merow et al. 2013) and there are several aspects of our approach to modeling habitat suitability at the local level which cumulatively addresses any issue of a bias effect in our presence data: (a) The points of local introduction are explicitly known for the various populations studied, and thus, clustered occurrence data are inevitable (and expected). (b) Our combined survey methods (citizen science, F I G U R E 1 Locations of UK breeding populations of Podarcis muralis (black) from which presence data informed models of climatic suitability (n = 25) and locations of confirmed sightings (red) arising from a citizen science campaign 2017-2019 (n = 8) canvasing, and visual surveys) ensured we could confidently survey from the core of the population out toward the extremities, giving confidence to our assessment of habitat use and local range extent. (c) Our models use as input occurrence data and manually generated background points from across all populations/study sites. The model was then projected to each of the study sites separately. Local projections, therefore, represented the suite of habitat types used (ranging from rural to highly urbanized) across all sitesoffsetting any local clustering around features (i.e., buildings and roads).

| Modeling local range expansion (IBM)
We used spatially explicit individual-based models to identify the functional responses of P. muralis to the habitat structure provided by MaxEnt models. These models were constructed using the platform RangeShifter v1.1 (Bocedi, Palmer, et al., 2014)

| Parameterization
Parameters of wall lizard demographics and behavioral attributes were based on empirical data in the published literature.
Parameterization was further calibrated 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 modeled as closely as possible the currently observed spatial extent of each study population (Fraser et al. 2015). Where published empirical data were not available, reasonable judgments and/or simplifying assumptions were made that were biologically realistic and justifiably reflect the functional biology of P. muralis (Table S1 in Appendix S1).

| Initialization
Simulations were initialized using known founder size where documented (Langham, 2019;Michaelides et al. 2015). 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 (probability every population (

| ANALYS IS
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 (R Core Team, 2017) using the package growthcurver (Wagner, 2016 McGarigal et al. 2002). An illustration of workflow in layer preparation from SDM to IBM to FRAGSTATS input is provided in Figure S2 Appendix S1.
We ran (simple stepwise) linear regression models with two FRAGSTAT metrics describing heterogeneity of suitable habitat patches within the landscape (Normalized Landscape Shape Index-a measure of patch aggregation where NLSI = 0 when the landscape consists of a maximally compact patch of the corresponding type and increases (to a maximum of 1) as the patch type becomes increasingly disaggregated; and Connectance-a measure of functional joinings of patches reported as a percentage of the maximum possible connectance given the number of patches) and average habitat quality as explanatory variables, and the growth rate parameters (K, r) and annual range expansion as response variables. We set the threshold distance within which patches are deemed "connected" to an arbitrary 100 m. Values for NLSI and Connectance for each site are provided in Table S2 Appendix S1. Spring radiance made the highest relative contribution to both models and minimum winter temperature the least contribution.

| National-scale climatic suitability
The models indicate suitable climatic conditions for P. muralis in  Figure S3 Appendix S2.

| Local scale habitat suitability
The model fit to the local study areas had an average test AUC of 0.88 (SD = 0.01) over the 10 areas and 50 replicated runs. The most important variable to the model was "habitat type," which made the highest relative contribution to the model (66%). Ten habitat classes out of 44 stood out as being influential to increased suitability for P. muralis; bare ground (1), residential garden (2), dense scrub (5), scattered scrub (6), rail track (17), road (18), introduced shrub (22), dry dwarf shrub (25), hard cliff (28), and quarry (37) (Figure 3).
Spring insolation had the second-highest percent contribution to the model (16%), where the amount of spring solar insolation had a positive influence on relative suitability (Figure 3). Relative suitability was also greater closer to buildings, railtrack, and roads. The response to NDVI is one of increasing suitability with an increase in vegetation from bare ground, followed by a rapid negative response past NDVI = 0.5. Maps indicating configuration of suitable habitat within local landscapes are presented in Figure 4 and Figure S4 Appendix S2.

| Individual-based models results
Patterns of range expansion from time since introduction to 2040, as determined by population dynamics and local landscape character, are presented in Figure 4 and Figure (Table S3). A positive relationship between habitat quality and carrying capacity (F (1,9) = 6.22, p = 0.03, R 2 = 43.74) was the only relationship observed between this growth parameter and the explanatory variables (Table 2). Annual range expansion was best explained by combined increases in NLSI and habitat quality (F (2, 9) = 29.65, p < 0.001, R 2 = 89.44) (Figure 5c), although habitat quality was not a significant predictor of annual dispersal distance on its own (F (1,9) = 1.21, p = 0.34, R 2 = 13.14). Greatest annual range expansion was predicted for the Eastbourne population

| D ISCUSS I ON
The predicted suitable climate for P. muralis in the UK is continuous along the southeast coast, the entire south coast through to the south coast of Wales, extending northward to a latitude of ~52°N-a latitudinal range most likely to reflect maritime climatic conditions found in the species' native origins. This northern limit  (Bomford et al. 2009;Mahoney et al. 2015;Pysek et al. 2010). Our climate suitability model output is similar to a previous SDM for P. muralis, which also highlighted suitable conditions in the UK up to ~53°N, (Wirga & Majtyka, 2015), despite a differing suite of climatic variables and species presence data informing models.
Our models show that amount of spring radiance, maximum winter temperature, and number of winter ice days were the most informative variables in predicting habitat suitability (where both native and allochthonous records were considered). The hibernation period is short in P. muralis, and individuals are often active in mid-winter during sunny mild spells, even in the northern extremes of their range, making them vulnerable to sudden or prolonged freezing (Claussen et al. 1990). Measurements of critical thermal minimum temperature in an introduced population of Podarcis sicula have been shown to be above temperatures likely experienced by some non-native populations in winter, suggesting individuals may need to find urban thermal retreats to survive winter conditions, or hibernate at a depth below soil freezing to survive (Burke et al. 2002;Liwanag et al. 2018). Interestingly, our UK only model accurately predicted the Greater London Urban Area as having Our fine-scale modeling of habitat suitability provides a detailed insight into local landscape structure and spatial pattern of available suitable habitat. The contribution of habitat classification and spring solar insolation to the model, and particularly the unimodal response observed toward vegetation cover (NDVI), is indicative of the species' affinities to disturbed habitats that provide resource for refugia (thermal and safety), egg deposition sites, and basking sites necessary for heliothermic temperature regulation (Bertram, 2004;Gherghel et al. 2009). It is possible that although we took great effort to assign habitat type in as much detail as practical, generalizations made during the construction of the habitat classification layer could possibly lead to overestimation of the extent of suitable F I G U R E 5 Relationship between growth rate (r) and (a) aggregation of suitable habitat (NLSI), time since introduction (b), relationship between NLSI and annual range expansion, and (c) in non-native population of Podarcis muralis in the UK habitat (e.g., not all habitat classed as residential garden would, in reality, be suitable to P. muralis). However, the combined effect of the NDVI variable would go some way to enhance fine-scale delineation between suitable and unsuitable habitat type.
The relative importance of railway line and introduced shrub habitat in the model can be explained by the number of presence records associated with those habitats in relation to the relative scarcity of those habitats in the landscape. Habitat associated with railway lines provides important habitat for P. muralis, facilitating both natural dispersal and accidental human movement of animals (CovaciuMarkov et al. 2006;Gherghel et al. 2009;Kühnis & Schmocker, 2008;Strugariu et al. 2008). Dispersal of the introduced P. muralis population in Ohio, Cincinnati, has been reported to be more rapid along the continuous hospitable terrain of rail embankments compared with the relatively slow spread through highly fragmented residential and commercial areas (Hedeen & Hedeen, 1999). Although our simulations of the West Worthing (trackside) population (see Appendix S1) did have relatively higher dispersal distance than most other populations, the pattern of spread did not indicate extensive natural dispersal along the railway, despite the core population being centered on, and around, disused sidings, and associated habitat. Instead, the simulated dispersal pattern is one of predominantly radial diffusion out into adjacent residential and commercial areas, where, although highly fragmented, the habitat was of suitable quality to facilitate this pattern of spread. Linear corridors may, therefore, only become important to natural dispersal when adjacent habitat is of low quality, or is less preferred, as is the case of invasive cane toads (Rhinella marina) selecting to use open roads for dispersal through less favorable vegetated habitat (Brown et al. 2006). The presence of other continuous, linear habitat features in our landscape models also increased rates of annual range expansion (e.g., vegetated cliff faces at Branksome and Canford; sea front garden along the promenade at Eastbourne), but this is likely a result of there being restrictions to radial dispersal as suitable habitat is bordered by inferior inland habitat and the shore line. Our findings are congruent with the theory that corridors may be most effective when they actively influence, direct, and channel dispersal rather than simply provide additional suitable habitat (Andrew & Ustin, 2010).
Growth curves derived from our predictive models suggest all the populations studied may be in the early stages of exponential growth, and have demonstrated (or are demonstrating) a lag before the onset of appreciable population growth that is often associated with such a growth trajectory (Sakai et al. 2001). The negative correlation we found between intrinsic growth rate and time since introduction, is to be expected as a function of logistic growth, where the longer-established populations approach local carrying capacity and density dependence constrains growth (Sibly & Hone, 2002).
Increasing densities in urban areas which present transport pathways may lead to an increase long-distance dispersal events and secondary introductions. A similar scenario has been described in the spread dynamics of invasive Pallas's squirrel (Callosciurus erythraeus), where a constant increase in the appearance of new populations occurred after a two-decade lag, and could be explained by increased vector activity (intentional translocations) as the population size at the initial introduction foci increased-causing a potentially exponential increase of translocation events when translocated populations start acting as a source themselves (Guichón et al., 2015).
Our models concur, however, that natural dispersal of P. muralis from points of introductions in the UK is likely to be slow (Foster, 2015), with annual population range expansion of between 5 and16 meters. Spread distances were particularly small for populations in areas of relatively continuous suitable habitat which allows for radial dispersal into suitable neighboring habitat with limited search effort (i.e., rural villages with interconnected gardens, quarries) (Baguette et al. 2013;Bonte et al. 2012). In such instances, it would appear that populations with limited opportunities/need for long-distance dispersal are increasing their numbers locally, but will be limited for establishing a population over a large area (Lustig et al. 2017). Increasing disaggregation of suitable habitat had a joint positive influence on dispersal rate and growth rate in our models.
We found this to be most apparent for the urban population of West Worthing, highlighting how the species' ability to exploit areas of human disturbance may facilitate overall invasion success (Marvier et al. 2004). Increasing abundance of discrete local patches of suitable habitat may provide opportunity for individuals to disperse more widely in the landscape, thus releasing density-dependent constraints on population growth that would be in effect when suitable habitat is more aggregated and compact. This pattern is in line with the theories of a percolation threshold, where invasive spread may occur most rapidly and extensively above a threshold level of disturbance (i.e., amount of habitat fragmentation) (With, 2002).
In addition, we found functional connectedness of suitable habitat patches had no relation to any of the growth parameters or rate  (Ingenloff et al. 2017), and invasive weeds (Bergelson et al. 1993).

| CON CLUS IONS
Extensive nationwide range expansion through secondary introduction is likely to be restricted by currently unsuitable climate beyond southern regions of the UK. However, exponential growth of local populations in habitats providing transport pathways (e.g., movement of aggregates, timber, plants, and general public) is likely to increase opportunities for long-distance dispersal and regional expansion.
Despite the fundamental physiological importance of sun exposure to diurnal reptiles, to the best of our knowledge, our models are the first to incorporate estimates of solar insolation into models of habitat suitability at this fine-scale (but see Garcia-Porta et al. (2019) for consideration of solar radiation as a variable in driving macroecological and macroevolutionary patterns in lizards). Our models demonstrate the inclusion of the variable at this scale, and indeed our entire approach to developing a fine-scale SDM could be very useful in other applications relating to ectotherm ecology (e.g., in developing habitat suitability indices, directing habitat management, and guiding survey effort for rare/cryptic species). Furthermore, while the use of SDMs and IBMs have become widely used to further understanding of mechanisms driving invasion dynamics (Fraser et al. 2015;Kadoya & Washitani, 2010;Suzuki-Ohno et al. 2017), the benefits of incorporating spatially explicit individual-based models into management plans for the control of invasive species have only recently been recognized (Day et al. 2018). In this regard, our models provide a best estimate for future expansion of P. muralis (under current climate conditions) at both the UK national and local scale, providing essential information (i.e., dispersal patterns, key habitat, current and projected population sizes) on which management decisions could be made. It is important to note that forecasts of climate warming may increase suitable area for range expansion and local carrying capacity. As our models show UK populations may be approaching an end to an inherent lag phase, there is an argument to suggest that a timely precautionary intervention may be justified to halt some populations before an abrupt end to the lag phase occurs.

ACK N OWLED G M ENTS
The authors would like to thank Network Rail, Natural England, and all those who granted site access for their invaluable help in facilitating this research. Thanks to all the media outlets who promoted our citizen science campaign, and to those who reported their lizard sightings. Thanks to Ellie Field, Joel Dixon, Hannah Grice, Miguel Martinez, Liz Cooper, and Jack Carney for assistance in the field.
Steve Langham provided additional site details for which we are very grateful. We thank Dr Steve Palmer and Dr Greta Bocedi for their guidance and technical support in using the RangeShifter software. This work was supported by a Natural Environment Research Council Studentship to RJW (NE/L002574/1).

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All spatial data including presence and background samples with