Modeling the ecology of parasitic plasmids

Single plasmid, single-population models

To understand the dynamics of parasitic plasmids in complex ecologies, we first need to understand their behavior in simple scenarios. In this section, we analyze the dynamics of plasmids spreading by different HGT mechanisms in single populations. We begin by modeling competition between plasmid-free cells and cells containing a conjugative plasmid. A nutrient, with concentration (C), is supplied to the system at rate (S). Cells grow at a rate proportional to (C) with proportionality constant (alpha) for plasmid-free cells or ((1 ,-, {Delta})alpha) for plasmid-containing cells. Since we are interested in parasitic plasmids, we assume that ({Delta} in (0,1)). Cells of both types die at a rate (delta). When a plasmid-containing cell divides there is a loss probability, (p_ell), for one of the daughter cells to contain no plasmids. As long as a daughter cell contains at least one plasmid, the original plasmid copy number (the number of copies of the plasmid maintained per cell) is regenerated (as depicted in Fig. 1A). Plasmids can spread horizontally by conjugation, as illustrated in Fig. 1B, wherein a plasmid-free cell and a plasmid-containing cell interact to produce two plasmid-containing cells. We model the rate of conjugation by a mass-action term with rate (gamma _{mathrm{c}}). The equations governing the dynamics of conjugation are therefore:

In this model, what are the conditions for a parasitic conjugative plasmid to be able to invade a plasmid-free population? Invasibility implies that the equilibrium containing only plasmid-free cells is locally unstable, which occurs when

where (rho ^ ast ,=, S/delta) is the steady-state abundance of the plasmid-free cells at the plasmid-free equilibrium. This invasibility condition has an intuitive physical interpretation: to invade, the rate of conjugation must overcome losses due to reduced host growth rate as well as plasmid loss during division. This condition is similar to those found in previous studies [15].

Given the condition for plasmid invasion in Eq. 4, what is the optimal behavior for a parasitic conjugative plasmid? The left-hand-side of the expression is linear in the plasmid-free population, meaning that it is more difficult for a plasmid to invade smaller populations. To favor invasion, the plasmid can minimize the right-hand-side of the equation. For a plasmid that relies on random segregation upon cell division, both the plasmid cost ({Delta}) and the loss probability (p_ell) are functions of plasmid copy number, (n_{mathrm{p}}), a property controlled by the plasmid itself. If the primary cost of a plasmid is its replication and its gene products, plasmid cost will scale with copy number such that ({Delta} ,=, {Delta}_{mathrm{p}}n_{mathrm{p}}), where ({Delta}_{mathrm{p}}) is the cost of an individual plasmid copy. The loss probability will be (p_ell ,=, 2^{1 ,-, n_{mathrm{p}}}), i.e., the probability that a daughter cell receives zero plasmids from random segregation. The right-hand-side of the invasion condition Eq. 4 is therefore (delta ({Delta}_{mathrm{p}}n_{mathrm{p}} ,+, 2^{1 ,-, n_{mathrm{p}}}(1 ,-, {Delta}_{mathrm{p}}n_{mathrm{p}}))), which has a minimum at finite (n_{mathrm{p}}). The minimum in the invasion boundary at finite (n_{mathrm{p}}) indicates that in our framework optimal conjugative plasmids have a moderate copy number.

What kinds of ecological dynamics does our model for a conjugative parasitic plasmid exhibit? To answer this question, we characterize the stability of the system’s equilibria (see SI Appendix 1 for details). For conjugative plasmids with the optimal copy number, the dominant form of loss will be from reduced host fitness (see SI Fig. S1), and thus we characterize the case of negligible loss rate (p_ell ,=, 0) (we consider the case of finite loss rates in SI Fig. S2 and find similar results). In Fig. 1C we show the phase diagram of possible ecological outcomes as a function of plasmid cost ({Delta}) and conjugation rate (gamma _{mathrm{c}}). For high values of plasmid cost and low values of conjugation rate, the plasmid is unable to invade and the plasmid-free equilibrium is the only stable state. As plasmid cost decreases or conjugation rate increases, plasmids are able to invade and there is a state of stable coexistence between plasmid-free and plasmid-containing cells. The range of conjugation rates permitting coexistence is larger for costlier plasmids. Once the plasmid cost is sufficiently low or the conjugation rate is sufficiently high, the unique stable state consists only of plasmid-containing cells (note that for finite values of loss rate (p_ell), this plasmid-only state will contain a small fraction of plasmid-free cells due to plasmid loss).

