Introduction

Plant invasions represent a threat to global biodiversity and affect almost all types of environment (Pimentel et al. 2005). Rapid urbanization may represent an important pathway of biological invasions, particularly when it occurs alongside other landscape-level transitions (Duguay et al. 2007; Giorgis et al. 2011). Montane environments that are accessible to large urban populations in adjacent lowlands may be vulnerable to expansion of settlements due to rising second home ownership, and associated increases in the introduction of non-native plants in gardens and hedgerows (Pauchard et al. 2009). The escape of these horticultural introductions into surrounding semi-natural habitats may then be promoted by simultaneous changes in land use and land management, such as declining rural populations and reductions in the density of domestic grazing animals (Pauchard et al. 2009). The complexity of these changes makes it very difficult to predict the rate and spatial pattern of species’ invasions, especially when long-term demographic data are lacking. In this paper we outline a methodology designed to assess the vulnerability of different habitats and inform control strategies in a setting where long-term data are unavailable and resources for monitoring are limited, and then apply this approach to predicting the range expansion of the invasive shrub Pyracantha angustifolia (hereon Pyracantha) in a rapidly changing montane landscape in northwestern Argentina.

Our approach relies on the use of spatial population models to provide valuable information about the relative speeds of invasion in different conditions. This is particularly useful in systems undergoing rapid invasions, especially with limited data on species distribution and the dynamics of spread (Buckley et al. 2003a; Jongejans et al. 2008; Harris et al. 2011; Travis et al. 2011). In general, modelling these spatial dynamics is challenging because range expansion is dependent on multiple processes including pollination, dispersal, germination, and seedling establishment, and these along with key demographic parameters related to fecundity, survival and dispersal can all be influenced by the local conditions (Harris et al. 2011; Damschen et al. 2014).

Quantifying how the demographic properties of the species (its age-specific fecundity, survivorship and dispersal ability) are influenced by, and interact with, the environment is central to understanding how an invasive species spreads in different conditions (Buckley et al. 2003a; Dawson et al. 2009; Harris et al. 2011). For example, habitat type influences the fecundity of the shrub Rhododendrum ponticum, leading to differences in range expansion and invasion density (Harris et al. 2011). Similarly, Hypericium perfiratum showed low amounts of fruit production in shaded sites compared to open habitats (Buckley et al. 2003a). Consequently, populations in shaded areas respond less to grazing, which is currently used as a control agent (Buckley et al. 2003b). These examples illustrate that understanding how demography is influenced by the local conditions may have important implications for management, as well as improving predictions of range expansion. However, there are few studies that link observations of local influences on demography with population-level impacts on plant invasive dynamics, and those that do, provide insights into tackling the spread of their study species (Buckley et al. 2003b; Harris et al. 2011).

Changes in propagule pressure due to changes in human activities have led to an increase of invasive woody species in mountainous areas across central and northern Argentina (Giorgis et al. 2011), and some of these species are having increasing impacts on native ecosystems. Many of these invasive species are bird-dispersed and their success is likely due to the attractiveness of the fleshy fruits to frugivorous birds facilitating dispersal (Gosper et al. 2005), which is enhanced by the availability of the fruit of some non-native species during seasons when the fruits of native species are scarce (Rouges and Blake 2001; Blendinger et al. 2012). Pyracantha is one such species, and it should be a major cause for concern, because it is currently in the early stages of invasion in some mountainous areas of northwestern Argentina. Studies on Pyracantha to date have focussed on its local impacts, ecology, and interactions with other species (Tecco et al. 2006; Giantomasi et al. 2008; Vergara-Tabares and Rojas 2016). There have been no studies on its demography and no predictions made for how rapidly it is likely to spread under different conditions.

In this study our aim was to gain new insight into how propagule pressure of Pyracantha is likely to depend upon both local habitat type and grazing pressure. Ultimately it will be important to consider impacts of habitat and grazing on multiple other potentially important processes including dispersal and establishment. However, our initial focus here is on determining how age-specific fecundity regimes are influenced by both habitat and grazing pressure and then using that new information to build models to assess how much this heterogeneity is likely to influence rates of range expansion.

