Local habitat scale suitability
The model fit to the local study areas had an average 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 probability of P. muralis presence; bare ground, residential garden, dense scrub, scattered scrub, rail track, road, introduced shrub, dry dwarf shrub, hard cliff, and quarry (Fig 4). Spring radiance had the second highest percent contribution to the model (15%), where the amount of spring solar insolation had a positive influence on probability of presence (Fig 4). Probability of occurrence was also greater closer to buildings, railtrack, and roads. The response to NDVI is one of increasing probability of presence with an increase in vegetation from bare ground, followed by a rapid negative response past NDVI = 0.25. Maps indicating configuration of suitable habitat within local landscapes are presented in Figure 5 and Appendix 2.
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 5 and Appendix 2. Growth curves for the 10 study populations are also presented in supplementary information (Appendix 3). Growth rates ranged from 0.07 (Shoreham) to 0.15 (Eastbourne). Following simple stepwise linear regression analysis (Table 2), growth rate (r ) was positively related to the NLSI (F(1, 9) = 8.39, p = 0.02,R2 = 51.13), and negatively related to time since introduction (F (1, 9) = 5.80, p = 0.04, R2 = 42.22) (Fig 6 a, b). Branksome and Canford – two populations on the Bournemouth coast – had the highest carrying capacity (10443 and 10315 individuals, respectively). Eastbourne had the lowest carrying capacity (1447). A positive relationship between habitat quality and carrying capacity (F(1,9) = 6.22, p = 0.03,R2 = 43.74) was the only relationship observed between this growth parameter and the explanatory variables (Table 2). Annual dispersal distance was best explained by combined increases in NLSI and habitat quality (F (2,9) = 29.65,p <0.001, R2 = 89.44) (Fig 6 C), although habitat quality was not a significant predictor of annual dispersal distance on its own (F (1, 9) = 1.21,p = 0.34, R2 = 13.14). Greatest annual dispersal was predicted for the Eastbourne population (16 m), whilst the Shoreham, Wembdon, and Newton Ferrers populations had similar low dispersal of ~4 m per year. Connectance between suitable habitat patches had no relationship with any of the dependent variables.
Discussion
The predicted suitable climate for P. muralis in the UK is contiguous along the southeast coast, the entire south coast through to the south coast of Wales, extending northwards to a latitude of ~52°N - a latitudinal range likely to reflect climatic conditions most akin to those found in the species’ native origins. This northern limit to suitable conditions is in keeping with climate matching being an important limiting factor in determining establishment success and range expansion of introduced species, particularly significant for reptiles (Bomford, Kraus, Barry, & Lawrence, 2009; Pysek et al., 2010; Mahoney et al., 2015). Our climate suitability model output is similar to a previous SDM for P. muralis , which also highlighted favourable 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 model shows that number of frost days and amount of annual sunshine were the most informative variables in predicting probability of occurrence. 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, Townsley, & Bausch, 1990). Measurements of critical thermal minimum temperature in an introduced population of P. 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, Haro, Callejas, Labib, & Pauly, 2018). Interestingly, our model accurately predicted the Greater London Urban Area as having relatively high habitat suitability, likely arising from matching to thermal characteristics associated with the “urban heat island” (UHI) effect (Trajer, Mlinarik, Juhasz, & Bede-Fazekas, 2014; Villalobos-Jimenez & Hassall, 2017). There are historic records of small, P. muralis populations persisting in this area (Langton, Atkins, & Herbert, 2011; Langham, 2019), and since we did not include these records in the input for the model (due to no recent confirmed sightings and no accurate location data), the predicted suitability in this area gives credence to the validity of the model and the theory of UHI in built environments facilitating overwintering for the species. Dependence on human structures to survive winter temperatures in northern extremes has been suspected for introduced populations of Mediterranean gecko (Hemidactylus tursicus ) (Locey & Stone, 2006). Microclimatic conditions close to human habitations may have also facilitated establishment of Argentine ant (Linepithema humilein ) in areas with otherwise unsuitable climate (Roura-Pascual et al., 2011). Such environments may, however, also act to shield populations from selective pressures that might lead to adaptive physiological responses that could facilitate more rapid diffusion and expansion across wider areas (Hulbert, Hall, Mitchell, & Warner, 2020).
Our fine scale modelling of probability of occurrence 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, Strugariu, Sahlean, & Zamfirescu, 2009). It is possible that although we took great effort to assign habitat type in as much detail as practical, generalisations made during the construction of the habitat classification layer could possibly lead to overestimation of the extent of suitable 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 is well documented as providing important habitat for P. muralis , facilitating both natural dispersal and accidental human movement of animals (Covaciu - Markov, Bogdan, & Ferenti, 2006; Kühnis & Schmocker, 2008; Strugariu, Gherghel, & Zamfirescu, 2008; Gherghel et al., 2009). 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 to 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 1) 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 centred 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 favourable vegetated habitat (G. P. Brown, Phillips, Webb, & Shine, 2006). The presence of other contiguous, 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).
Our models concur, however, that natural dispersal of P. muralisfrom points of introductions in the UK is likely to be slow (Foster, 2015), with annual population range expansion of between 5 -16 meters. Spread distances were particularly small for populations in areas of relatively contiguous suitable habitat which allows for radial dispersal into suitable neighbouring habitat with limited search effort (i.e., rural villages with interconnected gardens, quarries) (Bonte et al., 2012; Baguette, Blanchet, Legrand, Stevens, & Turlure, 2013). 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, Kareiva, & Neubert, 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 of spread, indicating that localised habitat fragments are acting as stepping stones to dispersal (Alharbi & Petrovskii, 2019). Similar effects of landscape heterogeneity on range expansion of invasive species have been observed in introduced populations of whistling frog (Eleutherodactylus johnstonei ) (Ernst, Massemin, & Kowarik, 2011), Eurasian collared dove (Streptopelia decaocto ) (Ingenloff et al., 2017), and invasive weeds (Bergelson, Newman, & Floresroux, 1993).