Simulations

The Miscanthus simulation model we used [18] allowed a founding population of an invader to begin in the center of the landscape and then grow through local reproduction and dispersal. Our simulations were similar to those in the original model in that the habitat type at each position in the landscape, which acts as a proxy for both abiotic habitat quality and the resistance to invasion by the extant community, affected the population growth rate and carrying capacity of the invader; however, we also allowed for the habitat type to influence the probability of dispersing seeds establishing in an empty patch. There was no reduction in the best habitat (grasslands) and the probability of establishment was reduced in other habitats proportionally to the relative values of their habitat expansion rates (r_i; the rate of population growth for an individual patch of a given habitat type).

We ran simulations on 1000 landscapes that were randomly selected from the METALAND database [29]. METALAND utilizes data from the US National Land Cover Database (NLCD; [30]) and splits the entire country into a grid of 6.48km x 6.48km landscape tiles. Each tile provides a landscape of 216 x 216 cells, with each cell being 30m x 30m. The only constraints on landscape selection were that the entire landscape was classified (i.e., landscape tiles at the edges of the NLCD coverage area that had empty pixels were excluded) and more than half of the landscape was terrestrial (i.e., not open water). While we used the same model structure as in our previous work, we allowed specific parameters to vary in order to represent hypothetical invaders with similar properties (such as fecundity and relative preference/suitability for different habitat types), but that encompassed a range of levels of certain traits associated with invasiveness. In particular, we varied values of 4 traits: dispersal ability, establishment ability, population growth rate, and ability to use dispersal corridors. To vary dispersal ability and population growth rate, we modified the mean values of the log-normal dispersal kernel and the habitat-specific expansion rates to range from one-third to double the values included in the original model (0.33, 0.5, 1, 1.5, 2). To vary establishment ability, we allowed the frequency with which accessible patches arose to be 0.005, 0.01, 0.025, or 0.05 (the original model used a value of 0.03). Finally, we modified corridor usage by changing the proportion of seeds that actually move along a corridor (given that they landed in corridor cells, here defined to be roads) from 0 to 100% (0, 25, 50, 75, 100%); this value was 50% in the original model.

Simulations for each unique parameter combination (500 total combinations) were run on all of the 1000 landscapes, with initial populations of 50 2^nd year plants in each of the 9 cells at the center of the landscape. Simulations were run with annual time steps to simulate 100 years of invasion dynamics. This provided 500,000 total simulations, each with 100 time steps. To capture different aspects of invasion dynamics, four separate response metrics were calculated. To evaluate short term, transient dynamics that reflect earlier stages of invasion, we recorded the yearly population growth rate averaged across the landscape and calculated the mean growth rate over the first 30 years of the simulation (growth rate). Additionally, to evaluate the rate of spread, we recorded the yearly increase in the proportion of landscape cells that had at least one adult plant and the mean across the first 30 years was calculated (expansion rate). Finally, to consider longer-term dynamics, we recorded the final population size averaged across the landscape at the end of the simulations (final population) and the proportion of the landscape with at least one adult plant (final extent); however, even at the end of 100 years, simulations did not reach equilibrium conditions.

Analysis of landscape characteristics

For each landscape we collected the habitat data from the original METALAND landscape but reclassified the map to account for habitat types that we combined [18] and calculated the proportion of the landscape in each of 11 categories (grassland, pasture, crop, deciduous forest, mixed forest, coniferous forest, wetland (riparian), shrubland, road, water, and unsuitable). Using these reclassified landscapes we calculated landscape-level spatial statistics for area, edge, shape, aggregation, and diversity metrics using FRAGSTATS and details of specific metrics can be found in the related documentation [24]. However, we did not calculate core area or contrast metrics because we had no strong a priori standards for edge depth or relative strength of contrast between habitat types. These metrics offer a detailed but generalizable description of many aspects of landscape structure including fragmentation, habitat heterogeneity and diversity, patch shape, and landscape connectivity [24,31].