Thus, our aim is to address two main questions (1) how do different habitats and grazing scenarios influence fruit production for an invasive population of Pyracantha and, (2) how do changes in Pyracantha fecundity influence population density and rates of expansion at the invasion front under different scenarios?

Materials and methods

Study site and species

Pyracantha angustifolia (Rosaceae) is a woody thorny shrub, native to China. In its native range it forms thickets on slopes and along roadsides at altitudes of 1600–3000 m (Cuizhi and Spongberg 2003). It has become widely spread in shrublands, grasslands, and in mountain areas across central (Tecco et al. 2006) and northwestern Argentina (Aráoz and Grau 2010). It produces fleshy fruits that reach maturity between May and October that are consumed and dispersed by birds (Vergara-Tabares and Rojas 2016).

Demographic data were collected at Tafi Del Valle, Tucuman, Argentina (26° 55′ S, 65° 41′ W) in May 2019 when Pyracantha was fruiting (Tecco et al., 2006). Tafi Del Valle is an intramountain valley, part of the Yungas ecoregion, where grasslands and shrublands intercalate with forest patches rich in Alnus acuminata (Betulaceae) (Aráoz and Grau 2010). Forest is located on south facing slopes and on moist sites at the bottom of the valley (Aráoz and Grau 2010). Patches of the tree Polylepis australis (Rosaceae) develop on upper slopes. Shrublands are commonly dominated by Acanthostyles buniifoliu (Acanthostyles). Species from the genera Chuquiraga and Adesmia (many of which are thorny) form dense shrublands that are vulnerable to invasion by Pyracantha. Native thorny species from the genera Berberis are the most similar shrub to Pyracantha, because they possess a similar fruit morphology and are also dispersed by birds, while non-native Cotoneaster sp. are potentially invasive in the valley in addition to Pyracantha. The grasslands are dominated by species of Festuca and Stipa. The altitude of the valley ranges from 1800 to 4000 masl (Busnelli et al. 2006) and our plots were located between 1840 and 2420 masl. The climate is monsoonal with mean annual rainfall of 400 mm and a mean annual temperature of 13.5 °C (Sampietro-Vattuone and Peña-Monné 2019). The winter is dry and relatively cold, with some freezing temperatures and occasional snowfall, the average winter temperature is 9 °C (Zinck and Sayago 1999). The summers are hot and humid; 82% of the annual rainfall occurs between November and March (Busnelli et al. 2006), the annual temperature of summer is 18 °C (Zinck and Sayago 1999).

Previously the valley was used primarily for cattle holdings but over the last 30 years the number of holdings has declined (Rainer and Gaitán 2014), which is contributing to an expansion of montane forest cover (Aráoz and Grau 2010). The area is also grazed by horses and sheep. The area has undergone rapid changes and urban expansion caused by an increase in second-home ownership (Rainer and Gaitán 2014), which has included the planting of Pyracantha as hedgerows and ornamental plants in private gardens, as in other intermountain valleys (Gurvich et al. 2005). The natural vegetation is strongly disturbed by grazing and consists mainly of open habitats dominated by grasses and shrubs (Busnelli et al. 2006).

We collected data from 10 sites across the valley (Fig. S1), chosen because they cover a range of environmental variables (habitat types, elevation, approximate distances to water sources and urban areas). We could not accurately measure these based on satellite images or in the field and therefore these data are not used in subsequent analysis. Where established reproductively mature Pyracantha occurred we chose sites where we observed that the population dynamics of Pyracantha were not affected by direct human activity (such as planting, pruning, removing etc.). Occupancy by Pyracantha was estimated based on preliminary observations of aerial images (Google Earth Pro) and field surveys. Preliminary field observations led us to define three common habitats where Pyracantha has invaded. These were categorized as grassland (agricultural pastures that are actively used for grazing cattle and horses), shrubland (characterized by the presence of native shrubs) and rocky (seasonally dry stream beds and boulder fields characterized by the presence of large boulders and low density of established vegetation).

Field methods

