Mechanistic insights into the role of large carnivores for ecosystem structure and functioning

findings also support previous studies suggesting that natural ecosystem dynamics have been severely modified and are still changing as a result of the widespread decline and extinction of large carnivores.


Introduction
Large terrestrial carnivores have experienced severe population declines and geographic range contractions since the Late Pleistocene (Beschta and Ripple 2009, Prugh et al. 2009, Dirzo et al. 2014, Ripple et al. 2014, Galetti et al. 2018. It is widely accepted that most of these declines have been driven by the rapid expansion of the human population followed by the loss or degradation of their habitat and depletion of their preferred prey (Larsen and Ripple 2003, Tittensor et al. 2014, Sandom et al. 2018, Enquist et al. 2020. Furthermore, historically large carnivores have been widely persecuted because they were perceived as a threat to both humans and livestock (Ripple et al. 2014). At present, large-bodied fauna is declining in abundance faster than smaller body sized organisms (Olff et al. 2002, Dirzo et al. 2014, and only 22 species of terrestrial large carnivores persist today (Wilman et al. 2014). These drastic declines have raised conservation concerns and trophic rewilding has been recently proposed as a method to restore the natural (top-down) control mechanisms and the associated trophic cascades (Svenning et al. 2016, Wolf andRipple 2018).
Large carnivores take the role of top predators in food webs, potentially exerting strong regulating feedbacks on all the lower trophic levels (Elmhagen et al. 2010, Estes et al. 2011, Ripple et al. 2014, Brose et al. 2019. Carnivores can exert control on their prey through predation as well as by non-lethal effects. The presence of top predators can force prey to avoid areas with high predation risk, altering the behavior of species further down the food chain (Le Roux et al. 2018, Pringle 2018, Gaynor et al. 2019, Sheriff et al. 2020. By changing herbivore densities and behavior, they can play an essential role in maintaining a high global plant biomass (i.e. 'green world hypothesis') (Bond 2005, Wilkinson andSherratt 2016). Furthermore, large carnivores can control the abundance of mesopredators through direct competition as well as predation (Elmhagen et al. 2010, Prugh andSivy 2020). In fact, the loss of large carnivores has been linked to rapid population growth of mesopredators (generally mammals), i.e. 'mesopredator release hypothesis' (Terborgh andWinter 1980, Soulé et al. 1988), which in several cases have caused catastrophic impacts on the diversity of the remaining animal and plant community (Brashares et al. 2010, Terborgh 2015. Ecosystems can respond in several ways to the loss of large carnivores, making the possible consequences difficult to anticipate (Terborgh 2015, Ford et al. 2017. For example, the strength of a possible response after the extinction of large carnivores can depend on the trophic structure and environmental conditions present in the pre-disturbed ecosystem (Pace et al. 1999, Terborgh 2015, Ford et al. 2017). These properties arise from the combination of the community composition, primary productivity, seasonality and other environmental factors (Hopcraft et al. 2010, Ripple andVan Valkenburgh 2010). The dynamics between trophic levels and interactions with local climate create unique equilibria between bottom-up or top-down forces at different levels in the food web, thereby shaping diversity, structure and functioning of entire ecosystems (Hopcraft et al. 2010, Doughty et al. 2016b, Galetti et al. 2018, Pringle 2018. The loss of regulating feedbacks exerted by megafauna species can, in many cases, be directly linked to drastic changes in ecosystem structure and functioning in combination with a loss in biodiversity (Miller et al. 2001, Prugh et al. 2009, Ripple et al. 2015.
Herbivore biomass is expected to increase with increased primary productivity (Pettorelli et al. 2009). According to the 'exploitation ecosystem hypothesis' (Oksanen et al. 1981, Letnic andCrowther 2013) this positive relationship can be weakened by the top-down effect of large carnivores that hold herbivore populations at low densities (Ripple and Beschta 2012, Letnic and Ripple 2017, Le Roux et al. 2019. In these predator-controlled (top-down regulated) ecosystems, the removal of large predators can lead to a growth of herbivore populations and a shift towards a system limited by resources (bottom-up controlled) (Owen- Smith 1987, Peterson et al. 2003, Ripple and Van Valkenburgh 2010, Letnic and Ripple 2017. According to exploitation ecosystem hypothesis, we would expect herbivores to respond more positively to the loss of large carnivores in high productive areas. These shifts are likely to have many additional cascading effects, changing both the structure and dynamics within lower trophic levels and can result in significant biodiversity losses and a turnover in species composition (Hebblewhite et al. 2005, Beschta andRipple 2009). Multiple studies concluded that in some ecosystems the biomass of large ungulates is only linked to primary productivity and experiences little to no effect from top to down regulation (Owen -Smith 1990, Fritz and Duncan 1994, Mduma et al. 1999. In these bottom-up regulation ecosystems predation may be concentrated on the 'doomed surplus' (Errington 1956), i.e. weak and old individuals, therefore exerting a form of compensatory mortality (Fritz et al. 2011, Ford et al. 2017. In these ecosystems, the extinction of large carnivores may have limited influence and cascading effects on lower trophic levels. Top-down control might also be weak in ecosystems with low bottom-up control because prey species can escape predation. Migratory ungulates, for example, are known to be little controlled by predation (Hebblewhite and Merrill 2009). Similarly, megaherbivores (> 1000 kg) are often too difficult or dangerous to subdue, therefore large carnivores are generally able to exert top-down control on (the smaller) large herbivores only (Estes et al. 2011, Fritz et al. 2011, Malhi et al. 2016, Le Roux et al. 2019. While a large body of evidence support these bottomup and top-down mechanisms, their relative importance is contingent on the specific conditions of the ecosystems studied. Therefore, general inference is largely constrained by the limited temporal and spatial scales at which studies are conducted (Beschta and Ripple 2009, Brashares et al. 2010, Navarro and Pereira 2017, thus limiting our predictive capacity under global change and for wildlife management. Generalizations of the top-down influence of large carnivores on ecosystems are further hampered by two factors. First, vegetation structures have been altered by extensive land use changes (Klein Goldewijk et al. 2017). Second, both large carnivores and large herbivores have, for the majority, declined in distribution and abundance simultaneously (Malhi et al. 2016). Quantifying and disentangling the effect of changes in the top-down regulating processes from the other simultaneous changes using empirical data sources and statistical approaches is therefore extremely challenging.
An alternative approach to investigate the top-down regulating processes of large carnivores for both the structure and functioning of ecosystems is through the use of general mechanistic ecosystem models (GEMs). By modelling fundamental ecosystem dynamics, we can estimate the role of large carnivores by focusing on the emergent properties of the ecosystem when carnivores are removed (Harfoot et al. 2014). This overcomes spatial and temporal scale limitations typical of empirical studies, while enabling full control over ecosystem conditions. Further, this approach can provide insights into many ecosystem components and processes that are generally hard to assess in experimental or observational studies. These methods obviously come with different limitations (as extensively discussed below), nonetheless, they may be key in underpinning ecological theory by providing additional and independent lines of evidence.
Here, we analyzed the influence of top-down control exerted by large carnivores on ecosystem structure by performing mechanistic ecosystem simulations using the Madingley model. The Madingley model is a mechanistic general ecosystem model that represents all (global) photoautotrophic and heterotrophic life with body mass larger than 10 mg (Harfoot et al. 2014, Newbold et al. 2018, Enquist et al. 2020. In our simulations, we focus specifically on the ecological role of extant large carnivores (> 21 kg) (Malhi et al. 2016), therefore estimating what would happen if they would further decline and go extinct. We only consider the direct role of predators in changing herbivore densities, possible behavioral effects, e.g. herbivores avoiding areas with high predation risk (Le Roux et al. 2018), are not incorporated in the simulation approach.
To examine how the removal of large carnivores (i.e. elimination of top-down forces) can influence the distribution of biomass across functional groups, we ran a global-scale simulation, and assessed how changes in biomass relate to gradients of productivity and seasonality. Additionally, to investigate how removal of large carnivores may alter food-web structure and ecosystem function, we simulated four local-scale removal experiments in ecosystems characterized by different environmental conditions, focusing on four illustrative examples: the Hustain Nuruu (Mongolia), Yosemite (United States), Białowieża (Poland) and the Serengeti (Kenya and Tanzania). Finally, we assessed whether the emergent properties of the simulation experiments are in line with the predictions of the 'green world hypothesis' (Bond 2005, Wilkinson andSherratt 2016), the 'exploitation ecosystem hypothesis' (Oksanen et al. 1981, Letnic andCrowther 2013) and the 'mesopredator release hypothesis' (Terborgh andWinter 1980, Soulé et al. 1988).

The Madingley model
The Madingley model simulates both photo-autotrophic and heterotrophic life by categorizing organisms into different categories with distinctive functional roles (Harfoot et al. 2014, Bartlett et al. 2016. Because of the vast number of organisms, especially on a global scale at a spatial resolution of 1 degree, the Madingley model groups organisms in cohorts, reducing the computational requirements (Harfoot et al. 2014). Each cohort describes a group of organisms with identical categorical and quantitative traits (e.g. same body mass, age and functional properties), all of which occur within the same spatial grid cell (Harfoot et al. 2014). In the terrestrial realm, heterotroph cohorts are divided into functional groups by: 1) trophic level (autotrophs, herbivores, omnivores and carnivores), 2) metabolic characteristics (endotherms and ectotherms) and 3) reproductive strategy (e.g. semelparity and iteroparity). In addition, cohorts are also defined by their juvenile, adult and current body mass. The current body mass characterizes most of their behaviors such as feeding preference, growth rate, dispersal capabilities, reproduction and metabolic needs (Pletcher 1999, Brown et al. 2004). Dynamics of autotrophs are simulated as a function of climate with a terrestrial carbon model (Smith et al. 2013). Environmental conditions are set per model grid cell using data supplied by the user, including from Behrenfeld and Falkowski (1997) and ISRIC-WISE (2012).
Simulations at the global scale led to emergent properties at the community and macro-ecological scales that closely resemble the distributions of plant and animal communities globally (Harfoot et al. 2014). The modelled processes are consistent across cells and change as a function of the cohorts present and the climate conditions. In other words, the difference in the emergent properties of different cells stem from the differences in environmental conditions, differences in the initialization parameters (e.g. the range in body weight of initialized cohorts) and stochasticity in ecological dynamics. Equations of the model and a more in-depth description of the various mechanisms can be found in (Harfoot et al. 2014).

Model improvements to study top-down control
In order to utilize the Madingley model for our purposes, we implemented a number of improvements to the original source code of the Madingley model first presented by Harfoot et al. (2014). We used a C++ translation of the code base, because this enabled demanding simulations to be run on a Linux cluster (Cartesius HPC at Surfsara), which was previously not possible with the C# version of the model.

4
The main improvements were: 1) the predation preference for all terrestrial, endothermic, carnivorous cohorts was updated, based on the allometric relationships described in Carbone et al. (1999). This changed the preferred predatorprey body size ratio in the Madingley model from 0.1 for all predators to 0.1 for predators smaller than 21 kg and 1.0 for predators larger or equal to 21 kg. An example how this alteration impacts predation among large mammals is provided in Supplementary material Appendix 1. It is important to note that probability of being preyed upon depends on the current body mass of the prey, which changes when cohorts are growing from juvenile to adult. 2) The Madingley initialization method was updated, allowing the user to set the maximum body mass of modelled cohorts spatially. 3) The model functionality of the Madingley model was expanded to allow the model state to be exported and used in the initialization of a new simulation, thereby extending the previous simulation run.
A full detailed description of all modifications implemented on the source code is presented Supplementary material Appendix 1 Table A1, Fig. A1. The updated C++ source code for the Madingley model as well as the climate input data can be found at: <https://github.com/SHoeks/ MadingleyCPP-openMP>.

Model initialization
The Madingley model was first initialized with cohorts, using the functional group properties described in Table 1. We initialized the simulation, setting a maximum body mass per functional group per cell that reflects the observed species distribution globally. For endotherms the spatial data was obtained by combining body mass information extracted from the EltonTraits 1.0 database (Wilman et al. 2014) and the species ranges from the IUCN Red List (2017) (Supplementary material Appendix 1). For ectotherms, we used data for reptiles and combined body mass data from the amniote life-history database (Myhrvold et al. 2015) with the species ranges from the GARD initiative (Roll et al. 2017) (Supplementary material Appendix 1). The categorical traits are those defined by the functional groups described in Table 1. Juvenile and adult mass are drawn in a probabilistic manner from the minimum to maximum mass range defined for the functional group and grid cell. After the initialization phase the model was run without any perturbations or manipulations for 500 yr to reach a dynamic equilibrium in the overall biomass (Harfoot et al. 2014, Newbold et al. 2018. The model state at the end of each (500 yr) spin-up run was exported and used as a baseline in the further topdown control scenario analyses. There was no anthropogenic land-use change implemented in the simulations.

Large-carnivore removal scenarios
To understand the influence of top-down control we compared a simulation in which all large (> 21 kg) endothermic carnivores were removed to a control (no removal) simulation. By setting the removal threshold to 21 kg, we excluded both large and megacarnivores as defined by Malhi et al. (2016) from the simulation. After the spin-up simulation of 500 yr, large carnivores were gradually removed over a period of 20 yr by starting from the heaviest cohorts (selected by adult body mass) until all cohorts of the targeted functional group were removed. During the removal phase, the removal threshold was consecutively lowered, starting at a threshold based on the heaviest cohort and ending at the desired threshold of 21 kg, making sure equal proportions of the summed biomass were removed at each model time step. The gradual removal procedure was not intended to mimic a realistic extinction scenario, but to minimize the shock induced on the running simulations. When removing large carnivores all at once, the ecosystem tends to fluctuate strongly, resulting in direct extinctions of various functional groups or specific body mass categories. After the removal phase, we run the simulation for an additional 500 years to ensure model conditions had stabilized. The same procedure was applied for the control simulation, except no cohorts were actively removed. We averaged the results across five simulation replicas, which may differ because of variation in initial conditions and stochasticity in ecological interactions. The simulations were done on a local and a global scale which is further explained below.

Global simulations
The Madingley model was run on a global scale to assess the resulting spatial patterns after removing large carnivores globally. The monthly outputs were summarized as yearly averages, in order to account for seasonal variations and limit the output size (70+ GB per replica). After the simulation, the exported yearly averages were summarized across a 20-yr period and mapped spatially. In order to provide insights into the removal of large carnivores and the responses of different ecosystems over a productivity gradient, we also evaluated functional group-specific biomass changes against the yearly average primary productivity extracted from individual grid cells. We further assessed the relative influence of NPP and its seasonality on the relative change in biomass of different categories of cohorts using a random forest algorithm (Liaw and Wiener 2002). Specifically, we considered changes in biomass of 1) autotrophs, 2) all endothermic herbivores, 3) large (> 100 kg) endothermic herbivores, 4) small (< 100 kg) endothermic herbivores, 5) endothermic omnivores and 6) medium-sized (10-21 kg) endothermic carnivores. The random forest models were ran using the relative difference in the biomass of a given functional group as response variable, and NPP and its seasonality as predictor variables. Each random forest model was ran using the 'randomForest' package in R with 1000 trees (see further information in Supplementary material Appendix 2). Contrasting to the food-web simulations, the global simulations were run with a maximum of 500 cohorts per spatial grid cell, in order to lower the total models' run time. These settings resulted in an average model runtime (per replica) of about 5 d on a cluster with two Xeon E5-2660 processors (total of 32 cores).

Small-scale simulations
We also conducted four small-scale simulation experiments to understand direct and indirect changes in food-web structure following the large endothermic carnivore removal scenarios under environmental conditions characteristic of four different locations and ecosystem types. The four locations vary widely in their net primary productivity and seasonality (Table 2). The first two locations (Hustain Nuruu and Yosemite) were selected for their low yearly average net primary productivity (NPP) (< 1 gC m −2 d −1 ), while the third (Białowieża) and fourth (Serengeti) were selected for their much higher yearly average NPP (> 1 gC m −2 d −1 ). Per NPP category, one location was characterized by relatively high seasonality (Hustain Nuruu and Białowieża) and one location was characterized by relatively low seasonality (Yosemite and Serengeti). We emphasize that these four locations were simulated with the local environmental characteristics (Table 2) and body mass ranges of local communities ( Table 2, Supplementary material Appendix 1 Fig. A1), disregarding ecosystem specific peculiarities. A more detailed description of procedure used for selecting the locations, can be found in Supplementary material Appendix 3. Each small-scale simulated consisted of an isolated area of 5 × 5 1-degree grid cells around a central coordinate ( Table 2). The biomass time series exported during the 500-yr spin-up simulation phase are presented in Supplementary material Appendix 4. We continued to export the basic time series results (biomass and abundance per functional group) for the 500-yr simulation phase post removal. In addition to the time series, we also exported the individual cohort properties, which were needed to create the food web. In order to simplify the food-web figures as well as to limit the amount of data exported to a maximum of 1 GB per modelled time step, food-web node biomass as well as the biomass fluxes between nodes were summarized on three levels: 1) endothermic or ectothermic, 2) feeding guild (carnivores, omnivores or herbivores) and 3) body mass bin (on a log 10 -scale). Moreover, food-web related results were solely exported at the very end of the simulation period (last 20 yr). We analyzed the influence of the removal scenarios on the level of 1) changes in total biomass of all cohorts within the same feeding guild, 2) changes on the total autotroph biomass, 3) changes in the biomasses between food-web nodes and 4) changes in biomass fluxes between food-web nodes. Foodweb simulations were run with a maximum of 1000 cohorts per spatial grid cell.

