1. Genetics and Genomics
  2. Microbiology and Infectious Disease

A high-throughput cytotoxicity screening platform reveals agr-independent mutations in bacteraemia-associated Staphylococcus aureus that promote intracellular persistence

  1. Abderrahman Hachani  Is a corresponding author
  2. Stefano G Giulieri
  3. Romain Guérillot
  4. Calum J Walsh
  5. Marion Herisse
  6. Ye Mon Soe
  7. Sarah L Baines
  8. David R Thomas
  9. Shane Doris Cheung
  10. Ashleigh S Hayes
  11. Ellie Cho
  12. Hayley J Newton
  13. Sacha Pidot
  14. Ruth C Massey
  15. Benjamin P Howden
  16. Timothy P Stinear
  1. Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Australia
  2. Infection and Immunity Program, Department of Microbiology and Biomedicine Discovery Institute, Monash University, Australia
  3. Biological Optical Microscopy Platform, University of Melbourne, Australia
  4. School of Microbiology, University College Cork, Ireland
  5. School of Medicine, University College Cork, Ireland
  6. APC Microbiome Ireland, University College Cork, Ireland
  7. School of Cellular and Molecular Medicine, University of Bristol, United Kingdom
  8. Microbiological Diagnostic Unit Public Health Laboratory, Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Australia
https://doi.org/10.7554/eLife.84778
Share this article
Cite this article
  1. Abderrahman Hachani
  2. Stefano G Giulieri
  3. Romain Guérillot
  4. Calum J Walsh
  5. Marion Herisse
  6. Ye Mon Soe
  7. Sarah L Baines
  8. David R Thomas
  9. Shane Doris Cheung
  10. Ashleigh S Hayes
  11. Ellie Cho
  12. Hayley J Newton
  13. Sacha Pidot
  14. Ruth C Massey
  15. Benjamin P Howden
  16. Timothy P Stinear
(2023)
A high-throughput cytotoxicity screening platform reveals agr-independent mutations in bacteraemia-associated Staphylococcus aureus that promote intracellular persistence
eLife 12:e84778.
https://doi.org/10.7554/eLife.84778
  2. Download BibTeX
  3. Download .RIS

Abstract

Staphylococcus aureus infections are associated with high mortality rates. Often considered an extracellular pathogen, S. aureus can persist and replicate within host cells, evading immune responses, and causing host cell death. Classical methods for assessing S. aureus cytotoxicity are limited by testing culture supernatants and endpoint measurements that do not capture the phenotypic diversity of intracellular bacteria. Using a well-established epithelial cell line model, we have developed a platform called InToxSa (intracellular toxicity of S. aureus) to quantify intracellular cytotoxic S. aureus phenotypes. Studying a panel of 387 S. aureus bacteraemia isolates, and combined with comparative, statistical, and functional genomics, our platform identified mutations in S. aureus clinical isolates that reduced bacterial cytotoxicity and promoted intracellular persistence. In addition to numerous convergent mutations in the Agr quorum sensing system, our approach detected mutations in other loci that also impacted cytotoxicity and intracellular persistence. We discovered that clinical mutations in ausA, encoding the aureusimine non-ribosomal peptide synthetase, reduced S. aureus cytotoxicity, and increased intracellular persistence. InToxSa is a versatile, high-throughput cell-based phenomics platform and we showcase its utility by identifying clinically relevant S. aureus pathoadaptive mutations that promote intracellular residency.

Editor's evaluation

This paper describes a new method to investigate Staphylococcus aureus intracellular virulence that has produced important insights into the mechanisms of staphylococcal pathogenesis. The results are convincing and the methodology is state-of-the-art. The authors have responded to the reviewer comments and resolved the issues identified during the review. This paper will be of interest to scientists studying microbial intracellular pathogenesis and cell biology.

https://doi.org/10.7554/eLife.84778.sa0

Introduction

Staphylococcus aureus is a leading cause of hospital-acquired infections, a problem exacerbated by increasing resistance to last-line antibiotics (Murray et al., 2022; Tong et al., 2015). S. aureus is traditionally considered an extracellular pathogen as it produces many secreted virulence factors, including superantigens, degradative enzymes, and cytolytic toxins (Tam and Torres, 2019). The potent pore-forming toxins (PFTs), including alpha-hemolysin and leukocidins, are among the major virulence determinants of S. aureus (Cheung et al., 2021). These toxins induce rapid host cell death, including death of the leukocytes and neutrophils recruited to remove bacteria from infected tissues (Surewaard et al., 2016; Thammavongsa et al., 2015).

Long-term asymptomatic S. aureus carriage is common but invasive infection is rare. Thus, understanding the changes enabling S. aureus to switch from a common coloniser of the nasal cavity to an invasive pathogen is a major research focus. A powerful discovery approach has been used to compare the cytolytic activities of secreted S. aureus virulence factors between colonising and bacteraemia isolates, followed by genomics to uncover the bacterial genetic loci linked with high/low cytolytic activity (Collins et al., 2008; Giulieri et al., 2018; Laabei et al., 2014; Laabei et al., 2021; McConville et al., 2022). These toxicity analyses use methods that monitor eukaryotic cell viability upon exposure to S. aureus culture supernatants (Dankoff et al., 2020; Das et al., 2016; Giulieri et al., 2018) and integrate these phenotypes within genome-wide association studies (GWAS) and other comparative bacterial population genomics techniques (Giulieri et al., 2018; Recker et al., 2017). However, such toxicity assays are limited in that they focus on exogenous virulence factors that have accumulated in the culture media during bacterial growth. Thus, such methods can be limiting as they report phenotypes that are temporally skewed and ignore the host cell-bacterium context under which the production of these factors is controlled during infection. When these endpoint toxicity assays are conducted at scale with many hundreds of bacterial isolates, additional issues caused by filter-sterilisation, freeze-storage, and other manipulations during the preparation of bacterial supernatants prior to incubation with host cells may increase assay variability (Giulieri et al., 2018; McConville et al., 2022).

Another restriction of toxicity assays is that they treat S. aureus as an obligate extracellular pathogen, whilst the literature has reported its capacity to adopt an intracellular behaviour, readily adhering to and invading various eukaryotic cells, including non-professional phagocytes (Flannagan et al., 2016; Foster et al., 2014; Soe et al., 2021). Upon internalisation, intracellular S. aureus initially reside in phagolysosomes, where low pH is a cue for bacterial replication (Flannagan et al., 2018; Lâm et al., 2010). S. aureus uses its arsenal of PFTs to escape from this degradative compartment, into the cytosol and cause the death of host cells (Moldovan and Fraunholz, 2019; Siegmund et al., 2021). This intracellular niche confers protection to S. aureus from antibiotics and immune responses (Peyrusson et al., 2020; Strobel et al., 2016). While guarding S. aureus from host attack, intracellular residency also enables the creation of a reservoir for the pathogen to persist in an infected host and could lead to bacterial transmigration into distal host tissues, from where the bacteria can cause protracted infections and more serious disease (Jorch et al., 2019; Surewaard et al., 2016). Toxicity and persistence of S. aureus in an intracellular context are critical pathogenesis traits. However, understanding these traits and their microevolution across diverse collections of clinical S. aureus strains, and correlating them with specific bacterial genetic variations has been hampered by the lack of a high-throughput platform for trait characterisation.

To address these issues, we took advantage of the capacity of S. aureus to invade cultured epithelial cell lines and established a 96-well assay format to accurately monitor over time the bacterial toxicity exerted from within host cells. HeLa cells are adherent epithelial cells and represent an amenable model to study the pathogenesis of most intracellular bacteria causing human disease, including S. aureus (Das et al., 2016; Stelzner et al., 2020a; Stelzner et al., 2020b). We modified an antibiotic/enzyme protection assay, using a combination of gentamicin and lysostaphin, to kill extracellular bacteria while preserving the viability of intracellular bacteria (Kim et al., 2019). Using propidium iodide (PI) as a marker of host cell death, the assay measured changes in fluorescence of HeLa cells over time caused by the presence of intracellular S. aureus. Combined with the power of bacterial genomics and evolutionary convergence analysis, we used the assay to screen a large collection of S. aureus bacteraemia isolates. Our large-scale pheno-genomics approach revealed known and previously undescribed loss-of-function mutations that were significantly associated with reduced intracellular cytotoxicity and increased intracellular persistence.

Results

InToxSa assay development

We set out to develop a high-throughput and continuous cell death assay that could measure the intracellular toxicity of S. aureus in a format we named InToxSa. We used adherent HeLa-CCL2 epithelial cells as an infection model in a 96-well format and infected them with S. aureus at a standardised multiplicity of infection (MOI). Following an infection period of 2 hr, infected cells were treated with an antibiotic/enzyme combination to specifically kill extracellular bacteria and prevent further reinfection by bacteria released by cells dying during the assay. HeLa cell viability was continuously monitored by measuring PI fluorescence (see methods). Reduced HeLa cell viability was indicated by increasing fluorescence over time (Figure 1A). We used regression to fit standardised curves to the PI uptake data and calculated seven kinetic parameters including the area under the curve (AUC) representing the total of PI uptake over time, peak PI uptake [μmax], the time to μmax [t(μmax)], the maximum rate in PI uptake [rmax], the time to rmax [t(rmax)], trough, and time to trough (Figure 1—figure supplement 1).

Figure 1 with 1 supplement see all
Download asset Open asset
Establishing the intracellular toxicity of S. aureus (InToxSa) assay.

(A) Overview of InToxSa assay. (B) Flow chart of the analytical pathway. (C) Plot of propidium iodide (PI) fluorescence uptake over 20 hr as a measure of S. aureus intracellular cytotoxicity in HeLa cells. Depicted curves are wild-type S. aureus JE2 (blue), the isogenic S. aureus JE2 agrA transposon mutant (red), and uninfected cells (yellow). The PI uptake curve for JE2 is annotated with five kinetic parameters. For each curve, the thick line represents the mean and the shading, the standard deviation. Curves are fitted with cubic smooth splines (see methods). To minimise batch effect, all kinetics data have been transformed using proportion of maximum scoring (POMS) using JE2 controls as reference minimum and maximum values (Little, 2013). x-axis is time and y-axis is PI uptake, represented as a proportion of maximal fluorescence in JE2-infected cells, where for every measured plate, a PI uptake value of 1 represents the maximum of JE2 PI uptake and zero its minimum. (D) Summary of five independent InToxSa experiments to assess assay and parameter variation. Violin plots represent the density distribution of all five replicates (based on Gaussian kernel density estimation) and the nested box plots (boxes represent median, first and third quartiles, whiskers represent are based on ×1.5 interquartile range) show the distribution of within plate replicates (3–5 technical replicates per plate replicate) for the three most discriminatory of seven parameters inferred from the PI uptake data (Figure 1—figure supplement 1).

We then mapped out a series of experiments to validate InToxSa and explored bacterial genetic factors linked to intracellular cytotoxicity (Figure 1B). To demonstrate method performance, we measured the intracellular toxicity of the wild-type S. aureus USA300 strain JE2 against an isogenic agrA mutant (Nebraska Transposon Library mutant NE1532; Fey et al., 2013), using non-infected cells as a baseline (Figure 1C). S. aureus JE2 caused a rapid and substantial increase in PI fluorescence over time, reflective of the known high cytotoxicity of this strain (Das et al., 2016; Grosz et al., 2014; McConville et al., 2022). Cells infected with the agrA mutant yielded significantly lower PI uptake (AUC) and slower (rmax), indicating HeLa cell viability during the infection course and is consistent with the reported low cytotoxicity of S. aureus agr mutants (Figure 1C; Laabei et al., 2021; McConville et al., 2022).

We then assessed the reproducibility and repeatability of InToxSa across five experimental replicates, each time using both biologically independent HeLa cell culture passages and independent S. aureus cultures with the same two comparator strains (JE2 wild type and the isogenic agrA transposon mutant) (Figure 1C). We measured seven PI uptake curve kinetic parameters (Figure 1C, Figure 1—figure supplement 1). We observed that the PI uptake AUC, peak PI uptake [μmax], and maximum PI uptake rate [rmax] for S. aureus JE2 compared to the agrA mutant and non-infected cells were significantly different, had very low intra-assay variation, and were among the most discriminatory and reproducible PI uptake curve parameters (Figure 1D, Table 1, Figure 1—figure supplement 1). At experimental endpoints, the acidity of the culture media had not changed, and no bacterial growth was observed after plating the media from infected 96-well plates, indicating that InToxSa assessed the cytotoxicity caused only by intracellular S. aureus (data not shown).

Table 1
Summary of intracellular toxicity of S. aureus (InToxSa) assay performance.
StrainNo. biological replicatesArea under the curve [AUC]Peak uptakeMax uptake
[μmax][rmax]
MeanStd. dev.CoVMeanStd. dev.CoVMeanStd. dev.CoV
S. aureus JE2 wild type2516214.60.0910.090.090.0030.00040.13
S. aureus JE2 agrA mutant1130.9140.450.170.060.350.00040.00020.52
Non-infected155.318.33.460.070.071.020.00010.000050.32

  1. Note: Std. dev. = standard deviation; CoV = coefficient of variation.

Confirmation of S. aureus internalisation in HeLa cells

We used confocal microscopy to confirm the presence of intracellular S. aureus (Figure 2A). HeLa cells grown on coverslips were infected with the same conditions used for InToxSa and analysed at 3, 8, and 24 hr post-infection. These timepoints were selected to reflect the key InToxSa timepoints highlighted during JE2 infection (start of PI uptake measurement, peak PI uptake [μmax], and experimental endpoint, respectively). We observed that at each timepoint, both JE2 and the agrA mutant co-localised with the lysosomal marker LAMP-1. However, the agrA mutant was detected in higher numbers within cells compared to wild type (Figure 2A, B). At later timepoints (8 and 24 hr) the number of JE2-infected cells decreased, with fewer detectable intracellular bacteria, suggesting that JE2-infected cells had died, consistent with bacterial cytotoxicity. In contrast, the number of cells infected with the agrA mutant did not vary significantly, indicating cell survival during the infection time course (Figure 2B). It also appeared that the number of intracellular bacteria increased for cells infected with the agrA mutant, suggesting bacterial replication over time. This latter observation was explored in more depth using high-content/high-throughput imaging (see later section). The microscopy results support the InToxSa assay outputs and indicate that non-cytotoxic S. aureus, such as the agrA mutant, can persist within cells without affecting cell viability.

