Authors: Heather A Hopkins, Christian Lopezguerra, Meng-Jia Lau, Kasie Raymann
Categories: Article, bacteria, evolution, exaptation, opportunistic pathogen, predation, virulence, AcademicSubjects/SCI01130, AcademicSubjects/SCI01140
Source: Genome Biology and Evolution
Doi: 10.1093/gbe/evae149
Opportunistic pathogens are environmental microbes that are generally harmless and only occasionally cause disease. Unlike obligate pathogens, the growth and survival of opportunistic pathogens do not rely on host infection or transmission. Their versatile lifestyles make it challenging to decipher how and why virulence has evolved in opportunistic pathogens. The coincidental evolution hypothesis postulates that virulence results from exaptation or pleiotropy, i.e. traits evolved for adaptation to living in one environment that have a different function in another. In particular, adaptation to avoid or survive protist predation has been suggested to contribute to the evolution of bacterial virulence (the training ground hypothesis). Here, we used experimental evolution to determine how the selective pressure imposed by a protist predator impacts the virulence and fitness of a ubiquitous environmental opportunistic bacterial pathogen that has acquired multidrug Serratia marcescens. To this aim, we evolved S. marcescens in the presence or absence of generalist protist predator, Tetrahymena thermophila. After 60 d of evolution, we evaluated genotypic and phenotypic changes by comparing evolved S. marcescens with the ancestral strain. Whole-genome shotgun sequencing of the entire evolved populations and individual isolates revealed numerous cases of parallel evolution, many more than statistically expected by chance, in genes associated with virulence. Our phenotypic assays suggested that evolution in the presence of a predator maintained virulence, whereas evolution in the absence of a predator resulted in attenuated virulence. We also found a significant correlation between virulence, biofilm formation, growth, and grazing resistance. Overall, our results provide evidence that bacterial virulence and virulence-related traits are maintained by selective pressures imposed by protist predation.
Keywords: evolution, exaptation, virulence, predation, bacteria, opportunistic pathogen
Opportunistic pathogens are environmental bacteria that typically coexist peacefully alongside hosts and only sometimes cause disease, e.g. in animal hosts with immunodeficiency or microbiome imbalance (Brown et al. 2012). Opportunistic pathogens are intriguing because they have diverse lifestyles and experience strong tradeoffs between adaptation to life in the environment versus life inside of a host (Sokurenko et al. 2006; Brown et al. 2012; Pandey et al. 2022). Unlike obligate pathogens, the proliferation of opportunistic pathogens does not depend on their ability to infect a they are fully capable of reproducing and living in the environment (Casadevall and Pirofski 2007; Brown et al. 2012). Thus, the factors driving virulence in opportunistic pathogens remain largely enigmatic from an evolutionary standpoint. Several hypotheses have been put forth, proposing that virulence is selected for, or maintained by, interactions that occur outside a host. The coincidental evolution hypothesis (CEH) postulates that virulence evolves indirectly due to selection that occurs in nonhost environments (Levin and Edén 1990; Levin 1996), and protists have been proposed to act as “Trojan horses” (Barker and Brown 1994) or training grounds (Molmeret et al. 2005) for the evolution of bacterial pathogens. These hypotheses stem from the fact that bacterial defense mechanisms against environmental predators (e.g. protists) such as adherence, biofilm formation, grazing resistance, intracellular toxin production, and outer membrane structures also play a role in host infection (Brown et al. 2012; Erken et al. 2013; Amaro and Martín-González 2021). Additionally, several studies have shown that pathogens challenged with microbial predators exhibited increased host virulence (Cirillo et al. 1999; Rasmussen et al. 2005; Adiba et al. 2010; Coombes et al. 2011; Rehfuss et al. 2011; Hosseinidoust et al. 2013).
Experimental evolution can help reveal how different environmental pressures impact the evolution of bacteria and which genes are involved in specific phenotypic traits (Lenski 2017; Cooper 2018; McDonald 2019). In recent years, experimental evolution has been used to investigate if eukaryotic predation increases the virulence of opportunistic bacterial pathogens, but results have been conflicting (Friman et al. 2009; Mikonranta et al. 2012; Friman and Buckling 2014; Zhang et al. 2014; Nair et al. 2019; Hoque et al. 2022; Leong et al. 2022). One major limitation of many previous studies has been the lack of thorough phenotypic and genotypic characterization. As not all populations will follow the same evolutionary trajectory (Blount et al. 2008; Tenaillon et al. 2012; Rodríguez-Verdugo et al. 2014; Santos-Lopez et al. 2019, 2021), it is critical to evaluate and compare genotypes and multiple phenotypes in experimental evolution studies. For example, even if predation generally selects for the maintenance or enhancement of pathogen virulence, it is possible that some traits that are adaptive for defense against predators have no impact on virulence or even result in attenuation of pathogenicity within hosts. Nevertheless, experimental evolution can be used to answer many questions about fundamental evolutionary processes (Lenski 2017; Cooper 2018; McDonald 2019) and help reveal genes and mechanisms associated with virulence and other phenotypes. Thus, the goal of this study was to experimentally evolve the opportunistic bacterial pathogen Serratia marcescens in the presence or absence of an environmental protist predator and characterize both the genotypes and phenotypes of the evolved bacteria to elucidate how predation, or the lack thereof, impacts opportunistic pathogen evolution. We hypothesized that virulence would be increased or at least maintained in the presence of a predator.
Serratia marcescens is an opportunistic pathogen that ubiquitously occurs in soil and water (Grimont and Grimont 1978) and has been used for several bacteria–predator experimental evolution studies (Friman et al. 2009; Mikonranta et al. 2012; Zhang et al. 2014). It is a host generalist that can colonize many different organisms, where it can be pathogenic and/or commensal (Grimont and Grimont 1978). In humans, S. marcescens is a common enteric bacterium that is generally harmless in the gastrointestinal tract (Almuneef et al. 2001). However, it was recently found to be capable of damaging gut epithelial cells in vitro (Ochieng et al. 2014) and is responsible for causing dangerous nosocomial infections worldwide, notably in the respiratory and urinary tracts (Hejazi and Falkiner 1997; Khanna 2013). In the last few decades, the frequency of nosocomial infections caused by S. marcescens has increased, several hospital-wide outbreaks have been reported, and most strains have acquired multidrug resistance (Hejazi and Falkiner 1997; Khanna 2013; Montagnani et al. 2015). Additionally, S. marcescens has been recognized as an opportunistic pathogen of many plants and other animals (Grimont and Grimont 1978). Thus, understanding the evolution of virulence in S. marcescens is important for the health of a wide range of organisms. In this study, we used S. marcescens KZ19 (Raymann et al. 2018), which was isolated from the gut of a honey bee (Apis mellifera) and was pathogenic to honey bees following microbiome disruption (Raymann et al. 2017; Motta et al. 2018; Powell et al. 2021; Steele et al. 2021) or if administered orally to bees at high doses (Raymann et al. 2018). This strain of S. marcescens resides in a monophyletic clade with human and plant isolates, with the closest sequenced relative being a nosocomial strain BIDMC 50 (Raymann et al. 2018).
Here, we serially passaged S. marcescens KZ19 every day for 60 d in the presence or absence of a generalist protist predator, Tetrahymena thermophila. Following 60 d of evolution, we performed whole-genome shotgun (WGS) metagenomic sequencing on the entire S. marcescens–evolved populations, sequenced individual evolved isolates, and characterized the mutations present in both evolved isolates and populations. We observed multiple cases of parallel evolution which, when combined with phenotypic assays of isolates, enabled us to predict genes and pathways involved in virulence-related phenotypes. Although all populations did not follow the same evolutionary trajectories, in general, we found that virulence and grazing resistance were attenuated in the absence of a predator and predation resulted in increased biofilm production and slightly elevated virulence. Additionally, we found strong correlations between growth, predation (grazing) resistance, biofilm production, and pathogenicity (in honey bees), lending support to the hypothesis that predation resistance and virulence are often governed by the same mechanisms. In conclusion, our study provides evidence that experimental evolution can be used to help elucidate how exposure to nonhost environments impacts opportunistic pathogen evolution.
We confirmed that T. thermophila grazes on the ancestral bacterial strain by visualizing a 24 h coculture of T. thermophila and a transconjugant S. marcescens KZ19 containing an E2 crimson fluorescent protein (Leonard et al. 2018; Raymann et al. 2018) via fluorescent microscopy (Fig. 1). We then created six media-evolved (ME) and six predator-evolved (PE) replicate lines (populations). We minimized predator/prey coevolution by passaging evolved S. marcescens to new axenic T. thermophila cultures each day. The same concentration of S. marcescens was transferred for each passage, and we confirmed successful passaging via plating on Luria–Bertani (LB) agar every day. Viability and population size of T. thermophila was evaluated before every passage via light microscopy. Following 60 d of evolution, the entire PE (n = 6) and ME (n = 6) populations were sequenced using WGS metagenomic sequencing. The ancestor (n = 1) and individual isolates from the three ME (Lines 1 to 3) and six PE (Lines 1 to 6) populations were sequenced using WGS sequencing (supplementary table S1, Supplementary Material online). We then performed phenotypic assays on the evolved isolates (ME = 3 and PE = 6) in order to associate mutations with phenotypes.
Fig. 1. Confirmation of grazing activity. a) Phase contrast image of T. thermophila after 24 h coculture with fluorescently labeled S. marcescens KZ19 in Neff media with 180 µg/mL spectinomycin, b) Cy5 filter image showing the E2 crimson–labeled S. marcescens KZ19 cells inside T. thermophila, and c) composite overlay of the phase contrast and fluorescent images. Images were taken at 20× on a Keyence BZ-X700 series all-in-one fluorescence microscope and overlayed in ImageJ.
All evolved populations were compared with the ancestral genome to identify mutations. First, the ancestral genome was resequenced using both long (Nanopore)- and short (Illumina)-read sequencing technologies and assembled using SPades v3.15.3 (Antipov et al. 2016) to generate a consensus ancestral genome (supplementary table S1, Supplementary Material online). We then corrected for false positives (i.e. sequencing errors) by mapping all the short reads (244× coverage) of the ancestral strain back to the new consensus ancestral genome using breseq v0.38.2 (Deatherage and Barrick 2014). All mutations predicted via breseq when the ancestral reads were mapped to the ancestral consensus genome were considered false positives if detected in the evolved populations (i.e. they most likely correspond to sequencing or assembly issues in the ancestral genome).
We first analyzed the mutation patterns across each population. The total number of mutations identified in each evolved population (MEp n = 6; PEp n = 6) ranged from 6 to 106 (minimum frequency cutoff of 0.05) with an average of 44 mutations per population (supplementary Dataset S1, Supplementary Material online). The population with the fewest mutations identified was MEp3, and the two with the largest number of mutations were PEp3 and MEp2. However, the average genome coverage was negatively correlated to the total number of mutations (P = 0.0001, R^2^ = 0.7805; supplementary fig. S1a, Supplementary Material online) and the number of low frequency mutations (frequency < 10%, P = 0.0002, R^2^ = 0.7575; frequency < 20%, P = 0.0079, R^2^ = 0.5227; supplementary fig. S1b and c, Supplementary Material online) detected. A negative correlation between the number of mutations (total and low frequency) and genome coverage indicates that our sequencing depth was sufficient to characterize the mutations present in the populations.
In all evolved populations, we identified nonsynonymous, synonymous, and intergenic mutations (Fig. 2a; supplementary Dataset S1, Supplementary Material online). To determine the percentage of nonsynonymous, synonymous, intergenic, nonsense, and RNA mutations expected at random, we simulated independent mutations in the ancestral S. marcescens KZ19 genome by performing 10,000 independent simulations using a Jukes and Cantor model (Jukes and Cantor 1969). We found that the percentage of nonsynonymous mutations detected in our populations (49% to 68%) was between 4% and 23% less than the random expectation of 71.5% (Fig. 2a), indicating that overall, both environments imposed negative (purifying) selective pressure on S. marcescens (i.e. deleterious nonsynonymous mutations are being purged).
Fig. 2. Characterization of the mutations present in the evolved populations. a) Percentage of nonsynonymous, synonymous, intergenic, nonsense, and RNA mutations detected in each evolved population (ME1 to ME6 and PE1 to PE6) and compared with the percentage of each mutation type expected at random (expected). b) Number of observed parallel (gene- and site-specific) mutation events observed in ME, in PE, and in both ME and PE (ME–PE) populations compared with random expectation (expected vs. observed). Significantly more gene-specific and site-specific parallel mutations were observed than would be expected to occur by chance (χ^2^ test, P < 0.0001). The percentage of KEGG (Kanehisa and Goto 2000) functional categories of c) all genes in the ancestral genome, d) all genes that contained mutations in each group (ME, PE, or ME–PE), and e) the genes in which parallel events occurred across populations. The functional categories of all genes that contained mutations in ME, PE, and ME–PE populations d) did not significantly differ from the representation of functional categories in the ancestral genome (χ^2^ tests, P > 0.5). The functional categories of genes that contained parallel mutations in ME, PE, and ME–PE populations e) significantly differed from the representation of functional categories in the ancestral genome (χ^2^ tests, P < 0.0001).
In addition to purifying selection, we also found evidence that our populations were undergoing positive selection based on the identification of numerous mutations either in the exact same position (site-specific parallel evolution) or in the same gene (gene-specific parallel evolution) across populations (Fig. 2b; supplementary Dataset S1, Supplementary Material online). Parallel evolution events were considered if a mutation was found in the same position or gene in independently evolved populations (MEp, PEp, or both). None of the parallel evolution events occurred across all populations, confirming that they were not present in the ancestral population. We observed a total of 118 parallel mutations; 57 (9 same-site) were shared between the MEp and PEp populations, 29 (11 same-site) were unique to PEp populations, and 32 (6 same-site) were unique to MEp populations (Fig. 2b; supplementary Dataset S1, Supplementary Material online).
To determine the number of parallel mutations expected to occur by chance, we simulated independent mutations in the ancestral S. marcescens KZ19 genome by performing 10,000 independent simulations using a Jukes and Cantor model (Jukes and Cantor 1969). The number of site-specific parallel mutations expected to occur by chance within or across all our evolved populations was less than one, indicating that we found overwhelmingly more site-specific parallel mutations than random expectation (Fig. 2b; MEp P = 2.2e − 16, PEp P = 2.2e − 16, MEp–PEp P = 2.2e − 16, χ^2^ test). The number of gene-specific parallel mutations observed across the MEp and PEp populations was also four times more frequent than expected at random (Fig. 2b; MEp P = 8.7e − 25, PEp P = 6.2e − 05, MEp–PEp P = 6.4e − 10, χ^2^ test). The presence of excessive convergent parallel mutations provides strong evidence that positive selection occurred in all populations to adapt to their respective environments (the media and the presence of a predator).
We compared the functional categories of the mutations identified in only PEp, only MEp, or in both (MEp–PEp) lines to the overall representation of functional categories in the ancestral genome (Fig. 2c and d). The percentage of functional categories of all the genes in which mutations were detected in the PEp, MEp, and MEp–PEp lines did not significantly differ from the distribution of functional categories in the ancestral genome (P > 0.5, χ^2^ tests). However, when we considered only the genes in which parallel evolution events occurred (32 genes were impacted and 8 different promoter regions), we found that the types of genes impacted by parallel evolution events significantly differed from the distribution of functional categories in the ancestral genome (Fig. 2d; χ^2^ tests, P < 0.000). Although parallel evolution events arose across and within PEp and MEp populations, both lines were impacted in different types of genes. A large portion of the parallel mutations that occurred across only MEp populations were in genes involved in environmental information processing, while the parallel evolution events unique to PE populations were concentrated in genes involved in signaling and cellular processes and metabolism. Overall, parallel mutations were identified in genes involved in genetic information and processing (PEp = 17%, MEp = 0%, and PEp + MEp = 6%), signaling and cellular processes (PEp = 33%, MEp = 11%, and PEp + MEp = 6%), environmental information processing (PEp = 0%, MEp = 44%, and PEp + MEp = 12%), and metabolism (PEp = 33%, MEp = 11%, and PEp + MEp = 18%) and unclassified genes (PEp = 17%, MEp = 33%, and PEp + MEp = 59%). Parallel mutations present in only ME or shared between PE and ME populations are likely conferring adaptation to the Neff media, whereas parallel mutations only present in the PE populations suggest that they are important for fitness in the presence of a predator. Since most of the genes in which we identified mutations have not been functionally confirmed in S. marcescens, all gene names and predicted functions discussed hereafter must be considered candidate genes/functions.
We plotted all the mutations identified in each MEp and PEp population and their observed frequencies and found that ten of the parallel mutations were fixed within one or more of the populations in which they occurred (Fig. 3a). Interestingly, of these ten fixed parallel mutations, seven occurred in genes within the same protein–protein interaction network as predicted via STRING v12 (Szklarczyk et al. 2023) in close relatives of S. marcescens, such as Serratia rubidaea and Yersinia enterocolitica (Fig. 3b). This protein network contains the BarA-UvrY two-component regulatory system, and the three core proteins of the Rcs phosphorelay system (rcsC, rcsD, and rcsB). Both BarA-UvrY and the Rcs regulatory systems have been shown to play a major role in environmental adaptation and the regulation of virulence in S. marcescens (Liu et al. 2023; Romanowski et al. 2021; Di Venanzio et al. 2014) or related enterobacteria (Pernestig et al. 2003; Teplitski et al. 2003; Herren et al. 2006; Tomenius et al. 2006; Palaniyandi et al. 2012). Mutations in the histidine kinase barA were found in multiple PEp (3/6) and MEp populations (5/6), whereas mutations in the response regulator uvrY were mostly restricted to MEp populations (Fig. 3c). Mutations in the response regulator of the Rcs system, rcsB, were identified at high frequency in virtually all MEp populations (5/6) but never in the PEp populations (Fig. 3c). One MEp population (MEp6) contained a fixed mutation in the transmembrane hybrid kinase rcsC and one PEp population (PEp3) also had a mutation in rcsC but at very low (<10%) frequency (Fig. 3c). Within this network, fixed (100% frequency) site-specific parallel mutations unique to MEp populations (MEp2 to MEp6) were identified in a promoter region upstream of a gene annotated as the biofilm dispersion protein bdlA (Fig. 3c). Four of the six PEp populations (PEp2, PEp4, PEp5, and PEp6) contained mutations in a gene annotated as the diguanylate cyclase yfiN, and of these four populations, three (PEp4 to PEp6) also contained mutations in a promoter upstream of a gene annotated as the outer membrane protein ompC (Fig. 3c). The genes bdlA, yfiN, and ompC have all been implicated in biofilm formation in other Enterobacteriaceae species (Katharios-Lanwermeyer et al. 2022; Morgan et al. 2006; Sanchez-Torres et al. 2011; Petrova and Sauer 2012; Yeom et al. 2012; Huertas et al. 2014). It is important to note that the six PEp (Lines 1 to 6) and three MEp (Lines 1 to 3) populations were evolved at the same time, but three of the MEp (Lines 4 to 6) populations were started over 1 year later, demonstrating the reproducibility of our results.
Fig. 3. Frequency of mutations present in the evolved populations. a) Plot of all the detected mutations and their frequencies in each population highlighting the parallel “fixed” mutations. b) STRING v12 (Szklarczyk et al. 2023) predicted protein–protein interaction network containing seven of the ten genes in which we identified “fixed” parallel mutations. c) Heatmap showing the frequency of mutations in or adjacent (indicated by *) to the seven genes from the same protein–interaction network within each evolved population.
To attempt to link specific mutations with phenotypes, a total of three MEi and six PEi isolates were sequenced from the evolved populations and compared with the ancestral genome to identify mutations (Fig. 4a; supplementary table S2, Supplementary Material online). At least one mutation was detected in each isolate, with a maximum of four and an average of two mutations per isolate (Fig. 4a). We identified a total of 24 mutations within or in upstream promoters of 14 genes, and 70% of these mutations (17/24) were present at a frequency of >90% in their respective population (supplementary table S2, Supplementary Material online). Additionally, the isolates possessed all the fixed (100% frequency) mutations present within the population from which they were isolated (Fig. 3a; supplementary Dataset S1, Supplementary Material online).
Fig. 4. Genotypes and phenotypes of evolved isolates. a) Mutations present in each isolate grouped by nonsynonymous, synonymous, and promoter (intergenic mutation in a promoter region), indicating the mutation type (SNV, nonsense mutation, duplication, or deletion) and the gene (or gene upstream) affected. b) Mean growth rate of PE
iand MEiisolates and the ancestral strain based on 600 nm OD readings every hour for 24 h in Neff media. Each point represents the mean of three biological replicates, and the error bars indicate the standard deviation. c) Mean number of colony-forming unit per milliliter of PEiand MEiisolates and the ancestor after 18 h of growth on Neff agar. Each point represents a biological replicate. d) Mean biofilm formation of MEiand PEiisolates and the ancestral strain based on 590 nm OD readings after 24 h of growth in Neff media. Four assays were performed in triplicate, and each point represents a biological replicate e) Mean grazing (predation) resistance of PEiand MEiisolates and the ancestral strain based on 600 nm OD readings of bacterial filtrate after 48 h in coculture with T. thermophila. Two assays were performed in triplicate, and each point represents a biological replicate. f) Percent death of honey bees 5 d after exposure to the ancestral, MEi, and PEiisolates (control group average percent death = 2.3%; see supplementary fig. S4, Supplementary Material online). Three replicate assays were performed with 5replicates of 20 bees each per isolate per assay (15 replicates total per isolate). Each data point represents a biological replicate of 20 bees. g) Percent cytotoxicity of the ancestor, MEi, and PEiisolates to murine macrophages (RAW264.7) based on LDH release. b to g) Significance was tested by comparing each evolved isolate with the others and the ancestor using one-way ANOVA with Dunnett's multiple comparisons test. *P < 0.01; **P < 0.001; ***P < 0.0001. Dashed lines represent the overall mean when combining all MEior PEiisolates.
We identified three instances of site-specific parallel evolution and three cases of gene-specific parallel evolution across the isolates (Fig. 4a), all of which occurred within coding regions or promoters upstream of genes within the same BarA-UvrY Rcs protein–protein interaction network (Fig. 3b). Two of the three site-specific parallel mutations were unique to MEi a 66 bp deletion in rcsB (MEi1 and MEi3) and an intergenic single nucleotide variant (SNV) in a promoter upstream of a gene annotated as bdlA (MEi2 and MEi3; Fig. 4a). The third site-specific mutation was unique to PEi an intergenic SNV in a promoter upstream of ompC (PEi4 and PEi6; Fig. 4a). Two of the gene-specific parallel evolution events occurred in both ME and PE barA (MEi2, MEv3, PEv2, and PEv4) and uvrY (MEi1 and PEi5). The third case of gene-specific parallel evolution arose in a gene annotated as the inner membrane protein yfiN in all PEi isolates except PEi1 and PEi3 (Fig. 4a). Overall, most of the mutations present in the isolates likely represent adaptive mutations as they were fixed in their respective populations.
To evaluate the phenotypes of the PEi and MEi isolates, we performed growth, biofilm production, predation resistance (population size after growth in the presence of T. thermophila), and virulence (mortality of honey bees and death of murine macrophage) assays. We chose to specifically analyze isolates, rather than populations, in order to correlate fixed mutations with phenotypic characteristics. To compare growth between the ancestral and evolved strains, we monitored optical density (OD) in fresh Neff media (Fig. 4b), spent Neff media (filtered media that had supported T. thermophila growth for 48 h; supplementary fig. S2a, Supplementary Material online), and LB media (supplementary fig. S2b, Supplementary Material online) every hour for 24 h. None of the MEi and PEi isolates significantly differed in growth in either media types when compared with the ancestor (Fig. 4b; supplementary fig. S2, Supplementary Material online; P > 0.05, analysis of variance [ANOVA] with Dunnett's multiple comparison test). We also counted the number of colony-forming unit per milliliter for each isolate after 18 h of growth on Neff agar and found that two of the ME isolates (MEi1 and MEi3) produced significantly less colony-forming units than the ancestor (Fig. 4c; MEi1 P < 0.0001, MEi3 P = 0.005, ANOVA with Dunnett's multiple comparison test). However, none of the PEi isolates yielded significantly less colony-forming units than the ancestor (Fig. 4c), and no significant difference in the number of colony-forming unit per milliliter was found when the mean of MEi or PEi isolates was compared with each other or the ancestral strain (supplementary fig. S3a, Supplementary Material online).
Because biofilm formation is a strategy for resistance to predation and is also associated with virulence, we measured biofilm production in our MEi and PEi isolates compared with the ancestor (Fig. 4d). None of the MEi isolates significantly differed in biofilm production when compared with the ancestor, whereas five of the six PE isolates produced significantly more biofilms than the ancestor (Fig. 4d; PEi1 P = 0.017, PEi3 P = 0.003, PEi4 P = 0.0001, PEi5 P < 0.0001, PEi6 P = 0.035, one-way ANOVA with Dunnett's multiple comparison test). When considered together, mean biofilm production of the all PEi isolates was significantly higher than the mean of all the MEi isolates and the ancestor (supplementary fig. S3b, Supplementary Material online; P < 0.01, with one-way ANOVA with Dunnett's multiple comparison test).
In order to determine if exposure to a predator resulted in increased predation resistance, we evaluated grazing resistance (i.e. ability to proliferate in coculture with a predator) of the MEi and PEi isolates. Each isolate was grown with T. thermophila for 48 h, and then, the population density of the bacteria was measured. All MEi isolates were found to be significantly more susceptible to predation than the ancestor (Fig. 4e; P < 0.0001, one-way ANOVA with Dunnett's multiple comparison test). In contrast, none of the PEi isolates significantly differed in resistance to predation relative to the ancestor (Fig. 4e). However, the mean grazing resistance of all PEi isolates combined was significantly higher than the mean of the MEi isolates (supplementary fig. S3c, Supplementary Material online; P < 0.01, with one-way ANOVA with Dunnett's multiple comparison test).
We further determined how evolution under predation pressure impacts the virulence of S. marcescens KZ19 by performing survival assays in a natural host, the honey bee. Honey bees were orally exposed to each individual isolate, the ancestral strain, or sucrose solution only (controls) and monitored every day for 5 d (supplementary fig. S4, Supplementary Material online). We evaluated the percentage of death, 5 d postexposure, and found that two of the three ME isolates (MEi1 and MEi3), both which contain large 66 bp deletion in rcsB, displayed significantly attenuated virulence when compared with the ancestral strain (Fig. 4f; ME1 P = 0.001, ME3 P < 0.0001, one-way ANOVA with Dunnett's multiple comparison test). Conversely, none of the PEi isolates significantly differed from the ancestor in terms of virulence (Fig. 4f). Based on probability of survival using a Kaplan–Meier method, we obtained the same results, with only MEi1 and MEi3 displaying significantly decreased virulence when compared with the ancestor (supplementary fig. S4, Supplementary Material online; P < 0.05, Mantel–Cox Log-rank test with Bonferroni correction). When comparing the mean percent death 5 d postexposure by combining all MEi and PEi isolates, the MEi isolates displayed significant attenuation of virulence when compared with the ancestor and PEi isolates, and the PEi isolates were significantly more virulent than the ancestor and the MEi isolates (supplementary fig. S3d, Supplementary Material online; P < 0.01, with one-way ANOVA with Dunnett's multiple comparison test).
Several clinical S. marcescens isolates have been shown to be cytotoxic to macrophages (Ishii et al. 2012; Krzymińska et al. 2012, 2010). Although S. marcescens KZ19 was isolated from the gut of a honey bee, we found that it is capable of killing RAW 264.7 murine macrophages (Fig. 4g). Thus, we evaluated the percent of murine macrophage killing by the ancestor and compared it with the MEi and PEi isolates based on the release of lactate dehydrogenase (LDH) after 2 h of coculture [multiplicity of infection (MOI) 1]. None of the MEi isolates differed from the ancestor in terms of cytotoxicity to macrophages. Conversely, three of the PEi isolates displayed significantly less cytotoxicity to macrophage when compared with the ancestral strain (Fig. 4g; PEi1 P = 0.0009, PEi5 P = 0.03, PEi6 P < 0.0001, one-way ANOVA with Dunnett's multiple comparison test). When all MEi and PEi isolates were combined, no difference was observed between the MEi isolates and the ancestor or the PEi and MEi isolates, but the PEi isolates were significantly less cytotoxic to macrophage cells than the ancestor (supplementary fig. S3e, Supplementary Material online). However, macrophage killing appears to be time and density dependent. In a separate time series experiment with a different starting ratio (MOI 10), we found that all evolved isolates (MEi and PEi) displayed overall less cytotoxicity to macrophage than the ancestor at 2 to 5 h postexposure, but after 6 h, they converged in their level of cytotoxicity with all macrophages being killed (lysed) by all isolates and the ancestral strain (supplementary fig. S5, Supplementary Material online). Therefore, macrophage cytotoxicity may not be an informative metric for gauging virulence in S. marcescens KZ19.
Finally, we tested whether there was a direct link between growth, pathogenicity in honey bees, cytotoxicity to murine macrophage, biofilm production, and grazing resistance (Fig. 5; supplementary fig. S6, Supplementary Material online). Despite our relatively small sample size, we found a significant positive correlation between biofilm production and grazing resistance (Fig. 5a; R^2^ = 0.56, P = 0.01, Pearson correlation) as well as between virulence and biofilm production (Fig. 5b; R^2^ = 0.5, P = 0.02, Pearson correlation) and virulence and grazing resistance (Fig. 5c; R^2^ = 0.68, P = 0.003, Pearson correlation coefficient). We also found a significant positive correlation between growth, based on the mean number of colony-forming unit per milliliter, and virulence in honey bees (Fig. 5d; R^2^ = 0.63, P = 0.04, Pearson correlation). No significant correlations were observed between cytotoxicity to macrophage and any of the other tested phenotypes (supplementary fig. S6a to d, Supplementary Material online) or between growth (colony-forming unit per milliliter) and biofilm production or grazing resistance (supplementary fig. S6e and f, Supplementary Material online).
Fig. 5. Correlations between mean a) biofilm production and grazing resistance, b) biofilm production and virulence in bees, c) grazing resistance and virulence in bees, and d) growth (colony-forming unit per milliliter) and virulence in bees. Dashed lines represent a simple linear regression, and the correlation strength and significance were tested using the Pearson correlation coefficient.
Here, we revealed that only 60 d of experimental evolution in the presence or absence of a predator resulted in overall consistent genotypic changes in S. marcescens, indicated by numerous parallel evolution events. The lower-than-expected number of nonsynonymous mutations coupled with the higher-than-expected number of parallel evolution events suggests that all evolved populations underwent both positive and purifying selection. We also found evidence of positive selection based on the presence of fixed mutations (sweeps) within the same genes, and sometimes in the same exact site, across populations. This is noteworthy as bottlenecks are inherent to experimental evolution studies that rely on serial passaging and could result in the fixation of random mutations due to drift (Wahl et al. 2002; Mahrt et al. 2021; Gamblin et al. 2023). The presence of fixed mutations in the same genes and positions strongly indicates that drift was not the main driver of the sweeps observed in our populations. Although parallel evolution events and sweeps suggest that the genes impacted are adaptive, it is possible that some of the observed high frequency mutations are the result of hitchhiking.
Despite overall variation in the types of genes impacted, virtually all the fixed mutations within and across PEp and MEp populations occurred in genes that are part of the same predicted protein–protein interaction network. Within this protein–protein interaction network, several of the proteins are key components in gene regulatory systems. For example, BarA and UvrY make up a two-component regulatory system (Liu et al. 2023), and RcsB and RcsC are two of the three major components in the modified two-component Rcs phosphorelay system (Pan et al. 2021), both of which have been shown to regulate the expression of numerous genes in S. marcescens (Pan et al. 2021; Liu et al. 2023; Trouillon et al. 2023) and have also been implicated in virulence mechanisms in multiple species of Enterobacteriaceae (Heeb and Haas 2001; Pernestig et al. 2003; Teplitski et al. 2003; Herren et al. 2006; Huang et al. 2006; Tomenius et al. 2006; Palaniyandi et al. 2012; Meng et al. 2021). Specifically, the Rcs system has been shown to regulate the expression of genes involved in motility, capsule biosynthesis, biofilm formation, and virulence (Meng et al. 2021), and the BarA-UvrY system has been implicated in regulating the production of toxins, quorum sensing, motility, and several metabolic functions (Lapouge et al. 2008). Thus, mutations in these genes likely have major impacts on multiple genes by upregulating, downregulating, or stopping expression, resulting in strong selective pressure to maintain beneficial mutations and purge “deleterious” mutations. However, further studies are needed to validate this hypothesis.
In order to link specific mutations with phenotypic characteristics, we randomly chose a single isolate from each population. It is important to note that for the MEp populations, we only evaluated isolates from Lines 1 to 3 because when we originally started this experiment, we only evolved three ME populations. We later realized the importance of having more “control” populations and subsequently evolved three more MEp lines, which were only analyzed at the population level in this study. However, we would like to point out that the three MEp lines (MEp4 to MEp6) that were evolved over a year after MEp Lines 1 to 3 presented mutations in many of the same genes, including several of the fixed mutations that were observed in the first three MEp populations. This finding underlines the reproducibility of our results and eliminates the probability that contamination could explain the parallel evolution events observed across populations. The genotypic analysis of our three MEi isolates from Lines 1 to 3 and six PEi isolates (one from each evolved population) revealed that all isolates contained the fixed mutations present in their respective population, but some isolates also contained mutations that were not fixed in the population. Moreover, only two isolates (PEi1 and PEi3) contained a single mutation, limiting our ability to directly correlate individual mutations with phenotypes. Based on our results, we predict that the sole mutations present in PEi1 and PEi3 in promoter regions upstream of the genes annotated as pdeH and a hypothetical protein, respectively, are responsible for the increased biofilm production observed in these two isolates. We also speculate that the large 66 bp deletion in rcsB present in the MEi isolates contributed to the attenuated virulence in bees, as seen in MEi1 and MEi3 isolates. Furthermore, based on a previous study that demonstrated a role of the BarA-UvrY two-component system in controlling the carbon storage regulatory system in E. coli (Pernestig et al. 2003), it is probable that the mutations we detected in barA and uvrY in the MEi and PEi isolates are the result of adaptation to the Neff media, which contains 5% glucose. We caution that these are only speculations that require further investigation to be validated.
The five phenotypes that we evaluated were growth, biofilm production, grazing resistance, virulence in honey bees, and cytotoxicity to murine macrophages. We found no change in population density (growth) based on OD readings over 24 h in the three media types tested. However, the number of colony-forming units present after 18 h of growth on Neff agar was significantly lower for two of the ME isolates, MEi1 and MEi3. If the proliferation of MEi1 and MEi3 isolates is also slower within the honey bee host, this could potentially explain the observed decreased virulence of these two isolates in honey bees. However, we did not quantify growth with the honey bee host in this study, so this hypothesis warrants further investigation. When evaluating biofilm production, we found that all PEi isolates displayed increased biofilm production on average when compared with the ancestor, although PEi2 was not found to significantly produce more biofilm. The MEi isolates produced less biofilms than the PEi and ancestor on average, but individually, none of them significantly produced less biofilm than the ancestor. As biofilm production is a natural defense mechanism against predators (Flemming and Wingender 2010; Wucher et al. 2021; Hoque et al. 2023), increasing biofilm production could be beneficial in the constant presence of a predator, whereas it is likely not essential in a nutrient-rich monoculture. Although we did not see a significant increase in grazing resistance in any of the PEi isolates, we did observe a decrease in all three MEi isolates when compared with the ancestor. It is important to note that in order to maintain T. thermophila in coculture with KZ19 for 24 h at the start of the experiment, we had to start with a ratio of only 1 (bacteria:protist) cell ratio, suggesting that the ancestral strain already possessed mechanisms to resist grazing and/or survive phagocytosis. Thus, it is possible that in addition to the pressure to survive predation, competition for resources was a major force impacting the evolution of S. marcescens in the presence of T. thermophila in our experiments. Evidence for this hypothesis is supported by the fact that the isolates that exhibited increased biofilm formation did not have significantly increased grazing resistance. Taken together, our results indicate that biofilm production promotes grazing resistance in S. marcescens KZ19 but also suggest that increasing biofilm production could be providing other functions in the presence of a predator aside from protection, such as the ability to better compete for resources (Oliveira et al. 2015; Rendueles and Ghigo 2015).
We found that growth, biofilm production, and grazing resistance were positively correlated with virulence in a natural host (the honey bee). The correlation between grazing resistance and bee pathogenicity was particularly strong, especially given our sample size, providing support for the idea that predation can indirectly select for host virulence (Levin and Edén 1990; Barker and Brown 1994; Levin 1996; Molmeret et al. 2005). Conversely, we found no correlation between cytotoxicity to murine macrophage cells (RAW 264.7) and growth, biofilm production, grazing resistance, or virulence in bees. In fact, three of the PE isolates were significantly less toxic to macrophage after 2 h of exposure, with PEi5 displaying virtually no cytotoxicity in this assay. Based on our time series experiment, we found that cytotoxicity to RAW 264.7 macrophages is highly time and density dependent. After 6 h of coculture, all isolates and the ancestral strain converged toward high cytotoxicity and had killed all macrophages in the culture. Other strains of S. marcescens have been shown to be highly toxic to macrophages (Ishii et al. 2012; Krzymińska et al. 2012, 2010), and our findings confirm this result for S. marcescens KZ19. The lack of correlation between macrophage toxicity and all the other phenotypes we tested suggests that the mutations present in our isolates have no impact on their ability to kill murine macrophage RAW 264.7 cells and, in some rare cases, might render them less toxic. These results indicate that although evolution in the presence of a predator may increase bacterial virulence, the increase in virulence could be host specific, likely due to the diversity of mechanisms involved in virulence. This underlines the importance of testing the virulence of PE bacteria in multiple systems.
A few other studies have used experimental evolution to directly test how predation by T. thermophila impacts the virulence of S. marcescens (Friman et al. 2009; Mikonranta et al. 2012; Zhang et al. 2014). Contradictory to our findings, all these studies found that predation attenuated virulence. Some reasons for inconsistencies between our study and previous studies could be due to (i) the use of different S. marcescens strains with different life histories, (ii) coevolving S. marcescens and T. thermophila together (Friman et al. 2009; Mikonranta et al. 2012; Zhang et al. 2014) rather than preventing the evolution of the predator (this study), (iii) using different hosts to test virulence, and (iv) differences in the methods used to evaluate virulence. Virulence of the PE S. marcescens in previous studies was evaluated based on testing either the entire populations (Friman et al. 2009) or the individual isolates that were then pooled to evaluate virulence (Mikonranta et al. 2012; Zhang et al. 2014). By testing a mixture of evolved isolates, it is likely that one will outcompete the others (and that the “winner” will be different depending on the environment). Combining or comparing data from assays performed with different isolates could strongly impact the results since each isolate might genotypically and phenotypically differ (as we observed here) and be more or less fit depending on the context. Moreover, none of these studies evaluated the genotypes of the evolved populations, limiting their ability to fully evaluate how predation impacts the evolution of S. marcescens. Thus, we believe that the differences in experimental design between the previous predator S. marcescens evolution studies and ours likely explains the discrepancy in the results.
Fundamentally, our study demonstrates the efficacy of using experimental evolution to identify genes and pathways involved in virulence-associated phenotypes, even over short evolutionary timescales. Despite some variability in the evolutionary outcomes in our experiment, our results provide overall support for the hypotheses (i.e. CEH and the training grounds) that propose that bacterial virulence is selected for or maintained by predation. Performing longer evolutionary experiments under different scenarios and with different opportunistic bacterial strains and species will provide a more accurate picture of the forces that drive the evolution of opportunistic pathogens and the mechanisms responsible for their virulence. Identifying virulence mechanisms and understanding how virulence evolves can help guide the development of more effective treatment strategies for combating bacterial infections.
To confirm that T. thermophila grazes on S. marcescens KZ19, we utilized transconjugant KZ19 containing an E2 crimson fluorescent protein with spectinomycin resistance (Leonard et al. 2018; Raymann et al. 2018). The transconjugant KZ19 strain was inoculated into 3 mL of LB broth containing 180 µg/mL spectinomycin and incubated at 30 °C for 24 h. After 24 h, 1OD of the transconjugant KZ19 was resuspended in fresh Neff media and 1 mL was pipetted into 3 mL of Neff media containing T. thermophila and 180 µg/mL spectinomycin. The culture was incubated for 24 h at 30 °C. One milliliter of the 24 h coculture was then centrifuged at 2000 × g for 10 m to collect the T. thermophila cells, the supernatant containing bacteria was discarded, and the T. thermophila were resuspended in 1 mL of axenic Neff media. The T. thermophila cells were then pipetted onto a microscope slide that was washed with 60% ethanol (EtOH; which immobilizes the cells without immediately lysing them), covered with a glass cover slip, and imaged on a Keyence BZ-X700 series all-in-one fluorescence microscope in brightfield and with a Cy5 filter.
Serratia marcescens strain KZ19 was evolved either in media alone (ME) or with the predator T. thermophila strain SB210 (PE). Tetrahymena thermophila strain SB210 was purchased from the Tetrayhema Stock Center located at the Cornell University. The ME experiment included three replicate lines, and the PE experiment included six replicate lines. First, T. thermophila and the ancestral S. marcescens KZ19 were cultured individually for 72 and 24 h, respectively, at 30 °C in Neff media. To start each evolved line, 1OD of the overnight KZ19 culture was resuspended in fresh Neff media and 1 mL was added to either 3 mL of Neff media alone (MEp lines) or 3 mL Neff media containing T. thermophila culture (PEp lines) at a ratio of 1 cells (bacteria:protist). The starting ratio of S. marcescens:T. thermophila was based on preliminary trials starting with various concentrations of S. marcescens and T. thermophila coinoculated in Neff media to determine the ratio that would allow both organisms to remain in coculture for 24 h.
The 3 mL experimental cultures were incubated at 30 °C without shaking. Every 24 h for 60 d, the cultures were vortexed, and 1% of each culture was sampled from the flasks and added to 3 mL of either (i) fresh Neff media or (ii) new axenic cultures of T. thermophila in Neff media. Every day, prior to passaging, the cultures were checked for T. thermophila viability and density via microscopy and the successful passaging and growth of S. marcescens presence were evaluated by plating on LB agar. Every day over the course of the 60 d evolution experiment, T. thermophila and S. marcescens were always detected. Although population size was not directly measured at every passage, based on daily plating (S. marcescens) and light microscopy (T. thermophila), we never observed a noticeable difference in population size. The number of generations reached for each S. marcescens evolved population over the course of 60 d is estimated to be approximately 396.
On day 60, six LB agar plates were streak inoculated from the culture flasks and incubated for 24 h at 30 °C. One bacterial colony was randomly selected from each of the six the PE populations (PE1 to PE6) and three of the ME populations (ME1 to ME3). Pure cultures of each isolate were created by inoculating in 3 mL LB broth with 24 h incubation at 30 °C. Samples (200 μL from each pure culture) were then frozen at −80 °C in 20% glycerol. These nine isolates were sequenced, analyzed, and used for phenotypic analysis (i.e. virulence in bees, growth in media, biofilm production, macrophage cytotoxicity, and grazing resistance).
Cultures of the evolved populations (ME n = 6, PE n = 6), the evolved isolates (ME n = 3, PE n = 6), and the ancestor were diluted to 1OD. DNA extractions were performed using the Zymo Quick-DNA Fungal/Bacterial Miniprep Kit (D6005). DNA was prepped for sequencing using the Oxford Nanopore Rapid Barcoding Sequencing gDNA Kit (SQK-RBK004) for long-read sequencing and the Illumina Nextera DNA Flex Library Prep Kit (20018704) for short-read sequencing. The ancestral genome was sequenced via long (Oxford Nanopore Minion)- and short (Illumina iSeq100)-read technologies. The evolved populations and isolates were sequenced via short-read technology on an Illumina iSeq100 with 2 × 150 paired end reads (see supplementary table S1, Supplementary Material online, for sequencing depth details).
The Illumina and Nanopore reads of the ancestral genome were then used to generate a hybrid assembly using hybridSPAdes v3.15.3 (Antipov et al. 2016). Assembled scaffolds were annotated using Prokka v1.14.5 (Seemann 2014) on the Department of Energy Systems Biology Knowledgebase (KBase) platform (Arkin et al. 2018). We also mapped all short reads of the ancestral strain back to the consensus ancestral genome using breseq v0.38.1 (Deatherage and Barrick 2014). Any mutations predicted by breseq when the ancestral reads were mapped to the ancestral consensus genome were considered sequencing or assembly issues in the ancestral genome and removed as candidate mutations in the evolved populations and isolates.
The sequencing reads of the evolved populations and isolates were trimmed using Trimmomatic v0.40 (Bolger et al. 2014). The trimmed reads were then mapped to the newly assembled ancestral genome using breseq v0.37.0 (Deatherage and Barrick 2014) with a base-quality cutoff PHRED score of 30 and minimum frequency cutoff of 0.05 to identify mutations. Promoters were predicted using iProEP (Lai et al. 2019). Functional characterization of genes (Fig. 2; supplementary table S2, Supplementary Material online) was predicted using KOALA KEGG Orthology and Links Annotation (Kanehisa and Goto 2000; Kanehisa et al. 2016). Protein–protein interaction network prediction was determined using the STRING v12 database multiple proteins search by querying the protein name of all ten of genes in which we identified fixed parallel mutations across the evolved populations (Szklarczyk et al. 2023) using S. marcescens, S. rubidaea, and Y. enterocolitica as reference organisms. Simulations to determine the percent mutations expected at random across intergenic, synonymous, and nonsynonymous positions and the probability of gene- and site-specific parallel mutations based on random expectation were conducted with a custom Python we conducted 10,000 independent simulations where 1 mutation was introduced at random under a Jukes and Cantor model in the reference genome of S. marcescens KZ19 (GCA_002915435.1). The annotations of the genome were used to infer whether the mutation was intergenic or not. Mutations affecting coding sequences were inferred as synonymous or nonsynonymous by translating the gene in silico before and after introducing the mutation.
Growth assays were performed in triplicate in among LB broth, fresh Neff media, and spent Neff media with T. thermophila filtered out (i.e. T. thermophila were grown in the media for 48 h before filtering). Pure 1OD cultures of the evolved isolates and the ancestor were diluted to 10^−5^, and 5 μL of the culture was pipetted into 200 μL of media in a 96-well microplate. Absorbance detection was performed at an OD of 600 nm at 30 °C with shaking every hour for 24 h on a BioTek Synergy two-plate reader.
Growth was also evaluated by counting colony-forming units. Pure 1OD cultures of each of the evolved isolates and the ancestor strain were diluted to 10^−6^, 10^−7^, and 10^−8^, and 1 mL of each dilution was plated in triplicate onto Neff agar plates using glass beads. The plates were incubated at 30C for 18 h, and then, colony-forming units were counted from the 10^−7^ dilution plates and multiplied by the dilution factor to obtain the number of colony-forming unit per milliliter.
Biofilm assays were performed in triplicate for each isolate. Pure 1OD cultures of the evolved isolates and the ancestor were diluted to 10^−5^, and 5 μL of the culture was pipetted into 200 μL of media in a 96-well microplate. The covered plate was incubated for 36 h at 30 °C without shaking. After 36 h, the supernatant was poured off, the wells were washed twice with 200 μL of sterile ddH2O, and the cells were stained with 150 μL of 0.04% crystal violet for 10 min. After 10 min, the excess crystal violet was removed, and the wells were washed again and stained once more for 10 min. The crystal violet was poured off, and cells were solubilized with 150 μL of 95% EtOH. OD readings were taken at absorbance reading 560 nm on the BioTek Synergy two-plate reader. A total of four assays were performed, each in triplicate.
Grazing resistance assays were performed in triplicate for each isolate. Cultures of the evolved isolates and the ancestor were grown overnight, normalized the 1OD, and 1 mL was added to a culture flask containing T. thermophila in 3 mL Neff media (10:1 ratio of bacteria to protist cells). Tetrahymena thermophila alone in 3 mL Neff media served as the negative control. Experimental culture flasks were incubated at 30 °C for 48 h. After 48 h, the cultures were mixed thoroughly and 1 mL was sampled and centrifuged at 5,000 rpm for 10 min. The supernatant was discarded, pellets were resuspended in 1 mL 1× phosphate-buffered saline (PBS), and the solution was filtered through 5.0 μm MilliporeSigma Millex-SV Sterile PVDF syringe filters to remove the T. thermophila cells. OD readings were taken of filtrates at 600 nm after blanking with the negative control filtrate. A total of two assays were performed, each in triplicate.
Pure 1OD cultures of the evolved isolates and the ancestor were pelleted via centrifugation and resuspended in a 1 sterile sugar syrup solution (SSS). Conventional adult honey bee workers (A. mellifera) were sampled from a single hive located on the Gateway North Research Campus in Browns Summit, North Carolina. Bees were immobilized at 4 °C, randomly distributed into groups, and exposed via the immersion method (Raymann et al. 2018) to ∼10 μL of one of three (i) sterile SSS only, (ii) the ancestor in sterile SSS, or (iii) each evolved isolate in sterile SSS. Bees were kept in cup cages (20 bees/cup; 100 bees/treatment) under hive 35 °C and 95% humidity. Mortality between the treatments was noted every day for 5 d. A total of 3 assays were performed with 5 replicates per assay equaling 300 bees tested in total per isolate. A survival curve (Kaplan–Meier) was created in GraphPad Prism v9.1.0.
Murine macrophage cells (Raw 264.7) were maintained in Dulbecco's Modified Eagle Medium (DMEM, ATCC: 30-2002) with 10% fetal bovine serum (FBS, ATCC: 30-2020) at 37 °C with 5% CO2. After reaching 90% confluency, cells were harvested, resuspended in 500 µL of DMEM media in a Corning Costar 24-well Clear TC-treated well plate (Fisher Scientific: 09-761-146), and incubated overnight at 37 °C with 5% CO2. After overnight incubation, the old DMEM media were discarded, the macrophage cells were washed with 1× PBS, and 500 μL of fresh DMEM media was added to each well. Overnight cultures of each bacterial isolate were normalized to 1OD prior to the experiment, and 10 mL of each 1OD bacterial culture (MOI 1) was added into a separate well containing macrophage and incubated at 37 °C with 5% CO2. For control wells (macrophage spontaneous release), 10 μL of 1× PBS was added instead of bacteria. The cytotoxicity was measured after 2 h using the CytoTox 96 Non-Radioactive Cytotoxicity Assay. In brief, 2 h after coculturing, the media from each well (including controls) were mixed two times and 60 μL was pipetted into a 1.7 mL tube and centrifuged at 10,000 rpm for 2 min to pellet bacterial cells. After centrifuging, 50 μL of the supernatant was transferred into a new 96-well plate (Fisher Scientific: 07-200-90) to measure LDH release (OD490). For control wells, the macrophage lysis buffer release was measured by adding 50 μL 10× lysis buffer to macrophage samples 45 min prior to LDH measurement, and the negative control was measured without adding lysis buffer. The percent cytotoxicity was calculated using the equation below. The ratio of 55/51 was used to correct the volume change caused by the lysis buffer. A single assay was performed with six replicates.
Since virtually all isolates were highly cytotoxic to macrophage after 2 h at MOI of 1 cells (bacteria:macrophage), we performed a time series experiment to determine if macrophage cytotoxicity is time and density dependent. The Raw 264.7 macrophages and bacterial cultures were prepared as described above, except this time the assay started with a ratio of 10 (MOI) bacteria to macrophage cells and cytotoxicity via LDH release (CytoTox 96 Non-Radioactive Cytotoxicity Assay) was measured at five time points following coculture (2, 3, 4, 5, and 6 h). Cytotoxicity in this assay was based on OD at OD490 rather than percent cytotoxicity. At 6 h, all macrophages cocultured with the ancestral and evolved isolates were dead as indicated by cell lysis observed under an inverted light microscope. One assay was performed in triplicate.
We would like to thank Dr. Carlos Goller for assisting us with the Nanopore sequencing, Dr. Stephanie Shames for providing us with the Raw 264.7 macrophages, and Dr. Louis-Marie Bobay for helping us with the simulations and providing constructive and helpful feedback on the manuscript.
Heather A Hopkins, Department of Plant and Microbial Biology, North Carolina State University, Raleigh, NC, USA; Department of Biology, University of North Carolina Greensboro, Greensboro, NC, USA.
Christian Lopezguerra, Department of Plant and Microbial Biology, North Carolina State University, Raleigh, NC, USA; Department of Biology, University of North Carolina Greensboro, Greensboro, NC, USA.
Meng-Jia Lau, Department of Plant and Microbial Biology, North Carolina State University, Raleigh, NC, USA.
Kasie Raymann, Department of Plant and Microbial Biology, North Carolina State University, Raleigh, NC, USA; Department of Biology, University of North Carolina Greensboro, Greensboro, NC, USA.
Supplementary material is available at Genome Biology and Evolution online.
H.A.H., C.L., M.-J.L., and K.R. performed the experiments and collected the data. K.R. designed and funded the research. All authors contributed to the data analysis and manuscript writing and editing.
This work was supported by the National Science Foundation under grant DEB-2344788 (to K.R.), the National Institutes of Health under grant 7R01GM145747-02 (to K.R.), National Science Foundation Graduate Research Fellowship under grant DGE-2137100 to C.L., the University of North Carolina Greensboro (UNCG) College of Arts and Sciences Faculty First Award (to K.R.), and the UNCG Department of Biology Graduate Student Support grants (to H.A.H and C.L.).
The data underlying this article are available in the NCBI Sequence Read Archive (SRA) at https://www.ncbi.nlm.nih.gov/, and can be accessed with BioProject identifiers PRJNA432218 and PRJNA838621.
The data underlying this article are available in the NCBI Sequence Read Archive (SRA) at https://www.ncbi.nlm.nih.gov/, and can be accessed with BioProject identifiers PRJNA432218 and PRJNA838621.