In each of the ten sites we randomly selected the middle point for three circular plots (20 m diameter) separated by a minimum distance of 20 m (30 plots representing a total area of 9426 m2). We recorded the coordinates and elevation at the centre of each plot with a gps device (Garmin eTrex Touch 35). We also recorded the habitat type and the number of saplings of Pyracantha (greater than 10 cm tall and non-reproductive) and reproductive adults in a plot. For each reproductive adult we measured the basal stem diameter of the main stem. We recorded whether the plant showed visible signs of damage by grazers, however it was not possible to measure grazing intensity. Finally, we counted all the fruits on segments comprising one eighth of the projected canopy area, based on a random compass bearing to select the segment for sampling. Fruit production per plant was determined by multiplying this value by eight for a total of 309 reproductive adults of known size.

We counted the number of seeds in 30 fruits from 10 different shrubs and consistently found five seeds per fruit. Raw fruit counts were used in the data analysis. To estimate the total fecundity for population models, the raw fruit counts were multiplied by eight (proportion of shrub counted) and then five (number of seeds per fruit).

Plant age estimation

To determine estimates of plant age we harvested 32 shrubs at ground level and removed a cylindrical block of stem from the base of the main stem (basal stem diameter ranged from 1.05 to 9.7 cm). One cross-section on each sample was sanded to produce a smooth surface. Growth rings were counted based on observation of denser and darker coloured latewood bands, which are assumed to represent annual growth rings. This was validated based on knowledge gained from satellite images of the approximate date of introduction of Pyracantha at the site compared with the projected age of the oldest shrubs. Counts were made along the radii with the most distinct rings and free of abnormalities or damage. Each sample was counted by three observers under a dissecting microscope with a light (Kyowa optical model: SDZ-TR-PL) at ×15 magnification and the basal stem diameter was measured.

Data analysis

All analyses were conducted in R version 3.6.1 (R Core Team 2019) using linear models fitted with the lm function, generalized additive models (GAMs) fitted with the “mgcv” package (Wood 2011) and graphs plotted with the “ggplot2” package (Wickham 2016). For all models the assumptions of homogeneity of variance and normality of residuals were checked using diagnostic plots (Appendix 2).

Age—size relationship model

The relationship between basal stem diameter and age was tested using a linear model with number of growth rings modelled as a function of basal stem diameter. This model was used to predict the age of all 309 shrubs in the data set and these ages were used in the subsequent fecundity modelling.

Fecundity models

Generalised additive models (GAMs) were used to estimate trends in fecundity with age for different treatments (elevation, habitat type, grazing and plant density).

All GAMs included a smooth function of age, so could support a non-linear relationship between plant age and fecundity. The estimated degrees of freedom (edf) were used to determine whether the relationship was non-linear. The response (raw fruit counts) was fitted on the natural log scale to account for the sharp increase in fecundity with age. All models assumed a Gaussian error and used an identity link function. To account for non-independence of plots within study sites, site was fitted as a random effect.

We tested different combinations of the fixed effects of habitat type, grazing, elevation and/or density of shrubs. Models were built where smooth functions could vary between different combinations of grazing regime (grazed or ungrazed) and/or habitat type to test whether the model supported the inclusion of interactions between age and the treatments. The model fit was assessed using the Akaike information Criterion (AIC, Akaike 1987). The final presented model had the lowest AIC value and for this model an assessment of the contribution of each variable was made using an Analysis of Variance (ANOVA).

Modelling invasion

There are a range of different approaches that can be applied to model the dynamics of spread (reviewed by Jongejans et al. 2008). Two common approaches are individual based models (IBMs) and analytical approaches, which couple the integrodifference equation (IDE) for dispersal with matrix models (Neubert and Caswell 2000; Travis et al. 2011). Since these approaches both have strengths and weaknesses, the two methodologies are applied to provide alternative predictions that were then evaluated by cross-validation (Travis et al. 2011).

We run the simulations in a modified version of RangeShifter V1 (Bocedi et al. 2014) and R version 3.6.1 (R Core Team 2019) for the IBM and analytical model respectively. To determine the influence of habitat type and grazing on invasion speed we set up a fully factorial design where outputs were derived for each of the six possible combinations of treatments. Survival and dispersal were constant for all models and only fecundity parameters changed.

Parameters

The regression equations from best fitting models determined as above were used to predict the age-specific fecundity (fecundity vs age relationship) for each habitat type for both grazed and ungrazed plants. The models were fitted using the raw data, so estimates of fruit counts were scaled up to the whole shrub. These data were used to construct the appropriate stage structure.

