Spatial heterogeneity promotes antagonistic evolutionary scenarios in microbial community explained by ecological stratification: a simulation study

Ecological Modelling 399 (2019) 66–76

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Spatial heterogeneity promotes antagonistic evolutionary scenarios in microbial community explained by ecological stratification: a simulation study

T

⁎

Alexandra I. Klimenkoa,b, , Yury G. Matushkina,b, Nikolay A. Kolchanova,b, Sergey A. Lashina,b a b

Institute of Cytology and Genetics, Siberian Branch of Russian Academy of Science, Lavrentiev avenue 10, Novosibirsk, 630090, Russian Federation Novosibirsk State University, Pirogova st. 1, Novosibirsk, 630090, Russian Federation

A R T I C LE I N FO

A B S T R A C T

Keywords: Microbial communities Individual-based modelling Ecological modelling Evolutionary modelling Spatial heterogeneity Stratification

There are two evolutionary trends in genome organization among microbes: towards either amplification or reduction. Which evolutionary scenario overcomes depends on environmental conditions and the complexity of gene networks determining phenotypic traits such as metabolic features of cells. In this simulation study, we have shown that the habitats characterized by nutrient gradients allow spatial subdivision of evolutionary trends depending on the distance to the nutrient source. We have considered interrelations between cell motility, metabolic complexity of dominant populations and ecological features of developing communities and have shown that the distribution of local dominant ecogroups follows clear patterns in chemotaxis-on and -off cases. Chemotaxis was shown to be a factor impeding introduction of new forms and decreasing total biomass of the community. Our simulations have shown that ecological patterns of self-organization of microbial communities cause sustainable different strategies underlying antagonistic evolutionary scenarios.

1. Introduction Microorganisms build up the world of large communities and compact genomes. However, depending on the environmental conditions, the evolution of prokaryotic genomes might result in genome reduction or genome amplification. Gene loss in vitro in bacteria (Mira et al., 2001), genome reduction in parasites (Moran, 2002) and genome size optimization under syntrophy (Kunin and Ouzounis, 2003; McCutcheon and Moran, 2011) exemplify the trend for reduction. On the other hand, the acquisition of antibiotic resistance genes (Rodrıguez-Martinez et al., 2006; Händel et al., 2016), acquisition of genes essential for surviving in extreme conditions (Nelson et al., 1999; Nesbo et al., 2001), and symbiogenesis (Margulis, 1970; Archibald, 2011; Cavalier-Smith, 2013; Provorov et al., 2016) present the trend towards genome amplification. Active processes of gene material exchange occur in microbial communities providing microorganisms with a capacity to adapt to a broad spectrum of environmental conditions without the necessity to include all the metabolic “toolboxes” in every particular genome (Barraclough et al., 2012; Niehus et al., 2015; Soucy et al., 2015). Therefore, if an environment changes, the organisms that have acquired

the adaptive genes might rapidly get an advantage over other species and then get rid of extra genes when environmental conditions alter again (Maslov et al., 2009; Pang and Maslov, 2011; Polz et al., 2013; Langridge et al., 2015). Generally, such processes involve so called accessory genome (Gogarten and Townsend, 2005). The major part of accessory genes with known function in Escherichia coli are known to be responsible for metabolism (Soucy et al., 2015; McInerney et al., 2017). However, a large genome size implies concomitant maintenance caused on the one hand by the increased time and resources spent on genome replication and on the other hand by increased probability of acquiring lethal mutations (Eigen, 1971; Weiße et al., 2015). That leads towards so-called Black Queen dynamics: certain costly biological functions (such as detoxification of hydrogen peroxide) tend to concentrate in a small subset of individuals in the community which provide this function for others as a “public good”. Consequently, this dynamics contribute to the development of symbiotrophic communities of mutually interdependent microbial species (Morris et al., 2012; Fullmer et al., 2015). Overall, the above factors constrain genome expansion in prokaryotes providing a stimulus towards reductionist evolutionary scenarios. Another factor of evolution determining the specifics of struggle for

⁎ Corresponding author at: Institute of Cytology and Genetics, Siberian Branch of Russian Academy of Science, Lavrentiev avenue 10, Novosibirsk, 630090, Russian Federation. E-mail address: [email protected] (A.I. Klimenko).

https://doi.org/10.1016/j.ecolmodel.2019.02.007 Received 20 September 2018; Received in revised form 17 January 2019; Accepted 19 February 2019 Available online 16 March 2019 0304-3800/ © 2019 Elsevier B.V. All rights reserved.

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

models) in our previous studies (Lashin et al., 2012). It was shown that in pessimal environmental conditions characterized by low inflow of nonspecific nutrient, populations of cells with “metabolically complete” or “metabolically rich” genomes (i.e. capable of utilizing all or almost all the nutrients that are present in the modelled environment respectively) prevail in the community. On the other hand, it was also demonstrated that suboptimal and optimal conditions manifest a steady trend towards genome reduction. The most primitive populations capable of only two metabolic activities (utilization of one nonspecific and one specific nutrients) displaced others. In this study we have applied the «Haploid evolutionary constructor» software package (HEC, for more details see Methods section) (Klimenko et al., 2015, 2016) to investigate the impact of spatial gradients of nutrients and environment structure on genome amplification and reduction trends in communities of microorganisms inhabiting a flow-through environment. The aim of the study is to gain insights out of modelling and to elucidate how these spatial gradients of ecological factors and cells’ motility affect the ecological and population structures of developing communities and evolutionary scenarios in local habitats. We turned on and off the ability of cells to migrate adaptively via chemotaxis. The following problems have been addressed:

existence in microbial ecosystems is adaptive migrations of microorganisms (Niehus et al., 2015; Bridier et al., 2017). Changes of environmental conditions make a population of microbes to face the following adaptive crossroads: the first option is the phenotypical adaptation (e.g. via adjusting its metabolism to alternative energy sources), the second is undergoing adaptive evolution that favors more efficient utilization of the available nutrient reserves, or the third way is to migrate pursuing optimal environmental conditions, intruding new habitats and engaging in competition with local species (Callahan et al., 2014; Chew et al., 2014; Liu et al., 2016). The impact of such adaptive individual behavior on a developmental scenario of the whole community is of current interest for acquiring more profound comprehension of interactions between various evolutionary mechanisms. There is a number of classifications of living organisms based on ecologically significant traits. The majority of them relies upon the ecological amplitude of the species projected on a certain ecological factor. For instance, microorganisms fall into the following ecophysiological groups (Zavarzin and Kolotilova, 2001) based on their preferred temperature of the environment: psychrophiles, mesophiles, thermophiles, extreme thermophiles and hyperthermophiles; based on their preferred pH: acidophiles (2 < pH < 6), neutrophiles (6 < pH < 8) and alkaliphiles (8.5 < pH < 11); based on oxygen conditions: aerobes and anaerobes, etc. Another classification is introduced based on the analysis of energy fluxes in an ecosystem proposing such categories as producers, first- second- and higher order consumers and decomposers (Pomeroy et al., 1988). There is another prominent classification based on pairwise interactions between species: commensalism, amensalism, mutualism, competition and neutralism (Lidicker, 1979; Faust and Raes, 2012). However, besides these classifications there is a typology based on ecological function that living organisms perform in their communities. For example, Henriette Poplavskaya (Poplavskaya, 1924) and Vladimir Sukachev (Sukachev, 1928) developed such a classification for plant communities. They proposed such phytosocial types as aedificators and assectators based on the role that species act upon in the community structure and its succession. The former act as community builders – either in natural conditions (autochthonous aedificators), or as a result of perturbations of the community caused by a short-term disturbance by a human being or an animal (degressive aedificators). The latter, on the other hand, act as community co-builders being, however, insignificant for its phytosocial milieu. Leonty Ramensky (Ramenskyi, 1935) and John Grime (Grime, 1974) independently proposed a system (triangular model of primary plant strategies, or CSR theory) classifying plants based on their adaptive strategies. The following two factors were regarded as foundational – stress intensity that restricts photosynthetic production and disturbance intensity that breaches the phytocenosis integrity caused by either biotic or abiotic factors. Three possible strategies determine such three major phytocenotypes as highly productive competitors, robust stress tolerators and ruderals, which are noncompetitive though rapidly colonizing recently disturbed areas in phytocenoses, and there is a number of intermediate types as well. The CSR theory has been recently adapted to microbial functional traits for methane oxidizing bacteria and in a generalized manner (Ho et al., 2013; Krause et al., 2014). Agent-based modelling (ABM) is an approach that has been extensively used for modelling of ecological (Grimm et al., 2017) and evolutionary (Liu et al., 2017) processes in the recent decade. ABM allows to develop flexible models taking into account multiple factors and individual interactions between the agents in the system. Applying agent-based modelling to microbial ecosystems shed the light on the processes of biofilm self-assembly, chemotaxis and other aspects of their functioning (Wimpenny and Colasanti, 1997; Kreft et al., 1998; Picioreanu et al., 2004; Emonet et al., 2005; Bihary et al., 2012; Ravikrishnan and Raman, 2018). The models of genome amplification and reduction in various ecological conditions were investigated for well-mixed environments (0D

• Is there any correlation of ecological conditions and metabolic complexity/genome size of the dominant population? • How does the chemotaxis ability affect species richness and total biomass of the community? • How does the chemotaxis ability influence the development of the ecological structure in evolving communities?

2. Methods 2.1. Model The models of functioning and evolution of microbial communities were simulated using the «Haploid evolutionary constructor» (HEC) software package (Lashin et al., 2007; Klimenko et al., 2015), which is based on the agent-based approach and super-individual concept (Scheffer et al., 1995). HEC meets the requirements for the next-generation ecological modelling (Grimm and Berger, 2016) since it is based on the first principles linking behavior and physiology to ecology, adopts trait-based approach to describe microbial functional types, adds microevolution, enables describing heterogeneous and dynamic abiotic environments and obtains modular multi-layer structure. The key simulation unit in HEC is a super-individual-like agent, which describes a genetically homogeneous population of cells that interacts with its environment and other populations that are genetically different from each other. All the present populations as a group constitute a microbial community. The following set of objects and processes is considered in this simulation approach. Objects:

• Environment. An aquatic flow-through environment is described •

67

by a lattice of well-mixed sub-habitats (nodes) populated with cells and substances; Nutrients. We regard any substance affecting the cell growth as a nutrient. Cells can either secrete nutrients that were synthesized inside them (in this case, we call such nutrients specific ones) or consume nutrients and use them for cell’s growth and reproduction. Those nutrients that come into the environment via the flow are called nonspecific ones. Nonspecific nutrients describe external energy source that nourishes microbial ecosystems that corresponds to dissolved organic matter, glucose or some other kind of substrate that is widely, nonspecifically used by members of microbial community for their growth and maintenance. The specific nutrients can either stimulate population growth or suppress it and even cause the

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

natural streams and brooks illustrate 1D-spb environments (Yao and Le Maguer, 1996; Gruber et al., 2011). The tubes are supplied with a nonspecific nutrient (which is consumed by all cells) from one side (source) – the «beginning» of the tube (Fig. 1c, left side of both tubes). The nonspecific nutrient as well as specific ones produced by cells of the community distribute further by diffusion and flow. The flow also transports the cells through the tube from source to sink. In the «end» of the tube (sink) nutrients and cells wash out of the environment (Fig. 1c, right side). Besides passive transport, we also consider chemotaxis ability of the cells, i.e. adaptive migrations towards more favorable environmental conditions in terms of nutrient saturation. We examined both cases of the community consisting of motile (chemotaxis-on case) and nonmotile (chemotaxis-off case) cells. Hereafter we call the community that has established at the end of a simulation run as the developed community.

drastic reduction of its population size;

• Populations. A population is characterized by a set of metabolic systems (model genes) corresponding to metabolic functions that determine which nutrients and how rapidly the cells of the population are able to consume and produce. Processes:

• Consumption of nutrients captured from the environment and secretion of products of cell metabolism; Nutrient exchange among the populations in the habitat mediated • by the fluid environment; • Reproduction in populations and natural selection; • Mutations. Mutations are described as a change of allele value of a certain gene in one or several cells of the population; • Gene loss. Cells might lose some genes generating a new population possessing a genome that is different from the original one; • Horizontal gene transfer. Cells that belong to different populations

2.1.1. Model M1: deterministic model of a «complete» microbial community We define the deterministic model of a «complete» microbial community exposed to the natural selection as model M1. A «complete» community according to our definition is such a community that consists of populations possessing all feasible metabolic combinations for the system with three specific and one nonspecific nutrient (see Fig. 1d). Initially those populations had the equal abundances (one million cells for each population, 105 cells per node of the tube for each population). Initial concentrations of nonspecific and specific nutrients (N1, S1, S2, S3) were set to 10−6 M in each node of the tube. We have examined the impact of cells’ motility on the final ecological structure that develops as a result of development of communities of a particular type. Thereby we have found typical distributions of ecogroup abundances independent of such stochastic factors as HGT/gene loss time and place of origin and rather determined by the fitness of the populations constituting corresponding ecogroups.

are able to exchange genes involved in metabolism functioning, resulting in an emergence of new populations.

These processes contributing into microbial population dynamics were described using the following equation (Lashin et al., 2012):

→ → ΔP = F (n 0 , S , r0, C , P , G ) =

r0 n 0 (P )

∑

ci si (P ) − kdeath P 2

i ∈ Iconsumed

−b

(G / k1 )n P (G / k1 )n − G n + 1

(1)

where ΔP is the change of population size (abundance, biomass) per iteration; n 0 – an amount of nonspecific nutrient consumed by the cells → of the population from the environment; S – numerical vector corre→ sponding to amounts of the specific nutrients (si are the elements of S ) consumed from the environment; r0 – genetically determined non→ specific nutrient utilization rate; C – vector of specific nutrient utilization rates (traits controlled by corresponding genes. ci are the ele→ ments of C ); P – population size (abundance, number of cells); G – number of metabolic systems of the population (characterizes the length of the corresponding genome); Iconsumed – set of indices of specific nutrients consumed by the cells of the population; kdeath – death coefficient; b – penalty coefficient; k1 – a parameter that determines the threshold value for increase of genome size penalty; n – a degree of nonlinearity of genome size penalty. Actually, n 0 and si are functions of P since the amount of consumed nutrients is proportional to the population size (see Supplementary materials for details). The first term refers to uptake of energy, the second term accounts for the negative density dependence and the last term describes a genome size/metabolic complexity penalty. The population growth rate is positively affected by an amount and number of consumed nutrients and negatively affected by both time spent on the replication (in our case it depends on the number of metabolic systems and, therefore, metabolic complexity can be considered as an indirect measure of the genome size) and time spent on utilizing the nutrients. We chose k1, n and b parameters (eq. 1) equal to 0.37, 10.0 and 0.01 respectively, which corresponds to a low genome size penalty according to the classification presented in (Lashin et al., 2012). In the case of motile cells, we regard 10% of the population as active migrants, which is in broad agreement with observed values for coastal seawater samples (Mitchell et al., 1995). We have carried out a number of simulations using the models of evolution and development of a microbial community inhabiting a spatially structured environment. We considered two types of 10-noded 1D-tubes: the first is the tube with non-permeable boundaries (Fig. 1c, 1D-simple), and the second is the tube with semi-permeable boundaries (Fig. 1c, 1D-spb). Both types of tubes have reasonable physical meaning. For instance, sewage pipes exemplify 1D-simple tubes, while

2.1.2. Model M2: stochastic model of evolutionary processes in a developing microbial community We define the stochastic model of evolutionary processes in a developing microbial community as model M2. Initially the community comprises a trophic ring of three symbiotic populations – all three initial populations of cells belong to a quasi-mutualist ecogroup. Each population produces a specific nutrient that is consumed by the following population in clock-wise direction (Fig. 1a). The initial conditions regarding to the population sizes and concentrations of nutrients were set the same as in the model M1. These model genes are meant to be essential for metabolic systems of the cells providing them with an ability either to utilize or to synthesize a specific nutrient. During the simulation course, the processes of horizontal gene transfer and/or gene loss occur with the probability equal to 10–7 per cell per generation. Thus, novel populations characterized by novel sets of metabolic functions emerge resulting in further evolution of the community (as it is previously described in Fig. 1b). 2.2. Ecological functional groups (ecogroups) For the convenience of further analysis, we have proposed the classification that takes into account, on the one hand, the role of organisms in their mutual interactions and, on the other hand, involves their function in the community structure estimating the complexity of their metabolism. Therefore, we classify all possible populations in the model with one nonspecific and three specific nutrients based on two features – number of specific metabolic systems (M, the range of metabolic activities) and skew (U = the number of metabolic systems of synthesis – the number of metabolic systems of utilization of specific nutrients) towards either utilization (negative semi-axis) or synthesis (positive semi-axis). Thereby, building our classification we have synthesized approaches based on pairwise interactions of species and those based on population 68

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

Fig. 1. The principal scheme of structures and processes described in the models of evolution and development of a microbial community inhabiting a spatially structured environment. a) Initial graph of the trophic relations in the community; b) Acquisition and loss of metabolic systems via horizontal transfer and loss of the genes that are responsible for particular metabolic traits; c) Two kinds of the spatial structure of the environment: the upper one is 1D-simple and the lower one is 1Dspb (semi-permeable boundaries). The nutrient supply, flow, diffusion and migration are also depicted (see the legend). The “beginning” is on the left side; d) Feasible set of metabolic systems combinations in the model with one nonspecific and three specific nutrients. Fig. 2. Ecological functional groups (ecogroups) classification scheme. X-axis represents the skew towards synthesis (positive semi-axis) or utilization (negative semi-axis) while y-axis represents total amount of specific metabolic systems. M]U – altruism line. M = -U – egoism line. U = 0 – balance line.

