Adaptive Resistance Mutations at Suprainhibitory Concentrations Independent of SOS Mutagenesis.

Emergence of resistant bacteria during antimicrobial treatment is one of the most critical and universal health threats. It is known that several stress-induced mutagenesis and heteroresistance mechanisms can enhance microbial adaptation to antibiotics. Here, we demonstrate that the pathogen Bartonella can undergo stress-induced mutagenesis despite the fact it lacks error-prone polymerases, the rpoS gene and functional UV-induced mutagenesis. We demonstrate that Bartonella acquire de novo single mutations during rifampicin exposure at suprainhibitory concentrations at a much higher rate than expected from spontaneous fluctuations. This is while exhibiting a minimal heteroresistance capacity. The emerged resistant mutants acquired a single rpoB mutation, whereas no other mutations were found in their whole genome. Interestingly, the emergence of resistance in Bartonella occurred only during gradual exposure to the antibiotic, indicating that Bartonella sense and react to the changing environment. Using a mathematical model, we demonstrated that, to reproduce the experimental results, mutation rates should be transiently increased over 1,000-folds, and a larger population size or greater heteroresistance capacity is required. RNA expression analysis suggests that the increased mutation rate is due to downregulation of key DNA repair genes (mutS, mutY, and recA), associated with DNA breaks caused by massive prophage inductions. These results provide new evidence of the hazard of antibiotic overuse in medicine and agriculture.

Introduction

Treatment failure of infectious diseases due to antimicrobial resistance is a major threat to human health (Ventola 2015). Antimicrobial resistance is generally the outcome of genetic modifications, obtained through advantageous mutations on existing genes or via gaining resistance-associated genes through horizontal gene transfer events (Blair et al. 2015; Piddock 2017). As the former mechanism is independent of external DNA source, it is of great interest to characterize the intrinsic capacity of bacteria to overcome antibiotic stress by the accumulation of resistance mutations, as may occur during therapy. Such a scenario is particularly important for those drugs in which resistance can be acquired through a single point mutation on an existing gene, for example, rifampicin resistance conferred by mutations on the RNA polymerase β-subunit (rpoB) gene (Goldstein 2014) or streptomycin resistance mediated by mutations in the ribosomal protein S12 (rpsL) and 16S rRNA (rrs) genes (Springer et al. 2001). Insights in microbiology indicate that microbes can respond to certain stresses by increasing their mutation rate, leading to an increased stress-responsive manner of evolution (Rosenberg 2001; Hersh et al. 2004; Cirz and Romesberg 2007; Foster 2007). In addition, some microorganisms can harbor subpopulations of cells that grow slowly in the presence of a stressor (El-Halfawy and Valvano 2015), potentially promoting the emergence and selection of beneficial mutations and thereby increasing the effective population size without the need to increase the mutation rate. The implications of such processes during antibiotic exposure could explain treatment failure and the high levels of drug resistance seen in bacteria of clinical importance (Friedman et al. 2016).

The in vitro exploration of microbial sensitivity and resistance profiles has led to major discoveries on the manner microorganisms respond to antibiotic stress. The emergence of microbial colonies within the inhibition zone of antibiotic diffusion tests, namely cells growing at apparent suprainhibitory antibiotic concentrations, has been reported in the literature since the mid-twentieth century (Gómez Lus 1959), and has been described for different antibiotic classes and many pathogenic species, including bacteria and fungi (Khan et al. 2008; Gomes et al. 2016; Sindeldecker et al. 2020). Potential explanations proposed that these colonies emerge due to: 1) spontaneous resistant mutations acquired prior to antibiotic exposure and their selection during the test (Macia et al. 2004; Drlica and Zhao 2007); 2) hetereoresistance, as a capacity of a pre-existing subpopulation to grow at a higher antibiotic concentration, usually through transient genetic or nongenetic physiological mechanisms other than the typical resistance mutations (El-Halfawy and Valvano 2015; Nicoloff et al. 2019); or 3) antibiotic degradation, allowing tolerant or persister cells to grow once subinhibitory concentrations are reached (Gefen et al. 2017). Whereas isolation of true genetic resistant clones could be explained with the first scenario, isolation of classical heteroresistant or persister cells do not, as they replicate the sensitivity profiles of the ancestral clone (El-Halfawy and Valvano 2015). Nevertheless, since wild-type bacteria have spontaneous mutation rates ranging from 10−8 to 10−11 bp/generation (Lynch 2010; Pan et al. 2021), a standardized preparation of antibiotic diffusion agar tests (with ∼106 cells plated per plate) (Institute CaLS 2017) reduces the chances of pre-existing mutants being selected, that is, adaptation from standing genetic variation, and raises the possibility that in some cases the emergence of antibiotic resistance may take place during antibiotic exposure on the test plates. Hence, while avoiding selection of preexisting mutants and delineating the antibiotic concentrations during incubation, the phenomenon of emergence of true resistant colonies within the inhibition zone of agar diffusion tests offers a unique model to explore the microbial capacity to respond and adapt to antibiotic therapy.

In case antibiotic degradation occurs during the incubation of agar diffusion tests, subinhibitory antibiotic concentrations could be reached, allowing surviving bacteria to react and potentially elicit mutagenic responses (Kohanski et al. 2010; Cairns et al. 2017). For example, sublethal antibiotic concentrations can promote the formation of mutations by stimulating the production of reactive oxygen species (ROS) that damage the DNA (Kohanski et al. 2010; Sebastian et al. 2017). However, persister bacteria that resume their growth within the inhibition zone on agar surfaces after replacing the antibiotic disks to nutrient disks, do not immediately evolve resistance despite the remaining sublethal drug conditions (Gefen et al. 2017). Contrastingly, if antibiotics are maintained at suprainhibitory concentrations, the evolution of resistance will have to occur in tolerant/surviving bacteria that have an exceptional capacity to overcome the stress and to accumulate mutations within the specific “resistance-conferring” locus. Such a scenario could be initially promoted by heteroresistant subpopulations, replicating under these conditions. Nevertheless, under conserved mutation rates (∼1–3×10−3 mutations per genome per generation) (Sung et al. 2012), the chances that a resistant mutation is acquired during the formation of a heteroresistant colony is extremely low, as a colony represents no more than 30 generations (Lee et al. 2012). Therefore, resistance emergence from sensitive cells during antibiotic exposure requires stress-associated mechanisms for mutagenesis.

It has been shown that de novo resistance mutations can be acquired during antibiotic stress (Cirz and Romesberg 2006). This acquisition requires specific stress responses that induce mutagenic mechanisms. In classical bacterial models (e.g., Escherichia coli and Bacillus subtilis), antibiotic stress can lead to mutagenic states, activated by the SOS mutagenic response (SOS response, led by RecA/LexA) or by other stress-induced mutagenesis mechanisms (RpoS-dependent) (Riesenfeld et al. 1997; Rosenberg 2001; Sung et al. 2003; Cirz and Romesberg 2006). Most of the mutagenic states described to date are characterized by the expression of Y family error-prone DNA polymerases (Foster 2007), which were first identified by their ability to perform mutagenic translesion DNA synthesis (Ohmori et al. 2001) and later to allow the replication to resume after double-strand DNA breaks (Ponder et al. 2005). Nevertheless, the availability of thousands of bacterial genomes has revealed that Y family error-prone polymerases are not universal, raising the question whether alternative pathways might facilitate evolution under stress in these bacteria. Alternative responses could include the derepression of essential DNA repair genes, as observed by the downregulation of MutS and MutH in an RpoS-dependent manner in E. coli (Feng et al. 1996; Tsui et al. 1997), or by inaccurate transcription or translation events, which include the production of altered proofreading subunits of the replicative polymerases (Slupska et al. 1998; Al Mamun et al. 2002; Balashov and Humayun 2002). Therefore, it is important to investigate the evolution of antibiotic resistance on unexplored bacterial phylogenies that lack typical mutagenic mechanisms.

Slow-growing bacteria can be intrinsically tolerant to typical bactericidal antibiotics as these drugs usually kill actively dividing cells (Kaiser et al. 2014; Brauner et al. 2016). Bartonella (Alphaproteobacteria) are slow-growing pathogenic bacteria, with a high degree of genetic diversity, distributed among mammalian and arthropod hosts (Harms and Dehio 2012; Okaro et al. 2017). The characterization of dozens of Bartonella genomes have shown that most of them lack many known genetic systems associated with stress response and mutagenesis, such as error-prone polymerases (Pol II, Pol IV or Pol V or Dnae2), RecBCD, RpoS stress-response regulator, among others (Saenz et al. 2007; Berglund et al. 2009; Gutiérrez et al. 2020). Recently, we have fully characterized the genomes of two Bartonella species, Bartonella krasnovii and Bartonella kosoyi, and estimated their spontaneous mutation rates under benign conditions (Gutiérrez et al. 2018, 2020). Although these Bartonella spp. exhibit conserved mutation rates, measured by fluctuation tests and mutation accumulation experiments (Gutiérrez et al. 2018), we present here evidence of their outstanding capacity to evolve resistance during antibiotic exposure, identified as colonies emerging within the inhibition zone of agar diffusion tests. We explored this phenomenon through experimental assays and a mathematical model, and delineated the conditions and bacterial responses leading to their ability to become resistant. Collectively, this study presents evidence of a microbial adaptive evolution in the absence of the classical stress-induced mutagenic complexes.

Results

Emergence of True Resistant Mutants within the Inhibition Zone of Antibiotic Diffusion Tests

