# Decomposing virulence to understand bacterial clearance in persistent infections | Panda Anku

### Fly population and maintenance

We used an outbred population of Drosophila melanogaster established from 160 Wolbachia-infected fertilised females collected in Azeitão, Portugal54, and given to us by Élio Sucena. For at least 13 generations prior to the start of the experiments the flies were maintained on standard sugar yeast agar medium (SYA medium: 970 ml water, 100 g brewer’s yeast, 50 g sugar, 15 g agar, 30 ml 10% Nipagin solution and 3 ml propionic acid; ref. 61), in a population cage containing at least 5000 flies, with non-overlapping generations of 15 days. They were maintained at 24.3 ± 0.2 °C, on a 12:12 h light-dark cycle, at 60–80 % relative humidity. The experimental flies were kept under the same conditions. No ethical approval or guidance is required for experiments with D. melanogaster.

### Bacterial species

We used the Gram positive Lactococcus lactis (gift from Brian Lazzaro), Gram negative Enterobacter cloacae subsp. dissolvens (hereafter called E. cloacae; German collection of microorganisms and cell cultures, DSMZ; type strain: DSM-16657), Providencia burhodogranariea strain B (gift from Brian Lazzaro, DSMZ; type strain: DSM-19968) and Pseudomonas entomophila (gift from Bruno Lemaitre). L. lactis43, Pr. burhodogranariea44 and Ps. entomophila45 were isolated from wild-collected D. melanogaster and can be considered as opportunistic pathogens. E. cloacae was isolated from a maize plant, but has been detected in the microbiota of D. melanogaster46. All bacterial species were stored in 34.4% glycerol at −80 °C and new cultures were grown freshly for each experimental replicate.

### Experimental design

For each bacterial species, flies were exposed to one of seven treatments: no injection (naïve), injection with Drosophila Ringer’s (injection control) or injection with one of five concentrations of bacteria ranging from 5 × 106 to 5 × 109 colony forming units (CFUs)/mL, corresponding to doses of approximately 92, 920, 1,840, 9200 and 92,000 CFUs per fly. The injections were done in a randomised block design by two people. Each bacterial species was tested in three independent experimental replicates. Per experimental replicate we treated 252 flies, giving a total of 756 flies per bacterium (including naïve and Ringer’s injection control flies). Per experimental replicate and treatment, 36 flies were checked daily for survival until all flies were dead. A sub-set of the dead flies were homogenised upon death to test whether the infection had been cleared before death or not. To evaluate bacterial load in living flies, per experimental replicate, four of the flies were homogenised per treatment, for each of nine time points: one, two, three, four, seven, 14, 21, 28- and 35-days post-injection.

### Infection assay

Bacterial preparation was performed as in Kutzer et al.24, except that we grew two overnight liquid cultures of bacteria per species, which were incubated overnight for approximately 15 h at 30 °C and 200 rpm. The overnight cultures were centrifuged at 2880 × g at 4 °C for 10 min and the supernatant removed. The bacteria were washed twice in 45 mL sterile Drosophila Ringer’s solution (182 mmol·L-1 KCl; 46 mol·L-1 NaCl; 3 mmol·L-1 CaCl2; 10 mmol·L-1 Tris·HCl; ref. 62) by centrifugation at 2880 × g at 4 °C for 10 min. The cultures from the two flasks were combined into a single bacterial solution and the optical density (OD) of 500 µL of the solution was measured in a Ultrospec 10 classic (Amersham) at 600 nm. The concentration of the solution was adjusted to that required for each injection dose, based on preliminary experiments where a range of ODs between 0.1 and 0.7 were serially diluted and plated to estimate the number of CFUs. Additionally, to confirm post hoc the concentration estimated by the OD, we serially diluted to 1:107 and plated the bacterial solution three times and counted the number of CFUs.

The experimental flies were reared at constant larval density for one generation prior to the start of the experiments. Grape juice agar plates (50 g agar, 600 mL red grape juice, 42 mL Nipagin [10% w/v solution] and 1.1 L water) were smeared with a thin layer of active yeast paste and placed inside the population cage for egg laying and removed 24 h later. The plates were incubated overnight then first instar larvae were collected and placed into plastic vials (95 × 25 mm) containing 7 ml of SYA medium. Each vial contained 100 larvae to maintain a constant density during development. One day after the start of adult eclosion, the flies were placed in fresh food vials in groups of five males and five females, after four days the females were randomly allocated to treatment groups and processed as described below.

