^{1}

^{*}

^{2}

^{3}

^{1}

^{4}

^{3}

^{1}

^{2}

^{3}

^{4}

Edited by: Thierry Tonon, Station Biologique de Roscoff, CNRS-UPMC, France

Reviewed by: Vangelis Simeonidis, University of Luxembourg, Luxembourg; Georg Basler, Consejo Superior de Investigaciones Cientificas, Spain

*Correspondence: Marianna Taffi, School of Biosciences and Veterinary Medicine, University of Camerino, Via Gentile III Da Varano, 62032 Camerino, Italy e-mail:

This article was submitted to Systems Biology, a section of the journal Frontiers in Genetics.

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

The pressure to search effective bioremediation methodologies for contaminated ecosystems has led to the large-scale identification of microbial species and metabolic degradation pathways. However, minor attention has been paid to the study of bioremediation in marine food webs and to the definition of integrated strategies for reducing bioaccumulation in species. We propose a novel computational framework for analysing the multiscale effects of bioremediation at the ecosystem level, based on coupling food web bioaccumulation models and metabolic models of degrading bacteria. The combination of techniques from synthetic biology and ecological network analysis allows the specification of arbitrary scenarios of contaminant removal and the evaluation of strategies based on natural or synthetic microbial strains. In this study, we derive a bioaccumulation model of polychlorinated biphenyls (PCBs) in the Adriatic food web, and we extend a metabolic reconstruction of

Aquatic ecosystems are subject to a mixture of synthetic organic chemicals, leading to adverse effects on organisms at different levels of biological organization and at all trophic levels of the food web. Over the last decades, many removal strategies have been proposed in order to reduce the bioavailability of persistent organic pollutants (POPs) and to limit the consequent bioaccumulation phenomena on species.

However, not all the living organisms in a polluted environment are prone to bioaccumulation. The sizeable variety of marine microbial life is metabolically involved in many transformation processes like biogeochemical cycles of elements, water quality conservation and biodegradation of many organic pollutants. Microbial communities are also an active compartment at the lower trophic levels of marine food webs. They interact with the grazing activities of planktonic groups and play a crucial role in the mineralization of organic matter through the complex trophic pathway known as the

Computational models and predictive tools have found wide applicability and usefulness both in ecotoxicological studies and in the reconstruction of genome-scale metabolic network of pollutant degrading bacteria. However, to the best of our knowledge, these techniques have never been considered for investigating, in a combined way, the multiscale effects of microbial bioremediation at the ecological level. In this work, we develop a computational framework that integrates bioaccumulation information at ecosystem level with genome-scale metabolic models of degrading bacteria. We apply it to the case study of the PCBs bioremediation in the Adriatic food web.

Specifically, we estimate the PCBs bioaccumulation model by using Linear Inverse Modeling, and we employ Flux Balance Analysis to extend the metabolic reconstruction of the toluene degrading bacteria