Figure 2
Download asset Open asset
Fluorescence confocal microscopy of intracellular S. aureus.

(A) HeLa cells were infected with S. aureus (wild-type JE2 or isogenic agrA mutant) and imaged at 3, 8, and 24 hr post-infection. Fixed cells were labelled with LAMP-1, S. aureus antibodies and 4′,6-diamidino-2-phenylindole (DAPI). (B) Manual quantification of confocal microscopy. Graph shows the percentage of cells infected with S. aureus at each of the three timepoints. At least 50 cells (n cells = 51–112) were counted in 3–5 fields of view, with at least 12 cells counted per field (n field = 12–40). Shown are all data points, mean, and standard deviation. Significance was assessed using two-way analysis of variance (ANOVA). Null hypothesis (no difference between means) rejected for adj p < 0.05. *p = 0.04, ***p = 0.007, ****p = <0.0001.

InToxSa benchmarking against trypan blue exclusion assay

In a previous study, we used a trypan blue exclusion assay with THP1 human macrophages exposed to culture supernatants to test within-host cytotoxicity variations from 51 clinical S. aureus isolated from 20 patients with bacteraemia (Giulieri et al., 2018). These 51 S. aureus isolates were originally selected because they were associated with phenotypic changes occurring during bacteraemia, such as infection relapse, antibiotic treatment failure, longer duration of bacteraemia, and augmented vancomycin MIC; phenotypes likely resulting from within-host evolution (Giulieri et al., 2018). Thus, the isolates were categorised by patient and by their sequential isolation, with the first isolate considered as ‘baseline’ and the subsequent isolates designated as ‘evolved’.

All 51 isolates were screened with InToxSa and compared with the original trypan blue assay data (Figure 3). Using trypan blue exclusion, only the evolved isolate from patient 50 (P_50) showed a significant phenotypic difference in cytotoxicity (Figure 3A). This difference was attributed to a loss-of-function mutation in agrA (T88M). InToxSa also identified a significant cytotoxicity difference for the P_50 isolate pair, validating the previous observation, despite the methodological and host cell type differences. However, InToxSa detected significantly reduced cytotoxicity for the evolved isolates from five more patient pairs compared to the original trypan blue screen (Figure 3B). These results support the higher sensitivity of InToxSa in uncovering S. aureus cytotoxicity variations resulting from the evolution of the bacterial population during bloodstream infection.

Figure 3
Download asset Open asset
Performance of intracellular toxicity of S. aureus (InToxSa) against trypan blue exclusion assay.

Comparative evaluation of S. aureus intracellular cytotoxicity with bacterial supernatants for 51 paired isolates from 20 patients with S. aureus bacteraemia. (A) Supernatant-based cytotoxicity on THP1 cells. Episode with significant difference in THP1 survival between baseline and evolved isolates is boxed in blue (p < 0.05). Assay performed in biological and technical duplicates. Bars represent mean percentage of dead cells; error bars show range between duplicates. Toxicity within isolate groups was compared using analysis of variance (ANOVA) with Bonferroni correction (B) propidium iodide (PI) uptake of infected HeLa cells. Values are mean area under the curve (AUC) and standard deviation. Episodes exhibiting significant phenotypic differences between baseline and evolved isolates are boxed in red (p < 0.05). Assay performed in biological and technical triplicates. PI uptake AUC within isolate groups was compared using ANOVA with Bonferroni correction.

Screening a large collection of clinical S. aureus to evaluate intracellular cytotoxicity

A major motivation for developing the InToxSa assay was to develop a pheno-genomics platform to efficiently measure the intracellular cytotoxicity profiles of a large collection of clinical S. aureus isolates and then use the power of comparative and statistical genomics to find bacterial genetic loci associated with intracellular cytotoxicity. To this end, we analysed 387 clinical S. aureus isolates, obtained from 298 episodes of bacteraemia and for which genome sequences were available (Giulieri et al., 2018; Holmes et al., 2011; Holmes et al., 2014; Holmes et al., 2018). A 164,449 single nucleotide polymorphism (SNP) core-genome phylogeny was inferred for this collection. The 387 isolates spanned 32 sequence types (STs) and were dominated by ST239, accounting for 30% of isolates (n:117), followed by ST22 (n:32, 8%), ST5 (n:34, 8%), ST45 (n:28, 7%), and ST1 (n:18, 5%). Fifty-three percent of the isolates were mecA positive (Figure 4A; Supplementary file 3).

Figure 4 with 4 supplements see all
Download asset Open asset
Intracellular cytotoxicity assessment of 387 bloodstream-associated clinical S. aureus isolates.

(A) Maximum likelihood phylogeny based on 164,449 core-genome SNPs for 387 S. aureus, showing sequence type (ST) and MRSA distribution. The heatmap depicts the mean area under the curve (AUC) of cytotoxicity based on intracellular toxicity of S. aureus (InToxSa) propidium iodide (PI) uptake assay. AUC values range from non-cytotoxic (score: 0, dark blue) to highly cytotoxic (score: 200, yellow). Adjacent to the heatmap (closed and open circles) are 28 pairs of genetically related, but phenotypically discordant isolates (see Figure 5). (B) ‘Variable importance plot’ showing different PI uptake metrics (features) in an unsupervised random forest (RF) machine learning model. The higher the value of ‘accuracy decrease’ or ‘Gini index decrease’, the higher the importance of the feature in the model. (C) Scatter plot of the two most discriminatory PI uptake kinetic metrics (AUC and maximum PI uptake rate). Dots are coloured based on the clustering obtained from the proximity matrix of the RF model. (D) Scatter plot showing the two principal components with the strongest association with PI uptake (lineage effect as measured using pyseer). Dots and ellipses are coloured based on the clustering obtained from the proximity matrix of the unsupervised machine learning model. (E) Manhattan plot of gene-burden genome-wide association study (GWAS) of cytotoxicity (PI uptake AUC) of 387 clinical isolates.

We assessed each of the 387 isolates by InToxSa, with the cytotoxicity phenotype for each isolate represented by mean PI uptake (AUC) and displayed as a heatmap, aligned with the core-genome phylogeny (Figure 4A). Several patterns were readily seen in the data. There was a wide range of cytotoxicity profiles across the 387 isolates, with notable variations within each ST suggesting frequent adaptation events affecting intracellular cytotoxicity levels. Two STs (ST239 and ST22) strongly associated with lower cytotoxicity. These variations are consistent with previous observations (Laabei et al., 2021).

To select the most suitable InToxSa parameters for our statistical genomics approach, we assessed their relative importance by fitting the data using unsupervised random forest (RF) machine learning (Mantero and Ishwaran, 2021). In this model, ‘observations’ were each of the PI uptake values for the 387 isolates and controls and ‘features’ were the seven PI uptake curve parameters (Supplementary file 4). We then tested RF feature importance and showed with variable importance plots that the PI uptake parameters, AUC, μmax, and rmax were the most informative features for the model (Figure 4B), consistent with our initial InToxSa assessment using JE2 and the agrA mutant (Figure 1C, D). Using the proximity matrix from the unsupervised RF model, we defined three main PI uptake clusters, corresponding with low, moderate, and high intracellular cytotoxicity categories. We labelled each of the 387 PI uptake data points with these three (low, moderate, and high) cytotoxicity categories and plotted the AUC and rmax values against each other (Figure 4C). As expected, these parameters were strongly, positively correlated, suggesting that the AUC alone is sufficient to capture intracellular cytotoxicity differences between S. aureus isolates. We used principal component analysis of the PI uptake data as an alternative unsupervised learning approach (Figure 4—figure supplement 1). When considering the first two components (67% of the variance explained), we observed a similar pattern where the same toxicity groups could be recognised within a cytotoxicity continuum among clinical isolates (Figure 4—figure supplement 1).

GWAS analysis using InToxSa outputs to identify S. aureus genes linked to intracellular cytotoxicity

We next used GWAS to identify genetic correlates of strain-level cytotoxicity, expressed as mean PI uptake AUC. The fraction of cytotoxicity variation explained by genetic variation (heritability: h2) was 49%, a figure lower than the ones obtained for other phenotypes such as vancomycin resistance (Giulieri et al., 2022b). A lower heritability could be resulting from the InToxSa assay variability or caused by differences in gene expression levels or due to epigenetic changes.

To assess the contribution of lineage effects relative to locus effects, we defined lineages using multidimensional scaling (MDS) of a pairwise genetic distance matrix generated by Mash, a tool that reduces genome content to a set of ‘sketches’ (hashed k-mers) (Ondov et al., 2016). Major MDS axes correlated with the most prevalent STs, for example ST239 was mainly defined by MDS1 (negative correlation) and MDS2 (positive correlation) (Figure 4—figure supplement 3). We then tested the association between the first 10 MDS axes (90% of the genetic variance explained) and the PI uptake phenotype in Pyseer (Earle et al., 2016; Giulieri et al., 2022b; Lees et al., 2018). In agreement with the initial observations based on the phylogeny and cytotoxicity heatmap (Figure 4A), we observed significant cytotoxicity-lineage associations represented by MDS3 and MDS4 (Figure 4D, Figure 4—figure supplement 2). Because of the ST-MDS lineages correlation, this is consistent with differences in cytotoxicity between clones (Figure 4—figure supplement 3B). Using the three cytotoxicity clusters defined by RF as categorical labels (Figure 4C), we plotted the 387 genomes along these two dimensions. While intracellular cytotoxicity was strongly associated with some S. aureus lineages, this analysis showed that lineage alone does not completely explain the phenotype, as indicated by the significant overlap between the three cytotoxicity clusters across MDS3 and MDS4 (Figure 4D). This pattern is consistent with other adaptive phenotypes (Earle et al., 2016; Giulieri et al., 2022b; Su et al., 2021) and suggests that locus effects from specific microevolutionary events modulate cytotoxicity, supporting the use of GWAS and convergent evolution approaches to identify these mutations.

Correcting for the observed population structure, we then used gene-burden GWAS to try and identify S. aureus loci significantly associated with intracellular cytotoxicity (PI uptake AUC) as a continuous variable. After correcting for multiple testing, only agrA reached the p < 0.05 significance threshold, supporting the important contribution of this locus to strain-level cytotoxicity (Figure 4E, Supplementary file 5). We also considered the highest ranking loci that did not reach genome-wide statistical significance. The second most significant gene, secA2 (p = 1.5 × 10−4) encodes the accessory ATPase to the Sec protein export system and is essential for the transport of SraP, a surface exposed and serine-rich staphylococcal protein which is associated with adhesion to, and invasion of epithelial cells and binding to human platelets (Siboo et al., 2005; Yang et al., 2014). Another high-ranking GWAS locus was ileS (p = 9.9 × 10−4), encoding an isoleucyl-tRNA synthetase linked with mupirocin resistance, and previously associated with S. aureus cytotoxicity (Yokoyama et al., 2018).

Identification of convergent mutations in genetic pairs with divergent InToxSa cytotoxicity profiles

Despite the relatively small sample size for this kind of analysis, the gene-burden GWAS detected the agrA locus with a high significance, but it did not have sufficient power to detect mutations other than those occurring in agr genes. This was expected with our sample size as demonstrated by our GWAS power calculations (Figure 4—figure supplement 4). Our power calculations have also showed that reducing the number of variants could substantially increase the yield of the analysis. This can be achieved by limiting the analysis to closely related isolates, with the added advantage of eliminating the bias related to the population structure. Therefore, we sought to identify rare mutations that might alter the intracellular cytotoxicity using comparative genomics approaches, a complementary strategy to microbial GWAS (Chen and Shapiro, 2021; Giulieri et al., 2022b; Saund and Snitkin, 2020). We used evolutionary convergence analysis to identify additional loci associated with intracellular cytotoxicity among the 387 S. aureus isolates. Our approach was to identify genetically related pairs of isolates with contrasting PI uptake AUC values from independent clades across the phylogeny and then search for homoplasic mutations between the pairs. We calculated genetic distances between all 387 genome-pairwise comparisons (149,769 combinations) and calculated a delta-PI uptake AUC value for each pair. We selected 28 phylogenetically independent S. aureus pairs (56 unique isolates) with a genetic distance <200 core-genome SNPs and a significant decrease in PI uptake AUC between reference (isolate-1) and control (isolate-2) (Wilcoxon rank-sum test) (Figure 5A). Variants within each pair (i.e., found in isolate-2 but not in isolate-1) were identified and annotated using a strategy that we have developed for S. aureus within-host evolution analysis (Giulieri et al., 2022a). We have previously shown that an SNP-calling approach using de novo assembly of one genome in a pair as a reference provided the most accurate estimate of the genetic distance (Higgs et al., 2022). There were between 0 and 206 mutations within the 28 pairs (Figure 5B). Mapping the genes in which these mutations were found back to a core-genome phylogeny constructed from the 56 paired S. aureus genomes showed convergent (homoplasic) mutations in agrA, supporting the GWAS findings. In addition, our analysis identified potentially convergent mutations in several other genes (6 with three independent mutations and 35 with two independent mutations) (Figure 5C–E, Supplementary file 6). However, because of the strong lineage effect and the paucity of representation for some S. aureus lineages (clonal complexes [CC] and STs), half of these mutations were only found in ST239, a well-represented lineage in our collection. In addition to target loci dedicated to the regulation of virulence factors such as the agr locus or involved in adhesion to host extracellular matrix proteins such as fibronectin and elastin (fnbA and ebpS), some of the convergent mutations were found in genes involved in metabolic processes (ribA, purF, sbnF, ilvB, lysA, and araB), associated with the cell wall (fmtB), devoted to the last step of the cell wall teichoic acid biosynthesis (tarL’), implicated in DNA repair (mfd), in protein transport (secA2), in solute transport (glcA, opp-3A, and thiW), in respiration (cydA) and found in a phage-associated locus (SAUSA300_1930). Aside from agrA and agrC genes (Giulieri et al., 2018; Laabei et al., 2014; Laabei et al., 2015; Mairpady Shambat et al., 2016), and those found in the promoter of the tar locus (Brignoli et al., 2022; Laabei et al., 2014), mutations in the other loci have not been associated with the reduction of cytotoxicity in clinical S. aureus isolates. Interestingly, homoplasic mutations were also found in the gene ausA, known to be involved in S. aureus escape from epithelial cell endosomes and the phagosome of phagocytic cells (Blättner et al., 2016).

