Biodiversity as a Firewall to Engineered Microbiomes for Restoration and Conservation

1. Introduction

Ecosystems are resilient entities capable of dealing with multiple sources of change, from environmental fluctuations to the undesirable outcomes of invader species. However, as with any other complex system involving positive feedbacks, they also involve breakpoints associated to changes in key parameters beyond given thresholds (Scheffer et al. 2001, Scheffer 2009, Lenton 2011, Barnovsky et al. 2012) . The increasing stressors associated with Anthropogenic drivers, from global warming to the increasing use of land associated with a growing human population. These factors already affect extant populations globally, regarding biodiversity declines and defaunation. Because of the accelerated pace of change, rapid transitions are expected to occur in the following decades (Sole and Levin 2022). Such changes will profoundly affect ecosystem services on an unprecedented scale.

Different strategies have been proposed to protect ecological communities across scales, from conservation to restoration. A commonality in all these cases is biodiversity conservation: biodiversity is a firewall to invasion and an indicator of ecosystem health. However, the time window for response is shrinking rapidly, and other paths might need to be taken. These include the engineering of ecosystems that target their fragilities (Solé, 2015; Solé et al., 2015; de Lorenzo, 2016; Piaggio et al., 2017). Such strategies include different microbial engineering strategies, from probiotics to the design of simple genetic circuits to perform specific functionalities. Can ecosystem bioengineering provide the right restoration effect while maintaining (or enhancing) biodiversity? What is the difference between using probiotic strains versus genetically modified cells? To answer these questions, dedicated efforts will be necessary to test them under controlled conditions to evaluate their effectiveness and safety. At this stage, models (both mathematical and computational) can provide great insight into potential implementations that fulfill our previous requirements. The aim of this paper is to provide such insight using a population dynamics approach.

Two potential examples for such an engineering approach are shown in Figure 1a–d. These involve drylands (a,b) and corals (c,d). Both have been extensively analysed and in both cases major declines (or even collapse) are likely to occur over the next decades (Solé and Levin 2022). In both cases too, increasing temperatures trigger community-level effects that jeopardize the stability of these communities. In both cases, the microbial compartment plays a crucial role. Coral reef shrinking is one of the most obvious manifestations of the effects of increasing temperatures, along with overfishing and eutrophication, but also due to pathogens. As a result, rapid declines have taken place in the last decades. What kind of interventions can be used? Within the context of coral reefs, they are known to be holobionts, including a microbial community (Figure 1c) involving a small number of taxa. The health of the corals is correlated with microbiome composition and treatments against undesirable pathogens or microbial diversity imbalances have been approached with probiotics following dysbiosis (Peixoto et al., 2021) but engineered bacteria have also been pointed out as a promising alternative (van Oppen et al., 2015; Van Oppen and Blackall, 2019).

Drylands on the other hand, despite the obvious differences in habitat and physical context, share some commonalities with the coral reef scenario. Here too, global warming can trigger community-level shifts as aridity thresholds are reached (Berdugo et al., 2017; 2020). The soil microbiome, along with soil organic carbon is known to be crucial to preserve diversity, which is a major driver of multifuncionality (Delgado-Baquerizo et al., 2016; 2017). Interventions have also been used in the past to restore degraded drylands (Maestre et al., 2006; Bowker 2007). Several studies involving inoculation of cyanobacteria have shown the potential for recovery, at least over some time scale (Abed et al., 2009, Sabarinathan et al., 2021). But a more directed approach using microbial bioengineering could provide a more efficient and durable effect (Maestre et al., 2017; de Lorenzo 2022). How reliable is this alternative and what can we expect from this kind of intervention? Are there unintended effects (such as extinction cascades) associated to such intervention? What kind of guidance can we get from a population-level study? Mathematical and computational models suggest that some well-defined “motifs” can be implemented to modify the location and nature of tipping points (Solé 2015a, 2015b; Solé et al., 2018; Vidiella et al., 2018, 2020). But these are low-dimensional models and cannot be used to test the potential effects of synthetic strains on a multi species system.

