Horizontal gene transfer increases microbiome stability

Genes encoding resistance to stressors, such as antibiotics or environmental pollutants, are widespread across microbiomes, often encoded on mobile genetic elements. Yet despite their prevalence, the impact of resistance genes and their mobility upon the dynamics of microbial communities remains largely unknown. Here we develop eco-evolutionary theory to explore how resistance genes alter the stability of diverse microbiomes. We show that adding resistance genes to a microbiome typically increases its overall stability, particularly for mobile resistance genes with high transfer rates that efficiently spread resistance throughout the community. However, the impact of resistance genes upon the stability of individual taxa varies depending upon the mobility of the resistance gene and the network of ecological interactions within the community. Non-mobile resistance genes can benefit susceptible taxa in cooperative communities, yet damage those in competitive communities. Moreover, whilst the transfer of mobile resistance genes generally increases the stability of previously susceptible recipient taxa to perturbation, it can, counterintuitively, decrease the stability of the originally resistant donor species. We confirmed these theoretical predictions experimentally using competitive soil microcosm communities. Here the stability of a susceptible microbial community to perturbation was increased by adding mobile resistance genes encoded on conjugative plasmids but was decreased when these same genes were encoded on the chromosome. Together these findings highlight the importance of horizontal gene transfer in driving the eco-evolutionary dynamics of diverse microbiomes.

impact of resistance genes upon the stability of individual taxa varies depending upon the mobility 23 of the resistance gene and the network of ecological interactions within the community. Non-24 mobile resistance genes can benefit susceptible taxa in cooperative communities, yet damage 25 those in competitive communities. Moreover, whilst the transfer of mobile resistance genes 26 generally increases the stability of previously susceptible recipient taxa to perturbation, it can, 27 counterintuitively, decrease the stability of the originally resistant donor species. We confirmed 28 these theoretical predictions experimentally using competitive soil microcosm communities. Here 29 the stability of a susceptible microbial community to perturbation was increased by adding mobile 30 resistance genes encoded on conjugative plasmids but was decreased when these same genes 31 were encoded on the chromosome. Together these findings highlight the importance of horizontal 32 gene transfer in driving the eco-evolutionary dynamics of diverse microbiomes. For each community, we calculated ∆R and ∆E for the whole microbiome (Fig 1E), the focal 107 species only (i.e., the species that originally carried the resistance gene), and the background 108 community only (i.e., all species except the focal species). Using this modelling framework, we 109 could then explore in depth the impact of resistance genes on microbiome stability. 110 111

Mobile resistance genes increase stability of non-interacting communities 112
We began by exploring how the presence and mobility of resistance genes influenced the stability 113 of microbiomes in the absence of between-species ecological interactions (that is, all aij = 0). To 114 do so, we generated a set of communities with and without resistance genes. We then 115 systematically varied the ability of these resistance genes to transfer within and between species 116 (Fig 2A), capturing all degrees of gene mobility from immobile (e.g., a chromosomally encoded 117 resistance gene) to highly mobile (e.g., a resistance gene encoded by a highly conjugative and 118 promiscuous plasmid). 119 120 Introducing any resistance gene increased overall microbiome stability, decreasing the average 121 drop in community abundances following the onset of the perturbation (Fig 2B, ∆R > 0). However, 122 examining background and focal species separately revealed that in most cases this increase 123 was driven solely by the increased stability of the focal species. That is, the presence of a 124 resistance gene in the focal species drove up the average stability of the community as a whole, 125 but in most cases the stability of background species remained effectively unchanged ( Fig 2B). 126 To substantially increase the stability of the background microbiome, we found that the novel 127 resistance gene must be highly mobile. This was because, prior to the perturbation, the resistance 128 gene did not confer any benefit, and thus could only spread into the background community when 129 its transfer rate exceeded its rate of decline caused by negative selection against its cost. 130 However, low-mobility resistance genes could increase background species stability provided 131 additional forces enhanced the spread of resistance prior to any perturbation. For example, prior 132 exposure to weak stressor selection introduced a weak benefit to harboring the resistance gene 133 prior to the perturbation, enabling it to spread into the background community even at lower rates 134 of gene mobility ( Fig 2C). As a consequence, prior exposure substantially increased the stability 135 of background species, even in communities harboring low-mobility resistance genes ( Fig 2D). 136 137 Interspecies interactions modulate the impact of resistance genes on microbiome stability 138 Having established these baseline properties of the system, we next examined how the impact of 139 resistance genes upon community stability is modulated by interspecies interactions. Specifically, 140 we allowed individual species within our simulated communities to interact with one another in a 141 variety of different ways, ranging from competition and ammensalism (-/-and -/0 interactions 142 respectively), through exploitation (+/-), to cooperation and commensalism (+/+ and +/0). We then 143 generated a range of different microbial communities, systematically varying the proportion of 144 each interaction type ( Fig 3A). As previously, we then simulated the introduction of resistance 145 genes into these communities, also systematically changing the mobility of the resistance gene. 146