The models were fitted with ten stages. Stages one and two were non-reproductive (seeds and saplings), the following eight stages were reproductive and structured to account for the initial rapid changes in fecundity with age (Figs. 1 and 3). Reproduction was set to start at year 6 and fecundity was capped at the value that was estimated for individuals of age 25, according to our observed data. Beyond this age, individuals may keep ageing, but their fecundity did not increase in our simulations.

Fig. 1
figure 1

Flow chart of stage structure used; numbers show the ages in each class. Stages 1 and 2 are non-reproductive

Full size image

Survival data were not available for Pyracantha therefore survival values were estimated based on studies of similar species. Seeds have a low probability of survival to seedlings. For example, for Acacia tortilis in arid and semi-arid habitats, annual seed survival was less than 0.001% (Hegazy and Elhag 2006). (Crisp 1978) found similarly low survival probabilities for three semi-desert shrubs. These species also have consistently high rates of mortality in the transition from seedlings to reproductive adults, but mortality declines and remains constant in the middle of the lifespan before rising again towards the end of the life span (Crisp 1978; Hegazy and Elhag 2006). Survival parameters were therefore estimated to be similar to those of the species examined by Crisp (1978) and Hegazy and Elhag, (2006). An elasticity analysis showed that a 20% reduction in annual survival probabilities for each age class always resulted in a < 9% change in wave-speed predicted using the analytical model (Fig. 6).

Seed dispersal was modelled with a 2Dt kernel where the probability density of seeds in a given location (x) is a function of dispersal distance with scale parameter (u) and shape parameter (p) (Clark et al. 1999):

$${f}_{(x)}= \frac{p}{\pi u{\left(1+\frac{{x}^{2}}{u}\right)}^{p+1}}$$
(1)

A dispersal kernel for Pyracantha was unavailable, therefore this model was parameterised using the scale and shape parameters derived for Ligustrum lucidum in Yungas forest (Powell and Aráoz 2018). These data provided information for incorporating rare long-distance dispersal events (Eq. 2).

$$f\left( x \right) = 0.9 \times \left( {\frac{{p_{1} }}{{\pi u_{1} \left( {1 + \frac{{x^{2} }}{{u_{1} }}} \right)^{{p_{1} + 1}} }}} \right) + 0.1 \times \left( {\frac{{p_{2} }}{{\pi u_{2} \left( {1 + \frac{{x^{2} }}{{u_{2} }}} \right)^{{p_{2} + 1}} }}} \right)$$
(2)

where parameters p1 = 5.58, p2 = 0.27, u1 = 11.34, and u2 = 7.59.

This is an unrelated species, but one with similar dispersal characteristics. In northern Argentina the two species often grow in the same habitat and produce fruits at the same time when few native fruits are available. They are therefore likely to attract the same species of birds (Tecco et al. 2006). We observed Turdus spp. eating fruits of Pyracantha, and species in the Turdidae family have been observed dispersing L. lucidium and Pyracantha when the two species grow together (Vergara-Tabares and Rojas 2016).

Individual based model (IBM)

We used a modified version of RangeShifter v1 (Bocedi et al. 2014) to build spatially explicit, stochastic IBMs to simulate the invasion of Pyracantha through homogenous landscapes. The RangeShifter platform was extended for this study to enable the inclusion of the 2Dt dispersal kernel.

A corridor of 20 cells (200 m) wide was set up and initially Pyracantha was introduced to the bottom five rows with the remaining cells empty. Carrying capacity was set at the average density of shrubs observed in the field and cells initialised at half the carrying capacity. Each simulation year an individual increased in age and size, and once the minimum age for reproduction was reached, its fecundity was determined by its age (as determined by the stage structured population projection matrix) and its seeds were dispersed according to the specified kernel (Eq. 2). Individuals transitioned to the next life stage only once the specified age was attained (Fig. 1). The model assumed no dispersal from outside the simulated environment, but seeds could disperse outside the simulated environment, in which case they died. This model is stochastic, so simulations were repeated 10 times for each parameter combination.

For each replicate and treatment, the invasion speed was calculated as the average distance moved per year by the invasive front at 50, 150 and 250 years from the beginning of the simulation. The density of the invasive population was calculated by dividing the total number of individuals by the area occupied after 150 years, at which point equilibrium had been reached (c.f. Travis et al., 2011).