Before injection, females were anesthetised with CO2 for a maximum of five minutes and injected in the lateral side of the thorax using a fine glass capillary (Ø 0.5 mm, Drummond), pulled to a fine tip with a Narishige PC-10, and then connected to a Nanoject II™ injector (Drummond). A volume of 18.4 nL of bacterial solution, or Drosophila Ringer’s solution as a control, was injected into each fly. Full controls, i.e., naïve flies, underwent the same procedure but without any injection. After being treated, flies were placed in groups of six into new vials containing SYA medium, and then transferred into new vials every 2–5 days. Maintaining flies in groups after infection is a standard method in experiments with D. melanogaster that examine survival and bacterial load (e.g. refs. 22, 63, 64). At the end of each experimental replicate, 50 µL of the aliquots of bacteria that had been used for injections were plated on LB agar to check for potential contamination. No bacteria grew from the Ringer’s solution and there was no evidence of contamination in any of the bacterial replicates. To confirm the concentration of the injected bacteria, serial dilutions were prepared and plated before and after the injections for each experimental replicate, and CFUs counted the following day.

### Bacterial load of living flies

Flies were randomly allocated to the day at which they would be homogenised. Prior to homogenisation, the flies were briefly anesthetised with CO2 and removed from their vial. Each individual was placed in a 1.5 mL microcentrifuge tube containing 100 µL of pre-chilled LB media and one stainless steel bead (Ø 3 mm, Retsch) on ice. The microcentrifuge tubes were placed in a holder that had previously been chilled in the fridge at 4 °C for at least 30 min to reduce further growth of the bacteria. The holders were placed in a Retsch Mill (MM300) and the flies homogenised at a frequency of 20 Hz for 45 s. Then, the tubes were centrifuged at 420 × g for one minute at 4 °C. After resuspending the solution, 80 µL of the homogenate from each fly was pipetted into a 96-well plate and then serially diluted 1:10 until 1:105. Per fly, three droplets of 5 μL of every dilution were plated onto LB agar. Our lower detection limit with this method was around seven colony-forming units per fly. We consider bacterial clearance by the host to be when no CFUs were visible in any of the droplets, although we note that clearance is indistinguishable from an infection that is below the detection limit. The plates were incubated at 28 °C and the numbers of CFUs were counted after ~20 h. Individual bacterial loads per fly were back calculated using the average of the three droplets from the lowest countable dilution in the plate, which was usually between 10 and 60 CFUs per droplet.

D. melanogaster microbiota does not easily grow under the above culturing conditions (e.g. ref. 42) Nonetheless we homogenised control flies (Ringer’s injected and naïve) as a control. We rarely retrieved foreign CFUs after homogenising Ringer’s injected or naïve flies (23 out of 642 cases, i.e., 3.6 %). We also rarely observed contamination in the bacteria-injected flies: except for homogenates from 27 out of 1223 flies (2.2 %), colony morphology and colour were always consistent with the injected bacteria (see methods of ref. 65). Twenty one of these 27 flies were excluded from further analyses given that the contamination made counts of the injected bacteria unreliable; the remaining six flies had only one or two foreign CFUs in the most concentrated homogenate dilution, therefore these flies were included in further analyses. For L. lactis (70 out of 321 flies), P. burhodogranaeria (7 out of 381 flies) and Ps. entomophila (1 out of 71 flies) there were too many CFUs to count at the highest dilution. For these cases, we denoted the flies as having the highest countable number of CFUs found in any fly for that bacterium and at the highest dilution23. This will lead to an underestimate of the bacterial load in these flies. Note that because the assay is destructive, bacterial loads were measured once per fly.

### Statistical analyses

Statistical analyses were performed with R version 4.2.166 in RStudio version 2022.2.3.49267. The following packages were used for visualising the data: “dplyr”68, “ggpubr”69, “gridExtra”70, “ggplot2”71, “plyr”72, “purr”73, “scales”74, “survival”75,76, “survminer”77, “tidyr”78 and “viridis”79, as well as Microsoft PowerPoint for Mac v16.60 and Inkscape for Mac v 1.0.2. Residuals diagnostics of the statistical models were carried out using “DHARMa”80, analysis of variance tables were produced using “car”81, and post-hoc tests were carried out with “emmeans”82. To include a factor as a random factor in a model it has been suggested that there should be more than five to six random-effect levels per random effect83, so that there are sufficient levels to base an estimate of the variance of the population of effects84. In our experimental designs, the low numbers of levels within the factors ‘experimental replicate’ (two to three levels) and ‘person’ (two levels), meant that we therefore fitted them as fixed, rather than random factors84. However, for the analysis of clearance (see below) we included species as a random effect because it was not possible to include it as a fixed effect because PPP is already a species-level predictor. Below we detail the statistical models that were run according to the questions posed. All statistical tests were two-sided.