To reduce the number of modeled landscape metrics in our analysis, we removed all metrics that showed little to no variation across the dataset; these were essentially trivial characteristics (e.g., total area, which was constant across all landscapes). Additionally, for characteristics that were based on a distribution of values from across the landscape (e.g., patch area) we limited the number of parameters used to describe the distribution to three. We only included median, area-weighted mean, and standard deviation metrics to capture information about the central tendency and variability of all patches, as well as characteristics of the most influential patches, while limiting parameters likely to be highly correlated (e.g. median and mean, or standard deviation and coefficient of variation). This led to a total of 50 landscape metrics, 11 of which reflected the composition of the landscape (the abundance of different habitat types) and 39 that described landscape-level spatial structure and patterning; these are hereafter referred to as composition and spatial structure metrics, respectively.

Statistical analyses

We conducted two levels of statistical analyses: 1) using a Bayesian lasso regression [32] we estimated the influence of each landscape characteristic on invasion for each unique combination of invader traits and response metric, and 2) using the parameter estimates from the lasso regressions, we estimated how the relative importance of landscape characteristics change with respect to the strength of specific invasive traits. The use of lasso regression provides an efficient method to deal with many highly correlated variables for which there are no clear expectations as to relative importance by allowing a shrinkage parameter to reduce the influence of uninformative variables [33].

We used a fully Bayesian approach to generate Markov Chain Monte Carlo (MCMC) samples of the posterior distributions of regression parameters. To fit the models, we used data from 900 of the simulations, and held the remaining 100 out for validation. Using a standard, generalized linear model framework with vague priors, the posterior of the regression model was: (Eq 1) where y is the values of the particular invasion metric of interest and β₁ to β₅₀ indicate coefficients for each of the 50 landscape characteristics which are quantified for each unique landscape in the matrix X. For each data set and response variable, we centered and scaled all independent variables and then ran three MCMC chains for 10,000 iterations after a burn-in of 5,000 iterations. We checked chains for convergence by visually checking the trace plots and calculating Gelman and Rubin’s convergence statistic. While sampling the posteriors for each model, we recorded the ability of the model to predict out-of-sample data (i.e., the results of simulations on the 10% of landscapes not included in the regression). This process was repeated for 24 different values of the rate parameter of the regression coefficients’ Laplace prior (τ ranging from 0 to 2,000); this rate parameter controls the degree to which parameters are shrunk towards 0. By plotting the predictive ability of models fit at each level of τ, we were able to determine the optimal degree of parameter regularization without imposing arbitrary model selection on our independent variables. All analyses were conducted with JAGS (Just Another Gibbs Sampler) using rjags [34] in the program R, version 3.1.2 [35].

After estimating β coefficients for each landscape characteristic at each unique combination of invader traits with the lasso regression we then evaluated how the importance of those landscape characteristics for predicting invasion dynamics changed with invader traits. We used the posterior means of parameters from the optimum lasso regressions as response variables in a standard Bayesian generalized linear model with the values of each of the invader traits as the independent variables in the regression. Here, our posterior was simply: (Eq 2) where is the estimated coefficients of given landscape characteristic from above, λ is the intercept for the linear predictor, and Φ is the coefficient for each of the invasiveness traits and Z is a matrix of the values of each of the invasiveness traits in each set of simulations. As with the lasso regression, we ran three MCMC chains for 10,000 iterations after a burn-in of 5,000 iterations and checked for convergence using standard methods. We calculated the mean and 95% credible intervals of the posterior distribution of each regression coefficient and considered Φ whose credible intervals did not include 0 to indicate a meaningful relationship between the plant invasiveness trait and the importance of the landscape characteristic on invasions.

Specific interactions with landscape characteristics for each plant invasiveness trait

The role of simulation models in evaluating invasion risks