Conjugation is the best studied mechanism of plasmid transmission, but plasmids can instead be transmitted by transformation, whereby plasmid-free cells are infected by free-floating plasmids, as illustrated in Fig. 1B. We therefore consider a model for plasmid-spread via transformation in which cell death results in release of free-floating plasmids which can then infect cells by mass action at rate (gamma _{mathrm{t}}). For every cell death, (n_{{mathrm{eff}}}) free-floating plasmids are released and these plasmids decay at a rate (delta _{mathrm{p}}). The dynamics of transformative plasmids are therefore:

What is the condition for transformative plasmid invasion? The plasmid-free equilibrium is unstable if

The left-hand-side of Eq. 9 is similar to the conjugative plasmid invasion condition, with the conjugation rate (gamma _{mathrm{c}}) replaced by the transformation rate (gamma _{mathrm{t}}). The numerator of the right-hand-side is also similar, with the cell death rate (delta) replaced with the plasmid decay rate (delta _{mathrm{p}}). The primary difference is in the denominator, which is the difference between the number of plasmids released on cell death, (n_{{mathrm{eff}}}), and the total replication deficit of plasmid-containing cells. If this denominator is negative, the inequality reverses and the plasmid-free equilibrium is always stable.

The invasion condition in Eq. 9 determines the optimal (n_{mathrm{p}}) of transformative plasmids: if each plasmid within a cell has a fixed probability of remaining viable after cell death, (p_{mathrm{v}}), then (n_{{mathrm{eff}}}) will scale linearly with (n_{mathrm{p}}) such that (n_{{mathrm{eff}}} ,=, p_{mathrm{v}}n_{mathrm{p}}). If the denominator of Eq. 9 is positive, the optimal copy number will be (n_{mathrm{p}} ,=, 1/{Delta}_{mathrm{p}}), the point at which the host’s growth rate is driven to zero and the plasmid relies entirely on horizontal transfer to survive. These results are substantially different than in the case of conjugation: instead of restricting itself to a limited portion of the host’s metabolic budget, a transformative parasite maximizes its spread by using as much of the host’s resources as possible. This is reminiscent of the behavior of phages—suggesting a possible evolutionary link between parasitic plasmids and phages.

As in the conjugation case, we now explore the ecological outcomes possible with transformative plasmids. We again consider the case of negligible loss rate (p_ell ,=, 0) and characterize the stability of the equilibria (see SI Appendix 1 for details). For (n_{{mathrm{eff}}} ,> , 1), the system has similar ecological outcomes to the conjugative case, with the system transitioning through no-plasmid, coexistence, and plasmid-only equilibria as ({Delta}) decreases and (gamma _{mathrm{t}}) increases. Interestingly, when (n_{{mathrm{eff}}} ,<, 1), there is now a regime in which both the plasmid-only and no-plasmid states are locally stable, leading to bistability (Fig. 1D). A community in this bistable regime is subject to dramatic composition shifts if subjected to sufficiently large perturbations.

Our simple models suggest that both conjugative and transformative parasitic plasmids can persist for sufficiently low fitness burdens and sufficiently high transfer rates. However, these different transfer mechanisms lead to significantly different optimal infection strategies and ecological outcomes. Conjugative plasmids in our framework have moderate copy optimal number, while transformative plasmids behave in a virus-like manner and produce as many copies as possible.

Multiplasmid, single-population models

Whether by invasion from other populations or mutation of existing plasmids, a population of cells is likely to encounter a wide variety of plasmids. To explore this scenario, we expand our model to allow for co-infection by multiple parasitic plasmids. We assume plasmids are from different incompatibility groups, meaning they do not interfere with each other’s replication, segregation, or other functions necessary for plasmid persistence and can coexist over long times [31]. We show a schematic of independently segregating plasmids in Fig. 1E. We also assume that plasmids transfer between cells independently. As shown in Fig. 1F, fitness costs are assumed to be multiplicative, such that the total cost for (m) plasmids is ({Delta}_{{mathrm{tot}}}(m) ,=, 1 ,-, (1 ,-, {Delta})^m). Note that assuming additive fitness costs does not qualitatively change the following results, and for small plasmid costs the two schemes are equivalent. For a system of (m) plasmids, this leads to (2^m ,+, 1) deterministic equations (we show an example model in SI Appendix 2). If the plasmids have identical properties and begin at equal abundances, this can be reduced to a system of (m ,+, 2) equations.