#### Do the bacterial species differ in virulence?

To test whether the bacterial species differed in virulence, we performed a linear model with the natural log of the maximum hazard as the dependent variable and bacterial species as a factor. Post-hoc multiple comparisons were performed using “emmeans”82 and “magrittr”85, using the default Tukey adjustment for multiple comparisons. Effect sizes given as Cohen’s d, were also calculated using “emmeans”, using the sigma value of 0.4342, as estimated by the package. The hazard function in survival analyses gives the instantaneous failure rate, and the maximum hazard gives the hazard at the point at which this rate is highest. We extracted maximum hazard values from time of death data for each bacterial species/dose/experimental replicate. Each maximum hazard per species/dose/experimental replicate was estimated from an average of 33 flies (a few flies were lost whilst being moved between vials etc.). To extract maximum hazard values we defined a function that used the “muhaz” package86 to generate a smooth hazard function and then output the maximum hazard in a defined time window, as well as the time at which this maximum is reached. To assess the appropriate amount of smoothing, we tested and visualised results for four values (1, 2, 3 and 5) of the smoothing parameter, b, which was specified using bw.grid87. We present the results from b = 2, but all of the other values gave qualitatively similar results (see Supplementary Table 2). We used bw.method = “global” to allow a constant smoothing parameter across all times. The defined time window was zero to 20 days post injection. We removed one replicate (92 CFU for E. cloacae infection) because there was no mortality in the first 20 days and therefore the maximum hazard could not be estimated. This gave final sizes of n = 14 for E. cloacae and n = 15 for each of the other three species.

$${{{{{rm{Model}}}}}},1:,{{log }}left({{{{{rm{maximum}}}}}},{{{{{rm{hazard}}}}}}right), sim ,{{{{{rm{bacterial}}}}}},{{{{{rm{species}}}}}}$$

#### Are virulence differences due to variation in pathogen exploitation or PPP?

To test whether the bacterial species vary in PPP, we performed a linear model with the natural log of the maximum hazard as the dependent variable, bacterial species as a factor, and the natural log of infection intensity as a covariate. We also included the interaction between bacterial species and infection intensity: a significant interaction would indicate variation in the reaction norms, i.e., variation in PPP. The package “emmeans”82 was used to test which of the reaction norms differed significantly from each other. We extracted maximum hazard values from time of death data for each bacterial species/dose/experimental replicate as described in section “Do the bacterial species differ in virulence?”. We also calculated the maximum hazard for the Ringer’s control groups, which gives the maximum hazard in the absence of infection (the y-intercept). We present the results from b = 2, but all of the other values gave qualitatively similar results (see results). We wanted to infer the causal effect of bacterial load upon host survival (and not the reverse), therefore we reasoned that the bacterial load measures should derive from flies homogenised before the maximum hazard had been reached. For E. cloacae, L. lactis, and Pr. burhodogranariea, for all smoothing parameter values, the maximum hazard was reached after two days post injection, although for smoothing parameter value 1, there were four incidences where it was reached between 1.8- and 2-days post injection. Per species/dose/experimental replicate we therefore calculated the geometric mean of infection intensity combined for days 1 and 2 post injection. In order to include flies with zero load, we added one to all load values before calculating the geometric mean. Geometric mean calculation was done using the R packages “dplyr”68, “EnvStats”88, “plyr”72 and “psych”89. Each mean was calculated from the bacterial load of eight flies, except for four mean values for E. cloacae, which derived from four flies each.

For Ps. entomophila the maximum hazard was consistently reached at around day one post injection, meaning that bacterial sampling happened at around the time of the maximum hazard, and we therefore excluded this bacterial species from the analysis. We removed two replicates (Ringer’s and 92 CFU for E. cloacae infection) because there was no mortality in the first 20 days and therefore the maximum hazard could not be estimated. One replicate was removed because the maximum hazard occurred before day 1 for all b values (92,000 CFU for E. cloacae) and six replicates were removed because there were no bacterial load data available for day one (experimental replicate three of L. lactis). This gave final sample sizes of n = 15 for E. cloacae and n = 12 for L. lactis, and n = 18 for Pr. burhodogranariea.

