The spatiotemporal profile and adaptation determine the joint effects and interactions of multiple stressors

Background Biodiversity is declining worldwide as ecosystems are increasingly threatened by multiple stressors associated with anthropogenic global change. Stressors frequently co-occur across scales spatially and temporally, resulting in joint effects that are additive or non-additive, i


Background
Worldwide terrestrial, marine and freshwater biodiversity is declining at an accelerating rate [1,2].Many ecosystems and their species are threatened by intensifying stressors associated with anthropogenic global change, including climate and land use change [3].Climate change threatens biodiversity by changes in mean climate parameters (e.g., annual temperatures or precipitation) but also in their variance, hence altering the frequency and intensity of climate extremes (e.g., floods or heatwaves).Land use change, typically through agriculture expansion, threatens biodiversity by habitat degradation and loss of connectivity, but also by agrochemical changes (e.g., pollution or eutrophication) [4].Furthermore, interactions between co-occurring climate and land use change related stressors are common (e.g., land cover affects local climate change effects or land structure influences species' dispersal in response to changing climatic conditions) [5].Nevertheless, the interactions between stressors are largely poorly understood, which hampers the development of appropriate conservation strategies in the face of future changes, despite potentially large consequences for biodiversity [4,5].
Global change related stressors frequently co-occur at different spatial and temporal scales, resulting in joint effects on ecosystems [6].Spatially, stressors may occur locally (e.g., pollution by point sources), regionally (e.g., droughts), or globally (e.g., ocean acidification) [7].Temporally, stressors may have many profiles ranging from discrete events (e.g., single climatic extremes) to continuous (e.g., constant discharge of toxic chemicals) [3,8].Multiple co-occurring stressors can interact in additive effects (joint effect = sum of individual effects), but also in non-additive effects, i.e., antagonisms (joint effect < sum of individual effects) or synergisms (joint effect > sum of individual effects) [4,9].However, research on multiple stressors has been criticized for emphasizing interaction classification but neglecting mechanistic understanding of joint effects [10,11].Although for toxins and environmental stressors for which a full dose-response relationship was available, effects were largely non-additive and could be explained by the stressor addition model [12], a qualitative meta-analysis of 15 meta-analyses of experimental studies with multiple stressors by Côté et al. [9] indeed yielded few consistent results.Partly, this owes to a bias in previous experimental studies toward simplified scenarios resulting in distorted, if not false, assessments of joint effects and interactions when extrapolated to real-world scenarios [6,13,14].
Non-additive effects from multiple co-occurring stressors can result from mechanistic stressor interactions, i.e., one stressor alters organisms' response to a second stressor [9], but also from other factors, like nonlinear stressor-response relationships, e.g., stressor alterations trigger a disproportionate response [12,15], the choice of the null model to predict the joint effect [16], or experimental duration [17].
The prevalence of additive and non-additive interactions is forecasted to change significantly in response to global anthropogenic changes [6,18].Notably, synergisms, considered to produce rapid, disproportionate losses of biodiversity and ecosystem functions [9,19,20], are forecasted to occur more frequently and with higher intensity in the future [21].Nevertheless, the capacity to forecast, and counteract the threats from intensifying stressor regimes requires a robust mechanistic understanding of joint effects, which is currently relatively low [9,22,23].
Process-based simulation models can simulate realistic spatiotemporal stressor profiles and, thereby, significantly contribute to a better mechanistic understanding of joint effects and related interactions across scales [6,24,25].In addition, they support the development of corresponding theory [11,26].Compared to data-based models, which rely on empirical information, often are very specific, and make fewer assumptions, process-based models have a high general predictive power, yet also structural uncertainty and associated risk of bias, e.g., potentially resulting from inappropriate assumptions [27,28].Process-based models support predictions beyond the range of observed stressors and allow for broader extrapolation, as they can simulate alternative developments and parts of real-world systems; they can also be applied to study hypothetical scenarios, whose design can be very complex with a high number of parameters [29] focusing on the identification of general principles [26,28,30].Nevertheless, process-based models in multiple stressor research are rare, although related studies can support future management and thereby optimize species conservation by assessing and predicting joint effects and related interactions [17,28,31].
By contrast, large-scale empirical experiments with multiple stressors are technically difficult to realize [11,32].Therefore, in experiments, multiple stressors have commonly been implemented with static profiles, where stressors are temporally constant or synchronous in pulses, and spatially homogeneous [6,33].However, at the landscape scale and beyond, multiple stressors typically occur discontinuously, i.e., intensities vary in space and time [25,34].Complex spatiotemporal stressor regimes [32,35] and thus associated complex patterns of joint effects and interactions emerge [6,20].For example, a variable spatial intensity of a continuous stressor such as habitat fragmentation together with another, regionally uniform discrete stressor such as drought events created locally dissimilar joint effects on insect populations [4].Yet, biotic exchange processes between patches within ecological networks, e.g., driven by the dispersal of organisms, can in turn compensate for local differences over time [36,37].Changes in the temporal stressor intensity, in the interval between or the frequency of stress events, can also modify the joint effects by e.g.altering the available recovery time [9,38].The majority of ecosystems are exposed to a variety of stressors occurring at different times, hence investigations of stressor dynamics are required to better understand overall impacts [39].For example, studies have shown that the sequence of occurrence of two discrete stressors can strongly modify their joint effects, which may turn from non-significant to non-additive [40,41].
Another shortcoming in multiple stressor research is that adaptation processes of organisms to stressors, likely modifying their effects and interactions, have been largely unconsidered in most previous experimental and modeling studies [35,42,43].Adaptation through phenotypic plasticity [44] or the evolution of life-history traits [43] can enable species to cope with changing stressor regimes.Consequently, future stress of similar intensity may have a lower effect on adapted individuals, populations, and communities, respectively [6,31].Moreover, adaptation can determine the response to additional stressors, for instance, through co-tolerance [45], tradeoffs [46], or maladaptation [47].
We aimed to identify how complex spatiotemporal stressor profiles and adaptation influence the joint effects and thereby interactions of multiple stressors at the landscape scale compared to simplified scenarios.Hereto, we used a spatially explicit, process-based meta-population model [48,49] to simulate hypothetical scenarios of two stressors considering potential adaptation to one stressor.Agricultural land use and climatic events were implemented as continuous and discrete stressors, respectively, as these are key global change stressors likely to have strong future impacts on ecosystems [50,51].Agricultural land use-related stress continuously affected patch qualities and was incorporated with multiple temporal profiles; i.e., steady ramping, and abrupt stepwise (positive or negative).Climate scenarios resulted in numerous, irregularly timed mortality events, with the frequency and intensity of events randomly drawn from three skewed normal probability distributions.Adaptation was considered as different levels of adaptability to climatic events based on previous events.Finally, we determined the response (i.e., joint effect and interaction size and type) of meta-populations to multiple combinations of stressors and adaptation levels.
Following recent conceptual studies [6,20], we hypothesize that scenarios incorporating complex spatiotemporal stressor profiles and adaptation profoundly change joint effects and interactions of multiple stressors compared to simplified scenarios, specifically altering the frequency of additive and non-additive effects.Given the lack of data, we had no specific expectations regarding the directional change of the joint effects concerning the multiple, interacting model parameters, which in isolation most likely either increase (the two stressors) or decrease (adaptation) the meta-population response.

Methods
We used a process-based, spatially explicit meta-population model for a generic hemimetabolous freshwater insect, parameterized based on traits of the European damselfly Coenagrion mercuriale [48,49], to simulate hypothetical scenarios incorporating complex, dynamic spatiotemporal profiles of a continuous (land use) and a discrete stressor (repetitive climatic events of mortality), considering adaptation to the latter.Hemimetabolous freshwater insects complete their larval development in streams but primarily disperse and mate on land [52].They the largest group of aquatic insects.However, metapopulation models for these insects are lacking, given most existing models are largely restricted to terrestrial or freshwater species.This limits our capacity for predicting the effects of changing environmental conditions on hemimetabolous species [49].
In the following, we provide an overview of the model utilized, i.e., the set-up of meta-population networks and the simulation process of population dynamics (i.e., reproduction and dispersal) within these, as well as introduce the stressor scenarios and the implementation of adaptation.The technical details of the model are described in Streib et al. [48], respectively, the softwareframework as well as supplementary data and material used for the present study is provided in a GitHub repository.

Meta-population networks and population dynamics Meta-population networks
Meta-population networks form the basis for the simulation process and were set up on a 12.5 km × 12.5 km landscape raster and consist of quality-assigned patches and interpatch connections.The landscape raster was extracted from real-world geospatial data, i.e., land cover data [53] around a stream network section [54] from South-West Germany.We classified the land cover data into three terrestrial types 'Urban' (LT1), 'Forestry' (LT2), and ' Agriculture' (LT3) and assigned the cells intersecting the stream network as the aquatic land cover type 'Stream' (LT4) (Fig. 1A).Each type was associated with (numeric) dispersal costs (LT1 = 100, LT2 = 75, LT3 = 50, LT 4 = 25), based on a literature review described in detail in Streib et al. [49], representing species-specific landscape permeability for Coenagrion mercuriale.
The meta-population patches were defined by randomly selecting 10% of LT4 cells from the inner 10 km × 10 km landscape area as suitable meta-population patches (Fig. 1B).The (most cost-efficient) interpatch connections (Fig. 1C) were identified based on the patches and the LT-specific dispersal costs in the landscape raster using least-cost path analysis [55], considering only connections below a maximum connectivity cost threshold [49].To minimize bias because of the random patch distributions along the stream network, we created a total of ten meta-population networks for simulation.
To capture the land use scenario, simulating different shifts of ' Agriculture' , starting from an initial 50/50 split of LT3 into intensive (LT3-i) and extensive (LT3-e) use, patch quality was determined time-step specific.Thereto, we considered the terrestrial land cover composition in the complete upstream catchment AP Total and set the specific quality for a patch Q P to its maximum 1.0 in the absence of 'Urban' and intensive ' Agriculture'; else, we linearly reduced the Q P as a function of proportion to a minimum value of 0.0 for 100% LT1 and LT3-i land cover in the catchment of a patch A P :

Population dynamics
Within each meta-population network, we simulated the following population dynamics: (1) Reproduction in each patch up to the maximum carrying capacity K P , and (2) inter-patch dispersal based on positive density-dependent emigration.K P depended on time-step specific patch quality and linearly increased with quality from 0 to 100 (1) individuals (i.e., K P = Q P × 100).If patches were extinct or not fully colonized, i.e., were below their maximum K P due to prior mortality by climatic events, recovery within the meta-population proceeds.Recovery was based on reproduction within a patch and on density-dependent dispersal within the network, i.e., the emigration of individuals into (directly) connected patches.The outcome of dispersal was determined by connection costs: dispersal mortality increased with cost, which represented higher risks and reduced energy reserves, and more dispersers were emitted to patches with lower costs when multiple patches were connected (based on higher colonization probability).

Agricultural land use scenarios
Agricultural land use (for simplicity hereafter referred to as 'land use') represented a continuous stressor.We assumed that intensification or extensification of agricultural land use in a catchment results in reduced or increased patch qualities, respectively, and thereby reduced or increased carrying capacities of patches [4].As we aimed to show general patterns, we implemented four different, hypothetical scenarios of agricultural land use change based on simple temporal profiles of a continuous stressor as suggested in Jackson et al. [6]: steady ramping shifts (Fig. 2), or abrupt stepwise shifts (Fig. 7) to 100% extensive agricultural land use (LT3-e), which decreases land use stress over time, or to 100% intensive agricultural land use (LT3-i), which increases land use stress over time, starting from 50% LT3-e and 50% LT3-i agricultural land use in the landscape raster.
For both profiles, results were highly similar at the end of the simulation.Therefore, we mainly present results for the 'ramping' scenarios, and differences for the 'stepwise' scenarios profile are discussed only briefly (see Appendix A-Figs.8 and 9).

Climate scenarios
We implemented climate scenarios as discrete stressor events of mortality resulting in population declines (i.e., reduction by an integer number of individuals) in a colonized meta-population patch, assuming the same effect is present at the regional scale, uniformly across all metapopulation patches.
We used normal probability distribution functions (PDFs; Figueiredo and Gomes [56]) with different skewness α to sample event mortality M in a range from 0 to 100 individuals for every time step.The PDFs provided probability weights for M sampling based on a logarithmic vector between log(0.01) and log (10) with e as the base and values rounded to integers.Three scenarios were simulated: (1) 'moderate' (α = 0), (2) 'severe' (α = − 2.5), and (3) 'intense' (α = − 5) (Fig. 3).Thus, event mortality increased on average with the scenarios.
To minimize bias, we sampled ten random sequences of climatic events (i.e., mortalities per time-step over the simulation) per climate scenario.

Adaptation
Adaptation comprises evolutionary or non-evolutionary processes and phenotypic plasticity by which species can possibly adapt to changing environmental conditions [6,57].We implemented adaptation in a greatly simplified form as a generic concept, excluding potential negative trade-offs (e.g., maladaptation).Based on a hypothetical approach, climatic event-induced mortality M for a time step t was buffered on the basis of the mean mortality of the last five time steps M t-5 in percentage terms, with a factor A controlling the power of the buffering: Overall, we simulated three levels of A: 'no' adaptation (A = 0), 'low' adaptation (A = 0.5), and 'high' adaptation (A = 1).

Simulation process
The simulation included the population dynamics under different land use and climate scenarios as well as adaptation to the latter (Fig. 4).75 time steps were simulated, starting from a fully colonized meta-population network, i.e., an initial population of 100% of the carrying capacity K of all patches.Land use scenarios and the resulting changes in patch qualities were simulated between time steps 25 and 50.Accordingly, patch qualities remained stable in the first and the last 25 time steps, so a stable initial and final state was achieved.The latter was done to capture long-term land use impacts, i.e., whether joint effects remain stable over time.Climate scenario-induced (2) mortality events were simulated over the entire simulation process and reduced the patch populations depending on the adaptation level.
In total, we ran 4500 simulations based on 10 metapopulation networks, 5 land use scenarios, 10 sequences for each of the 3 climate scenarios, and 3 adaptation levels.Accordingly, 100 runs (10 meta-population networks × 10 climate scenario sequences) were simulated for each combination of land use and climate scenario and adaptation, and the mean value was determined.

Quantification of joint effects
The joint effect E for a scenario combination S (i.e., land use × climatic events × adaptation) was calculated as the mean change in meta-population size relative to a baseline.The baseline B was defined as static land use (i.e., no land use change), excluding the simulation of climatic events and, thus, adaptation.This baseline provided an easily interpretable comparison between the discrete and continuous stressor.In addition, stressor levels were generic and defined to systematically study the effect levels of a discrete and continuous stressor rather than to represent a current or past climate scenario.
The meta-population size corresponded to the sum of individuals in all patches at the end of a simulated time step i.For all 100 simulation runs (i.e., 10 meta-population networks × 10 climatic events sequences) related to a scenario combination, we determined the deviation of the resulting mean population size N S from the mean population size of all baseline simulations N B at all simulated time steps i and calculated E(i) as:

Assessment of stressor interactions
The type and size of a stressor interaction for a scenario combination S were calculated using the multiplicative null model (also termed Response Addition, Effect (3) Addition, or Independent Action), given mortality as the ecological response [9,16].Hereto, we compared the simulated effect E(i) of a land use and climate scenario combination to predicted effects P(i), i.e., the probabilistic sum of their individual effects: where E c and E d of a scenario combination S are the individual effects of the continuous (i.e., simulated for one specific land use scenario without climatic events) and discrete (i.e., simulated for one specific climatic eventscenario and the 'static' land use scenario) stressor alone.
Following the concept of model deviation ratio MDR [58], we defined interactions I(i) by the ratio of predicted to simulated effect [i.e., (1 + P(i))/(1 + E(i))].A ratio of 1.0 corresponds to an additive, a ratio > 1.0 to an antagonistic, and a ratio < 1.0 to a synergistic interaction.Consequently, the magnitude of non-additive (antagonistic or synergistic) interactions increases with departure from 1.0.

Results
The present study intended to identify how complex spatiotemporal profiles of a continuous and a discrete stressor, with potential adaptation to one stressor influence joint effects and interactions on meta-populations.Agricultural land use represented the continuous stressor impacting meta-population patch quality and network connectivity, while climatic events of temporary mortality represented the discrete stressor; adaptation reduced time-specific mortality based on past climatic events.
Excluding adaptation, we found that the discrete stressor (i.e., climatic events) primarily dominated joint effects, whereas the continuous stressor (i.e., land use change) always dominated the interaction type.Adaptation lowered joint effects but changed interaction sizes and, thereby, classification inconsistently across land use and climate scenarios; for decreasing land use stress, adaptation (partly considerably) reduced ( 4) (See figure on next page.)Fig. 4 Flowchart of the stylized simulation process for one meta-population network.START: the initial population size N per meta-population patch P is set to its carrying capacity K P .SIMULATION: population dynamics and dispersal, discrete mortality events based on one climate scenario are simulated over 75 time steps t or until an entire meta-population is extinct (i.e., ΣN t = = 0); between time steps 25 and 49, land use scenario related shifts in agricultural land 'LT3' use are simulated.The event mortality M At of a time step is buffered over the respective adaptation level A and subtracted from a patch population N Pt-1 .Based on the resulting N Pt (logistic) population growth is simulated; in the time steps concerned, the land use scenario modifies the carrying capacity KP via the patch quality Q P , determined as a function of the proportion of 'urban' A P LT1 and intensive 'agricultural' A P LT3-i land use in a catchment A P Total.Dispersal-driven processes are next simulated and result in changes in N Pt based on the difference of immigrants N EMt from directly connected patches and emigrants N IMt ; for clarity, the processes to calculate immigration and emigration are not shown here-these are described in detail in Streib et al. [49].STOP: Store and export the results for subsequent data analysis Fig. 4 (See legend on previous page.)meta-population declines, yet had little effect for increasing land use stress.Details are presented below.

Joint effects
In the absence of adaptation, severe and intense climate stress dominated the joint effects in any related scenario combination, with major declines or extinctions of the meta-population (Fig. 5).Even for decreasing land use stress (i.e., decreasing intensive agriculture), the joint effects were only slightly lower at 0.8 (i.e., 80% meta-population size reduction) than for increasing land use stress (i.e., increasing intensive agriculture) at 0.95.However, for moderate climate stress, the land use scenario dominated the joint effects.
Adaptation had negligible influence on the trajectory of joint effects in the moderate scenario, but significant influence in the severe and intense climate scenarios (Fig. 5).Here, adaptation slowed down meta-population decline for increasing land use stress, yet declines over time were still severe, virtually resulting in extinction.For decreasing land use stress, however, adaptation had a pronounced positive influence; high adaptation even yielded a slight recovery close to the baseline in the moderate climate scenario.
The land use scenarios of the 'stepwise' profile showed qualitatively similar patterns in the long term, resulting in comparable levels of meta-population sizes at the end of the simulation (Fig. 7).

Interactions
In general, we found that non-additive interactions emerged over time triggered by land use change (Fig. 6); however, no and low adaptation resulted in additive interactions for the intense climatic scenario or close to additive interactions for the severe climatic scenario.Antagonistic interactions progressively developed with increasing land use stress and synergistic interactions with decreasing land use stress.For any given stressor combination, adaptation changed the interaction sizes, but changes depended on the climate scenario.The interactions inconsistently decreased in the moderate climate scenario with adaptation levels but increased in the severe and intense.As a result, for severe climate stress low adaptation resulted in a similarly strong antagonism as high adaptation and severe climate stress at decreasing land use stress, but in much lower synergism at increasing land use stress (Fig. 6, I5).By contrast, in low climate scenarios, the interaction size of synergism and antagonism decreased consistently across land use scenarios with adaptation levels.

Discussion
We simulated two stressors (i.e., agricultural land use and climate) with complex spatiotemporal profiles in varying scenarios and evaluated their effect on metapopulations of a freshwater insect.Compared to a static baseline scenario, we found that joint effects and interactions developed fairly differently over time, depending on the stressor levels and dynamics.Adaptation to climatic stress reduced the realized joint effects, but partly increased the strength of interactions, which resulted in a change in stressor classification from additive to nonadditive, supporting previous calls to place less emphasis on stressor classification.
To our knowledge, this is the first modeling study to show that interactions depend on spatiotemporal stressor profiles and adaptation.We conclude that static stressor scenarios over short periods, as often used in experiments, without including potential adaptation, may be insufficient to reliably predict the joint effects and interactions of multiple stressors under real-world conditions.
Below we discuss our results in detail.

Relevance of discrete and continuous stressors under dynamic scenarios in the absence of adaptation
We found that moderate climate stress had no or only negligible long-term impact on meta-populations, yet became dominant for severe and intense stress.For moderate climate stress, land use dominated the joint effects and a relatively strong antagonism developed with increasing land use stress, whereas decreasing land use stress resulted in synergism.Here, the two temporal profiles of the continuous stressor (i.e., 'steady ramping' and 'stepwise') were clearly reflected in the result, however, similar joint effects develop over time.For severe and intense climate stress, the joint effects were largely or completely disjointed from land use (regardless of the simulated profile), thus only minor non-additive ('severe' scenario) or additive ('intense' scenario) interactions emerged.Côté et al. [9] showed that management of a regional continuous stressor to seagrass populations in multiple stressor environments lowered mortality and thereby led to higher populations, in turn resulting in antagonism.Consequently, higher continuous stress likely can have the opposite effect, i.e., synergism via higher mortality and thus lower population sizes.Furthermore, optimizing land use can mitigate long-term impacts from more severe extreme events and thereby strong joint effects, as shown in a review of climate and land use change feedbacks [4].However, positive benefits from land use optimization are only to expect if climatic stress remains below a certain level and sufficient time for recovery is provided, as shown in a study on dry coniferous forests in western North America under rapid climate change and altered disturbance regimes [38].Otherwise, the discrete stressor dominates, resulting in additive interactions as the benefits of managing a continuous stressor become minimal [9].This matches also a recent review on abrupt changes in ecological systems, indicating that in multiple stressor environments, a strong increase of a discrete stressor likely results in additive interactions [19].
Under moderate climate stress, the dominant role of land use for the joint effects is explained by intensive agriculture in the catchment, determining the mean patch quality in the meta-population networks.Positive land use change resulted in antagonism, whereas negative land use change resulted in synergism, as higher patch qualities produced higher meta-populations based on increased carrying capacities, while lower qualities produced lower meta-populations.Comparable to the results of a laboratory experiment by Bible et al. [39] which showed that recovery can neutralize synergistic effects of multiple stressors on the survival of Olympia oysters, we found that extreme climatic events resulted only in temporal population declines.Here, mean event mortality was on average low and, simultaneously, metapopulations had sufficient time to fully recover between events (Bruder et al. [36]).For details on underlying general processes, see Streib et al. [48].Accordingly, the two temporal land use profiles are clearly reflected in the results during the simulation period (i.e., between time steps 25 and 50); however, as 100% of extensive agricultural land use (LT3-e) or intensive agricultural land use (LT3-i), respectively, is reached thereafter, the results are very similar.For severe and intense climate stress recovery after individual events may be incomplete and insufficient to survive subsequent events, given that with more intense climatic extreme events: (1) mean event mortality increased and (2) time to the next intense events decreased; notably, a random event sequence, compared to a synchronous, non-realistic sequence, can result in very short intervals [6].This produced strong joint effects that were predominantly or entirely driven by climate stress, as events became so severe that even improved patch quality with decreasing land use stress failed to prevent the extinction of most or all patches over time (Figs. 10 and 11).Consequently, weak non-additive interactions emerged for severe climate stress and additive interactions emerged for intense climate stress, as metapopulations became extinct or were reduced to almost zero.

Role of adaptation on joint stressor effects and interactions under dynamic scenarios
Adaptation reduced the joint effects for all scenario combinations.Consequently, non-additive interactions decreased with adaptation for the moderate climate scenario, as the dominant impact of land use became even stronger.However, non-additive interactions increased or emerged for the severe and intense scenarios and joint effects were strongly reduced, because under adaptation climate stress no longer dominated; so the temporal land use profile is also reflected in the results during the simulation period (i.e., between time steps 25 and 50).Here, antagonistic interactions based on decreasing land use stress increased more than synergistic interactions based on increasing land use stress, indicating that land use optimization may provide more positive outcomes than to expect.
Our findings support the idea that adaptation needs to be considered in predicting joint effects and interactions of multiple stressors [34,43].As theoretically outlined by Bush et al. [44] or Jackson et al. [6], reduced joint effects by adaptation are expected when ecosystems or species have sufficient time to adapt to single or multiple changed environmental stressors.In the best case, stressor effects are even reduced despite more intense events.Hughes et al. [34] showed empirically that adaptation alleviated the effects of marine heat waves on coral reefs.Adapted corals were significantly less affected by a second heat wave, albeit this was more intense.The effect of joint stressors may be overestimated if not considering adaptation as demonstrated in a modeling study on fruit flies [44].Furthermore, beneficial impacts by management, e.g., land use optimization, may be missed.Our second major finding, that adaptation produces higher interactions in intensified stressor regimes despite strongly reduced joint effects, is initially surprising yet supported by a recent study by Orr et al. [43].Using data from an evolutionary experiment with the rotifer Brachionus calyciflorus, they showed that higher interaction sizes can emerge if adaptation reduces both individual and joint effects compared to control.A change in interactions from antagonism to synergism followed when the reduction in individual effects was greater than the reduction in joint effects.For management, this implies that actions should focus on reducing joint effects based on a mechanistic, predictive understanding of interactions [28], rather than automatically seeking to prevent synergies [16,43].However, adaptation to new stressor regimes is difficult to predict [44,59], making predictions of joint stressor effects and thus interactions challenging [44,59].While adaptation may occur in many species, the speed and strength of adaptation remain unclear [60].Additionally, adaptation to one specific stressor can result in trade-offs in terms of increasing sensitivity to additional stressors [4].For example, pesticide-adapted populations of Gammarus pulex were more sensitive to temperature increases, possibly reflecting the fitness cost of genetic adaptation to pesticides [46].However, species may also respond positively to an increasing stressor if the environment shifts toward their niche optimum, or if competition or predation is reduced by removing more vulnerable species, thereby increasing their robustness to other stressors [45,61].However, local adaptation may prove detrimental if environments change rapidly, i.e., maladaptation, which likely to become more frequent in the future with global change [62,63].For example, maladaptation was demonstrated for the rotifer Brachionus calyciflorus, where adaptation to stressors shifted populations' environmental optimum, resulting in a higher fitness costs when returning to the original environment [47].
The changes observed in joint effects and interactions are explained by the reduced extreme event mortality associated with adaptation.Notably under severe and intense climate stress, meta-populations were more resilient as single patches survived short intervals of high mortality more frequently and in higher numbers when adaptation is present allowing them and, subsequently, the meta-population to recover.Compared to the absence of adaptation, it turns out that the recovery potential is highly dependent on land use, as improved patch quality now prevents more patches from extinction over time with decreasing (Fig. 10) than with increasing land use stress (Fig. 11).Consequently, (both under 'low' and 'high') adaptation produced relatively stronger reductions in joint effects and thereby stronger non-additive interactions for reduced land use stress.

Conclusions
We provide theoretical evidence that scenarios of complex spatiotemporal dynamics and adaptation are critical to understanding how species respond to modified multiple stressor regimes under global change.Albeit the approach is primarily hypothetical we are confident that this work contributes to the mechanistic understanding of how multiple stressors in real-world environments act across space and time by demonstrating and explaining general principles.
Analysis of simplified static regimes at local scales over short periods is likely insufficient for reliable prediction of future joint effects and interactions, notably since recovery processes are not sufficiently considered.Moreover, adaptation likely reduces joint effects and, thereby, can alter interactions with inconsistent direction and size.We expect that these findings could be tested in an experimental setting with moderate effort.
Consequently, regarding management, actions should focus on reducing the strongest individual or joint effects, rather than placing too much emphasis on interactions.Not considering adaptation can result in overestimating joint effects and potentially missing beneficial management outcomes.In A the red line shows the results for N with 'no' adaptation, the orange line with 'low' adaptation, and the green line with 'high' adaptation.In B the red area shows the rate of colonized patches in the range between 0 and 1 with 'no' adaptation, the orange area with 'low' adaptation, and the green area with 'high' adaptation; note that the red area is displayed in front of the orange area, which is displayed in front of the green area Fig. 11 Exemplary representation of A total population size N and B rate of colonized patches (y-axis) over time i (x-axis; i.e., 75 time steps) for one meta-population network simulated for one sequence of the severe climate scenario and the 'ramping' land use scenario with increasing intensive agriculture ↘.Event mortality M (secondary y-axis; i.e., population decline per patch in the range of 0 to 100) per time step is represented by black lines vertically downward from the top.In A the red line shows the results for N with 'no' adaptation, the orange line with 'low' adaptation, and the green line with 'high' adaptation.In B the red area shows the rate of colonized patches in the range between 0 and 1 with 'no' adaptation, the orange area with 'low' adaptation, and the green area with 'high' adaptation; note that the red area is displayed in front of the orange area, which is displayed in front of the green area

Notes on the graphs shown in Figs. 10, 11
Successive mortality events with low to medium mortality in the first time steps result in only small changes in the meta-population size (~ − 80%) and rate of colonized patches (~ − 50%).However, without adaptation, a high mortality event in time step 19 sharply reduces meta-population size and the rate of colonized patches, with levels remaining relatively constant in the following; with adaptation, reduction in meta-population size and the number of colonized patches is less severe, with a partial recovery following (more pronounced for 'high' relative to 'low').The effects of the land use scenarios simulated from time step 25 to time step 50 are only pronounced with adaptation.Here, despite a high mortality event on time step 40, decreasing intensive agriculture (Fig. 10) generally results in progressively increasing meta-population size and a constant colonization rate (both, 'low' and 'high' adaptation), while increasing intensive agriculture (Fig. 11) results in decreasing meta-population size and colonization rate (both, 'low' and 'high' adaptation).Consequently, without adaptation population sizes and colonization rates are comparable for both land use scenarios, but with adaptation they are clearly higher under decreasing land use stress.After time step 50, the long-term effects of the land use scenarios become clear, i.e., high mortality events in time steps 64 and 73 result in metapopulation extinction under increasing land use stress, whereas the meta-populations survive and recover subsequently under decreasing land use stress.

Fig. 1
Fig. 1 Flow of the set-up of meta-population networks-A Landscape raster extracted from real-world geospatial data, with following land cover types: grey = 'Urban' (LT1), green = 'Forestry' (LT2), white = 'Agriculture' (LT3), and blue = 'Stream' (LT4).B Random patch arrangement scenarios along LT4 of the landscape raster, patches are represented as 'black' points.C Meta-population networks determined via least-cost path analysis based on (A) and (B); networks consist of patches, represented by 'black' points, and the most cost-efficient connections, represented as 'black' lines

Fig. 2 Fig. 3
Fig. 2 Temporal profile associated with the land use scenario for the 'ramping' scenario.Here, the red lines show the trend of intensive agriculture (LT3-i), and the green lines show the trend of extensive agriculture (LT3-e)

Fig. 5
Fig. 5 Development of joint stressor effects (y-axis) compared to the baseline B (black line) over the simulation period (x-axis) for all stressor scenarios-adaptation level combinations, with lines corresponding to the mean over 100 simulation runs per combination.Joint effects sizes > 0.0 imply a lower mean population size compared to B, and < 0.0 imply a higher mean population size compared to B. The climate scenario is labeled at the top of the figure, and the land use scenario at the right, where ↗ represents 'ramping' increasing extensive agricultural land use, and ↘ 'ramping' decreasing extensive agricultural land use.The adaptation level is coded by the line color: red is 'no' , orange is 'low' adaptation, and green is 'high' adaptation

Fig. 6
Fig.6 Interaction size (y-axis) over the simulation period (x-axis) for all stressor scenarios-adaptation level combinations, with lines corresponding to the mean over 100 simulation runs per combination.Interaction sizes > 1.0 indicate antagonism, and < 1.0 synergism.The climate scenario is labeled at the top of the figure, and the adaptation level is at the right and color-coded: red is 'no' , orange is 'low' adaptation, and green is 'high' adaptation.For the land use scenario, ↗ represents 'ramping' increasing extensive agriculture (i.e., decreasing land use stress), and ↘ 'ramping' decreasing extensive agriculture (i.e., increasing land use stress)

Fig. 7 Fig. 8 Fig. 9
Fig. 7 Temporal profile associated with the land use scenario for the 'stepwise' scenario.Here, the red lines show the trend of intensive agriculture (LT3-i), and the green lines show the trend of extensive agriculture (LT3-e)