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


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.



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.


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
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
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
  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.

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.

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
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
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, HFigure 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
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).


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).


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)
(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.


  1. Book
    1. Little TD
    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.
    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.

Article and author information

Author details

  1. Abderrahman Hachani

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    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
    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
    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
    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
    Competing interests
    No competing interests declared
  5. Marion Herisse

    Department of Microbiology and Immunology, Doherty Institute, University of Melbourne, Melbourne, Australia
    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
    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
    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
    Formal analysis, Investigation
    Competing interests
    No competing interests declared
  9. Shane Doris Cheung

    Biological Optical Microscopy Platform, University of Melbourne, Melbourne, Australia
    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
    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
    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
    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
    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
    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
    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
    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


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.

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)


© 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.


  • 1,529
  • 109
  • 0

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

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
A high-throughput cytotoxicity screening platform reveals agr-independent mutations in bacteraemia-associated Staphylococcus aureus that promote intracellular persistence
eLife 12:e84778.

Share this article


Further reading

    1. Computational and Systems Biology
    2. Genetics and Genomics
    Ardalan Naseri, Degui Zhi, Shaojie Zhang
    Research Article Updated

    Runs-of-homozygosity (ROH) segments, contiguous homozygous regions in a genome were traditionally linked to families and inbred populations. However, a growing literature suggests that ROHs are ubiquitous in outbred populations. Still, most existing genetic studies of ROH in populations are limited to aggregated ROH content across the genome, which does not offer the resolution for mapping causal loci. This limitation is mainly due to a lack of methods for the efficient identification of shared ROH diplotypes. Here, we present a new method, ROH-DICE (runs-of-homozygous diplotype cluster enumerator), to find large ROH diplotype clusters, sufficiently long ROHs shared by a sufficient number of individuals, in large cohorts. ROH-DICE identified over 1 million ROH diplotypes that span over 100 single nucleotide polymorphisms (SNPs) and are shared by more than 100 UK Biobank participants. Moreover, we found significant associations of clustered ROH diplotypes across the genome with various self-reported diseases, with the strongest associations found between the extended human leukocyte antigen (HLA) region and autoimmune disorders. We found an association between a diplotype covering the homeostatic iron regulator (HFE) gene and hemochromatosis, even though the well-known causal SNP was not directly genotyped or imputed. Using a genome-wide scan, we identified a putative association between carriers of an ROH diplotype in chromosome 4 and an increase in mortality among COVID-19 patients (p-value = 1.82 × 10−11). In summary, our ROH-DICE method, by calling out large ROH diplotypes in a large outbred population, enables further population genetics into the demographic history of large populations. More importantly, our method enables a new genome-wide mapping approach for finding disease-causing loci with multi-marker recessive effects at a population scale.

    1. Chromosomes and Gene Expression
    2. Genetics and Genomics
    Lisa Baumgartner, Jonathan J Ipsaro ... Julius Brennecke
    Research Advance

    Members of the diverse heterochromatin protein 1 (HP1) family play crucial roles in heterochromatin formation and maintenance. Despite the similar affinities of their chromodomains for di- and tri-methylated histone H3 lysine 9 (H3K9me2/3), different HP1 proteins exhibit distinct chromatin-binding patterns, likely due to interactions with various specificity factors. Previously, we showed that the chromatin-binding pattern of the HP1 protein Rhino, a crucial factor of the Drosophila PIWI-interacting RNA (piRNA) pathway, is largely defined by a DNA sequence-specific C2H2 zinc finger protein named Kipferl (Baumgartner et al., 2022). Here, we elucidate the molecular basis of the interaction between Rhino and its guidance factor Kipferl. Through phylogenetic analyses, structure prediction, and in vivo genetics, we identify a single amino acid change within Rhino’s chromodomain, G31D, that does not affect H3K9me2/3 binding but disrupts the interaction between Rhino and Kipferl. Flies carrying the rhinoG31D mutation phenocopy kipferl mutant flies, with Rhino redistributing from piRNA clusters to satellite repeats, causing pronounced changes in the ovarian piRNA profile of rhinoG31D flies. Thus, Rhino’s chromodomain functions as a dual-specificity module, facilitating interactions with both a histone mark and a DNA-binding protein.