• Minimal genome – the only one population that is able to utilize only nonspecific nutrient and does not produce anything. • Producers are positioned along the altruism line – a class of organ-

function in the community organization. We have called the resulting classes “ecogroups” and given them names according to their ecological role performed by their representatives in a community (Fig. 2). The number of specific metabolic systems M can be used to distinguish populations of microorganisms possessing relatively small and large genomes and respective ranges of metabolic activities. Different criteria are to be applied for each of these groups considering their ecological functional status. Given M ≤ 3 (populations possessing small genomes):

• 69

isms generating metabolic impasses that can be utilized by other organisms as an energy source, structural component, electron donor or acceptor, etc. Quasi-commensals are positioned along the egoism line – a class of organisms consuming products of others. In contrast with classic commensals, quasi-commensals compete for a nonspecific nutrient

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

populations perform in the microbial community; 4 Ecogroup abundance dynamics analysis: associate transformations of ecological structure of a community with its population dynamics for both motile and nonmotile cells; 5 Principal component analysis (PCA) of the ecogroup abundance data: revealing clusters that correspond to certain ecological structures and estimating significance of the revealed clusters by Welchtest for means; 6 Analysis of population characteristics: species richness and total biomass.

with their breadwinners.

• Lying between egoism and altruism lines symmetrically with respect

to balance line, quasi-mutualists constitute a class of organisms possessing a small and relatively balanced genome (with respect to utilization/synthesis) and, therefore, potentially capable of making symbiotic trophic rings with other similar organisms. In contrast with classic mutualists, quasi-mutualists compete for a nonspecific nutrient with their symbionts. Given 4 ≤ M ≤ 6 (populations possessing large genomes):

• Given

•

|U| = 2, a population belongs either to quasi-commensal ecogroup or to producer ecogroup depending on the U sign (see Fig. 6). Since we consider in our model only those ecogroups that can be distinguished based on three kinds of specific nutrients, none of the populations can have more than three metabolic systems of the same type (either utilization or synthesis) that inevitably implies, in case of large genomes, a necessity to have metabolic systems of another type as well. Therefore, the populations possessing a large genome skewed significantly towards either utilization or synthesis also fall into these classes respectively. All other populations belong to aedificator class – the populations possessing large though balanced genome and characterized by a wide range of metabolic activities and, thereby, transforming the environment significantly thus affecting many other populations in the community.

3. Results 3.1. Genetic and ecological features of developing communities We have carried out a number of simulation runs using models of microbial communities’ evolution (see the model M2). The analysis of the dominant metabolic complexity revealed the following spatial localization patterns of species with various metabolic complexity: in chemotaxis-off case, metabolically rich populations with large genomes tend to dominate closer to the source of the nonspecific nutrient whereas the fraction of dominant populations with metabolically poor short genomes increases in line with the distance from the source. 3.1.1. Metabolic complexity level A mean metabolic complexity level calculated for every sub-habitat elucidates the issue of temporal stability of the previously revealed trend:

2.3. Analysis workflow We have carried out a number of simulations using the models of evolution and development of a microbial community inhabiting a spatially structured environment. There is a feasible set of metabolic systems combinations in the model that can be reached as new forms emerge due to horizontal gene transfer (HGT) and gene loss (Fig.1d). The populations with metabolically rich genome gather on one pole whereas metabolically poor and minimal genomes meet on another one. The populations of cells characterized by a particular combination of metabolic systems may be divided into five categories according to their functional role in the ecosystem (we refer to these groups as ecological functional groups or ecogroups further): aedificators, quasicommensals, quasi-mutualists, producers and minimal genome (for more details see the Methods section). The ecogroups have been used to analyze the processes of ecosystem development in the model. We regard a population in the community possessing the maximum size among all of the populations as a dominant population. We defined a threshold distinguishing the dominant population from non-dominant populations as a number of cells that discriminate between the sizes of the dominant population with maximum size and the non-dominant population with maximum size among all non-dominant populations (see Fig. S0 in Supplementary materials). So that, multiple dominants are possible. We have focused our analysis primarily on the case of the tubes with semi-permeable boundaries (1D-spb) since this case exhibit more diverse patterns than the case of the tube with non-permeable boundaries (1D-simple) does. However, we discuss the influence of the structure of the environment in detail in the 3.4 section. The results of each simulation run were analyzed as follows:

MCL =

∑i pi * gsi ∑i pi

where pi – the i-th population abundance, gsi – number of metabolic systems for the i-th population. In the case of nonmotile cells metabolic complexity level decreases as the distance from the source increases (see Fig. S1 in Supplementary materials). In the chemotaxis-on case, we have not found any patterns, which indicate the fact that the developmental processes in such systems can not be explained by trivial steady trends. 3.1.2. Locally dominant ecological groups The metabolic complexity level alone can not say too much about the evolutionary direction of the community without taking into account spatially heterogeneous micro-conditions affecting the adaptation. Along these lines, we have made use of analytical tools that reveal the ecological meaning of the processes manifesting in the community. We have analyzed a spatial distribution of dominant ecological functional groups. This analysis unveiled the mechanisms underlying the detected patterns of the metabolic complexity distribution. Analyzing the results of the natural selection acting in the model of the community with a full set of metabolic systems combinations (model M1) is very helpful before tackling the investigation of evolutionary trends in microorganisms towards metabolism augmentation and simplification. The initial composition of populations consists of populations representing all feasible set of metabolic systems combinations in equal proportions (see Model M1 section) and no new population emerge during the course of simulation, though numerous populations extinct. Our simulations have demonstrated that in the deterministic model of the «complete» community (M1 model) the natural selection acts differently upon motile and nonmotile organisms (Fig. 3a,c): in the chemotaxis-off case ecosystems develop the structure of aedificators and quasi-commensals, notably biomasses of both ecogroups are almost equal in the sub-habitats that are closer to the source of the nonspecific nutrient, while in the distant sub-habitats quasi-commensals tend to prevail. In the chemotaxis-on case after 3000 generations, only the

1 Local dominant metabolic complexity analysis: what locations of the environment are dominated by the populations with generalized and specialized metabolism; 2 Evaluation of the metabolic complexity level: which of the trends detected on the 1st step are true not only for the dominant populations but for the whole community as well; 3 The analysis of ecological function that locally dominant 70

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

Fig. 3. The results of the analysis of locally dominant ecological groups. Left side: averaged by time biomass fractions of ecogroups for chemotaxis-off (a) and chemotaxis-on (b) cases in the sub-habitats of the «complete» community (deterministic model M1). Right side: dominant ecogroup distribution for chemotaxis-off (c) and chemotaxis-on (d) cases in the stochastic model of evolutionary processes in a developing microbial community (M2 model). Color depicts corresponding dominant ecogroup in a particular sub-habitat. Bar size of a certain color corresponds to a number of simulation runs which demonstrate predominance of the ecogroup corresponding to the bar color. We regard an ecogroup as a dominant one in a sub-habitat if total amount of cells that belong to the ecogroup outnumbers any other ecogroup in the system. The description of ecogroups is presented in the «2. Methods» section.

exhibits one more stable distribution of ecogroup abundances (quasicommensals and quasi-mutualist symbiosis as a core of the community, see Fig. 4d). Notably, that mode is feasible not only by developing communities on their early phase but also by relatively established systems (e.g., after 10 000 or even 18 000 generations).

metabolically complete population survives. Moreover, the closer subhabitat to the source of the nonspecific nutrient, the faster all the populations except for the metabolically complete population extinct. Thus, in this case the natural selection favors the population that has self-sufficient metabolism. On the other hand, the investigation of stochastic model of evolutionary processes in a developing microbial community (M2 model) sheds the light on ecological structural patterns of evolving communities. A comparison of local dominant ecogroup distribution in the chemotaxis-on and -off cases (see Fig. 3b, d) shows distinct patterns. In the first case, dominant ecogroups of final communities are typically aedificators or quasi-commensals. In the second case, nonmotile cells compose spatially stratified communities where local dominant ecogroups are specifically distributed depending on their proximity to the nutrient source. Aedificators tend to dominate near the source while quasi-commensals dominate closer to sink. It is worth mentioning that the population dynamics grouped by ecogroups (see Supplementary files, “EcoGroupAbundanceDynamics” folder) demonstrate that there is a biomass parity of these two ecogroups (see Fig. 4a for a typical example) although the portions of space under control are unequal. Thereby, it has been shown that the nonspecific nutrient gradient can cause spatial stratification of dominant ecogroups in the populations of nonmotile cells in contrast with the well-mixed case (0D-model) (Lashin et al., 2012), where either aedificators dominate in the ecosystem or quasi-commensals do as it happens in the chemotaxis-on case. Thus, spatial heterogeneity facilitates local maintenance of antagonistic evolutionary scenarios resulting in amplification or reduction of microbial genomes. In the case of chemotactically active cells, in general, there are two kinds of outcomes observed: either the population with self-sufficient metabolism (aedificator, see Fig. 4b) dominates or specialized quasicommensals do due to subdominant quasi-mutualists who feed them. It is worth noting that unlike deterministic model of the «complete» community (M1 model) with a full set of combinations of metabolic systems where the selection favors aedificators, stochastic model of evolutionary processes in a developing microbial community (model M2) considering genetic variation processes such as HGT and gene loss