While evaluating the minimum inhibition concentration (MIC) of B. krasnovii wild-type strain using rifampicin (RIF) Etest (BioMériux, Marcy-l’Étoile, FR), a discrete number of bacterial colonies emerged within the inhibition zone and near the growth border (supra-MIC area) of chocolate agar plates after prolonged incubation (two to seven days after the culture lawn was visually evident, fig. 1). Sequential isolation of the colonies that grew within the inhibition zone and re-evaluation of the MIC by Etest and agar dilution methods revealed that the isolates had MIC values ranging from 0.064 to 512 µg/ml (supplementary table 1, Supplementary Material online) which is 3- to 22,260-fold higher than the original RIF MIC (MICoriginal = 0.023 µg/ml; fig. 1). The same phenomenon was recorded using in house RIF diffusion disks (25 µg) and on alternative culture media, refuting the possibility of a test-related effect (fig. 1). Notably, most of the colonies (∼95%) isolated from the inhibition zones had identical incubation times as compared with the original clones (average 4 days to reach a >1 mm colony), and only in a few cases was a delay in the growth recorded (up to 7 days). Moreover, the resistant phenotypes and their MICs were fully stable, as demonstrated after five sequential passages on RIF-free agar plates and MIC re-evaluation by Etest, indicating a heritable trait. The emergence of resistant colonies within antibiotic diffusion tests was later confirmed with other Bartonella spp. and strains, including B. kosoyi, B. elizabethae, B. grahamii, and B. henselae.

Fig. 1.

Rifampicin resistant mutants are generated within the inhibition zone of antibiotic diffusion tests after prolonged incubation. (A) Bartonella krasnovii wild-type strain showed the emergence of colonies (indicated by yellow arrows) within the inhibition zone produced by rifampicin (RIF) Etest after prolonged incubation. (B) The colonies that grew within the inhibition zone presented higher resistance phenotype after isolation and re-evaluation of their minimal inhibition concentration (MIC) on subsequent Etests. (C) RIF resistant colonies emerged on chocolate and/or blood agar plates containing in house RIF disks (RD, 25 µg RIF). (D) Resistant colonies are not preexisting and only emerge on RIF diffusion disk plates (RD), whereas plates with homogeneous RIF concentrations inhibited the emergence of resistant colonies.

As single mutations on the RIF-resistant determinant region (RRDR) of the RNA polymerase β subunit gene (rpoB) are known to confer >95% of the resistance to RIF (Telenti et al. 1993), we sequenced the RRDR locus from 63 colonies from different Bartonella spp. isolated from the inhibition zone of Etests or disk-diffusion tests. The latter resulted in the identification of 15 different mutation types within 12 positions of the RRDR (table 1). All RIF-resistant colonies had one mutation in the RRDR. Remarkably, different mutation types were commonly identified from colonies isolated from the same agar plate, suggesting that independent mutational events had occurred. Finally, in two cases, we identified sequential mutations at the RRDR from colonies that emerged from consecutive RIF diffusion tests (supplementary table 2, Supplementary Material online).

Table 1.

Nonsynonymous Mutations on the rpoB Gene and Their Associated MIC Recorded from Different Resistant Bartonella spp. Colonies Isolated from the Inhibition Zone of Antibiotic Diffusion Tests.

rpoB Position	Reference	Mutation	AA Effect (position)	MICmutant (µg/ml)a	Times Recorded	Species
1574	T	C	L to S (525)	0.125	2	B. krasnovii b
1577	C	T	S to L (526)	0.032	1	B. henselae c
1579	C	A	Q to K (527)	128	1	B. krasnovii
1588	G	A	D to N (530)	0.75	11	B. krasnovii
1618	C	T	H to Y (540)	≥32	11	B. krasnovii, B. elizabethae d
1627	C	T	R to C (543)	0.5–1	5	B. krasnovii
1628	G	A	R to H (543)	64–128	2	B. krasnovii
1634	C	A	S to Y (545)	256	3	B. krasnovi, B. grahamii e
1634	C	T	S to F (545)	256	6	B. krasnovii, B. elizabethae, B. grahamii
1642	G	T	G to C (548)	0.125–1	7	B. krasnovii, B. kosoyi f
1733	C	T	P to L (578)	0.75–1	5	B. krasnovii, B. elizabethae
1733	C	A	P to H (578)	0.25	1	B. elizabethae
1758	T	G	I to M (572)	0.094	1	B. kosoyi
1763	C	A	S to Y (588)	0.064	1	B. krasnovii
1763	C	T	S to F (588)	0.38	6	B. krasnovii

MIC estimated by Etest and/or agar dilution methods.

Bartonella krasnovii MICoriginal=0.023 µg/ml.

Bartonella henselae MICoriginal=0.006 µg/ml.

Bartonella elizabethae MICoriginal=0.008 µg/ml.

Bartonella grahamii MICoriginal=0.008 µg/ml.

Bartonella kosoyi MICoriginal=0.008 µg/ml.

The emergence of colonies within the inhibition zone was further studied on B. krasnovii, as this species was recently characterized pheno- and geno-typically, and its spontaneous mutation rates were estimated by fluctuation tests and mutation accumulation experiments (Gutiérrez et al. 2018, 2020). First, the phenomenon was explored using additional antibiotics (supplementary fig. 1A, Supplementary Material online). We tested antibiotics for which resistance can be conferred by point mutations on existing genes (streptomycin [Springer et al. 2001], enrofloxacin and ciprofloxacin [Hooper et al. 1987]), antibiotics for which resistance involves predominately the acquisition of specific genes (chloramphenicol and doxycycline; Schwarz et al. 2004; Grossman 2016) as negative controls, and mitomycin C (a mutagenic drug) as a positive control (Tsukamura and Tsukamura 1962). Results showed that only streptomycin and mitomycin C disks resulted in clear isolated colonies within the inhibition zone (supplementary fig. 1B and C, Supplementary Material online). All colonies that emerged from the inhibition zone of streptomycin disks showed MICs higher than the starting strain in subsequent tests, and mutations on the rpsL gene were identified (supplementary table 3, Supplementary Material online). All mitomycin C-derived colonies showed incremental increases in resistance relative to the original population, as visualized by the smaller inhibition zone detected upon retesting (supplementary fig. 1D, Supplementary Material online). In the case of ciprofloxacin, colonies were observed in two occasions at the edge of the culture lawn of tested plates, but sequential isolation of these colonies resulted in identical MICs as the original populations. Colonies within the inhibition zone were never observed when enrofloxacin, chloramphenicol, or doxycycline disks were used. Thus, the emergence of resistant clones during antibiotic exposure occurs only with certain antibiotics, as observed with RIF and mitomycin C, although the frequency of colonies within the inhibition zone was the highest with RIF relative to the other antibiotics tested.

Resistant Mutants Are Generated during Antibiotic Exposure in an SOS-Independent Manner

To evaluate whether the resistant colonies that emerged within the inhibition zone of RIF diffusion tests resulted from selection of pre-existing spontaneous mutants (standing genetic variation) or were the outcome of de novo mutations generated during antibiotic exposure, we ran several experiments. First, we reduced the possibility of carrying spontaneous rpoB mutants by working with an extremely diluted fresh liquid culture (Ni = 4,200 ± 275 cells/100 µl) as a starting inoculum (fig. 1). Then, we split this culture into 100-µl subcultures and incubated 83 parallel subcultures for a short period to slightly increase the bacterial loads (Nf = 1.9 ± 0.44×106 cells/100 µl), which corresponded to 8.8 ± 0.38 cell divisions. Following the incubation, 20 µl of the 100-µl cultures were spread on chocolate agar plates (average 3.7 ± 0.87×105 cells/plate) with RIF diffusion disks (RD) or chocolate agar plates containing homogeneous concentrations of RIF, ranging from 0.04 to 0.1 µg/ml (representing 2- to 5-folds the MICoriginal). As shown in figure 1, on 95% (19/20) of the RD plates colonies grew within the inhibition zone, whereas zero colonies were recorded on the homogeneous-RIF containing plates (pairwise Wilcoxon test, P = 8×10−8). Despite an additional 10 days of incubation (20 days in total), plates with homogeneous RIF concentrations (0.04–0.1 µg/ml) showed no colonies. By contrast, out of 62 isolated colonies that grew within the inhibition zone of the RD plates, 92–97% were successfully isolated on 0.02–0.1 µg/ml RIF plates, respectively (supplementary table 4, Supplementary Material online). These results demonstrated that the studied phenomenon is associated with the diffusion of the antibiotic, promoting the generation of resistance postantibiotic exposure, whereas a direct exposure to suprainhibitory RIF inhibited or prevented the emergence of resistant mutants.

To understand the extent of this phenomenon, we serially diluted (10-fold dilutions) fresh Bartonella liquid cultures in PBS and seeded each dilution on five RD plates. After incubation, the number of colonies within the inhibition zone were recorded in all groups (supplementary fig. 2A, Supplementary Material online), including those inoculated with the lowest bacterial loads (6,800 cells per plate). Additionally, the number of colonies within the inhibition zone did not accumulate exponentially across dilutions (pairwise Wilcoxon rank sum test with Bonferroni correction, all P > 0.05; supplementary fig. 2B, Supplementary Material online), contrary to the expectation for pre-existing mutants in the population. The mutation frequencies, calculated as the number of resistant colonies per inoculated bacteria within each inhibition zone (as the area of this zone increases inversely proportional to the inoculum size; supplementary fig. 2A, Supplementary Material online) ranged from 1.3×10−6 to 6.2×10−4 mutants/bacteria, and showed a tendency to increase as the inoculum size decreases (supplementary fig. 2C, Supplementary Material online), resulting in 70 to over 34,000-fold higher than the estimated mutation rate (µ = 1.4×10−8rpoB mutations/generation) (Gutiérrez et al. 2018). These results suggested that an external factor must be driving the increased mutation rates.