$${{{{{rm{Model}}}}}},2 :,{{log }}({{{{{rm{maximum}}}}}},{{{{{rm{hazard}}}}}}), sim ,{{log }}({{{{{rm{geometric}}}}}},{{{{{rm{mean}}}}}},{{{{{rm{bacterial}}}}}},{{{{{rm{load}}}}}}),\ times ,{{{{{rm{bacterial}}}}}},{{{{{rm{species}}}}}}$$

To test whether there is variation in pathogen exploitation (infection intensity measured as bacterial load), we performed a linear model with the natural log of infection intensity as the dependent variable and bacterial species as a factor. Similar to the previous model, we used the geometric mean of infection intensity combined for days 1 and 2 post injection, for each bacterial species/dose/experimental replicate. The uninfected Ringer’s replicates were not included in this model. Post-hoc multiple comparisons were performed using “emmeans”, using the default Tukey adjustment for multiple comparisons. Effect sizes given as Cohen’s d, were also calculated using “emmeans”, using the sigma value of 2.327, as estimated by the package. Ps. entomophila was excluded for the reason given above. The sample sizes per bacterial species were: n = 13 for E. cloacae, n = 10 for L. lactis and n = 15 for Pr. burhodogranariea.

$${{{{{rm{Model}}}}}},3:,{{log }}({{{{{rm{geometric}}}}}},{{{{{rm{mean}}}}}},{{{{{rm{bacterial}}}}}},{{{{{rm{load}}}}}}), sim ,{{{{{rm{bacterial}}}}}},{{{{{rm{species}}}}}}$$

#### Are persistent infection loads dose-dependent?

We tested whether initial injection dose is a predictor of bacterial load at seven days post injection22,25. We removed all flies that had a bacterial load that was below the detection limit as they are not informative for this analysis. The response variable was natural log transformed bacterial load at seven days post-injection and the covariate was natural log transformed injection dose, except for P. burhodogranariea, where the response variable and the covariate were log-log transformed. Separate models were carried out for each bacterial species. Experimental replicate and person were fitted as fixed factors. By day seven none of the flies injected with 92,000 CFU of L. lactis were alive. The analysis was not possible for Ps. entomophila infected flies because all flies were dead by seven days post injection.

$${{{{{rm{Model}}}}}},4:,{{log }}({{{{{rm{day}}}}}},7,{{{{{rm{bacterial}}}}}},{{{{{rm{load}}}}}}), sim ,{{log }}({{{{{rm{injection}}}}}},{{{{{rm{dose}}}}}}),+,{{{{{rm{replicate}}}}}},+,{{{{{rm{person}}}}}}$$

#### Calculation of clearance indices

To facilitate the analyses of clearance we calculated clearance indices, which aggregate information about clearance into a single value for each bacterial species/dose/experimental replicate. All indices were based on the estimated proportion of cleared infections (defined as samples with a bacterial load that was below the detection limit) of the whole initial population. For this purpose, we first used data on bacterial load in living flies to calculate the daily proportion of cleared infections in live flies for the days that we sampled. Then we used the data on fly survival to calculate the daily proportion of flies that were still alive. By multiplying the daily proportion of cleared flies in living flies with the proportion of flies that were still alive, we obtained the proportion of cleared infections of the whole initial population – for each day on which bacterial load was measured. We then used these data to calculate two different clearance indices, which we used for different analyses. For each index we calculated the mean clearance across several days. Specifically, the first index was calculated across days three and four post injection (clearance index3,4), and the second index was calculated from days seven, 14 and 21 (clearance index7,14,21).

#### Do the bacterial species differ in clearance?

To test whether the bacterial species differed in clearance, we used clearance index3,4, which is the latest timeframe for which we could calculate this index for all four species: due to the high virulence of Ps. entomophila we were not able to assess bacterial load and thus clearance for later days. The distribution of clearance values did not conform to the assumptions of a linear model. We therefore used a Kruskal-Wallis test with pairwise Mann-Whitney-U post hoc tests. Note that the Kruskal-Wallis test uses a Chi-square distribution for approximating the H test statistic. To control for multiple testing we corrected the p-values of the post hoc tests using the method proposed by Benjamini and Hochberg90 that is implemented in the R function pairwise.wilcox.test.

$${{{{{rm{Model}}}}}},5:,{{{{{{rm{clearance}}}}}},{{{{{rm{index}}}}}}}_{3,4}, sim ,{{{{{rm{bacterial}}}}}},{{{{{rm{species}}}}}}$$

