Design and optimization of microbial communities

Universidade do Minho Escola de Engenharia Sophia Torres Santos Design and Optimization of Microbial Communities abril de 2024 UMinho | 2024Sophia Torres Santos Design and Optimization of Microbial Communities Universidade do Minho Escola de Engenharia SOPHIA TORRES SANTOS DESIGN AND OPTIMIZATION OF MICROBIAL COMMUNITIES TESE DE DOUTORAMENTO DOUTORAMENTO EM BIOENGENHARIA TRABALHO REALIZADO SOB ORIENTAÇÃO DE: PROFESSORA DOUTORA ISABEL CRISTINA DE ALMEIDA PEREIRA DA ROCHA E DE DOUTOR OSCAR MANUEL LIMA DIAS abril 2024 ii ESTE É UM TRABALHO ACADÉMICO QUE PODE SER UTILIZADO POR TERCEIROS DESDE QUE RESPEITADAS AS REGRAS E BOAS PRÁTICAS INTERNACIONALMENTE ACEITES, NO QUE CONCERNE AOS DIREITOS DE AUTOR E DIREITOS CONEXOS. ASSIM, O PRESENTE TRABALHO PODE SER UTILIZADO NOS TERMOS PREVISTOS NA LICENÇA ABAIXO INDICADA. CASO O UTILIZADOR NECESSITE DE PERMISSÃO PARA PODER FAZER UM USO DO TRABALHO EM CONDIÇÕES NÃO PREVISTAS NO LICENCIAMENTO INDICADO, DEVERÁ CONTACTAR O AUTOR, ATRAVÉS DO REPOSITÓRIUM DA UNIVERSIDADE DO MINHO. Licença concedida aos utilizadores deste trabalho Atribuição CC BY https://creativecommons.org/licenses/by/4.0/ iii 1. AGRADECIMENTOS E que aventura fantástica está prestes a terminar! Mas antes de terminar quero agradecer a todos aqueles que tornaram este trabalho possível. Agradeço à Professora Doutora Isabel Rocha. Obrigada por todos os conhecimentos científicos e pessoais transmitidos, pela constante disponibilidade e acompanhamento. Pelo incentivo, motivação e confiança depositada. Obrigada por me ajudar a crescer e a abrir horizontes. Será sempre um exemplo para mim. Agradeço ao Doutor Oscar. Por ter aceitado o desafio de me orientar, de me aconselhar, de me motivar. Por me estimular a perseguir outros desafios e a dar sempre mais de mim em cada tarefa proposta. Agradeço à instituição de acolhimento, Centro de Engenharia Biológica da Universidade do Minho, por me ter proporcionado as condições necessárias à realização deste trabalho. Agradeço também à Fundação para a Ciência e Tecnologia pelo financiamento atribuído através da bolsa de doutoramento SFRH/BD/121695/2016. Agradeço em especial a todos os parceiros do Instituto de Tecnologia Química e Biológica António Xavier, Instituto Gulbenkian e Ciência e Universidade dos Açores, em especial à Ana e ao Ricardo, que sem o seu trabalho e ajuda este trabalho não seria possível. Agradeço a todos os elementos do grupo BioSystems, presentes e passados. Obrigada pelo bom ambiente de trabalho que sempre proporcionaram, pelas discussões científicas e partilha de ideias. Um agradecimento ao Vítor Pereira pela ajuda na parte da implementação dos métodos de otimização no contexto de comunidade microbianas. À Joana que me acompanhou desde o início nesta caminhada doutoral. Obrigada pela partilha de conhecimentos, pelas longas horas de conversa, pelas gargalhadas, pelas conversas sérias e pela amizade conquistada. À Ana e à Sara e “às frizes” que bebemos e que espero que continuemos a beber. À minha família e amigos. Por estarem sempre presentes e serem um presente na minha vida. À minha mãe, a quem devo todos os meus valores e educação. Por sempre me encorajar a ser melhor. Ao Dinho, pelo amor, dedicação e paciência. Por ser o meu porto seguro e aconchego. Pelo incentivo a me superar todos os dias. Por nunca me deixar desistir. Ao Tomé, que renovou a minha vida! iv 1. STATEMENT OF INTEGRITY I hereby declare having conducted this academic work with integrity. I confirm that I have not used plagiarism or any form of undue use of information or falsification of results along the process leading to its elaboration. I further declare that I have fully acknowledged the Code of Ethical Conduct of the University of Minho. University of Minho, April 15th, 2024 _____________________________________ Sophia Torres Santos v 3. ABSTRACT Design and Optimization of Microbial Communities Microbial communities directly affect surrounding environments and are an important biological process, with potential applications in a variety of fields, such as biotechnology, environmental, and human health. However, the overall understanding of interactions and dynamics in microbial communities remains a challenge. Synergies between computational methods and genome-scale metabolic models have been explored in the last years, as a way to unravel community interactions and behavior, as demonstrated by the numerous simulation methods developed for application in the context of microbial communities. The available simulation methods, with application to microbial communities, were here evaluated and revealed good predictions for phenotypic behavior. However, few studies are available in terms of optimization tools in the community context. Hence, this work describes the implementation of algorithms for the optimization of minimal medium composition, as well as genes/reactions for the production of target compounds. These tools were implemented in MEWpy to transform it into an integrative Python workbench for metabolic engineering to explore constraint-based models of microbial communities. Five hydrothermal samples from the São Miguel Island, Azores, were analyzed to determine prokaryotic community composition to further reconstruct individual and community genome-scale metabolic models, and through simulation and design methods try to unveil possible routes to produce compounds with industrial interest. The first manually curated genome-scale metabolic model for the thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 was developed, uncovering the details of its metabolic capabilities and suggesting for the first time that S. azorense Az-Fu1 may have metabolic potential for bacterial cellulose production. Moreover, the microbial communities of the different samples were modeled, and co-culture optimization was performed using the implemented methods. Among other results, it was shown that S. azorense Az-Fu1 can enhance its cellulose production capabilities when fed with acetate produced by another organism. Keywords: Design, extremophile environments, genome-scale metabolic modeling, microbial communities, optimization. vi 4. RESUMO Design e Optimização de Comunidades Microbianas As comunidades microbianas afetam diretamente os ambientes circundantes e são um processo biológico de enorme relevância, com aplicações potenciais em vários campos, como biotecnologia, meio ambiente e saúde humana. No entanto, a compreensão geral das interações e dinâmicas nas comunidades microbianas continua um desafio. Têm sido exploradas nos últimos anos diversas sinergias entre métodos computacionais e modelos metabólicos à escala genómica como forma de desvendar interações e comportamentos de comunidades, tal como demonstrado pelos inúmeros métodos de simulação desenvolvidos para aplicação neste contexto. Os métodos de simulação disponíveis com aplicação a comunidades microbianas foram avaliados e revelaram boa capacidade preditiva do comportamento fenotípico. No entanto, poucos estudos estão disponíveis em termos de ferramentas de otimização no contexto de comunidades microbianas. Assim, este trabalho descreve a implementação de algoritmos para a otimização da composição mínima do meio, bem como genes/reações para a produção de compostos alvo. Essas ferramentas foram implementadas no MEWpy de forma a transformá-lo em um ambiente de trabalho Python integrado para engenharia metabólica de comunidades microbianas. Cinco amostras hidrotermais da ilha de São Miguel, Açores, foram analisadas para determinar a composição das comunidades procarióticas de forma a reconstruir modelos metabólicos individuais e comunitários à escala genómica e, através de métodos de simulação e design, tentar desvendar possíveis formas de produzir compostos com interesse industrial. Foi desenvolvido o primeiro modelo metabólico à escala genómica manualmente curado para a bactéria termofílica Sulfurihydrogenibium azorense AzFu1, revelando detalhes das suas capacidades metabólicas e sugerindo, pela primeira vez, que S. azorense Az-Fu1 pode ter potencial metabólico para produção de celulose bacteriana. Além disso, foram modeladas as comunidades microbianas das diferentes amostras e a otimização de co-culturas foi realizada usando os métodos implementados. Entre outros resultados, foi demonstrado que S. azorense Az-Fu1 pode aumentar a sua capacidade de produção de celulose quando suplementado com acetato produzido por outro organism. Palavras-chave: Ambientes extremófilos, comunidades microbianas, modelação metabólica à escala genómica, otimização. vii 5. TABLE OF CONTENTS AGRADECIMENTOS .............................................................................................................. iii STATEMENT OF INTEGRITY .................................................................................................... iv ABSTRACT .......................................................................................................................... v RESUMO ........................................................................................................................... vi TABLE OF CONTENTS ......................................................................................................... vii LIST OF FIGURES ................................................................................................................ xi LIST OF TABLES ................................................................................................................ xiv Motivation, Objectives, and Thesis Outline ...................................................... 1 Context and Motivation ................................................................................................ 2 Research Aims ............................................................................................................ 3 Outline of the Thesis .................................................................................................... 4 Scientific Output .......................................................................................................... 5 Chapter 1 Systems Biology in Microbial Communities .................................... 6 1.1 Introduction .................................................................................................... 7 1.2 Systems Biology.............................................................................................. 8 1.3 Genome-Scale Metabolic Modelling ................................................................. 9 1.3.1 Genome-Scale Metabolic Models Reconstruction ..................................... 10 1.3.2 Simulation with Genome-Scale Metabolic Models ..................................... 17 1.3.3 Optimization using Genome-Scale Metabolic Models ................................ 18 1.4 Microbial communities ................................................................................. 19 1.4.1 Historical Perspective ............................................................................. 20 1.4.2 Microbial Community Organisms’ Identification ....................................... 22 1.4.3 Microbial interactions .............................................................................. 26 xiv 7. LIST OF TABLES Table 1.1 Main online data sources used for the reconstruction of genome-scale metabolic models. .................................................................................................................................. 11 Table 1.2. Types of extremophiles based on their habitat. Adapted from (Durvasula et al., 2018). .............................................................................................................................................. 28 Table 2.1 Physical and chemical characteristics of Azorean hydrothermal springs sample sites. .............................................................................................................................................. 39 Table 2.2. Bacterial and Archaeal 16s rRNA sample detection. Low - +, Medium - ++, and High - +++. Samples CV – Caldeira Velha, NP – Nascente Poente, ESG – Esguicho de Maio, FCRG – Fumarola Caldeira da Ribeira Grande, PCRG – Piscina Caldeira da Ribeira Grande. ................................ 43 Table 2.3. Obtained number of reads, contigs, and scaffolds for each sample analyzed. Samples CV – Caldeira Velha, NP – Nascente Poente, ESG – Esguicho de Maio, FCRG – Fumarola Caldeira da Ribeira Grande, PCRG – Piscina Caldeira da Ribeira Grande. ................................................... 44 Table 2.4. Binning algorithm prediction results. Number of total predicted bins, bins with high completeness (>97%) and low contamination (<10%) values, and respective species identification. Samples CV – Caldeira Velha, NP – Nascente Poente, ESG – Esguicho de Maio, FCRG – Fumarola Caldeira da Ribeira Grande, PCRG – Piscina Caldeira da Ribeira Grande. ................................ 45 Table 2.5. Operational Taxonomic Units algorithm results. For each sample only organisms with respective abundance higher than 1% are presented. Shaded organisms belong to the Domain Bacteria, and non-shaded organisms belong to the Domain Archaea. Samples CV – Caldeira Velha, NP – Nascente Poente, ESG – Esguicho de Maio, FCRG – Fumarola Caldeira da Ribeira Grande, PCRG – Piscina Caldeira da Ribeira Grande. Organisms: Acidimicrobium ferrooxidans DSM 10331 ( A. ferrooxidans DSM 10331), Acidithiobacillus caldus SM-1 ( A. caldus SM-1), Desulfurococcus amylolyticus DSM 16532 ( D. amylolyticus DSM 16532), Pyrobaculum aerophilum str. IM2 ( P. aerophilum str. IM2), Sulfurihydrogenibium azorense Az-Fu1 ( S. azorense Az-Fu1), Thermodesulfovibrio yellowstonii DSM 11347 ( T. yellowstonii DSM 11347), Thermofilum adornatus 1505 ( T. adornatus 1505), Thermoplasma acidophilum DSM 1728 ( T. acidophilum DSM 1728), Thermorudis peleae ( T. peleae ), Thermus antranikianii DSM 12462 ( T. antranikianii DSM 12462), Thermus scotoductus SA-01 ( T. scotoductus SA01), Thiomonas intermedia K12 ( T. intermedia K12). .............................................................. 48 xv Table 3.1. Chemolithoautotrophic organisms with published genome-scale metabolic models. CBB - Calvin-Benson-Bassham cycle, rTCA – reverse Tricarboxylic acid cycle, 3-HP/4-HB - 3hydroxypropionate/4-hydroxybutyrate cycle, DC/4-HB - dicarboxylate/4-hydroxybutyrate cycle, WL - Wood-Ljungdahl pathway. ....................................................................................................... 56 Table 3.2. List of phylogenetic similar organisms/genus to S. azorense given to the automatic workflow feature in merlin . ...................................................................................................... 58 Table 3.3. S. azorense genome annotation automatic workflow results taking into account phylogenetically related genus. ............................................................................................... 63 Table 3.4. Biomass macromolecular composition of the S. azorense Az-Fu1 model. .... 66 Table 3.5. S. azorense Az-Fu1 genome information ..................................................... 67 Table 3.6. i SS352 metabolic model information .......................................................... 67 Table 3.7. Minimal medium composition for each condition tested: chemolithoautotrophic and heterotrophic growth. Oxygen and Ferrous iron (highlighted in grey) were only supplied under microaerophilic conditions. ..................................................................................................... 71 Table 3.8. i SS352 model validation against experimental conditions from literature. .... 72 Table 3.9. FVA analysis of i SS352 model cellulose production capabilities. Total consumption of the carbon source was imposed, and the specific growth rate was set to at least 10% of the specific growth rate obtained with pFBA simulation under chemolithoautotrophic conditions. Uptake values are represented as negative and production values are presented as positive values. ..................... 74 Table 3.10. iSS352 prediction of byproduct production under N-limiting conditions and restriction of bicarbonate production. ...................................................................................................... 74 Table 4.1. List of tools used for analyzing microbial community models and their main features. All tools were run under their Python implementation. ............................................................. 90 Table 4.2. Minimal medium composition used to run single organism pFBA simulations using GSMMs of N. europaea and N. vulgaris , respectively. Uptake rates are shown in mmolgDW-1h-1. .. 91 Table 4.3. Minimal medium composition used to run microbial community simulations. Uptake rates are shown in mmolgDW-1h-1. ............................................................................................... 92 Table 4.4. N. europaea model validation using pFBA as simulation method. Specific growth rate (h-1), uptake (Consumption), and export (Production) rates (mmolgDW-1h-1) are presented. ........... 93 xvi Table 4.5. N. vulgaris model validation using pFBA as simulation method, after manual curation of the oxidative phosphorylation pathway. Specific growth rate (h-1), uptake (Consumption), and export (Production) (mmolgDW-1h-1) rates are presented. ........................................................................ 93 Table 4.6. Microbial community, composed of N. europaea and N. vulgaris , simulation using pFBA, SMETANA, SteadyCom, MICOM, and OptCom as simulation methods. When available, specific growth rate (h-1), uptake (Consumption), export (Production) (mmolgDW-1h-1), and interaction rates are presented. Consumed metabolite rates are represented as negative and produced metabolite rates are represented as positive. N.vu. – Nitrobacter vulgaris , N.eu. – Nitrosomonas europaea , Com – Community. n.a. – Data not available. *Abundancy value used as input for the simulation method. ............. 95 Table 5.1 FVA analysis of the co-culture composed by E. coli and S. cerevisiae L-tyrosine production capabilities. The specific growth rate was set to at least 10% of the specific growth rate obtained with the pFBA simulation. ....................................................................................... 112 Table 5.2. pFBA analysis of the co-culture composed by E. coli and S. cerevisiae for L-tyrosine production using the MEWpy deletion and/or overexpression reactions predictions. Consumed metabolite rates are represented as negative and produced metabolite rates are represented as positive. 114 Table 5.3. Microbial community case studies used for the Minimal Medium Optimization and information on the published and manually curated GSMMs used for each community case study. ............................................................................................................................................ 115 Table 6.1. Genome’s information used to import genome files into merlin . ................ 125 Table 6.2. Caldeira Velha microbial community composition and organism physicochemical properties. Organisms’ abundances were recalculated maintaining the relative abundances of the original data. .................................................................................................................................... 133 Table 6.3. Nascente Poente microbial community composition and organism physicochemical properties. Organisms’ abundances were recalculated maintaining the relative abundances of the original data. .................................................................................................................................... 134 Table 6.4. Esguicho de Maio microbial community composition and organism physicochemical properties. Organisms’ abundances were recalculated maintaining the relative abundances of the original data. .................................................................................................................................... 135 Table 6.5. Biomass macromolecular composition according to the organism domain and gram staining. ............................................................................................................................... 136 xvii Table 6.6. Metabolic information of the final Genome-Scale Metabolic Reconstructions. Compartments are divided into c – cytosol, p – periplasm, and e – extracellular. Archaea organisms include a pseudo-compartment to simulate the proton motive force (PMF). ............................ 137 Table 6.7. Minimal medium composition for each condition tested: chemolithoautotrophic and heterotrophic growth. Oxygen and Ferrous iron (highlighted in grey) were only supplied under aerophilic conditions. ........................................................................................................................... 139 Table 6.8. Models’ validation against experimental conditions from literature. ............ 141 Table 6.9. Pairwise interactions predicted by the SMETANA simulation method. Only inter-species interactions with values greater than 0.5 were considered. .................................................... 143 Table 6.10. FVA analysis of S. azorense cellulose production capabilities using a 2 organisms community. Total consumption of the carbon source by S. azorense was imposed, and the specific growth rate was set to at least 10% of the specific growth rate obtained with pFBA simulation under nitrogenlimiting conditions. ............................................................................................................... 147 Table 6.11. Reaction KO analysis of the community formed by S. azorense and E. coli for cellulose production capabilities. Total consumption of the carbon source by S. azorense was imposed, and the specific growth rate was set to at least 10% of the specific growth rate obtained with pFBA simulation under nitrogen-limiting conditions. pFBA simulations were performed under nitrogen-limiting conditions. ACALD - acetaldehyde dehydrogenase (acetylating), HEX7 – hexokinase, PGI - glucose-6-phosphate isomerase, PYK - pyruvate kinase, TPI - triose-phosphate isomerase. ..................................... 149 |1 Introduction Chapter 0 Motivation, Objectives, and Thesis Outline “Our unity is our strength, and diversity is our power” Kamala Harris Microbial communities are widespread in nature and have been recognized to be more adequate in industrial settings than pure cultures, especially regarding the robustness towards contaminations. Moreover, the metabolic capabilities can be greatly extended when compared with individual species. However, very few applications have been described besides the production of food products using natural communities. In order to use fully the industrial potential of microbial communities it is important to be able to customize their behavior by optimizing the composition both in terms of species and genes, as well as optimizing the environmental conditions. This thesis aims to contribute to developing optimized microbial communities for the production of a target compound by improving both modeling and simulation methods and developing metabolic engineering tools. The in silico results will be applied to the optimization of the production of cellulose by Sulfuhydrogenibium azorense Az-Fu in a co-culture system. Introduction 2| PhD Thesis | Sophia Santos Context and Motivation Microbial communities play pivotal roles throughout all environments showing unique biological properties in maintaining Earth’s biosphere (Fierer, 2017; Gilbert et al., 2018; Sunagawa et al., 2015). Microbial communities in environments that once were described as unhabitable (Rothschild et al., 2001) are the ones that nowadays are gaining huge interest due to their remarkable diversity of morphology, biochemistry, genomics, and biosynthesis of many adaptive compounds (Durvasula et al., 2018). Moreover, microbial communities in high-temperature environments are generally less diverse (Inskeep et al., 2013), making hydrothermal habitats an ideal model system for studying principles of community structure and function (Sahm et al., 2013). Genome-scale metabolic models (GSMM) are a relevant example with proven results in accurately predicting metabolic states for an increasing number of single organisms (Feist et al., 2009) as well as for microbial communities (Colarusso et al., 2021), reducing the time and cost implicated in experimental tasks. Together with the growing availability of ‘omics’ data and multiple algorithms to process and make knowledge from metagenomic raw data, these metabolic models are an indispensable systems biology tool to try to understand function, interaction, and dynamics within these microbial consortia (Zaramela et al., 2021). GSMMs are usually used to compute metabolic phenotypes of an organism in response to environmental and genetic perturbations (Feist et al., 2010; Rocha et al., 2008) most commonly through simulation methods, such as Flux Balance Analysis (FBA) (Varma & Palsson, 1994). However, the construction of GSMMs with good phenotypic predictions needs manual curation and experimental validation (Lieven et al., 2020), which is still a laborious task even with available semi-automatic reconstruction pipelines (Arkin et al., 2018; Machado et al., 2018). But more challenging than that is trying to define biologically relevant objective functions (García-Jiménez et al., 2021) in a microbial community context. A myriad of extension methods for FBA with application to community GSMMs is available (Chan et al., 2017; Gomez et al., 2014; Mahadevan et al., 2002; Zelezniak et al., 2015), with proven reliability when predicting growth, response to nutrients, and gene essentiality in single organisms and even microbial communities (Chng et al., 2020; Machado et al., 2021; Nayfach et al., 2020). However, not much research has been made on computational optimization methods to rationally design microbial Introduction PhD Thesis | Sophia Santos |3 communities using GSMMs, which is one of the most promising features for the use of community models for the production of desired compounds in industrially relevant amounts (García-Jiménez et al., 2021; Tsoi et al., 2019; Zaramela et al., 2021). This is particularly relevant given that microbial communities are often easier to use than pure cultures in industrial settings. Research Aims The present thesis aims at reconstructing extremophilic microbial communities' metabolic models and performing simulations with the ultimate objective of using those models to design communities better fitted for industrial aims. Microbial communities are often more adequate than pure cultures for industrial applications and have, at the same time, great potential in terms of metabolic capabilities by combining individual capacities of single strains. The specific objectives are: - Identifying, through metagenomic taxonomic profiling approaches, extremophile organisms on hydrothermal samples from São Miguel, Azores; - Developing genome-scale metabolic models of the extremophile organisms identified to further combine individual models to obtain an accurate representation of the capabilities of a microbial community; - Evaluating existing tools and methods with an emphasis on Flux Balance Analysis based techniques, and using the most promising ones to perform simulations using as a case study the well-established nitrification bioprocess catalyzed by Nitrosomonas euroapaea and Nitrobacter vulgaris . - Adapting metaheuristics algorithms previously developed at the Biosystems group to be applied to the design of improved communities for industrial aims, by allowing the manipulation of the genome of untargeted and target species, and the environmental conditions; - Use the selected simulation tools and methods and metaheuristics algorithms developed to find design strategies for the production of cellulose, using extremophile community models. Introduction 4| PhD Thesis | Sophia Santos Outline of the Thesis To address the above-mentioned objectives, this thesis has been structured into 7 Chapters: ▪ In Chapter 1 a comprehensive review of systems biology approaches on microbial communities was conducted, trying to give an overall view of the metagenomic analysis process and how systems biology using genome-scale metabolic reconstruction, simulation, and optimization can improve knowledge on function, interaction, and dynamics within these microbial consortia. ▪ In Chapter 2 a metagenomic study of five hydrothermal samples from São Miguel, Azores, was performed to determine their prokaryotic community compositions. Taxonomic profiling using assembly-based and read-based analysis was combined and the most abundant organisms were predicted for each sample. Sulfurihydrogenibium azorense Az-Fu1 was one of the most abundant predicted organisms found in two of the hydrothermal samples. ▪ In Chapter 3 the first manually curated genome-scale metabolic model for the thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 is presented. The reconstruction revealed the presence of the main components of the bacterial cellulose operon and its regulators, suggesting that Sulfurihydrogenibium azorense Az-Fu1 may have metabolic potential for cellulose production. ▪ In Chapter 4 a systematic evaluation of different steady-state simulation methods applied to microbial communities has been used to model the well-established nitrification bioprocess catalyzed by Nitrosomonas euroapaea and Nitrobacter vulgaris . Methods performances were compared to assess which ones should be used in a specific community-level context. ▪ In Chapter 5 The implementation of metabolic engineering methods in MEWPy to explore constraint-based models of microbial communities allowing the optimization of microbial communities using Evolutionary Algorithms is described. MEWpy also allows the use of the simulation methods for microbial communities included in the REFRAMED library. ▪ In Chapter 6 the microbial communities of the different samples were modeled and both untargeted and targeted co-culture optimization was performed using the methods described in the previous chapter to evaluate Sulfurihydrogenibium azorense Az-Fu1 capabilities of cellulose production. ▪ Finally, in Chapter 7 the main conclusions of this thesis are summarized and some perspectives on future work based on the questions raised throughout this work are presented. Introduction PhD Thesis | Sophia Santos |5 Scientific Output The main scientific outputs of this thesis are listed below and include peer-reviewed publications and oral presentations at international conferences. Peer-reviewed Publications Santos, S., Dias, O., Rocha, I. Genome-scale metabolic model of thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 (in preparation). Santos, S., Pereira, V., Dias, O., Rocha, M., Rocha, I. Optimization of Microbial Communities using MEWpy (in preparation). Santos, S., Dias, O., Rocha, I. Modeling and Design of Microbial Communities from Extremophilic Environments in the Azores (in preparation). van den Berg, N.I., Machado, D., Santos, S. et al. Ecological modelling approaches for predicting emergent properties in microbial communities . Nat Ecol Evol (2022). Oral Presentations S. Santos, S. Correia, I. Rocha. Inferring optimal minimal medium on genome-scale metabolic models using evolutionary algorithms. – Metabolic Pathway Analysis Conference, Riga, Latvia, 12-16 August 2019. PhD Thesis | Sophia Santos |6 Chapter 1 Chapter 1 Systems Biology in Microbial Communities “Coming together is a beginning. Keeping together is progress. Working together is success” Henry Ford Microbial communities are widespread in nature and have been recognized to be more adequate in industrial settings than pure cultures, especially regarding the robustness towards contaminations. Moreover, the metabolic capabilities can be greatly extended when compared with individual species. However, very few applications have been described besides the production of food products using natural communities. It is important to be able to customize the microbial communities’ behavior by optimizing the composition both in terms of species and genes, as well as optimizing the environmental conditions to use fully the industrial potential. This chapter will review systems biology approaches applied to microbial communities giving insights into how can genome-scale metabolic reconstruction, simulation, and optimization improve knowledge of function, interaction, and dynamics within these microbial consortia. Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |13 unique identifiers such as Enzyme Commission (EC) (Barrett, 1997) and Transporter Classification (TC) (Saier et al., 2021) numbers. Genes involved in regulatory or signaling processes are not included in GSMM reconstructions. Public repositories of genomic data have available annotated genomes, such as NCBI and KEGG. However, in specific cases annotation of genes may be incorrect or even missing; therefore, a re-annotation is always recommended. Information from phylogenetically closely related organisms can also be used to improve genome annotation. 1.3.1.2 Assembling the Metabolic Network The assembly of a metabolic network begins with the identification and collection of all reactions present in the organism in study, from biological databases such as KEGG and BiGG. All reactions catalyzed by enzymes encoded by genes annotated with EC numbers in the previous step must be incorporated into the model, as well as spontaneous reactions. Reactions should also be balanced (generic metabolites and/or metabolites without formula must be curated), and reversibility must be confirmed to avoid mispredictions. Although these steps are automated in most reconstruction tools, manual curation based on curated information from literature or databases, such as MetaCyc, BRENDA, and KEGG is highly recommended. Compartmentalization The next step in the reconstruction process is the compartmentation of the reactions, identifying all the organelles in which enzymes can operate. In prokaryotic organisms, compartments are typically limited to extracellular space, cytosol, and periplasm (in gram-negative bacteria). However, for instance, in Fungi, reactions can occur in up to sixteen compartments (Lu et al., 2019), and, for higher eukaryotes, reactions should be differentiated between tissues. Information on the location of reactions is available in the literature, but several bioinformatics tools have been developed to predict enzymatic location from protein sequences (Jiang et al., 2021) Genes, Proteins, and Reactions A high-quality GSMM should include curated Gene-Proetin-Reactions relationships (GPRs), which allows the accurate prediction of the effect of genetic modifications. These associations are usually defined according to databases and literature (Rocha et al., 2008; Thiele et al., 2010). Most GSMM Chapter 2 14| PhD Thesis | Sophia Santos reconstruction tools, such as merlin (Capela et al., 2021), ModelSEED (Seaver et al., 2021), and CarveMe (Machado et al., 2018) do automatically include GPRs rules in the model. Transport Reactions After compartmentalization, reactions that transport metabolites across the inner and outer membranes must be added to the model. Again, literature information and databases contain these reactions. However, tools such as TranSyT (Lagoa et al., 2021) can use genome information, to generate predicted system-specific transport reactions, providing the associated gene-protein-reaction (GPR) rule. 1.3.1.3 Converting the Metabolic Network to a Stoichiometric Model In this step, the metabolic reaction set is converted into a stoichiometric matrix, and other constraints are added to the model. These constraints include the definition of an abstraction of the biomass in the form of a mathematical equation that represents the drain of macromolecules to generate a new unit (gram) of biomass, energetic requirements through the inclusion of equations that represent the depletion of ATP for cell growth and maintenance, and constraints that represent the environmental conditions. Biomass Composition The composition of the biomass should be experimentally determined through chemostat experiments, or alternatively during the log phase of the cell’s growth. However, in the absence of organism-specific experimental information, data from genome information (particularly nucleotides, deoxynucleotides, and amino acids) can be used or adapted from phylogenetically related organisms. The importance of an accurate biomass composition determination has been reported (Santos et al., 2016), showing that even small differences in biomass content coefficients may considerably impact GSMMs’ predictions. Therefore, a well-defined biomass equation is a crucial step of the GSMM reconstruction process. Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |15 For n biomass constituents, the biomass equation can be formulated as: ∑𝑐𝑖𝑋𝑖→𝐵𝑖𝑜𝑚𝑎𝑠𝑠 𝑛 𝑖=1 where 𝑐𝑖 is the coefficient of each macromolecule or building block, 𝑋𝑖, considered in the biomass. The units of all the coefficients are defined in mmol per gram of dry weight (mmol/gDW) and the biomass units are defined per hour (h-1). The formulation of the biomass objective function can be obtained at different levels of detail: basic level (define the macromolecular content on the cell, i.e., protein, RNA, DNA, lipids), intermediate level (calculate the necessary biosynthetic energy), and advanced level (detailing the necessary vitamins, elements, and cofactors) (Feist et al., 2010). If a biomass precursor is not considered in the biomass reaction, synthesis reactions may not be required for growth, as well as associated genes, which play an important role in in silico gene deletion experiments (Thiele et al., 2010). Growth and Maintenance ATP Requirements Growth-associated ATP maintenance (GAM) reaction accounts for the energy required to replicate a cell, including the synthesis of macromolecules (e.g., Proteins, DNA, and RNA). GAM should be experimentally determined through chemostat experiments. When no experimental data is available, GAM can be estimated by calculating the number of phosphate bonds for macromolecular synthesis (Neidhardt et al., 1990). Non-growth associated ATP (NGAM) maintenance reaction represents the ATP requirements to the cell to maintain for example its membrane leakage (Feist et al., 2007), represented by n ATP hydrolysis reaction, which transforms ATP into ADP and phosphate. The following reaction should then be included in the model’s reaction set: 𝑥 𝐴𝑇𝑃+ 𝑥 𝐻2𝑂 → 𝑥 𝐴𝐷𝑃+ 𝑥 𝑃𝑖+𝑥 𝐻+ where x is the number of phosphate bounds and the NGAM value is set as its reaction rate. As for GAM, the value for the NGAM reaction rate should be experimentally determined or can either be found in literature or estimated by fitting the mode results to experimental data (Rocha et al., 2008). Chapter 2 16| PhD Thesis | Sophia Santos Medium Constraints Transport fluxes for nutrients in the medium should be constrained between zero and infinity. These constraints allow replicating the organism’s physiological behavior under limited substrate availability or maximal uptake rates. The limiting substrate must be constrained to a specific rate value while the remaining ones are usually left unconstrained. Nutrients not available in the medium should be constrained to zero (Rocha et al., 2008). 1.3.1.4 Metabolic Model Curation and Validation The model curation step is an iterative process that continues until simulation results match experimental data in the literature. Tools, such as merlin (Dias, Rocha, Ferreira, & Rocha, 2018), provide a graphical user interface that allows to easily perform re-annotation, correct reactions stoichiometric balance and directionality, include and/or exclude a reaction from the model, and ultimately to export the model in a standard and computational friendly format, such as in Systems Biology Markup Language (SMBL) format. Gap-filling The metabolic network must be screened for possible gaps. The presence of gaps can compromise the synthesis of biomass components and therefore the predictability of the model. Literature and databases (KEGG, MetaCyc, BRENDA, for instance) information should be used to assist the gap-filling process. Several tools can automate the network gap-filling (Mendoza et al., 2019). However, manual curation should follow such procedure. This is an iterative process that is repeated until all biomass precursors and other relevant compounds can be synthesized, and a feasible model is obtained. Metabolic Model Consistency Tests and Comparison to Experimental Data Once gap-filling is performed, consistency tests can verify the actual phenotypic prediction accuracy of the GSMM. Firstly, the metabolic model must be able to compute growth under specific conditions and the predicted rates assessed to published characterization studies (including growth, secretion, and uptake rates). Growth under a defined growth medium is strongly advised whenever possible (Feist et al., Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |17 2008). Specific growth conditions can include growth in limiting substrates, aerobic and anaerobic conditions, and different growth conditions reported in the literature. Also, the analysis of active pathways under specific growth conditions can be performed for model validation. If GSMM predictions are not in accordance with the experimental results, the model should be examined, potentially missing reactions included, and incorrect reactions removed. Another approach for validating GSMMs is assessing the simulation results of experiments performed with deletion mutants. This approach can provide valuable insights into the predictive capabilities of the model and a good training set may be of great value for model debugging. After this validation step, the high-quality GSMM should be exported in a standard format, such as Systems Biology Markup Language (SBML) to be used on GSMM simulation platforms. A high-quality standardized test of a GSMM can be accessed using the MEMOTE test suit (Lieven et al., 2020), which tests fundamental SBML semantic and conceptual requirements. A variety of platforms for simulation with GSMMs are available, such as the user-friendly java tool Optflux (Rocha et al., 2010), Matlab® through the COBRA toolbox (Heirendt et al., 2019)), and Python through COBRApy (Ebrahim et al., 2013) or Reframed packages. 1.3.2 Simulation with Genome-Scale Metabolic Models The most common GSMM simulation approach, Flux Balance Analysis (FBA) (Varma & Palsson, 1994), calculates intracellular reaction rates (also known as flux distributions), under a steady-state assumption (no accumulation of internal metabolites), along with nutrient uptake and secretion, and the cellular growth rate (Orth et al., 2010). Since the mathematical representation of a GSMM forms an underdetermined system of linear equations, the determination of the most biologically plausible solution requires the specification of an objective function. The metabolic network can be represented by a stoichiometric matrix 𝑆 , of dimensions 𝑚 × 𝑛, where 𝑚 corresponds to the total number of metabolites and 𝑛 to the total number of reactions in the network. The coefficients of the 𝑆 matrix define the relationship between the reactions and compounds of the metabolic network. An optimal solution, consistent with the known constraints and the maximization or minimization of a given objective function ( 𝑍 ), can be obtained by solving the linear problem: Chapter 2 18| PhD Thesis | Sophia Santos 𝑀𝑎𝑥𝑖𝑚𝑖𝑧𝑒 𝑍 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑆𝑣 =0 ∝𝑗≤𝑣𝑗≤𝛽𝑗, 𝑗 = 1,…𝑁 where 𝑣 is a vector of fluxes of each reaction, and ∝𝑗 and 𝛽𝑗 are the lower and upper limits for the fluxes, respectively. These limits are used to set the reversibility of the reactions, to limit uptake and secretion rates, and eventually to specify measured fluxes (Rocha et al., 2008). The solution to a FBA problem with respect to a given objective function is a set of optimal metabolic flux distributions. Although the objective value found is unique, flux distributions rarely are. To tackle this problem a second optimization criterion can be added such as found in the FBA variants of Parsimonious enzyme usage Flux Balance Analysis (pFBA) (Lewis et al., 2010), which attempts to find a more biologically viable flux distribution, by trying to minimize the absolute value of the sum of all fluxes through all reaction in the network while maintaining the reached optimum or Flux Variability Analysis (FVA) (Mahadevan et al., 2003), which sets a constraint that requires the objective flux to be equal to its optimal value, and assess the robustness of a flux distribution regarding its production capability of the target compounds. Robustness is analyzed by maximizing and minimizing each reaction flux. High robustness is verified when the predicted maximum and minimum flux values differ slightly. The major challenge in FBA is the definition of an objective function, with biological relevance (Gianchandani et al., 2010). A variety of objective functions have been used to define an FBA problem. However, in the case of microbes, the maximization of growth rate is the most commonly used assumption based on selection pressure (Feist et al., 2010). Other objective functions can be used, such as minimizing ATP production or maximizing the desired compound. FBA has a wide range of applications, such as the optimization of bio-processes in industries or the identification of drug targets (Raman et al., 2009). 1.3.3 Optimization using Genome-Scale Metabolic Models Genetic manipulation of microbes for industrial purposes is nowadays widely used to produce biobased, sustainable, environmentally friendly, and viable compounds (Stephanopoulos, 1999). GSMMs are in that matter a commonly used tool for rational strain design by biotechnology companies with Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |19 successful applications (Jullesson et al., 2015) aided by the development of computational strain optimization methods. These methods try to identify a set of possible genetic modifications towards a desired phenotypical trait. Typically, the objective is to optimize the production of the desired compound, maintaining cell viability. Based on this bi-level objective, OptKnock (Burgard et al., 2003) established the foundation for the development of many other computational strain optimization methods. A variety of computational strain design optimization methods have been developed to search for non-intuitive genetic designs in more efficient and scalable ways. Today over 30 different computational strain optimization methods are published and can be divided into 3 main branches: bi-level mixed-integer programming, metaheuristics, and elementary-mode analysis-based methods. These methods can use different options for optimization targets (genes, reactions) or tasks (deletions, insertions, up/down-regulation) (Maia et al., 2016). Most of the mixed-integer programming and elementary-mode analysis methods guarantee to reach a global optimal solution. However, such approaches do not scale well with larger models or a high number of perturbations. Alternatively, metaheuristic-based methods, such as Evolutionary Algorithms and Simulated Annealing, not only usually have good scalability but are more flexible to the specification of objective functions (Rocha et al., 2008). However, these methods do not guarantee that a globally optimal solution is found. Accessing and using these methods is not always straightforward as the implementation of some methods is not available (Maia et al., 2016). Platforms such as Opflux (Rocha et al., 2010), and CAMEO (Cardoso et al., 2018) incorporate metaheuristic (OptGene (Patil et al., 2005) family methods) and deterministic methods (OptKnock) (Burgard et al., 2003), but are currently restricted to the use of GSMM containing only metabolic information. The recently developed python package MEWpy (Pereira et al., 2021), offers a practical interface to several optimization heuristics allowing to model and optimize microbial production on GSMMs with metabolic, transcriptional, and translational information. 1.4 Microbial communities Microbial communities are immensely abundant throughout all environments and their composition is equally diverse. Microorganisms within their community share not only the same ecosystem but also the same resources. From extreme environments to host symbionts, these communities play pivotal roles in a multitude of processes and have unique biological properties Chapter 2 20| PhD Thesis | Sophia Santos maintaining Earth’s biosphere and contributing to plant and animal physiology and health (Fierer, 2017; Gilbert et al., 2018; Sunagawa et al., 2015). Thus, the possibility of controlling and engineering natural and synthetic microbial communities would be of huge interest to the environmental, health, and industrial processes (García-Jiménez et al., 2021). 1.4.1 Historical Perspective The first attempt to “measure the invisible” microorganisms was made by Antoine Van Leeuwenhoek in 1676, using an indigenous microscope to observe the oral microbiome (Hamarneh, 1960). But only in the late 1800s, Robert Koch was able to isolate and cultivate microorganisms in a solid phase. His studies not only helped to understand microbial physiology but also to establish causality between microorganisms and disease (Blevins & Bronze, 2010). Since then, and until the emergence of the microbial ecology concept, based on Winogradsky's works in the 1920s and 1930s to foment the growth of symbiotic microbes (Dworkin & Gutnick, 2012), studies were focused on pure cultures. In fact, the study of microorganisms and their environmental roles within communities, initiated by Winogradsky, revolutionized microbiology and these culture-based approaches were widely applied to microbial communities from all environments (Ackert, 2013). Identification of microorganisms was the next step in microbiological history, and indeed a vast number of methods were developed over time, such as the gram-staining technique, the agar selective and differential media, and biochemical tests. However, these culture-dependent methods besides being time-consuming and labor-intensive are also inadequate for identifying phenotypically similar species (Buszewski et al., 2017). Perhaps the breakthroughs in these almost 400 years of history (Figure 1.2) were the works of Carl Woese, Fred Sanger, and Kary Mullis. The emergence of 16s rRNA subunits as molecular taxonomy markers (Woese et al., 1977) and the development of automated DNA sequencing methods (Sanger et al., 1977) and techniques, such as polymerase chain reaction (PCR) (Mullis et al., 1986) allowed not only the identification and classification of microorganisms but also to unveiled microbial communities unculturable species. Supported by the development of these techniques, metagenomics as a research field surfaced at the beginning of the 1990s, with the first description of a method for the identification of species in a natural microbial community using 16s rRNA libraries (Ward et al., 1990). However, the term Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |21 “Metagenomics” was only proposed later that decade, disclosing all the industrial potential of accessing genomes of unculturable organisms (Handelsman et al., 1998). Although the amount of data that was generated by these techniques and instruments was already vast, limited throughput and the high costs of sequencing were still barriers identified by the Human Genome Project development (Goodwin et al., 2016). Next-generation sequencing, presented in the mid-2000s, is a truly high-throughput sequencing platform, greatly reducing the necessary reaction volume while dramatically extending the number of sequencing reactions, and reducing the time, costs, and complexity required to sequence large amounts of DNA (Schuster, 2007). Figure 1.2. Metagenomics timeline and milestones. Timeline showing advances in microbial communities’ studies. Adapted from (Escobar-Zepeda et al., 2015) These new achievements contributed to the exponential growth in the size of sequencing repositories and launched metagenomics as a new scientific field. Since then, a global effort has been Chapter 2 22| PhD Thesis | Sophia Santos made to catalog the uncultured microbial diversity through the Global Ocean Sampling Expedition (Rusch et al., 2007), the Human Microbiome Project (Turnbaugh et al., 2007), and more recently the Earth Microbiome Project (Gilbert et al., 2014). Metagenomics has reshaped our view of the tree of life (Hug et al., 2016; Parks et al., 2018), led to the identification of deeply rooted and metabolically diverse organisms (Jaffe et al., 2020; Wilson et al., 2014), and leveraged the characterization of new biosynthetic pathways and products (Crits-Christoph et al., 2018; Freeman et al., 2017). 1.4.2 Microbial Community Organisms’ Identification One of the primary objectives of studying microbial communities is determining their composition at the species or strain level, ideally with quantitative information on occurrences (Frioux et al., 2020). The emergence of high-throughput sequencing technologies and new computational pipelines, combined with shotgun metagenomics methods allowed exploring genomic information of all organisms in a microbial community, not only the culturable ones (Zorrilla et al., 2021). In fact, metagenomics has enabled large-scale investigation of complex microbiomes and unveiled species that carry out complete nitrification of ammonia (Daims et al., 2015), or the widespread presence of antibiotic genes in commensal gut bacteria (Donia et al., 2014), for instance. However, the process of identifying microorganisms has some limitations due to the yet large amount of unculturable organisms without a reference genome available (Zorrilla et al., 2021). The identification of organisms and extraction of microbial abundances from the raw data is usually a complex and multiple-step procedure. The process can be summarized in four steps: (i) Experimental pipeline; (ii) Pre-processing analysis; (iii) Sequence analysis and (iv) Post-processing analysis (Quince et al., 2017). For each of these steps, numerous experimental and computational approaches are available, each of which with limitations and challenges. A short description of each step is available below. 1.4.2.1 Experimental pipeline The experimental pipeline consists of collecting, processing, and sequencing the metagenome samples. The importance of this step is often underestimated in metagenomics. Sample collection and preservation protocols are extremely important steps once any miscalculation can affect both the quality and the accuracy of metagenomics data. Moreover, careful preliminary work is often required to ensure Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |29 evolve, including structural and biochemical adaptations, which allowed growth and survival under conditions that imposed a lack of nutrients and energy (Ando et al., 2021). Often these adaptations are a consequence of key changes in the organism’s enzymes amino acid sequences, which are translated into variations in the structure, flexibility, charge, and/or hydrophobicity (Sarmiento et al., 2015). Therefore, extremozymes can offer exciting industrial opportunities when compared with organic synthesis and even mesophilic biocatalysts, allying better chemical precision, sustainability, and costeffectiveness, to harsh enzymatic conditions capable of withstanding industrial demands (Van den Burg, 2003). Indeed, extremophiles' enzymes that have been isolated and functionally characterized, led to the identification of the thermostable DNA polymerase used in the polymerase chain reaction (PCR) (Ishino et al., 2014) and enzymes used in the biofuel (Barnard et al., 2010), food, paper and cosmetic industries (Sarmiento et al., 2015), for instance. Until now, hyperthermophilic and thermophilic extremophiles/extremozymes have attracted the most attention. Prokaryotic diversity in hydrothermal ecosystems has been extensively studied (Ando et al., 2021). Moreover, the temperature endured by thermophilic enzymes, besides increasing reaction rates during high-temperature processing, reduces microbial contamination, lowers substrate viscosity, and increases the solubility of many polymeric substrates (Van den Burg, 2003). 1.4.5 Genome-scale Metabolic Modelling of Microbial Communities The increased interest in studying microbial communities due to their key roles in human and environmental health was highly supported by the recent advances in high-throughput multi-omic technologies (metagenomics, metatranscriptomics, metaproteomics) (Zaramela et al., 2021). However, the use of these techniques to understand the individual role of microbes and their interactions, with each other and the community, is still a challenge. The use of GSMMs through computational and mathematical modeling approaches is nowadays a widely used platform to generate testable hypotheses on microbial community behavior at the taxon and community levels (Colarusso et al., 2021). Here a brief description of the current existent methods and tools for reconstruction, simulation, and optimization of microbial communities using GSMMs will be performed, including an assessment and a discussion of current limitations and challenges. Chapter 2 30| PhD Thesis | Sophia Santos 1.4.5.1 Microbial Community Genome-Scale Metabolic Model Reconstruction Microbial community GSMMs’ reconstruction is based on the reconstruction process of a single species. Two main approaches are used for these reconstructions: the mixed-bag and the compartmentalized network approaches (Faria et al., 2017). The mixed-bag network approach considers the microbial community a single supra-organism (Abubucker et al., 2012), constructing a model composed of one cytosolic and one extracellular compartment, ignoring species boundaries. In this case, the resulting model will be analogous to a singlespecies prokaryotic network. Species-level resolution is ignored, and these networks are usually used to analyze community—environment interactions (Henry et al., 2016). The reconstruction of a mixed-bag network requires either the annotated genome sequences of all species present in the community or the whole deep metagenome sequence with annotated reads. From this stage, the network reconstruction is performed similarly to a single species network (see section 1.3.1). In the final model, redundant reactions (reactions included in at least two of the community species) are ignored. The compartmentalized network approach considers the different species of the community in distinct compartments. For each organism in the community, a single species metabolic model is constructed based on its genome sequence. The community model will be composed of as many compartments as the number of organisms in the community, plus an extracellular pool to allow metabolite exchange. Species-level resolution is obtained and therefore cross-species metabolic interactions can be predicted and optimized. Metabolic capabilities can be associated with the respective species within the community, which will allow microbial community design (Faust, 2018). However, in this approach, microbial community composition in terms of species must be known. Hence, reconstructing models from metagenomic samples, with undefined composition, demands an initial step of species identification from the metagenome(see section 1.4.2). The selection of an approach to reconstruct a microbial community network will depend on data and resource availability, as well as on the purpose. The limitations and challenges of each approach will be detailed in section 1.4.5.4. Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |31 1.4.5.2 Simulation of Microbial Communities using Genome-Scale Metabolic Models The application of GSMMs to the simulation of microbial communities was pioneered by Stolyar et al, who have applied them in the simulation of a two-species methanogenic community (Stolyar et al., 2007). Several other simulation methods have since been published, but a comprehensive and systematic analysis will be detailed in Chapter 4 of this thesis. Here, pioneering methods for microbial community simulation will be discussed. As demonstrated before, multiple GSMMs can be easily merged into a compartmentalized multi-species GSMM, where organisms are allowed to exchange compounds through a common extracellular pool. However, formulating a biologically and ecologically meaningful objective function for a microbial community model is not trivial, and is one of the main differentiating features of different simulation methods. Two other main differentiating features are the inclusion of temporal and spatial components. The simplest class of methods is used for the simulation of steady-state flux distributions in microbial communities, allowing the prediction of individual growth rates and cross-feeding interactions in community equilibrium. OptCom (Zomorrodi et al., 2012) implements a bi-level objective function where the inner objective is to maximize the growth rate of individual species, and the outer objective is to maximize the total community biomass. The underlying assumption is that the community will optimize the management of available resources while still considering the individualistic objectives of its members. Other methods, such as cFBA (Khandelwal et al., 2013), SteadyCom (Chan et al., 2017), and MICOM (Diener et al., 2020), implement different objective functions following similar assumptions. One advantage of steady-state methods is that they do not require parameterization. However, they are characterized by a large uncertainty in the space of optimal solutions. The second class of methods considers the temporal component when performing time-course simulations. These are generally implemented as multi-species extensions of dynamic FBA (dFBA) (Zhuang et al., 2011), where time is divided into discrete steps and an FBA simulation is performed for each species at each time step. The d-OptCom (Zomorrodi et al., 2014) method combines dFBA with OptCom to implement dynamic multi-objective optimization. Recently, a faster implementation of multispecies dFBA, µbialSim (Popp et al., 2020), was able to scale up to communities with hundreds of species. These methods account for the temporal variation in metabolite concentrations and species abundance, which can play a role in community assembly, due to sequential species colonization, and allow simulating the response to changes in initial conditions and other environmental perturbations. Chapter 2 32| PhD Thesis | Sophia Santos However, every organism requires the characterization of the substrate’s uptake kinetics for all metabolites that are consumed (either from the growth medium or through cross-feeding). The impractical nature of such a massive in vitro characterization results in the adoption of default parameter values for all species and compounds, which limits the predictive ability of these methods. The third class of methods adds the spatial component as well, resulting in spatiotemporal simulations. These methods are especially relevant in situations of heterogeneous or spatially segregated environments, such as biofilm colonies or the tract along the human gut, where nutrient diffusion and access to resources play a role in community assembly and the potential for cross-feeding interactions. The implementation of spatiotemporal methods mostly differs in the way the spatial component is modeled. COMETS (Harcombe et al., 2014) extends dFBA with a 2D-grid where each square represents a population of cells at a given point in time. BacArena (Bauer et al., 2017) also uses a 2D-grid structure but follows an individual-based modeling (IBM) approach, where each cell represents one individual, allowing for a more fine-grained resolution of inter and intra-species interactions. IndiMeSH (Borer et al., 2019) also uses an IBM approach, but the spatial structure is based on angular pore networks, which allows the creation of more complex structured environments. In addition to the parameters on uptake kinetics mentioned before, these methods also require the specification of diffusion rates for cells and metabolites. These three main classes of methods represent increasing levels of simulation detail but also require more experimental data for model setup and, naturally, more computational power. Therefore, the selection of the most suitable kind of method depends on the complexity of the microbial community under study, as well as the amount of available experimental data. Regardless of this choice, GSMMbased simulation allows for a mechanistic interpretation of observed community phenotypes and of the competitive and cooperative behaviors contained therein. Multi-omics data (transcriptomics, proteomics, metabolomics) can be used to constrain the solution space of metabolic models further highlighting the metabolic pathways active in each condition, and analyzing how different kinds of perturbations will elicit a coordinated response comprising the molecular, cellular, and population levels. 1.4.5.3 Microbial Community Genome-Scale Metabolic Model Optimization The biological properties of microbial communities are closely related to the organisms present within them and resource availability (Johns et al., 2016). If the community composition (species or Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |33 medium) is disturbed, its functional capabilities can be affected. Such findings challenged the scientific community to manipulate microbiomes to improve crop productivity (O’Connell et al., 1996), to deal with contaminated groundwater (Löffler et al., 2006), or even to recover valuable resources from wastewater (McCarty et al., 2011), showing the huge potential of engineering microbial communities. Two main approaches have been currently used for engineering microbial communities: the top-down approach (Widder et al., 2016), which aims to control metabolic processes for stabilizing complex and natural communities; and the bottom-up approach (Lindemann et al., 2016), which aims at designing defined, synthetic communities with desired functionalities (Zerfaß et al., 2018). When combined, these approaches can offer complementary strategies to design defined, synthetic communities used to impact and engineer the behavior of complex communities (Lawson et al., 2019). Given the success of computational strain optimization methods to rationally design single organisms using GSMMs, already addressed in section 1.3.3, it is reasonable to assume that the same strategies may lead to success using microbial communities. Microbes rarely live in isolation, thus engineering microbial communities is a strategy closer to natural and physiologic behavior, than for isolated organisms. Moreover, the engineered product will be a cooperative effort from a structured community able to carry out complex interactions and synthesize complex molecules (Bosi et al., 2017; García-Jiménez et al., 2021). Until now, the number of studies that aim at developing algorithms for designing microbial communities using GSMMs has been low. Most of the tools developed use strategies to optimize medium composition and/or best community configuration for a given objective. Among them, the SIM algorithm (Klitgord et al., 2010), developed over the MATLAB® platform, uses a mixed-integer linear programming approach, to optimize a medium composition that supports the growth of multispecies co-cultures in a specific condition. An initial minimal set of potential exchanged metabolites that allows a positive growth rate of both organisms is identified. The core function of the SIM algorithm identifies all potential metabolites that can be switched in such a medium and the initially available source for a given element (carbon, nitrogen, sulfur, or phosphate). Hence, the chemical formula of each metabolite present in the GSMMs must be available. The algorithm is based on the identification of a symbiosis-inducing media, assuming that each species will secrete what the other needs. Similarly, FLYCOP (García-Jiménez et al., 2018), supports the design and engineering of microbial communities by selecting a consortium configuration that optimizes a given goal (community growth rate, Chapter 2 34| PhD Thesis | Sophia Santos stability, medium composition, etc). Hence, FLYCOP selects and evaluates candidate consortium configurations, through an iterated local search algorithm (SMAC (Hutter et al., 2011)). For strain design, individual strain engineering must be performed beforehand, optimizing only the conditions for fermentation. Using a different approach, DOLMN (Thommes et al., 2019) explores the space of feasible singlestrain or multistrain metabolic networks, by systematically limiting the number of intracellular and transport reactions (performing reaction knockouts) in each metabolic model, using a MILP formulation. However, optimal solutions result in deleting a large number of reactions, which is not applicable for practical purposes. The increasing interest in designing microbial communities demands that model-guided microbial community engineering should trend toward the development of technologies capable of predicting potential genetic modifications at the community level, like individual-level design platforms, such as OptKnock. These new implementations of microbial community design methods may focus not only on optimizing exchange reactions (species cross-feeding interactions, medium composition) but also on untargeted and targeted (directed to specific species in the community) genes/reactions, as well as identifying the best community composition for a given objective, using top-down and bottom-up approaches. 1.4.5.4 Genome-Scale Metabolic Modelling of Microbial Communities - Limitations and Challenges Although microbial community modeling is widely established, there are known limitations and challenges. A good GSMM must be manually curated and validated with experimental data to accurately predict phenotypic behavior. The process of reconstructing a community GSMM can be similar to that of a single metabolic model. In this approach, one would start by annotating all genes in a metagenomics sample (mixed bag approach). Alternatively, one could initially perform the identification of all organisms present in the community and then reconstruct the corresponding GSMMs (compartmentalized approach) (Frioux et al., 2020). In the latter approach, depending on the number of organisms and the availability of curated GSMMs in databases, such as the BiGG knowledgebase (King et al., 2016a), the time it takes to process, collect, and eventually reconstruct GSMMs can increase exponentially. Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |35 In the last years, pipelines that semi-automatically reconstruct prokaryotic organisms’ GSMMs in a matter of minutes, such as KBase (Arkin et al., 2018), the python package CarveMe (Machado et al., 2018) or the AGORA knowledgebase which is focused in studying the human microbiota or host-pathogen interactions (Almut Heinken et al., 2020), have been published. These models have already shown their reliability in predicting growth, response to nutrients, and gene essentiality in single organisms and microbial communities (Chng et al., 2020; Machado et al., 2021; Nayfach et al., 2020). However, when the objective is the rational design of microbial communities, these models should be used with caution, as they often have some inaccuracies (Lieven et al., 2020), requiring further manual curation. To accurately study a microbial community, it is essential to correctly identify the available species. Besides wet-lab limitations, the lack of a complete reference genome and missing functional annotation of microbial genes in current databases (Bharti et al., 2021; Quince et al., 2017) impair the reconstruction of community GSMM. When reconstructing microbial community’s GSMMs, the struggle starts with how to compare simulation results. Several methods for simulating microbial communities using GSMMs have been published lately, although the full potential of those methods has not been achieved. The current experimental techniques are designed for individual organisms, with still limited applicability to microbial communities, as not enough experimental data on abiotic and biotic interactions, and perturbations is yet available (Machado et al., 2018). Although one of the main objectives of studying microbial communities is the possibility of controlling and engineering natural and synthetic microbial communities, few studies are available in terms of optimization tools capable of predicting potential genetic modifications at the community level (García-Jiménez et al., 2021). The challenge here is the development of sophisticated and integrative platforms that support different levels of community optimization (e.g. medium, interactions, community composition, targeted organism design within the community) to unveil the full potential of microbial communities (Eng et al., 2019). PhD Thesis | Sophia Santos |36 Chapter 2 Chapter 2 Metagenomic study of thermophilic and hyperthermophilic environments from Azores “You are never alone. You are externally connected with everyone” Amit Ray Five hydrothermal samples from hydrothermal vents at São Miguel, Azores, were analyzed to determine prokaryotic community composition. Taxonomic profiling of all samples was performed using assembly-based and read-based analysis showing combined results specifically for the predicted most abundant organisms in each sample. Samples showed to be abundant in members of the Aquificales and Crenarchaeota orders, in specific, Sulfurihydrogenibium azorense Az-Fu1, Desulfurococcus amylolyticus DSM 16532, Pyrobaculum islandicum DSM 4184/ Pyrobaculum aerophilum str. IM2 and Thermofilum adornatus 1505. S. azorense Az-Fu1 was first isolated at São Miguel, Azores but in a different location (Furnas). Thus, the metabolic capabilities of S. azorense Az-Fu1 in isolation and within its microbial community, hence a systems biology approach, is worthy to be investigated. Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |37 Sophia Santos1, Ana Alão Freitas2, Ricardo B. Leite3, Duarte N. Toubarro4, José B. Leal5, Nelson Simões4, Oscar Dias1, Isabel Rocha2 1 Centre of Biological Engineering of University of Minho, Braga, Portugal 2 Instituto de Tecnologia Química e Biológica António Xavier, Universidade Nova de Lisboa (ITQB-NOVA), 3 Instituto Gulbenkian de Ciência, Oeiras, Portugal 4 Centro de Biotecnologia dos Açores, Departamento de Biologia, Universidade dos Açores, Ponta Delgada, Portugal 5 Ophiomics, Lisboa, Portugal Autors’ contributions Sophia Santos conceived the study, carried out the read-based taxonomic profiling, analyzed the results, and drafted this chapter. Duarte N. Toubarro carried out the sample collection and DNA extraction. Ana Alão Freitas and Ricardo B. Leite carried out the metagenomics and bioinformatics pipeline. Nelson Simões and José B. Pereira Leal participated in the design of the study. Oscar Dias and Isabel Rocha conceived the study participated in its design and coordination and helped to draft this chapter. Chapter 2 38| PhD Thesis | Sophia Santos 2.1 Introduction From extreme environments to host-symbionts, microbial communities play pivotal roles in a multitude of processes and have unique biological properties maintaining Earth’s biosphere, contributing to plant and animal physiology, as well as human health (Fierer, 2017; Gilbert et al., 2018; Sunagawa et al., 2015). The appearance of Metagenomics as a field of study has transformed the way these microbial communities are perceived, being nowadays possible to control and engineer natural and synthetic microbial communities, with applications on environmental, health, and industrial processes (GarcíaJiménez et al., 2021). The study of microbial communities in environments that once were described as unhabitable (Rothschild et al., 2001) – extremophiles – which are capable of sustaining biological life under one or more environmental stressors and are known for their remarkable diversity of morphology, biochemistry, genomics, and biosynthesis of many adaptive compounds (Durvasula et al., 2018) has nowadays gained huge importance. In harsh environmental conditions, microorganisms evolved several structural and biochemical adaptations that allowed them to grow and survive under conditions that impose, among others, a lack of nutrients and energy (Ando et al., 2021; Sarmiento et al., 2015). Therefore, the study of extreme microorganisms through their outstanding metabolic capabilities or their extremozymes can offer exciting industrial opportunities (Van den Burg, 2003). When studying microbiomes via metagenomics, one of the primary objectives is the identification of the microorganisms, ideally with quantitative information on their occurrences (Frioux et al., 2020). Despite the existence of multiple algorithms, the identification of organisms and extraction of microbial abundances from the raw data is still a complex and challenging procedure. However, microbial communities in high-temperature environments are generally less diverse (Inskeep et al., 2013) making hydrothermal habitats an ideal model system to study principles of community structure and function (Sahm et al., 2013). In fact, prokaryotic diversity in hydrothermal ecosystems has been extensively studied at hydrothermal sites, such as Yellowstone National Park (Inskeep et al., 2013; Jay et al., 2016), repeatedly revealing characteristic taxonomic groups, like the main presence of Aquificales and Crenarcheota members (Strazzulli et al., 2017). Being widespread in nature, in marine and terrestrial geothermal systems, members of both phyla are essentially thermophilic and hyperthermophilic microorganisms, almost exclusively obligate or facultative autotrophs (Hedlund et al., 2015) and are also Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |45 In fact, all samples appear to have high prokaryotic diversity, as shown by the total number of bins predicted for each sample. However, Kraken results are hugely dependent on the availability of marker genes and reference genomes in databases, making it difficult to profile organisms present in low quantities in samples, as well as uncharacterized organisms that are typically present in the soil, marine, and freshwater environments (Breitwieser et al., 2019), such as the ones being analyzed in this work. Table 2.4. Binning algorithm prediction results. Number of total predicted bins, bins with high completeness (>97%) and low contamination (<10%) values, and respective species identification. Samples CV – Caldeira Velha, NP – Nascente Poente, ESG – Esguicho de Maio, FCRG – Fumarola Caldeira da Ribeira Grande, PCRG – Piscina Caldeira da Ribeira Grande. Sample Number of bins Species Identification Domain Total High Completeness Low Contamination CV 17 2 Sulfurihydrogenibium azorense Az-Fu1 Thiomonas intermedia K12 Bacteria Bacteria NP 12 2 Pyrobaculum islandicum DSM 4184 Thermus scotoductus SA-01 Archaea Bacteria ESG 10 2 Sulfurihydrogenibium azorense Az-Fu1 Pyrobaculum islandicum DSM 4184 Bacteria Archaea FCRG 11 1 Thermoplasma acidophilum DSM 1728 Archaea PCRG 10 2 Acidithiobacillus caldus SM-1 Acidimicrobium ferrooxidans DSM 10331 Bacteria Bacteria As mentioned before, only a small number of genomes could be identified using this approach. Given the methodology of this taxonomic profiling approach, the identified genomes have a fair probability of being the most abundant in each sample (Hug et al., 2016). Indeed, all samples except sample FCRG, which is also one of the samples that revealed higher archaeal diversity, show bacterial diversity in accordance with the results obtained by 16S rRNA gene amplicon sequencing (section 2.3.1). Binning was not able to predict any species identification of archaeal diversity in CV and PCRG samples, although species identification of around 80% of the predicted bins remains inconclusive in both samples. Moreover, the results are also in accordance with previous studies that report the high abundance of Aquificales ( Sulfurihydrogenibium azorense Az-Fu1) and Crenarchaeota members ( Pyrobaculum islandicum DSM 4184) in hydrothermal sites (Strazzulli et al., 2017). Chapter 2 46| PhD Thesis | Sophia Santos 2.3.3.2 Read-based analysis - Operational taxonomic units (OTUs) Read-base taxonomic profiling of metagenomes more than identifying microbial species in a metagenome also estimates their abundance. Based on simply mapping reads to reference genomes, these methods are ideal for profiling low-abundance organisms (Nielsen et al., 2014). OTUs generated with the mOTUs profiler v.2.0.0 (Milanese et al., 2019) for each sample are presented in Table 2.5. Abundance percentage values presented were normalized for the species with an abundance percentage higher than 1%. Complete mOTUs results are presented in Supplementary Table S6. Accounting for all samples, a total of 12 organisms were identified. Figure 2.2. Domain-level composition of microbial communities based on mOTUs algorithm organism identification and respective abundance prediction. The obtained results agree with the ones obtained in section 2.3.3.1, for both prokaryotic diversity and organism identification. As predicted by assembly-based analysis, all samples seem to have bacterial and archaeal diversity, exception made again for sample PCRG that is composed only of bacteria (Figure 2.2). Although mOTUs algorithm only predicted the presence of two organisms in the PCRG sample, the Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |47 binning algorithm predicted 10 bins, that can or not be mapped to actual different genomes (Hug et al., 2016). The identification procedure in the taxonomic profile process of both approaches is highly dependent on reference genomes in databases, so archaeal diversity detected using 16S rRNA gene amplicon sequencing (section 2.3.1) could be resultant of unknown archaeal organisms (Zorrilla et al., 2021). CV reveals to be the sample with higher prokaryotic diversity (composed of 9 different microorganisms), showing higher bacterial diversity than archaeal, which is in accordance with all previous results. Considering organism identification, all samples show that the most abundant organisms identified by the read-based approach are also the ones identified by the binning approach. However, Thiomonas intermedia K12 and Thermus scotoductus SA-01 in samples CV and NP, respectively, were predicted by both approaches but are not included in the most abundant by the mOTUs algorithm. Moreover, results for both approaches exhibit remarkable similarity for samples with lower prokaryotic diversity (ESG, FCRG, and PCRG), showing that indeed read-based and assembly-based analysis approaches validate and complement each other. A unique divergence in the results has been spotted: the two approaches have identified in samples CV and ESG a different member of the Pyrobaculum genus, Pyrobaculum islandicum DSM 4184 by the binning approach and Pyrobaculum aerophilum str. IM2 by the read-based approach. In fact, the two organisms are genetically and metabolically similar (Feinberg et al., 2008) so these differences can be due to the use of different reference genome databases (mOTUs uses an internal reference genomes database while BinSanity uses NCBI RefSeq database) or even to incorrect fragment assignment which can prevent achieving accurate contigs of high-abundant genomes, before organism identification (Breitwieser et al., 2019). Chapter 2 48| PhD Thesis | Sophia Santos Table 2.5. Operational Taxonomic Units algorithm results. For each sample only organisms with respective abundance higher than 1% are presented. Shaded organisms belong to the Domain Bacteria, and non-shaded organisms belong to the Domain Archaea. Samples CV – Caldeira Velha, NP – Nascente Poente, ESG – Esguicho de Maio, FCRG – Fumarola Caldeira da Ribeira Grande, PCRG – Piscina Caldeira da Ribeira Grande. Organisms: Acidimicrobium ferrooxidans DSM 10331 ( A. ferrooxidans DSM 10331), Acidithiobacillus caldus SM-1 ( A. caldus SM-1), Desulfurococcus amylolyticus DSM 16532 ( D. amylolyticus DSM 16532), Pyrobaculum aerophilum str. IM2 ( P. aerophilum str. IM2), Sulfurihydrogenibium azorense Az-Fu1 ( S. azorense Az-Fu1), Thermodesulfovibrio yellowstonii DSM 11347 ( T. yellowstonii DSM 11347), Thermofilum adornatus 1505 ( T. adornatus 1505), Thermoplasma acidophilum DSM 1728 ( T. acidophilum DSM 1728), Thermorudis peleae ( T. peleae ), Thermus antranikianii DSM 12462 ( T. antranikianii DSM 12462), Thermus scotoductus SA-01 ( T. scotoductus SA-01), Thiomonas intermedia K12 ( T. intermedia K12). CV NP ESG FCRG PCRG Organism % Organism % Organism % Organism % Organism % S. azorense Az-Fu1 58.6 P. aerophilum str. IM2 79.9 P. aerophilum str. IM2 96.0 T. acidophilum DSM 1728 98.7 A. caldus SM-1 64.9 T. adornatus 1505 9.9 T. antranikianii DSM 12462 12.4 S. azorense Az-Fu1 1.8 A. caldus SM-1 1.1 A. ferrooxidans DSM 10331 35.1 T. yellowstonii DSM 11347 7.2 T. scotoductus SA-01 5.3 D. amylolyticus DSM 16532 1.3 A. caldus SM-1 7.0 D. amylolyticus DSM 16532 1.3 A. ferrooxidans DSM 10331 4.9 T. adornatus 1505 1.1 T. acidophilum DSM 1728 4.8 P. aerophilum str. IM2 2.9 T. peleae 2.5 T. intermedia K12 2.1 Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |49 Once more, mTOUs algorithm, in accordance with the binning algorithm, revealed a high abundance of the hydrothermal characteristic taxonomic groups Aquificales ( Sulfurihydrogenibium azorense Az-Fu1) and Crenarchaeota ( Desulfurococcus amylolyticus DSM 16532, Pyrobaculum aerophilum str. IM2 and Thermofilum adornatus 1505) members (Strazzulli et al., 2017). In particular, Sulfurihydrogenibium azorense Az-Fu1 was identified in the samples CV and ESG, being the most abundant organism in the CV sample. Curiously, S. azorense Az-Fu1 was first isolated in Furnas, São Miguel, Azores (Aguiar et al., 2004) and observed in other metagenomics studies of the same site (Sahm et al., 2013). Therefore, the prediction of the presence of this bacterium on other hydrothermal sites suggests that this bacterium is well adapted to the physical and chemical environmental conditions of São Miguel Island. Being a representative of the Aquificales order, the metabolic capabilities of S. azorense Az-Fu1 in isolation, as well as within its microbial communities is therefore of huge interest. Chapter 2 50| PhD Thesis | Sophia Santos 2.4 Conclusions Since the mid-2000s, with the genesis of next-generation sequencing, metagenomic studies have transformed the understanding of microbial communities' function and dynamics through the global effort to catalog the uncultured microbial diversity throughout the most diverse environments (Gilbert et al., 2014; Rusch et al., 2007; Turnbaugh et al., 2007). More recently, the scientific community has focused on the study of environments that once were described as unhabitable (Rothschild et al., 2001) where a large number of organisms, known as extremophiles (Nicolaus et al., 2010) were revealed, presenting remarkable genomic and metabolic attributes (Durvasula et al., 2018) with industrial applications. In this work, the prokaryotic diversity of five samples from hydrothermal vents at São Miguel, Azores, was characterized, using metagenomic approaches. The taxonomic profiling of all samples was performed using assembly-based and read-based analysis algorithms. Although both approaches are usually used separately and selected depending on the origin of the sample in the study, the use of both is usually recommended (Quince et al., 2017) whenever possible. Indeed, both taxonomic profiling approaches presented very similar results demonstrating that the two approaches validate and complement each other. Differences were only spotted on the assignment of different members of the Pyrobaculum genus depending on the profiling algorithm, which can be justified by the use of different reference genome databases. High abundant organisms were identified in all samples and domain-level composition was also predicted. The results of both taxonomic profiling approaches showed a similar profile to the experimental 16S rRNA gene amplicon sequencing profile obtained. Differences were only encountered in sample PCRG, which showed a considerable detection of archaeal 16s RNA, however, no archaeal organisms were predicted to be present when applying the taxonomic profiling algorithms. No conclusion can be made on the presence or absence of archaeal diversity in the PCRG sample, since the taxonomic profiling algorithms available have indeed some limitations, due to the still big amount of unculturable organisms without a reference genome available in databases (Zorrilla, Patil, et al., 2021). Aquificales and Crenarchaeota members, regularly present in hydrothermal sites (Strazzulli et al., 2017), were predicted to be present in high abundance in samples CV, NP, and ESG, by both taxonomic profiling approaches. Specifically, Sulfurihydrogenibium azorense Az-Fu1 and Pyrobaculum islandicum Metagenomic study of thermophilic and hyperthermophilic environments from Azores PhD Thesis | Sophia Santos |51 DSM 4184/ Pyrobaculum aerophlium str. IM2, Aquificales, and Crenarchaeota members, respectively, appear to be abundantly present in various samples. Environments rich in Aquificales members have recently received scientific interest for being believed to be the earliest lineage within the domain Bacteria (Takacs-Vesbach et al., 2013) and for being a common presence in thermophilic and hyperthermophilic environments. Curiously, S. azorense Az-Fu1 was first isolated in January 2001 from terrestrial hot springs at Furnas, São Miguel Island, Azores, Portugal (Aguiar et al., 2004) and its re-occurrence at different sites with similar environmental conditions is worthy to be investigated. A systems biology approach to study S. azorense Az-Fu1 metabolic capabilities in isolation, as well as part of these communities, is suggested. Chapter 2 52| PhD Thesis | Sophia Santos 2.5 Supplemental Material Additional file in Excel format: Chapter2_Supplementary_Material.xlsx Link: DesignOptimizationMicrobialCommunities/Data Table S1 FastQC Report Summary Table S2 FastP Report Summary Table S3 MetaQuast Quality Results Table S4 BinSanity Complete Results Table S5 Kraken Results for the bins with high completeness (>97%) and low contamination (<10%) values and respective species identification Table S6 mOTUs Complete Results PhD Thesis | Sophia Santos |53 Chapter 3 Chapter 3 Genome-scale metabolic model of thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 “All for one, one for all” Alexandre Dumas i SS352 is the first manually curated genome-scale metabolic model for the thermophilic bacterium S. azorense Az-Fu1, comprising 352 genes and 772 reactions. The construction of this model involved performing a reannotation, which revealed the presence of the main components of the bacterial cellulose operon and its regulators, suggesting that S. azorense Az-Fu1 may have the metabolic potential for cellulose production. The model also clarifies this chemolithoautotrophic organism's carbon fixation route, central carbon, and sulfur metabolisms. The i SS352 metabolic model will serve the ongoing fundamental research of chemolithoautotrophic metabolism in extreme environments, provide clues regarding new extremophilic enzymes, and studies of interactions with the identified environmental microbial community in which it was discovered. Chapter 3 54| PhD Thesis | Sophia Santos 3.1 Introduction The Aquilicales order is of significant interest as it is believed to be the earliest-lineage within the domain Bacteria (Takacs-Vesbach et al., 2013). Being widespread in marine and terrestrial geothermal systems, members of Aquificales are thermophilic and hyperthermophilic bacteria almost exclusively obligate or facultative autotrophs (Hedlund et al., 2015). Considered to be primary producers, members of this order can also use a wide range of energy sources, such as H2, elemental sulfur, and thiosulfate (Takacs-Vesbach et al ., 2013, Yamanaka, 2008). Moreover, the study of extremophiles and their enzymes has increased in the last years (Counts et al., 2017; Han et al., 2019; Kumar et al., 2019) due to the need for the industry to have more stable biocatalysts at higher temperatures and pressures or extreme pHs to resist to harsh industrial settings (Atalah et al., 2019) and also for using less costly purification technologies, such as heat treatment that does not involve chemical reagents, protein tags, nor adsorbents (Sun et al., 2012). Sulfurihydrogenibium azorense Az-Fu1, a representative of the Aquificales order, commonly found as dominant taxa in hot spring vent communities, is a gram-negative, thermophilic, chemolithoautotrophic, and microaerophilic bacterium (Lalonde et al., 2005). The bacterium was isolated in January 2001 from terrestrial hot springs at Furnas, São Miguel Island, Azores, Portugal (Aguiar et al., 2004). S. azorense grows optimally at 68ºC, pH 6, and at low concentrations of NaCl and can also grow heterotrophically (Nakagawa et al., 2005) and use elemental sulfur, thiosulfate, hydrogen, and ferrous iron as energy sources (Aguiar et al., 2004). As an Aquificales member, S. azorense is believed to fix CO2 via reductive Tricarboxylic Acid Cycle (rTCA) (Hügler et al., 2007) to generate acetyl-CoA as an end product (Gupta et al., 2013), like its closely related relative Sulfurihydrogenibium subterraneum . In fact, the rTCA is reported to be the most plausible candidate for the first autotrophic metabolism in the earliest life (Kitadai et al., 2017). Also, (Lalonde et al., 2005) reported that, under stress conditions, this bacterium produces chemolithoautotrophically sufficient amounts of exopolysaccharides (EPS). EPS are polymeric structures of repeated sugar units of the same or different types that microorganisms usually produce under unfavorable growth conditions. These structures have a high range of physiological roles, such as preventing the entering of harmful substances within the cell, contributing to providing nutrients and/or water or promoting adherence to surfaces (biofilm formation) (Rana & Upadhyay, 2020). In extremophiles, EPS accumulation functions as a relevant adaptation strategy of cell protection and membrane stabilization (Kambourova et al., 2016). Due to their adaptation Genome-scale metabolic model of thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 PhD Thesis | Sophia Santos |61 stoichiometric balance and directionality, include or exclude reactions from the model, and finally, export models in Systems Biology Markup Language (SMBL) format to be used in simulation using platforms such as Optflux (I. Rocha et al., 2010), Matlab® or python (COBRApy (Ebrahim et al., 2013), REFRAMED - https://github.com/cdanielmachado/reframed) Gap-filling Before being ready for simulations, the metabolic network must be screened for possible gaps. The presence of gaps can deeply compromise the synthesis of biomass components and other relevant compounds. To assist in this process, merlin ’s “BioISO” (Capela et al., 2021) was used to trace back the network to identify gaps that can be originated from errors in genome annotation, absence of enzymatic, transport, or exchange reactions, or incorrect reaction irreversibility or direction. Other features included in merlin , such “Blocked reactions”, which identifies reactions that contain dead-end metabolites, and “Draw in Browser” which opens on a web browser a selected KEGG pathway, showing specifically highlighted enzymes and reactions present in the model, were used to facilitate the detection of gaps in the network (Capela et al., 2021). Literature and databases (KEGG, MetaCyc, BRENDA, for instance) were also used to assist the gap-filling process. This is an iterative process that is repeated until all biomass precursors, and other essential compounds, can be synthesized, and a feasible model is obtained. Microaerobic and anaerobic metabolism S. azorense is a facultative chemolithoautotrophic (Nakagawa et al., 2005), microaerophilic bacterium (Aguiar et al., 2004). Unlike plants, cyanobacteria, and other chemolithoautotrophic Proteobacteria that fix CO2 through the Calvin-Benson-Bassham cycle, S. azorense metabolism was assessed to determine if it fixates CO2 through the rTCA, such as its closely related organism S. subterraneum (Hügler et al., 2007) and most of Aquificales members. Anaerobic and microaerobic growth (specific growth rate of 0.28 h-1) was also tested as reported in the literature (Aguiar et al., 2004; Nakagawa et al., 2005). EPS production (Lalonde et al., 2005) was also screened. Chapter 3 62| PhD Thesis | Sophia Santos Carbon source utilization S. azorense can grow heterotrophically using yeast extract, bactopeptone, trypticase peptone, and casamino acids (Nakagawa et al., 2005). Except for casamino acids, all other medium components reported as carbon sources are not chemically defined; thus, only casamino acids were tested for heterotrophic growth. 3.2.2.5 Simulations Software All simulations were performed using the python package REFRAMED, version 1.1.0, using CPLEX 12.8.0 as a solver, through the PyCharm Integrated Development Environment (IDE). REFRAMED provides tools for phenotype simulation such as Flux Balance Analysis (FBA) (Varma & Palsson, 1994) Flux variability analysis Quantitative evaluation of the model was performed using FVA (Mahadevan et al., 2003) to understand S. azorense EPS production capabilities. The analysis included setting the specific growth rate to at least 10% of the specific growth rate obtained with pFBA (Lewis et al., 2010) in the respective reference flux distribution. 3.3 Results and Discussion 3.3.1 Genome annotation Assisted on the phylogenic tree developed in section 3.2.2.1 and merlin ’s feature automatic workflow, genome re-annotation identified 774 metabolic genes (Table 3.3) out of 1657 candidates (Table 3.5), representing 47% of the whole S. azorense Az-Fu1 genome. The complete list of genes reviewed is available in Supplementary Table S1. Genome-scale metabolic model of thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 PhD Thesis | Sophia Santos |63 Manual curation was performed and the number of genes in the model decreased to 352 (21%) due, for instance, to the removal of pseudo and truncated genes, incomplete EC numbers, and blocked reactions with encoded genes. Table 3.3. S. azorense genome annotation automatic workflow results taking into account phylogenetically related genus. Organism Gene count Sulfurihydrogenibium azorense Az-Fu1 6 genus Sulfurihydrogenibium 101 genus Persephonella 13 genus Aquifex 202 genus Hydrogenobaculum 25 genus Thermus 21 genus Thermotoga 5 genus Deinococcus 12 genus Desulfobacterium 3 Default annotation 385 The analysis of the annotation results shows the presence of bacterial cellulose synthase (Bcs) subunits A (BcsA – SULAZ_1378) and B (BcsB – SULAZ_1376). As described elsewhere (Lalonde et al., 2005), S. azorense Az-Fu1 can produce sufficient amounts of EPS, whose composition is still unknown, under stress conditions. Hence, the presence of annotated subunits of cellulose synthase suggests that cellulose can be a part of S. azorense EPS. However, cellulose production is dependent not only on the catalytic subunits BcsA and BcsB but also on several additional BCS components and gene regulators (Krasteva et al., 2017; Römling, 2002). Hence, S. azorense ‘s genome was screened to identify other cellulose operon subunits missing from the current annotation. Chapter 3 64| PhD Thesis | Sophia Santos The cellulose operon has three main types (Römling et al., 2015), as shown in Figure 3.1. All contain the catalytic subunits BcsA (cytosol) and BcsB (periplasm). These two subunits are essential for in vitro cellulose production. Subunits BcsC, present in type I and II operons, and BcsK, in type III are involved in cellulose transport, whereas subunit BcsD is present only in type I operon and is required for the Bcs complex arrangement. In type II and III operons, the BcsZ subunit regulates cellulose production. In type II operon, subunits BcsE and bcsG are required for optimal cellulose production, and subunit BcsQ seems to determine the cellular localization of the respective Bcs complexes (Römling et al., 2015). Other Bcs subunits have already been identified, but their function has not yet been unveiled. The regulation of the cellulose production process has also been described (Chang et al., 2001; Romling et al., 2013), showing the importance of bis-(3’,5’) cyclic diguanylic acid (c-di-GMP) molecule as an allosteric activator of subunit BcsA. In Escherichia coli , c-di-GMP is produced by GGDEF domain proteins such as AdrA or YedQ, and its degradation and shut-off of cellulose biosynthesis is achieved by specific phosphodiesterases, such as YhjH (Monteiro et al., 2009). Figure 3.1. Comparison of cellulose operon types I, II, and III, respectively, (A), (B), and (C), with (D) - predicted S. azorense Az-Fu1 cellulose operon. Colors indicate that a match was found as a result of NCBI tblastn searches. Red – BcsA unit, Dark gray – BcsB unit, Dark blue – BcsZ, Green – BcsC, Yellow – BcsQ. (C) – coverage; (I) – identity. As no information was found reporting the presence of any type of cellulose operon in S. azorense , gene accession numbers of all known units of the different BCS types and their regulators (Supplementary Genome-scale metabolic model of thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 PhD Thesis | Sophia Santos |65 Table 3.1) were analyzed. This analysis involved performing BLAST searches (NCBI tblastn) against S. azorense Az-Fu1’s whole genome to identify the presence of BCS units and the respective BCS type, if possible. As a result, besides the already described presence of BcsA and BcsB, the BcsC and BcsZ units were found in the same order type II BCS operon, suggesting that if present, S. azorense ’s cellulose operon will be a type II. Although all units were found to have adequate identity percentages (Figure 3.1), BcsB and BcsC show low coverage values. A gene sequence similar to the BcsQ unit was also found but not at the same position as in type II operons, which renders its identification inconclusive. Nevertheless, the re-annotation of genes SULAZ_1375 and SULAZ_1374 should be updated as bacterial cellulose synthase Z unit and bacterial cellulose synthase C unit, respectively. Regulators AdrA, YedQ, and YhjH were also identified with acceptable coverage and identity percentages (AdrA - C.50.3%, I. 35.2%; YedQ - C.42%, I.35.9%; YhjH – C.64.0%, I.28:7%). Although the main Bcs units and regulator factors were identified in the S. azorense Az-Fu1 genome, the bacterium's actual capability of producing cellulose must be in vivo experimentally assessed. 3.3.2 Biomass Composition The biomass macromolecular composition, adapted from E. coli (Feist et al., 2007) is presented in Table 3.4. The detailed biomass composition is available in Supplementary Table S2. As mentioned before, the amino acid, deoxynucleotide, and nucleotide composition were calculated based on S. azorense Az-Fu1 genome’s information using the e-Biomass feature in merlin . This feature also automatically includes cofactor composition based on the study of universal essential cofactors in prokaryotes (Xavier et al., 2017). Usually, ubiquinone is included in the cofactors pool. However, most enzymes required for ubiquinone biosynthesis are not available in the genome annotation, and the defined medium (Aguiar et al., 2004) does not include ubiquinone; thus, this cofactor was omitted in the biomass formulation. The lipids composition was adapted from the phylogenetically closely related organism H. thermophilus (Yoshino et al., 2001). However, S. azorense Az-Fu1 does not exhibit the enzymes responsible for producing phosphatidylserine, phosphatidylethanolamine, and phosphatidylcholine in the model. Hence, these compounds were excluded from biomass composition. The coefficients of the remaining compounds were recalculated, maintaining the relative abundances of the original data. The Chapter 3 66| PhD Thesis | Sophia Santos average fatty acid composition was adopted from S. subterraneum (Takai et al., 2003) (Supplementary Table S3). Table 3.4. Biomass macromolecular composition of the S. azorense Az-Fu1 model. Biomass Composition g gDW-1 (%) Protein 53.3 DNA 2.7 RNA 13.6 Lipids 2.9 Carbohydrates 10.7 Cell Wall 6.8 Cofactors 10 Total 100 Cell wall and carbohydrate components were adapted from E. coli (Feist et al., 2007). Specifically for cell wall components, E. coli ’s composition was reconciled between the KEGG reactions assigned through merlin ’s annotation to its biosynthesis pathway (Supplementary Table S4). The growth-associated energy (GAM) and non-growth-associated energy (NGAM) requirements for S. azorense are not been experimentally determined yet. Therefore the GAM requirements of 56.64 mmol ATP gDW−1 h−1 were estimated according to (Thiele et al., 2010) and based on data for E.coli (Neidhardt et al., 1990). The NGAM requirements of 8.39 mmol ATP gDW−1 h−1 were adopted from E. coli (Feist et al., 2007). 3.3.3 Metabolic Model The genome-scale metabolic reconstruction of S. azorense Az-Fu1, i SS352, was generated through a bottom-up approach. The final reconstruction contains 352 genes, 772 reactions (62 exchange, 141 Genome-scale metabolic model of thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 PhD Thesis | Sophia Santos |67 transport, and 25 gap-fill reactions), and 642 metabolites (Table 3.6), distributed over 70 subsystems and three compartments: extracellular, periplasm, and cytoplasm (Supplementary Tables S6, S7, and S8). The gene ratio is about 21%, a value in accordance with other published GSMMs for chemolithoautotrophic organisms (Supplementary Table S5), and 70% of all reactions have a GPR associated (Table 3.6). Reactions without GPR associations include spontaneous, exchange, outer membrane transport, and diffusion transport reactions of metabolites, such as CO2, water, or O2. The final GSMM, identified as i SS352, is included in Supplementary Material ( i SS352.xml) in the SBML level 3 version 2 format. The model was able to score over 97% on all consistency tests on the MEMOTE (Lieven et al., 2020) test suit. Table 3.5. S. azorense Az-Fu1 genome information Genome Information Genome length (bp) 1640877 G+C content (%) 32.8 No. of ORFs 1657 Table 3.6. i SS352 metabolic model information Metabolic Model Reactions 772 Metabolites 642 Genes 352 Gene Rules 546 3.3.4 Metabolism of S. azorense as represented in the i SS352 model 3.3.4.1 Carbon metabolism Autotrophic CO2 fixation is one of the most critical biosynthetic processes in nature. The CalvinBenson-Bassham cycle is one of the most recognized pathways, mainly used by plants and cyanobacteria. Additional CO2 fixation mechanisms that motivate the scientific community to explore possible alternative CO2 fixation routes for reducing the human carbon footprint (Bar-Even et al., 2010; Liu et al., 2020) include five natural autotrophic pathways, namely the rTCA, dicarboxylate/4-hydroxybutyrate cycle, 3hydroxypropionate cycle, 3-hydroxypropionate/4-hydroxybutyrate, and Wood-Ljungdahl cycles. Chapter 3 68| PhD Thesis | Sophia Santos Figure 3.2. Sulfurihydrogenibium azorense Az-Fu1 proposed central carbon metabolism under chemolithoautotrophic growth and a potential route for cellulose production. rTCA – reverse Tricarboxylic Citrate Acid cycle; PPP – Pentose Phosphate Pathway; EMP - Embden-Meyerhof-Parnas (EMP) glycolytic pathway. Genome-scale metabolic model of thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 PhD Thesis | Sophia Santos |69 As reported for the closely related organism S. subterraneum (Hügler et al., 2007) and as most members of the Aquificales order (Gupta et al., 2013), S. azorense Az-Fu1 seems to fix CO2 through the rTCA and produce all biomass components (Figure 3.2). The presence of three specific enzymes characterizes this fixation pathway: fumarate reductase (SULAZ_0630/R00408), ferredoxin-dependent 2-oxoglutarate synthase (SULAZ_0641/R01197), and ATP-citrate lyase (SULAZ_0527/R00352), all present in the S. azorense Az-Fu1 genome. Other CO2 fixation pathways have been screened. However, several vital enzymes were missing from the organism’s genome. Moreover, given the fact that S. azorense Az-Fu1 is a thermophilic organism that grows under microaerobic conditions, which is an energylimiting environment, it makes sense that the rTCA cycle is the preferred route of CO2 fixation for this organism. Compared to the Calvin-Benson-Bassham cycle, the rTCA cycle requires significantly less energy as ATP and fewer reducing equivalents to synthesize a three-carbon unit (Dahl et al., 2008). Surprisingly, S. azorense has a reduced central carbon metabolism, lacking all the oxidative part of the Pentose Phosphate Pathway (PPP) and the enzyme fructose-bisphosphate aldolase, originating an incomplete Embden-Meyerhof-Parnas (EMP) glycolytic pathway (Figure 3.2). Initially described as an obligate chemolithoautotroph (Aguiar et al., 2004), able to use a variety of electron donors and acceptors, S. azorense Az-Fu1 was later considered a bacterium with the ability to grow heterotrophically in yeast extract, bactopeptone, trypticase peptone, and casamino acids (Nakagawa et al., 2005). In fact, amino acid transporters were predicted to be present in S. azorense Az-Fu1’s genome, suggesting these as the only carbon source besides CO2. Moreover, no sugar transporters other than cellulose were predicted in this work. As mentioned in section 3.2.2.1, S. azorense Az-Fu1 has the main cellulose biosynthesis operon subunits and respective regulators in its genome, suggesting that cellulose production seems metabolically plausible. 3.3.4.2 Sulfur metabolism As a sulfur-oxidizing bacterium, the study of S. azorense Az-Fu1 sulfur metabolism is of significant interest. Sulfur chemolithotrophy is regarded as the earliest self-sustaining metabolism, as it is found in several extremophiles of the deepest phylogenetic branches of Archaea and Bacteria (Ghosh et al., 2009). In fact, sulfur occurs in nature in oxidation states ranging from +2 to -6, and is thus used by prokaryotes to build cell constituents and as an energy source (Dahl et al., 2008). Over the last years, several studies have attempted to identify the sulfur-oxidizing pathways. Neutrophilic sulfur-oxidizing bacteria, such as S. Chapter 3 70| PhD Thesis | Sophia Santos azorense , seem to use two pathways for sulfur oxidation: the multienzyme complex catalyzing the complete oxidation of reduced sulfur compound to sulfate (Sox pathway) and a pathway having sulfite and sulfur as important intermediates (via sulfide:quinone reductase) (Friedrich et al., 2001). According to the annotation, S. azorense Az-Fu1 seems to have the Sox enzyme system. This system is present in the periplasm and is responsible for the oxidation of thiosulfate to sulfate. When other sulfur sources are used a first step of conversion into thiosulfate is needed. The Sox gene cluster soxXYZABCD encodes four periplasmic proteins, SoxXA, SoxYZ, SoxB, and Sox(CD)2 (sulfur dehydrogenase). In S. azorense Az-Fu1 the protein Sox(CD)2 is absent, indicating that this bacterium is likely to have a truncated Sox system, similar to the closely related organism Aquifex aeolicus (Friedrich et al., 2005). Figure 3.3. Sulfurihydrogenibium azorense Az-Fu1 proposed sulfur metabolism through a truncated sulfur-oxidizing (Sox) system. Here elemental sulfur is being used as the main sulfur source. The proposed process is similar when hydrogen sulfide, sulfite, or thiosulfate are used as a sulfur source. Studies in the purple sulfur bacterium Allochromatium vinosum , which also has a truncated Sox system lacking Sox(CD)2, indicate that elemental sulfur is produced from thiosulfate throughout the sulfuroxidizing process and transferred to the cytoplasm (Figure 3.3) or growing sulfur globules. The mechanism for this process is still unresolved, although it may involve a rhodanese-like enzyme (SoxL) (Dahl & Friedrich, 2008). The truncated Sox enzyme system is less energy efficient once it yields only two Genome-scale metabolic model of thermophilic bacterium Sulfurihydrogenibium azorense Az-Fu1 PhD Thesis | Sophia Santos |77 Supplementary Figure 3.1. Phylogenetic tree of Sulfurihydrogenibium azorense related genus. This tree was built using the EMBL-EBI Clustal OMEGA multiple sequence alignment tool. Numbers in front of each genus represent the branch lengths to each node generated automatically by the tool using the Neighbour-joining method. Supplementary Table 3.1. Cellulose operon units for each main type and respective gene accession numbers. Gene Type I Komagataeibacter xylinus Type II Escherichia coli Type III Methylobacterium extorquens bcsA / YhjO AAA21884.1 NP_417990.4 WP_080518407.1 bscB / YhjN AAA21885.1 NP_417989.1 WP_012252997.1 bcsC/ YhjL AAA21886.1 YP_026226.4 --------- bcsD AAA21887.1 --------- --------- bcsQ/ YhjQ --------- WP_011310329.1 --------- bcsZ / YhjM --------- NP_417988.1 WP_012252998.1 bcsE/ YhjS --------- NP_417993.1 --------- bcsF/ YhjT --------- NP_417994.2 --------- bcsG/ YhjU --------- NP_417995.1 --------- bcsR/ YhjR --------- NP_417992.1 --------- bcsH/ CcpA AAA16970.1 --------- --------- bcsK --------- --------- WP_012252999.1 bcsN --------- --------- WP_012253000.1 bcsS --------- --------- WP_003606498.1 Chapter 3 78| PhD Thesis | Sophia Santos Supplementary Table 3.2. Simulation results under chemolithoautotrophic and microaerophilic conditions. Simulation Results Specific growth rate (µ): 0.2636 h-1 Consumption Production Metabolite Uptake Value (mmol gDW-1 h-1) Metabolite Production Value (mmol gDW-1 h-1) CO2 12.00 Sulfate 12.23 H2O 4.286 Sulfur 12.27 Fe2+ 0.004437 Ammonia 2.674 O2 9.518 H3PO418.90 Supplementary Table 3.3. Simulation results under mixotrophic and anaerobic conditions. Simulation Results Specific growth rate (µ): 0.1069 h-1 Consumption Production Metabolite Uptake Value (mmol gDW-1 h-1) Metabolite Production Value (mmol gDW-1 h-1) Amino Acids < 0.1 H+ 4.190 Sulfur 0.001790 H2O 1.660 Fe2+ 0.001790 Ammonia 0.4677 H3PO40.1106 PhD Thesis | Sophia Santos |79 Chapter 4 Chapter 4 Microbial Community Simulation Methods “The lone wolf dies but the pack survives.” George R. R. Martin Microbial communities directly affect surrounding environments and are an important part of all biological processes. Thus, the study of microbial communities’ behavior and composition has been proven to be useful in areas such as biotechnology, environmental, and human health. However, the overall understanding of microbial communities’ interactions and dynamics remains a challenge. Synergies between computational methods and genome-scale metabolic models have been explored in the last years, as a way to unravel community interactions and behavior, as demonstrated by the numerous simulation methods developed for application in the context of microbial communities. Here, different steady-state simulation methods applied to microbial communities have been used to model the well-established nitrification bioprocess catalyzed by Nitrosomonas euroapaea and Nitrobacter vulgaris . The different methods' performances were compared to assess which method(s) should be used in a specific community-level context. The available simulation methods, with application to microbial communities, revealed good phenotypic behavior predictions. Each of the simulation methods exhibited strengths and weaknesses. Hence, to better predict the communities’ behavior, it is recommended the use of various simulation methods, whenever possible. Chapter 4 80| PhD Thesis | Sophia Santos 4.1 Introduction During the last decades, genome-scale metabolic models (GSMMs) have proven to be a valuable tool in Systems Biology, with outstanding biotechnological applications (Gu et al., 2019), predicting cellular behavior under different genetic and environmental conditions (Woolston et al., 2013), and reducing time and costs implicated in experimental tasks. In the last decade, the development of high-throughput techniques and the scientific desire to understand microbial communities’ behavior and capabilities, through computational and mathematical modeling approaches, turned communities’ GSMMs into the new frontier in the field. The emergence of numerous simulation methods developed to be applied in a microbial community context (Colarusso et al., 2021), with proven applicability in industrial, environmental biotechnology, and human health environments (García-Jiménez et al., 2021) support this claim. The currently available simulation methods, developed to perform simulations with microbial communities GSMMs, can be classified as steady-state, time-course, and spatial-temporal (Figure 4.1). Steady-state methods allow for predicting individual growth rates and cross-feeding interactions in a situation of community equilibrium. These methods are based on the well-established Flux Balance Analysis (FBA) (Varma et al., 1994) method for individual species. This constraint-based approach is usually applied to biochemical networks, allowing the prediction of biological phenotypes. The steadystate reaction flux distribution in a network is determined by Linear optimization of biologically and ecologically meaningful objective(s) function(s), subject to a set of underlying constraints (Raman et al., 2009). However, these methods are characterized by a large uncertainty in the space of optimal solutions. Time-course methods, based on dynamic flux balance analysis (dFBA) (Mahadevan et al., 2002), account for the temporal variation in metabolite concentrations and species abundance, allowing for simulation of the response to changes over initial conditions and other environmental perturbations. However, this approach requires the characterization of the substrate uptake kinetics for consumed metabolites (either from the growth medium or through cross-feeding). The impractical nature of such a massive in vitro characterization results in the adoption of default parameter values for all species and compounds, which limits the predictive ability of these methods. Microbial Community Simulation Methods PhD Thesis | Sophia Santos |81 Figure 4.1 Summary of the main features of the current simulation methods with application to microbial communities. Spatiotemporal methods add a spatial component to the temporal one. These methods are especially relevant in situations of non-heterogeneous or spatially segregated environments such as biofilm colonies or the tract along the human gut, where nutrient diffusion and access to resources play a role in community assembly and potential for cross-feeding interactions. In addition to the parameters on uptake kinetics mentioned before, these methods also require the specification of diffusion rates for cells and metabolites. Besides having common grounds, as all are based on the FBA method for single organisms, these methods differ in the number and kind of inputs needed and have different levels of simulation detail and Chapter 4 82| PhD Thesis | Sophia Santos objective functions. Hence, each method can produce a variety of output information, enabling the handling of microbial communities’ GSMMs under a wide range of different applications. However, the increasing level of simulation detail simultaneously requires more experimental data for model setup and, naturally, more computational power, which is only possible for small communities composed of wellcharacterized microorganisms. Steady-state simulation methods benefit from not requiring parameterization, thus being scalable for large microbial communities. Also, they are compatible with flux variability analysis (FVA) (Gudmundsson et al., 2010), flux coupling analysis (David et al., 2011), and random flux sampling (Schellenberger et al., 2009). Hence, developing steady-state methods for the analysis of microbial communities (Figure 4.2) (Colarusso et al., 2021) has been the main focus of research. Having a great variety of methods, the most difficult task is understanding which method(s) should be applied in a specific context, considering available experimental data, microbial community complexity, and particular final purpose(s). Understanding the differentiating features allows selecting the best strategy for a specific microbial community. Figure 4.2. Lineage of steady-state simulation methods with application to microbial communities. Each branch of the tree represents different simulation method assumptions. Microbial Community Simulation Methods PhD Thesis | Sophia Santos |83 Here, a systematic analysis of the performance of the most used steady-state simulation tools for predicting microbial communities’ behavior is presented. This analysis will cover key features that allow properly selecting the best method, that fits not only the purpose of the analysis but also available experimental data. For that, the most promising simulation methods have been installed and applied to the well-known case of the biological nitrification process catalyzed by the symbiotic interaction between Nitrosomonas europaea and Nitrobacter vulgaris (Ilgrande et al., 2018). The GSMMs for these bacteria have been reconstructed and manually curated with experimental data within our research group (Cruz, 2018; Raposo, 2017). 4.1.1 Current state of steady-state simulation methods with application to microbial communities A brief description of the current steady-state simulation methods used in a microbial community context is presented hereafter, highlighting the main features of each one. Flux Balance Analysis Although its formulation was developed for predicting the phenotypic behavior of single organisms, Stolyar et al pioneered the use of FBA to simulate the behavior of microbial communities, namely in studying a community of two methanogenic species (Stolyar et al., 2007). The steady-state model must be created before the FBA simulation by merging the stoichiometric matrices of any number of individual organisms. FBA can determine steady-state flux distribution throughout the community network, restricted to a defined medium and corresponding uptake rates of the exchanged metabolites. Hence, by analyzing flux distribution across the network, community, and individual growth rates, as well as cross-feeding interactions in the community can be easily established. Although the definition of an objective function for microbial communities is still ambiguous, FBA and the variant parsimonious enzyme-usage flux balance analysis (pFBA) (Lewis et al., 2010) have been successfully used using growth maximization of the entire microbial community as the objective function. Examples of applications include understanding microbe-microbe and host-microbe interactions (Bordbar et al., 2010; Heinken et al., 2013) and the production of biocompounds in diverse settings, such as the production of vitamin C (Ye et al., 2014) or 1,3-propanediol (Bizukojc et al., 2010). Chapter 4 84| PhD Thesis | Sophia Santos FBA can be used under diverse platforms, such as MATLAB® through the COBRA toolbox (Heirendt et al., 2019), the python packages COBRApy (Ebrahim et al., 2013) and REFRAMED, or user-friendly frameworks such as Optflux (I. Rocha et al., 2010) or KBase (Arkin et al., 2018). OptCom OptCom (Zomorrodi et al., 2012) is an FBA framework that relies on a bi-level optimization structure. The inner level maximizes the growth of individual organisms, while the outer level maximizes total community fitness. The underlying assumption is that the community will optimize the usage of available resources while not compromising the individualistic objectives of its members. If more constraints are provided, OptCom is also able to compute cases where organisms do not present their maximum individual growth, exhibiting cooperative behavior instead. This approach can use any number of GSMMs and returns strain ratio, substrate uptakes, and secretion rates showing applicability in inferring the syntrophic association through hydrogen between Desulfovibrio vulgaris and Methanococcus maripaludis (Zomorrodi et al., 2012) and the production of butyrate for cancer prevention (El-Semman et al., 2014). OptCom has been implemented in the COBRA toolbox (Heirendt et al., 2019) and integrated into the MICOM Python framework (Diener et al., 2020). cFBA Community Flux Balance Analysis (cFBA) (Khandelwal et al., 2013) is a computational method whose primary focus is elucidating the metabolic capabilities of a community and understanding metabolic interactions. This approach is a direct translation of FBA to microbial communities requiring balanced growth and postulating a single objective. Balanced growth implies that the entire community is in steady-state and thereby all metabolites (intraand extracellular) are at a steady-state level as well. Although the problem definition of cFBA is non-linear, when biomass fractions are considered a variable, by fixing the individual biomass values, a linear programming problem arises, thus providing a result that identifies the optimal specific flux values. Notice that the community growth rate may not match the maximum growth rate of any individual organism. However, when considering a mutualistic interaction between growing organisms, to preserve the cross-feeding metabolites production and consumption rates equilibrium, organisms must grow at an equal pace. Along with the abundance of all Microbial Community Simulation Methods PhD Thesis | Sophia Santos |85 species present in the community, cFBA predicts fluxes distribution, growth rates, and exchange fluxes between the microbes and the environment. Hence, the GSMMs of the organisms present in the community must be merged before simulation. The cFBA method is applied essentially to microbial communities exposed to constant environmental conditions, such as the ones involved in bioremediation, wastewater treatment, or in laboratory settings. SMETANA Species Metabolic Interaction Analysis (SMETANA) (Zelezniak et al., 2015), a Python-based command line tool, uses mixed-integer linear programming to estimate the interaction potential of the species in a microbial community and returns the probability of inter-species metabolite exchange capacity. SMETANA can be applied with few microbial community information, such as species identity and genome sequence. GSMMs for each species are automatically reconstructed and assembled into a community model. This step can be ignored if GSMMs of all species in the community are available. SMETANA calculates two scores to predict the species’ interaction potential: the Metabolic Resource Overlap (MRO) and the Metabolic Interaction Potential (MIP). The MRO is an intrinsic property of any community and is defined by the maximum possible overlap between the minimal nutritional requirements of all member species. MIP represents the tendency of community members to exchange metabolites and is given by the maximum number of essential nutritional components that a community can provide for itself, through the interspecies metabolic exchange. Thus, the higher its value, the higher the probability of community members benefiting from metabolite production from other members. MRO and MIP represent opposite circumstances. While the first indicates the predisposition to competition between the organisms, as both require the same metabolite(s) from the environment, MIP represents the tendency of the community’s organisms to rely on each other, not being able to grow on their own. SMETANA calculates three other scores to predict potential inter-species interactions, Species Coupling Score (SCS), Metabolic Uptake Score (MUS), and Metabolite Production Score (MPS). The SCS score measures the growing dependence of species A in the presence of species B, whereas the MUS score measures the growing dependence of species A in the presence of metabolite m , provided by the other members of the community. MPS is a binary score that reflects the ability of a given species to produce metabolite m . Finally, the SMETANA score that evaluates the growth dependency of species A on metabolite m produced by species B is calculated as a product of the scores SCS, MUS, and MPS. Chapter 4 86| PhD Thesis | Sophia Santos SMETANA has been successfully applied to communities composed of a high number of species, spanning habitats as diverse as soil, water, and human gut (Machado et al., 2021; Zelezniak et al., 2015; Zorrilla et al., 2021) as well as in a three-species stable synthetic community, with nitrogen overflow capabilities (Ponomarova et al., 2017). OptDeg OptDeg (Koch et al., 2016), uses as input a community model for predicting growth and flux rates. This approach applies the previously described concept of balanced growth, in which all organisms in the community grow at the same growth rate. For predicting the community composition, a hierarchical optimization approach that consists of two objective functions is employed. Firstly, the maximization of the community’s growth rate is optimized, followed by the maximization of each organism's growth rate under optimal utilization of substrates, which is defined by the authors, as biomass yield. For the quantification of the overall biomass yield, an optimality degree constant (OptDeg) was introduced. This constant is calculated as the quotient of the maximal community growth rate and the individual minimal expected growth rate if the organism uses its substrates optimally. Thus, for OptDeg = 1, all the community species grow at the maximum specific growth rate and maximum biomass yields. OptDeg has been implemented in CellNetAnalyzer , a MATLAB package for structural and functional analysis of metabolic and signaling networks (Klamt et al., 2007). SteadyCom SteadyCom (Chan et al., 2017) is an FBA-based method that predicts the metabolic flux distribution and relative abundance of each species in a community, assuring that a steady state is imposed. The novelty of this implementation relies on the imposition of a time-average constant growth rate to ensure the co-existence and stability of all members of the community. For instance, when a single organism is considered, the biomass flux is normalized by the organism’s rates of consumption or production. However, when multiple organisms are growing, there is not a constant growth rate for all microbes. Therefore, the fastest-growing organism can outgrow the rest of the community members. To avoid such situations, SteadyCom imposes a steady-state condition, that includes a restriction to force zero flux through an organism with zero abundance. As the sub-problems to be solved are independent of Microbial Community Simulation Methods PhD Thesis | Sophia Santos |93 with the experimental data reported by Cruz (2018). Results of the simulation of N. europaea GSMMs (Table 4.4) were in agreement with the information available in the literature (Raposo, 2017). Table 4.4. N. europaea model validation using pFBA as simulation method. Specific growth rate (h-1), uptake (Consumption), and export (Production) rates (mmolgDW-1h-1) are presented. Nitrosomonas europaea Specific growth rate: 0.007528 h-1 Consumption Production CO2 0.3000 H+ 0.6654 Orthophosphate 0.002598 H2O 1.182 SO420.001556 NO20.6606 Ammonia 0.7265 Urea 0.001338 O2 0.6624 Fe2+ 0.00001359 Table 4.5. N. vulgaris model validation using pFBA as simulation method, after manual curation of the oxidative phosphorylation pathway. Specific growth rate (h-1), uptake (Consumption), and export (Production) (mmolgDW-1h-1) rates are presented. Nitrobacter vulgaris Specific growth rate: 0.03179 h-1 Consumption Production CO2 0.0100 H+ 0.8133 Orthophosphate 0.0005715 NO30.2612 SO420.0004815 NO20.3134 O2 0.06714 Fe2+ 0.0002978 Ethanol 0.250 N. europaea was able to use ammonia and almost totally convert it into NO2- (91% of nitrogen conversion) as reported in the literature, whereas N. vulgaris showed a slightly lower proportion of conversion of NO2into NO3- (83% of nitrogen conversion), yet consistent with experimental results (Cruz, 2018) when no ammonia is supplemented to the medium and NO2is also used as a nitrogen source for growth. Chapter 4 94| PhD Thesis | Sophia Santos 4.3.3 Steady-state simulations of community Each method's features were analyzed according to method-specific outputs. These outputs include 1) the organism’s individual growth rate within the community; 2) carbon source and ammonia consumption; and 3) cross-feeding metabolites, namely NH3+, NO2-, and NO3-. A summary of all simulations is provided in Table 4.6. 4.3.3.1 FBA Initially, a community model was created using the REFRAMED Python package. When applying REFRAMED’s Community function, providing the individual GSMMs for each organism as input, a compartmentalized model is returned. Given the fact that the biomass equation on the community model is defined as the sum of the individual biomass equations at the same proportion when applying the FBA method to microbial communities, all organisms contribute equally to the community biomass. Thus, to maintain the steady-state balance, all specific growth rates are equal. In this specific case (Table 4.6), the results show that the specific growth rate of each organism in the co-culture (0.0077 h-1) mimics the slower-growing organism, N. europaea (0.007528 h-1). Microbial Community Simulation Methods PhD Thesis | Sophia Santos |95 Table 4.6. Microbial community, composed of N. europaea and N. vulgaris , simulation using pFBA, SMETANA, SteadyCom, MICOM, and OptCom as simulation methods. When available, specific growth rate (h-1), uptake (Consumption), export (Production) (mmolgDW-1h-1), and interaction rates are presented. Consumed metabolite rates are represented as negative and produced metabolite rates are represented as positive. N.vu. – Nitrobacter vulgaris , N.eu. – Nitrosomonas europaea , Com – Community. n.a. – Data not available. *Abundancy value used as input for the simulation method. Experimental pFBA SMETANA SteadyCom MICOM OptCom N.vu. N.eu. Com N.vu. N.eu. Com N.vu. N.eu. Com N.vu. N.eu. Com N.vu. N.eu. Com N.vu. N.eu. Com Specific Growth Rate n.a. n.a. 0.0147 0.0077 0.0077 0.0154 n.a. n.a. n.a. 0.0243 0.000199 0.0244 0.0150 0.0075 0.0151 0.0138 0.0075 0.0139 Metabolite Consumption/Production CO2 n.a. n.a. n.a. -0.124 -0.297 -0.421 n.a. -0.403 -0.0177 n.a. -0.195 -0.0441 -0.245 n.a. n.a. -0.239 NH4 n.a. n.a. n.a. -0.0308 - 1.22 -1.25 Producing Consuming n.a. -1.701 -0.100 n.a. -0.824 -0.236 -1.06 n.a. n.a. -1.06 NO2Consuming Producing n.a. -1.15 1.15 0 Consuming Producing n.a. 1.700 -1.700 n.a. 0 0.999 0.999 n.a. n.a. 0.999 Urea n.a. n.a. n.a. 0 0.00128 0.00128 n.a. 2.05 x 10-5 0 n.a. 0.145 0 0.145 n.a. n.a. 0.143 NO3n.a. n.a. n.a. 1.15 0 1.15 n.a. 0 1.700 n.a. 0 0 0 n.a. n.a. 0 Abundancy n.a. n.a. 0.5* 0.5* n.a. n.a. n.a. 0.99 0.01 0.99* 0.01* 0.99* 0.01* Chapter 4 96| PhD Thesis | Sophia Santos The analysis of the simulation results for the exchange fluxes demonstrates that, as reported in the literature, the microbial community converts NH3+ into NO3-. Moreover, the proportion of conversion is similar to the individual organisms (92% of nitrogen conversion), using CO2 as a unique carbon source. Furthermore, the interacting reactions’ fluxes (Table 4.6), show that the uptake of CO2 is analogous to the uptake rate of carbon source for each organism, where N. europaea is consuming 0.31 mmolgDW1h-1 of CO2 growing in isolation (Table 4.4) and -0.297 mmolgDW-1h-1 of CO2 growing in co-culture with N. vulgaris (Table 4.6). Although the total CO2 uptake is similar, the growth rate of N. vulgaris within the community is 10-fold smaller, which was already reported in the literature (R. Cruz, 2018) which states that N. vulgaris has smaller growth rates when CO2 is the sole carbon source. As shown in Table 4.6 and the literature (Cruz et al., 2018; Grunditz et al., 2001; Mellbye et al., 2018), when NH3+ is used by N. vulgaris as nitrogen source, all NO2excreted by N. europaea is consumed by N. vulgaris , demonstrating the biological efficiency of the nitrification process catalyzed by these organisms. 4.3.3.2 SMETANA SMETANA estimates the interaction potential of the species in a microbial community and returns the probability of inter-species metabolite exchange. No growth and flux rates are directly obtained through SMETANA. For the case study of the community composed of N. europaea and N. vulgaris , a pre-simulation step was required. The individual GSMMs were reconstructed using merlin (Dias et al., 2018), which relies on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2016), for metabolites and reactions identifiers. However, SMETANA requires BiGG metabolite identifiers, thus those exchange metabolites identifiers were converted. The SMETANA potential inter-species interactions are available in Supplementary Tables S2 and S3. Only inter-species interactions with values greater than 0.4 were considered. While analyzing possible inter-species interactions, as expected, and reported by literature, N. europaea acts within the community as a donor of nitrogen sources. Surprisingly, NH3+ seems to have a high probability of being the nitrogen source shared by the organisms. Indeed, when simulating individual growth, when NH3+ is added to the medium, N. vulgaris rather uses it as a nitrogen source to produce Microbial Community Simulation Methods PhD Thesis | Sophia Santos |97 biomass components (R. Cruz, 2018). SMETANA also predicted a high dependence on glycine, with a SMETANA score of 0.68 (data not shown), which can act as a carbon source and has a relevant role in assimilating one-carbon compounds through the reductive glycine pathway of chemolithoautotrophic organisms, as discussed in previous studies (Sánchez-Andrea et al., 2020; Yishai et al., 2018). However, experimental validation is required to confirm this possibility. When analyzing the species interaction potential, SMETANA’s results show no growth dependence on each other, indicating that both organisms can grow independently (MIP score = N/A). However, the resulting MRO score of 0.6 indicates the tendency to competition between the organisms for the main elemental sources, namely the nitrogen sources. 4.3.3.3 SteadyCom SteadyCom predicts the metabolic flux distribution and relative abundance of each species in a community. Using the same environmental conditions as before (Table 4.3), SteadyCom predicted the growth of the two organisms, with a community growth rate of 0.02441 h-1. As shown in Table 4.6, despite small numerical differences, the metabolites’ interaction profile and respective ratios are quite similar to the experimental data. Importantly, this prediction shows, as predicted by the pFBA method and as reported in the literature (Cruz et al., 2018; Grunditz et al., 2001; Mellbye et al., 2018), that all NO2excreted by N. europaea is consumed by N. vulgaris . SteadyCom also calculates relative abundances for all organisms in the microbial community. Regarding this case study, SteadyCom's results are corroborated by individual growth experiments, which show that N. europaea grows at a much lower rate than N. vulgaris resulting in relative abundancies of 0.01 and 0.99, respectively. The predicted growth rates using this method are in the same range for N. vulgaris when compared with isolated growth, and for N. europaea 10-fold lower when compared to isolated growth under the same conditions. Within a community, microorganisms' growth rates’ perturbations are expected when organism abundances are significantly distinct, as seems to be the case for this community. However, these results cannot be assessed as additional experimental trials are required. Chapter 4 98| PhD Thesis | Sophia Santos 4.3.3.4 MICOM MICOM mathematical formulation also considers two different classifications for growth rates: individual growth rates and community growth rates. However, MICOM’s implementation does not assume a fixed community growth, thus requiring the community ‘s taxa abundance. The first analysis was performed using the same environmental conditions as before and setting identical abundancies for both organisms, such as defined to run the pFBA simulation. MICOM’s simulation results showed that in these conditions N. europaea is unable to grow (data not shown). MICOM predicts a given microbial community’s growth by fixing the minimum for all community members (equal or different for each community member). Therefore, when fixing the minimum growth rate to the slowest-growing organism in this community ( N. europaea - a growth rate of 0.0075 h-1), MICOM predicts the growth of both organisms (data not shown). However, the growth rates, as well as uptake and withdrawal rates, are 10-fold and 100-fold greater than the ones predicted with the previous methods. These results are probably associated with the fact that MICOM implementation scales the whole community’s internal exchange fluxes to the respective organism’s abundance. To prove this hypothesis, the organism's input abundances were adjusted to the abundance values predicted by SteadyCom (0.01 for N. europaea and 0.99 for N. vulgaris ). The results corroborate the suggested hypothesis, as the growth rates and uptake/production rates are in the same range of the pFBA and SteadyCom methods’ predictions (Table 4.6) However, in both predictions, the nitrogen flow through the networks (Supplementary Table S4) does not match the previous predictions, nor the literature on the nitrification process. Although NH3+ is consumed and NO2is produced by N. europaea , N. vulgaris does not consume the produced NO2and consequently no NO3conversion is achieved by the community. 4.3.3.5 OptCom As mentioned in section 4.2.2 the Python implementation within the MICOM framework was used for assessing OptCom. Hence, taxa abundancies were a mandatory input for simulations. The organism abundances predicted by SteadyCom were used (0.01 for N. europaea and 0.99 for N. vulgaris ) and a minimal growth rate of 0.0075 h-1 for both organisms was also fixed. Microbial Community Simulation Methods PhD Thesis | Sophia Santos |99 Growth and uptake/production rates differ slightly from the ones obtained using the MICOM simulation method, although are in the same range of values. Likewise, the nitrogen flow prediction does not match the previous predictions, showing that NH3+ is consumed and NO2is produced by N. europaea , but no conversion of NO2to NO3is achieved by the community. No interspecies metabolite interactions were obtained. 4.3.4 General assessment overview All simulation methods assessed showed different capabilities, exhibiting strengths and weaknesses, as shown in Figure 4.4. For instance, on one hand, SteadyCom allowed the prediction of Organisms Abundancies , while such abundancies are not predicted on any other method, being instead an obligatory input parameter for all other methods. Feature Simulation Methods FBA SMETANA SteadyCom MICOM OptCom (MICOM) Software availability Scalable for large communities Customizable Fast Individual Growth Rate Community Growth Rate Flux Rates Species Interactions Organisms Abundances Unsatisfactory Poor Satisfactory Good Outstanding Figure 4.4. Qualitative assessment of the studied steady-state simulation methods, with application to microbial communities. We evaluated each method from an unsatisfactory (red) to an outstanding performance (dark green). Chapter 4 100| PhD Thesis | Sophia Santos The evaluation of poor and satisfactory in the FBA and MICOM methods, respectively, are associated with the option of adjusting initial taxa abundances (altering the SBML file or as required input, respectively) to perform simulations. FBA, SteadyCom, and MICOM are all methods that allowed to inspect Individual and Community Growth and Flux Rates. Figure 4.5. Summary of the main output capabilities of each simulation method with application to microbial communities. SMETANA performed unsatisfactorily in Individual Growth Rate , Community Growth Rate and Flux Rates features and poor on the Customizable feature because of its interface's inability to adjust organism abundances and to provide access to growth or flux rates. On the other hand, SMETANA performed as outstanding when checking the ability to predict metabolite Species Interactions between organisms. Although FBA and SteadyCom allow inspecting intra-organisms metabolic flux rates easily, and correctly Microbial Community Simulation Methods PhD Thesis | Sophia Santos |101 predict the nitrogen flow throughout the network, no tool other than SMETANA provided insights on potential metabolite exchanges. MICOM does not allow to easily obtain exchange flux distributions. The retrieval of internal Flux Rates is not straightforward, needing a conversion step taking into account species abundancies, resulting in an evaluation of this feature of good. OptCom’s original Matlab® version was not available for download (thus evaluated as poor). However, all Python methods evaluated are freely available for download, with straightforward information for installation and scalable for large communities. As a summary of the main capabilities of simulation methods with application to microbial communities (Figure 4.5) and taking into account the strengths and weaknesses of the respective method, SMETANA should be used in case of screening a microbial community for pairwise interspecies interactions and species interaction potential, and SteadyCom should be used in case of the determination of species abundancies in a community, although in this case, no experimental data is available to validate the abundance predictions obtained. Exception made for SMETANA, all other methods can predict, to some extent, flux and growth rates. According to the results using the case study of the interaction of N. europaea and N. vulgaris , FBA was the method with prediction values closer to the experimental ones reported in the literature. Chapter 4 102| PhD Thesis | Sophia Santos 4.4 Conclusions In this work, the performances of different steady-state simulation methods with application to microbial communities were analyzed. A quantitative comparison between simulation methods was not performed due to differences in method-specific definitions of objective functions, which by itself contributes to differences in flux distributions (Santos et al., 2016; Schuetz et al., 2007). Comparing the simulation methods' performance to predict the nitrification process catalyzed by the bacteria N. europaea and N. vulgaris , pFBA, and SteadyCom seem to provide the best results. SteadyCom only overcame pFBA by providing more information on the optimum organism abundances, regarding the specific environmental conditions defined. Thus, although its formulation has been developed for predicting single-organism phenotypic behavior, pFBA seems relevant for understanding microbial community interactions. SMETANA proved to be outstanding at identifying known interactions that are characteristic of the nitrification process and potential interactions between both organisms. However, no information on growth and uptake/withdrawal rates is obtained. MICOM, and the MICOM implementation version of OptCom, revealed limitations regarding the correct prediction of interactions between the organisms in the nitrification process, as well as the range of flux distribution when an equal organism abundance was assumed. However, when more accurate organism abundances are provided, the obtained fluxes’ distribution is comparable to the pFBA and SteadyCom methods. To summarize, the available simulation methods with application to microbial communities could, to some extent, predict the phenotypic behavior that characterizes the nitrification process catalyzed by N. europaea and N. vulgaris . Each one of the simulation methods showed strengths and weaknesses and the success of either approach depends on the microbial community composition and complexity. Consequently, the use of more than one simulation method is recommended whenever possible as these showed to complement and validate each other. Systems biology approaches allow to understand how different species interact and affect their environment. However, these methods are limited to Matlab or Python, which restricts their manipulation to bioinformaticians or other expert researchers. Hence, the implementation of these simulation methods in frameworks with user-friendly interfaces such as Kbase (Arkin et al., 2018) or Optflux (I. Rocha et al., Designing Microbial Communities with MEWpy PhD Thesis | Sophia Santos |109 framework, the modification strategy (deletion, under/overexpression) as well as the target product and the environmental conditions. In the phenotype simulation layer all the constraint-based methods already used to evaluate wildtype and mutant strains such as Flux Balance Analysis (FBA) (Varma & Palsson, 1994) or Minimization of Metabolic Adjustment (MOMA) (Segrè et al., 2002), as well as Flux Variability Analysis (FVA), can be selected. The different phenotype simulation methods are provided by COBRApy (Ebrahim et al., 2013) and REFRAMED libraries. In the optimization layer, all the optimization heuristics used for strain design and respective objective functions can be found. The candidate solution fitness of each optimization iteration is evaluated by running the respective phenotype simulation required by the objective function. Here an objective function was included to minimize the molecular weight of the compounds of exchange reactions when using the minimal medium optimization. The EAs in MEWpy are implemented by the Inspyred (Tonda, 2020) and JMetalPy (Benítez-Hidalgo et al., 2019) Python libraries. MEWpy requires a compatible linear programming solver (CPLEX, GUROBI, or GLPK), with installed Python dependencies. Currently, in specific for microbial communities, MEWpy allows the optimization of 4 different scenarios from a community perspective (Figure 5.1): (i) Minimal medium optimization, including the objective function of minimizing the molecular weights ( MolecularWeight fitness evaluation function) of the exchange reactions compounds (possible when chemical formulas of compounds are included in GSMMs); (ii) Interactions of metabolites among species; (iii) Identification of reactions/enzymes/genes in the community to manipulate targeting the production of desired compounds, and (iv) Identification of organism-specific reactions/enzymes/genes to manipulate targeting the production of desired compounds. MEWpy can be installed from PyPi ( pip install mewpy ) and the source code is available at https://github.com/BioSystemsUM/mewpy under the GPL license. Chapter 5 110| PhD Thesis | Sophia Santos Figure 5.1 Summary of the different strategies to optimize and design microbial communities using MEWpy. A – Medium optimization. B – Intraspecies metabolite exchange optimization. C – Community reaction/gene optimization. D – Organismspecific reaction/gene optimization. 5.3 Usage examples To illustrate the applicability of the evolutionary algorithms to GSMMs of microbial communities, only two specific cases were tested, namely regarding Minimal Medium Optimization and Organism Specific Reaction Optimization due to the difficulty of finding suitable case studies for the other optimization scenarios. For each scenario, different case studies were evaluated. All simulations were performed using CPLEX 12.8.0 as a solver, through the PyCharm Integrated Development Environment (IDE). All scripts, input/output auxiliary files, and GSMMs used in this work are available on GitHub at SophiaSantos/DesignOptimizationMicrobialCommunities. 5.3.1 Optimization of target compound production through reaction manipulation 5.3.1.1 Case study The selected case study shown here considers a synergistic co-culture of model organisms Escherichia coli ( E. coli ) and Saccharomyces cerevisiae ( S. cerevisiae ) developed for the heterologous Designing Microbial Communities with MEWpy PhD Thesis | Sophia Santos |111 production of the high-value flavonoid naringenin from xylose, by performing further genetic manipulation in only one of the co-culture organisms, in this case, E. coli (Zhang et al., 2017). The genetic manipulation performed used gene deletions ( Δpyk and ΔpheA ) and gene over-expressions ( tktA, ppsA aroGfbr , aroE, and tyrAfbr ) as represented in Figure 5.2. Naringenin is widely produced by citrus plants, especially in grapefruit and bitter orange. The heterologous route to produce naringenin was introduced in S. cerevisiae , which is more suitable than prokaryotic organisms to express compounds produced by plants. Figure 5.2. Schematic representation of the metabolic pathway for biosynthesis of naringenin via the co-culture of S. cerevisiae and E. coli with the experimental metabolic engineering strategy used. Gene deletions of Δ pyk and Δ pheA ; Gene over-expression of tktA, ppsA aroGfbr , aroE, and tyrAfbr . Adapted from (Zhang et al., 2017). The heterologous pathway for the synthesis of naringenin from L-tyrosine is composed of four enzymes (in yellow). Abbreviations: PEP - phosphoenolpyruvate, E4P - erythrose-4-phosphate, DAHP - 3-deoxy-D-arabino-heptulosonate-7-phosphate, DHS - 3-dehydroshikimic acid, SHK - shikimic acid, CHA - chorismic acid, 4HPP - 4-hydroxyphenylpyruvic acid, L-Phe - L-phenylalanine, L-Tyr - Ltyrosine, EPSP - 5-enolpyruvylshikimate-3-phosphate, p-CA - p-coumaric acid, p-CA-CoA - p-coumaroyl-CoA, NC - naringenin chalcone. Chapter 5 112| PhD Thesis | Sophia Santos 5.3.1.2 GSMMs and parameter definition For the validation of Reaction Optimization by MEWpy, the GSMMs for E. coli (iAF1260) developed by (Feist et al., 2007) and for S. cerevisiae (iMM904) developed by (Mo et al., 2009) were used. Both models were retrieved from the BiGG database (Schellenberger et al., 2010), having the same metabolite identifiers, which is required for the construction of the community model. Biomass equations identifiers of both GSMMs were switched to match the same identifier (‘R_BIOMASS’) for community model construction purposes. The defined medium was inferred from information in the literature (Zhang et al., 2017). Xylose was set as E. coli’s carbon source and acetate produced by E. coli was set as S. cerevisiae's carbon source using reaction community identifiers. As shown in Figure 5.2 the metabolic engineering strategy in this case study includes overexpression and deletion of E. coli genes. Both strategies were tested for reactions (ROUProblem and RKOProblem), using MEWpy. Target reactions for the deletion or under-over expression optimization were defined as the reactions assigned with iAF1260 ( E. coli GSMM identifier), removing from those all the essential and exchange reactions. Although the case study objective is the production of naringenin by S. cerevisiae , all the metabolic engineering performed in E. coli had the objective to produce L – tyrosine, and for that reason, the objective of all optimizations will be the production of L – tyrosine by E. coli . All optimizations were run using 100 as the maximum number of generations and the candidate maximum size of two was used for the RKOProblem (maximum number of reaction knock-outs) and five for the ROUProblem (maximum number of over/under-expression of reactions), as shown in the case study culture (Zhang et al., 2017). A Jupyter Notebook with the whole process is available in Supplementary Material – naringenin_optimization.ipynb. Table 5.1 FVA analysis of the co-culture composed by E. coli and S. cerevisiae L-tyrosine production capabilities. The specific growth rate was set to at least 10% of the specific growth rate obtained with the pFBA simulation. L – tyrosine Production (mmol gDW-1 h-1) E. coli [-0.0572, 0.463] S. cerevisiae [-0.463, 0.0572] An FVA was performed before the optimization process to try to verify the community model's ability to produce L-tyrosine (Table 5.1). The FVA analysis showed that under the environmental conditions defined, L–tyrosine is being exchanged between the two organisms. Moreover, higher production of L– tyrosine is obtained by E. coli (0.463 mmol gDW-1 h-1). Designing Microbial Communities with MEWpy PhD Thesis | Sophia Santos |113 5.3.1.3 Reaction Deletion Optimization (RKOProblem) When performing the reaction deletion optimization (RKOProblem), although a maximum of two modifications was set, the best solutions included only 1 modification, Δpyk which is one of the two deletions tested in the case study (Figure 5.3). Figure 5.3. Schematic representation of the metabolic pathway for biosynthesis of naringenin via the co-culture of S. cerevisiae and E. coli with the results of MEWpy reaction knock-out optimization and the pFBA prediction values of the effect of deletion of the Δpyk corresponding reaction. Adapted from (Zhang et al., 2017). The heterologous pathway for the synthesis of naringenin from L-tyrosine is composed of four enzymes (in yellow). Abbreviations: PEP - phosphoenolpyruvate, E4P - erythrose-4-phosphate, DAHP - 3-deoxy-D-arabino-heptulosonate-7-phosphate, DHS - 3-dehydroshikimic acid, SHK - shikimic acid, CHA - chorismic acid, 4HPP - 4-hydroxyphenylpyruvic acid, L-Phe - L-phenylalanine, L-Tyr - Ltyrosine, EPSP - 5enolpyruvylshikimate-3-phosphate, p-CA - p-coumaric acid, p-CA-CoA - p-coumaroyl-CoA, NC - naringenin chalcone. Even in solutions with two modifications, Δpyk is included in all of them. ΔpheA appears in solutions when the number of modifications is higher than 2 jointly with the deletion of reactions in the Chapter 5 114| PhD Thesis | Sophia Santos L–serine biosynthesis and in mannose metabolism (Supplementary Material – naringenin_optimization.ipynb) but never together with Δpyk . When performing a pFBA simulation using the knock-out of pyk and pheA, the predicted value for the production of L-tyrosine by E. coli and consequent consumption by S. cerevisiae is very similar to the sole deletion of pyk (Table 5.2). In summary, the results of the reaction deletion optimization suggest that indeed Δpyk seems to have a relevant role in silico for the E. coli production of L-tyrosine under the co-culture with S. cerevisiae , showing also the MEWpy applicability on reaction deletion optimization in the context of microbial communities. 5.3.1.4 Reaction Over/Under expression Optimization (ROUProblem) When performing the reaction over/under expression optimization (ROUProblem), although a maximum of five modifications was set, the best solutions included only one modification, a 2-fold overexpression of tktA which is one of the five over-expression genes used in the case study (Figure 5.2). Even in solutions with more modifications, tktA is included in all of them jointly with over-expressions in reactions on sugar metabolism (Supplementary Material – naringenin_optimization.ipynb). Table 5.2. pFBA analysis of the co-culture composed by E. coli and S. cerevisiae for L-tyrosine production using the MEWpy deletion and/or overexpression reactions predictions. Consumed metabolite rates are represented as negative and produced metabolite rates are represented as positive. E. coli L–tyrosine Production in co-culture (mmol gDW-1 h-1) Wild-type -0.0571 Δpyk 0.0416 Δpyk and ΔpheA 0.0415 tktA 2-fold overexpression 0.0419 Δpyk and tktA 2-fold overexpression 0.0415 Indeed, when performing a pFBA simulation constraining the tktA reaction to a 2-fold flux value, Ltyrosine is predicted to be produced by E. coli (Table 5.2) in the same range as with pyk knock-out. However, when adding to the pFBA simulation constraints the knock-out of Δpyk , there is no effect on the production of L-tyrosine. Nevertheless, MEWpy over/under expression optimization showed significant Designing Microbial Communities with MEWpy PhD Thesis | Sophia Santos |115 information on which reactions to try to over/under express for this co-culture showing once more, its applicability in optimizing microbial communities. 5.3.2 Minimal Medium Optimization For the validation of the Minimal Medium Optimization by MEWpy, eight published and manually curated GSMMs of single prokaryotic and eukaryotic organisms were used that are part of one or more of the five selected natural and synthetic microbial communities (Table 5.3). Table 5.3. Microbial community case studies used for the Minimal Medium Optimization and information on the published and manually curated GSMMs used for each community case study. Genome-scale Metabolic Models Microbial Community ID Organism Reference ID Reference iAO358 Lactococcus lactis (Oliveira et al., 2005) 1 (Ponomarova et al., 2017) iBT721 Lactobacillus plantarum (Teusink et al., 2006) iMM904 Saccharomyces cerevisiae (Mo et al., 2009) iBif452 Bifidobacterium adolescentis (El-Semman et al., 2014) 2 (El-Semman et al., 2014) iFap484 Faecalibacterium prausnitzii (El-Semman et al., 2014) iJN1411 Pseudomonas putida (Nogales et al., 2017) 3 (Sgobba et al., 2018) iAF1260 Escherichia coli (Feist et al., 2007) iYS432 Corynebacterium glutamicum (Shinfuku et al., 2009) 4 (Sgobba et al., 2018) iAF1260 Escherichia coli (Feist et al., 2007) iAF1260 Escherichia coli (Feist et al., 2007) 5 (Zhang et al., 2017) iMM904 Saccharomyces cerevisiae (Mo et al., 2009) GSMMs were screened for the same external metabolite identifiers, and the BiGG (Schellenberger et al., 2010) annotation was adopted when different annotations were found. Biomass equations Chapter 5 116| PhD Thesis | Sophia Santos identifiers of all GSMMs were switched to match the same identifier (‘R_BIOMASS’) for community model construction purposes, and compound chemical formulas were included in all GSMMs that lacked that information. The experimentally defined medium of each microbial community was determined by information in the literature (Table 5.3) - Supplementary Material – Table S1. Figure 5.4. Number of exchange reactions defined on the Experimental medium (green line), on the optimized medium by MEWpy (blue line) without using Molecular Weight minimization, and the number of exchange reactions that are common in both mediums (pink line). For the optimization of the minimal medium for each microbial community, a carbon source, defined in the literature, was set as a constraint, as well as all exchange reactions that are essential for each community GSMM have growth. Minimal Medium Optimization was tested with and without minimizing exchange compound Molecular Weight. The best results were compared and validated with Designing Microbial Communities with MEWpy PhD Thesis | Sophia Santos |117 experimental data from literature (Figure 5.4) - Supplementary Material – Table S2, Table S3, Table S4, Table S5, Table S6. When analyzing the obtained results, major differences are exhibited in communities ID 1 (coculture Lactococcus lactis , Lactobacillus plantarum, and S. cerevisiae ) and ID 2 (co-culture Bifidobacterium adolescentis and Faecalibacterium prausnitzii ). Indeed, in these two cases, a rich medium is used in the experimental setup, which resulted in an Experimental defined medium with a high number of metabolites. In community ID 1 most of the Experimental medium components are ions, vitamins, and cofactors, as well as amino acids which are essential to Lactococcus lactis. Whereas in community ID 2 Experimental medium components are mainly amino acids, which are an important part of the gut ecosystem. In both cases, the majority of the amino acids in the Experimental medium were not included in the Optimized medium , indicating that individual amino acid auxotrophies are being bridged by other organisms within the community. Importantly, a large part of the Optimized medium components is similar to the Experimental medium . Differences between experimental and predicted mediums are higher in the community ID 3 composed of Pseudomonas putida and E. coli (Sgobba et al., 2018), where the predictions revealed that vitamins (for instance, riboflavin, thymidine, and pantothenate) are not required for the community growth. These differences cannot be directly translated for the definition of a new experimental minimal medium because most of them are associated with the lack of metabolite requirement information included in the GSMM. In all other cases, the differences happen because MEWpy often selects complex molecules that can, at the same time, be nitrogen, phosphate, and/or sulfur sources (e. g. guanosine 3phosphate, adenosine 3, 5–biphosphate or L-cysteine), as the objective is to determine a minimal number of exchange reactions that allow the growth of a specific community. Also, the best solutions generally differed among them on the nitrogen, phosphate, sulfur, or iron sources when the GSMMs included more than one suitable replacement (for instance, sulfate, thiosulfate, hydrogen sulfide, or L-cysteine). These results suggest that in cases where the exchange metabolite chemical formula is available in the GSMM, the minimization of the sum of the Molecular Weight of the metabolites in the solution would be helpful to eliminate mainly complex molecules that can function as nitrogen, phosphate, and/or sulfur sources. This hypothesis was corroborated in part using the minimization of exchange metabolite Molecular Weight, which eliminated the best solutions that had complex molecules, such as the replacement of Chapter 5 118| PhD Thesis | Sophia Santos sulfate (molecular weight of 96.06 g mol-1) for hydrogen sulfide (molecular weight of 34.08 g mol-1). However, in none of the case studies the number of metabolites in the best solutions changed. Overall, the results obtained for Minimal Medium Optimization using MEWpy showed biological significance, correctly predicting a minimal medium for each community. However, results should be carefully analyzed considering differences between the information included in each organism GSMM within the community and the experimental data available. Modeling and Design of Microbial Communities from Extremophilic Environments in the Azores PhD Thesis | Sophia Santos |125 reversibility of reactions workflow, PSORTb 3.0 (Yu et al., 2010) was used to predict reactions compartments, the Transporters Classification Database (TCDB) (Saier et al., 2021) was used to predict the model transport reactions through merlin ’s plug-in TranSyT. Table 6.1. Genome’s information used to import genome files into merlin . Organism Taxonomy ID Assembly Accession Number Desulfurococcus amylolyticus DSM 16532 NCBI: txid768672 ASM23101v3 Pyrobaculum aerophilum str. IM2 NCBI: txid178306 ASM722v1 Thermodesulfovibrio yellowstonii DSM 11347 NCBI: txid289376 ASM2098v1 Thermofilum adornatus 1910b NCBI: txid697581 ASM81324v1 Thermoplasma acidophilum DSM 1728 NCBI: txid273075 ASM19591v1 Thermus scotoductus SA-01 NCBI: txid743525 ASM18700v1 6.2.3 Metabolic Models Reconstruction merlin (Capela et al., 2021) is a user-friendly framework that allows performing several steps of the reconstruction process semi-automatically, downloading relevant information from several databases, and was used to assist in the reconstruction of all GSMMs. Moreover, it has a graphical interface that facilitates GSMM information reviewing and manual curation. The main steps of each GSMM model reconstruction process are hereafter described. 6.2.3.1 Genome Annotation merlin allows performing the functional annotation of a genome, using as similarity search engines the Basic Local Alignment Search Tool (BLAST) (Altschul et al., 1990) and Diamond (Buchfink et al., 2021) against databases that contain reviewed (such as UniProt/Swiss-Prot) and unreviewed (such as UniProt/TrEMBL) enzyme information. The EC numbers and enzymatic functions assigned to each gene Chapter 6 126| PhD Thesis | Sophia Santos are scored based on the taxonomy and frequency of similar sequences, as described elsewhere (Dias et al., 2018). The genome annotation is a crucial step of the GSMM reconstruction process, as incorrect EC numbers and enzymatic function assignments can significantly impact the model performance. Once the similarity search is complete, merlin ’s automatic annotation workflow feature (Capela et al., 2021) can prioritize the obtained gene products and EC numbers annotations. This operation considers a list of organisms ordered by phylogenetic similarity provided by the user and defines a confidence level (A to I) for each gene annotation when a match is found. For all organisms within this study , due to the lack of available information on closely related organisms and to maximize retrieval of reviewed information, the automatic workflow feature was configured to use genus instead of species as input. The ranked list of the closely related phylogenetic genus of all organisms and the associated phylogenetic trees (Supplementary Material) were constructed from 16S RNA sequences reference organisms of each genus using the EMBL-EBI Clustal OMEGA multiple sequence alignment tool (Sievers et al., 2011). 6.2.3.2 Assembling the Metabolic Network The assembly of a metabolic network starts with gathering all reactions present in the organism. merlin retrieves reactions by importing them from KEGG, based on the annotated EC numbers from the previous step and spontaneous reactions. Reactions should also be balanced (generic and metabolites without formula must be curated), and reversibility must be confirmed to avoid mispredictions of the model. To assist these steps, merlin checks whether a reaction is balanced and also corrects the reactions’ reversibility. Although these steps are automated in merlin , we performed manual curation based on curated information from literature or databases, such as MetaCyc (Caspi et al., 2014), as this step is strongly encouraged. Compartmentalization The compartmentalization of each model was based on results obtained from PSORTb 3.0 (Yu et al., 2010) choosing the specific organism type and gram stain before online submission. The “Long format” report generated was imported, and merlin ’s compartments feature assigned each reaction to its specific compartment (Capela et al., 2021). Modeling and Design of Microbial Communities from Extremophilic Environments in the Azores PhD Thesis | Sophia Santos |127 Transport Reactions After compartmentalizing the model, we defined transport reactions between compartments and also with the exterior. merlin ’s TranSyT (Lagoa et al., 2021) was used to generate system-specific transport reactions associated with Gene-Protein-Reactions (GPRs) rules that are automatically integrated into each model. Genes, proteins, and reactions A high-quality GSMM requires GPRs’ rules to predict genetic modifications accurately. These associations are usually defined according to databases and literature (Rocha et al., 2008; Thiele et al., 2010). merlin ’s “Gene-Protein-Reaction rules” feature automatically adds GPRs rules to the model. The algorithm used by merlin to implement these rules is described elsewhere (Dias et al., 2015). 6.2.3.3 Converting the Metabolic Network to a Stoichiometric Model Biomass equation Biomass composition must be experimentally determined in cells growing in the log phase before being included in the model. However, in the absence of experimental information organism-specific, data from genome information (particularly nucleotides, deoxynucleotides, and amino acids) can be used or adapted from phylogenetically related organisms. The importance of an accurate biomass composition determination has been reported (Santos et al., 2016), showing that even minor differences in biomass component coefficients can significantly impact simulations’ predictions. Therefore, a well-defined biomass equation is crucial for the GSMM reconstruction process. Information regarding the biomass composition of the organisms present in this study is scarce in the literature. Hence, for organisms of the domain Bacteria, the macromolecular composition was adapted from the gram-negative bacterium Escherichia coli ( E. coli ) (Feist et al., 2007) or the grampositive bacterium Bacillus subtilis ( B. subtilis ) (Dauner & Sauer, 2001). For organisms for the domain Archaea, the macromolecular composition was adapted from the archaeon Methanosarcina fusaro ( M. fusaro ) (Goyal et al., 2014). The composition of amino acids, nucleotides, and deoxynucleotides was estimated from genome information through merlin ’s “e-Biomass Equation” feature. This feature also Chapter 6 128| PhD Thesis | Sophia Santos automatically includes cofactor composition based on a study of universal essential cofactors in prokaryotes (Xavier et al., 2017). The fatty acids composition was taken from experimental data from the literature for all organisms. The lipid, carbohydrates, and cell wall composition (when required) compositions were adapted from the information of reference organisms considering enzyme annotation to include or exclude specific elements. When required, new coefficients were calculated, maintaining the relative abundances of the original data. Alternative biomass equations were defined for simulations under aerobic and anaerobic conditions. Since in the latter environmental conditions, Heme is not required, the compound was removed from the Cofactor composition, and all other coefficients were recalculated as mentioned before. Growth and maintenance ATP requirements No information was found for any of these organisms for growth and ATP maintenance. Therefore, such data was adapted from experimental data for E. coli (Feist et al., 2007; Neidhardt et al., 1990). 6.2.3.4 Metabolic Model Curation and Validation The model curation is an iterative process that stops when simulation results match experimental data in the literature. merlin ’s interface allows the user to efficiently perform re-annotations, correct reactions stoichiometric balance and directionality, include or exclude reactions from the model, and finally, export models in Systems Biology Markup Language (SMBL) format to be used in simulation using platforms such as Optflux (I. Rocha et al., 2010), Matlab®’s COBRA toolbox (Heirendt et al., 2019), COBRApy (Ebrahim et al., 2013) or REFRAMED. Gap-filling Before being ready for simulations, the metabolic network must be screened for possible gaps. The presence of gaps can deeply compromise the synthesis of biomass components and other relevant compounds. To assist in this process, merlin ’s “BioISO” (Cruz et al., 2021) helps trace back the network to identify gaps that can be originated from errors in genome annotation, absence of enzymatic, transport, or exchange reactions, or incorrect reaction irreversibility or direction. Other features included in merlin , Modeling and Design of Microbial Communities from Extremophilic Environments in the Azores PhD Thesis | Sophia Santos |129 such “find Blocked reactions”, which identifies reactions that contain dead-end metabolites, and “Draw in Browser” which opens on a web browser a selected KEGG pathway, showing specifically highlighted enzymes and reactions present in the model, facilitate the detection of gaps in the network (Capela et al., 2021). Literature and databases (KEGG, MetaCyc, BRENDA, for instance) should be used to assist the gap-filling process. This is an iterative process, repeated until all biomass precursors, and other essential compounds, can be synthesized, and a feasible model is obtained. Aerobic and anaerobic metabolism and Carbon Source utilization Experimental data from the literature was consulted to identify the growth conditions for each organism. A defined minimal medium was favored when available. Whenever an organism was reported to have aerobic and anaerobic growth, both conditions were tested. Growth under different carbon sources was also tested and the presence of the specific transport reactions was screened. 6.2.3.5 Simulations and Strain Optimizations Software SMETANA (Zelezniak et al., 2015) was used to predict pairwise interactions. All other simulations and community strain optimizations were performed using the Python package MEWpy (Pereira et al., 2021) version 0.1.28, CPLEX 12.8.0 as a solver, through the PyCharm Integrated Development Environment (IDE). MEWpy allows for phenotype simulation using constraint-based methods provided by COBRApy (Ebrahim et al., 2013) and REFRAMED libraries, such as Flux Balance Analysis (FBA) (Varma & Palsson, 1994), or SteadyCom (Chan et al., 2017), specific for simulating microbial communities. Quantitative evaluation of the model was performed using Flux Variability Analysis (FVA) (Mahadevan et al., 2003) to determine S. azorense’s cellulose production capabilities in co-culture. The analysis included setting the specific growth rate to at least 10% of the specific growth rate obtained with parsimonious FBA (pFBA) (Lewis et al., 2010) in the respective reference flux distribution. All Python scripts, input/output auxiliary files, and GSMMs used in this work are available on GitHub at SophiaSantos/DesignOptimizationMicrobialCommunities. Chapter 6 130| PhD Thesis | Sophia Santos 6.3 Results and Discussion 6.3.1 Samples Characterization Taxonomic profiling of all samples from hydrothermal vents in São Miguel, Azores was performed using assembly-based and read-based analysis algorithms (Chapter 2). Joint information from both approaches revealed the potential presence of 12 different microorganisms. Several of these organisms were predicted to be present in various samples. Physiological and metabolic characteristics of all identified organisms were screened to determine if their presence in such environments is known or viable, based on each sample site's physiochemical properties. Two of the five samples, FCRG and PCRG, were excluded from this characterization, as both are composed of only two organisms present in the CV sample and therefore their interactions will be screened in that sample. From all organisms identified through the profiling algorithms nine were selected for the GSMM reconstruction process: Acidithiobacillus caldus SM-1 Acidithiobacillus caldus SM-1 ( A. caldus ) is a gram-negative, moderately thermophilic, acidophilic (pH range 1.0 - 3.5), chemolithoautotrophic bacterium. Growth temperature ranges between 32 and 52ºC. Mixotrophic growth is obtained with tetrathionate and glucose or yeast extract (Hallberg et al., 1994). A. caldus is capable of oxidizing elemental sulfur and other reduced inorganic sulfur compounds (Hallberg et al., 1996) and therefore is one of the dominant sulfur-oxidizing bacteria in bioleaching reactors together with iron-oxidizing bacteria (L. Chen et al., 2012; Watkin et al., 2009; Watling et al., 2014). The GSMM of this organism was reconstructed as part of the master thesis entitled “Genomescale metabolic modeling of an extremophile microbial community” by Rui Barros Nunes. Acidimicrobium ferrooxidans DSM 10331 Acidimicrobium ferrooxidans DSM 10331 ( A. ferrooxidans ) is a gram-positive, moderately thermophilic, acidophilic bacterium. Growth has been observed between temperatures of 45 and 50ºC of temperature and 1 and 3.5 pH values. Autotrophic growth is possible on ferrous iron and heterotrophic growth is possible on yeast extract (Cleaver et al., 2007). A. ferrooxidans is usually found in mixed cultures of thermophilic microorganisms in bioleaching processes, to enhance metal extraction in heaps (Watling Modeling and Design of Microbial Communities from Extremophilic Environments in the Azores PhD Thesis | Sophia Santos |131 et al., 2014). The GSMM of this organism was reconstructed as part of the master thesis entitled “Genome-scale metabolic modeling of an extremophile microbial community” by Rui Barros Nunes. Desulfurococcus amylolyticus DSM 16532 Desulfurococcus amylolyticus DSM 16532 ( D. amylolyticus ) is an obligately anaerobic and hyperthermophilic archaeon. Growth is optimal at pH values between 6.0 and 6.5 and temperature values between 80 and 92ºC (Perevalova et al., 2016). Can grow on a broad range of carbon sources (sugars, polysaccharides, and amino acids) but is one of the few Crenarchaeota able to grow in cellulose (Reischl et al., 2018). Pyrobaculum aerophilum str. IM2 Pyrobaculum aerophilum str. IM2 ( P. aerophilum ), in contrast with most species of the genus Pyrobaculum, is an aerobic, hyperthermophilic archaeon. Organic and inorganic compounds serve as substrates under both aerobic and anaerobic respiration. Growth temperatures range between 75 and 105ºC and pH values between 5.8 and 9.0. Autotrophic growth is achieved by the oxidation of hydrogen or thiosulfate (Volkl et al., 1993). Sulfurihydrogenibium azorense Az-Fu1 Sulfurihydrogenibium azorense Az-Fu1 ( S. azorense ) is a gram-negative, thermophilic, chemolithoautotrophic, and microaerophilic bacterium (Lalonde et al., 2005). The bacterium grows optimally at 68ºC, pH 6, and at low concentrations of NaCl and can also grow heterotrophically (Nakagawa et al., 2005) and use elemental sulfur, thiosulfate, hydrogen, and ferrous iron as energy sources (Aguiar et al., 2004). As an Aquificales member, S. azorense fixes CO2 via reductive Tricarboxylic Acid Cycle (rTCA) (Hügler et al., 2007) to generate acetyl-CoA as an end product (Gupta et al., 2013). The GSMM of this organism was reconstructed in Chapter 3 of this thesis. Chapter 6 132| PhD Thesis | Sophia Santos Thermodesulfovibrio yellowstonii DSM 11347 Thermodesulfovibrio yellowstonii DSM 11347 ( T. yellowstonii ) is a thermophilic sulfate-reducing, strictly anaerobic gram-negative bacterium. Growth was observed for temperatures between 40ºC and 70ºC and pH between 5.5 and 8.5. T. yellowstonii can use sulfate, thiosulfate, and sulfite terminal electron acceptors (Bhatnagar et al., 2015). In the presence of sulfate, growth was only observed with lactate, pyruvate, hydrogen plus acetate, or formate plus acetate. Pyruvate and lactate are oxidized to acetate (Henry et al., 1994). Thermofilum adornatus strain 1910b Thermofilum adornatus strain 1910b ( T. adornatus ) is a strict anaerobe archaeon. As moderately acidophilic and hyperthermophilic crenarchaea, T. adornatus requires various components of cells or culture broths of other Crenarchaeota, specifically culture broth filtrate of Desulfurococcus and Pyrobaculum . Growth was observed for temperatures between 50ºC and 95ºC, and pH between 5.3 and 8.5. Peptone, yeast extract, cellulose, starch, glucose, lactose, mannose, and pyruvate can be used as carbon sources. However, T. adornatus can utilize cellulose as the sole carbon and energy source (Zayulina et al., 2020). Thermoplasma acidophilum DSM 1728 Thermoplasma acidophilum DSM 1728 ( T. acidophilum ) is a thermophilic, acidophilic archaeon lacking a cell wall. Growth takes place over a range of temperatures from 45ºC to 62ºC and a pH between 1 and 3.5 (Darland et al., 1970). Anaerobic growth is highly supported by the addition of sulfur which results in the accumulation of thiosulfate. Under fully aerobic conditions, sulfur compounds show no influence on growth. Yeast extract is required for growth in all conditions (Segerer et al., 1988). Thermus scotoductus SA-01 Thermus scotoductus SA-01 ( T. scotoductus ) is a gram-negative aerobic thermophilic bacterium, which can grow at temperatures between 42ºC and 73ºC and pH values between 5.0 and 11. Is a Modeling and Design of Microbial Communities from Extremophilic Environments in the Azores PhD Thesis | Sophia Santos |133 facultative anaerobe, reducing nitrate to nitrite anaerobically. Growth was observed in proline, glutamate, pyruvate, and sucrose (Kristjánsson et al., 1994). 6.3.1.1 Caldeira Velha (CV) The metagenomic approaches used to determine the microbial composition predicted that nine microorganisms were present in this sample. However, there was experimental evidence of the presence of the genus (Sahm et al., 2013) of only seven of such microorganisms, which were therefore selected to proceed to the GSMM reconstruction process (Table 6.2). Nevertheless, it should be noted that inconsistencies in terms of optimal growth temperatures and pH values can be spotted concerning the sample global physicochemical properties for some of the predicted organisms. Table 6.2. Caldeira Velha microbial community composition and organism physicochemical properties. Organisms’ abundances were recalculated maintaining the relative abundances of the original data. Organism Domain Temperature (ºC) pH Abundance Ambiental conditions 98 2.23 S. azorense Az-Fu1 Bacteria (Gram -) 50 - 73°C 5.8 - 9.0 61.4% T. adornatus 1910b Archaea 70 - 90°C 5.5 - 7.0 10.4% T. yellowstonii DSM 11347 Bacteria (Gram -) 40 - 70°C 5.5 - 8.5 7.6% A. caldus SM-1 Bacteria (Gram -) 32 - 52°C 1.0 - 3.5 7.3% A. ferrooxidans DSM 10331 Bacteria (Gram +) 45 - 50°C 1.0 - 3.5 5.1% T. acidophilum DSM 1728 Archaea 45 - 63°C 0.8 - 4.0 5.0% P. aerophilum str. IM2 Archaea 75 - 104°C 5.0 - 7.0 3.1% The taxonomic profiling algorithms varied only on the assignment of different members of the Pyrobaculum genus ( Pyrobaculum aerophlilum str. IM2 and Pyrobaculum islandicus DSM 4184). These differences can be justified by the different reference genome databases in each approach. However, Pyrobaculum aerophlilum str. IM2, besides having more growth and metabolic information available (188 Chapter 6 134| PhD Thesis | Sophia Santos versus 64 hits on the Web of Science database), has more predicted genes (2518 versus 2043) and therefore it was selected for the GSMM reconstruction process. Thermorodis peleae and Thiomonas intermedia K12 were excluded from the study, as there is no evidence of neither being present in such environments (Sahm et al., 2013). 6.3.1.2 Nascente Poente (NP) The taxonomic profiling approaches predicted the presence of five microorganisms for this sample. Among them, are two different Thermus strains, Thermus antranikianii DSM 12462 and Thermus scotoductus SA-01. The two strains are phylogenetically very close (Chung et al., 2000; Lapierre et al., 2006), and therefore only one was selected for the GSMM reconstruction process. Thermus scotoductus SA-01 was selected (Table 6.3) to proceed with this study because it is reported to be one of the Thermus species isolated in the hot springs of the Azores (Santos et al., 1989). Table 6.3. Nascente Poente microbial community composition and organism physicochemical properties. Organisms’ abundances were recalculated maintaining the relative abundances of the original data. Organism Domain Temperature (ºC) pH Abundance Ambiental conditions 95.8 6.88 P. aerophilum str. IM2 Archaea 75 - 104°C 5.0 - 7.0 79.9% T. scotoductus SA-01 Bacteria (Gram -) 65 - 80°C 5.0 - 11 17.7% D. amylolyticus DSM 16532 Archaea 68 - 97°C 5.7 - 7.5 1.28% T. adornatus 1910b Archaea 70 - 90°C 5.5 - 7.0 1.14% Note that T. adornatus is characterized by a growth dependence on various components of cells or culture broths of other Crenarchaeota, specifically, culture broth filtrate of Desulfurococcus and Pyrobaculum species (Zayulina et al., 2020), which curiously were also predicted as present in this samples.