The apparent increase in the observed mutation rates could be explained by an SOS-like induced mutagenesis. Although B. krasnovii does not encode the error-prone polymerases Umuc or DinB, it does carry the recA and lexA genes, which are required for SOS response. Therefore, we wondered whether UV exposure, a known SOS response trigger (Aksenov 1999) would promote the generation of RIF-resistant colonies at higher levels. Results from 30 parallel cultures exposed to the minimal lethal dose of UV (two seconds) and no-exposure controls did not reveal any significant difference in the number of colonies present within the inhibition zone of RD plates (supplementary fig. 3, Supplementary Material online). These results suggest that SOS mutagenesis does not play a role in the emergence of resistance during RIF exposure.

Bartonella Shows a Minimal Heteroresistance Capacity to RIF

Since the emergence of a discrete number of colonies within the inhibition zone of antibiotic diffusion tests has been considered as a trademark for heteroresistance (El-Halfawy and Valvano 2015), we evaluated the RIF-heteroresistance capacity of B. krasnovii by the gold standard method, the population analysis profile (PAP) test, and followed El-Halfawy and Valvano criterion (2015) for heteroresistance identification (see Materials and Methods). Starting with high-load B. krasnovii cultures (∼107–108 cells/ml), results from three independent experiments showed a marked decline in the capacity of B. krasnovii to generate colonies in agar plates with RIF concentrations ≥0.02 µg/ml, resulting in a narrow window of heteroresistance (only 6- to 8-fold, fig. 2). The RIF concentration required to inhibit the 99% of the population (MIC99) was 0.02 µg/ml, confirming estimations by Etests and agar dilution methods. At supra-MIC levels, B. krasnovii was only able to generate colonies on plates containing 0.04 µg/ml RIF (two times higher the MIC), representing a small heteroresistant subpopulation of about 0.0002% (two per 106 cells, fig. 2). The later was confirmed in subsequent tests of independent cultures seeded on 0.04 µg/ml plates, which gave colonies only when >106 cells were loaded (supplementary table 5, Supplementary Material online). Colonies (N = 4) that grew on the PAP 0.02 and 0.04 µg/ml plates, exhibited transient heteroresistant phenotype, reproducing the original MIC (MIC = 0.023 µg/ml) in a subsequent Etest (supplementary fig. 4, Supplementary Material online). According to El-Halfawy and Valvano criterion (2015), B. krasnovii present an intermediate heteroresistance capacity to RIF.

Fig. 2.

Bartonella krasnovii phenotypic response to RIF exposure. (A) Bartonella krasnovii population analysis profile (PAP) for RIF, indicating the minimal inhibition concentration for the 99% of the population (MIC99) and the heteroresistance window from the lowest antibiotic concentration giving maximum growth inhibition (0.005 µg/ml) to the highest noninhibitory concentration (0.04 µg/ml). (B) Dye-dilution test of B. krasnovii liquid cultures marked with eFluor450. After 5 days of incubation with different RIF concentrations at 37 °C (control at 4 °C), the bacteria were labelled with live/dead markers, Syto9 and propidium iodine (PI), and analyzed on ImageStream. eFluor450 intensities correspond to the survival population (Syto9+, PI−). (C) Death curve on RIF (0.2 µg/ml; 10X-MIC), indicating the minimal duration for killing the 99% of the population (MDK99=8 days). Analysis based on liquid cultures of B. krasnovii monitored for live bacteria every 24 h on antibiotic-free plates. Error bars indicate SD from the mean.

Since the PAP method relies on the ability of a bacterium to produce a macroscopic colony forming unit (CFU) at supra-MIC, which represents the outcome of 20–30 cell divisions from a single cell (Gutiérrez et al. 2018), we explored the possibility that the B. krasnovii population could undergo undetectable cell divisions at supra-MICs (less than required to generate a visual CFU). Accordingly, we performed a dye-dilution flow cytometry assay to measure cell divisions in individual Bartonella cells exposed to RIF (see Materials and Methods). This method is based on labeling cells with a fluorescent dye that binds intracellular amines (eFluor450, eBioscience, Invitrogen) that is diluted with each cell division (namely its final fluorescence will be inversely proportional to the cell divisions, supplementary fig. 5A and B, Supplementary Material online). As shown in figure 2, after 5 days of RIF exposure (0.02, 0.04, and 0.1 µg/ml), the surviving cells (Syto9-positive, as live marker) had average eFluor450 intensities 1.7- to 2.0-fold lower than the nonmultiplying control (parallel cultures incubated at 4 °C; t test with Bonferroni correction, P < 0.002; supplementary fig. 5C, Supplementary Material online). These results indicated that the B. krasnovii survivors underwent an average of one cell division during RIF exposure (halving the eFluor450 fluorescence level). Therefore, B. krasnovii exhibits a minimal population heteroresistance, with similar capacity on a range of supra-MIC RIF levels.

One of most frequent heteroresistance mechanisms known is gene amplification (Nicoloff et al. 2019). To evaluate if amplification of the rpoB gene played any role in the minimal RIF heteroresistance capacity of B. krasnovii, we quantified the relative copy number of rpoB gene (in relation to the citrate synthase gene) on heteroresistant colonies that multiplied on agar media with RIF (≥MIC levels), using qPCR assays. Results showed that rpoB relative copy numbers on pools of colonies exposed to RIF plates did not increase in comparison to the control cultures grown on RIF-free plates (2−ΔΔ = 1.004; supplementary table 6, Supplementary Material online). To confirm the latter, we excised single colonies from 0.02 and 0.04 µg/ml RIF plates and compared their rpoB-Cq values with single colonies excised from RIF-free plates. The number of rpoB copies in the colonies exposed to RIF concentrations was not higher than that of the control colonies (pairwise Wilcoxon rank sum test, all P > 0.05; supplementary fig. 6, Supplementary Material online), indicating that rpoB gene amplification was not responsible for B. krasnovii heteroresistance upon RIF exposure. We conclude that heteroresistance plays a minimal role on the studied phenomenon based on small subpopulation capable of producing colonies at supra-MICs, reduced population heteroresistance capacity and no rpoB amplification identified.

Bartonella krasnovii Exhibits High Tolerance to RIF

To further evaluate the B. krasnovii response to RIF exposure, we explored the survival rate of B. krasnovii by measuring the minimum duration to killing on liquid media containing RIF (0.2 µg/ml). The population of B. krasnovii exhibited a slow killing rate, with the minimal duration for killing of 99% of the population (MDK99) requiring approximately eight days (fig. 2). In fact, surviving cells were obtained even after 11-day postincubation on 0.2 µg/ml RIF liquid cultures. Notably, all the surviving clones maintained the MICoriginal phenotype after isolation (supplementary table 7, Supplementary Material online). Our results suggest that B. krasnovii presents a population-level tolerance response to RIF exposure with a large proportion of bacteria surviving supra-MIC of RIF during the first days of the agar experiment.

RIF Readily Accumulates and Is Maintained at Supra-MICs on the Inhibition Zone of Antibiotic Diffusion Plates

To test whether antibiotic degradation was playing any role on the emergence of colonies within the inhibition zone, RIF concentrations were measured by liquid chromatography tandem mass spectrometry (LC-MS/MS) directly from excised agar plugs from the inhibition zones (fixed distance of 2.0 cm from the center of the disk) over time (fig. 3). As shown in figure 3, the RIF concentrations reached MIC levels within 24 hours (average at 24 h = 0.029 ± 0.014 µg/ml), increased during the first four days, and then declined but remained at supra-MIC levels until the end of the experiment (average concentration of 0.07 ± 0.01 µg/ml after 10 days). We also identified the oxidized RIF metabolite, RIF quinone, during the LC-MS/MS analyses. Notably, RIF quinone is known to maintain its antimicrobial activity (Staudinger et al. 2014), a characteristic confirmed with RIF quinone disk assays in B. krasnovii (supplementary fig. 7, Supplementary Material online). To confirm the accumulation dynamics of RIF on these plates, fresh B. krasnovii cultures were exposed to RD sterilely “aged” plates (see Materials and Methods). In all time points, controls of fresh RD media were prepared in parallel with the same B. krasnovii inoculum, and after reincubation with the bacteria the radius of each plate inhibition zone was measured (fig. 3). If RIF degradation occurred during incubation, it was expected that the bacteria will produce smaller inhibition zones (smaller radius). Contrary to this assumption, in all cases the inhibition zone of the “aged” plates was larger than the fresh media controls, resulting in fold change of 1.1–1.5 of the controls radius size (fig. 3), and reaching a plateau after the third day. Interestingly, colonies emerged within the inhibition zone of “aged” plates, even after being sterilely incubated with RIF for 48 and 72 h, but always in the proximity of the culture lawn (supplementary fig. 8, Supplementary Material online). Altogether, these results confirmed that at the inhibition zone, where RIF resistant colonies emerge, RIF concentrations reaches supra-MIC after 24 h, accumulates during the first 3–4 days and is maintained at supra-MIC levels during the course of the study.