#### Do exploitation or PPP predict variation in clearance?

To assess whether exploitation or PPP predict variation in clearance we performed separate analyses for clearance index3,4 and clearance index7,14,21. As discussed above, this precluded analysing Ps. entomophila. For each of the two indices we fitted a linear mixed effects model with the clearance index as the response variable. As fixed effects predictors we used the replicate-specific geometric mean log bacterial load and the species-specific PPP. In addition, we included species as a random effect.

In our analysis we faced the challenge that many measured clearance values were at, or very close to zero. In addition, clearance values below zero do not make conceptual sense. To appropriately account for this issue, we used a logit link function (with Gaussian errors) in our model, which restricts the predicted clearance values to an interval between zero and one. Initial inspections of residuals indicated violations of the model assumption of homogenously distributed errors. To account for this problem, we included the log bacterial load and PPP as predictors of the error variance, which means that we used a model in which we relaxed the standard assumption of homogenous errors and account for heterogenous errors by fitting a function of how errors vary. For this purpose, we used the option dispformula when fitting the models with the function glmmTMB91.

$${{{{{rm{Model}}}}}},6 :,{{{{{{rm{clearance}}}}}},{{{{{rm{index}}}}}}}_{3,4},{{{{{rm{or}}}}}},{{{{{{rm{clearance}}}}}},{{{{{rm{index}}}}}}}_{7,14,21}, \ sim ,{{log }}({{{{{rm{geometric}}}}}},{{{{{rm{mean}}}}}},{{{{{rm{bacterial}}}}}},{{{{{rm{load}}}}}}),+,{{{{{rm{PPP}}}}}}+{{{{{{rm{bacterial}}}}}},{{{{{rm{species}}}}}}}_{{{{{{rm{random}}}}}}}$$

#### Does longer-term clearance depend upon the injection dose?

In contrast to the analyses described above, we additionally aimed to assess the long-term dynamics of clearance based on the infection status of dead flies collected between 14 and 35 days and 56 to 78 days after injection. Using binomial logistic regressions, we tested whether initial injection dose affected the propensity for flies to clear an infection with E. cloacae or Pr. burhodogranariea before they died. The response variable was binary whereby 0 denoted that no CFUs grew from the homogenate and 1 denoted that CFUs did grow from the homogenate. Log-log transformed injection dose was included as a covariate as well as its interaction with the natural log of day post injection, and person was fitted as a fixed factor. Replicate was included in the Pr. burhodogranariea analysis only, because of unequal sampling across replicates for E. cloacae. L. lactis injected flies were not analysed because only 4 out of 39 (10.3%) cleared the infection. Ps. entomophila infected flies were not statistically analysed because of a low sample size (n = 12). The two bacterial species were analysed separately.

$${{{{{rm{Model}}}}}},7 :,{{{{{rm{CFU}}}}}},{{{{{{rm{presence}}}}}}/{{{{{rm{absence}}}}}}}_{{{{{{rm{dead}}}}}}}, sim ,{{log }}({{log }}({{{{{rm{injection}}}}}},{{{{{rm{dose}}}}}})),\ times ,{{log }}({{{{{rm{day}}}}}},{{{{{rm{post}}}}}},{{{{{rm{injection}}}}}}),+,{{{{{rm{replicate}}}}}},+,{{{{{rm{person}}}}}}$$

To test whether the patterns of clearance were similar for live and dead flies we tested whether the proportion of live uninfected flies was a predictor of the proportion of dead uninfected flies. We separately summed up the numbers of uninfected and infected flies for each bacterial species and dose, giving us a total sample size of n = 20 (four species × five doses). For live and for dead homogenised flies we had a two-vector (proportion infected and proportion uninfected) response variable, which was bound into a single object using cbind. The predictor was live flies, and the response variable was dead flies, and it was analysed using a generalized linear model with family = quasibinomial.

$${{{{{rm{Model}}}}}},8:,{{{{{rm{cbind}}}}}}({{{{{rm{dead}}}}}},{{{{{rm{uninfected}}}}}},,{{{{{rm{dead}}}}}},{{{{{rm{infected}}}}}}), sim ,{{{{{rm{cbind}}}}}}({{{{{rm{live}}}}}},{{{{{rm{uninfected}}}}}},,{{{{{rm{live}}}}}},{{{{{rm{infected}}}}}})$$

### Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.