In a previous work (Maull and Solé, 2022) we explored the problem of engineering a given species from a resident microbiome in such a way that the modified strain would be reintroduced afterwards. This was a first-order approach to the more general problems outlined above. By design, we considered a modification of the resident community where the synthetic strain is obtained from a small change of one of its members. As a result, the synthetic organism would inherit most (if not all) the links of the wild type used as a target. Such a requirement is relevant because we start from a stable community with an associated ecological network. In this context, the introduced strain would be expected to be almost redundant and remain in place or replace the wild type. Such scenarios were shown to occur under very general conditions within a multispecies competition model (Maull and Solé, 2022). However, the engineered strain was not performing a specific functionality. In this paper we explore this situation using a resource-consumer multi species model, where we target one species that will be modified to reduce the loss of one of the resources. As will be shown below, this approach does not produce diversity losses while an improvement in the abundance of the key resource is obtained.

2. Methods

2.1. Consumer-Resource Model

This study builds upon the established framework of consumer-resource models, specifically drawing from the foundational work of MacArthur in the 1970s (Macarthur, 1970, Chesson, 1990). However, we have made slight modifications to the model to accommodate our synthetic intervention requirements. We propose a set of differential equations that govern the model, consisting of one equation for resource abundance and another for population density over time. The consumer-resource equations can be expressed as follows:

$$\frac{d{R}_k}{dt}={\rho }_k-\left(\sum _i^n{\xi }_{ki}{N}_i\right)R_k-F(R_k,N_{syn})$$

(1)

$$\frac{d{N}_i}{dt}={\eta }_i{N}_i\left(\sum _k^m{\xi }_{ik}{R}_k\right)-{\sigma }_i{N}_i+ϵ$$

(2)

In the proposed model, the resource abundance is denoted by $R_k$, where $k=1,2,...m$, and the population abundance of each species is denoted by $N_i$, where $i=1,2,...n$. The growth rates for each resource and species are represented by $\rho \in (0,1)$ and $\eta \in (0,1)$, respectively. Notably, $\rho$ is not resource-dependent and is treated as a general growth rate for each inorganic resource. In contrast, $\eta$ represents classic replication. The growth of species is also dependent on resource uptake, which is mediated by a matrix ${\xi }_{ki}$ that weighs the strength of the interaction between each resource and its consumer. The terms of the resource-consumer matrix $\xi$ come up from a uniform rectangular distribution U [0, 1], with a connectivity $C=0.3$.

Species decay rate is represented by $\sigma \in (0,1)$, and to prevent unrealistically low biomass, a small immigration factor of $ϵ=0.01$ is introduced. Resource decay is driven by resource uptake by the species, mediated by the $\xi$ matrix, and a decay rate function $F(R_k,N_s)$.

$$F(R_k,N_s)=\frac{{\beta }_k{R}_k}{1+{\delta }_{k\alpha }\frac{N_s}{\lambda }}$$

(3)

The function $F(R_k,N_{syn})$ introduces a non-linear modulation of the resource degradation process. Specifically, the presence of a synthetic species in the community triggers a damping effect on the degradation of a specific resource $\alpha$. The maximum decay rate for a given resource k is represented by $\beta \in (0,1)$, while the constant $\lambda \in (0,1)$ represents the rate of inhibition of the degradation process ($\beta$) due to the presence of the synthetic strain, $N_s$. Lastly, ${\delta }_{k\alpha }$ is the Dirac’s delta function, i. e. ${\delta }_{k\alpha }=1$ when $k=\alpha$ and zero otherwise. Therefore, the damping process will only affect the targeted resource.

2.2. Synthetic Invasion

The core idea is to numerically test the effect of a synthetic invader, under a fundamental control framework, if such organism can have an influence in the extant environment. We take as a fundamental control framework the ideas particularly outlined in (Maull and Solé, 2022). The underlying idea, in contrast with standard microbial invasions (Mallon et al., 2015; Vila et al., 2019) we propose a designed invasion scenario that considers both cell-level features as well as higher-order phenomena connected with community diversity. In order to modify ecosystems by means of introducing new functions using synthetic biology, there are key aspects worth implementing in the design process. The use of an extant species (wild type) as a chasis for the synthetic strain (that will be introduced) can provide of control and effective colonisation outcome. The synthetic, under this specific layout, will have the same or almost the same community interactions than the wild type, therefore the spread potential will be constrained by its niche. At the same time the chances of successful establishment can increase substantially.