Figure 5 with 3 supplements see all
Download asset Open asset
Evolutionary convergence analysis to identify S.aureus genes linked with intracellular cytotoxicity.

(A) Distribution of genetic distance determined by pairwise comparisons using MASH distances among the 387 S. aureus genomes, against the difference in propidium iodide (PI) uptake area under the curve (AUC) between each pair. The shaded circles denote the 95% multivariate t-distribution (blue: pairs included in the convergence analysis; red: pairs excluded from the analysis). (B) Ranked distribution of the difference in PI uptake AUC between the 28 pairs. Heatmap shows reduction in AUC values. (C) Core-genome phylogeny for the 28 pairs of isolates. Tree tips are coloured by PI uptake AUC. Aligned with the phylogeny, the 10 first genes targeted by convergent mutations are shown. (D) Number of mutations detected for each of the 20 genes, coloured by S. aureus ST. (E) Location of convergent mutations in each gene (non-synonymous in orange, truncating in maroon), purple triangles indicate the position of transposon inserted in tested transposon mutant. (F) Effect of loss of function for each of the 20 genes on intracellular cytotoxicity measured by intracellular toxicity of S. aureus (InToxSa), using mutants from the Nebraska transposon library. Dotted line shows mean PI uptake AUC of positive control strain JE2. Depicted are mean (dot) and standard deviation (SD; bar) of biological triplicates. Mutants causing significantly lower PI uptake AUC to JE2 are depicted in light blue, non-significant changes are in dark blue (Wilcoxon rank-sum test, corrected for multiple testing). (G, H) Operetta high-content imaging analysis for each of the 20 Nebraska transposon mutants and JE2 positive control. Heatmaps show the percentage of HeLa cells infected with each transposon mutant (blue) and the number of bacteria per infected cells at 3 and 24 hr post-infection (red).

We focussed our analysis on genetic pairs with loss of cytotoxicity, to be consistent with the hypothesis that adaptive mutations arising during infection would promote a low-toxicity phenotype as previously shown (Giulieri et al., 2018; Laabei et al., 2014). This hypothesis is supported by our gene-burden GWAS, where 10 out of 10 most significant hits are associated with loss of cytotoxicity (mean normalised PI AUC decrease = −0.8, Supplementary file 5 and Figure 4E). To explore whether including pairs with toxicity gain would alter the analysis, we selected four additional phylogenetically independent pairs with gain of toxicity. Adding these pairs did not reveal new genes with convergent mutations, however, this analysis revealed an additional homoplasic mutation in ausA (Figure 5—figure supplement 1).

To further confirm the roles of genes bearing homoplasic mutations in pairs with discordant phenotype, we compared PI uptake between pairs with and without homoplasic mutations. We ran the analysis in 100 replicates of independent pairs combinations and assessed the associations between within-pairs mutations (aggregated at the gene level) and PI uptake difference within the pair using linear regression. For 6 out of the 20 genes originally identified in the genetic pairs analysis (including agrA, glcA, ribA, fmtB, sbnF, and tarL’), this analysis showed that for >95% of replicates, the regression coefficient was above 0, thus providing a robust statistical support for our findings (Figure 5—figure supplement 2). For the gene ausA, 86% of replicates had a beta value above 0, however 78% of agrC regression replicates had a beta value below 0.

Functional assessment of genes with convergent mutations

To assess the functional consequences of the convergent mutations (caused by at least two homoplasic mutations per gene), we again turned to the Nebraska transposon library and selected transposon mutants for 20 genes we had identified (Figure 5D). We used InToxSa to assess the effect of gene disruption on the intracellular cytotoxicity phenotype for each mutant compared to the JE2 wild type (Figure 5F). Over the 20 transposon mutants tested, six showed a statistically significant reduction in cytotoxicity, namely NE1532 (agrA), NE119 (ausA), NE188 (mfd), NE873 (agrC), NE1140 (lysA), and NE117 (cydA). Strains with transposon insertions in agrA, agrC, and ausA showed a highly significant reduction in PI uptake AUC (adjusted p = 5.4 × 10−4, Wilcoxon rank-sum test), confirming their reported roles in affecting bacterial cytotoxicity and validating our convergence analysis (Figure 5F; Blättner et al., 2016; Das et al., 2016; Laabei et al., 2021; Mairpady Shambat et al., 2016). We extended this analysis and used high-content, high-throughput microscopy to observe and quantify in an unbiased manner the impact of each mutation on the S. aureus infectivity and intracellular persistence (see methods). There was an inverse relationship between InToxSa and high-content imaging outputs, with the three mutants most reduced in cytotoxicity showing both a higher percentage of infected cells recovered after 24 hr of infection, and a high number of bacteria per infected cell at 24 hr post-infection as compared to the wild-type control JE2 (Figure 5G, H – Figure 5—figure supplement 3).

Functional assessment of specific convergent mutations

To further assess the impact of specific convergent mutations on intracellular cytotoxicity, we used site-directed mutagenesis in S. aureus BPH3370 (ST239) to recreate a subset of the convergent mutations. We selected isolate BPH3370 for these experiments as it displayed high InToxSa PI uptake AUC (comparable to JE2, Figure 4C) cytotoxicity without bearing any of the convergent mutations we intended to introduce. We focussed our attention on mutations likely to affect protein function and based on the attenuation in cytotoxicity of the cognate transposon mutants. We selected six mutations, previously not documented nor characterised, including non-synonymous mutations leading to residue substitution (agrA E7K and cydA R390C [a reversion of C390R]), frameshifts leading to truncated proteins (agrC G310 frameshift, ausA K2308 frameshift, and lysA K354 frameshift), and introduction of a stop codon (mfd W568 stop codon) in the sequences of convergent genes (Figure 6A and Figure 6—figure supplement 1). We then used InToxSa to assess the cytotoxicity of each targeted mutant, compared to BPH3370 wild type, JE2 and the corresponding JE2 Nebraska transposon mutant, for each of the six loci (Figure 6B, Figure 6—figure supplement 1). We observed that recreation of the E7K agrA mutation, the agrC G310 frameshift mutation and the ausA K2308 frameshift mutations lead to a significant reduction in intracellular cytotoxicity in BPH3370 (Figure 6B). However, the W568 stop codon mfd mutation, the K354 frameshift lysA mutation, and the cydA R390C mutation had no significant effect on the cytotoxicity of the BPH3370 strain (Figure 6—figure supplement 1). It is noteworthy that transposon insertions in these three genes also had a less pronounced effect on the phenotype of JE2 strain as compared to the agr and ausA loci (Figure 5F).

Figure 6 with 1 supplement see all
Download asset Open asset
Introduction of convergent agrA, agrC, and ausA mutations in the clinical isolate BPH3370 reduces its intracellular cytotoxicity while ausA mutation affects aureusimine B production.

(A) Position and nature of convergent mutations identified in the genes agrA, agrC, and ausA. For each gene, the amino acid position affected by mutations is shown on the x-axis for each gene. Convergent mutations causing a significant contrasting propidium iodide (PI) uptake phenotype are coloured according to their consequence on protein function: non-synonymous (orange), truncating characterised by the introduction of a frameshift (fs) or a stop codon (*) (maroon). (B) Effect of convergent mutations on the intracellular cytotoxicity of the clinical isolate BPH3370. The PI uptake area under the curve (AUC) values for JE2, the cognate Nebraska transposon mutants of convergent genes, BPH3370 wild type and BPH3370 bearing the mutations affecting agrA, agrC, and ausA. The crossbar represents mean and standard deviation (p < 0.0001). (C) Predicted impact of K2308 frameshift mutation (K2308fs) on aureusimines. (D) HPLC analysis of S. aureus ethyl-acetate extracts for aureusimines compared to an Aureusimine B synthetic standard.

We predicted that the ausA K2308 frameshift causing a 11 base pair deletion mutation in BPH3370 would lead to a loss of aureusimine biosynthesis. This was because the frameshift occurred within the ausA reducing domain and would thus prevent the release of the dipeptide L-Val-L-TyrT2 to form the intermediate amino aldehyde, with no cyclisation to form the imine (Figure 6C; Zimmermann and Fischbach, 2010). As expected, high-performance liquid chromatography (HPLC) analysis confirmed the absence of aureusimines in the BPH3370 K2308 frameshift mutant, similar to the transposon ausA mutant (NE119) (Figure 6D).

Discussion

Recurrent and persistent staphylococcal infections have been proposed to result from within-host selective pressures leading to the evolution of adaptive traits by the bacteria, a process also observed in other human bacterial pathogens (Didelot et al., 2016; Gatt and Margalit, 2021; Giulieri et al., 2022a). The emergence of mutations affecting regulators controlling toxin production has been proposed as a mechanism enabling S. aureus to adapt to its host while evading cellular immune responses (Giulieri et al., 2022a; Young et al., 2012; Young et al., 2017). Identifying the molecular signatures supporting the pathoadaptation of S. aureus at the host cell interface is important for understanding how S. aureus can cause persistent, difficult-to-treat infections lasting many months (Gao et al., 2015).

Several studies of S. aureus clinical isolates have attempted to identify such signatures by assessing the cytotoxicity of bacterial supernatants applied onto host cells, in an ex-cellulo fashion (Giulieri et al., 2018; McConville et al., 2022; Recker et al., 2017). Such assessments would be adequate if S. aureus was an extracellular pathogen exerting its cytotoxicity from without host cells (Soe et al., 2021), but S. aureus is a facultative intracellular pathogen able to invade and persist in a wide range of eukaryotic cells (Al Kindi et al., 2019; Krauss et al., 2019; Luqman et al., 2019; Musilova et al., 2019; Sinha and Fraunholz, 2010). We developed the InToxSa cytotoxicity assay to address this shortcoming and to try and identify S. aureus pathoadaptive mutations that support a S. aureus intracellular lifestyle. Moreover, the InToxSa assay provides a maximum throughput of 7× 96-well plates per week, thus corresponding to 98 distinct clinical strains testable per week (encompassing 6 individual replicates, each tested across 2 different days/plates) which could be upscaled 4-fold with the adoption of a 384-well plate format. By harnessing the power of comparative and statistical bacterial genomics with InToxSa readouts for a large collection of bacteraemia-associated S. aureus isolates, we identified mutations in S. aureus that reduced the intracellular cytotoxicity and increased intracellular persistence.

We showed the performance and sensitivity of InToxSa with the identification of cytotoxicity differences between S. aureus isolates that had not previously been detected by ex-cellulo methods (Figure 3; Giulieri et al., 2018). The difference in phenotypic outputs for both methods may be in part explained by the different cell types exploited for the readout (THP1 macrophages for trypan blue exclusion assay versus HeLa-CCL2 epithelial cells for InToxSa) and the bacterial fraction examined (culture supernatants versus bacterial cells) (Figure 3). The capacity of InToxSa to detect subtle phenotypes missed by gross cytotoxicity assessments is also conferred by its temporally granular and objective measurements of PI uptake as a marker of host cell viability. InToxSa assesses the S. aureus toxicity caused by bacterial virulence factors produced in response to the intracellular environment and is proportional to a defined bacterial load. This approach contrasts with methods relying on the presence of toxins accumulating over time in bacterial supernatants and whose production relies almost solely on the functionality of the Agr quorum sensing system (Altman et al., 2018; Giulieri et al., 2018; McConville et al., 2022).

We further showed the performance of this analytical pipeline by readily identifying mutations in loci such as the Agr quorum sensing system, which is well known to control S. aureus cytolytic activity (Giulieri et al., 2018; Giulieri et al., 2022a; Laabei et al., 2021; Mairpady Shambat et al., 2016; Recker et al., 2017). However, our approach also enabled discovery of mutations in less characterised systems, including changes in ausA, that reduced S. aureus cytotoxicity and increased intracellular persistence of clinical isolates. AusA is a non-ribosomal peptide synthetase responsible for production of aureusimines, pyrazinone secondary metabolites. Our observations are consistent with previous reports showing that aureusimines contribute to the phagosomal escape of S. aureus JE2 to the cytosol (Blättner et al., 2016; Wilson et al., 2013).

Interestingly, S. aureus mutants that were most affected in cytotoxicity also had a propensity to persist intracellularly (Figure 5). Infected host cells have been proposed as Trojan horses for intracellular S. aureus, increasing the risks of systemic dissemination to organs, such as the liver and kidneys, following bacteraemia and contribute to infection persistence (Jorch et al., 2019; Surewaard et al., 2016; Thwaites and Gant, 2011). Whilst agr-dysfunctional isolates were associated with persistent infections (Fowler et al., 2004; Schweizer et al., 2011), intracellular persistence caused by mutations in agr loci could possibly constitute a population bottleneck (Pollitt et al., 2018; Spaan et al., 2013). However, such population bottlenecks may be transient as it has been suggested that mutations arising in agr defective pathoadapted clinical isolates could possibly compensate for the loss of agr functionality and restore S. aureus virulence, suggesting a stepwise within-host evolution of clinical isolates (Altman et al., 2018; Giulieri et al., 2022a).

Current statistical genomics strategies in human genetics support combining allele-counting methods (GWAS), for the detection of common variants, with comparative genomics approaches to identify rare variants (Singh et al., 2022; Trubetskoy et al., 2022). In microbial genomics, this strategy is best achieved by combining microbial GWAS and convergent evolution studies (Chen and Shapiro, 2021; Giulieri et al., 2022b; Giulieri et al., 2018;(Guérillot et al., 2018) Saund and Snitkin, 2020). Whilst our GWAS approach only identified agrA as significantly associated with low cytotoxicity (Figure 4E), our evolutionary convergence analysis on genetic pairs among our 387 bacteraemia isolates identified mutations in several S. aureus genes that led to reduced cytotoxicity (Figure 5). However, only convergent mutations occurring in agrA, agrC, and ausA were confirmed to affect the cytotoxicity and intracellular persistence phenotypes when introduced into a clinical isolate (Figure 6). This may be due to epistatic effects or combinations of mutations within a specific S. aureus strain may be acting in concert to control the expression of the numerous bacterial cytolytic determinants, underscoring the need to functionally confirm the findings of the convergence analysis. We acknowledge as a limitation of the convergent evolution analysis its reliance on a set of phylogenetic independent pairs, selected by the investigator. Thus, the choice of a different set of pairs may lead to a different list of homoplasic changes. To overcome this limitation, we have performed a repeated pair sampling and assessed the association between homoplasic mutations in PI changes. This approach provided robust support to 6 out of 20 findings of the convergent analysis (including agrA and ausA).