Our simulation modeling approach allowed us to address the role of landscapes on invasion dynamics. We were able to simulate invasions on real, remotely sensed, landscapes where the invader’s population dynamics varied among habitats; this provided the basis for evaluating complex spatial dynamics. As a stochastic variation of an integro-difference equation (IDE) model in a complex landscape, rather than an analytical model (i.e., classic IDE models), the influence of landscape structure on invasions can be seen as an emergent property of the landscape. By using a model-based approach to understanding invasion dynamics, we were able to 1) run simulations across a wide range of landscapes to produce large datasets and 2) replicate invasions in the same landscapes with invaders that have slightly different traits. This combination allowed us the opportunity to look at both how different landscape characteristics influence invasions and how the importance of particular landscape characteristics in predicting invasions changes for invaders with different ecological traits. We were also able to do this in ways that would not have been possible with experimental or observational methods. Additionally, by using a Bayesian model regularization approach [32], we were able to overcome many of the challenges of managing datasets of numerous landscape metrics, particularly those relating to highly correlated data.

Our use of dynamic simulations was also critical to identifying the influence of landscape structure on invasions as this was more clearly seen in invasion dynamics rather than long-term patterns. Much of the recent effort in predictive modeling of invasive species has been focused on identifying and then quantifying the abundance of high-quality habitat, particularly with species distribution models (SDM) or habitat suitability models (HSM) [36,37]. To a degree, our analysis supports this focus as we also found habitat quality to be the most informative predictor of invasion dynamics. However, the fact that spatial structure characteristics were also informative and that relative importance changed with ecological traits of invaders suggests that a sole focus on habitat is a limited approach for predicting invasion dynamics. Our ability to identify these more complex influences of landscapes on invasion dynamics is in part because we explicitly considered both short-term and long-term dynamics. There has been a growing appreciation that invasive species violate the equilibrium assumption embedded in SDM and HSM techniques, and approaches to adapt to this issue are being developed [37,38]. Yet, even these efforts only focus on the lack of equilibrium when determining habitat suitability; predictions made by most models still assume complete penetration of invaders and thus spatial or temporal dynamics inherently become inconsequential. For most management and decision making purposes, this is not a reasonable assumption. There is a need to consider short-term dynamics because they can have important implications for risk analysis and design of management or best practices strategies [39,40]. Short-term dynamics may be the most important determinants of whether spread can be prevented and as such may be paramount concerns in the case of potentially invasive biomass crops, as society may be unwilling to tolerate any substantial risk of widespread establishment of such species.

Invasions are complex emergent processes with multiple drivers

The variability between multiple response metrics shows that invasions are complex processes, which cannot be described or compared by a single measure. Invasion dynamics are the emergent outcome of numerous interacting processes. This complexity makes them difficult to predict, but suggests that there may be a variety of options for management and intervention, and that a focus should be placed on different targets depending on the particular goal [41–43]. While the long-term outcome of most successful invasions is extensive spread across the entire landscape, patterns over short-term timescales might be more practically relevant to questions of management. Over short timescales, invasions will certainly be less spatially extensive, but are also more constrained within certain habitat types; identifying those differences can inform where monitoring and eradication efforts should be focused. Similarly, concerted efforts on certain landscape features (e.g., roads that act as dispersal corridors or habitat islands that can act as reservoirs) can have an outsized influence on the speed and extent of spread [44]—this insight has also influenced management strategies for marine and aquatic invasive species [45,46]. Depending on the particular ecology and/or impacts of an invader, the importance of managing for spatial spread versus population size may be very different. For example, managing dispersal of species that have little impact at low abundance may be less important than focusing control efforts on the habitat patches where populations can grow to harmful levels. In contrast, the most effective strategy for species that have strong local impacts but have low natural dispersal capability could be to prevent access to dispersal corridors. Relatedly, our results indicate that changing the strength of invasiveness traits does not influence all invasion metrics equally. As might be expected, increased dispersal ability tends to influence spatial extent of invasions more than total population size (though both increase together, Fig 1). Increasing population growth rates display the opposite pattern with stronger effects on population sizes.