In our framework we focus on the case of PCBs bioaccumulation in the Adriatic sea, a semi-enclosed basin characterized by high biodiversity (Coll et al.,

We assume that organic chemicals follow the same paths as biomasses, moving via feeding link through the trophic structure of the food web, which is a common approach in the field of ecotoxicological modeling (Hendriks et al., ^{1}

Figure ^{−1} wet weight-based. Biomasses are measured in t km^{−2} wet weight organic matter, and biomass flows in t km^{−2} year^{−1}. PCBs flow rates are thus expressed in mg km^{−2} year^{−1}. We denote the contaminant flow from prey _{i→j}, and the PCBs concentration in _{i}_{i→j}) are known quantities and are estimated as reported in (Taffi et al.,

_{j} c_{j→i} − ∑_{j} c_{i→j} = 0 |

The bioconcentration of a generic group ^{a} |

_{i} ⋈ |

where _{i} |

_{j→i} = _{j→i} · _{j} |

The contaminant flow from group _{j→i} and the PCBs concentration in the source _{j→Export}); respiration flows (_{j→Respiration}), which account for part of the unassimilated fraction of ingested biomass; or removal due to fishing activity, which can be directed to the landings (_{j→Landing}) or to the discards (_{j→Discard}). The latter enters back the biomass cycle and is modeled as a mass-balanced group, with its own bioconcentration value. |

_{Import→i} = _{Import→i} · _{i} |

This class of constraints describes generic imports of PCBs coming from external contaminant inflows (e.g., immigration), which we group in the Import compartment. In this case, the PCBs concentration in the biomass imported into group |

_{Water→i} = _{i}_{Water} |

where _{i}_{Water} is the concentration in water^{b} |

_{i} ≥ 0 |

_{Water} is unknown, the constraint becomes non-linear and _{i}_{Water→i} is treated as a single unknown. We assume null _{i}

Various environmental and biological factors limit the natural PCBs degradation process, among which the high selectivity of bacteria for specific PCBs congeners. Higher chlorinated congeners typically tend to accumulate in marine sediments, where anaerobic bacteria use these compounds by reductive dechlorination as alternative electron acceptors in their respiration processes, thus making PCBs less chlorinated and more aerobically degradable. This step is generally slow but crucial in the whole detoxification process, and various PCBs-dechlorinating bacteria, mainly belonging to the phylum

In this work, we construct a synthetic model of PCBs degrading bacteria using the FBA approach (see Section 2.3), by extending the metabolic reconstruction of

Starting from biochemical reactions and stoichiometric coefficients, the Flux Balance Analysis (FBA) framework is based on the assumption of a metabolic steady state (Orth et al.,

Formally, let _{h}_{k}_{h}_{hk}_{hk}_{hk}

Typically, there are more reactions than metabolites, thus equation ^{⊥}_{k} ≤ _{k}^{⊤}_{k}, where ^{⊥}_{k} and ^{⊤}_{k} are the minimum and maximum flux rates for the ^{n}_{k = 1} _{k}v_{k}_{k}

When two objective functions are taken into account, an FBA problem can be formulated as a bilevel linear programming problem (e.g., for optimizing growth and product yield Burgard et al.,

where _{k1} (e.g., biomass production) and the synthetic objective _{k2} (e.g., contaminant uptake), we set _{k1} = 1 and _{k2} = 1. The solution of the bilevel problem (1) is a pair indicating the maximum natural objective (inner problem) allowed by the constraints ^{⊥}_{k} ≤ _{k}^{⊤}_{k}, and the maximum synthetic objective allowed in the flux distribution that maximizes the natural objective. The bilevel problem can be converted to a single-level problem using the duality theory applied to the inner problem, which is replaced by additional constraints for the outer problem.

We introduce a method for integrating the PCBs bioaccumulation network with the FBA-based metabolic reconstruction of

The basic idea is encoding each link

where the rate of a reaction (_{i,j}) is upper bounded by the original corresponding flow rate _{i→j}. This formulation admits a space of solutions with potentially reduced (even zeroed) contaminant flows, which is required in order to reproduce the contaminant removal by the bacterial metabolism.

Any admissible vector of fluxes for the reactions in _{FW}_{i,j} for any group _{i→j} > 0) but _{i}

Additionally, we consider the following set of exchange reactions for expressing the external inputs and outputs of the food web:

The set _{FW}_{FW}_{j} c_{i→j}). Note that it is sufficient to constrain the import reactions in order to obtain a consistent FBA encoding of the bioaccumulation model. Indeed, by mass-balance, throughflow values are conserved by the encoding and it can be shown that the food webs entailed by the reactions in _{FW}_{FW}_{FW}

Finally, for a generic group _{i}

An effective way to accomplish this task is adding a dummy metabolite

We define the ^{2}

where _{FW}_{FW}_{FW}

Evidently, not all integrations are ecologically and biologically plausible. In our model, we consider two bioremediation scenarios, as also shown in Figure

The integrated models have been obtained after converting PCBs flows (mg km^{−2} year^{−1}) to the flux units used in FBA (mmol h^{−1} gDW^{−1}). The conversion factor is ^{−1}; 4-Chlorobiphenyl: 188.6529 mol g^{−1}); ^{2}). ^{−3} gDW, enough to import the totality of the connected PCBs flows into the