Here we use a resource-consumer model to make explicit the relationship between the biotic community and the abiotic resource environment. The species compete for the resource pool through a bipartite network $\xi$, as specified in Equations (1) and (2). The synthetic species is thus deployed within a resident community ${\mathcal{C}}_R$, therefore creating a new, synthetic community ${\mathcal{C}}_S$. Because of our designed intervention, the synthetic species will essentially share the same interactions and similar growth rate than the wild type strain. The simplest formulation is to assume that ${\xi }_{sk}={\xi }_{wk}$, but it seems reasonable to consider that a genetic modification process on a given strain can potentially change the metabolic balances and therefore alter their R-C interactions. To incorporate these deviations, we randomly shift the matrix elements by adding a noise term, i. e. ${\xi }_{sk}={\xi }_{wk}\pm \Delta \xi$, with $\Delta \xi \in (0,0.1)$ a random number with a uniform distribution. The functional impact of the synthetic species will take place on a randomly chosen resource $R_{\alpha }$. This effect, as discussed above, involves a dampening in the resource degradation rate, as specified in Equation (3). Additionally, the engineered strain is likely to have a fitness that might be different from the wild type. It can have a direct positive feedback effect on the metabolic performance or conversely, it can imply a metabolic burden. To incorporate these two options we modify the efficiencies using ${\eta }_s={\eta }_w\pm 0.1{\eta }_w$. The neutral scenario (${\eta }_s={\eta }_w$) will also be considered.

The intervention approach thus starts with an m-resources, n-species ecosystem. Each population will start the simulation with a small initial value of $N_i\left(0\right)=0.01$ and an initial resource amount of $R_k\left(0\right)=0.01$. The community state ${\mathcal{C}}_R$ is followed over a transient of ${\tau }_1=800$ time steps. At that moment the community is considered stable and is therefore invaded with an extra species, $N_s$, which will be based on a randomly chosen wild type $N_w$. The chosen mold species must have a population $N_w\left({\tau }_1\right)>0.2$. We invade with a small inoculation of synthetic strain as $N_s\left({\tau }_1\right)=0.01$. It is important to notice that the set of species that compose the ${\mathcal{C}}_R$ does not change in terms of diversity because there is an immigration factor of $ϵ$ that prevents for total extinctions. Once the ${\mathcal{C}}_R$ is stable, the invasion takes place, leading the community to the synthetic state ${\mathcal{C}}_S$, after time ${\tau }_2=2000$ time steps, until the system reaches stability again. The R-C interaction matrix goes from:

$${\mathcal{X}}_R=\left(\begin{array}{cccc}{\xi }_{1,1}& {\xi }_{1,2}& \cdots & {\xi }_{1,n}\\ {\xi }_{2,1}& {\xi }_{2,2}& \cdots & {\xi }_{2,n}\\ ⋮& ⋮& \ddots & ⋮\\ {\xi }_{m,1}& {\xi }_{m,2}& \cdots & {\xi }_{m,n}\end{array}\right)$$

That is replaced by the “synthetic” R-C interaction matrix, namely

$${\mathcal{X}}_R=\left(\begin{array}{cccc}{\xi }_{1,1}& {\xi }_{1,2}& \cdots & {\xi }_{1,n+1}\\ {\xi }_{2,1}& {\xi }_{2,2}& \cdots & {\xi }_{2,n+1}\\ ⋮& ⋮& \ddots & ⋮\\ {\xi }_{m,1}& {\xi }_{m,2}& \cdots & {\xi }_{m,n+1}\end{array}\right)$$

Alternatively, it has been also considered the control scenario where the invader is a random species. Meaning by random that it does not come from a wild type species present in the community. Therefore, such introduced strain will have unknown random connections in the $\xi$ matrix and a random $\eta$ growth rate. We considered two interesting scenarios under this approach. First, the random invader do positively affect a randomly chosen resource, just at it is depicted in Equation (3). This scenario mimics the introduction of a enginered microorganism (not based on an extant species) in order to perform the restoration task. Second, the random invader do not have a specific positive interaction with any resource.