Analytical model

We used the analytical approach described by Neubart and Caswell (2000) to model population spread (Eq. 3), which accounts for population growth and dispersal. Analysis was conducted in R version 3.5.2 using code based on the case study discussed by Neubart and Caswell (2000).

$$n\left(x, t+1\right)= {\int }_{-\infty }^{\infty }[K(x-y)\circ {B}_{n}]n\left(y,t\right)dy.$$
(3)

The population density (n) at each location x changes with time t and is estimated by integrating the process for all locations y. Here K(x–y) represents a probability density function (known as a dispersal kernel) that defines the probability of an individual dispersing from x to y (Eq. 4). This is expressed as a matrix of dispersal kernels for each stage in an age-structured population projection matrix (Bn,). Here dispersal only occurs for seeds dispersing from adults (represented as the 2Dt kernel—Eq. 2). When there is no dispersal associated with a stage the Dirac delta function (δ(x)) is used as this behaves as a probability density function and integrates to 1 so these individuals remain at their location and do not disperse.

$$k(x)=\left[\begin{array}{cccccc}\delta \left(x\right)& \delta \left(x\right)& 2Dt& \cdots & 2Dt& 2Dt\\ \delta \left(x\right)& \delta \left(x\right)& \delta \left(x\right)& \cdots & \delta \left(x\right)& \delta \left(x\right)\\ \delta \left(x\right)& \delta \left(x\right)& \delta \left(x\right)& \cdots & \delta \left(x\right)& \delta \left(x\right)\\ \delta \left(x\right)& \delta \left(x\right)& \delta \left(x\right)& \cdots & \delta \left(x\right)& \delta \left(x\right)\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \delta \left(x\right)& \delta \left(x\right)& \delta \left(x\right)& \cdots & \delta \left(x\right)& \delta \left(x\right)\end{array}\right]$$
(4)

Since the solution is impacted by the way in which the number of stages is specified (Travis et al. 2011), the demographic data used in the IBM is used to structure the population projection matrix (Eq. 5) in the general form for n age classes. This is based on the range of age classes specified for the 10 stages used in the IBM (Fig. 1).

$${B}_{n}=\left[\begin{array}{cccccc}{f}_{0}& {f}_{1}& {f}_{2}& \cdots & {f}_{n-1}& {f}_{n}\\ {s}_{0}& 0& 0& \cdots & 0& 0\\ 0& {s}_{1}& 0& \cdots & 0& 0\\ 0& 0& {s}_{2}& \cdots & 0& 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0& 0& 0& \cdots & {s}_{n-1}& {s}_{n}\end{array}\right]$$
(5)

Here Parameters fn are the per-capita fecundities of plants and sn are the survival rates of plants of each age (See Appendix 3 for values used for each model treatment).

Population spread is simulated from the population projection matrix (Eq. 5) across a domain from left to right and forms a wave of constant shape. From this the equilibrium wavespeed (m year−1) is derived, which is directly comparable with the invasion speed calculated from the IBM.

Elasticity analysis

We used elasticity analysis to estimate the proportional change in wavespeed for a proportional change in a vital rate (fecundity and survival). We reduced the fecundity and survival parameters for each age class and treatment by 20% and ran the analytical model (as before) and recorded the wavespeed.

Results

Plant age estimation

Statistical models established that plant age could be adequately estimated as a simple linear function of basal stem diameter (Fig. 2; lm: F1,29 = 76.2, R2 = 0.7149, P < 0.0001) and the resulting model was used to predict age for all individuals in the data set. The resulting predictions range from 5.8 to 25.7 years with an average age (± SE) of 11.07 ± 0.21 years. These estimates are consistent with suggestions that Pyracantha was introduced roughly 30 years ago immediately following the start of urbanisation linked to tourism (Ricardo Grau, pers. comm.).

Fig. 2
figure 2

Relationship between basal stem diameter and number of rings—used as a proxy for age (years), the line shows the fitted model (y = 1.61x + 3.62)

Full size image

Environmental impacts