Fig. 3.

RIF degradation test during incubation. (A) Sterile chocolate agar plates with RIF diffusion disks (RD) were incubated at 37 °C + 5% CO2. At every time point (from 0.5 to 10 days), three plates were recovered, and three agar plugs were cut and processed per plate. RIF was extracted with acetonitrile and measured by HPLC. (B) RIF concentrations surpassed the MICoriginal (0.02 µg/ml, indicated by a blue-dash line) after 24 h, and maintained at supra-MIC for the rest of the experiment. (C) Confirmation assay of accumulation of RIF during time. Five sterile RD plates were incubated, recovered at different time points, and inoculated with fresh Bartonella krasnovii liquid cultures. Fresh RD plates were prepared in parallel and used as controls (N = 5). (D) After reincubation of aged and fresh inoculated RD, the inhibition zones were measured and compared. Fold change of the radius measured showed that aged plates always produced higher inhibition zones (above the FC = 1, red dash line). (E) Representative RD plates from the previous experiment (C and D) showing the increment of the inhibition zone in aged plates.

Whole-Genome Sequencing Analysis Reveals Single rpoB Mutations in Resistant Clones and Rejects Any Tolerance or Heteroresistance Epistasis Effect

To evaluate if multiple mutations accumulate in the colonies that emerge from the RIF inhibition zone, and if this was associated with a tolerant or heteroresistant genetic background, we performed the following controlled experiment. Four ancestral/original clones (confirmed as RIF sensitive by Etest) were fully sequenced along with seven derived resistant colonies isolated from the inhibition zones of the original RIF diffusion plates, seven tolerant clones recovered after 5 and 11 days from liquid cultures of the original clones exposed to RIF at 10-fold higher the MICoriginal and five heteroresistant clones isolated from PAP experimental plates (from 0.02 and 0.04 µg/ml RIF plates). In comparison to the original clones, nucleotide polymorphisms (SNPs) were only identified in the colonies isolated from the inhibition zone of RD plates (resistant clones). Each resistant isolate acquired only a single point mutation per genome, and it was always within the rpoB RRDR locus (table 2). By contrast, no mutations were identified in the tolerant and heteroresistant genomes. These results demonstrated that B. krasnovii tolerance and heteroresistance are not due to mutations, but most probably due to a physiological response, and that resistant clones acquired a single rpoB mutation during drug exposure without accumulating other nonrelated mutations.

Table 2.

Single-Nucleotide Polymorphisms (SNP) Identified on Bartonella krasnovii Genomes from Different Experimental Setups, in Comparison to Their Original Clone Genomic Data.

Sample Description	Etest MIC (µg/ml)	Chromosome Position	Type	Reference	Alternative	Evidence	Gene	Strand	rpoB Position	Effect
Resistance clone from RD	1	613444	SNP	G	A	A:165 G:0	rpoB	+	1588	D to N (530)
Resistance clone from RD	0.75	613444	SNP	G	A	A:352 G:0	rpoB	+	1588	D to N (530)
Resistance clone from RD	>32	613474	SNP	C	T	T:125 C:0	rpoB	+	1618	H to Y (540)
Resistance clone from RD	>32	613484	SNP	G	A	A:126 G:0	rpoB	+	1628	R to H (543)
Resistance clone from RD	>32	613484	SNP	G	A	A:89 G:0	rpoB	+	1628	R to H (543)
Resistance clone from RD	>32	613490	SNP	C	T	T:174 C:0	rpoB	+	1634	S to F (545)
Resistance clone from RD	0.064–0.094	613498	SNP	G	T	T:304 G:0	rpoB	+	1642	G to C (548)
Tolerant clone from 10X-MIC (5 days)	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Tolerant clone from 10X-MIC (5 days)	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Tolerant clone from 10X-MIC (5 days)	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Tolerant clone from 10X-MIC (5 days)	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Tolerant clone from 10X-MIC (5 days)	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Tolerant clone from 10X-MIC (11 days)	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Tolerant clone from 10X-MIC (11 days)	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Heteroresistance clone from 1X-MIC	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Heteroresistance clone from 1X-MIC	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Heteroresistance clone from 1X-MIC	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Heteroresistance clone from 2X-MIC	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND
Heteroresistance clone from 2X-MIC	0.023	ND	ND	ND	ND	ND	ND	ND	ND	ND

Note.—RD, RIF diffusion test plate; ND, not detected.

A Mathematical Model Outlines Hypothetical Conditions Required to Reproduce the Emergence of Resistant Colonies during RIF Exposure

Using a mathematical model (fig. 4), we estimated the number of resistant colonies that could originate from a susceptible population during drug exposure, as observed in the RIF diffusion tests. The model assumes logistic growth and its performance were compared with an experimentally controlled procedure (see Materials and Methods). We tested different model settings to identify the potential scenarios that could replicate the in vitro results. The experimental controlled procedure resulted in the emergence of 1–3 confirmed resistant colonies per inhibition zone of 95% (19/20) on the tested RD plates, whereas no colonies were obtained from 20 plates containing RIF at 0.1 µg/ml seeded in parallel. Approximately 1.3 ± 0.03×103 bacteria (N0 exp = 1,300) were seeded within the edge of the inhibition zone of the RD plates (area composed of 0.5 cm from the growth lawn to the disk). In the model, we first tested different mutation rates (µ) and phenotypic switching rates from sensitive to heteroresistant (v), maintaining the remaining parameters under the experimentally estimated values. Results indicated that at the normal µ and v values the model could not replicate the empirical results (fig. 4). Instead, resistant (true mutants) and heteroresistant (not mutants) colonies were observed in the model only with higher v and µ values (fig. 4). Particularly, true mutant colonies emerged in the model when µ was 1,000-fold higher and v was 10- or 50-fold higher, without involving an accumulation of heteroresistant colonies (fig. 4). Secondly, as a 24 h sub-MIC window was identified by the LC-MS/MS assay, we tested if this period could have contributed to the emergence of mutants on the RD plates by potentially increasing the population size. Accordingly, we evaluated the model with an altered inoculum size (N0 alt = 20,800) representing the increment of N0 if no inhibition was suffered during the sub-MIC period (see Materials and Methods). We found that even with a larger population size, the model could not replicate the empirical results without significantly increasing the mutation rates (fig. 4). In fact, resistant colonies were only obtained when the mutation rates from the heteroresistant subpopulation (µH) was ≥100-fold higher than the normal rates. Notably, under these hypothetical increased population size and mutation rates, altered switching rates were not required. Third, we tested if resistant mutations could be generated during the sub-MIC period, potentially allowing selection during supra-MICs. Starting with the N0 exp, we recorded the emergence of mutated individual cells during the first 24h. Results showed that mutations occurred only when µS are ≥100-fold higher than the normal mutation rates and the growth rate of sensitive cells is maintained at 100% capacity (no inhibition) (fig. 4). This demonstrates that in order to replicate the empirical results, B. krasnovii should undergo a hypermutator state, together with an increased capacity to replicate during RIF exposure (higher heteroresistance capacity) or a maximum capacity to multiply during the sub-MIC period (first 24h of incubation).

Fig. 4.

A mathematical model for microbial growth and adaptation in the presence of a drug due to heteroresistence and resistance. (A) Illustration of the ordinary differential equation system of the model, logistic growth is assumed but not shown. The model was run at experimentally defined parameters (i.e., normal rates) versus hypothetical rates representing increments of µ, mutation, and v, sensitive-to-heteroresistance switching rates, and inoculum size (all N = 20 runs per group). (B) Resulting number of hypothetical resistant (top panel) and heteroresistant colonies (bottom panel) obtained in the model (altering µ and/or and v) versus the experimental data. Model parameters that best fitted the experimental results (purple) are highlighted in yellow. (C) Resulting number of hypothetical resistant colonies obtained by increasing the initial population size (N0 alt) and the heteroresistant mutation rates (µH). (D) Resulting number of mutated cells during a 24 h period obtained per parameter set, altering heteroresistant, µH, or general population mutation rates, µS, and allowing maximum growth rate of sensitive cells during incubation (rS=rR).

The Emergence of Resistant Colonies Is Inhibited at Microaerophilic Conditions

Since it has been reported that nonlethal concentrations of antibiotics could lead to elevated mutagenesis in bacteria, in association with oxidative stress and the generation of reactive oxygen species (ROS) (Kohanski et al. 2010; Sebastian et al. 2017), we evaluated the emergence of RIF mutants on RD plates incubated at different oxygen levels, namely capnophilic (17% O2 with 5% CO2) and microaerophilic (5% O2 with 10% CO2). B. krasnovii bacteria grow optimally under both conditions, showing efficiency of plating (EOP) of 99.1% (microaerophilic/capnophilic; data from four independent biological replicates). Results from 60 RD plates, prepared in parallel and incubated either at capnophilic or microaerophilic conditions (30 paired replicates), showed that the number of colonies within the inhibition zone was significantly higher under the capnophilic conditions (t test, P = 5.2×10−8; fig. 5). Notably, the MIC level was not found to be significantly different between these two conditions (supplementary fig. 9A, Supplementary Material online; Wilcoxon test, P = 0.734). Notably, the level of heteroresistance was augmented under microaerophilic conditions in comparison to capnophilic conditions (supplementary fig. 9B, Supplementary Material online). These results prompted the hypothesis that the emergence of resistant colonies is associated with oxidative stress during incubation but not via an increased heteroresistance capacity alone.

Fig. 5.