A statistical analysis of the synthetic invasion impact is performed by determining the population change of the community from ${\tau }_1$ to ${\tau }_2$. Population change in terms statistical biotic increase of the present species. Both the wild type and the synthetic species are excluded from the change counting. We generate 500 random stable resident communities ${\mathcal{C}}_R$ of $m=10$ resources and $n=20$ species. Each ${\mathcal{C}}_R$ community is randomly generated. Once each community is stable, they are subsequently invaded. To simulate the dynamics of the resource-consumer model, we used a fourth-order Runge-Kutta algorithm.

3. Results

In this section we will study the general properties of the model described by Equations (1) and (2). The first step is to consider the simplest model that captures the main assumptions, but with some extra assumptions that allow to obtain a one-dimensional version. The general model is then considered using numerical experiments where the synthetic strain is introduced within a stable RC community and the statistical patterns of population change are analysed and compared with those associated to a random invader performing the same functionality.

3.1. Two Resources and Two Consumers

Before we proceed to the general analysis of the multidimensional model, let us consider a simpler scenario that allows to perform some analytics and define the bounds required for a synthetic strain to persist. As defined above, this engineered microorganism is designed to enhance the stability of one of the resources by reducing its loss. The basic diagram of this model is summarised in Figure 2. In this section we exploit the simplicity of this low-dimensional scenario along with an extra constraint that allows to define basic inequalities required for the synthetic strain to be successful.

For simplicity, we fix ${\xi }_k=1$ and assume a constant population constraint for the consumer populations, i.e., $N_s+N_w$ is a fixed quantity. This choice allows for reducing the initial system to a single-equation model. The equations now read:

$$\begin{array}{c}\hfill \frac{d{R}_1}{dt}={\rho }_1-(N_w+N_s)R_1-\beta R_1\end{array}$$

(4)

$$\begin{array}{c}\hfill \frac{d{R}_2}{dt}={\rho }_2-(N_w+N_s)R_2-\frac{\beta R_2}{1+\frac{N_s}{\lambda }}\end{array}$$

(5)

for the resource dynamics and we also consider a modified set for the microbiome dynamics, namely:

$$\begin{array}{c}\hfill \frac{d{N}_w}{dt}={\eta }_w{N}_w(R_1+R_2)-\sigma N_w-N_w\Phi \end{array}$$

(6)

$$\begin{array}{c}\hfill \frac{d{N}_s}{dt}={\eta }_s{N}_s(R_1+R_2)-\sigma N_s-N_s\Phi \end{array}$$

(7)

where we indicate $\Phi =\Phi (N_w,N_s)$.

A general analysis is given in the SM without the previous assumptions, but leading to the same basic result. Using this simplified model, we can further reduce its dimensionality by assuming that the dynamics of the resources is fast and thus we can consider that

$$\left(\frac{d{R}_1}{dt}\right)=\left(\frac{d{R}_2}{dt}\right)\approx 0$$

(8)

and thus write the resource terms as a function of microbiome populations, i.e., after some algebra we obtain

$$R_1=\frac{{\rho }_1}{1+\beta }$$

(9)

for the first resource and

$$R_2\left(N_s\right)=\frac{{\rho }_2(\lambda +N_s)}{(1+\beta )\lambda +N_s}$$

(10)

Now, because of the assumption that the total population N remains the same, we consider for convenience that $N=1$ and, from the condition $dN/dt=d(N_w+N_s)/dt=0$, it is possible to show that:

$$\Phi (N_s,N_w)=({\eta }_s{N}_s+{\eta }_w{N}_w)(R_1+R_2)-\sigma$$

(11)

Our main goal in this context is to guarantee that the synthetic strain can persist and define the conditions for such persistence. Using $N_w=1-N_s$ and the previous result, we have now, after rearranging some terms we obtain:

$$\frac{d{N}_s}{dt}=(R_1+R_2)({\eta }_s-{\eta }_w)(1-N_s)N_s$$