In general, we observed segregation in response patterns between landscape metrics that described the composition of habitat types versus metrics related to the spatial structure. Broadly, the effects of landscape composition characteristics were generally stronger than those of spatial structure metrics, similar to patterns that have been identified by others [47]. However, we also observed that for transient dynamics, the importance of spatial structure increased as invasiveness increased (Fig 4), while that of composition decreased, this suggests that analyses without consideration of such an interaction could be misleading. Overall, this pattern indicates that at low invasiveness, landscape pattern has less influence than composition on the dynamics of invasions. This is not because of an overall cessation of invasion, as invaders are still spreading. Rather, when invasions are slow the species may be unable to explore much of the landscape and only the habitat directly adjacent to the founding colony matters. As invasiveness increases and spread is greater, the composition and arrangement of habitat begins to play a larger role; however the particular means by which invasiveness increases (e.g., through higher propagule pressure from rapid population growth, greater dispersal distances, or higher seed survival) have greater impacts in certain types of landscapes versus others. This perspective is bolstered by the fact that these patterns are reversed when long-term metrics are considered. At these time scales, invaders likely have explored the entire landscape so compositional characteristics again take primacy. This question of landscape structure and its relationship to invasions has been considered before, generally showing that heterogeneity increases invasibility [15] and that landscape characteristics that are associated with disturbance can also be indicators of invasion [48]. However, this work has been limited by the relatively few types of landscape characteristics considered and because it generally has not considered how that relationship might change for invaders with different natural histories (but see [19]).

The basis for an interaction between spatial structure and specific ecological traits of invaders is clearer when there is an explicit focus on dynamics rather than equilibrial outcomes. Growth rates, fecundity, and dispersal do not impact carrying capacities, but they can influence how efficiently invaders spread through a landscape [49]. This is the aspect of invasions that spatial structure can most strongly influence; however, we observed strong contingencies in the effects of spatial characteristics, many of which had strong effects only in concert with particular invasiveness traits, indicating that the influence of spatial context is related to each trait. For example, we saw qualitatively different patterns between invasiveness traits. The spatial characteristics whose importance changed in concert with dispersal ability were largely related to fragmentation, as might be expected. Changing establishment ability most strongly influenced the importance of characteristics related to the presence of large, regularly shaped patches. This suggests that when an invader is very good at establishing, the most effective approach for spreading is from high propagule pressure produced by large local populations. In contrast, the ability to take advantage of dispersal corridors reduces the need for high propagule pressure to support spread, provided corridors are available. Thus, the importance of high-quality habitat can be relatively lower while patch characteristics indicative of corridors (abundance of road habitat, as well as parameters indicating long, linear patch shapes) have a larger effect.

Implications for management of invasive species

Inferences of management implications from an exploratory modeling study must be cautious. However, we highlight several implications that emerge from our results, though much additional analysis would be needed before use in practice, and these are best regarded as indications of how our approach to modeling landscape effects on invasion could provide substantial insights for management. Overall it seems apparent that the presence of high suitability habitat is the central risk factor for the spread of Miscanthus invasions, but spatial structure can limit—or exacerbate—the early stages of invasion, particularly in highly-invasible habitats. If preventing invasion in such habitats is important, then manipulations of landscape spatial structure can provide some measure of control. By evaluating invasion risks in individual landscapes we see that grasslands with facilitative spatial structures present the worst-case scenarios for Miscanthus invasion, particularly in terms of the rate of spread. We observed that many of these landscapes have very regular and consistent habitats, which may indicate that Miscanthus can move more efficiently through large single patches of habitat. Also, those landscapes often contain roads or other linear elements that can act as dispersal corridors. Conversely, many landscapes where spatial structure had relatively strong negative effects on rates of spread had complex spatial patterning featuring many small patches of habitat, suggesting that complex fragmented landscapes have slower spread rates. This finding suggests that management practices that create such landscapes might be useful interventions in particularly favorable landscapes for Miscanthus invasion, such as grasslands with large, rectangular habitat features. In practice, fragmentation might be introduced with long, narrow patches of a contrasting habitat type, e.g., so-called prairie strips that are being used to add habitat diversity to annual row-crop landscapes in Midwest US [50]. Of course, these management implications only apply specifically to Miscanthus and the invasiveness variants considered in this modeling study. Our results also argue that the influences of landscape structure vary with respect to specific traits of the invader and thus different landscape patterns may be better optimized for limiting the spread of other invasive species. Thus, it is important to recognize that manipulation of spatial structure, as a management option for limiting initial spread of invasive species, is highly tentative and requires examination across invasiveness traits and habits.