The emergence of resistant colonies within the inhibition zone of RD is inhibited at microaerophilic conditions. (A) Comparison of parallel cultures seeded on RD plates and incubated at either capnophilic or microaerophilic conditions (N = 30 per group). Right panel illustrating representative RD plates incubated at the different oxygen atmospheres (yellow arrows highlighted the emergence of CFUs on the inhibition zone of capnophilic incubated plates). ****P < 0.00001. (B) Mutation spectra obtained from colonies isolated from the inhibition zone of RD plates versus colonies recovered from typical fluctuation test RIF plates. χ2 test for homogeneity is indicated.

To further explore the oxidative stress hypothesis, we compared the mutation type frequencies between resistant colonies that emerged from the inhibitory zone of RD plates (postantibiotic resistance development) and pre-existing resistant mutants selected through fluctuation tests (preantibiotic exposure development). If oxidative damage is playing a role on the emergence of colonies within the inhibition zone of RD plates, a difference in the frequency of oxidative mutations would be expected, such as C-T (ROS induced cytosine deamination) and G-T (by oxo-8-G) mutations (Kreutzer and Essigmann 1998; Foti et al. 2012; Degtyareva et al. 2013). As shown in figure 5, the mutation spectra obtained from 51 RD colonies and 74 fluctuation test colonies revealed significant differences (χ2 for homogeneity, χ2 = 15.323, P = 0.012; with simulated data of 2000 replicates). Notably, the oxidative mutations C-T and G-T were higher from colonies of RD plates. Thus, these results support a role for oxidative stress in the onset of the rpoB mutations.

The Bartonella Response to RIF Stress Is Characterized by a Mutator State and the Upregulation of All Prophages and Plasmid Genes

To study the effect of RIF on the gene expression of B. krasnovii, we performed RNA sequencing (RNA-seq) analysis from cultures incubated under capnophilic conditions for 18 hours (18h) in the presence or absence of RIF. For accurate determination of the expression changes, four replicates of time-zero cultures (starting cultures) were included and compared with four replicates of 18h RIF exposure (at 0.02 µg/ml) and four replicates of 18h antibiotic-free cultures. Out of the 1,826 protein coding genes annotated for the B. krasnovii genome (NCBI accession nos. CP031844.2 and CP042965.1), RNA products from 1,760 genes were identified, including 1,759 protein coding genes and one ncRNA (6S gene). Both 18h treatments showed a significant differential expression compared with the time zero cultures (fig. 6), but the RIF treatment elicited major expression changes, with 71.4% of the genes differentially expressed, whereas only 29.6% were found significantly different between the 18h antibiotic-free and the time zero cultures (fig. 6supplementary table 8, Supplementary Material online). Remarkably, the comparison between 18h antibiotic-free and 18-h RIF cultures showed 69.4% of the genes differentially expressed (supplementary table 8, Supplementary Material online). Regarding KEGG pathways, ribosome and metabolic processes were significantly downregulated to a large degree on 18h RIF in comparison to the controls (supplementary table 9, Supplementary Material online).

Fig. 6.

Differential RNA expression of Bartonella krasnovii in response to RIF stress. (A) Heatmap of the RNAseq results from the seeding cultures (0h) and the cultures incubated for 18h with or without RIF (0.02 µg/ml). (B) Volcano plots of the pairwise comparison between the treatments.

The DE analysis detected four major changes upon exposure to RIF: 1) deregulation of key DNA repair genes, including significant downregulation of mutS, mutY, and rmuC genes (with ≥3.2 fold difference, FDR < 2.5×10−12), and the upregulation of the mfd, recO, mutM, and abbB genes (with >2.4 fold difference, FDR < 3.2×10−21; fig. 7); 2) upregulation of 90% (104/121) of all prophage genes, representing all the predicted prophage loci (fig. 7); 3) upregulation of all plasmid genes (FDR < 1.3×10−9; fig. 7); and 4) deregulation of transcriptional regulators, including the upregulation of rpoE (2.0-folds difference, FDR < 4.4×10−6) and downregulation of rpoH and dksA genes (>3.3-folds difference, FDR < 3.9×10−19; fig. 7). Additionally, changes on the expression of many virulence factor genes were recorded, including the upregulation of badA, filamentous hemagglutinin, toxin–antitoxin genes, as well as downregulation of virB and hemin binding protein genes supplementary fig. 10, Supplementary Material online). Furthermore, the upregulation of lexA gene (fold change = 1.7, FDR = 0.0003), and the downregulation of recA (fold change = 0.42, FDR = 5.3×10−29, fig. 7) during RIF exposure, confirmed that RIF did not elicit an SOS response in B. krasnovii. All these results were supported by the comparison of 18h RIF versus 0- and 18h free versus 0 h (supplementary fig. 11, Supplementary Material online).

Fig. 7.

Major differences in the RNA expression of Bartonella in response to RIF stress. Fold change (log2) of the RNA expression between 18h RIF versus 18h free cultures of: (A) Repair genes; (B) Phage genes (prophage nomenclature based on phaster.ca); (C) Plasmid genes; (D) Transcriptional regulators. Red dash lines indicate the threshold of FC of 1.5 applied. Abbreviations: GTA, gene transfer agent loci; ROR, run-off replication loci; BER, base excision repair; DSBR, double-strand break repair; HR, homologous recombination; MMR, DNA mismatch repair; NER, nucleotide excision repair; T4SS, type IV secretion system. Error bars indicated SD from the mean.

The exposure to RIF did not cause a clear trend in the expression of the DNA polymerases (supplementary fig. 12, Supplementary Material online). Moreover, except for rpoZ, all the RNA polymerase subunits were downregulated during both 18h cultures (with and without RIF), including the RIF gene target, rpoB (supplementary fig. 12, Supplementary Material online).

We also compared the expression of 18h cultures with and without RIF under microaerophilic conditions. Remarkably, the expression under these conditions mirrored the capnophilic incubation (supplementary fig. 13, Supplementary Material online). Indeed, only 11 genes were found significantly differently expressed between of capnophilic and microaerophilic cultures exposed to RIF.

Discussion

The term adaptive mutation was originally intended to describe mutations formed in response to a stressful environment (Rosenberg 2001). This concept challenged the notion that mutations arise spontaneously and randomly during replication, as demonstrated by Luria and Delbruck’s pioneering study (Luria and Delbruck 1943). Thus, adaptive mutation is a process that is thought to involve the formation of mutations in nonactively dividing cells, namely replication-independent mutations. Cairns et al. (1988) showed that such a scenario is possible when cells are exposed to nonlethal stress conditions. Although these studies were controversial, because of a lack of distinction between stress conditions and selection for mutations, follow up studies, using different conditions for stress and selection were consistent with the hypothesis that cells under stress can generate mutations (Rosenberg 2001). Later, the discovery of error-prone polymerases expressed during stress (Sommer et al. 1993; Foster 2005) shed light on the potential mechanisms promoting transient mutagenic states. Here, we present a phenomenon in which bacteria undergo adaptive mutations despite the presence of an antibiotic at lethal concentrations. In a series of microbial experiments (figs. 1–3) supported by a mathematical model (fig. 4), we found that Bartonella spp. exposed to gradually increasing antibiotic concentrations, elicit a mutagenic state that yields resistant mutations, even in the absence of known error-prone polymerases and the conventional SOS response. Furthermore, this occurred within a population having a limited degree of heteroresistance (fig. 3).

One of the most important and interesting characteristics that we identified was the extraordinary antibiotic dose-dependence of the phenomenon: Direct exposure of bacteria to supra-MIC RIF plates prevented the emergence of de novo mutants. Our data suggest that gradual exposure to RIF concentrations, which includes a 24h sublethal period before supra-MICs are reached, is critical for the acquisition of resistant mutations (fig. 1). Since B. krasnovii are slow-growing bacteria, able to produce a >1-mm CFU in antibiotic-free media in ≥4 days (Gutiérrez et al. 2018, 2020), the first 24h may represent a period of adaptation and response to the emerging stress, but under a limited growing capacity. Accordingly, this period could allow the stimulation of a mutagenic response in tolerant/surviving cells, leading to the acquisition of resistant mutations at supra-MIC concentrations. Remarkably, RIF exposure at MIC levels during 18h elicited the upregulation of rpoE sigma factor, whereas rpoD and rpoH were downregulated (fig. 7). In the absence of the rpoS general stress regulator in Bartonella, RpoE has been implicated as the main sigma factor during environmental-related Bartonella stresses, such as exposure to arthropod-associated temperature and hemin concentrations (Abromaitis and Koehler 2013; Abromaitis et al. 2013). In E. coli, RpoE has also been associated with stress response mutagenesis (Gibson et al. 2010). Additionally, expression and experimental data demonstrated that the emergence of resistant colonies is independent of the SOS-response, as recA was downregulated and lexA was upregulated during RIF exposure, and UV exposure did not promote the emergence of resistance in RD plates (fig. 7 and supplementary fig. 3, Supplementary Material online). Hence, our results suggest that rpoE drives the SOS-independent Bartonella response to gradual RIF stress, like it does for cold shock and hemin toxicity responses.

