Biological traits and threats interact to drive extinctions in a simulation study

How a particular threat influences extinction risk may depend on biological traits. Empirical studies relating threats and traits are needed, but data are scarce, making simulations useful. We implemented an eco-evolutionary model to analyse how five threat types influence the extinction risk of virtual organisms differing in body size, maturity age, fecundity, and dispersal ability. Results show that direct killing mostly affected slow-living and low dispersal organisms. Habitat loss and fragmentation both affected larger and less fecund organisms, but drove contrasting responses according to dispersal ability. Habitat degradation and the introduction of invasive competitors had similar effects, mostly affecting large, fast-living, and highly fecund organisms. Many of the reported results confirm previous studies, while others were never tested, creating new hypotheses for future empirical work. Statement of authorship FC, LC and PC designed the study, FC implemented the model and ran the statistical analyses. FC and PC wrote the first draft, and all authors contributed substantially to further revisions.


Introduction
Species are going extinct at unprecedented rates in human history (Vos et al. 2015). Understanding what drives species to extinction is crucial if further extinctions are to be minimized approaching natural, pre-human, levels. Yet, the mechanisms leading to extinction, both intrinsic (species traits) and extrinsic (environmental threats), and the way they interact remain obscure (Chichorro et al. 2020).
Much progress has been made recently, with over 150 studies trying to identify similarities in traits between threatened species in the last two decades (Chichorro et al. 2019). These studies were both comparative and correlative in the sense that they compared traits of threatened species with those of non-threatened species, often using the IUCN red list categories (IUCN 2020) as proxies of extinction risk (Purvis et al. 2000;Cardillo et al. 2005;González-Suárez et al. 2013;Marco et al. 2015).
A recent study (Chichorro et al. 2020) tested, for the first time, multiple traits related to extinction risk across a large part of the Tree of Life and the entire world, minimizing functional, phylogenetic, and spatial biases of previous studies. Using data from close to 900 species of vertebrates, invertebrates, and plants, it was found that five traits, namely habitat breadth, fecundity, dispersal ability, altitudinal range, and degree of human influence were universal predictors of extinction risk.
The response of five other traits, i.e. body size, offspring size, generation length, microhabitat, and change in the degree of human influence across the species' range was dependent on the taxon, often on the threat type affecting each taxon. It showed large promise for future studies that the mechanisms driving their similarities and differences across taxa and threat types could be elucidated.
Previous correlative studies indicate that external factors interact in important ways with the relation between intrinsic traits and extinction risk (Olden et al. 2007;González-Suárez et al. 2013;Murray et al. 2014). Habitat destruction, pollution, invasive species, overexploitation, and climate change (Secretariat of the Convention on Biological Diversity 2014) all drive extinctions, and depending on the driver, some traits might be particularly relevant in a species trajectory towards extinction.
In a "statistically ideal" world, we would have full data on species traits, external drivers, and potential mechanisms linking these to extinction. However, as models become more complex, the amount of data required to test hypotheses increases. Additionally, because observations are rarely fully independent (dispersal ability might be connected with body size, fecundity with offspring size, etc.), it is often difficult to partition this complexity in independent components.
Conceptual models have been suggested to complement correlative models in order to better understand mechanistic drivers of diversity patterns (Murray et al. 2014). Within these, individual-based models (or agent-based models, henceforth referred to as ABMs) are increasingly used to aid understanding of the functioning of species and ecosystems (DeAngelis & Grimm 2014).
Their power derives from the possibility of replicating patterns found at the population and community levels from phenomena occurring at the individual level (the level at which the basic deciding biological entity exists). Despite difficulties in specifying and implementing these models, today tools exist that facilitate both specification (Grimm et al. 2006(Grimm et al. , 2010 and implementation (Wilensky 1999) of ABMs.
In this study, we develop a conceptual ABM to evaluate the influence of different threats, namely direct killing, habitat loss, habitat fragmentation, habitat degradation, and invaders, on the risk of extinction of species presenting differing body size, maturity age, fecundity, and dispersal ability. The landscape includes organisms competing for resources, which they use to invest in survival, growth, and reproduction. After creating virtual species adapted to given simulated worlds, we apply the threats and observe how the average trait values of populations respond to each threat. These results from in silico experiments are related to real-world knowledge on different organisms and reveal the mechanisms behind species extinction. Additionally, we hypothesize mechanisms to be 4 tested in the future, which can guide data collection, and can be confirmed or disproved as more empirical data become available.

5
Materials and methods The model The purpose of the model was to analyse the shift in mean trait values of an evolving population of organisms competing for resources, when an external threat event was introduced into the virtual world where they lived (Fig. 1). We implemented it as an ABM that simulated a landscape where organisms moved in the environment in order to obtain the required resources to survive and reproduce. At initialization, organisms presented a random set of traits. The population then converged towards normally distributed trait values around a stable mean. Each stable population (after convergence) was then treated as an experimental unit. We subjected each experimental unit to a range of threats at several intensities. As the population converged to a new set of trait value distributions, we registered the variation in mean trait values occurring under each threat type and intensity.
The ABM model was developed in Netlogo 6.1.1 (Wilensky 1999). A full description of the model following the Overview Design Concepts and Details (ODD) protocol (Grimm et al. 2006(Grimm et al. , 2010 is available in Supplementary Materials S1.

The environment
The environment consisted of a continuous 2D space in which organisms foraged in the form of an underlying toroidal matrix of 33x33 patches ( Fig. 2A). Patches were classified as either resource patches or bare patches. The resource patches generated resources linearly every timestep until they reached a maximum amount of resources. The quantity and spatial distribution of resource patches were controlled at the beginning of the simulation (Supplementary materials S1). 6 The virtual organism Organisms were mobile entities that reproduced asexually. They were characterized by an energetic currency, their age, their position in the continuous environment, and their focal traits. Organisms competed directly, accessing resources that fuel their maturation, reproduction, and survival needs ( Fig. 3; Sibly et al. 2013;van der Vaart et al. 2016;Brown et al. 2018). Excess energy was stored in reserves up to a maximum determined by body size . An organism dispersed whenever the intake of resources at the current timestep was lower than its needs for maintenance and maturation/reproduction. If the amount of energy in the reserves reached a certain ratio of max-energy ( energy-to-reproduce ) and the organism was an adult, it reproduced, and part of its energy was allocated to its offspring.
Organisms died whenever the energy in reserves fell below 0, if they were killed by the perturbation direct killing , or due to old age. The focal traits influenced the fitness of organisms by determining the order in which they foraged ( body-size ), the volume of energetic intakes and expenditures of organisms ( body-size and dispersal-ability ), how quickly they reached maturity and died ( maturity-age ), how quickly, and thus how far, they could move in the environment ( dispersal-ability ), and how efficient they were in leaving progeny and spreading their traits in the population ( fecundity ).

Body size
Being larger than other organisms competing for resources is a competitive advantage for both animals and plants through access to more resources (Kingsolver & Pfennig 2004;Hone & Benton 2005;Kingsolver & Huey 2008). On the other hand, energetic expenditure and acquisition scale positively with the body size of organisms (Brown et al. 2004;Glazier 2005Glazier , 2010. Therefore, in this model organisms executed their functions in an energy-sorted 7 order: organisms with more energy, which was positively correlated with body size, had higher chances of foraging first. Larger organisms intook and used energy at a faster rate than smaller organisms (Supplementary materials S1).

Maturity age
The ability to reach maturity earlier is accompanied by a reduced ability to correct mistakes when replicating DNA, which has multiplicative negative effects on cell functioning, limiting their ability to live longer (Kirkwood 1977;Dowling & Simmons 2009;Selman et al. 2012). In the model, longevity was proportional to the maturity age of an organism.

Fecundity
Investment in many offspring has a negative effect on each offspring's survival shortly after birth (Fox & Czesak 2000). In the model, organisms with higher fecundity generated offspring with lower initial energy.

Dispersal ability
Higher dispersal ability may impose energetic, risk, and opportunity costs on organisms that may be expended during the pre-dispersal phase, transfer period, or resettlement phase (Bonte et al. 2012

Threat phase
The threat phase started by subjecting a stable population to a threat. The threats in this model were built to mimic the main biodiversity extinction drivers (Table 1): direct killing, habitat loss, habitat fragmentation, habitat degradation, and invaders. The simulation ended when the mean trait values stabilized again (the stopping conditions test was called again). Each threat was applied following a gradient of intensity (only one per population, Fig. 1B). We set the maximum threat value as the value at which half of the populations per value went extinct before stabilizing in new trait values. We then created a gradient ranging from no threat (0) to the maximum value of each threat type in 20 steps. In total, we subjected 64 populations to each threat type and intensity.

Statistical analysis
The relative change in trait value ( ) was calculated using the following equation: where is the mean value of the trait k in the population when the simulation ended, V kf and the mean value of trait k of the initial population (before adding the threat). To V ki calculate the signal of the relationship across all replicate populations and to see whether this signal was significantly different from zero (thus positive or negative), we performed linear mixed-effects models for each (threat, trait) pair. Since each replicate population was run at different threat intensities (see subsection "Burn-in phase" above), we modelled the logRC within a linear mixed-effects model framework: where logRC ij is the j th observation of population replicate i , and population_replicate i being the random intercept and slope, assumed to be normally distributed with mean 0 and variance 2 . We forced the intercept through the origin because under threat intensity = 0, logRC is expected to be 0. In those (perturbation, trait) pairs where the relationship was significant, we recorded the signal of the relationship.
As parameter values of the standard analyses (e.g. initial number of resource patches, energetic intake of organisms, and mutation rates of traits), we chose values that generated both realistic scenarios, and that allowed trait values to increase or decrease after adding threats (values presented in the Supplementary Materials S2). In addition to the standard analyses, we ran robustness tests (Supplementary Materials S2) to determine whether 10 changing parameter values had a qualitative influence on the model outputs. We considered that the choice of values of a given parameter influenced the results if changing its value qualitatively modified the results, i.e., if we observed that the relationship between logRC and perturbation intensity was negative under a robustness experiment, while under the standard experiment values the relationship was positive, and vice-versa.
Thus, we considered that qualitative changes occurred if there was a reversal in the sign of the relationship between the logRC and perturbation intensity. The scripts to run the Netlogo model and statistical analyses are available in Supplementary Materials S3. 11 Results Body size and maturity age decreased and fecundity and dispersal ability increased with increasing levels of direct killing of individuals (Fig. 4). Results were consistent across all robustness analyses for maturity age and dispersal ability (i.e. there were no robustness analyses results that led to opposite effects on traits, Table S3.1), but not for body size and fecundity (Supplementary Materials S2).
As a response to habitat loss intensity, body size, maturity age, and dispersal ability decreased, while fecundity increased. All changes were consistent independently of simulation parameters, except for maturity age, which showed a positive signal under a range of robustness tests (Table S2.1).
Habitat fragmentation led to mostly similar outputs as direct killing (Fig. 4). Body size and maturity age decreased, and dispersal ability and fecundity increased. All patterns were consistent independently of simulation parameters, except for maturity age, which under certain values showed a positive response (Supplementary Materials S2). Although this threat also implied habitat loss, the change in trait values was greater under habitat fragmentation, revealing the direct effects of fragmentation beyond those of loss.
Maturity age and dispersal ability increased with increasing levels of habitat degradation, although there was an inflection point at extreme levels of perturbation (which did not qualitatively have an impact on the direction of change). Body size and fecundity decreased, but again with inflection points in the intensity gradient (Fig. 4). All patterns were robust to parameter values, except dispersal ability, which decreased under certain scenarios (Table   S2.1). 12 Adding competing invasive organisms with traits similar to the native led to increasing maturity age, and to decreasing body size and fecundity (Fig. 4). No change was detected in dispersal ability. All results were robust to parameter changes (Table S2.1). 13

Discussion
Our work confirms in simulated settings that different threat types act differently on organisms depending on their traits. The way that organisms respond to threats or the probability with which species can go extinct may therefore depend on the specific threats to which they are subjected. Most patterns found were robust to different parameter values, meaning that the underlying mechanism should be consistent across many ecological systems in nature. We should note that, being a single-species experiment, our simulations did not test species extinction directly. They did test the mortality of individuals with different traits that make the species evolve, in this way indirectly testing extinction.

Direct killing
Direct killing negatively affected large, slow-living, poorly fecund, and poorly dispersive organisms. Large sizes were removed from the environment because as organisms were killed at random, competition for resources diminished, which then reduced the need for the competitive advantage that larger size brings. Under some parameter value choices, however, the small organisms were instead disfavoured (Supplementary Materials S2). This seems to be due to the indirect effects of direct killing on resource availability. As organisms die, resources become more abundant. Abundant resources thus allow organisms to become larger under some parameter values.
Direct killing also heavily impacted organisms with slower life-cycles. Organisms taking longer to reach maturity faced greater probability of being killed before they were able to reproduce, and hence, the simulations became dominated by fast-living organisms.
14 Organisms with poor fecundity were disfavoured for most of the parameter values tested.
As direct killing affected all organisms with the same probability, individual survival became less relevant than producing more offspring to compensate for the mortality. In simulations where organisms had a high longevity-to-maturity age ratio (due to low maturity age), they tended to evolve towards low fecundity, and hence, highly energetic offspring. These offspring quickly reproduce, reducing the chances of being killed before reproduction.
Our results also suggest that direct killing disfavours organisms with poor dispersal ability.
This effect was again seemingly due to the increased resources available at each patch; as organisms died more frequently, more resources became available. These resources were better exploited by those organisms that found them first, or those with more dispersal ability.
The widely reported extinction of the megafauna before the 18th century has already shown that large-sized, slow-living, low fecund, and less mobile fauna are particularly susceptible, although it must be remarked that human hunting and fishing have long been selective towards large species (Diamond 1989). The same dual pattern seems to be occurring in recent times, at least for vertebrate taxa. Among fish and mammal species threatened by overexploitation, larger sizes tend to be more at risk (Olden et al. 2007;González-Suárez et al. 2013), a fact that is typically associated with their slow life-cycles.
Additionally, threatened mammals facing overexploitation show lower fecundity independently of size effects than those not threatened (Strauss et al. 2006;Warzecha et al. 2016). Other sources of direct killing, such as the presence of an alien predator species (Matsuzaki et al. 2011) or disease (Murray & Hose 2005), could be affecting organisms with the same traits that we found, but thus far studies examining interactions between threats 15 and traits have been limited (Murray et al. 2014). Interestingly, the increased vulnerability of small-sized individuals also emerged in a model evaluating the impact of predation on the size and maturity age of organisms, showing that prey size often increased when predation had an indirect effect on resource availability of prey (Abrams & Rowe 1996).

Habitat loss
Our model suggests that habitat destruction has a negative impact on organisms with larger sizes, slow life-cycles, low fecundity, and high dispersal ability.
The decrease in body size could be explained by the reduction of total resource area in the environment. As habitat loss decreased the number of suitable habitat patches, it decreased the area of continuous suitable habitat, threatening larger size individuals that require abundant resources. Lower population sizes due to resource scarcity also incur greater probability of local stochastic accidents (Pimm et al. 1988). Edge effects further exacerbate this due to the high spatial resource unpredictability seen through the organism's point of view.
Long-living organisms were also negatively affected by habitat loss. The higher spatial heterogeneity may negatively affect slow life-cycles, because these organisms attain reproductive age more slowly, and thus, reproduce less often, not compensating for mortality losses due to movement into the matrix of the unsuitable habitat. Under some robustness scenarios, though, habitat loss was detrimental to fast life-cycles. Typically, slow-lived populations survive local extinction better at lower population densities than 16 fast-lived populations (Pimm et al. 1988), and with less stochastic fluctuations, thus possibly conferring higher resilience in spatially heterogeneous landscapes.
Habitat loss intensity was detrimental to organisms with lower fecundity. High fecundity of organisms in the model came at the cost of decreased survival of individual offspring. When faced with a reduction of habitat, however, individual survival is offset by increased edge effects and spatial unpredictability.
Organisms with larger dispersal ability were negatively affected by habitat loss. As habitat decreases in size without alternatives to which to move, individuals that invest in long-distance dispersal incur higher mortality.
The impacts of habitat loss on species' traits have not been studied extensively.
Conceptually, larger sizes are thought to be the most susceptible to decreases in fragment area (Ewers & Didham 2006). However, empirical evidence supporting the higher risk faced by larger organisms is very limited, due to the existence of many confounding factors under the umbrella of body size (Ewers & Didham 2006); large size may correlate with other traits that increase (e.g. dispersal ability (Warzecha et al. 2016)) or decrease (e.g. trophic level (Henle et al. 2004)) the organism's capability of surviving in areas with increased habitat loss. Likewise, the effects of habitat loss on lifespan have not been conclusive. Plants with greater longevity were hypothesized to be better at coping with habitat isolation (Lindborg 2007), as they face less fluctuations in their population levels, but fast life-cycles have been shown to survive well in isolated fragments, likely because fast life-cycles were also characterized by being very fecund (Lindborg et al. 2012;Marini et al. 2012). In a grassland that has faced extensive habitat loss, surviving plants in isolated habitat patches had greater seed output than the original species population prior to the habitat loss event (Saar et al. 17 2012), which authors associated with their probable better ability to find microsites than species with lower seed output. A meta-analysis indicates that butterflies with poorer fecundity are more affected by habitat fragmentation (Öckinger et al. 2010). As for dispersal ability, evidence shows that when patches of remaining habitat are very isolated, organisms with reduced dispersal ability often cope better (Henle et al. 2004;Ewers & Didham 2006;Saar et al. 2012). However, isolation and no dispersal may have consequences for genetic drift and erosion for organisms reproducing sexually (Zambrano et al. 2019) and for species' ability to respond and shift their geographical ranges in the face of climate change (Hodgson et al. 2011).

Habitat fragmentation
Like habitat loss, habitat fragmentation had a negative effect on large, slow-living, and poorly fecund organisms, but a negative effect on poor dispersers. At the same proportion of habitat removal, the effect of habitat fragmentation on traits was stronger than that of habitat loss, revealing the importance of these two effects combined to drive population declines. Slow-moving organisms, being favoured under the pure habitat loss scenario were disadvantaged under habitat fragmentation, as they could not move among a complex matrix of suitable and unsuitable patches.
With increased fragmentation, organisms with higher dispersal that could move to new, relatively easy-to-find patches, were likely to benefit, also avoiding the dangers of low population size and stochasticity in very small fragments.
Increasing fragmentation has been shown to benefit good dispersers. For example, species of butterflies and moths in fragmented landscapes were characterized by greater dispersal ability (Öckinger et al. 2010). Differences in dispersal ability among landscapes can be found 18 even within species, whereas isolated patches in fragmented landscapes were more frequently visited by larger, more mobile wild bees (Steffan-Dewenter & Tscharntke 1999;Warzecha et al. 2016).
Studying the impact that habitat loss and fragmentation have on species traits is challenging. First, because often traits correlate with each other. Second, because the literature on the effects of habitat loss and habitat fragmentation on species' traits is fuzzy, since these terms are often used interchangeably (Tscharntke et al. 2012;Fahrig 2017;Zambrano et al. 2019). Third, because habitat loss and fragmentation can be subdivided into five components, each of which is hypothesized to have different effects on species' life-histories. These components comprise a decrease in fragment area, distance from the edge, shape complexity, isolation of the remaining habitat, and contrast between suitable and non-suitable habitats (Ewers & Didham 2006). For example, the decrease in fragment area is thought to negatively affect mostly the intermediate dispersers, whereas the long-distance dispersers can find suitable patches elsewhere and poorly dispersers do not incur unnecessary dispersal events that lead to increased mortality (Henle et al. 2004;Ewers & Didham 2006).

Habitat degradation
Habitat degradation in our model largely equates to resource scarcity. Larger organisms require more resources, therefore being especially vulnerable to reductions in resource quantity. In the model, organisms that mature more quickly survive less time as adults.
Under these settings, the chances that they may be able to accumulate enough energy to reproduce are smaller. Therefore, greater longevity is favoured, as it gives a chance to accumulate the energy required for reproduction. In the same way, by concentrating more energy onto only a few offspring, organisms increase the odds that their offspring will survive famine during the juvenile stage. Organisms that are poorer dispersers are disfavoured. As resources become less abundant, the necessity to rapidly find these elsewhere increases.
The effect of habitat quality as related to resource availability on biodiversity variables at a landscape level has been poorly studied, despite the fact that along with habitat loss and fragmentation, it is one of the major drivers of species extinction (Mortelliti et al. 2010 (Skogland 1985;Ferguson 2002). For instance, when food is scarce, pregnant wild deer abandon their fetuses in preference for their own survival (Skogland 1985). Plants living in nutrient-poor scenarios, usually have traits that confer low growth rates, such as high root-to-stem ratios and low specific leaf area (Grime 1977). Likewise, cave organisms survive extreme resource-depleted scenarios.
They are invariably long-lived, with low fecundity, and large organisms simply cannot survive entirely in caves due to lack of resources. All of these traits were favoured by our simulations, which replicate similar circumstances to those frequently occurring in degraded habitats. Our understanding of the mechanisms through which habitat degradation affects species is far from clear since, as we do not understand the types of traits that are selected (Mortelliti et al. 2010), even though habitat quality may typically outweigh the importance 20 of habitat loss and fragmentation to species extinctions, as suggested by other agent-based simulations (Heinrichs et al. 2016).

Invaders
The effect of invasives was somewhat similar to habitat degradation: a reduction in body size and fecundity, and an increase in maturity age (but no effect on dispersal ability).
The invasive organisms added competition for the same resources, reducing the number of resources available and causing similar outcomes for most traits. Organisms become smaller as a direct response to decreasing resource availability, compromising their ability to access resources first due to their smaller size than invasive organisms. To compensate for lower resource intake, they become longer-lived so as to accumulate resources over a longer period of time before reproducing.
Additionally, they maximize their offspring's survival by putting more energy into a few, rather than many, offspring. Invasive species had, however, no effect on the dispersal ability of natives. Under habitat degradation, the population density of organisms is reduced.
Therefore, higher dispersal may be useful to find other patches of suitable habitat. Such is no longer the case with the addition of an invasive species. The density of native and invasive organisms combined does not change much in relation to the pre-invasive event, and therefore, if an organism disperses to another patch it may find exactly the same density of organisms, curtailing the advantage of moving fast.
The impact of invasive competitors on the traits of native species competing for the same resources has not been examined thoroughly. A comparative study of extinction risk compared the traits of freshwater invasive, threatened native, and non-threatened native 21 fish species, and found that threatened natives had significantly faster life-cycles and higher fecundity than those of non-threatened native species (Liu et al. 2017). The authors suggested that the reduced fecundity and faster lives could be related to higher environmental stochasticity, and hence reduced resilience in the face of environmental change. River systems are especially threatened by the presence of invasive species; in many regions the number of invasive fish species surpasses a quarter of the total species richness (Leprieur et al. 2008), and in this system in particular they could be the number one cause of extinctions (Light & Marchetti 2007). Therefore, the faster life cycle of threatened fish species could be a consequence of the presence of alien competitors. A meta-analysis of trait differences between invasive and non-invasive plant species has shown that natives are on average smaller and have lower growth rates than invaders (Kleunen et al. 2010). Trait values of successful native plant species were however similar to successful invaders in another study of invasiveness across many habitats in Central Europe (Loiola et al. 2018), and reflect the fast growth capabilities needed to deal with increased disturbance in their environment (Leishman et al. 2010).

Conclusions
In this study, we explored how widely the examined drivers of extinction can differentially   28 Tables   Table 1: Types of threat studied, their description, and real-world equivalents.

Threat Description Examples
Direct

Purpose and patterns
The model is designed for theoretical exploration of biological trait value changes in response to five perturbation types: direct killing, habitat loss, habitat fragmentation, habitat degradation, and invasives. Specifically, we will explore changes in the values of four traits (body size, maturity age, fecundity, and dispersal ability) of virtual organisms moving and competing in virtual landscapes. From the resulting trait changes implications for extinction risk will be discussed. The generic patterns used to claim that the model is internally consistent and realistic enough for its purpose are previously published trends in trait-extinction risk relationships.

Entities, state variables, and scales
The model's entities are organisms, patches, and the environment. Their state variables are listed in Table S1.1. The model parameters are listed in Table S1.  One timestep has no real unit of time equivalent, but can be delineated by considering a "minimum resolution" which is the necessary amount of time for an organism to perform any of these functions: feeding, maturation (if juvenile) or reproduction (if adult), and movement. The total extent of time of a simulation is very variable, because it depends on how quickly the mean trait values of the population converge before, and after a perturbation. Convergence of trait values can be attained after only 500 timesteps, but typically simulations last between 1000-10000 timesteps.
The model world consists of 33x33 patches, which is large enough to generate maps with realistic spatial autocorrelation. We confirmed the robustness of our finding with both smaller and larger world sizes (Supplementary Materials S2).

Process overview and scheduling
Every timestep each organism executes the following submodels that are explained in detail in ODD section 7, "Submodels": maintenance , in which it deducts energetic costs of maintaining its body, and die if its energetic level becomes lower than zero; feeding , in 35 which it obtains energy from the patch where it is; movement , in which it moves if the energy intake of the current timestep is lower than the predicted costs of maturation and maintenance; if they have not yet reached maturity age, maturation is executed, in which it spends energy in maturation; if maturity age was reached and enough energy is available, reproduction , in which the organism reproduces asexually; then mortality is executed, in which the organism may die due to the direct killing perturbation (during the perturbation phase), or due to old age; and then age ( ageing ). Each time step, the order by which the organisms execute their sequence is ordered according to the organisms' current energy level, with some noise ( sort-organisms ), this way favoring organisms with higher energy levels on average. When an offspring is born, its body size , fecundity , maturity age , and dispersal ability are inherited from the parent, but change due to mutation. Under organism-initialization high level variables are calculated, such as energy budgets considering the agent's body size and system level energetic parameters and the longevity dependent on maturity age. 36

Design concepts
The sections Learning, Objectives, Prediction and Collectives do not apply to this model.

Basic principles
Four core assumptions govern the trade-offs between the focal traits (Table S1.2). set of traits depends highly on the availability of resources and spatial arrangement of resources in the world, as well as on the intensity and type of the perturbation.

Adaptation and Sensing
Organisms perceive their internal energy/reserves level and energy-income and decide to move if their energy-income was lower than the costs of maintenance and growth.

Interaction
Interactions between individuals are indirect, via competition for a limited resource.

Initialization
Initialization begins with generating a map that determines the spatial distribution of

Submodels
To ensure that each submodel was working properly we made use of the Netlogo features such as inspection of agent's values of state and other variables, the use of monitors to see the evolution of the trait values, and other useful information, such as the number of organisms, total amount of resources available, mean age at reproduction, etc, as well as the visualization tool incorporated, that allows the user to inspect each agent's action.
Under each of the submodels we provide some examples of problems we faced while implementing some of the original code.

generate-map
This submodel creates a landscape with dimensions given by the parameters w-width and w-length , and distributes resources to patches according to the parameters nr-resource-clusters and resource-patch-fraction . An initial number of patches, given by nr-resource-clusters , is assigned as resource patches, their max-resources assigned to 1, and their resource-regen initialized as regen-rate . The patches around each resource-patch are then also assigned as resource patches, until the number of resource patches in the environment reaches the total number of patches * resource-patch-fraction . The remaining non-assigned patches are then initialized as bare-patches, with max-resources , 42 resource-regen , and regen-rate as 0. The parameter map-seed allows the creation of landscapes that are replicable by setting its value to a number larger than 0.

Rationale:
As resource-patch-fraction increases, the number of patches with resources increases. The more nr-resource-clusters , the less spatial aggregation between resource patches. If nr-of-resource-clusters is 1 resource patches will tend to aggregate in a single area, whereas when nr-resource-clusters is very high their spatial distribution is completely random and non-correlated ( Fig. S1.1). Similar algorithms are used throughout the literature to simulate landscapes with different degrees of patchiness of resources (Fahrig 1998;Milles et al. 2020).  (Table S1.1).
Rationale: Derived variables are not strictly necessary in the model, in the sense that they are just the product of other variables and parameters. But their calculation once in a simulation makes the code faster.
The maximum amount of energy of organisms ( max-energy ) is initialized as body-size ^ metabolic-allometric-exponent . Energetic balances, representing the maximum intake of energy during a timestep ( max-energy-intake ), maintenance and maturation costs per timestep ( maintenance-cost , maturation-cost ), and costs of dispersal per unit of distance traveled ( dispersal-cost ) are calculated in relation to their max-energy : max-energy-intake = max-energy * max-energy-intake-coef , and so forth (Table S1.1). Thresholds of reproduction, such as the minimum energy to reproduce ( energy-to-reproduce ), and the energy remaining in the organism after reproducing ( energy-after-reprod ), are calculated in the same way, by multiplying energy-to-reproduce-coef and energy-after-reprod-coef to maximum-energy . The longevity of organisms is also calculated as the product of longevity-maturity-coef with maturity-age . 45 update-patches Patches regrow resources . Resources are incremented by resource-regen * max-resources , but not exceeding max-resources (in this case resources become max-resources ).
Rationale: In the real world, resources typically regenerate until reaching a maximum level.
The rate at which resources regenerate may depend greatly on the system being studied.
When resources are organisms, such as plants, their growth is typically modelled with logistic growth curves. Some systems are modelled with adapted logistic growth curves, such as vegetation under predation effects from herbivores. In such a case, the existence of below-ground reserves enhances growth at low aboveground biomass (Mortensen et al. 2018) . Linear growth of resources may be modelled for systems with constant input of nutrients, such as the constant flow of detritus in lotic systems, used by detritivores. We opted for modelling the linear resources as it allows one parameter less than logistic growth, and it is also realistic.

sort-organisms
Each organism draws a number, being sorted in descending order of the outcome. This number is randomly sampled from a normal distribution, with mean energy , and standard deviation energy * stand-dev-to-body-size . The energy intake of an organism ( energy-intake ) takes the minimum value of three possible quantities: max-energy-intake , their storage capacity ( max-energy -energy ) or the amount of resources in the patch they are on ( resources ). The energy-intake is then subtracted from the patch's resource s.
Rationale: Organisms will always try to feed as much as they can, limited by max-resource-intake . However, they might not be able to store all the incoming energy if their energetic income surpasses their max-energy . Additionally, they may only acquire as many resources as the environment allows ( resources ).
If the organism is not an alien, energy-intake is added to its energy .  (Hardin 1960) .

movement
After feeding , the organism moves to another patch if the energy-intake is lower than maintenance-cost plus maturation-cost .
Rationale: If an organism's energy-intake is lower than the requirements for maintenance and maturation it might not be incorporating enough energy for its costs. In the case of adults, maturation-cost accounts for the accumulation of energy for reproduction.
If their energy-intake is not enough to cover those costs , organisms should be compelled to search for resources elsewhere. Otherwise there is no reason for the organism to move.
The organism moves by selecting a random direction ( heading is randomly assigned to an integer between 0 and 359, covering the entire spectrum of directions in a Netlogo model), and a random distance between 0 and disp-ability . If the organism is not an alien , they pay for the cost of movement, which is proportional to the distance dispersed ( energy = energydispersal-cost * dispersal distance).

Rationale: The dispersal algorithm models a simple form of dispersal, which is essentially a
Brownian motion (Brown 1828), or isotropic random walk (Codling Edward A et al. 2008).

Nonetheless, it is the form of choice of many organisms, like ballooning spiders and
wind-dispersed plants.
reproduction 48 If age is larger than maturity-age the organism is considered to be mature, and may reproduce (asexually). A mature organism may only reproduce if its energy exceeds energy-to-reproduce . The energy that is allocated to offspring is the difference between the energy level and energy-after-reprod . Then, each offspring receives an equal fraction of the allocated energy (equal to the allocated-energy divided by fecundity ). Each offspring inherits the traits of the parent with a slight mutation , and its age set to 0. maturation If the organism is a juvenile, and not an alien , its energy is subtracted by maturation-cost .
Rationale: maturation simulates costs of developing and assembling the body structures for organisms to be able to reproduce. In many agent-based models, the attaining of maturation comes with attaining a given size at maturation (Sibly et al. 2013

Direct killing
Check mortality above.
Rationale: Direct killing simulates background mortality rates.

Habitat loss process
A random patch is selected and its resources , resource-regen , and max-resources set to 0.
Then, neighbouring patches also have their state variables set to 0, until the maximum number of patches affected reaches the maximum, which is the total number of patches * habitat_loss / 100.

Rationale:
The growth of affected patches follows the same algorithm as the growth of resource patches around initially marked resource patches. Habitat loss simulates the loss of suitable habitat aggregated in space.
Habitat fragmentation process 55 The resources , resource-regen , and max-resources of a random number of patches is set to 0. The number of random patches selected is equal to total number of patches * ( habitat_fragmentation / 100).
Rationale: Habitat fragmentation simulates the reduction in suitable habitat, when the reduction of habitat happens randomly in the landscape.

Habitat degradation process
The resources , resource-regen , and max-resources of all patches is multiplied by (1 -habitat_degradation ).
Rationale: Habitat degradation simulates the reduction of resource availability.

Invasion process
A random number of organisms creates a copy of themselves ( aliens ), inheriting all their state variables. Aliens do not reproduce (check reproduction ), do not gain or lose energy (check feeding , movement , maturation , and maintenance ), and do not die. The number of aliens added is equal to number of organisms * ( invasives / 100).
Rationale: In the real world, alien invasive species often are successful in relation to the natives because they have no antagonists (predators, diseases, parasites, etc. (Table S1.3). In most cases, extreme parameter values were selected as being three times lower and higher than standard trait values (thus capturing about a magnitude change in trait values). For parameters that are bounded on one side (e.g. energetic costs, resource clusters), we either selected the extreme low and high values for that parameter (maximum resource patch fraction is 1), or we selected a value close to the extreme (maintenance costs cannot be higher than 0.2, otherwise the energy intake of organisms is lower than the energy expenses). In some cases, selecting either an extreme low or extreme high did not make sense as the standard value was already in such extreme low or high value (relative energy intake, energy to reproduce, energy-after-reprod). We varied the parameters locally, i.e., varying only one parameter at a time. In the case of the metabolic allometric exponent, we only tested the typical values found in the literature: 0.66 to 1. indicates positive change in trait value. "-" Indicates negative change. "*" indicates cases in which the signal was inspected visually, as the GLMM did not converge. In orange, cases where the signal is opposite to that of the standard runs. 62