3.2. Investigating ecological structures of developed communities with the principal component analysis 3.2.1. Whole-environment analysis Despite the apparent fruitfulness of the aforesaid approach, the data on dominant ecogroups per se is insufficient to acquire comprehensive knowledge about the ecological structure of the community. Therefore, we have applied the principal component analysis (PCA) to our simulation results (model M2) using the information about total biomass of each ecogroup (ecogroup abundance) in particular for the whole environment comprised by all the nodes and for the local habitats (particular nodes). Thereby, we have obtained an opportunity to reveal the clusters corresponding to particular ecological structures. The simulation results of stochastic modelling have been used as observations for the PCA whereas we regard a number of cells falling to a specific ecogroup as a variable. This number was measured both at simulation finish and in average. Thus, if x is an observation corresponding to a simulation run then x = (aedificator total biomass, producer total biomass, quasi-commensal total biomass, quasi-mutualist total biomass, minimal genome total biomass). Resulting picture (based on the data obtained at the end of simulation runs) reflects the biomass distribution of ecogroups that is typical for climax communities while average biomass analysis (based on the data averaged by time) allows revealing those ecogroups whose participation is essential for the development of corresponding communities. The analysis has shown that the chemotaxis-on and -off cases comprise distinct clusters and do not fuse (Fig. 5a, c). It gives us the evidence that an ability to migrate adaptively by the chemotaxis brings about distinction between the ecological structures of developed communities. The significance of such differences has been approved by the Welch’s t-test for aedificator and quasi-mutualist biomasses (p71

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

Fig. 4. Illustrative dynamics of ecogroups abundances for chemotaxis-off (left side: a,c) and chemotaxis-on (right side: b,d) cases in the model M2. x-axis – number of generations, yaxis – ecogroup abundances; a – established parity of aedificators with quasi-commensals; b – the dominance of aedificators; c – a transient phase characterized by a gradual increase of quasi-commensal fraction in the community (ecological succession); d – the dominance of quasi-commensals. Quasi-mutualists act as subdominants.

corresponding variable with the first two principal components calculated for averaged data, the contribution of producers into the final structure of the community is absolutely negligible (the absolute value of corresponding PCA loading < 1e-4 – see the Supplementary materials, “PCA loadings” folder). The same reason lies behind the fact that minimal genome does not play a significant role in the explanation of observed diversity. The simulation results for the chemotaxis-on case fall into two clusters (see Fig. 5a) that correspond to two final ecological configurations: 1) Dominance of aedificators accompanied by negligible

value < 0.05) using both final and averaged observations data (the results of performed statistical tests are presented in the Supplementary files, see t_tests.R). The conclusion to be drawn states: on average, aedificator biomass is higher in ecosystems consisting of nonmotile organisms (incapable of chemotaxis) whereas quasi-mutualist biomass is higher in the case of motile organisms. For biomasses of producers, quasi-commensals and minimal genome bearers, the differences are insignificant (p-value > 0.05). Despite the fact that producers do play a certain role in the community’s development as one can see from the correlations between

Fig. 5. Principal component analysis applied to the data on total biomasses of ecogroups for all the simulation runs (on the left side) or to the local observations from the chemotaxis-off case runs (on the right side). The projection on the first two principal components is presented. On the left side, circles depict chemotaxis-off case runs, triangles – chemotaxis-on case runs. On the right side, different habitats where the biomasses of ecogroups were calculated are depicted using different symbols. Variable axes are shown with the arrows: aedificator (AED) – aedificator total biomass; producer (PRD) – producer total biomass; quasi.commensal (QC) – quasicommensal total biomass; quasi.mutualist (QM) – quasi-mutualist total biomass; minimal genome (MIN) – minimal genome total biomass. The top (а, b) figures represent the final biomasses distributed by ecogroups; the bottom (c, d) figures represent average biomasses distributed by ecogroups (i.e. averaged by time). 72

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

Fig. 6. Final species richness (top) and total biomass (bottom) distributed by the habitats for chemotaxis-off (black) and chemotaxis-on (gray) cases (mean values for a sample of simulation runs of the model M2, error bars depict standard errors of the means). Standard errors of the mean for the biomass values is generally too low to be displayed in the figure. See these values in Supplementary materials ("Species richness and biomass" spreadsheets).

through x-axis appears to be explained by the distance towards the nonspecific nutrient source (notably, PCA loadings of the variables show that each ecogroup except for the minimal genome tends to prevail near the source), which causes the development of biomass gradient, whereas y-axis is characterized by the greatest and opposite contribution of aedificator and quasi-commensal biomasses. The biomass ratio between these two ecogroups mostly explain the variation through this axis (see the Supplementary materials, “PCA loadings” folder). In the chemotaxis-on case, there are two major directions of variability that can be distinguished – in the total biomass of quasi-commensals together with quasi-mutualists and in the aedificator biomass. Remarkably, the former direction corresponds to 95.9% of explained variance and high biomass values are achieved only near the nutrient source.

presence of other ecogroups (see Fig. 4b); 2) Dominance of quasicommensals accompanied by quasi-mutualists under almost complete absence of aedificators in the system (see Fig. 4d). The only exception is one simulation run where quasi-mutualists dominate while aedificators act as subdominants, however, during 20 000 iterations the model did not reached the steady state expressing irregular dynamics. The chemotaxis-off cases are clustered less distinctly into three groups; each of them shares common structure: dominating aedificators and quasi-commensals that was mentioned above, and the differences concern only the quantitative balance between these two ecogroups. 3.2.2. Local-habitat analysis We performed local-habitat analysis in order to understand how spatially heterogeneous local ecological conditions of the biotope affect the structure of the ecosystem developing in this biotope (M2 model). In the previous subsection, we have presented an approach to analyze the ecological structure of the community on a whole-environment level. Here, we apply a similar approach at the level of local habitats to identify local patterns. For each individual simulation run an observation is represented by the data obtained from a local habitat rather than from the whole environment. Since we possess a priori knowledge about where the observations were made, we are able to analyze the spatial patterns of biomass distribution of ecogroups without prior putting our knowledge into the analysis. Considering the fact that the total biomass decreases in line with the distance to the nonspecific nutrient source (see the next section), the null hypothesis states that each ecogroup preserves this trend. However, the analysis has shown that it is not the case (see details in the Supplementary materials, “Linear regression models” folder) for motile cells – aedificator biomass significantly increases in line with the distance to the nonspecific nutrient source (p-value < 0.01 for both averaged and final observations). The trend also is not abided for the cells with the minimal genome in the chemotaxis-off case (p-value < 0.05, for averaged data only). The PCA results for local biomasses of the ecogroups in the chemotaxis-off case are presented in Fig. 5b, d. The variation observed

3.3. Adaptive migration impact on species richness and total biomass in a community As we have mentioned above adaptive migrations play a significant role in the evolution. Therefore, their impact on species richness and total biomass, which includes cells of all ecogroups put together, in a community becomes a reasonable object of scientific scrutiny. To estimate this impact, let’s define two values:

¯ 〉 = < species richness>chemoff < species richness> 〈SR chemon

¯ 〉 = < total biomass>chemoff < total biomass> 〈TB chemon where < species richness>chemoff – the species richness in the whole environment averaged by the number of chemotaxis-off simulation runs, < species richness>chemon – the species richness in the whole environment averaged by the number of chemotaxis-on simulation runs, < total biomass>chemoff – the total biomass averaged by the number of chemotaxis-off simulation runs, < total biomass>chemon – the total biomass averaged by the number of chemotaxis-on simulation runs. 73

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

¯ 〉 in a developed It turned out that mean number of species 〈SR ¯ 〉 is 2.02 fold community is 2.84 fold lower and mean total biomass 〈TB lower, when the cells are able to migrate towards habitats with better conditions (see Fig. 6 for the average results across simulations). We interpret this result this way (see Fig. S2 in Supplementary materials): high motility allows populations utilizing the nutrient that is present in all parts of the environment with the maximum effect; it results in a rapid growth of populations utilizing available nutrients. As a result, well-established populations of microorganisms, which exhibit maximum resource utilization efficacy, prevent new populations from ecological niche realization though the latter carry a potential for transforming the environment to create new ecological licenses (Osche, 1966; Schmitt, 1987;Levchenko and Starobogatov, 1990). Thus, we have demonstrated in this simulation study that active migrations act as an impediment to the introduction of new forms possessing novel combinations of metabolic systems and as a factor decreasing total biomass of the community that is caused by a rapid expansion of wellestablished populations.