The fact that emerging resistant mutants accumulated a single mutation only at the rpoB gene and no other mutation in the whole genome is remarkable (table 2). Notably, as RIF stress was maintained at supra-MICs from the 24th hour to the end of the experiment, the potential appearance of additional mutations may represent a fitness cost for the bacteria, promoting their loss during the experiment (negative epistasis) (Wong 2017). These results led to the question: how is it possible that from a small number of bacterial cells seeded on the RIF plates, an rpoB mutation is generated and selected before the bacteria is killed by the antibiotic? Using a mathematical model, we showed that to reproduce this phenomenon, the mutation rates and heteroresistance capacity should be elevated by several orders of magnitude compared with the normal estimated values (fig. 4). Interestingly, Krasovec et al. (2014) showed that mutation rates to RIF can be plastic and inversely proportional to cell density. We observed a similar relationship between the number of bacteria seeded per inhibition zone and the mutation rate (supplementary fig. 2C, Supplementary Material online). In line with these results, the Bartonella expression response demonstrated a hypermutator state during RIF treatment (discussed in detail below). On the other hand, experimental data demonstrated a minimal heteroresistant capacity of B. krasnovii, with a population multiplication level of 0.2 cell-division/day, active heteroresistant subpopulations of 0.0002%, and no evidence of rpoB gene-amplification or additional mutations on heteroresistant clones (fig. 2, table 2, and supplementary fig. 6, Supplementary Material online). Therefore, these results suggest that the model-predicted increment in heteroresistance capacity must occur in a rare stressed subpopulation of cells, undetected by our methods, during the gradual accumulation of RIF. Interestingly, if the sub-MIC period of the RD plates could allow a few cell divisions on the survival cells, increasing the effective population size, the model showed that the occurrence of resistant colonies would be exclusively dependent on the hypermutator state, without the need of increasing the heteoresistance capacity. In summary, the model reveals three scenarios that could replicate the experimental data: 1) synergy of increased mutation rates and heteroresistant capacity; 2) increased population size during the first 24h together with increased heteroresistant mutation rates; and 3) maximum cell multiplication capacity during sub-MICs with increased mutation rates of the entire population. Nevertheless, since PAP experiments showed inhibition under sub-MICs (ranging from 3% to up to 35%; fig. 2), and resistant colonies emerged with a delayed growth and even in “aged” RD plates, the second and third scenarios become highly unlikely. Therefore, our model results advocate for an integrative and transient hypermutator and hyperheteroresistant states during gradual and supra-MIC.

Three major factors are proposed as the potential enhancers for the generation of mutations in B. krasnovii during RIF exposure. First, RIF provokes a mutator state, which includes the downregulation of mutS, a mismatch repair gene and mutY, a DNA oxidative-damage repair gene (fig. 7). Consistent with this, Long et al. (2016) found that in E. coli mutS and mutY knockouts, sublethal antibiotic exposure (to norfloxacin) induced an increased mutation rate independent of error-prone polymerases. This suggests a role for the high-fidelity polymerases, Pol I and III, in the generation of mutations. Additionally, mfd, a transcription-coupled repair gene, was upregulated during RIF exposure. mfd has been associated with elevated mutagenesis in Bacillus subtilis in nongrowing stressed cells (Ross et al. 2006), particularly in MutY-deficient backgrounds (Leyva-Sanchez et al. 2020). Second, RIF induces the expression of 90% of all phage genes encoded on the B. krasnovii genome, representing all predicted prophages (fig. 7). In other bacteria, prophages are induced after exposure to sub-MIC antibiotic concentrations (Goerke et al. 2006; Long et al. 2016; Molina-Mora et al. 2020). These results suggest that stressed cells could undergo DNA breaks (single- or double-strand breaks) due to the excision of induced phages (Craig and Nash 1983; Ramisetty and Sudhakari 2019), all of this under a limited DNA repair capacity (fig. 7). Third, high oxygen levels were significantly associated with the appearance of mutations in RIF-exposed cells, suggesting that oxidative damage plays a role in the emergence of resistant colonies (fig. 5). Support for this idea comes from the observation that RIF induces ROS formation in bacteria that were exposed to the drug in liquid cultures (Kohanski et al. 2010; Piccaro et al. 2014; Sebastian et al. 2017), and it was proposed that on agar surfaces, DNA damage by oxidative agents also has a large impact (Foster 2007). Therefore, we propose that the parallel onset of an error-prone response, DNA breaks due to prophage inductions, and oxidative damage during RIF exposure, are sources of the mutagenic state described here.

In this study, the emergence of resistant mutants within the inhibition zone of antibiotic diffusion tests was only recorded with specific antibiotics (RIF, streptomycin, and mitomycin C; supplementary fig. 1, Supplementary Material online). Interestingly, only with mitomycin C, a known mutagenic drug (Tsukamura and Tsukamura 1962; Maccubbin et al. 1997), could quantitatively replicate the RIF phenomenon. Similarly, several researchers have recorded a high frequency of rpoB mutations conferring RIF resistance even in the absence of RIF selection (Taddei et al. 1995, 1997; Bjedov et al. 2003; Foster 2007; Wrande et al. 2008; Katz and Hershberg 2013). Among these studies, the “resting organisms in a structured environment” (ROSE) phenomenon showed an elevation of RIF-resistant mutation frequencies in a manner dependent on SOS-response but independent of error-prone polymerases (Taddei et al. 1995). Later, it was demonstrated that rpoB mutations confer an adaptative advantage in aged colony microenvironments, independent of a mutator state (Katz and Hershberg 2013; Hershberg 2017). In contrast to those studies, here the mutations were generated during RIF exposure in a SOS-independent manner. The fact that other resistant-associated loci, associated with point mutations, such as with quinolone resistance determinants, do not result in a similar outcome as for RIF, remains to be elucidated, but raises the query “which factors confer a mutational hotspot to rpoB in Bartonella?”

In a clinical setting, the pharmacokinetics of antibiotics is expected to mirror the diffusion dynamics before reaching the infection site (Levison and Levison 2009), raising the possibility that the studied phenomenon may also occur in vivo. In fact, treatment failure and relapses have been reported during treatment of Bartonella infections (Koehler and Tappero 1993; Maurin and Raoult 1993; Rolain et al. 2004). Hence, the demonstrated capacity of Bartonella to sense and respond to the gradual accumulation of antibiotic stress, promotes the exploration of this phenomenon and its implication in a clinical setup.

In conclusion, this study uncovers a novel and alternative mutagenic state expressed in a slow-growing pathogenic bacterium under antibiotic stress. Our findings demonstrate that the gradual accumulation of antibiotic stress induces a response in surviving cells that promotes the downregulation of key DNA repair genes in nondividing cells in parallel with an overall prophage induction and oxidative stress. These changes provide an alternative drug resistance adaptation mechanism unveiled for a slow-growing and pathogenic bacterial genus. We raise the possibility that these findings may have clinical implications, as similar pharmacokinetics of antibiotics may also occur during antibiotic therapy.

Materials and Methods

Bartonella Strains

Bartonella species and strains included in this study are: B. krasnovii and B. kosoyi type strains (Gutiérrez et al. 2020), B. elizabethae (ATCC 49927), B. grahamii (ATCC 700132), and B. henselae (strain 54 A) (Gutiérrez et al. 2013).

Media and Incubation Conditions

Chocolate agar media from two providers (Novamed, Jerusalem, IS and Hylabs, Rehovot, IS), blood agar (Hylabs, Rehovot, IS) and chocolate agar media containing RIF at different concentrations (Hylabs and in house preparation) were used in this study. All incubations for solid media were performed at 37 °C with 5% CO2 in Water-Jacketed CO2 incubators (ThermoFisher Scientific, MA) or under microaerophilic conditions, generated by CampyGen atmospheric generation systems (Oxoid, ThermoFisher Scientific, MA) in anaerobic jars. For liquid cultures, Schneider’s insect medium (Sigma-Aldric, MO) supplemented with 10% of fetal bovine serum (Biological Industries, Beit HaEmek, IS) and 5% sucrose, sterilized by filtration (0.2 µm) was used. Liquid culture incubations were performed at 37 °C, in sealed anaerobic jars under humid atmosphere and constant shaking (150 rpm), with or without CampyGen bags.

Antibiotic Diffusion Tests

Two types of diffusion tests were employed: 1) Etest (BioMériux, Marcy-l’Étoile, FR); and 2) in house disk diffusion tests. For the latter, Whatman filter papers (No. 2) were cut in circles with a hole-puncher (6.0 ± 0.5 mm diameter) of and sterilized in an autoclave, following previous described methods (Gefen et al. 2017). Then, the antibiotic was added to the paper disks, at the desired concentration, and was placed on the culture lawns. All antibiotics were bought from Sigma-Aldrich (Rehovot, IS), and prepared and maintained according to previously described methods (Andrews 2001). Culture lawns were prepared from fresh bacterial solutions in sterile PBS homogenized with sterile glass beads of 1 mm. All antibiotic diffusion plates employed were incubated for at least 10 days (unless indicated otherwise).

Isolation of Colonies within the Inhibition Zone of Antibiotic Diffusion Plates

Colonies that emerged within the inhibition zone of antibiotic diffusion plates were picked with a bacteriological loop and seeded on fresh chocolate agar media without antibiotic. After isolated colonies were obtained from the secondary plate, five colonies were suspended on sterile PBS with a swab and plated in a new antibiotic diffusion test to confirm the resistant phenotype. In parallel, a bacterial suspension was prepared on LB broth supplemented with 20% of glycerol and stored at −80 °C.

Agar Dilution Method