In order to assess the effects of bioremediation on our contaminated food web, we combine the evaluation of bioconcentrations with the study of ecological network indices. Typically, global indices (Kones et al.,

FBC gives the topological importance of a species in maintaining the flows among other groups. The FBC of a group _{i}

where _{G} c_{j→k} is the maximum flow between _{G\i} _{j→k} is the maximum flow between

We employ

where

The approach for the estimation of the PCBs bioaccumulation model and for the analysis of network indices was implemented in R (using packages

By applying bilevel FBA, the growth rate of ^{−1}) for PCBs uptake rate up to 9.8 mmol h^{−1} gDW^{−1} (Figure ^{−1} gDW^{−1}, since the rate of PCBs uptake stays constant for upper bounds greater than this value. Therefore, optimal growth is maintained until the rate of PCBs uptake is almost at its maximum, while for uptake rates greater than 9.8 mmol h^{−1} gDW^{−1}, biomass production drops to 71% of its optimal value (1 h^{−1}). Further, the addition of the PCBs bioremediation pathways to the

^{−1} gDW^{−1}. The maximum PCBs uptake rate is 10 mmol h^{−1} gDW^{−1}, and the optimal growth rate is thus achieved for almost the whole range of PCBs uptake. ^{−1} gDW^{−1}. The symmetric bilevel problem (with toluene limited from 0 to 20 mmol h^{−1} gDW^{−1}) gives the same linear front. This tradeoff delineates two phenotypes in the PhPP analysis (L2: biomass, L1: toluene+PCBs uptakes): in the lower half (green region), we have optimal growth; in the upper half (blue region), growth is limited to 71% of the optimal growth.

In order to investigate the relationship between growth rate and oxygen uptake, and between PCBs uptake and oxygen uptake, we apply a single-level FBA analysis. In Figure ^{−1} gDW^{−1}, and then remains stable even at high oxygen uptake. The

We also analyse the interdependence between the PCBs degradation pathways (introduced in this work) and that of toluene (in the original reconstruction). We derive an optimality front between them by solving two bilevel problems. In the first, we evaluate the maximum toluene uptake when PCBs imports are favored, while in the second problem, we consider the symmetric objectives. Both problems identify the identical linear trade-off (red dashed line in Figure ^{−1} gDW^{−1}, toluene flux > 19.7 mmol h^{−1} gDW^{−1}), as also seen for the PCBs case in Figure

We analyse the integrated models obtained by applying the two scenarios introduced in Section 2.4. In Scenario 1, microbial degradation pathways reduce contaminant concentrations through outflows from natural detritus and fishing discards (functional groups 38 and 39, trophic level=1), simulating a bioremediation at the level of the microbial loop. In Scenario 2, PCBs bioremediation is assumed to act in the water compartment by reducing simultaneously all the PCBs uptakes in each functional group. The following results are obtained by solving the bioremediation problem (Section 2.4) and computing bioconcentrations and network indices on the resulting (the entailed) bioaccumulation networks.

In Figure

In Scenario 1 [plot (b)], we clearly notice a simpler pattern of PCBs contamination among functional groups, due to a considerable reduction of feeding links active in the transport of PCBs. Specifically, the redirection of outflows from detritus and discard out of the food web causes the inactivation of several PCBs flows, and subsequent presence of groups with null PCBs concentration. Therefore, these groups (not plotted) are disconnected from the bioaccumulation network but still active in the biomass network. They include detritus (38) and discard (39); detritivores (6,7,8); group 30, which only feed on planktonic groups; and groups feeding on those mentioned so far (33, 34, 35). Other variations are detectable in groups 5, 27, and 28 (detritivores and planktivores) that no longer acquire PCBs from food. Finally, we can observe that groups 12 and 27 gain in this scenario a central role in the contaminant diffusion, becoming the preferential source of most of their predators, while in Scenario 0 their contribution appears less relevant.

The bioaccumulation network under Scenario 2 [plot (c)] exhibits a similar structure to Scenario 1. Detritus groups (38,39) and detritivores (6,7,8) are no longer connected to the rest of the food web, showing that bioremediation of the water compartment tends to disrupt the pathways of contaminant uptake at the lowest trophic levels. Another similarity is the promotion of group 12 as a central node in the acquisition of PCBs by its predators. On the other hand, group 27 has no outgoing PCBs flows, while in plot (b) the opposite situation (no inflows) is observed for the same group. In general, we notice a lower number of active links with respect to Scenario 1, especially in species at higher trophic levels.

Another kind of analysis enabled by our framework is the study of the networks obtained by solving the bioremediation problem at increasing efficiencies, limiting the amount of PCBs flow allowed into the bacterial metabolism. In both scenarios, we analyse the variations in PCBs bioconcentrations (Figures ^{−2} year^{−1} and 3312 mg km^{−2} year^{−1}, respectively), which mainly depends on the structure of the network. Applying the conversion factor in Section 2.4, the maximum remediated flow in Scenario 1 corresponds to a PCBs uptake rate of 3.15 mmol h^{−1} gDW^{−1} (31.5% of the maximum uptake), while in Scenario 2 to 2.45 mmol h^{−1} gDW^{−1} (24.5% of the maximum uptake).

As regards Scenario 1, no remarkable reductions in bioconcentrations are observed in the entire food web [see plot (a)], apart from the natural detritus and the discard groups, whose PCBs values are zeroed at 89 and 16%, respectively, of the maximum bioremediation efficiency. We register only minor drops in a number of groups at TL 4 (14, 16, 23, 24) and in group 37 (feeding on discard). This tendency is also visible in the sum of PCBs, which is practically constant.

On the other hand, Scenario 2 gives a considerable decrease in the bioconcentrations of all groups [plot (b)]. This is explained by the fact that the estimated uptakes from water constitute an important fraction of imported contaminant, whose degradation also mitigates, indirectly, uptakes from food. The only exceptions are natural detritus and discard (38, 39), which have null PCBs flow from water (see Table

In Scenario 1, the analysis of the FBC index [plot (c)] highlights the topological importance of natural detritus in the bioaccumulation network, which derives from the fact that every group in the food web contributes (via natural death and unassimilated food) to its contaminant uptake. Indeed, this group maintains its central role up to 89% of bioremediation efficiency. After this point, a structural disruption occurs, related to the detritus becoming disconnected from the network (i.e., no incident flows). This leads to cascade effects also in the centrality of groups 13, 15, 16, 19, and 32. Apart from this case, FBC exhibits quite a robust pattern, showing a number of groups (2, 21, 22, 25) with unchanged centralities regardless of the amount of bioremediated flux. This structural robustness is evidenced also by the link density values, indicating that, globally, the number of links active in the PCBs diffusion is relatively constant.

On the contrary, Scenario 2 [plot (d)] produces prominent changes in the centrality of most species. Here, natural detritus loses its dominant role in the network at 10% of maximum bioremediation. Moreover, at 34% of efficiency, we observe a sudden fall in the FBC of group 24, as also registered on its bioconcentration values [see plot (b)]. Only functional groups 2 and 25 show robust topological importance, in agreement with Scenario 1. The evolution of the link density index also evidences the high sensitivity of the network structure. Indeed, the index reaches an average of 3.1429 active links per group, 36% lower than the initial value.

Recent biotechnological advances and novel discovery tools in marine metagenomics are paving the way for new integrated solutions in environmental bioengineering, turning empirical hypotheses into practical methods. In this context, we presented a computational framework for the analysis of contaminated ecosystems and for the evaluation of different hypothetical bioremediation scenarios. We considered the case of PCBs bioaccumulation in the Adriatic food web and PCBs degradation mediated by

Our computational experiments indicated that the extended

To the best of our knowledge, this is the first computational method linking genome-scale reconstructions of bacterial metabolism with food web bioaccumulation models for designing and analysing bioremediation strategies. Approaches based on high dimensional omics data and network inference methods (Perkins et al.,

In this work we wanted to stress a different view on marine ecosystems, regarding them not just as ensembles of macro-species, but as complex multiscale networks linking classical food webs and microbial groups, toward a new perspective of “eco-metabolic” networks. We believe that bringing the study of microbial metabolic activity into the field of ecotoxicological modeling can highlight bottlenecks and advantages of different bioremediation approaches and shed light on the ecological role of marine microbial life. Furthermore, PCBs degrading bacteria live in communities structurally organized in biofilms (Abraham et al.,

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

^{1}Being a linear operation, the mean of valid solutions to a system of linear constraints is in turn a valid solution to the system (see also van Oevelen et al.,

^{2}Due to the mass balance assumption of FBA, influxes could have been equivalently considered.