biomass of the community in comparison to the communities of nonmotile organisms. We explain this phenomenon with a concept of chain reduction of available ecological licenses caused by the decrease in diversity (see Fig. S2 in Supplementary materials). Despite the processes’ of genetic variation being inherently stochastic, there are robust self-organizing patterns on ecological structure level. Along with that, there are multiple basins of attraction for evolution of such systems. Thus, the space of evolutionary trajectories is channeled. Independently of the origin and exact metabolic sets of the populations, particular configurations of ecological functional types develop, and only in terms of these configurations we achieve plausible interpretations for the logic underlying the evolutionary processes in microbial communities. It alludes us to the ideas of George Zavarzin who pointed out that «in search for knowledge on the history of biosphere evolution one has to rely on the functional traits rather than on phylogeny of contemporary bacteria, i.e. on ecotype rather than ribotype» (Zavarzin, 2006). The importance of migration for evolutionary success of horizontal gene transfer events was shown in the work of Niehus et al. (Niehus et al., 2015). The migration provides a permanent gene flow of sufficient density that allows mutants to gain a critical biomass, which is necessary for the fixation in the community. The fixation of a population which obtains biomass exceeding a critical value, in its turn, enables such a rare event as horizontal gene transfer to determine the evolution of a microbial community overcoming the constraints imposed by the harsh competitive milieu of the established community. Thereby, the existence of such critical thresholds of biomass can lead to bi- or even to multistability of the systems of motile populations at the level of the distribution of ecological groups. Besides things described above, the energy cost of motility is a factor that might influence an outcome of evolution processes in a community. Though both theoretical estimates (Purcell, 1977; Mitchell, 1991) and the evidence (Stocker and Seymour, 2012) indicate that for large and slowly swimming bacteria chemotaxis constitute very small fraction (˜0.1%) of their energy budgets. However, bacteria living in the ocean tend to exhibit much higher velocities and energy cost of swimming is much higher for them comprising up to 50% of their total energy expenditure (Stocker and Seymour, 2012). Therefore, the investigation of its role is very promising and will certainly shed the light on the evolution in the communities of marine bacteria in the future. Using the bottom-up approach at the level of microbial communities allows us to elucidate the system properties underlying self-organizing patterns that we observe both in our model and in other theoretical and experimental studies (Matsuyama and Matsushita, 1993; Tsyganov et al., 1999; Bihary et al., 2012; Ghoul and Mitri, 2016; Goldschmidt et al., 2017). These patterns are of a great interest for solving the problems of ecological engineering (Mitsch and Jørgensen, 2003) since that knowledge allows to predict the structure of a climax community based on its initial composition. Taking into account genetic variation processes inducing changes in ecological structure of a microbial community proves to be extremely important for unravelling the principles of ecosystem’s self-design. Thus, in this simulation study we have shown that given even onedimensional nutrient gradient is presented in the environment, feasible outcomes in the terms of dominance of particular ecogroups are distributed throughout habitats. We have found that local dominant ecogroups distribution distinguishes in the chemotaxis-on and -off cases. Motile cells with a few exceptions exhibit the dominance of either aedificators or quasi-commensals throughout the whole environment, whereas nonmotile cells exhibit spatial stratification of dominant ecogroups in line with the distance towards nonspecific nutrient source. Aedificators dominate closer to the source while quasi-commensals do farther from it. It is worth mentioning that the analysis of abundance dynamics of the ecogroups show that there is a biomass parity trend among these two groups despite unequal portions of space under control. Generally, aedificator biomass is higher in the ecosystems that

3.4. The influence of the structure of the environment Investigating the influence of the structure of the environment in the models of evolution and development of a microbial community we have arrived to a number of conclusions. As we stated above different types of tubes correspond to different physical structures of the environment and can be considered as an abiotic factor influencing the formation of a community. We have analyzed the patterns of dynamics of ecogroups abundances to establish the differences distinguishing both cases. Along these lines, in the deterministic model of the «complete» community (M1 model) switching the environmental structure from the tube with semi-permeable boundaries (1D-spb) to the tube with nonpermeable boundaries (1D-simple) affects the patterns significantly only for the communities of nonmotile organisms (see Fig. S3 in Supplementary materials). It results in adjusting the ratio between the abundances of two dominant ecogroups (quasi-commensals and aedificators). As one can see from the dynamics, the effect of eliminating additional effluxes targets specifically quasi-commensals increasing the total biomass of this ecogroup. It can be explained by the fact that nonpermeable tube favors the accumulation of specific nutrients boosting the growth of quasi-commensals in the distant nodes. In the stochastic model of evolutionary processes in a developing microbial community (M2 model), switching from 1D-spb to 1D-simple affects both communities of motile and nonmotile organisms. In the chemotaxis-off case, the pattern of aedificator domination is lost for the same reasons as in M1 model since additional effluxes are no longer present (see Fig. S4 in Supplementary materials). In the chemotaxis-on case, the patterns exhibiting oscillating or switching dynamics of ecogroups abundances no longer observed (see Fig. S5 in Supplementary materials) in 1D-simple environments in contrast with the tubes with semi-permeable boundaries (1D-spb) (see Fig. S6 in Supplementary materials). 3.5. Discussion The computer modelling we have carried out shows that the spatial organization of the environment shapes the landscape of selection pressure, mediating thereby the evolutionary processes in the ecosystem. It leads to a stratification on ecosystem level resulting in spatially specific realization of various evolutionary scenarios that is in accordance with the results obtained in other investigations both theoretical (Goldschmidt et al., 2017) and experimental (Callahan et al., 2014; Goldschmidt et al., 2017) ones. Migrations as a process impeding expression of this stratification cause the decrease in community diversity and expansion of several fitted populations. However, it does not produce an increase in total 74

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

consist of nonmotile organisms (incapable of chemotaxis) while quasimutualist biomass is higher for ecosystems of motile organisms. Thus, different evolutionary trends sustain in habitats with contrasting ecological conditions due to the nutrient gradients that structure the environment spatially. The development of microbial community under natural selection in the deterministic model of the «complete» community (M1 model), given start sizes of all populations are equal, results in the dominance of the aedificator with a complete set of metabolic systems. However, the processes of horizontal gene transfer and gene loss in communities of motile microorganisms cause a novel attractor to emerge toward which the system tends to evolve in terms of biomass distribution of the ecogroups: not only the developing communities on their early phase attract to that mode but also the established communities do. Moreover, active migrations act as a factor impeding introduction of new forms and decreasing total biomass of the community that is caused by a rapid expansion of well-established populations resulting in the chain reduction of available ecological licenses. The structure of the environment is also a significant abiotic factor influencing the formation of a community. Our modelling shows that semi-permeability of environmental boundaries while providing additional effluxes for the substances and cells permits more diverse ecological patterns than those reached in the tubes with non-permeable boundaries. Along these lines, the context of interactions among populations and spatially heterogeneous ecological conditions affecting the adaptation determine the selection pressure at the ecosystem level favoring those trophic configurations that appear to be the most robust in these particular conditions.