The confirmation of the RIF MIC was performed through agar dilution assay, following methods described elsewhere (Andrews 2001). In brief, fresh chocolate agar media was prepared, and RIF was added to the media to final concentrations representing 2-fold increments (0.01–512 µg/ml) before purring into sterile plastic petri plates. Then, 2 µl drops of a fresh culture suspended in PBS at a concentration of McFarland 1, were deposited in triplicates into the plates. Plate controls (without RIF) were run in parallel. Plates were checked every 24h to identify growth at the different RIF concentrations.

DNA Extraction from Bacterial Cultures

Two DNA extraction protocols were performed, depending on the goal of the experiment. For conventional PCR assays, clonal bacterial suspensions DNA was extracted by a thermal protocol, as follows: 2-ml Eppendorf tubes were cap-sealed and incubated at 95 °C for 12 min, followed by a centrifugation at 8,500 rpm for 5 min at 4 °C. DNA was then collected from the supernatant and stored at −20 °C until further analyzed. For real-time qPCRs and whole-genome sequencing preparations, QIAGEN DNeasy blood and tissue kits (Hilden, DE) were used following the manufacturers protocols.

Conventional PCR and Sanger Sequencing

All PCR reactions were performed in 25 μl final volume of PCR-Ready High Specificity ready mix (Syntezza Bioscience Ltd, Jerusalem, IS) containing 1 μl of 10 μM solution of each primer, 21 μl ultrapure water (UPW; Biological Industries), and 2 μl of each extracted isolate-DNA. The amplification products were obtained by the following protocol: 5 min at 95 °C, followed by 35 cycles of 30 s at 95 °C, 30 s at 55 °C, 1 min at 72 °C, and a final step of 5 min at 72 °C. Finally, the PCR products were run on 1.0% agarose gel, stained with ethidium bromide, and visualized under UV light.

For the identification of RRDR of the rpoB, the following primers were employed: 1400 F (CGCATTGGCTTACTTCGTATG) and 2300 R (GTAGACTGATTAGAACGCTG) for an 800-bp rpoB fragment (Renesto et al. 2001). For identification of mutations conferring streptomycin resistance, a locus of the ribosomal protein S12 (rpsL) gene was amplified with rpsL-F (TACCCGATTGGGTGAGACGA) and rpsL-R (TGACGACGGGACATTGACTT) primers (this study), that amplify a 450-bp fragment of the rpsL gene.

All positive PCR products were purified and cleaned by NEB Exo-SAP PCR purification kit (New England Biolabs, Inc., MA) and subsequently sequenced with sense and antisense primers using BigDye Terminator v3.1 cycle sequencing chemistry (Applied Biosystems, CA) in an ABI 3730xl DNA Analyzer and the ABI’s Data collection and Sequence Analysis software. Sequence analyses were performed in MEGA version X (Kumar et al. 2018).

Quantification of Bacterial Inoculum by CFU Counting

All initial and final bacterial loads were quantified by preparing 10-fold dilutions in sterile PBS and plating them in chocolate agar media. Accordingly, 5 μl drops (four to five replicates) per dilution were seeded sequentially in fresh plates. Plates were incubated and CFUs were counted after 48–72h of incubation under a stereoscope.

Pre-Existing RIF Mutants Test

A fresh Bartonella culture was extremely diluted (∼103 cells), to minimize the changes to carry rpoBR spontaneous mutants. Then, this inoculum was seeded on 83 parallel cultures (100 µl each) and shortly incubated at 37 °C ± 10% CO2 and 5% O2 (in Campygen atmospheric bag, Oxoid) overnight, to increase the bacterial load. Three parallel cultures were used to calculate the average final obtained number of bacteria per well by CFU counting. Eighty of the parallel cultures were seeded on the following plates (with 20 replicates per group): 1) chocolate agar plates with a 25 µg RIF diffusion disk (RD); 2) RIF-chocolate agar plates containing 0.04 µg/ml (∼2-fold higher the MICoriginal), 3) RIF-chocolate agar plates containing 0.08 µg/ml (∼4-fold higher the MICoriginal); and 4) RIF-chocolate agar plates containing 0.1 µg/ml (∼5-fold higher the MICoriginal). After 10 days of incubation, the number of colonies that emerged within the inhibition zone of the RD plates or colonies growing on the RIF-containing plates (0.04–0.1 µg/ml) were recorded.

Population Analysis Profile

PAP tests were performed following El-Halfawy and Valvano method (El-Halfawy and Valvano 2015). In short, chocolate agar plates containing 2-fold increments of RIF were prepared. Then, three heavy inocula (∼107 to 108 cells/ml) of fresh Bartonella cultures were prepared from isolated colonies and 10-fold dilutions (10−1 to 10−7) were seeded on the RIF chocolate agar plates (5–10 µl drops, in triplicates). Plates were incubated at 37 °C with 5% CO2 for up to 20 days. Once bacterial growth was identified, CFUs were counted under a stereoscope. Heteroresistance was defined according to El-Halfawy and Valvano criterion (El-Halfawy and Valvano 2015), as the presence of subpopulations that are able to multiply at antibiotic concentrations >8-fold higher than the maximum noninhibitory concentration for the general population.

qPCR Assays

All qPCR reactions were performed in StepOnePlus Real-Time PCR System (ThermoFisher Scientific) following methods and reagent recipes described earlier (Gutiérrez et al. 2018). For absolute quantification, a standard curve was prepared from 10-fold dilutions of a fresh Bartonella culture with loads determined by CFU plating. Primers qgltA-F2 (AGCCCATAAGGCGAAAAGGA) and qgltA-R2 (CATGGTGGTGCCAATGAAGC) (Gutiérrez et al. 2018), targeting a 117-bp fragment of the citrate synthase gene (gltA), were used.

For gene amplification assay, the relative copy number of rpoB gene relative to the gltA gene copy number (loci at 110,928 bp distance) was evaluated. For this, B. krasnovii cultures were seeded on five chocolate agar plates with RIF (0.04 µg/ml) and five chocolate agar plates without RIF. Once colonies were grown on both type plates, five colonies per plate were suspended in sterile PBS and DNA extracted. For qPCR amplification, a 109-bp rpoB fragment was targeted with qrpoB-F (GATCGGCCAGTAACGGTAGG) and qrpoB-R (GCTGCTGCGTAACAAGTGAG) primers (this study), and gltA was targeted with the qgltA-F2 and qgltA-R2, described above. Primer efficiencies were evaluated with a standard curve of known CFU values. Finally, delta-delta Cq analysis was performed (Livak and Schmittgen 2001). As a confirmatory assay, the rpoB Cq values were compared between five single 1-mm CFUs (one CFU/sample as normalization) excised from RIF-free, 0.02 and 0.04 µg/ml RIF plates.

Dye-Dilution Proliferation Assay

A dye-dilution flow cytometry assay was developed to evaluate cell divisions in Bartonella cells after incubation in liquid media containing RIF, following protocols described elsewhere (Ueckert et al. 1997; Flannagan and Heinrichs 2018). Following manufacturers’ protocols, we first confirmed that the cell proliferation dye eFluor450 (eBioscience Cell Proliferation Dye eFluor450, Invitrogen, CA) could serve as a marker of bacterial division prior antibiotic exposure. Accordingly, the use of this dye was confirmed under fluorescence microscopy (supplementary fig. 5A, Supplementary Material online). Then, as a preliminary test, labelled B. krasnovii cells (sensitive and resistant) were incubated with RIF (0.2 µg/ml) and without RIF for 5 days, live and dead cells were marked with a viability kit (LIVE/DEAD BacLight Bacterial Viability Kit Protocol, Molecular probes, Invitrogen) after incubation, and fluorescence levels of eFluor450 were measured on the surviving cells (gated by Syto9; supplementary fig. 5B, Supplementary Material online). Following standardization of the method, we prepared fresh 1-ml B. krasnovii liquid cultures (Ni = 3.7×105 cells) with different RIF concentrations (0.02, 0.04, and 0.1 µg/ml) previously labelled with eFluor450 and measured the fluorescence levels of eFluor450 on the surviving cells, after 5 days of incubation. Results were compared with a nonmultiplication control incubated in parallel at 4 °C without RIF (supplementary fig. 5C, Supplementary Material online).

Death Curves

Survival rates of B. krasnovii were measured through death curves on 10-ml liquid media containing 0.2 µg/ml RIF (10-fold higher the MIC) by plating 100-µl aliquots on chocolate agar plates every 24 h, following Balaban et al. (2019) recommendations. CFUs were counted from the plates after 3–4 days of incubation.

Whole-Genome Sequencing and Analysis

Whole genome sequencing was performed for the following isolates: 1) seven resistant colonies isolated from the inhibition zone of RIF diffusion plates; 2) seven tolerant clones recovered after 5 and 11 days from liquid cultures exposed to RIF at 10-folds the MIC levels; 3) five confirmed heteroresistant clones isolated from PAP experimental plates (from 0.02 and 0.04 µg/ml RIF plates); 4) four original clones, obtained after at least three population bottlenecks (to ensure clonality) as the starting cultures used on these experiments. Accordingly, all clones were grown on chocolate agar media, DNA was extracted, and the MIC was confirmed by Etest. Nextera libraries (Illumina, CA) were prepared and Illumina MiSeq 150-bp paired-end sequencing was performed in Genomics Applications Laboratory of The Hebrew University of Jerusalem, Israel. In average, 1.78×106 sequences/read was obtained, giving an average genome coverage of 243×.