Spatial consequences of large-carnivore removal
Following large-carnivore removal, the overall carnivore biomass increased in the tropical regions (Fig. 1a) and decreased in the temperate regions (Fig. 1a). The increase in overall carnivore biomass in tropical regions resulted from overcompensation in abundance of smaller to medium-sized carnivores. Herbivore biomass increased at all locations, but the increase was stronger in tropical areas (Fig. 1b). Omnivores showed a similar pattern to carnivores (Fig. 1c), with strong increases in tropical regions and decreases mostly confined to the temperate regions. Due to herbivore biomass increase, the autotroph biomass decreased globally, especially in tropical regions (Fig. 1d).

Top-down and bottom-up shifts
In response to the removal of large carnivores, the positive relationships between NPP and the biomass of carnivores, omnivores and herbivores, increased in steepness (Fig. 2a-c). This effect was particularly pronounced in herbivores ( Fig. 2b), in line with the prediction of the 'exploitation ecosystem hypothesis'. In contrast, the relationship between NPP and autotroph biomass decreased in steepness after the removal of large carnivores (Fig. 2d), which is consistent with the 'green world hypothesis' that predicts a decrease in autotroph biomass when herbivores are no longer topdown controlled. Furthermore, the removal of large carnivores increased the spread of data points around the fitted line (i.e. the variation between the biomasses extracted from the individual spatial grid cells) suggesting that the transition from top-down to bottom-up regulation reduces ecosystem stability. The coefficient of variation calculated over all grid cells increased from 0.56 to 1.12 in carnivores, from 0.36 to 0.66 in herbivores, and 1.17 to 1.45 in omnivores. Both NPP and its seasonality were important drivers for the geographic pattern in relative change in biomass following large carnivore removal, with NPP being the most important variable in explaining global differences. The relative increase in the biomass of different functional groups was stronger at higher levels of NPP (Supplementary material Appendix 2 Fig. A2), but weaker in more seasonal environments (Supplementary material Appendix 2 Fig A3). The effects were consistent across all functional groups, suggesting that large carnivore removal may have stronger consequences in highly productive and low seasonal environments, and milder effects in low productive and high seasonal environments. A more detailed description of the relative effect of NPP and seasonality on different functional groups is provided in Supplementary material Appendix 4.