The best fitting model included habitat type and grazing as fixed effects (Table 1). Grazed plants had a lower fecundity than ungrazed plants (GAM main effect of grazing: F1,18.08 = 25.3, p < 0.005). Habitat type also influenced fecundity (GAM main effect of habitat: F2,18.08 = 3.2, P = 0.04); plants in shrubland had a higher fecundity for their age than those in grassland and rocky habitats (Fig. 3, Table 2). Including plant density and an interaction with habitat type was not supported by the AIC and therefore did not improve predictions of fecundity (Table 1). Fecundity increased with plant age and this relationship was nonlinear (grazed: edf = 4.66, p < 0.001, and ungrazed: edf = 3.20, p < 0.001): plants increased in fecundity steeply between ages of 5 years and approximately 10 years in both grazed and ungrazed environments, but then declined in fecundity with increasing age with a more marked decline in grazed environments (Fig. 3). Interestingly, the best fitting model included a smoothing term for plants in both grazed and ungrazed environments so supports the hypothesis that grazing changes the shape of the relationship between plant age and fecundity (Table 1). The fitted line for plants that have been grazed was asymptotic with age while for ungrazed plants there was no evidence that an asymptote in fecundity has yet been reached, therefore the largest differences in fecundity between grazing and non-grazing regimes are between older individuals (Fig. 3). Based on the fitted models for fecundity and estimated survival data a life table was constructed for Pyracantha (Table 2).

Table 1 Description of fitted models with associated degrees of freedom and AIC values
Full size table
Fig. 3
figure 3

The relationship between age and number of fruits for each treatment, the panels show each of the 3 habitat types and black represent grazed and grey not grazed. Fitted lines show the predicted lines from the best fitting model

Full size image
Table 2 Life table showing the fecundity values estimated for each treatment for each stage (ages within stages also shown) and the estimated annual survival rate
Full size table

The habitats varied in specie plant density, grasslands had the greatest density of both adult shrubs and saplings (Table 3). The average plant age was consistent across habitats (Table 3).

Table 3 Average plant age, number of adults and number of saplings within each plot for the 3 habitat types
Full size table

Dynamics of invasion

The equilibrium wave-speed calculated from both the analytical and IBM models consistently predicted faster rates of spread in ungrazed habitats (Fig. 4). The IBM predicted a wave-speed that increased with time from the start of the simulations (Fig. 4). At all timepoints grazed habitats had a lower invasion speed. For example, the simulation predicted a range in wave speeds from 4.38 to 5.00 m year−1 in grazed habitats and from 5.06–5.47 m year−1 in ungrazed habitats after 150 years (Fig. 4). In addition, shrublands consistently had a higher rate of spread than in grassland and rocky areas. The analytical model predicted values that ranged from 3.99 to 5.22 m year−1 in ungrazed and 2.65 to 3.54 m year−1 in grazed habitats. The analytical model produced consistently larger differences between treatments than at all timepoints in the IBM, but the wave speed in grazed habitats was consistently higher in the absence of grazing in both methods.

Fig. 4
figure 4

Wavespeed (m year−1) for each treatment estimated by the analytical model and the IBM. Error bars for the IBM represent standard errors derived from 20 replicates

Full size image

Elasticity analysis showed that a reduction in parameter values by 20% resulted in small changes in the wavespeed (Fig. S6), ranging from 0 to 9%. Survival of young non-reproductive individuals has the biggest influence on wavespeed, but these effects are small (Fig. 5).

Fig. 5
figure 5

The mean proportional change in wavespeed (elasticity) averaged across treatments for individuals ranging from ages 0–25 responding to proportional changes in either the fecundity or survival parameters

Full size image

Invasion in the absence of grazing generated a higher equilibrium population density (after 150 years) than invasion in the presence of grazing, with a difference of approximately 1 individual per m2 (Fig. 5). Equilibrium population density was highest in the shrubland habitat (Fig. 6).

Fig. 6
figure 6

Average plant density (individual per m2) for each treatment estimated with the IBM. Error bars for the IBM represent standard errors derived from 20 replicates

Full size image

Discussion