Faust, K., Raes, J., 2012. Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10, 538–550. Fullmer, M.S., Soucy, S.M., Gogarten, J.P., 2015. The pan-genome as a shared genomic resource: mutual cheating, cooperation and the black queen hypothesis. Front. Microbiol. 6. Ghoul, M., Mitri, S., 2016. The ecology and evolution of microbial competition. Trends Microbiol. 24, 833–845. Gogarten, J.P., Townsend, J.P., 2005. Horizontal gene transfer, genome innovation and evolution. Nat. Rev. Microbiol. 3, 679–687. Goldschmidt, F., Regoes, R.R., Johnson, D.R., 2017. Successive range expansion promotes diversity and accelerates evolution in spatially structured microbial populations. ISME J. 11, 2112–2123. Grime, J.P., 1974. Vegetation classification by reference to strategies. Nature 250, 26–31. Grimm, V., Berger, U., 2016. Structural realism, emergence, and predictions in nextgeneration ecological modelling: synthesis from a special issue. Ecol. Modell. 326, 177–187. Grimm, V., Ayllón, D., Railsback, S.F., 2017. Next-generation individual-based models integrate biodiversity and ecosystems: yes we can, and yes we must. Ecosystems 20, 229–236. Gruber, M.F., Johnson, C.J., Tang, C.Y., Jensen, M.H., Yde, L., Hélix-Nielsen, C., 2011. Computational fluid dynamics simulations of flow and concentration polarization in forward osmosis membrane systems. J. Memb. Sci. 379, 488–495. Händel, N., Hoeksema, M., Freijo Mata, M., Brul, S., Ter Kuile, B.H., 2016. Effects of stress, ROS and the SOS response on de novo acquisition of antibiotic resistance in Escherichia coli. Antimicrob. Agents Chemother. 60, 1319–1327. Ho, A., Kerckhof, F.M., Luke, C., Reim, A., Krause, S., Boon, N., Bodelier, P.L.E., 2013. Conceptualizing functional traits and ecological characteristics of methane-oxidizing bacteria as life strategies. Environ. Microbiol. Rep. 5, 335–345. Klimenko, A.I., Matushkin, Y.G., Kolchanov, N.A., Lashin, S.A., 2015. Modeling evolution of spatially distributed bacterial communities : a simulation with the haploid evolutionary constructor. BMC Evol. Biol. 15, S3. Klimenko, A.I., Matushkin, Y.G., Kolchanov, N.A., Lashin, S.A., 2016. Bacteriophages affect evolution of bacterial communities in spatially distributed habitats: a simulation study. BMC Microbiol. 16, S10. Krause, S., Le Roux, X., Niklaus, Pa., van Bodegom, P.M., Lennon, T.J.T., Bertilsson, S., et al., 2014. Trait-based approaches for understanding microbial biodiversity and ecosystem functioning. Front. Microbiol. 5, 1–10. Kreft, J.-U., Booth, G., Wimpenny, J.W.T., 1998. BacSim, a simulator for individual-based modelling of bacterial colony growth. Microbiology 144, 3275–3287. Kunin, V., Ouzounis, Ca, 2003. The balance of driving forces during genome evolution in prokaryotes the balance of driving forces during genome evolution in prokaryotes. Genome Res. 13, 1589–1594. Langridge, G.C., Fookes, M., Connor, T.R., Feltwell, T., Feasey, N., Parsons, B.N., et al., 2015. Patterns of genome evolution that have accompanied host adaptation in Salmonella. Proc. Natl. Acad. Sci. 112, 863–868. Lashin, S.A., Suslov, V.V., Kolchanov, N.A., Matushkin, Y.G., 2007. Simulation of coevolution in community by using the “Evolutionary Constructor” program. In Silico Biol. (Gedrukt) 7, 261–275. Lashin, S.A., Matushkin, Y.G., Suslov, V.V., Kolchanov, N.A., 2012. Computer modeling of genome complexity variation trends in prokaryotic communities under varying habitat conditions. Ecol. Modell. 224, 124–129. Levchenko, V.F., Starobogatov, Y.I., 1990. Succession changes and evolution of ecosystem. Zh. Obshch. Biol. 51, 619–631. Lidicker, W.Z., 1979. A clarification of interactions in ecological systems. Bioscience 29, 475–477. Liu, W., Røder, H.L., Madsen, J.S., Bjarnsholt, T., Sørensen, S.J., Burmølle, M., 2016. Interspecific bacterial interactions are reflected in multispecies biofilm spatial organization. Front. Microbiol. 7. Liu, C., Bridges, M.E., Kaundun, S.S., Glasgow, L., Owen, M.D.K., Neve, P., 2017. A generalised individual-based algorithm for modelling the evolution of quantitative herbicide resistance in arable weed populations. Pest Manage. Sci. 73, 462–474. Margulis, L., 1970. Origin of Eukaryotic Cells: Evidence and Research Implications for a Theory of the Origin and Evolution of Microbial, Plant, and Animal Cells on the Precambrian Earth. Yale University Press, New Haven. Maslov, S., Krishna, S., Pang, T.Y., Sneppen, K., 2009. Toolbox model of evolution of prokaryotic metabolic networks and their regulation. Proc. Natl. Acad. Sci. 106, 9743–9748. Matsuyama, T., Matsushita, M., 1993. Fractal morphogenesis by a bacterial cell population. Crit. Rev. Microbiol. 19, 117–135. McCutcheon, J.P., Moran, Na., 2011. Extreme genome reduction in symbiotic bacteria. Nat. Rev. Microbiol. 10, 13–26. McInerney, J.O., McNally, A., O’Connell, M.J., 2017. Why prokaryotes have pangenomes. Nat. Microbiol. 2, 17040. Mira, A., Ochman, H., Moran, N.A., 2001. Deletional bias and the evolution of bacterial genomes. Trends Genet. 17, 589–596. Mitchell, J.G., 1991. The influence of cell size on marine bacterial motility and energetics. Microb. Ecol. 22, 227–238. Mitchell, J.G., Pearson, L., Bonazinga, A., Dillon, S., Khouri, H., Paxinos, R., 1995. Long lag times and high velocities in the motility of natural assemblages of marine bacteria. Appl. Environ. Microbiol. 61, 877–882. Mitsch, W.J., Jørgensen, S.E., 2003. Ecological engineering: a field whose time has come. Ecol. Eng. 20, 363–377. Moran, N.A., 2002. Microbial minimalism: genome reduction in bacterial pathogens. Cell 108, 583–586. Morris, J.J., Lenski, R.E., Zinser, E.R., 2012. The black queen hypothesis : evolution of dependencies through adaptative gene loss. MBio 3, e00036–12.