Sequencing reads were processed in the following manner: quality was checked using fastqc (Andrews 2010). Adapters and bad quality reads (q < 30 and length <100) were filtered out with TrimGalore (software code available at: https://github.com/FelixKrueger/TrimGalore; last accessed March 29, 2021). SNP variant calling analysis was performed on reads mapped to the wild-type B. krasnovii genome (NCBI accession nos. CP031844.2 and CP042965.1) using SNIPPY (Seemann 2015). Then, all SNPs present in the original clones or with a coverage <100 bp were filtered out from the analysis.

RIF Concentration Determination from Agar Plugs by LC-MS/MS Assay

RIF concentration was measured on agar plugs excised at a fixed distance from the disk (at 2.0 cm) from chocolate agar plates incubated for 12h to 10 days, using organic extraction and LC-MS/MS. In brief, 33 sterile chocolate agar plates with RD (25 µg) were prepared and incubated at 37 °C with 5% CO2. At every time point, nine cylindrical agar plugs of 0.4 cm×0.7 cm, from three chocolate agar plates, were cut, weighted, and processed. The organic extraction was performed with 400 µl of acetronitrile (HPLC grade), sonicated for 20 min. Then, samples were centrifuged (13,500 rpm for 10 min) and the supernatant was dried at 60 °C for 10 min. Finally, the pellet was resuspended in 300 µl of a solution of acetonitrile 1:1 double distilled water, recentrifuged, and the supernatant analyzed by LC-MS/MS. A standard curve of RIF was run in parallel. As negative controls (triplicates), incubated RIF-free plates were run in parallel, which resulted in no detection of RIF. For standardization of the extraction method, nine agar plugs from RIF-free plates were inoculated with 10 µl of 5 µg/ml (three plugs), 10 µl with 0.5 µg/ml (three plugs), and 10 µl with 0.05 µg/ml. This method resulted in good correlation (R2 = 0.9964; supplementary fig. 14, Supplementary Material online) between expected and observed measurements.

RIF Accumulation Test during Incubation

To confirm the accumulation of RIF on the plates, we performed the following experiment. Twenty-five RD plates were incubated sterilely for different periods (1, 2, 3, 5, and 10 days). At every time point, five plates were recovered from the incubator and inoculated with Bartonella fresh cultures and placed back to the incubator. In all time points, controls of fresh RD media were prepared in parallel with the same B. krasnovii inoculum. After reincubation with the bacteria, and once the culture lawn was formed, the inhibition zone radius of each plate inhibition zone was measured by picturing the plates and analyzing them in ImageJ (Schneider et al. 2012).

Mathematical Model

We developed a stochastic mathematical model to estimate the number of resistant colonies that could originate from a sensitive population during drug exposure, as the outcome of a phenotype switch from sensitive phenotype to heteroresistant phenotype followed by a mutation to resistant phenotype, or a mutation from sensitive phenotype to resistant phenotypes (fig. 4). The model assumes logistic growth and considers the following variables (eq. 1): 1) number of inoculated sensitive bacteria, N0; 2) the growth rates of sensitive, rS, heteroresistant, rH, and resistant, rR, bacteria during drug exposure; 3) the rate of phenotype switching from sensitive to heteroresistant, v; 4) the mutation rate of the sensitive, µs, and heteroresistant, µH, phenotypes; 5) the death rate of sensitive bacteria during drug exposure, δ; and 6) the length of the experiment. The stochastic model was implemented in Python (Van Rossum and Drake 2009) using the Gillespie algorithm and the τ-leap method (Gillespie 2001) and it can be summarized by the following differential equation system,

The model was then compared with an experimental procedure prepared from an extremely diluted fresh culture as a starting point (148.3 ± 27.5 bacteria/ml), split it on 200 µl cultures and shortly incubated overnight (reaching ∼8.9×103 bacteria/well). Then, the entire volume of each parallel culture was spread on 20 RD plates. As controls for pre-existing mutants or mutants generated overnight, other 20 cultures were spread on agar plates containing 0.1 µg/ml (5-fold higher than MICoriginal). After incubation, the average inhibition zone was calculated from RD plates using ImageJ, resulting in 25.1 ± 1.2 cm2. As most resistant colonies emerged close to the edge area (∼0.5 cm from the growth lawn), we calculated the number of bacteria inoculated within this area. We estimated 1.3 ± 0.03×103 bacteria seeded within the edge of the inhibition zone. Colonies that emerge from the inhibition zone were counted, and one or two colonies per plate were confirmed as resistant by subsequent Etests and Sanger sequencing of the RRDR rpoB locus (supplementary table 10, Supplementary Material online). Following this, we examined the parameter values under which the model replicates the empirical results. We evaluated the model under the normal experimental estimated values: N0 exp = 1,300, the inoculum size estimated on the experimental procedure; rS= 0.008 cell division/hour, the growth rate of sensitive cells at supra-MICs, estimated by the flow cytometry dye-dilution method; δ = 0.04 death/hour; estimated by the death curve at 0.2 µg/ml; rH = 0.09, the growth rate of heteroresistant colonies from PAP tests; rR = 0.24, the growth rate of RIF resistant colonies on RIF agar plates; µ = 1.4×10−8 mutations/rpoB locus/generation, the estimated mutation rate obtained by RIF fluctuation tests; v = 2×10−6, the rate of phenotype switching from sensitive to heteroresistant phenotype, obtained by PAP tests; 10 days, the usual incubation limit of the plates. We also evaluated the model in three modes: 1) increasing both v and µ, from normal rates to 1,000-fold higher; 2) increasing the inoculum size, N0 alt = 20,800 sensitive cells, under the scenario that all inoculated cells could have multiplied at a normal rate during the first 24 h (N0 alt = N0 exp×2, being that n = 4 is the estimated number of generations per day in a Bartonella colony) (Gutiérrez et al. 2018) and different v and µ values (either for sensitive, µS, and heteroresistant, µH); and 3) screening for individual bacterial mutations during the first 24h, setting rS as normal rates in antibiotic-free media (rS= 0.24 cell division/hour) and different µS and µH values. Every parameter combination was evaluated 20 times (as the model is stochastic). Colonies obtained from each model run were classified as RH, RS, or H. True RH were obtained from a mutated H colony and were identified when the resistant cells reach at least 105 (for a visual CFU).

Characterization of Spontaneous rpoB Mutations through Fluctuation Tests

Classical fluctuation tests, using RIF as selective marker, were run following recommendations of standardized methods for Bartonella spp. (Gutiérrez et al. 2018). This method consists of seeding high inoculum cultures on plates containing homogenous supra-MIC RIF concentrations, selecting only pre-existing rpoB mutants (spontaneously generated during the RIF-free liquid culture). This assumption is supported by the fact that RIF-resistant colonies only appeared on plates containing homogenous RIF concentrations when high inocula (>107) are seeded, and as demonstrated in this study, plates containing homogeneous RIF concentrations obliged the studied phenomenon.

UV Mutagenesis Test

The effect of UV light on the emergence of RIF mutants within the inhibition zone of RD plates was evaluated as follows. First, UV death curve was performed to identify the minimal lethal UV dosage (MLUV, supplementary fig. 15, Supplementary Material online). Then, a fresh B. krasnovii culture was prepared, seeded on 30 RD plates, and treated in the following manner: 1) no UV exposure (controls, ten replicates); 2) exposed to the MLUV dosage at time zero (T0); and 3) exposed to the MLUV dosage after 72h postincubation (once the inhibition zone is formed). All plates were incubated in parallel at the conditions previously mentioned. After 10-day postincubation, plates were recovered and colonies within the inhibition zone were counted.

RNAseq Analysis

RNA was extracted from Bartonella liquid cultures following previously described methods (Quebatte et al. 2010; Omasits et al. 2013). Accordingly, liquid cultures were incubated at 37 °C under capnophilic and microaerophilic conditions with constant shaking (110 rpm), with or without RIF (0.02 µg/ml) and during 18h. An 18h period was chosen as a reasonable time for the population to respond to the RIF stress without causing significant death and 0.02 µg/ml as a concentration readily exposed in the first day on RD plates. Four replicates of the time-zero cultures (the seeding cultures for this experiment) were also processed to identify the RNA expression baseline. rRNA were depleted using NEBNext rRNA depletion kit (bacteria) and purified with RNAClean beads (NEB) quantified with Qubit RNA BR assay kit (Invitrogen) and quality checked (pre- and post-rRNA depletion) with TapeStation Agilent 4200 (Agilent Technologies, CA). The Library preparation was done using KAPA Stranded mRNA-Seq Kit Illumina Platforms. Sequencing was done on Nextseq 500 Illumina, in single read of 78 bp. Obtained RNA reads were analyzed through SPARTA pipeline (Johnson et al. 2016) to identify the differential expression (DE) between treatments. We further analyzed the DE by the Kyoto Encyclopedia of Genes and Genomes (KEGG; https:https://www.kegg.jp/kegg/kegg2.html; last accessed March 29, 2021), using KEGGREST conducted in R (R Development Core Team 2020; Tenenbaum and Maintainer 2020). We defined as significant different expressed genes, those genes with a fold change ≥1.5 (upregulated) or ≤0.67 (downregulated) and a FDR ≤ 0.01.

Data and Statistical Analyses

All analyses presented in this study were produced in R (R Development Core Team 2020), using the following packages: ggplot2, ggsignif, ggpubr, tidyverse, dplyr, forcats, scales, rstatix. Other analyses were performed in Microsoft Office Excel 2016. Illustrative figures were created using Biorender (http://BioRender.com; last accessed March 29, 2021).

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.