Temporal trends in biomass
The relative carnivore biomass trends fluctuated right after the start of the large-carnivore removal phase, after these initial fluctuations, the biomass of carnivores stabilized at lower proportions of the control scenario (Fig. 3a). Herbivore biomass increased for all selected locations (Fig. 3b), with Serengeti showing the stronger increase (5.07 times the biomass observed in the control scenario) due to the relatively high NPP (>1 gC m −2 d −1 ) and low seasonality (Table 2). Nevertheless, the herbivore biomass in other locations also showed large increases, stabilizing at 1.91 (Białowieża), 1.84 (Yosemite) to 1.46 (Hustain Nuruu) times the biomass observed in the control scenario. In general, omnivores benefited from the loss of large carnivores, with the exception of the Hustain Nuruu, which is characterized by very low NPP (0.23 gC m −2 d −1 ) and high seasonality when compared to the other locations ( Table 2). The average yearly autotroph biomass decreased the most in the Serengeti (Fig. 3d; 0.43 times), due to the strong increase in overall herbivore biomass (5.07 times). In the other three locations the decrease in average yearly autotroph biomass was less pronounced (0.71-0.93 times).

Alteration of the food-web structure
Large endothermic herbivores with a body mass between 100 kg and 1000 kg increased in biomass after the removal of large carnivores in all four locations (Fig. 4). Megaherbivores (≥ 1000 kg) showed a similar but weaker response as probability of predation events were rare above a certain size (Fig. 4d). We found that the increase in large herbivores was more prominent in the Serengeti than in lower productivity (Hustain Nuruu and Yosemite) or higher seasonality areas (Białowieża). In all four locations, the removal of large carnivores was followed by an increase in medium-sized (10-21 kg) endothermic carnivores. In the Białowieża and Yosemite locations, the biomass of large (≥ 100 kg) endothermic omnivores decreased to 0.58 and 0.36 times the biomass observed in the control simulation, respectively. This decrease stemmed from the increased competition with medium-sized carnivores and endothermic herbivores.