For the case of highly infectious plasmids, we can infer the behavior of the multiplasmid model from the single-plasmid dynamics. Consider a plasmid that is sufficiently infectious such that when it invades a population, every plasmid-free cell is infected. Now consider the case of another similarly infectious plasmid invading this population. Since every cell contains the first plasmid, the second plasmid has a similar (slightly smaller in the multiplicative case) effective fitness cost as it would have had invading a plasmid-free population. Thus, the second plasmid is also able to invade and infect all hosts. The second invasion deforms the fitness landscape further, meaning that a third plasmid can then invade. This process of repeated invasions will lead to plasmids reducing the host fitness to arbitrarily low levels. These runaway invasions can be viewed as a plasmid “tragedy of the commons”. In the classical tragedy of the commons, individuals acting in their own self-interest deplete a shared resource until it is unusable for all individuals. In this case, the host cell is a common resource for plasmids. Individual plasmids only use part of the cell’s resources, but together they can completely deplete the resource, ultimately destroying the host.

What about the case of moderately infectious plasmids that do not infect every plasmid-free cell? When multiple moderately infectious plasmids invade a population, the result is a distribution of cells with varying numbers of unique plasmid types per cell. If a single plasmid stably infects a fraction (p) of the population, then if another identical plasmid invades and is completely independent of this first plasmid, it will also infect a fraction (p) of the population. Thus, if plasmids are identical and independent, the plasmids will be distributed among cells as ({mathrm{Binomial}}(m,p)) where (m) is the number of unique plasmid types. If plasmids are independent but not identical, they will be distributed as a Poisson binomial. While real plasmids are not independent of each other (they are coupled by their fitness costs to the host), the distribution of plasmids with low fitness cost will be well-approximated by a Poisson binomial. The discrepancy between the true distribution and a Poisson binomial will grow as fitness cost increases (see SI Fig. S3). In Fig. 1G, we show an example of this distribution for the case of eight identical moderately conjugative infectious plasmids. Below the critical conjugation rate for a single plasmid, (gamma _{mathrm{c}}^ ast), plasmids do not persist. As the conjugation rate increases past (gamma _{mathrm{c}}^ ast), the distribution of unique plasmid types per cell becomes increasingly skewed towards the total number of unique plasmid types. Since the distribution of plasmid types is binomial-like, as the total number of unique plasmids in the system increases, the mean number of unique plasmid types per cell increases as well (for the exact binomial this mean is (mp), where (m) is the number of unique plasmid types). Thus, even invasion by moderately infectious plasmids can lower host fitness to arbitrarily low levels. At the single-population level, there is no mechanism to stop this process, except perhaps extinction of the population, which may occur with finite population sizes (see SI Appendix for details). However, even if hosts with extremely high numbers of plasmids went extinct, this repeated invasion scenario would still imply that hosts with large numbers of unique plasmid types should be common in nature. In striking contrast to this scenario, prior analyses have found that most genomes contain only 0–3 plasmids [12]. What mechanisms might be limiting plasmid invasion in nature?

Muliplasmid, multipopulation models

In single populations, runaway plasmid invasion occurs because each invading plasmids further deforms the fitness landscape. However, the natural world does not consist of a single well-mixed population. In nature, different populations are separated by physical and genetic barriers that limit HGT. These barriers limit the ability of plasmids to spread, possibly stopping runaway plasmid invasion.

To explore the impact of HGT barriers, we use a Wright–Fisher framework to model a metapopulation of cells. We model a set of (N) populations living in (N) isolated demes. At every time period (“epoch”), a new resident of each deme is selected from a weighted distribution of the current (N) populations. The probability of being selected is proportional to the average fitness of the population, with a population containing (i) plasmids having fitness of (w_i). During every epoch, each population may also be invaded by a new plasmid type with probability (q). We assume the invading plasmid is compatible with the existing plasmids within the cell (we explore case of significant incompatibility interactions in the SI Appendix). We assume that adaptation of an invading plasmid to a new population occurs on a timescale much faster than the length of a epoch, such that successful adaptation is accounted for within the invasion probability. We illustrate three epochs of this process in Fig. 2A. Formally, each epoch can be described as (N) draws from the following multinomial distribution:

A population’s fitness is dependent on the number of unique plasmid types it contains. Thus far, we have considered a simple multiplicative model. However, it has been demonstrated that plasmid–plasmid interactions can modulate plasmid properties. For example, one study found that the presence of a plasmid can reduce the fitness cost of an invading plasmid [12]. To account for this epistasis between plasmids, we also consider fitness costs that increase sub-multiplicatively (positive epistasis) or super-multiplicatively (negative epistasis). We show examples of positive epistasis, negative epistasis, and no epistasis in Fig. 2B.