Our approach to integrate GWAS with convergence analyses underscores an important point regarding the choice of samples in genotype–phenotype studies. For future studies, it would be ideal to enrich for genetically related isolates, for example, isolated from the same patient or from a clonal outbreak, if logistically feasible. This would reduce the genetic distance and could therefore increase the yield of the analysis.

Our study also shows that intracellular cytotoxicity levels vary between STs. Despite causing bacteraemia, the ST22 and ST239 isolates were overall less cytotoxic than the ST8 isolates in our collection (as shown on the heatmap in Figure 4A), further corroborated by the direct cytotoxicity comparison between strains JE2 (ST8) and BPH3370 (ST239) (Figure 6B). The evolution of reduced Agr functionality (and thus cytotoxicity) in hospital-acquired ST239 and ST22 isolates has already been reported by our group and others and is confirmed by the InToxSa outputs (Baines et al., 2015; Collins et al., 2008; Giulieri et al., 2018; Laabei et al., 2021; Li et al., 2016). Consistent with their reduced cytotoxicity and with our hypothesis of inverse correlation between toxicity and intracellular replication, ST239 isolates caused higher degrees of bacterial persistence in infected animal models (Baines et al., 2015; Li et al., 2016) and showed increased intracellular persistence in osteoblasts (Bongiorno et al., 2021). Within the limits of our experimental settings, the relatively lower cytotoxicity of ST239 and ST22 isolates indicates that the amplitude of this phenotype should probably be considered within a genetic lineage. The inclusion of representative isolates per lineage, with defined cytotoxicity levels, would identify cytotoxicity thresholds and perhaps allowing identification of more subtle genomic changes affecting phenotypes. Moreover, some of the loci detected by the convergence and GWAS analyses may also have more pronounced effects in some lineages than in others. For instance, mutations affecting tarL and secA2 may affect the export and secretion of virulence factors that are only present in a subset of lineages, thus explaining the absence of effect on cytotoxicity caused by the cognate transposon mutants (Figure 5F).

We developed InToxSa using HeLa cells, a well-defined, adherent, and non-phagocytic cellular model (Das et al., 2016; Stelzner et al., 2020a; Stelzner et al., 2020b). We used adherent epithelial cells because they can be maintained for extended infection periods and so allow the acquisition of useful kinetic measurements of cytotoxicity. However, we also acknowledge the limitation in using these cells in that they do not have the bactericidal modalities of the phagocytes encountered by S. aureus in the bloodstream (Brinkmann et al., 2004; Chow et al., 2020; Krause et al., 2019). Neutrophils are among the first immune cells to engage S. aureus during bacteraemia (Brinkmann et al., 2004). However, neutrophils have a relatively short in vitro lifespan following their purification from blood and would not be well suited to an InToxSa-style assay format (Ge et al., 2020; Rosales, 2020; Tak et al., 2013; Zwack et al., 2022). Polymorphonuclear cell lines such as HL-60, exploited in other ex-cellulo assays, may represent an alternative to primary neutrophils (McConville et al., 2022; Rose et al., 2015). These cells display some of the same important biological functions as neutrophils, including neutrophil extracellular traps (Scieszka et al., 2020), critical in the clearing of S. aureus (Brinkmann et al., 2004; Greenlee-Wacker et al., 2014; Zwack et al., 2022). While PI uptake by these cells could be used as a readout of their viability, HL-60 cells also do not cover all the bactericidal enzymatic activities of primary neutrophils, a potential limitation for their use (Nordenfelt et al., 2009; Yaseen et al., 2017).

We used InToxSa to identify S. aureus pathoadaptive mutations, enriched in bacterial populations that are associated with human disease (e.g., upon transit from colonising to invasive). We hypothesised that these mutations would support an intracellular persistence for S. aureus. Our future research will focus on understanding how these genetic changes might be allowing the bacterium to avoid cell-intrinsic surveillance systems, such as lytic programmed cell death; the self-destructive processes restricting systemic progression of intracellular bacterial pathogens (Wanford et al., 2022). Unlike well-described intracellular Gram-negative bacteria, S. aureus does not have effector proteins to block lytic programmed cell death (Soe et al., 2021). Pathoadaptive mutations such as those arising in the agr locus might prevent cellular injuries caused by S. aureus toxins under Agr control, that would be sensed by cell-intrinsic surveillance platforms such as the inflammasomes and trigger cell death (Krause et al., 2019; Soe et al., 2021). Loss-of-function mutations in ausA, preventing the biosynthesis of aureusimines might be confining S. aureus to a lysosomal compartment where the bacteria have the potential to replicate, and conceivably evade host surveillance mechanisms (Blättner et al., 2016; Flannagan et al., 2016; Grosz et al., 2014; Moldovan and Fraunholz, 2019).

Conclusion

Current large-scale comparative genomics of S. aureus bacteraemia isolates can be further refined by including underexplored pathogenicity traits such as the capacity of S. aureus to invade and survive in host cells. We have addressed this poorly characterised trait of S. aureus pathogenicity by creating the InToxSa assay that measures the intracellular cytotoxicity of many hundreds of S. aureus clinical isolates at scale. We showcase the robustness and reproducibility of phenotypic outputs which, in combination with comparative and statistical genomics, have confidently identified convergent mutations arising in agr and ausA genes that reduced the intracellular cytotoxicity and increased the intracellular persistence of bacteraemia isolates during infection. The adoption of the InToxSa methodology in future pheno-genomics studies would improve the detection of pathoadaptive mutations supporting the persistence and relapse of S. aureus infections.

Materials and methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Gene (S. aureus)BPH2947GenBankAccession GCF_900620245.1
Strain, strain background (Escherichia coli)E. coli strain IM08BMonk et al., 2015IM08BElectrocompetent cells
Cell line (Homo sapiens)Epithelial cellsATCCHeLa-CCL2
Chemical compound, drugAntibioticBaxterGentamicin80 and 40 mg/ml
Recombinant DNA reagentAntibioticAmbiLysostaphin10 mg/ml
Chemical compound, drugNucleic markerSigmaPropidium iodide1 mg/ml
AntibodyAnti-S. aureus (Rabbit polyclonal)WEHI-antibody technology platformCustomisedIF (1:1000)
AntibodyAnti-rabbit
(Donkey polyclonal) coupled to Alexa Fluor 488		Thermo Fisher-InvitrogenCat#: A-21206IF (1:2000)
OtherPhalloidin-TRITCSigmaP1951IF (1:4000)
OtherDAPISigma-Merck10236276001IF (1:2000)

S. aureus isolates

Request a detailed protocol

Clinical isolates were selected from a combined collection of 843 clinical isolates of S. aureus bacteraemia (Giulieri et al., 2018) that was obtained from the vancomycin substudies of the Australian and New Zealand Cooperative on Outcome in Staphylococcal Sepsis (ANZCOSS) study (Holmes et al., 2011) and the Vancomycin Efficacy in Staphylococcal Sepsis in Australasia (VANESSA) study (Holmes et al., 2018). We selected 387 isolates to maximise the likelihood to detect phenotype–genotype associations by sampling different lineages and enriching for episodes where multiple isolates per patient were available (see Supplementary file 3).

Whole-genome sequencing

Request a detailed protocol