Discussion
Our results support the idea that large carnivores play key roles in ecosystem structure and functioning (Ripple et al. 2014, Wolf andRipple 2018). First of all, the results of our global simulations showed that the presence of large carnivores suppresses the positive relationship between NPP and herbivore biomass density (Fig. 2b), indicating that ecosystems become resource limited after the removal of large carnivores (Letnic and Ripple 2017). This is consistent with the predictions of the 'exploitation ecosystem hypothesis' (Oksanen et al. 1981) and supports previous correlative studies testing this hypothesis Beschta 2012, Letnic andRipple 2017). Second, our model simulation suggests that the complete removal of large carnivores can trigger trophic cascades that alter the entire structure and functioning of ecosystems (Fig. 4). In line with the 'green world hypothesis' (Paine 1980, Beschta and Ripple 2009, Sherratt and Wilkinson 2009), our findings imply that the decline of large carnivores can result in a sharp increase in herbivore biomass, followed by a strong decrease in autotroph biomass ( Fig. 1-3). This supports the idea that the control exerted by large carnivores on herbivore populations can promote standing stocks of plant biomass. However, we must consider that the cascading effect of carnivores on vegetation biomass following the removal of large carnivores may be overestimated by the model. For example, various species-specific plant defense mechanisms (e.g. secondary plant compounds) can limit herbivore populations from consuming all plant biomass (Ford et al. 2014).
It is also possible that the impact of top-down control may have been overestimated by the Madingley model as the population carrying capacity is determined through the  Table 2. limitation and competition for resources (prey consumption) only. This ignores that competition for space (territories) may limit maximum predator densities to lower levels.
It is important to note that our simulations are only able to mimic the role of predators in changing herbivore densities, we did not consider behavioral influences, such as the avoidance of areas in which the risk of predation is large (Le Roux et al. 2018). Furthermore, although the updated predator-prey ratios resulted in more realistic predator prey interactions (Supplementary material Appendix 1), the Madingley model still simplifies predator-prey interactions and disregards mechanisms to escape predation other than through body size. For example, migration can be considered as an important escape mechanism and is highly relevant under certain environmental conditions (Hebblewhite and Merrill 2009). Another supplementary factor to consider in top-down dynamics concerns disease dynamics and its role in controlling regime shifts (Getz 2009). Diseases in wildlife populations can be responsible for keeping large herbivore populations at low numbers (Mcnaughton 1992, Holdo et al. 2009), thereby fulfilling a complementary role to large predators. Our results showed that the stability in biomass across grid cells decreased in the global simulations after the removal of   (Table 2). (a) Hustai Nuruu (low NPP and high seasonality), (b) Yosemite (low NPP and low seasonality), (c) Białowieża (high NPP and high seasonality) and (d) Serengeti (high NPP and low seasonality). Node size depicts relative yearly average standing stocks (complete biomass) of each of the log 10 -binned cohort groups. Note that the node representing endothermic carnivores with a body mass of 10-100 kg is split into two sub-nodes, the left sub-node depicts the biomass of endothermic carnivores with a body mass of 10-21 kg and the right node represents the endothermic carnivores with a body mass of 21-100 kg, matching the threshold used for the removal experiment. The color and width of the node border represent the strength of increase (blue) and the strength of decrease (red) in node biomass as a consequence of removing large carnivores. Numbers within nodes display the relative change in biomass calculated as a ratio (removal scenario/ control scenario). Lines connecting the nodes represent the biomass flows. Red lines indicate decreases and blue lines show increases in biomass flows. Line width provides an indication of the relative importance of the biomass flowing between the connected nodes, comparing all flows corresponding to a single consuming node. Dotted grey lines represent flows that were lost after the removal of large carnivores.
large carnivores. This is also in line with previous studies (Miller et al. 2001, Ripple et al. 2014, Van Valkenburgh et al. 2016, demonstrating that top-down regulating forces promote long-term ecosystem stability through maintaining the trophic structure. The widespread decline of large carnivores and the demise of megacarnivores (Ripple et al. 2014) may therefore have lowered the stability of ecosystems globally. In addition to the impacts concerning herbivore and vegetation biomass, we were also able to identify secondary ecosystem effects. The small-scale simulation results show that the removal of large carnivores triggered an increase in the number and total biomass of medium-sized (10-21 kg) carnivores (Fig. 4), as expected under the 'mesopredator release hypothesis' (Brashares et al. 2010, Elmhagen et al. 2010, Ripple et al. 2014). This is a result of the loss in predation pressure previously exerted by large carnivores and the loss of competition between medium-sized (10-21 kg) and large carnivores. Interestingly, the global results indicate that the strength of the mesopredator release effect increased with NPP and decreased with seasonality (Supplementary material Appendix 2 Fig. A2c, f ). In contrast to the strong response of mesopredators, large omnivores (100-1000 kg) decreases in total biomass after the removal of large carnivores (Fig. 4b-c). The decrease in large omnivores as a result of large carnivores removal was not expected and most likely reflects the simplified feeding model incorporated in Madingley, which assumes that 1) medium-sized carnivores and medium to large omnivores can directly compete for the same prey and 2) all herbivores and omnivores rely on similar herbivore feeding mechanics, thereby competing for the same resources.
Our global simulations also provide novel insights into the context-dependency of ecosystem responses to large carnivore removal. In highly seasonal environments, the effects of large-carnivore removal were less pronounced compared to environments with low seasonality (Supplementary material Appendix 2 Fig. A2, A3). The minimum yearly NPP creates a bottleneck for herbivore populations (Norrdahl et al. 2002), thereby attenuating the increase in herbivores after the removal of large carnivores. It follows that ecosystem consequences of large and megacarnivore decline are expected to be more severe in tropical non-seasonal environments that, historically, have suffered less anthropogenic impacts overall than temperate areas (Rapacciuolo et al. 2017, Santini et al. 2017. The context-dependence of the ecosystem response can also have implications for trophic rewilding actions (Olff et al. 2002, Svenning et al. 2016). The effect of large carnivore reintroduction and recolonizations (Chapron et al. 2014) in ecosystems is expected to vary across biomes, with stronger effects in tropical non-seasonal environments than temperate and seasonal environments.
General ecosystem models, like the Madingley model, are particularly valuable to understand complex ecosystem dynamics -such as the ecological role of large carnivores in ecosystems -overcoming spatial and temporal limitations of empirical studies (Newbold et al. 2018, Enquist et al. 2020. Our results highlight the importance of top-down control exerted by large carnivores on the structure, functions and stability of ecosystems, and open up several possible lines of future research. Our results can have implications for setting baselines and expectations for trophic rewilding approaches. However, future research including field and modeling studies are needed to explore the effects of the reintroductions of previously extinct functional groups in ecosystems (Svenning et al. 2016). Our findings also indicate that the control exerted by large carnivores on herbivore populations globally promotes standing stocks of plant biomass, therefore serving as a possible climate mitigation tool by reducing the amount of atmospheric carbon (Wilmers et al. 2012, Ripple et al. 2014, Enquist et al. 2020, although this generally is not recommended due to the negative effects on biodiversity (Veldman et al. 2015, Bond et al. 2019).
However, our study focuses on the role of large carnivores present in current ecosystems without considering anthropogenic pressures. The effect of the absence of large predators may have been compensated by humans' over-exploitation on large herbivores (Darimont et al. 2015), thereby substituting the role of large carnivores in the ecosystem -although often this effect has been, and still is, far beyond, leading to strong population reductions or losses of the herbivore species (Ripple et al. 2015, Doughty et al. 2016a, Malhi et al. 2016). At the same time, humans may to some extent have replaced the ecological functioning of megaherbivores through agricultural practices (Bocherens 2018), albeit again the effect of humans has overcompensated these processes with detrimental effects (Ellis et al. 2013, Ellis 2015. Finally, our results support previous studies suggesting that the widespread decline and extirpation of large carnivores globally may have substantially altered the global ecosystem structure, patterns and dynamics (Faurby andSvenning 2015, Santini et al. 2017). Future research, potentially based on method and principles shown here, may specifically aim to estimate how global ecosystems have plausibly changed as result of the end-Pleistocene and Holocene megafauna extinctions.

Data availability statement
The C++ source code utilized for our simulation experiments is openly available from Github: <https://github. com/SHoeks/MadingleyCPP-openMP>. The summarized simulation outputs are accessible on a Figshare repository: <https://figshare.com/articles/dataset/ Madingley_Outputs_Global_rasters/12807068>. considers this work a contribution to his Carlsberg Foundation Semper Ardens project MegaPast2Future (grant CF16-0005) and to his VILLUM Investigator project 'Biodiversity Dynamics in a Changing World' funded by VILLUM FONDEN (grant 16549). MH acknowledges the support of the KR and Hempel Foundations through the project 'Towards a brighter future for humanity and biodiversity' to which this research is a contribution. Author contributions -LS conceived the project. All (co-)authors contributed to the research design. SH performed the simulation experiment and interpreted the results. SH wrote the manuscript with assistance from all co-authors. Conflicts of interest -The authors declare no conflicts of interest.