The results show that both habitat and grazing influence the relationship between fecundity and age in Pyracantha. By quantifying this relationship and incorporating these changes into population models we have shown that projected rates of invasive spread and population density would vary across habitats and grazing regimes. We provide a first indication as to how these would vary (in form and magnitude). Most notably, reducing grazing pressure can substantially increase the speed of invasion. Further we identify that shrublands are more vulnerable to invasion than rocky or grassland habitats. This study highlights the importance of quantifying environmental and land use influences on demographic variables for predicting invasion speed. These results are useful for designing optimal control strategies by identifying vulnerable habitats where resources may be focussed to effectively target invasion by Pyracantha.

It must be noted that we estimated seed production at one point in time during the fruiting season, so the results provide a relative index of fruit production in local conditions rather than estimates of total lifetime fecundity. Nevertheless, we show that grazing reduces the realized fecundity of Pyracantha which has consequences for the wavespeed and density. This could be due to a reduction in plant photosynthetic area and plant vigour due to browsing damage or the direct consumption of fruits (Strauss and Agrawal 1999; Barton and Koricheva 2010). For example, cattle were regularly observed consuming the fruits of adult Pyracantha despite their robust thorns. Our results suggest that browsing may have a larger impact on the fecundity of older plants. However, grazing on younger plants may play a critical role in slowing their transition to larger plants that produce more fruits. Further work targeted at better understanding the mechanisms by which herbivory reduces fecundity in Pyracantha would be informative and provide the opportunity to explore alternative grazing control scenarios. For example, it would be useful to know whether a focus on moving cattle to areas in the valley invaded by the plant during the fruiting period would be more beneficial than having the grazing in those areas spread more evenly throughout the year. In addition, grazing may also be effective at reducing survival of Pyracantha saplings. The elasticity analysis showed that invasion speed was most vulnerable to changes in survival of young individuals, which suggests that increased browsing pressure of young Pyracantha may help to reduce wave speed in grazed areas.

The two models employed in this study provided similar estimates of relative fecundity between treatments. However, the wave speed predicted from the IBM increased with time. This suggests that the spread of the invasion front is, in part, driven by seeds from the core; as time increases the number of older individuals in the population increases, and the chance of rarer long-distance dispersal events increases (Harris et al. 2011). These older and larger individuals are likely to be found towards the invasion core and this may drive accelerating expansion of the population. In contrast, the analytical model only used information about the dynamics at the front of the invasion. Effectively, it implicitly assumes a pulled front whereas the IBM can incorporate pushed front dynamics too (see Gandhi et al. 2016 and Lewis et al. 2015 for discussion on pushed and pulled invasion dynamics). However, the conditions, modelled by the IBM assumed infinite habitat and space for spread with no competition from other vegetation. In reality, the invasion speed would reach an equilibrium as the population becomes limited by elevation, available suitable habitat, herbivory and competition. For example, Pyracantha seedlings survival is limited above 2400 m (which is consistent with the maximum elevation of our plots), therefore future models in real landscapes could include elevation as a limiting factor for seedling survival. Notably, the results of the two models are similar after 150 years of simulation; this is also when the population density in the IBM reaches an equilibrium.

As we noted in the introduction, there are multiple reasons why spread rates of invasive species are likely to differ between habitats. Here, our focus has been on understanding the contribution that changes in fecundity in different habitats are likely to make to spread rates in Pyracantha. Future work would usefully extend this to consider how dispersal and establishment success vary across habitats. In terms of dispersal, several studies have found that dispersal, including of invasive plants, varies across habitats, particularly for bird-dispersed species (Medellin and Gaona 1999; Carlo et al. 2013; Powell and Aráoz 2018). The presence of perches influences the distance travelled by birds (Medellin and Gaona 1999; Powell and Aráoz 2018). Thus, we might anticipate that dispersal of seeds into the shrubland habitat would be accelerated relative to the other habitat types due to the increased availability of perches. Additionally, the relative absence of perches at the front of an invasive shrub species spreading into a grassland habitat might be limited by the absence of perches: a form of Allee effect (Medellin and Gaona 1999; Powell and Aráoz 2018). The likelihood that individual seeds are consumed by the bird and dispersed may depend upon the abundance of seeds on a plant and in the neighbourhood. It is plausible that birds focus their foraging on areas of high fruit abundance thus resulting in strong heterogeneities in the likelihood that individual seeds are dispersed as well as where they are moved to. A further level of interest is provided in the system by the presence of other fruiting plants that also influences the behaviour of frugivorous birds. Notably, during the winter months in this region of Argentina the fleshy-fruits provided by exotic species, including Pyracantha, are the only food-source for frugivorous native birds (Rouges and Blake 2001; Blendinger et al. 2012). This contributes to increased inputs of seeds in invaded areas and could explain the observed association among exotic plant species dispersed by birds (Tecco et al. 2006; Vergara-Tabares and Rojas 2016). There is obvious potential for the invasion of one exotic species to facilitate the spread of others by providing additional fruit resource to support a population of frugivores, by providing additional perches from which seeds may be deposited and by acting as nurse plants and providing potential safe sites for seed germination and establishment. Therefore, the impact of other exotic species on the behaviour of frugivores and the consequences for dispersal would be interesting to quantify and incorporate into models of invasion spread in heterogeneous landscapes for Pyracantha and other invasive species that occupy the same habitats.