Software availability The HEC software is available at http://evol-constructor.bionet.nsc. ru/?page_id=34&lang=enpage. Conflict of interest The authors declare no conflict of interest. Acknowledgement The study was supported by the Budget Project 0324-2019-0040. Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.ecolmodel.2019.02. 007. References Archibald, J.M., 2011. Origin of eukaryotic cells: 40 years on. Symbiosis 54, 69–86. Barraclough, T.G., Balbi, K.J., Ellis, R.J., 2012. Evolving concepts of bacterial species. Evol. Biol. 39, 148–157. Bihary, D., Kerenyi, A., Gelencser, Z., Netotea, S., Kertesz-Farkas, A., Venturi, V., Pongor, S., 2012. Simulation of communication and cooperation in multispecies bacterial communities with an agent based model. Scalable Comput. Pract. Exp. 13. Bridier, A., Piard, J.-C., Pandin, C., Labarthe, S., Dubois-Brissonnet, F., Briandet, R., 2017. Spatial organization plasticity as an adaptive driver of surface microbial communities. Front. Microbiol. 8. Callahan, B.J., Fukami, T., Fisher, D.S., 2014. Rapid evolution of adaptive niche construction in experimental microbial populations. Evolution (N. Y). 68, 3307–3316. Cavalier-Smith, T., 2013. Symbiogenesis: mechanisms, evolutionary consequences, and systematic implications. Annu. Rev. Ecol. Evol. Syst. 44, 145–172. Chew, S.C., Kundukad, B., Seviour, T., van der Maarel, J.R.C., Yang, L., Rice, S.A., et al., 2014. Dynamic remodeling of microbial biofilms by functionally distinct exopolysaccharides. MBio 5 e01536-14-e01536-14. Eigen, M., 1971. Selforganization of matter and the evolution of biological macromolecules. Naturwissenschaften 58, 465–523. Emonet, T., Macal, C.M., North, M.J., Wickersham, C.E., Cluzel, P., 2005. AgentCell: a digital single-cell assay for bacterial chemotaxis. Bioinformatics 21, 2714–2721.

75

Ecological Modelling 399 (2019) 66–76

A.I. Klimenko, et al.

Theory and Practice. CRC Press, Boca Raton. Rodrıguez-Martinez, J.M., Pascual, A., Rodriguez-Martinez, J.M., Pascual, A., 2006. Antimicrobial resistance in bacterial biofilms. Rev. Med. Microbiol. 17, 65–75. Scheffer, M., Baveco, J.M., DeAngelis, D.L., Rose, K.A., van Nes, E.H., 1995. Super-individuals a simple solution for modelling large populations on an individual basis. Ecol. Modell. 80, 161–170. Schmitt, M., 1987. Ecological niche’ sensu Guenther and “ecological licence” sensu Osche–two valuable but poorly appreciated explanatory concepts. Zool. Beiträge 31, 49–60. Soucy, S.M., Huang, J., Gogarten, J.P., 2015. Horizontal gene transfer: building the web of life. Nat. Rev. Genet. 16, 472–482. Stocker, R., Seymour, J.R., 2012. Ecology and physics of bacterial chemotaxis in the ocean. Microbiol. Mol. Biol. Rev. 76, 792–812. Sukachev, V.N., 1928. Rastitelnye Soobshchestva: Vvedenie V Fitosotsiologiyu (Plant Communities: Introduction to Phytosociology), 4th ed. “BOOK” Leningrad-Moscow, Leningrad. Tsyganov, M.A., Kresteva, I.B., Aslanidi, G.V., Aslanidi, K.B., Deev, A.A., Ivanitsky, G.R., 1999. The mechanism of fractal-like structure formation by bacterial populations. J. Biol. Phys. 25, 165–176. Weiße, A.Y., Oyarzún, D.A., Danos, V., Swain, P.S., 2015. Mechanistic links between cellular trade-offs, gene expression, and growth. Proc. Natl. Acad. Sci. 112, E1038–E1047. Wimpenny, J.W.T., Colasanti, R., 1997. A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models. FEMS Microbiol. Ecol. 22, 1–16. Yao, Z., Le Maguer, M., 1996. Mathematical modelling and simulation of mass transfer in osmotic dehydration processes. Part I: conceptual and mathematical models. J. Food Eng. 29, 349–360. Zavarzin, G.A., 2006. Does evolution make the essence of biology? Her. Russ. Acad. Sci. 76, 292–302. Zavarzin, G.A., Kolotilova, N.N., 2001. Introduction to Environmental Microbiology. University Book House, Moscow.

Nelson, K.E., Clayton, R.A., Gill, S.R., Gwinn, M.L., Dodson, R.J., Haft, D.H., et al., 1999. Evidence for lateral gene transfer between Archaea and Bacteria from genome sequence of Thermotoga maritima. Nature 399, 323–329. Nesbo, C.L., L’Haridon, S., Stetter, K.O., Doolittle, W.F., 2001. Phylogenetic analyses of two “archaeal” genes in thermotoga maritima reveal multiple transfers between archaea and bacteria. Mol. Biol. Evol. 18, 362–375. Niehus, R., Mitri, S., Fletcher, A.G., Foster, K.R., 2015. Migration and horizontal gene transfer divide microbial genomes into multiple niches. Nat. Commun. 6, 8924. Osche, G., 1966. Grundzüge der allgemeinen phylogenetik. In: Gessner, F. (Ed.), Handbuch Der Biologie III. Athenaion, pp. 817–906 Frankfurt/M. Pang, T.Y., Maslov, S., 2011. A toolbox model of evolution of metabolic pathways on networks of arbitrary topology. PLoS Comput. Biol. 7, e1001137. Picioreanu, C., Kreft, J., Van, M.C.M., Van Loosdrecht, M.C.M., 2004. Particle-based multidimensional multispecies biofilm model particle-based multidimensional multispecies. Biofilm Model. 70. Polz, M.F., Alm, E.J., Hanage, W.P., 2013. Horizontal gene transfer and the evolution of bacterial and archaeal population structure. Trends Genet. 29, 170–175. Pomeroy, L.R., Hargrove, E.C., Alberts, J.J., 1988. In: Pomeroy, L.R., Alberts, J.J. (Eds.), Concepts of Ecosystem Ecology: a Comparative View. Springer-Verlag, New York, New York. Poplavskaya, G.I., 1924. Opyt fitosotsiologicheskogo analiza rastitelnosti tselinnoy zapovednoy stepi Askaniya-Nova (an essay on phytosociological analysis of vegetation of the virgin protected steppe of Askania-Nova). Zhurn. Rus. Bot. 9, 125–146 (in Russian). Provorov, N.A., Tikhonovich, I.A., Vorobyov, N.I., 2016. Symbiogenesis as a model for reconstructing the early stages of genome evolution. Russ. J. Genet. 52, 117–124. Purcell, E.M., 1977. Life At low reynolds number. Pdf. Am. J. Phys. 45, 3–11. Ramenskyi, L.G., 1935. O principialnyh ustanovkah, osnovnyh ponyatiyah i terminah proizvodstvennoi tipologii zemel, geobotaniki i ekologii (On the principal approaches, main concepts and terms of the economic typology of land). Sov. Bot. 4, 25–42 (In Russian.). Ravikrishnan, A., Raman, K., 2018. Systems-Level Modelling of Microbial Communities:

76