(12)

Since $R_1+R_2>0$, this single equation provides a threshold condition for growth of the synthetic strain ($d{N}_s/dt>0$), which in this case is simply namely

$${\eta }_s>{\eta }_w.$$

(13)

This result, under all the simplifications involved, suggests that a synthetic microorganism able to reduce the loss of a common good will likely establish itself provided that its efficiency of resource use is large enough. Our derivation seems far from the multispecies problem, but -as shown below- this result turns out to be robust.

A more general case can be considered if deviations in the interaction matrix between the wild and synthetic components occur (as expected) when a genetic modification is introduced. Therefore upgrading the coarse-grained model with ${\xi }_{wk}\ne {\xi }_{sk}$ states an scenario worth exploring. Therefore, the quartet of equations is the following:

$$\begin{array}{c}\hfill \frac{d{R}_1}{dt}={\rho }_1-({\xi }_{wk}{N}_w+{\xi }_{sk}{N}_s)R_1-\beta R_1\end{array}$$

(14)

$$\begin{array}{c}\hfill \frac{d{R}_2}{dt}={\rho }_2-({\xi }_{wk}{N}_w+{\xi }_{sk}{N}_s)R_2-\frac{\beta R_2}{1+\frac{N_s}{\lambda }}\end{array}$$

(15)

Rearranging $\psi =\eta \xi$:

$$\begin{array}{c}\hfill \frac{d{N}_w}{dt}={\psi }_{wk}{N}_w(R_1+R_2)-\sigma N_w-N_w\Phi (N_w,N_s)\end{array}$$

(16)

$$\begin{array}{c}\hfill \frac{d{N}_s}{dt}={\psi }_{sk}{N}_s(R_1+R_2)-\sigma N_s-N_s\Phi (N_w,N_s)\end{array}$$

(17)

We can reduce again the dimensionality by using the fast dynamics condition for resources(8) along with the constant population constraint. It can be shown that:

$$\Phi (N_s,N_w)=({\psi }_{sk}{N}_s+{\psi }_{wk}{N}_w)(R_1+R_2)-\sigma$$

(18)

And finally rearranging the synthetic strain differential equation to the following expression. For the sake of clarity we recover $\eta \xi =\psi$:

$$\frac{d{N}_s}{dt}=(R_1+R_2)({\eta }_s{\xi }_{sk}-{\eta }_w{\xi }_{wk})(1-N_s)N_s$$

(19)

Notice that now $(R_1+R_2)$, assuming again fast dynamics, reads:

$$\frac{{\rho }_1}{\beta -N_s({\xi }_{wk}-{\xi }_{sk})+{\xi }_{wk}}+\frac{{\rho }_2}{\frac{\beta \lambda }{N_s+\lambda }-N_s({\xi }_{wk}-{\xi }_{sk})+{\xi }_{wk}}$$

(20)

Even though it changes in relation to (9,10), since populations are normalized to one, the sum still positive for all scenarios. Which, overall, provides a new threshold condition, namely:

$${\eta }_s>\left(\frac{{\xi }_{wk}}{{\xi }_{sk}}\right){\eta }_w$$

(21)

Here, the $\xi$ term plays a multiplicative role with the intrinsic growth rate. In the multispecies extrapolation, in general terms, it can be stated that the pairwise dynamics between the wild type and the synthetic will fall on the same space of the possible: competitive exclusion, failure and coexistence (when they are dynamically the same species), depending on the $\eta \xi$ relation.

3.2. Multiple Species and Resources

Consider now the full model, as described by the set of Equations (1) and (2), along with the functional modulation of one resource given by (3). This section aims to study, under general conditions, the effect of different strains (engineered or not) that invade an extant community. Two types of interventions are considered, namely synthetic and random. We define as synthetic a genetically modified strain obtained from a wild-type species present in the resident community. The synthetic invaders can display different growth rate performances in relation to the wild-type organism. They can present a metabolic boost, a metabolic burden or remain unchanged. On the contrary, we define those strains as random, genetically modified or not, and thus not based on any previous existing wild-type species of the extant community. Therefore all the attributes are entirely random.