147
As in non-interacting communities, introducing any resistance gene typically increased overall 148 community stability, and this increase was higher for more mobile resistance genes ( Fig 3B, S1-149 5). However, this beneficial effect of resistance genes varied with interaction type and was far 150 stronger in microbiomes with a high proportion of cooperative interactions. In these cooperative 151 communities, individual species benefited both directly from acquiring resistance genes, and 152 indirectly from their cooperative partners acquiring resistance, which in turn helped to buffer the 153 negative impact of any perturbation. The principal effect of prior weak stressor exposure was once 154 again to reduce the level of gene mobility required for resistance genes to spread within the 155 community. And, as a consequence, this prior stressor exposure increased the stabilizing effect 156 of mobile resistance genes on overall community stability, across all interaction types (Fig 3E, S1-157 5, although, notably this effect was strongest for resistance genes with intermediate mobility). 158 While mobile resistance genes were beneficial for overall microbiome stability, examining focal 160 and background species separately revealed surprisingly different dynamics -with background 161 and focal species showing markedly different responses to resistance genes depending upon the 162 precise manner in which species were interacting and the mobility of the resistance gene ( Fig 3C,  163 D, S1-5). In cooperative microbiomes, background species benefited from the introduction of a 164 resistance gene regardless of its mobility. This occurred because, by promoting the survival of 165 species with whom a susceptible species cooperates, resistance genes aid recovery of the 166 susceptible species regardless of whether they have access to the resistance gene through HGT. 167 However, in competitive communities highly mobile resistance genes increased background 168 community stability while low-mobility genes reduced background community stability ( Fig 3C). 169 That is, in highly competitive communities most species were less stable when another member 170 of the community harbored an immobile or low mobility resistance gene than when all species 171 were susceptible. What drove this phenomenon? In fully susceptible competitive communities, 172 during a perturbation every species experienced a reduction in their net growth rate, typically 173 resulting in a decrease in their overall abundance. However, as a consequence each species also 174 benefited from some competitive release -that is, the negative impact of competitors was 175 reduced as these competitors also decreased in abundance. In contrast, if one species acquired 176 an immobile or low-mobility resistance gene then this focal species remained at a high density 177 during the perturbation -and as such, susceptible species not only suffered from the stressor-178 mediated inhibition, but also from continued strong competitive inhibition by the focal species. 179 180 This impact of competitive release also modulated the impact of prior selection -again with very 181 different impacts on background and focal species. Background species benefited from prior 182 exposure to low-level stressors regardless of community context (Fig 3F), as this initial selection 183 promoted the spread of mobile resistance into the community. Moreover, this spread of mobile resistance also stabilized cooperative focal species, as these species benefited from their 185 cooperative partners acquiring resistance genes and thus remaining at high densities during 186 perturbations ( Fig 3G). However, in certain competitive communities prior selection could, 187 counterintuitively, slightly reduce the stability of the focal species (∆E<0, Fig 3G)  pQBR103 and pQBR57. Whereas the chromosomally encoded resistance is non-mobile, the 202 conjugative plasmids can transfer mercury resistance to other species [23][24][25] . Having assembled 203 these communities, we allowed them to equilibrate for ten serial transfers (40 days) either with or 204 without weak mercury selection (Fig 4A). After equilibrating, communities were perturbed with a 205 pulse of high concentration mercury, then propagated without mercury for two additional serial 206 transfers, mirroring the mode of perturbation used in our modelling framework. We determined 207 the composition of bacterial communities by 16S rDNA gene amplicon sequencing before and 208 after the perturbation, and estimated the abundances of the total community by colony counts. As 209 in our models, we then quantified community stability as the change in species abundances 210 between the pre-and post-perturbation samples. Pseudomonas genera 28-31 , suggesting this is likely to have been a competitive community. As 219 such, our model predicts that an introduced resistance gene would be likely to increase overall 220 microbiome stability, but that background species would only benefit from plasmid-encoded 221 resistance, and might suffer from the introduction of chromosomally encoded genes. Moreover, 222 our modelling also predicts that prior selection might, counterintuitively, slightly decrease the 223 stability of the focal species when it harbored plasmid-encoded mercury resistance, due to 224 transfer of mobile resistance limiting the benefits of competitive release. 225 226 As predicted, the presence of any resistance gene in the focal species SBW25 significantly 227 increased overall microbiome stability, strongly reducing the decline in abundance caused by the 228 mercury pulse (Fig 4C, effect  Moreover, they also underline the importance of considering exactly how properties such as 302 microbiome stability are quantified. Relying solely on coarse, whole-community metrics such as 303 overall community abundances or diversity may mask striking differences between individual 304 community members, and risks obscuring important eco-evolutionary dynamics. 305

306
To identify broad patterns in the impact of mobile resistance genes upon microbiomes here we 307 used relatively simple ecological models. The advantage of these simple models is that we can 308 analyze large numbers of microbiomes in high-throughput. However, as a consequence, there 309 are ecological and evolutionary features that we have not explicitly accounted for. For example, 310 previous studies have suggested that plasmid transfer is more likely between closer phylogenetic relatives, due to constraints on plasmid host-range or spatial structuring 35 . Indeed, in our 312 experiments plasmid-mediated transfer of Hg R from P. fluorescens SBW25 appeared to be limited 313 to congenerics -with more distantly related members of the community such as B. megaterium 314 apparently unable to gain the resistance genes, suggesting they were unable to acquire or 315 maintain either plasmid. Similarly, our model does not explicitly include de novo mutations or 316 phenotypic changes that may modulate how individual taxa respond to stressors without requiring 317 the acquisition of resistance genes from other community members. Indeed, in our experiments 318 we observed an increase in the stability of background species after prior exposure to mercury 319 within the chromosomally-encoded resistance treatment (Fig. 4D), suggesting mechanisms other 320 than plasmid transfer may have driven this increased stability 36 . Each of these phenomena has 321 the potential to further modulate the impact of horizontal gene transfer upon microbiome stability, 322 and disentangling exactly how will be critical if we wish to fully understand and ultimately predict 323 microbiome dynamics. 324 325 326

Materials and Methods 327
Underlying microbiome model 328 In line with previous work 18,19 , we model each microbiome as a network, where each node 329 represents a species and each edge captures the interaction between them. However, now we 330 extend this model to include two populations for each species, one with plasmids and one 331 without 25 . For each species i, plasmid-free cells grow at a rate ri and plasmid-bearing cells grow 332 at a rate ri -c, where c indicates the cost of plasmid carriage. Plasmids can transfer within and 333 between species, with the per cell rate of plasmid transfer from species j into species i defined as 334 γ !" = γ $ + ϵ !" , where γ $ is the average plasmid transfer rate, and ϵ !" is a noise term drawn from a 335 Normal distribution to introduce variability in plasmid transfer rates between species (note, in instances where γ !" < 0 we set γ !" = 0). Finally, we introduce an inhibition term −β , where D 337 defines the level of antibiotic or toxin in the environment, and β the susceptibility of the population 338 in question. More specifically, we set β # = 1 for the susceptible population, and β $ = 0.1 for the 339 resistant population. Together, this enables us to define the growth rate of the plasmid-negative, 340 # , and plasmid-positive, $ , populations of a given species as, 341 community. Notably, this model can also be extended to incorporate the phenomenon whereby 347 plasmids can be lost during segregation at a rate δ, however, this does not qualitatively alter the 348 results ( Fig S5). 349

350
Having established this new model, we generated a series of microbiomes, each composed of 351 N=10 species. Within any given microbiome, each species i interacts with species j with probability 352 C. To investigate how interspecies interactions modulate the effect of resistance genes, we 353 systematically alter the proportion of individual interaction terms, aij, that are positive, P + , such 354 that when P + = 0 the community is purely competitive, when P + = 0.5 the community contains 355 a mixture of all interaction types, from competitive and ammensal, through exploitation, to 356 cooperative and commensal, and when P + = 1 the community is purely cooperative. Finally, we 357 choose the magnitude of each aij is from a half-normal distribution with standard deviation σ = 358 0.05, and set the intrinsic growth rates rij such that the community has a linearly asymptotically stable equilibrium ! * ≈ 1∀ ∈ 1, … , in the absence of any resistance gene. Importantly, we find 360 similar results when varying each of these key microbiome parameters (Figs S3-5). 361 362 Quantifying microbiome stability 363 We calculate the stability of any given community by simulating its dynamics in response to a 364 perturbation. Specifically, we allow the community to equilibrate, then briefly perturb it by setting To test the accuracy of our predictions, we assembled experimental model microbiomes, 380 composed of a defined background community augmented with a predetermined focal species. 381 For our focal species we used Pseudomonas fluorescens SBW25 labelled with a gentamicin 382 resistance marker using the mini-Tn7 transposon system 25,36-38 . More specifically, we generated 383 independent P. fluorescens strains that were either susceptible to mercury, or harbored harboured while horizontal gene transfer enables resistance genes to spread within and between species 586 (black arrows). B. Schematic illustrating representative microbiome dynamics, capturing species 587 dynamics before, during, and after an external perturbation (lightning bolt). We calculate each 588 species's abundance immediately before, A b , and after, A a , the perturbation period, then define 589 stability as the average logged fold change for each species (mean log10(A a /A b )) C. Schematic 590 illustrating our four modelling scenarios: communities with and without resistance genes, and with 591 and without prior exposure to low-level stressors. D. Comparing the four scenarios allows us to 592 calculate the change in microbiome stability that results from the introduction of a resistance gene 593 (∆R), and the change in stability that results from prior exposure to weak low-level selection (∆E). 594

E.
To disentangle the impact of resistance and selection on different taxa we calculate ∆R and 595 ∆E for the total community, the background community only and the focal species alone. 596