What is the distribution of unique plasmid types across populations in our model with HGT barriers? We derive the stationary distribution of this model for the three different epistasis functions in Fig. 2B and plot them in Fig. 2C (see SI Appendix 3 for details). For the case of no epistasis, the stationary distribution is Poisson-like. Positive epistasis favors carriage of multiple plasmids and results in an exponential-like distribution with a long tail. Negative epistasis has the opposite effect: it penalizes carriage of multiple plasmids and results in a sub-Poissonian distribution with a reduced tail. Importantly, in all cases the runaway invasion of plasmids is stopped. While there is nothing stopping individual populations from being overrun by invading plasmids, these populations are more likely to be out-competed by populations with fewer plasmids. Thus, the single-population “tragedy of the commons” is counteracted at a higher level of selection.

Analysis of natural genomes

How does our predicted distribution of unique plasmid types per cell compare to that in natural genomes? To make this comparison, we downloaded all complete bacterial genomes from NCBI (a total of 17,725 genomes) and analyzed their plasmid content. In Fig. 3A, we show the overall distribution of unique plasmid types per genome and corresponding model fits for both positive and no epistasis cases (see “Methods” for fitting details). The natural distribution is exponential-like and is well-fit by a model with positive epistasis. The model fit with no epistasis has too short a tail to be able to fit the data, and this problem becomes even more severe for negative epistasis. Thus, interestingly, we find that the distribution of unique plasmid types in real-world genomes is consistent with parasitic plasmids that ameliorate each other’s fitness costs. The degree of positive epistasis suggested by the data is quite strong—the distribution is nearly a pure exponential. In our model, this corresponds to the case in which the cost of all plasmids beyond the first is zero, such that for (m ,> , 1) the parameters controlling both population replication and plasmid invasion are independent of plasmid number. This means that the ratio between consecutive elements of the distribution is constant, yielding an exponential tail. In order to determine whether our conclusions are influenced by oversampling of clinically relevant species, we excluded 91 genera known to be clinically relevant or human-associated and repeated our analysis. The remaining dataset contains nearly 5000 genomes and still shows clear exponential behavior (see SI Fig. S4). We also analyzed whether the presence of engineered strains within the NCBI database influences our results. We found that there are only a small number of these engineered strains and that removing them had negligible impact on our results (see SI Fig. S5).

Can our model capture variation within smaller, related groups of genomes? In Fig. 3B we show the distribution of unique plasmid types per cell within the genus Escherichia. As can be seen, the data is very well fit by a model of parasitic plasmids with positive epistasis. However, our model was not able to capture some of the within-genus distributions we encountered. A notable exception is the distribution of unique plasmid types per cell in the genus Klebsiella, shown in Fig. 3C. In this genus, there is a substantial discontinuity between the zero-plasmid class and the rest of the distribution. While our simple Wright–Fisher model with some positive epistasis can capture the tail of the distribution, it then fails to capture the first few classes. Despite such exceptions, we find that the positive epistasis model is generally able to capture the overall trends in plasmid distributions over the bulk of natural genomes (see SI Fig. S6).

It should be noted that our current model of constant plasmid invasion probability and strong positive epistasis is not the only Wright–Fisher model that can produce an exponential distribution matching the data. We analyzed a more general form of the Wright–Fisher model in which the invasion probability and total fitness cost are arbitrary functions of unique plasmid number (see SI Appendix). We find that the general condition to yield an exponential is that the plasmid invasion probability and total fitness cost must be comparable regardless of the number of plasmids in the cell. These results indicate that even if there is no epistasis in fitness cost, an exponential can still result if there is positive epistasis in the invasion probability (i.e., if existing plasmids make it more likely for a new plasmid to successfully invade).

HGT barriers are not the only mechanism that can plausibly limit runaway plasmid invasion. Cells also have specialized systems to defend against foreign DNA, notably the CRISPR-Cas system [32]. To explore whether CRISPR-Cas is responsible for limiting plasmid invasion in natural genomes, we searched for cas genes within the NCBI complete bacterial genomes using HMMER (see “Methods” for details). We expect that if CRISPR-Cas plays a major role in limiting the spread of plasmids, the distribution of unique plasmid types per cell would be shifted towards lower plasmid numbers in genomes containing cas genes versus those lacking cas genes. In Fig. 3D, we show the distribution of unique plasmid types per genome in genomes containing at least one cas gene and those not containing any cas genes. The distributions are very similar, with no large differences between them. These results suggest that CRISPR-Cas is not a major mechanism limiting the spread of plasmids in bacteria. There are additional defense systems that may also influence plasmid carriage. However, a prior bioinformatics study found results similar to ours for restriction-modification (RM) systems, another defense system that protects against foreign DNA; the study examined the distribution of RM systems in bacterial genomes and found almost no relationship between the number of RM systems a genome encodes and the presence of plasmids (in one subset of data the authors actually found a positive relation) [33].

Source: Ecology - nature.com