Two different families of results are presented. We numerically generate the time series of the deployment scenarios. Those time series highlight the pairwise dynamics between the wild type and the synthetic, and the possible impacts of a random invasion on the community. By analysing multiple runs of these time series, a statistical analysis is performed by measuring the frequency distributions of population shifts. Specifically, the population change of the j-th species after inoculation is computed as:

$$\Delta \left(N_j\right)=N_j\left({\tau }_1\right)-N_j\left({\tau }_2\right)$$

(22)

associated to the reconfiguration event post invasion, accumulated through all the statistical trials. Here ${\tau }_1$ is the last time step before the introduction of the synthetic, once the resident population has reached its stable state, and the final population is computed after ${\tau }_2$ steps (see methods).

In Figure 3a–c we display the three outcomes emerging from our intervention scenario as described above. As predicted by the two-species, two-resource scenario, the efficiency of the introduced strain largely determines the result of the bioengineering. However, now we can also evaluate its impact on community diversity: is the synthetic strain triggering community shifts? The answer seems to be negative, as illustrated by the time series. Here we use thick orange and blue to represent wild type and synthetic populations, respectively, while all other species abundances are shown in green. The resources are draw in black and the enhanced resource ($R_{\alpha }$) in light blue. In (a) we show a typical time series for the symmetric scenario (${\eta }_{wt}={\eta }_{syn}$). As we appreciate, the two populations converge to the same population level. In contrast, if ${\eta }_{syn}={\eta }_{wt}+0.1{\eta }_{wt}$ (b), a switch occurs. Finally, in (c) the growth rate of the synthetic is a $10\%$ smaller than the wild type due the extra metabolic load (I.e. ${\eta }_{syn}={\eta }_{wt}-0.1{\eta }_{wt}$). In this case, the synthetic fails to establish itself. In each run, a random wild type species has been chosen from the pool with a population greater than $0.2$ and the initial $N_s\left({\tau }_1\right)=0.1$. Here the connectivity is $C=0.3$. Notice that in both scenarios b and c, the species that end up failing, are not complete washed out due to the (small) immigration term $ϵ$ in Equation (2).

What is the statistical nature of the previous patterns? Specifically, we want to analyse the expected impact of the synthetic “invader” on community structure. Exploring this question is strongly tied to the problem of whether such intervention can trigger “unintended consequences”. In Figure 4, we display the distribution of changes in population density under the three different growth performances under the two previous matrix choices. The plots are represented in linear-log scale, thus helping to detect asymmetries. But they also point into a first important result: the vast majority of population shifts are very small, thus indicating that the introduced strains have little impact at this level. In all plots we have indicated the domain of $\Delta \left(N_j\right)\le -1$. These are associated to those (rare) cases where population shifts push the given species to extinction. However, since we use a constant immigration rate $ϵ$, this corresponds to a turnover event. In any case, the fraction of events under this threshold gives us a measure of the frequency of community changes that involve a marked reduction of a component at time ${\tau }_1$.

In the first row, a, b and c correspond to the (neutral) unchanged growth rate, metabolic boost and metabolic burden respectively. All of them entail symmetric interactions between the wild type and the synthetic (${\xi }_{sk}={\xi }_{wk}$). In the same figure, d, e and f correspond to the same scenarios regarding the intrinsic growth rate, but the R-C interactions between the species is slightly randomized. Interestingly, there is a robust statistical skewness towards the positives in all scenarios, which indicates a statistical biotic increase as a result of the synthetic intervention. In other words, the introduced strain triggers population improvements over population decreases. This skewness is more dramatic in the symmetric case (a, b, c).

The positive statistical biotic increase gets highlighted when compared against a random invader, see Figure 5. Here we find two scenarios which correspond to two interesting cases. The first one (a) corresponds to a random invader in all regards. However, the newcomer still possesses the ability to reduce a randomly chosen resource decay rate (here indicated by an arrow with a positive impact). Therefore representing a potentially engineered species, i.e., synthetic, yet not a copy of someone in the community. Secondly (b), we display the results for a random invader with no effect on resources. In both cases, the results show a dramatic difference. The skewness in the distribution of changes drives to the negatives, depicting a statistically biotic decrease resulting from the invasion.