After subculturing strains twice from −80°C glycerol stock, DNA was extracted using the Janus automated workstation (PerkinElmer) or manual extraction kits (Invitrogen PureLink genomic DNA kit or the Sigma GenElute kit). Normalised DNA (at a concentration of 0.2 ng/ml) was prepared for sequencing using Nextera XT DNA (Illumina) and sequencing was performed on Illumina MiSeq and NextSeq platforms. Reads quality was assessed based on mean read depth and percentage of S. aureus reads as computed using Kraken2 (Wood et al., 2019). Reads were assembled using Shovill, an assembly pipeline that optimises the Spades assembler (https://github.com/tseemann/shovill) (Bankevich et al., 2012). Annotation was performed using Prokka, with a minimal contig size of 500 bp (Seemann, 2014). Assembly and annotation metrics were used for further quality control of the reads. Genetic distance between clinical isolates was calculated using Mash with a sketch size of 10,000 (Ondov et al., 2016). We used the distance matrix generated by Mash to perform MDS using the function ‘cmdscale()’ in base R. Multi-locus ST were inferred from the assemblies using mlst (https://github.com/tseemann/mlst) (Seemann, 2020a). We assessed the correlation between the most prevalent ST and the MDS axes using the get_correlations() function in the R package bugwas (Earle et al., 2016).

Variants calling: single reference

Request a detailed protocol

Clinical isolates’ reads were mapped to internal reference BPH2947 (accession GCF_900620245.1), an ST 239 reference genome that was generated from the collection. We used snippy, v4.6.0 for mapping and variant calling, with default settings (https://github.com/tseemann/snippy) (Seemann, 2020b). The core-genome alignment was constructed using Snippy-core. We defined core genome as positions where at least 90% of the sequences had a minimum coverage of 10 reads and used Goalign v0.3.4 and SNP-sites v2.5.1 to extract core-genome positions. To infer a maximum likelihood phylogenetic tree of the clinical isolates collection we ran IQ-TREE v2.0.3 using a GTR-G4 model. We used HomoplasyFinder (Crispell et al., 2019) to identify homoplasic sites based on the consistency index. The consistency index was calculated with the following formula: (Number of nucleotides at site − 1)/Minimum number of changes on tree.

Construction of mutants by allelic exchange

Request a detailed protocol

Engineering of convergent mutations in agrA, agrC, ausA, mfd, lysA, and cydA genes in the strain BPH3370 was performed by allelic exchange as described previously (Monk and Stinear, 2021) using oligonucleotides described in Brinkmann et al., 2004; Supplementary file 7 (outlining residues modified by convergent mutations). Upstream and downstream regions of each mutation were PCR amplified and gel extracted and then a splice by overlap extension PCR was performed to generate each insert. Each insert was cloned into linearised pIMAY-Z vector by Seamless Ligation Cloning Extract (SLiCE) cloning (Zhang et al., 2012) to generate six plasmids. Each plasmid was separately transformed into E. coli strain IM08B (Monk et al., 2015) confirmed by colony PCR, then purified and transformed into S. aureus strain BPH3370 by electroporation. Mutant candidates were screened by Sanger Sequencing (Australian Genome Research Facility, Melbourne, VIC, AUS) and positive clones were validated by whole-genome sequencing on an Illumina Miseq or NextSeq550 platforms (Illumina, San Diego, CA, USA) to confirm their genotype. The resultant reads were mapped to the BPH3370 reference genome and mutations were identified using snippy (v4.6.0, https://github.com/tseemann/snippy) (Crispell et al., 2019).

Clinical isolates library preparation

Request a detailed protocol

The collection of clinical isolates was prepared to be readily inoculated from 96-well microtitre plates. Clinical isolates were grown in 10 ml Brain Heart Infusion (BHI) broth (BD Bacto) from single colonies to stationary phase. Briefly, a volume corresponding to 1-unit OD600 for each culture was centrifuged at 10,000 × g for 5 min. The bacterial pellets were washed once with 500 μl of fresh BHI and centrifuged again. The washed bacterial pellets were resuspended in 600 μl of storage media (BHI containing 40% glycerol), vortexed briefly and 200 μl were distributed across 96-well microtitre plates. To prevent operator and plate effect biases, the 387 isolates were randomly distributed with each plate to include 29 distinct isolates, represented in non-contiguous technical triplicates. Built-in controls for cytotoxicity were included in each plate. The wild-type JE2 strain was used as positive cytotoxicity control and the BPH3757 strain, an ST239 isolate bearing the T88M agrA mutation described in Giulieri et al., 2018, as a non-cytotoxic control. Six wells were kept empty to monitor the viability of non-infected controls and account for residual PI uptake. Plates were stored at −80°C. Each plate was at least tested in three biological replicates.

Tissue culture

Request a detailed protocol

HeLa-CCL2 cells (sourced from the ATCC) were maintained and propagated in Dubelcco’s modified Eagle medium (DMEM) + GlutaMAX (4.5 g/l D-glucose and 110 mg/l sodium pyruvate) supplemented with heat-inactivated 10% foetal bovine serum (Gibco) and in absence of antibiotics. The cells were regularly tested for the presence of mycoplasma.

InToxSa assay

Request a detailed protocol

S. aureus isolates were inoculated directly from stabbed frozen parsed plates stock into 100 μl of BHI broth dispatched in flat bottom 96-well microtitre plates. Inoculated plates were incubated for 16 hr in a heat-controlled plate reader (CLARIOstar plate reader, BMG Labtech) set at 37°C. Bacterial growth was assessed by OD600 measurement every 10 min. The endpoint optical densities of cultures were used to infer bacterial density (1-unit OD600 corresponding to 5 × 108 bacteria/ml). Bacterial cultures were standardised and serially diluted in DMEM to reach an MOI 10. 100 μl of bacterial suspension were added to infect 40,000 HeLa-CCL2 cells grown (70% confluence per well) in 96-well black plates, clear bottom (Sigma). Infection was synchronised by centrifugation at 500 × g for 10 min (Eppendorf 5810R) at room temperature. Infected plates were incubated 2 hr at 37°C and 5% CO2 to allow for S. aureus internalisation. The infective media was then discarded, and cells washed once with sterile phosphate-buffered saline (PBS) and further incubated 1 hr with 100 μl DMEM containing cell impermeable antibiotics (80 μg/ml gentamicin [Baxter] and 10 μg/ml of lysostaphin [Ambi]) at 37°C and 5% CO2 (Kim et al., 2019). This first step of antibiotic-protection assay was followed by another using a lower antibiotic concentration (40 μg/ml gentamicin and 10 μg/ml lysostaphin), in media supplemented with 5% fetal bovine serum (Gibco), and 1 μg/ml PI, a live cell-impermeant nucleic acid dye (Sigma). Plates were then incubated in the CLARIOstar Plus plate reader (BMG Labtech) set at 37°C and 5% CO2 throughout the infection (up to 20 hr post-infection). The fluorescence signal emitted by PI-positive cells was acquired every 6 min from each well (excitation at 535 nm, emission at 617 nm, using the spiral well scanning mode with 50 flashes per well). Non-infected control cells were permeabilised with 0.1% Triton X-100 to determine the maximum level PI uptake and HeLa cell death.

High-content imaging

Request a detailed protocol

The Operetta high-content microscope (PerkinElmer) was employed to accurately quantify and analyse intracellular persistence at the single-cell resolution. HeLa-CCL2 cells were seeded in Cell Carrier-96 black and optically clear bottom plates (PerkinElmer) to reach a density of 15,000 cells per well at the day of infection. HeLa cells were infected as described in the above section.

Post-infection, cells were washed twice with sterile PBS and fixed with 40 μl of freshly prepared 4% paraformaldehyde (Thermo Fisher Scientific) for 10 min. Fixed cells were further washed five times and stored at 4°C in PBS. Fixed cells were first permeabilised with 40 μl of 0.2% Triton X-100 for 3 min, washed thrice with PBS, and incubated 1 hr in 40 μl of blocking solution (PBS-bovine serum albumin [BSA] 3%). Bacteria were detected with polyclonal antibodies raised in rabbits against whole fixed cells of S. aureus USA300 strains, JE2::spa, BPH2919, and BPH3672 (WEHI antibody technology platform, https://www.wehi.edu.au/research/research-technologies/antibody-technologies). Sera were used at 1:1000, diluted in PBS-BSA 3%, Tween 0.05% for 5 hr at room temperature. Wells were then washed thrice with PBS and incubated 45 min with a secondary antibody (donkey anti-rabbit coupled to Alexa 488, 1:2000 dilution, Invitrogen) in PBS-BSA 3% containing 0.05% Tween-20 (Sigma) and 10% normal donkey serum (Abcam). Wells were washed thrice with PBS and incubated with Phalloidin-TRITC (1:4000) and DAPI (1:4000) (Sigma) in PBS for 15 min. Finally, wells were washed five times with PBS and covered with 200 μl of PBS. Plates were covered with aluminium foil and stored at 4°C until image acquisition on the Operetta microscope.

Confocal microscopy

Request a detailed protocol

HeLa-CCL2 grown on coverslips were infected using the same conditions described above. Coverslips were treated with PBS supplemented with 1% BSA and 0.2% Triton X-100 for 20 min at room temperature to permeabilise cells and incubated overnight in a blocking buffer (PBS supplemented with 1% BSA and 0.1% Tween-20). Coverslips were then probed 1 hr at room temperature with an anti-LAMP1 monoclonal antibody (1:250, clone H4A3 [mouse], Developmental Studies Hybridoma Bank) and 1:1000 polyclonal anti-S. aureus diluted in blocking buffer supplemented with 10% normal goat serum (Abcam). Coverslips were washed thrice with PBS then incubated overnight at 4°C with 1:2000 anti-rabbit (488), anti-mouse (647) secondary antibodies diluted in blocking buffer supplemented with 10% normal goat serum. Coverslips were then incubated 7 min with DAPI (1:5000), washed five times and mounted in Prolong Gold antifade (Thermo Fisher Scientific). Samples were imaged on the Zeiss LSM780 confocal microscope.

High-content imaging acquisition and analysis

Request a detailed protocol

Cells were analysed using the Operetta CLS high-content analysis system (PerkinElmer). For each well, images were acquired in a single plane at 11 non-overlapping fields of view (675 × 508 μm/1360 × 1024 pixels in size) using a ×20 PLAN long working distance objective (NA 0.45). DAPI fluorescence (HeLa cell nuclei) was imaged with the filter set: excitation = 360–400 nm, emission: 410–480 nm; 50 ms exposure. A488 fluorescence (S. aureus) was imaged with the filter set: excitation = 460–490 nm, emission = 500–550 nm; exposure = 200 ms. A594 fluorescence (HeLa actin stained by Phalloidin-TRITC) was imaged with the filter set ‘StdOrange1/Cy3’ filter set (excitation: 520–550 nm, emission: 560–630 nm; 0.5-s exposure).

Image processing and analysis were performed using the PhenoLOGIC machine-learning option in the Harmony software (PerkinElmer, v4.1). Nuclei were segmented from the DAPI channel using the ‘Find Nuclei algorithm’ (Method B, Area filter >40 μm2, Common Threshold of 0.4). Cells were segmented from the A594 channels using the Find Cells algorithm (Method C, Area filter >100 μm2, Common Threshold of 0.5). The A488 signals corresponding to S. aureus were further processed using a sliding parabola (curvature, 50 pixels) and Gaussian filter (filter width, 1 pixel) to remove noise and improve the signal-to-noise ratio. S. aureus were segmented by applying the Find Spots algorithm (Method A, Relative Spot Intensity 0.280, Splitting Coefficient 0.5).

Processing of PI fluorescence signals

Request a detailed protocol

For every 96-well plate, the PI uptake data for each well at each timepoint were standardised to the JE2 strain control using proportion of maximum scoring (POMS): (PI uptake − min(Pi uptake [JE2]))/range(PI uptake [JE2]). Experiments with less than two JE2 replicates available per plate were excluded from our analyses. Standardised data were used to fit a cubic smoothing spline (Little, 2013) using the R function smooth.spline(). Technical replicates within each plate were classified as outliers and excluded if >10% of their timepoint values differed by more than 1.5 times the interquartile range (Tukey method), between the fitted value and the mean for a given isolate. After excluding outlier replicates, fitted data were used to calculate the following PI uptake parameters: AUC, maximum PI uptake rate (rmax), peak PI uptake (μmax), time to maximum PI uptake rate (t(rmax)), time to peak PI uptake (t(μmax)), trough PI uptake, time of trough, and final PI uptake.

Dimensionality reduction of PI uptake data

Request a detailed protocol

Principal component analysis was performed using the ‘dudi.pca()’ function in the R package ‘adegenet’ and the randomForest package in R was used for fitting an unsupervised random forest model. We used the similarity matrix generated by the model to define similarity cluster of PI uptake. We used the output of the RF model to calculate the importance of each PI uptake parameter defined as mean decrease in Gini index and mean decrease in accuracy (Breiman, 2001).

PI uptake GWAS

Request a detailed protocol

We transformed the mean PI uptake AUC data using the automated normalisation package bestNormalize in R. A GWAS using the normalised AUC data and the 158,169 core-genome variants (all positions where at least 90% of the strains had at least 10 reads coverage, see above) obtained after mapping isolates reads to reference genome BPH2947. All 387 isolates were included in the GWAS, including serial isolates from the same patient, as this outcome is bacterial and not expected to be affected by the host (Tonkin-Hill et al., 2022). To correct for the population structure, we used the factored spectrally transformed linear mixed models (FaST-LMM) implemented in pyseer v1.3.6 (Lees et al., 2018). Random effects in Fast-LMM were computed from a kinship matrix based on the core-genome SNPs generated by Gemma v0.98.1 (Zhou and Stephens, 2012). The Bonferroni method was used to correct p values for multiple testing. We performed the GWAS using single variants and the gene-burden test implemented in pyseer. We excluded synonymous mutations for single variants and gene-burden GWAS. For the single variants GWAS, only mutations with a minimum allele fraction (MAF) of 0.01 (as suggested in the pyseer documentation [https://pyseer.readthedocs.io/en/master/index.html]) were kept and at least two independent acquisitions across the phylogeny. For the gene-burden GWAS, we aggregated rare mutations (MAF <0.01), with 3072 coding sequences of the reference and modelled the association of the PI AUC with the presence/absence of any protein-altering mutations in each BPH2947 gene. For consistent annotation of mutations, we identified BPH2947 genes homologs using BLASTP and annotated FPR3757 genes using aureowiki (Fuchs et al., 2018) and Microbesonline (Alm et al., 2005). The GWAS analysis was run using a customised in-house pipeline (https://github.com/stefanogg/CLOGEN).

GWAS power calculation

Request a detailed protocol

We used BacGWASim (Saber and Shapiro, 2020) to simulate a genotype matrix with 1000 and 150,000 variants for a range of sample size between 300 and 9600. We set mutation rate at 0.06 and recombination rate at 0.01. We simulated 10 phenotype replicates with 16 true causal variants and a heritability of 0.5. We then run pyseer on each replicated simulation using the same pipeline described above. As in Denamur et al., 2022, power was calculated as the fraction of identified variants based on a Bonferroni corrected p value below 0.05.

Determination of genetic pairs with contrasting PI uptake

Request a detailed protocol

Isolate pairs for the convergent evolution analysis were identified by screening pairs separated by less than 200 mutations distance for statistically significant differences in the PI uptake AUC (Mann–Whitney test), wherein an isolate causing low PI uptake (isolate-2) was compared to a reference isolate causing higher PI uptake (isolate-1). The genetic distance between closely related isolates was calculated using Snippy and is based on the number of variants identified when mapping the reads of isolate-2 on the draft assembly of isolate-1 (Higgs et al., 2022). To avoid biases related to assembly errors and uneven reads coverage between the two isolates, variant calls were filtered as previously described (https://github.com/stefanogg/staph_adaptation_paper; Giulieri et al., 2022a; Giulieri, 2022). Non-redundant and phylogenetically independent genetic pairs were identified by manual inspection of the phylogenetic tree.

Genetic pairs analysis

Request a detailed protocol

Mutations identified in genetic pairs and filtered as described above were further characterised using a multilayered annotation strategy as previously described (Giulieri et al., 2022a). Firstly, mutated coding regions (amino acid sequences) across draft genomes were clustered using CD-HIT. We then used BLASTP to identify homologs of each cluster within the S. aureus USA300 FPR3757 reference genome that was annotated using the AureoWiki repository (Fuchs et al., 2018), with a 90% identity and 50% coverage threshold. As genetic pairs were phylogenetically independent and non-redundant, emergence of the same mutation or mutations in the same locus in multiple pairs indicated convergent evolution and was suggestive of positive selection. Based on this, we ranked USA300 FPR3757 homologs according to the number of pairs with mutations.

Comparison of PI uptake AUC between phylogenetically independent pairs with or without homoplasic mutations

Request a detailed protocol

A set of phylogenetically independent pairs was sampled from the tree using the R package ggtree. The impact of the presence/absence of homoplasic mutations on the PI uptake AUC difference within pairs was modelled using linear regression. The sampling and regression were repeated in 100 replicates. As this approach is not suitable for p-value statistical significance, the distribution of the regression coefficient beta of the replicates was used to infer the statistical support for the associations of the homoplasic mutations with the PI uptake phenotype.

Code and data availability

Request a detailed protocol

Scripts to process PI uptake data and to perform genomic analyses are available on github at https://github.com/stefanogg/InToxSa (copy archived at Giulieri and Guerillot, 2023). The code for genomic analyses is available on https://github.com/stefanogg/CLOGEN (copy archived at Giulieri, 2023) (GWAS analysis), and on https://github.com/stefanogg/staph_adaptation_paper (Giulieri, 2022) (comparative genomics of genetic pairs).

Whole-genome sequences of the 387 clinical strains are available in the European Nucleotide Archive under Bioproject accession number PRJEB27932.

Aureusimine B identification

Request a detailed protocol

Bacterial extracts were isolated from 30 ml cultures grown in TSB at 37°C overnight under agitation. Bacterial cells were pelleted by centrifugation at 4000 × g during 30 min and the culture supernatants were sterilised by passage through a 0.22 mM filter. For each strain, 10 ml of supernatant were added to an equal volume of ethyl acetate in glass tubes, vortexed and allowed to extract at room temperature overnight. Ethyl acetate extracts were dried in vacuo. Dried ethyl acetate extracts were resuspended in 100 ml methanol and 2 ml of each sample was analysed by HPLC using the Shimadzu Prominence HPLC system coupled to a SPD-M20A diode array detector. The column oven (CTO-20A) was set to 40°C and aureusimine B was separated on Kinetex C18, 75 × 3 mm, 2.6 μm column (Phenomenex). Purified aureusimine B was used as reference standard (Bioaustralis Fine Chemicals). All used chemicals were of analytical grade.

Samples were run with water, 0.1% trifluoroacetic acid (solvent A) and acetonitrile (solvent B). The gradient elution was performed on the HPLC at a flow rate of 0.5 ml/min as follows: 10% B for 3.5 min, 10–100% B over 12.5 min, 100–10% B over 1 min, then 10% B for 7 min (total run time, 24 min).

Data availability

-The genome sequences of Staphylococcus aureus clinical isolates are deposited on the European Nucleotide Archive under Bioproject accession number PRJEB27932.- All data generated during this study are included in the supplementary files of the manuscript. The generated codes are publicly available on: https://github.com/stefanogg/InToxSa (copy archived at Giulieri and Guerillot, 2023).

The following previously published data sets were used
    1. Giulieri SG
    (2018) European Nucleotide archive
    ID PRJEB27932. Genomic Exploration of Sequential Clinical Isolates Reveals a Distinctive Molecular Signature of Persistent Staphylococcus aureus Bacteraemia.
    https://www.ebi.ac.uk/ena/browser/view/PRJEB27932

References

  53. Book
    1. Little TD
    (2013)
    Longitudinal Structural Equation Modeling
    The Guilford Press.
    1. Murray CJL
    2. Ikuta KS
    3. Sharara F
    4. Swetschinski L
    5. Robles Aguilar G
    6. Gray A
    7. Han C
    8. Bisignano C
    9. Rao P
    10. Wool E
    11. Johnson SC
    12. Browne AJ
    13. Chipeta MG
    14. Fell F
    15. Hackett S
    16. Haines-Woodhouse G
    17. Kashef Hamadani BH
    18. Kumaran EAP
    19. McManigal B
    20. Achalapong S
    21. Agarwal R
    22. Akech S
    23. Albertson S
    24. Amuasi J
    25. Andrews J
    26. Aravkin A
    27. Ashley E
    28. Babin FX
    29. Bailey F
    30. Baker S
    31. Basnyat B
    32. Bekker A
    33. Bender R
    34. Berkley JA
    35. Bethou A
    36. Bielicki J
    37. Boonkasidecha S
    38. Bukosia J
    39. Carvalheiro C
    40. Castañeda-Orjuela C
    41. Chansamouth V
    42. Chaurasia S
    43. Chiurchiù S
    44. Chowdhury F
    45. Clotaire Donatien R
    46. Cook AJ
    47. Cooper B
    48. Cressey TR
    49. Criollo-Mora E
    50. Cunningham M
    51. Darboe S
    52. Day NPJ
    53. De Luca M
    54. Dokova K
    55. Dramowski A
    56. Dunachie SJ
    57. Duong Bich T
    58. Eckmanns T
    59. Eibach D
    60. Emami A
    61. Feasey N
    62. Fisher-Pearson N
    63. Forrest K
    64. Garcia C
    65. Garrett D
    66. Gastmeier P
    67. Giref AZ
    68. Greer RC
    69. Gupta V
    70. Haller S
    71. Haselbeck A
    72. Hay SI
    73. Holm M
    74. Hopkins S
    75. Hsia Y
    76. Iregbu KC
    77. Jacobs J
    78. Jarovsky D
    79. Javanmardi F
    80. Jenney AWJ
    81. Khorana M
    82. Khusuwan S
    83. Kissoon N
    84. Kobeissi E
    85. Kostyanev T
    86. Krapp F
    87. Krumkamp R
    88. Kumar A
    89. Kyu HH
    90. Lim C
    91. Lim K
    92. Limmathurotsakul D
    93. Loftus MJ
    94. Lunn M
    95. Ma J
    96. Manoharan A
    97. Marks F
    98. May J
    99. Mayxay M
    100. Mturi N
    101. Munera-Huertas T
    102. Musicha P
    103. Musila LA
    104. Mussi-Pinhata MM
    105. Naidu RN
    106. Nakamura T
    107. Nanavati R
    108. Nangia S
    109. Newton P
    110. Ngoun C
    111. Novotney A
    112. Nwakanma D
    113. Obiero CW
    114. Ochoa TJ
    115. Olivas-Martinez A
    116. Olliaro P
    117. Ooko E
    118. Ortiz-Brizuela E
    119. Ounchanum P
    120. Pak GD
    121. Paredes JL
    122. Peleg AY
    123. Perrone C
    124. Phe T
    125. Phommasone K
    126. Plakkal N
    127. Ponce-de-Leon A
    128. Raad M
    129. Ramdin T
    130. Rattanavong S
    131. Riddell A
    132. Roberts T
    133. Robotham JV
    134. Roca A
    135. Rosenthal VD
    136. Rudd KE
    137. Russell N
    138. Sader HS
    139. Saengchan W
    140. Schnall J
    141. Scott JAG
    142. Seekaew S
    143. Sharland M
    144. Shivamallappa M
    145. Sifuentes-Osornio J
    146. Simpson AJ
    147. Steenkeste N
    148. Stewardson AJ
    149. Stoeva T
    150. Tasak N
    151. Thaiprakong A
    152. Thwaites G
    153. Tigoi C
    154. Turner C
    155. Turner P
    156. van Doorn HR
    157. Velaphi S
    158. Vongpradith A
    159. Vongsouvath M
    160. Vu H
    161. Walsh T
    162. Walson JL
    163. Waner S
    164. Wangrangsimakul T
    165. Wannapinij P
    166. Wozniak T
    167. Young Sharma T
    168. Yu KC
    169. Zheng P
    170. Sartorius B
    171. Lopez AD
    172. Stergachis A
    173. Moore C
    174. Dolecek C
    175. Naghavi M
    (2022) Global burden of bacterial antimicrobial resistance in 2019: A systematic analysis
    The Lancet 399:629–655.
    https://doi.org/10.1016/S0140-6736(21)02724-0
    1. Trubetskoy V
    2. Pardiñas AF
    3. Qi T
    4. Panagiotaropoulou G
    5. Awasthi S
    6. Bigdeli TB
    7. Bryois J
    8. Chen C-Y
    9. Dennison CA
    10. Hall LS
    11. Lam M
    12. Watanabe K
    13. Frei O
    14. Ge T
    15. Harwood JC
    16. Koopmans F
    17. Magnusson S
    18. Richards AL
    19. Sidorenko J
    20. Wu Y
    21. Zeng J
    22. Grove J
    23. Kim M
    24. Li Z
    25. Voloudakis G
    26. Zhang W
    27. Adams M
    28. Agartz I
    29. Atkinson EG
    30. Agerbo E
    31. Al Eissa M
    32. Albus M
    33. Alexander M
    34. Alizadeh BZ
    35. Alptekin K
    36. Als TD
    37. Amin F
    38. Arolt V
    39. Arrojo M
    40. Athanasiu L
    41. Azevedo MH
    42. Bacanu SA
    43. Bass NJ
    44. Begemann M
    45. Belliveau RA
    46. Bene J
    47. Benyamin B
    48. Bergen SE
    49. Blasi G
    50. Bobes J
    51. Bonassi S
    52. Braun A
    53. Bressan RA
    54. Bromet EJ
    55. Bruggeman R
    56. Buckley PF
    57. Buckner RL
    58. Bybjerg-Grauholm J
    59. Cahn W
    60. Cairns MJ
    61. Calkins ME
    62. Carr VJ
    63. Castle D
    64. Catts SV
    65. Chambert KD
    66. Chan RCK
    67. Chaumette B
    68. Cheng W
    69. Cheung EFC
    70. Chong SA
    71. Cohen D
    72. Consoli A
    73. Cordeiro Q
    74. Costas J
    75. Curtis C
    76. Davidson M
    77. Davis KL
    78. de Haan L
    79. Degenhardt F
    80. DeLisi LE
    81. Demontis D
    82. Dickerson F
    83. Dikeos D
    84. Dinan T
    85. Djurovic S
    86. Duan J
    87. Ducci G
    88. Dudbridge F
    89. Eriksson JG
    90. Fañanás L
    91. Faraone SV
    92. Fiorentino A
    93. Forstner A
    94. Frank J
    95. Freimer NB
    96. Fromer M
    97. Frustaci A
    98. Gadelha A
    99. Genovese G
    100. Gershon ES
    101. Giannitelli M
    102. Giegling I
    103. Giusti-Rodríguez P
    104. Godard S
    105. Goldstein JI
    106. González Peñas J
    107. González-Pinto A
    108. Gopal S
    109. Gratten J
    110. Green MF
    111. Greenwood TA
    112. Guillin O
    113. Gülöksüz S
    114. Gur RE
    115. Gur RC
    116. Gutiérrez B
    117. Hahn E
    118. Hakonarson H
    119. Haroutunian V
    120. Hartmann AM
    121. Harvey C
    122. Hayward C
    123. Henskens FA
    124. Herms S
    125. Hoffmann P
    126. Howrigan DP
    127. Ikeda M
    128. Iyegbe C
    129. Joa I
    130. Julià A
    131. Kähler AK
    132. Kam-Thong T
    133. Kamatani Y
    134. Karachanak-Yankova S
    135. Kebir O
    136. Keller MC
    137. Kelly BJ
    138. Khrunin A
    139. Kim S-W
    140. Klovins J
    141. Kondratiev N
    142. Konte B
    143. Kraft J
    144. Kubo M
    145. Kučinskas V
    146. Kučinskiene ZA
    147. Kusumawardhani A
    148. Kuzelova-Ptackova H
    149. Landi S
    150. Lazzeroni LC
    151. Lee PH
    152. Legge SE
    153. Lehrer DS
    154. Lencer R
    155. Lerer B
    156. Li M
    157. Lieberman J
    158. Light GA
    159. Limborska S
    160. Liu C-M
    161. Lönnqvist J
    162. Loughland CM
    163. Lubinski J
    164. Luykx JJ
    165. Lynham A
    166. Macek M
    167. Mackinnon A
    168. Magnusson PKE
    169. Maher BS
    170. Maier W
    171. Malaspina D
    172. Mallet J
    173. Marder SR
    174. Marsal S
    175. Martin AR
    176. Martorell L
    177. Mattheisen M
    178. McCarley RW
    179. McDonald C
    180. McGrath JJ
    181. Medeiros H
    182. Meier S
    183. Melegh B
    184. Melle I
    185. Mesholam-Gately RI
    186. Metspalu A
    187. Michie PT
    188. Milani L
    189. Milanova V
    190. Mitjans M
    191. Molden E
    192. Molina E
    193. Molto MD
    194. Mondelli V
    195. Moreno C
    196. Morley CP
    197. Muntané G
    198. Murphy KC
    199. Myin-Germeys I
    200. Nenadić I
    201. Nestadt G
    202. Nikitina-Zake L
    203. Noto C
    204. Nuechterlein KH
    205. O’Brien NL
    206. O’Neill FA
    207. Oh S-Y
    208. Olincy A
    209. Ota VK
    210. Pantelis C
    211. Papadimitriou GN
    212. Parellada M
    213. Paunio T
    214. Pellegrino R
    215. Periyasamy S
    216. Perkins DO
    217. Pfuhlmann B
    218. Pietiläinen O
    219. Pimm J
    220. Porteous D
    221. Powell J
    222. Quattrone D
    223. Quested D
    224. Radant AD
    225. Rampino A
    226. Rapaport MH
    227. Rautanen A
    228. Reichenberg A
    229. Roe C
    230. Roffman JL
    231. Roth J
    232. Rothermundt M
    233. Rutten BPF
    234. Saker-Delye S
    235. Salomaa V
    236. Sanjuan J
    237. Santoro ML
    238. Savitz A
    239. Schall U
    240. Scott RJ
    241. Seidman LJ
    242. Sharp SI
    243. Shi J
    244. Siever LJ
    245. Sigurdsson E
    246. Sim K
    247. Skarabis N
    248. Slominsky P
    249. So H-C
    250. Sobell JL
    251. Söderman E
    252. Stain HJ
    253. Steen NE
    254. Steixner-Kumar AA
    255. Stögmann E
    256. Stone WS
    257. Straub RE
    258. Streit F
    259. Strengman E
    260. Stroup TS
    261. Subramaniam M
    262. Sugar CA
    263. Suvisaari J
    264. Svrakic DM
    265. Swerdlow NR
    266. Szatkiewicz JP
    267. Ta TMT
    268. Takahashi A
    269. Terao C
    270. Thibaut F
    271. Toncheva D
    272. Tooney PA
    273. Torretta S
    274. Tosato S
    275. Tura GB
    276. Turetsky BI
    277. Üçok A
    278. Vaaler A
    279. van Amelsvoort T
    280. van Winkel R
    281. Veijola J
    282. Waddington J
    283. Walter H
    284. Waterreus A
    285. Webb BT
    286. Weiser M
    287. Williams NM
    288. Witt SH
    289. Wormley BK
    290. Wu JQ
    291. Xu Z
    292. Yolken R
    293. Zai CC
    294. Zhou W
    295. Zhu F
    296. Zimprich F
    297. Atbaşoğlu EC
    298. Ayub M
    299. Benner C
    300. Bertolino A
    301. Black DW
    302. Bray NJ
    303. Breen G
    304. Buccola NG
    305. Byerley WF
    306. Chen WJ
    307. Cloninger CR
    308. Crespo-Facorro B
    309. Donohoe G
    310. Freedman R
    311. Galletly C
    312. Gandal MJ
    313. Gennarelli M
    314. Hougaard DM
    315. Hwu H-G
    316. Jablensky AV
    317. McCarroll SA
    318. Moran JL
    319. Mors O
    320. Mortensen PB
    321. Müller-Myhsok B
    322. Neil AL
    323. Nordentoft M
    324. Pato MT
    325. Petryshen TL
    326. Pirinen M
    327. Pulver AE
    328. Schulze TG
    329. Silverman JM
    330. Smoller JW
    331. Stahl EA
    332. Tsuang DW
    333. Vilella E
    334. Wang S-H
    335. Xu S
    336. Adolfsson R
    337. Arango C
    338. Baune BT
    339. Belangero SI
    340. Børglum AD
    341. Braff D
    342. Bramon E
    343. Buxbaum JD
    344. Campion D
    345. Cervilla JA
    346. Cichon S
    347. Collier DA
    348. Corvin A
    349. Curtis D
    350. Forti MD
    351. Domenici E
    352. Ehrenreich H
    353. Escott-Price V
    354. Esko T
    355. Fanous AH
    356. Gareeva A
    357. Gawlik M
    358. Gejman PV
    359. Gill M
    360. Glatt SJ
    361. Golimbet V
    362. Hong KS
    363. Hultman CM
    364. Hyman SE
    365. Iwata N
    366. Jönsson EG
    367. Kahn RS
    368. Kennedy JL
    369. Khusnutdinova E
    370. Kirov G
    371. Knowles JA
    372. Krebs M-O
    373. Laurent-Levinson C
    374. Lee J
    375. Lencz T
    376. Levinson DF
    377. Li QS
    378. Liu J
    379. Malhotra AK
    380. Malhotra D
    381. McIntosh A
    382. McQuillin A
    383. Menezes PR
    384. Morgan VA
    385. Morris DW
    386. Mowry BJ
    387. Murray RM
    388. Nimgaonkar V
    389. Nöthen MM
    390. Ophoff RA
    391. Paciga SA
    392. Palotie A
    393. Pato CN
    394. Qin S
    395. Rietschel M
    396. Riley BP
    397. Rivera M
    398. Rujescu D
    399. Saka MC
    400. Sanders AR
    401. Schwab SG
    402. Serretti A
    403. Sham PC
    404. Shi Y
    405. St Clair D
    406. Stefánsson H
    407. Stefansson K
    408. Tsuang MT
    409. van Os J
    410. Vawter MP
    411. Weinberger DR
    412. Werge T
    413. Wildenauer DB
    414. Yu X
    415. Yue W
    416. Holmans PA
    417. Pocklington AJ
    418. Roussos P
    419. Vassos E
    420. Verhage M
    421. Visscher PM
    422. Yang J
    423. Posthuma D
    424. Andreassen OA
    425. Kendler KS
    426. Owen MJ
    427. Wray NR
    428. Daly MJ
    429. Huang H
    430. Neale BM
    431. Sullivan PF
    432. Ripke S
    433. Walters JTR
    434. O’Donovan MC
    435. Indonesia Schizophrenia Consortium
    436. PsychENCODE
    437. Psychosis Endophenotypes International Consortium
    438. SynGO Consortium
    439. Schizophrenia Working Group of the Psychiatric Genomics Consortium
    (2022) Mapping Genomic Loci Implicates genes and synaptic biology in schizophrenia
    Nature 604:502–508.
    https://doi.org/10.1038/s41586-022-04434-5

Decision letter

  1. Arturo Casadevall
    Senior and Reviewing Editor; Johns Hopkins Bloomberg School of Public Health, United States
  2. Matthew Culyba
    Reviewer; University of Pittsburgh, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "A high-throughput cytotoxicity screening platform reveals agr-independent mutations in bacteraemia-associated Staphylococcus aureus that promote intracellular persistence" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Arturo Casadevall as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Matthew Culyba (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

Reviewer #1 (Recommendations for the authors):

Consider depositing the long read sequence of the HeLa strain, or at least chromosome mapping/ structure variant mapping. These could be important when other investigators attempt to reproduce results using their version of the cell line.

Line 145 – What specifically is "very low intra assay variation"?

Line 283 – Unclear if these are the only 28 pairs that met the criterion, or if there were others but that weren't selected.

Line 290 – Were homoplasic mutations in these loci found between other pairs of strains (ie where there were no significant differences)

Line 324-328 – This correlation between π AUC and intracellular bacterial titer should be shown as a scatterplot for clarity

Reviewer #2 (Recommendations for the authors):

1. The InToxSa assay was devised to reflect the toxicity of intracellular S. aureus, however, it is not made clear in the manuscript that the intracellular nature of this assay is actually an important feature that distinguishes it from what would be observed by testing the same strains in the equivalent extracellular toxicity assay. For example, the interpretation of the InToxSa versus Tryptan blue exclusion assay data was that InToxSa was more sensitive. However, two key variables differed in this comparison: (1) extracellular toxicity vs intracellular toxicity and (2) cell type used. Is InToxSa more sensitive in this comparison because of the intracellular nature of the assay, or because it employs HeLa cells instead of THP1 macrophages? Maybe HeLa cells are just more sensitive regardless? What would happen if you repeated the Tryptan blue exclusion assay with HeLa cells?

2. High throughput screening assays often report a 'reconfirmation rate' for the assay. This helps address day-to-day assay signal variation from the actual experimental samples (instead of controls). If you just rerun the 56 isolates that make up the 28 pairs from the convergence analysis, how well do the data replicate in terms of the difference in π uptake? I'm wondering how much of the low validation/confirmation rate with the NTLM and allelic exchange strains could be due to day-to-day assay noise in the primary screen. Maybe some of these 'hits' were simply false positives?

3. The analysis focuses on negative π difference values. What is the interpretation of a positive value for π differences? Higher toxicity? Were any of these identified by GWAS or the convergence analysis? Are these of biological interest?

4. I'm wondering how the authors envision others deploying this system to screen clinical isolates since the 'gene burden GWAS' only found agr, which is a well-described virulence locus. How many isolates would have to be screened and analyzed in this manner to find additional biologically relevant mutant alleles? The statistical power of this approach is related to the number of isolates (N), their genetic relatedness, the magnitude of the effect size (in this case, π uptake signal), and the gene burden. A power analysis could be a useful way how to think about deploying this assay using the 'phenomics' approach. Would it be possible to model these parameters in the gene burden GWAS, and provide an estimate of the number of isolates (N) that would need to be screened using the assay to find mutant alleles as a function of their effect size (i.e. PI-uptake), gene burden, and the genetic relatedness of the sample population? Perhaps you will find that the assay is best deployed in settings where the strains being screened are all very closely genetically related (e.g. serial isolates from the same patient or transmission outbreak). Outside of these situations, I wonder if N might be too large for this to be a feasible approach to finding novel genes involved in intracellular toxicity. Modeling this as a power analysis would describe this quantitatively and likely point to the best types of clinical strain samples to pursue given the practical constraint that screening more than a few thousand isolates in this system would not be feasible. It seems you already took steps to enrich for serial isolates from the same patient and then the convergence analysis was even further enriched for these samples. These are clearly critical steps to get good information output. Modeling the statistical power of this would be a nice complement I think.

5. All of the 'hits' from the gene burden GWAS and convergence analyses were attempted to be validated using the corresponding strains from the NTML. This begs the question: why not just screen the entire NTML as a starting point? It seems this approach would have not only offered better validation for the utility of the assay but also provided maximal statistical power for genotype-phenotype correlations given the isogenic background within this strain library. With this in mind, can you elaborate further on your choice of approach?

6. The manuscript describes using HeLa cells as an in vitro model system for professional phagocytic cells that are thought to be important for S. aureus clearance in vivo. The biological relevance of the HeLa cell system would be further supported by also studying some mutants in a system that employs professional phagocytes to study intracellular S. aureus persistence. For example, it would be of interest to know if the ausA mutants have a persistence phenotype in a macrophage or not.

7. The statistical 'gene burden GWAS' procedure only identified agr as having a significant P-value. In contrast, a manual 'convergence' analysis and hunt for homoplastic mutations for π uptake seemed to identify additional plausible candidates, but this manual procedure was not supported by statistical analysis. Can you explain why the convergence analysis would find loci that the 'gene burden GWAS' failed to identify? Can the mutations identified by the convergence analysis be further supported within a formal statistical framework? In the 'convergence' analysis, you only consider pairs with <200 SNP differences, which effectively increases the chance (i.e. statistical power) that any given mutation could be responsible for the observed π uptake effect size. What would happen if you re-ran the 'gene burden GWAS' using only a subset of the more genetically related strain pairs? Would that increase sensitivity enough to identify the homoplastic mutations you found in the 'convergence' analysis?

https://doi.org/10.7554/eLife.84778.sa1

Author response

Essential revisions:

Reviewer #1 (Recommendations for the authors):

Consider depositing the long read sequence of the HeLa strain, or at least chromosome mapping/ structure variant mapping. These could be important when other investigators attempt to reproduce results using their version of the cell line.

We refer interested investigators to the American Type Culture Collection (ATCC) from where the HeLa CCL2 were obtained and refer them to the following manuscript describing the phenotypes obtained during bacterial infection of various HeLa cell types (Tang L. et al., DOI: 10.1038/s41592-019-0375-1). The haplotype-resolved genome and epigenome of the aneuploid HeLa CCL2 cell line used in this study can be found in the following publication: DOI: 10.1038/nature12064.

Line 145 – What specifically is "very low intra assay variation"?

We substantiate “very low intra assay variation” by providing the coefficient of variation (CoV) and standard deviation (Std. Dev.) values in Table 1. The CoV values are lower than those recommended in the ISO standards for diagnostic tests. We now also provide standard deviation values (Std. Dev.) in Table 1 as requested.

Line 283 – Unclear if these are the only 28 pairs that met the criterion, or if there were others but that weren't selected.

Pairs with contrasting π uptake were selected according to 3 criteria:

1) < 200 mutations distance between isolate 1 and isolate 2,

2) significantly lower π uptake in isolate 2 compared to isolate 1 (Mann-Whitney test), and

3) phylogenetic independence (i.e pairs belonging to non-overlapping clades) (Methods lines 691-701).

While the total number of pairs satisfying criteria 1 and 2 is high, the need for phylogenetic independence dramatically restricts the number of eligible pairs. Without this independence filter, the occurrence of homoplasic mutations would be unduly inflated according to the size of the clade and meaningless for detecting and comparing the strength of the genetic association signals (i.e representing redundant genetic comparisons). Previous studies have successfully used the same approach (Pidot et al., 2018, DOI: 10.1126/scitranslmed.aar6115) to address this issue.

Line 290 – Were homoplasic mutations in these loci found between other pairs of strains (ie where there were no significant differences)

This is an interesting question that we have now addressed by comparing the difference in π uptake between independent genetic pairs that carried mutations in the most convergent genes and those that didn’t. Because multiple combinations of phylogenetically independent pairs are possible, we repeated the analysis in 100 replicates and calculated the median of absolute difference in π uptake for each group and fitted 100 replicated linear regressions to assess the strength of the association between the presence/absence of homoplasic mutations and the π uptake difference.

For six genes (agrA, glcA, ribA, fmtB, sbnF, tarL’), this supplementary analysis provided full support of the association with the cytotoxicity phenotype (100% of the linear regression replicates had a coefficient above 0). Similarly, 84% replicates supported the association between mutations occurring in the gene ausA and the π uptake. This new analysis is reported in the results (lines 319-339), methods (lines 780-787) and illustrated in Figure 5 —figure supplementary 2.

Together, our analyses show that most genes with homoplasic mutations in phenotypically divergent pairs were not mutated in phenotypically similar pairs, with the gene agrC being the only exception. While the agrC transposon mutant has reduced cytotoxicity, confirming the expected role of agrC loss-of-function mutations in reduced cytotoxicity, it remains possible that the mutations identified in our clinical strains could be compensatory and therefore not directly associated with the phenotype (see the following references reporting such phenomena, DOI: 10.1128/IAI.003331-18 and DOI: 10.7554/eLife.77195) or, that the null effect on cytotoxicity of aggregated agrC mutations could be due to the presence of both cytotoxicity-increasing and decreasing mutations within the same isolate. We have now included these supporting references in the text (Lines 444-447).

Line 324-328 – This correlation between π AUC and intracellular bacterial titer should be shown as a scatterplot for clarity

We now provide two scatterplots showing the statistically significant inverse Pearson correlations of: (1) π AUC with the percentage of infected cells; and (2) π AUC with the number of S. aureus per infected cell at 24 hours post infection. These two plots are cited in the main text and available in Figure 5 –Supplementary figure 5.

Reviewer #2 (Recommendations for the authors):

1. The InToxSa assay was devised to reflect the toxicity of intracellular S. aureus, however, it is not made clear in the manuscript that the intracellular nature of this assay is actually an important feature that distinguishes it from what would be observed by testing the same strains in the equivalent extracellular toxicity assay. For example, the interpretation of the InToxSa versus Tryptan blue exclusion assay data was that InToxSa was more sensitive. However, two key variables differed in this comparison: (1) extracellular toxicity vs intracellular toxicity and (2) cell type used. Is InToxSa more sensitive in this comparison because of the intracellular nature of the assay, or because it employs HeLa cells instead of THP1 macrophages? Maybe HeLa cells are just more sensitive regardless? What would happen if you repeated the Tryptan blue exclusion assay with HeLa cells?

We are aware of the potential limitations of both Trypan blue exclusion and InToxSa methods, and these have been presented in the introduction (Lines 65-97) and in the discussion of the manuscript (Lines 403-428). In Figure 3, the InToxSa assay is more sensitive than the trypan blue exclusion assay in the sense that InToxSa detects more differences between isolates and thus increases the power of our analytical framework to identify genetic differences between these isolates. It would be possible to measure the sensitivity of HeLa cells to extracellular, secreted staphylococcal toxins using Trypan Blue, but such experiments are outside the aims of the current study. Furthermore, such assays are low-throughput in nature and are not continuous (end-point only), requiring extensive trial-and-error experimentation to establish consistent timepoints that represent thresholds of cytotoxicity between isolates.

2. High throughput screening assays often report a 'reconfirmation rate' for the assay. This helps address day-to-day assay signal variation from the actual experimental samples (instead of controls). If you just rerun the 56 isolates that make up the 28 pairs from the convergence analysis, how well do the data replicate in terms of the difference in π uptake? I'm wondering how much of the low validation/confirmation rate with the NTLM and allelic exchange strains could be due to day-to-day assay noise in the primary screen. Maybe some of these 'hits' were simply false positives?

The 56 isolates were assessed on at least two different plates on different days, with three biological replicates on each plate (2 to 5 plates/days per isolate; 6 to 15 replicates in total). Day-to-day assay noise is not responsible for the low reproduction of phenotype in mutants. We believe that this low confirmation rate is most likely due to epistatic interactions where the functional impact of the detected genomic signatures is conditioned upon the presence of other interacting gene(s)/mutation(s) in a specific genetic background. Unfortunately, the Nebraska transposon mutant library has been made in a ST8 background but none of the convergent mutations identified in our study have emerged in this genetic background.

3. The analysis focuses on negative π difference values. What is the interpretation of a positive value for π differences? Higher toxicity? Were any of these identified by GWAS or the convergence analysis? Are these of biological interest?

We focussed on pairs with loss of cytotoxicity to be consistent with the hypothesis that loss of toxicity (rather than gain) is an adaptive phenotype during infection. This is supported by our gene-burden GWAS, where 10/10 most significant hits are associated with loss of cytotoxicity (mean normalised π AUC decrease = –0.8, Supplementary file 5). However, we agree that there could also be an interest in investigating gain of cytotoxicity, as this phenotype might be advantageous for the bacteria in certain clinical settings (e.g. abscess formation). To further explore this, we selected 7 genetic pairs that had a significant π uptake increase from iso1 to iso2 and were phylogenetically independent from each other and from the previously selected pairs. We repeated the convergence analysis and found that adding these pairs further increased the number of independent (homoplasic) ausA and fmtB mutations. This reinforces the importance of ausA as an agr-independent cytotoxicity locus (Result: lines 320-329; Figure 5- Supplementary Figure 5).

4. I'm wondering how the authors envision others deploying this system to screen clinical isolates since the 'gene burden GWAS' only found agr, which is a well-described virulence locus. How many isolates would have to be screened and analyzed in this manner to find additional biologically relevant mutant alleles? The statistical power of this approach is related to the number of isolates (N), their genetic relatedness, the magnitude of the effect size (in this case, π uptake signal), and the gene burden. A power analysis could be a useful way how to think about deploying this assay using the 'phenomics' approach.

To address this issue, we have now used the bacterial GWAS power calculation approach suggested in Denamur et al., PLoS Genetics 2022 (this reference is included in the revised manuscript and in the methods) and used the simulator BacGWASim to simulate genome datasets of increasing size (from 300 samples to 9600) for a genome length of 1 million bp, a mutation rate of 0.06 and recombination rate of 0.01. For each dataset 16 causal variants were distributed in 10 genomes. We then ran pyseer using the same approach as in our GWAS analysis and calculated power as the proportion of causal variants that were above the Bonferroni-corrected threshold. We found that the power of the mutations GWAS was overall low, with 0% mutations called for sample size < 2,400 and a 20% power for the largest sample. As commented by the reviewer and consistent with the BacGWASim paper (Saber, Morteza and Shapiro, Microb Genom 2020), bacterial GWAS has limited power when applied to low frequency recombining pathogen like S. aureus. This supports our strategy to use convergent evolution analysis in our exploration of mutations associated with cytotoxicity. Despite not being specifically designed for convergence evolution analysis, our simulations for a low variants GWAS (n=1,000) further confirms our strategy (see Figure 4 – Supplementary Figure 4).

Note however, that this power calculation approach has inherent limitations. Firstly, the magnitude of the functional impact of mutations is not accounted for. Secondly, it relies heavily on rough estimates of core genome sizes, recombination and mutation rates that would be dataset dependent. Finally, it doesn’t simulate rare variants and the impact of aggregating them in a gene burden test. We agree with the reviewer that a more sophisticated bacterial GWAS simulation framework would provide more accurate information to plan bacterial phenotype-genotype studies. However, we feel that a more detailed power analysis is outside the aims of the current study and thus should be the topic of a specific statistical genomics methods paper.

5. Would it be possible to model these parameters in the gene burden GWAS, and provide an estimate of the number of isolates (N) that would need to be screened using the assay to find mutant alleles as a function of their effect size (i.e. PI-uptake), gene burden, and the genetic relatedness of the sample population?

We agree with the reviewer that modelling these parameters in the gene burden GWAS would be a useful and complementary approach to model the statistical power of such studies. However, such analysis is outside the aims of the current study and would represent a standalone statistical genomics methods paper.

6. All of the 'hits' from the gene burden GWAS and convergence analyses were attempted to be validated using the corresponding strains from the NTML. This begs the question: why not just screen the entire NTML as a starting point? It seems this approach would have not only offered better validation for the utility of the assay but also provided maximal statistical power for genotype-phenotype correlations given the isogenic background within this strain library. With this in mind, can you elaborate further on your choice of approach?

We acknowledge the reviewer’s point. A previous study has successfully identified mutations associated with intracellular persistence using libraries of transposon mutants (https://doi.org/10.1073/pnas.1520255113). However, the focus of our study was to identify naturally occurring and potentially clinically important adaptive mutations directly from a cohort of clinical isolates.

7. The manuscript describes using HeLa cells as an in vitro model system for professional phagocytic cells that are thought to be important for S. aureus clearance in vivo. The biological relevance of the HeLa cell system would be further supported by also studying some mutants in a system that employs professional phagocytes to study intracellular S. aureus persistence. For example, it would be of interest to know if the ausA mutants have a persistence phenotype in a macrophage or not.

The manuscript does not describe HeLa cells as an in vitro model for professional phagocytes. We agree with the reviewer regarding the verification of the phenotypic outputs identified using HeLa cells in professional phagocytes, which is the focus of future work.

8. The statistical 'gene burden GWAS' procedure only identified agr as having a significant P-value. In contrast, a manual 'convergence' analysis and hunt for homoplastic mutations for π uptake seemed to identify additional plausible candidates, but this manual procedure was not supported by statistical analysis. Can you explain why the convergence analysis would find loci that the 'gene burden GWAS' failed to identify? Can the mutations identified by the convergence analysis be further supported within a formal statistical framework?

We believe that our convergence analysis provides a complementary approach to allele-counting GWAS methods like that implemented in pyseer. The main advantage of convergence-based methods is the elimination of biases due to the underlying population structure. Pyseer uses a linear mixed model to correct for population structure, whereby the genetic relatedness of the isolates is expressed as kinship matrix and modelled as random effect. This method is widely used in human GWAS and might be sufficient for highly recombinogenic bacteria such as Streptococcus pneumoniae (see for example Lees et al., Nat Commun 2019). However, it is possible that this indirect correction for population structure is not complete in S. aureus, preventing the discovery of new genetic loci when performing a GWAS. By contrast, convergent evolution-based analysis eliminates biases due to the populations structure, with lower risk of identifying false positive. This is demonstrated by a large number of studies performed by both our group (Howden, PLoS Pathogen 2011; Pidot, Sci Transl Med 2018) and others (Das, PNAS 2016), who were able to identify biologically meaningful mutations from genetically close but-phenotypically divergent isolates. These studies were able to identify new phenotype-genotype associations without a statistical framework. However, we agree with the reviewer that providing statistical support would improve the confidence in the strength of the association and thus help prioritising genes for further investigation. In our response to Reviewer 1, we have now provided solid support for the contributions of 7 out of 20 genes (including agrA and ausA) in affecting intracellular cytotoxicity (see Figure 5 and Figure 5 Supplementary figure 1).

9. In the 'convergence' analysis, you only consider pairs with <200 SNP differences, which effectively increases the chance (i.e. statistical power) that any given mutation could be responsible for the observed π uptake effect size. What would happen if you re-ran the 'gene burden GWAS' using only a subset of the more genetically related strain pairs? Would that increase sensitivity enough to identify the homoplastic mutations you found in the 'convergence' analysis?

Whilst we agree with the importance of providing statistical support for our convergent evolution analysis, we feel that using an allele-counting approach like the one implemented in pyseer or other similar packages wouldn’t be appropriate for a small dataset of 56 paired closely related strains. Firstly, such approach would lead to drastic loss of power given the inclusion of 7-times less strains that in the original GWAS. Secondly the correction for population structure implemented in pyseer is not designed for pairs of closely related isolates. As mentioned above, we now provide a robust framework to support the associations between homoplasic mutations and cytotoxicity changes.

https://doi.org/10.7554/eLife.84778.sa2

Article and author information

Author details

  1. Abderrahman Hachani

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Conceptualization, Resources, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Stefano G Giulieri and Romain Guérillot
    For correspondence
    abderrahman.hachani@unimelb.edu.au
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8032-2154

  2. Stefano G Giulieri

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Conceptualization, Resources, Data curation, Formal analysis, Validation, Investigation, Methodology, Writing – original draft, Writing – review and editing
    Contributed equally with
    Abderrahman Hachani and Romain Guérillot
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5366-1943

  3. Romain Guérillot

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Conceptualization, Resources, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – review and editing, Supervision
    Contributed equally with
    Abderrahman Hachani and Stefano G Giulieri
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9915-1420

  4. Calum J Walsh

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  5. Marion Herisse

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Investigation, Methodology, Writing – original draft
    Competing interests
    No competing interests declared

  6. Ye Mon Soe

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4776-2523

  7. Sarah L Baines

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Resources, Data curation, Writing – original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0557-0518

  8. David R Thomas

    1. Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    2. Infection and Immunity Program, Department of Microbiology and Biomedicine Discovery Institute, Monash University, Clayton, Australia
    Contribution
    Formal analysis, Investigation
    Competing interests
    No competing interests declared

  9. Shane Doris Cheung

    Biological Optical Microscopy Platform, University of Melbourne, Melbourne, Australia
    Contribution
    Data curation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8733-4756

  10. Ashleigh S Hayes

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2926-7561

  11. Ellie Cho

    Biological Optical Microscopy Platform, University of Melbourne, Melbourne, Australia
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2932-1890

  12. Hayley J Newton

    1. Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    2. Infection and Immunity Program, Department of Microbiology and Biomedicine Discovery Institute, Monash University, Clayton, Australia
    Contribution
    Investigation, Writing – original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9240-2001

  13. Sacha Pidot

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Investigation, Methodology, Writing – original draft
    Competing interests
    No competing interests declared

  14. Ruth C Massey

    1. School of Microbiology, University College Cork, Cork, Ireland
    2. School of Medicine, University College Cork, Cork, Ireland
    3. APC Microbiome Ireland, University College Cork, Cork, Ireland
    4. School of Cellular and Molecular Medicine, University of Bristol, Bristol, United Kingdom
    Contribution
    Resources, Data curation
    Competing interests
    No competing interests declared

  15. Benjamin P Howden

    1. Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    2. Microbiological Diagnostic Unit Public Health Laboratory, Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Resources, Supervision, Funding acquisition, Writing – review and editing, Joint senior author
    Contributed equally with
    Timothy P Stinear
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0237-1473

  16. Timothy P Stinear

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    Contribution
    Resources, Supervision, Funding acquisition, Visualization, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Benjamin P Howden
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0150-123X

Funding

National Health and Medical Research Council (GNT2018880)

  • Abderrahman Hachani
  • Stefano G Giulieri
  • Romain Guérillot

National Health and Medical Research Council (GNT1105525)

  • Timothy P Stinear

National Health and Medical Research Council (GNT1145631)

  • Timothy P Stinear

National Health and Medical Research Council (GNT1196103)

  • Benjamin P Howden

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

Senior and Reviewing Editor

  1. Arturo Casadevall, Johns Hopkins Bloomberg School of Public Health, United States

Reviewer

  1. Matthew Culyba, University of Pittsburgh, United States

Version history

  1. Received: November 8, 2022
  2. Preprint posted: December 13, 2022 (view preprint)
  3. Accepted: May 23, 2023
  4. Version of Record published: June 8, 2023 (version 1)

Copyright

© 2023, Hachani, Giulieri, Guérillot et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 953
    Page views
  • 56
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Abderrahman Hachani
  2. Stefano G Giulieri
  3. Romain Guérillot
  4. Calum J Walsh
  5. Marion Herisse
  6. Ye Mon Soe
  7. Sarah L Baines
  8. David R Thomas
  9. Shane Doris Cheung
  10. Ashleigh S Hayes
  11. Ellie Cho
  12. Hayley J Newton
  13. Sacha Pidot
  14. Ruth C Massey
  15. Benjamin P Howden
  16. Timothy P Stinear
(2023)
A high-throughput cytotoxicity screening platform reveals agr-independent mutations in bacteraemia-associated Staphylococcus aureus that promote intracellular persistence
eLife 12:e84778.
https://doi.org/10.7554/eLife.84778

Categories and tags

Further reading

    1. Genetics and Genomics

    The DBL-1/TGF-β signaling pathway tailors behavioral and molecular host responses to a variety of bacteria in Caenorhabditis elegans

    Bhoomi Madhu, Mohammed Farhan Lakdawala, Tina L Gumienny
    Research Article

    Generating specific, robust protective responses to different bacteria is vital for animal survival. Here, we address the role of transforming growth factor β (TGF-β) member DBL-1 in regulating signature host defense responses in Caenorhabditis elegans to human opportunistic Gram-negative and Gram-positive pathogens. Canonical DBL-1 signaling is required to suppress avoidance behavior in response to Gram-negative, but not Gram-positive bacteria. We propose that in the absence of DBL-1, animals perceive some bacteria as more harmful. Animals activate DBL-1 pathway activity in response to Gram-negative bacteria and strongly repress it in response to select Gram-positive bacteria, demonstrating bacteria-responsive regulation of DBL-1 signaling. DBL-1 signaling differentially regulates expression of target innate immunity genes depending on the bacterial exposure. These findings highlight a central role for TGF-β in tailoring a suite of bacteria-specific host defenses.

    1. Epidemiology and Global Health
    2. Genetics and Genomics

    Within-host diversity improves phylogenetic and transmission reconstruction of SARS-CoV-2 outbreaks

    Arturo Torres Ortiz, Michelle Kendall ... Louis Grandjean
    Research Article

    Accurate inference of who infected whom in an infectious disease outbreak is critical for the delivery of effective infection prevention and control. The increased resolution of pathogen whole-genome sequencing has significantly improved our ability to infer transmission events. Despite this, transmission inference often remains limited by the lack of genomic variation between the source case and infected contacts. Although within-host genetic diversity is common among a wide variety of pathogens, conventional whole-genome sequencing phylogenetic approaches exclusively use consensus sequences, which consider only the most prevalent nucleotide at each position and therefore fail to capture low frequency variation within samples. We hypothesized that including within-sample variation in a phylogenetic model would help to identify who infected whom in instances in which this was previously impossible. Using whole-genome sequences from SARS-CoV-2 multi-institutional outbreaks as an example, we show how within-sample diversity is partially maintained among repeated serial samples from the same host, it can transmitted between those cases with known epidemiological links, and how this improves phylogenetic inference and our understanding of who infected whom. Our technique is applicable to other infectious diseases and has immediate clinical utility in infection prevention and control.

    1. Cancer Biology
    2. Genetics and Genomics

    Complementary CRISPR screen highlights the contrasting role of membrane-bound and soluble ICAM-1 in regulating antigen specific tumor cell killing by cytotoxic T cells

    Ann-Kathrin Herzfeldt, Marta Puig Gamez ... Lee Kim Swee
    Research Article

    Cytotoxic CD8+ T lymphocytes (CTLs) are key players of adaptive anti-tumor immunity based on their ability to specifically recognize and destroy tumor cells. Many cancer immunotherapies rely on unleashing CTL function. However, tumors can evade killing through strategies which are not yet fully elucidated. To provide deeper insight into tumor evasion mechanisms in an antigen-dependent manner, we established a human co-culture system composed of tumor and primary immune cells. Using this system, we systematically investigated intrinsic regulators of tumor resistance by conducting a complementary CRISPR screen approach. By harnessing CRISPR activation (CRISPRa) and CRISPR knockout (KO) technology in parallel, we investigated gene gain-of-function as well as loss-of-function across genes with annotated function in a colon carcinoma cell line. CRISPRa and CRISPR KO screens uncovered 187 and 704 hits respectively, with 60 gene hits overlapping between both. These data confirmed the role of interferon‑γ (IFN-γ), tumor necrosis factor α (TNF-α) and autophagy pathways and uncovered novel genes implicated in tumor resistance to killing. Notably, we discovered that ILKAP encoding the integrin-linked kinase-associated serine/threonine phosphatase 2C, a gene previously unknown to play a role in antigen specific CTL-mediated killing, mediate tumor resistance independently from regulating antigen presentation, IFN-γ or TNF-α responsiveness. Moreover, our work describes the contrasting role of soluble and membrane-bound ICAM-1 in regulating tumor cell killing. The deficiency of membrane-bound ICAM-1 (mICAM-1) or the overexpression of soluble ICAM-1 (sICAM-1) induced resistance to CTL killing, whereas PD-L1 overexpression had no impact. These results highlight the essential role of ICAM-1 at the immunological synapse between tumor and CTL and the antagonist function of sICAM-1.