Grazing by domestic livestock may have complex effects on the dynamics of invasion. In some cases, grazing may open established vegetation to colonisation by invading species (Relva et al. 2010) and reduce the survival of seedlings of more palatable species (Marcora et al. 2013). If the seed is not damaged, cattle can act as a dispersal agent and disperse seeds for long distances (Chuong et al. 2016). For example, in California seeds of Pyracantha are viable after ingestion by coyote (Canis latrans) or cattle (Silverstein 2005). Our models only considered dispersal by birds, but this highlights the potential for mammals to act as dispersers. However, grazing may also influence the dynamics of invasions if the grazers consume the seeds or fruits of non-native species. The influence of grazers on the population dynamics of invasive plants may be in part determined by whether seeds survive ingestion by grazers, although this is a relatively unexplored topic in invasion biology. Our results highlight that any change in grazing regimes, driven by changing patterns of human settlement or demography, may influence the dynamics of invasions.

Furthermore, urbanisation may be associated with an increase in planting of Pyracantha, which may enhance propagule pressure and lead to an increase in invasiveness especially in areas near urban areas (Busnelli et al. 2006; Dawson et al. 2015). In addition, Pyracantha resprouts with a high vigour after fires which could have important consequences in ecoregions where wildfires are frequent or where burning is used in management (Herrero et al. 2016). However, the response to fire may depend on the habitat. The presence of large rocks results in an increase in the post-fire resprout size of Pyracantha (Herrero et al. 2016). There may also be an interaction between fire, grazing and urbanisation that could have complex impacts of the invasion dynamics of Pyracantha. For example, in grasslands the influence of fire on woody species is reduced where there is grazing resulting in an overall increase in woody species (Briggs et al. 2002). Therefore, if future studies aimed to use this approach in ecoregions where fires are frequent it may be important to quantify this interaction and incorporate this into models.

Our results highlight that ignoring habitat-dependent heterogeneities in age-specific fecundity regimes can lead to inaccurate estimates of the potential for spread. We show that expansions occur more rapidly and give rise to denser populations in a shrubland habitat than in grasslands or rocky habitat. They suggest that the shrubland areas represent the most vulnerable fronts for invasion. Therefore, our findings are useful for designing monitoring surveys and management plans by identifying conditions where rapid and dense invasions are more likely. However, it is unclear what the distribution of effort (plant removals, managed grazing, etc.) between the different habitat types should be. Therefore, models could be extended to examine the optimal distribution of management effort between the habitat types in order to minimise total area occupied by the invasive species, total abundance of the invasive species, or both. Such work could also potentially account for the value of the different habitats in terms of biodiversity, ecosystem service provision or farm income.

Furthermore, in a management context, we highlight that relaxed grazing pressure due to reduction in the number of cattle farms (Busnelli et al. 2006) could increase the projected spread rate by up to 53% due to the impact on fecundity alone. We emphasise that fecundity is just one route by which heterogeneity in habitat or grazing regime can influence spread rate. Future studies are required to gather the data to ascertain how dispersal and establishment vary according to local conditions and then modelling is needed to determine the relative importance of habitat-specific fecundity, dispersal and establishment on spread rates, including considering the potential for interactions between them. With this information at hand, it will be possible to use models of Pyracantha in complex real landscapes to provide robust advice on optimal management options for Pyracantha in the region under the continuing context of upland urbanisation.