The comparison between the synthetic and random scenarios reveals important insights. In the first set of cases, it is observed that the synthetic invader consistently thrives. Whether the growth rate remains unchanged or experiences a boost, the synthetic invader successfully establishes itself by either coexisting with the wild type or by outcompeting it. On the contrary, the random invader fails more frequently. The frequency distributions reveal a high kurtosis value can be expected, since the distributions peak around zero (notice that the frequency is in logarithmic scale). The synthetic intervention scenarios exhibit right-skewed distributions, indicated by positive skewness values, while the random scenarios are left-skewed.

Notably, the symmetric scenarios (4a-b) exhibit higher skewness values than the asymmetric scenarios (4d-e). However, among the asymmetric scenarios, when the synthetic experiences a metabolic burden (Figure 4f), it demonstrates the highest positive skewness compared to all cases. This finding raises intriguing questions regarding the potential impact of positive and transient interventions. Can temporary synthetic invasions enhance the ecosystem at a broader scale and subsequently disappear, leaving behind a reorganised and healthier ecosystem? Furthermore, does this approach align with concerns regarding the control of organism release in natural environments? These problems are addressed again in the discussion below.

4. Discussion

Microbial bioengineering might be a necessary step to improve the resilience of some ecosystems that face major challenges associated to climate change. Can this strategy improve ecosystem reliability? Will instead engineered strains trigger undesirable, community-level effects that prevent their use within natural communities? The truth is that, despite the sometimes heated debate concerning the use of engineered strains in conservation biology, no real test nor serious theoretical attempts have been promoted to truly address the previous questions. As a consequence, the potential outcomes are poorly understood. Moreover, there is a chasm between most of the research approaches used within synthetic biology (dominated by the molecular and cellular scales) and those ecosystem-level studies that consider the community dynamics on meso- or global scales.

In this paper we have introduced a model of resource-consumer interactions aimed at testing the potential population-level changes when an engineered microbial strain has been introduced. In contrast with a previous paper (Maull and Solé, 2022), where the new strain simply competed with the resident community, here a specific functional improvement of the resource pool has been included. Specifically, for a given resource, a synthetic strain designed from a member of the resident community acts by reducing its decay or spillover. Our results indicate that the impact is typically positive, dominated by small population changes, consistently with our previous findings. This suggests that, against the idea of unintended consequences, a properly designed bioengineering can fulfil to important goals: modify what is needed while preserving diversity. In this context, this is an engineering strategy that exploits emergent properties: by maintaining (and even increasing) biodiversity, multifunctionality will also be protected. By contrast with standard ecological engineering (Odum and Odum, 2003), here a global property (diversity) is also the target. In other words, our strategy includes a specific top-down controlled feature along with a bottom-up, emergent property. This is, we believe, an important conceptual novelty of the proposed intervention that is grounded in a systems-level view of engineering where emergent dynamics is at the center (Holling 1996; Solé 2015; Krakauer, 2019; Gorochowski et al., 2020).

We consider these results a starting point to a more ambitious roadmap to further ways to terraform ecosystems. Many open problems remain. On the theoretical side, evolutionary, spatial and network dynamics and the architecture of the environment have been neglected. Space, for example, can play a key role as a biodiversity enhancer. One key result from the early days of spatial ecology is that complex spatial dynamics (not considered here) can be crucial to maintain ecological diversity (Solé et al., 1992; Bascompte and Solé 1995; see also Pierce et al., 2020). In this context, the trade offs between resource and species dynamics in both space and time will need to be integrated in a more rich context where stochastic dynamics might also play a relevant role. Network architecture, on the other hand, is known to influence the propagation of changes across trophic or mutualistic webs (Montoya et al., 2011). Further exploration needs to be done in relation with realistic web architectures. Similarly, we have reduced our bioengineering design to a single species change, but potential multi species consortia might be of relevance to achieve the desired effects, perhaps targeting different resources. Above all, experimental tests using micro- and mesocosm frameworks are much needed to put to test the theoretical ideas outlined here.