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).