In order to study the common and varying mechanisms of biofilm formation in MRSA clones, three MRSA strains representing three pandemic MRSA clones were used: (i) UAS391 [USA300, PVL+/ACME+, abscess/wound isolate, 2591 coding sequences (CDS); mobile content: 5.12%] belonging to ST8, harboring SCCmec type IV, previously identified as a prolific biofilm former; (ii) HEMRSA-15 (EMRSA-15; unknown source, 2562 CDS; mobile content: 7.75%) belonging to ST22, harboring SCCmec type IV; and (iii) ST239-16 (Hungarian/Brazilian ST239, respiratory tract isolate, 2785 CDS, mobile content: 12.18%) belonging to ST239, harboring SCCmec type III. These strains have been shown previously by us to have strikingly different biofilm phenotypes in terms of biomass or cell-to-matrix ratios (Vanhommerig et al., 2014).

First, we studied growth rates of each strain to understand if these differences might be due to differences in replication rates. Next, biofilm formation was assessed in a static biofilm model with a high-nutrient environment, followed by quantification after 24 h. Subsequently, these strains were allowed to form static biofilms for 24, 48, and 72 h with planktonic cultures generated in parallel. Samples for total RNA extraction were obtained at the given time points for both phenotypes and used for transcriptomic analyses. Genes found to be differentially expressed (at least 2-fold) between (i) the same strain’s biofilm and planktonic phenotypes at the same time point, or (ii) the same strain’s consecutive biofilm time points were studied further using metabolic pathway analysis. In parallel, in silico co-expression networks were constructed to observe inter- and intrastrain patterns.

Last, since it has been previously shown that presence of SCCmec influences MRSA biofilm formation (Pozzi et al., 2012), we investigated this further in the pandemic USA300 clone. Also, the arginine catabolic mobile element (ACME), which is unique to USA300, is shown to aid in skin colonization. Both SCCmec and ACME were excised to study their role in biofilm formation by USA300. Subsequently, the effect of single gene disruptions in ACME on biofilm formation was also investigated using Bursa aurealis transposon mutants.

All strains were grown either on 5% horse blood Columbia agar plates (Becton, Dickinson & Company, New York, NY, United States) or in brain heart infusion (BHI; Becton, Dickinson & Company, New York, NY, United States) broth supplemented with 0.1% D(+)-glucose monohydrate (Merck Millipore, Boston, MA, United States) as previously described (De Backer et al., 2018). Bursa aurealis transposon mutants in S. aureus USA300-JE2 were provided by the Network on Antimicrobial Resistance in Staphylococcus Aureus (NARSA) for distribution by BEI Resources (Nebraska Transposon Mutant Library (NTML),¹ Supplementary Table 1; Fey et al., 2013).

Growth rates of USA391, HEMRSA-15, and ST239-16 cultures were determined as previously described (De Backer et al., 2018). At each time point, average optical densities were calculated from quadruplicates. Additionally, biofilms of all three strains, as well as the NARSA transposon mutants, were formed on pre-sterilized 96-well flat-bottom polystyrene microtiter plates (CELLSTAR^®, Greiner Bio-One, Kremsmünster, Austria) in triplicate to measure biofilm biomass using the crystal violet assay as described previously (De Backer et al., 2018).

All three strains, UAS391, HEMRSA-15, and ST239-16, stored at −80°C underwent two subsequent subcultures on blood agar (see section “Bacterial growth conditions” for media). A single colony was utilized to start an overnight broth and agar subculture from the same plate. The overnight broth culture was utilized to start and collect, in the exponential growth phase, 0.05 McF of each strain that was inoculated in a 136-mm polystyrene Petri dish (Greiner Bio-One, Austria), and grown under static conditions for 24, 48, or 72 h before rinsing (PBS, Thermo Fisher Scientific Inc., Boston, MA, United States) and scraping off adherent bacteria. Simultaneously, one colony from the blood agar subcultures was re-suspended in medium and incubated for 24, 48, or 72 h with shaking (250 rpm, 37°C) for each strain. The growth medium was refreshed every 24 h by spinning down the culture, removing the supernatant, and adding the fresh medium. After incubation, bacteria were collected by centrifugation at 4,600 × g for 20 min. Subsequently, approximately 0.200 g of both biofilm and planktonic bacterial suspensions were placed in RNAprotect™ Bacteria Reagent (Qiagen, Hilden, Germany) and incubated for 5 min at room temperature to stabilize the RNA. Bacteria were then harvested by centrifugation at 5,000 × g for 10 min at 4°C. Then two biological repeats were analyzed for each strain, phenotype, and time point.

Total RNA was purified from harvested bacteria using the Masterpure™ Complete DNA and RNA Purification kit (Epicentre^®, Biotechnologies, Madison, WI, United States) according to the manufacturer’s protocol with modifications. Briefly, cells were mechanically disrupted using Lysing Matrix B 0.1 mm silica sphere bead beating (FastPrep^® 24 instrument, MP Biomedicals Inc., Irvine, CA, United States) before extraction. Each RNA sample was suspended in 30 μL of RNA storage solution, and the quality was determined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Germany). Ribosomal RNA (rRNA) was depleted using a Ribo depletion kit (Thermo Fisher Scientific Inc., Netherlands), followed by preparation of a stranded TruSeq library and 1 × 75 bp sequencing using the NextSeq 500 system, v2 (Illumina Inc., San Diego, CA, United States).

Raw reads were trimmed by Trimmomatic (v.0.39) (Bolger et al., 2014), and FastQC (v.0.11.8)² was used for quality assessment. All processed data were normalized to the same sequencing depth (4,464,784 reads per sample) and then used for alignment by BBMap (v.38.47) against their reference genomes with default parameters to obtain transcripts per million (TPM) values (Bushnell, 2014). For UAS391 and HEMRSA-15, reference genomes (GenBank: CP007690.1 and CP007659.1, respectively) were obtained from the National Center for Biotechnology Information (NCBI). ST239-16 was whole-genome-sequenced in-house, trimmed as described previously, assembled by SPAdes (Bankevich et al., 2012), and annotated using Prokka (Seemann, 2014). DEGs (up- or downregulated) identified between (i) the biofilm and planktonic phenotype for each clone and (ii) in early (24 h), intermediate (48 h), and late (72 h) stages during biofilm formation were mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database of S. aureus. Expression levels of housekeeping and three randomly selected genes were validated using RT-PCR. Selected genes and used primers are listed in Supplementary Table 2.

For each strain under biofilm or planktonic phenotypes, six samples representing the three time points (24, 48, and 72 h) were used to study temporal changes in gene (co-)expression by generating a network based on Spearman correlations with p_adj < 0.01 and |r| > 0.6. In total, six networks were generated for both the biofilm and planktonic phenotypes of the three strains. Networks were visualized by Gephi (v0.9.2) (Bastian et al., 2009), and their topological properties were calculated using the igraph package in R (v3.6) (Csárdi and Nepusz, 2006; R Core Team, 2020). In total, 1,000 Erdös–Réyni random networks, with equal numbers of nodes and edges as the corresponding observed network, were created using the igraph R package (Csárdi and Nepusz, 2006).

In UAS391, both ACME and SCCmec or single elements were cured by electroporation (20°C, 100 Ω, 25 μF, and 2.3 kV) of pSR2, a thermosensitive plasmid containing the ccrAB2 gene complex and a tetracycline resistance marker, isolated from Escherichia coli SF8300 (ZymoPURE plasmid miniprep kit, Zymo Research, Irvine, CA, United States), using the method described by Katayama et al. (2000) and Diep et al. (2008). Here, the pSR2 plasmid was first transferred into the restriction-deficient intermediate S. aureus cloning host RN4220 to adapt the plasmid DNA to S. aureus modifications. RN4220 and UAS391 transformants were grown for three consecutive days at 30°C in BHI broth supplemented with tetracycline. Single colonies were picked and grown for 24 h in the BHI broth at 42°C to promote loss of pSR2. The deletion of ACME, SCCmec, or ACME and SCCmec was confirmed by PCR using X1-X5, X2-X3, and X1-X3 PCR primers (Diep et al., 2008). Additionally, transformants were screened for the presence of arc and opp3 gene clusters by PCR-based assays, using the primer pairs AIPS.27 and AIPS.28 for arcA, and AIPS.45 and AIPS.46 for opp3AB (Diep et al., 2008). SCCmec typing was based on ccr recombinase and mec PCRs (Kondo et al., 2007).

For RNASeq data, DESeq2 (v.1.14.1) was used to determine significantly differently expressed genes (DEGs) in each pairwise comparison based on TPM values of every gene (Love et al., 2014). Specifically, genes with p < 0.05 and fold change (FC) ≥ 2 were extracted, and their p-values were adjusted by using the false discovery rate (FDR) approach (p_adj). Finally, genes with p_adj < 0.05 and FC ≥ 2 in each pairwise comparison were considered as DEGs. One-way ANOVA was conducted in R (v3.6) to evaluate the differences in biomass and growth rates at the three time points between the three tested strains. Transcriptomic profiles across both phenotypes and time points were compared using the permutational multivariate analysis of variance (PERMANOVA) using the Adonis function with 999 permutations in the vegan package³ in R.

In the static biofilm assay, HEMRSA-15 (OD492 = 0.856) and ST239-16 (OD492 = 1.366) showed higher biomass after 24 h of growth (+16% and +83%, respectively; p ≤ 0.03) than UAS391 (OD492 = 0.746, Supplementary Figure 1A). However, in the growth assay (OD620 measured after 24 h of growth at 37°C), the high biomass-forming HEMRSA-15 (0.149 min^–1, doubling time: 4.65 min) and ST239-16 (0.160 min^–1, doubling time: 4.31 min) showed a lower growth rate than UAS391 (0.164 min^–1, doubling time: 4.21 min; p ≤ 0.046, Supplementary Figure 1B).

Principal component analysis in each of the three strains revealed substantial differences in transcriptional profiles between biofilm and planktonic phenotypes, as well as across the three time points (Figure 1A and Supplementary Tables 3, 4). After 24 h of growth, there were 120, 238, and 43 differentially expressed genes (DEGs) (p_adj ≤ 0.05 and FC ≥ 2) identified between the biofilm and planktonic phenotypes for UAS391, HEMRSA-15, and ST239-16, respectively (Supplementary Figure 2 and Supplementary Tables 5–7). This corresponded to 4.63, 9.29, and 1.54% of all CDSs, respectively. Among these genes, many hypothetical and phage-related genes were observed. As a control, the expression of five housekeeping genes (gyrB, groE, glpF, gmk, and yqiL) did not show significant differences between the biofilm and planktonic phenotypes in any of the three strains (p_adj > 0.05, Supplementary Table 8). Additionally, expression levels of randomly selected DEGs in the biofilm phenotype were confirmed by RT-PCR (Supplementary Figure 3).

All identified DEGs between the biofilm and planktonic phenotypes after 24 h of growth were mapped against the KEGG database with 56 (46.67%), 107 (44.96%), and 13 (30.23%) DEGs assigned KO terms for UAS391, HEMRSA-15, and ST239-16, respectively (Figure 1B). When comparing the gene expression of the top 10 mapped up and downregulated genes in biofilm with the planktonic phenotype in UAS391, HEMRSA-15, and ST239-16 at 24, 48, and 72 h, common genes were observed between UAS391 and HEMRSA-15 after 24 h of growth (Table 1). At 48 h, common genes were observed between UAS391 and ST239-16 but showing opposite expression patterns, while at 72 h, only one common gene between UAS391 and HEMRSA-15 was identified (Table 1). Genes involved in the KEGG pathways “Staphylococcus aureus infection,” “Transporters,” “Transcription,” and “Peptidases” were upregulated in all three biofilm phenotypes compared to their corresponding planktonic phenotypes. However, except these four KEGG pathways, transcriptomic signatures differed across the three strains. The up- and downregulated genes in UAS391 were involved in 17 and 10 KEGG pathways, respectively, in HEMRSA-15 in 34 and 24 KEGG pathways, respectively, and in ST239-16 in 10 and 1 KEGG pathways, respectively (Figure 1). UAS391 and HEMRSA-15 shared eight and five KEGG pathways with up- and downregulated genes, respectively, while HEMRSA-15 also showed 19 unique KEGG pathways, each with up- and downregulated genes. The total number of affected pathways in ST239-16 (n = 11) was remarkably lower than that in UAS391 (n = 27) and HEMRSA-15 (n = 58). Interestingly, the majority of the affected KEGG pathways were not linked to the accessory genome, indicating that expression of core genome pathways differed substantially in 24-h biofilms between the studied strains (Figure 1B).

Of note, in the “Staphylococcus aureus infection” pathway, the only DEG that was commonly upregulated in all three strains in 24-h biofilms was clfA (2.20-, UAS391; 2.07-, HEMRSA-15; and 2.66-fold, ST239-16; Table 2). The clfA gene encodes clumping factor A, an important virulence factor also involved in initial cell adherence. Other virulence genes such as those encoding the infection-related factors CHIPS (chemotaxis-inhibiting protein of S. aureus, chs), staphylokinase (sak), and a complement inhibitor SCIN family protein (scn) were downregulated in 24-h biofilms of UAS391 (−2.36-, −1.66-, and −2.21-fold, respectively) and HEMRSA-15 (−3.33-, −2.92-, and −3.20-fold, respectively) compared to planktonic bacteria (Table 2). Additionally, staphylococcal protein A (lysM) and gamma hemolysin subunit B (hlgB) showed downregulation in HEMRSA-15 (−2.11- and −2.33-fold, respectively, Table 2). In ST239-16 biofilms, scn and lysM were not differentially expressed compared to the planktonic phenotype (p_adj > 0.05), but virulence factors chs, sak, and hlgB were absent in the ST239-16 genome.

Since we previously showed the influence of the arginine-synthesizing, urea, and tricarboxylic acid (TCA) cycles on USA300 biofilm formation (De Backer et al., 2018), arginine-related pathways were further investigated here. Under anoxic conditions, found in biofilms, arginine can provide energy in the form of ATP through catabolism with the formation of ammonia through the arginine deiminase (ADI) pathway, or be utilized as a substrate in the urea cycle with additional ammonia generation (Cunin et al., 1986). In 24-h biofilms of UAS391 and HEMRSA-15, the ADI pathway (arcA, arcB, arcC2_1, arcD, and arcR) was upregulated at least 2- and 3-fold, respectively, compared to the planktonic phenotype (Table 2). Among these, arcA and arcC2_1 were most highly upregulated (2.98- and 3.95-fold, respectively) in UAS391 and HEMRSA-15 24-h biofilms. In line with these results, transposon inactivation of arcA in USA300-JE2 proved to have a significant impact on produced biofilm, compared to the wild type (OD492 = 0.372, 57%; p ≤ 0.001, Supplementary Figure 4). Moreover, arcC_2 inactivation decreased the biofilm production to an equal level, as observed in the arcA:Tn mutant (OD492 = 0.373, 57%; p ≤ 0.001, Supplementary Figure 4). Contrasting the UAS391 and HEMRSA-15 biofilms, in ST239-16 biofilm, the upregulation of the ADI pathway was only observed after 48 h of biofilm growth with all elements, except arcB, showing a minimum of 2-fold upregulation compared to planktonic cells (Table 2).

Arginine can be synthesized in the arginine biosynthesis pathway (argBCJ) coupled with the urea cycle (argFGH, rocF) (Cunin et al., 1986). Expression levels were stable in 24-h biofilm phenotypes of all three strains with no DEGs detected in both pathways when comparing to respective planktonic phenotypes (Table 2). On the other hand, in 24-h biofilms, differential expression was observed for genes related to urease (ureABCDEFG), which catabolizes urea formed by arginase (rocF) in the urea cycle. In UAS391 24-h biofilms, all urease-related genes showed FCs ranging from 1.99 to 2.37 (Table 2). Another arginine-synthesizing pathway, the pyrimidine nucleotide biosynthesis pathway (pyrRPBC, carAB, and pyrFE), did not show significant differential expression in 24-h biofilms of the three strains (p > 0.05).

In addition to arginine catabolism-mediated ammonia production, the high-affinity K⁺-specific transport system (kdpABC and kdpDE) also contributes to pH homeostasis through the exchange of K⁺ for H⁺ (Cotter and Hill, 2003). Also here, 24-h biofilms showed differential expression when compared to its planktonic counterpart: in UAS391, kdpABC was upregulated (2.25-, 2.23-, and 2.29-fold, respectively), and in HEMRSA-15, both kdpABC and kdpDE were upregulated (2.61-, 3.17-, 3.01-, 1.08-, and 1.94-fold, respectively). By contrast, in the ST239-16 24 h biofilm, this was not observed (Table 2).

As with 24-h biofilms, gene expression patterns in 48- and 72-h biofilms showed substantial differences from their planktonic counterparts, as shown by PCA (Figure 1A). Of all identified DEGs, KEGG mapping could be performed for 595 (65.7%), 70 (25.4%), and 250 (49.3%) of the identified DEGs for the 48-h time point for UAS391, HEMRSA-15, and ST239-16, respectively (Figure 2). After 72 h of growth, the number of mapped DEGs were limited to 72 (82.8%), 2 (20%), and 11 (40.7%) for the three strains, respectively (Supplementary Figure 5). As with the 24-h time point, transcriptomic signatures differed remarkably between strains; especially, UAS391 showed a significant amount of affected KEGG pathways, which were unaffected in the two other strains (Figure 2). In contrast to the 24-h time point, more affected pathways were linked to the accessory genome. These pathways were also affected in two or all three strains, indicating this change in gene expression was common across all three strains. Also noticeable were a number of core genome-linked KEGG pathways, which were downregulated in UAS391 but upregulated in ST239-16, again indicating differences in gene expression patterns between the strains (Figure 2).

Comparing temporal biofilm time points within strains, differential expression analysis showed high numbers of DEGs: 48-h biofilm to 24-h biofilm [653 (UAS391), 563 (HEMRSA-15), and 422 (ST239-16)] and 72-h biofilm to 48-h biofilm [3 (UAS391), 198 (HEMRSA-15), and 406 (ST239-16), Supplementary Figure 6]. Housekeeping genes did not show significant differences between time points, and RT-PCR confirmed expression levels of random DEGs (Supplementary Table 9 and Supplementary Figure 3). Furthermore, null model testing showed that the global gene expression patterns were non-stochastic and therefore could not be attributed to random fluctuations at the single-gene level (Supplementary Table 10).

However, as with expression patterns in 24- and 48-h biofilms compared to the planktonic phenotype, ST239-16 showed striking differences during biofilm maturation compared to the other strains (Supplementary Figure 6). In the 24- to 48-h biofilm timeframe, expression levels of several metabolism pathways showing downregulation in UAS391 and HEMRSA-15, such as glycolysis/gluconeogenesis and oxidative phosphorylation, showed upregulation in ST239-16 (Supplementary Figure 6). Also, genes encoding ribosomal proteins and the biogenesis of ribosomes showed this trend with 65 and 56 downregulated genes between 24- and 48-h biofilms for UAS391 and HEMRSA-15, respectively, but 41 upregulated genes in ST239-16 (Supplementary Figure 6).

The ADI pathway, showing significant upregulation in 24-h biofilms of UAS391, ST239-16, and HEMRSA-15 compared to planktonic cells, did not show differential expression in biofilms of 48 h and 72 h for both strains (Table 2). However, when comparing 48-h UAS391 biofilm to 48-h planktonic cells, all components were downregulated in the biofilm phenotype with FC ranging from −1.29 to −2.8-fold (Table 2). In HEMRSA-15, this pattern was not observed. Similar results were observed for urease-related genes, upregulated in 24-h UAS391 biofilm, with most components downregulated comparing 48-h UAS391 biofilm to 24-h biofilm (FC: −1.49 to −2.08-fold) and 48-h UAS391 planktonic (FC: −1.2 to −1.94-fold, Table 3), suggesting a distinct role in early UAS391 biofilm formation.

Interestingly, the pyrimidine nucleotide biosynthesis pathway, unaffected in 24-h biofilms, showed significant downregulation of all components (except pyrP) in 48-h UAS391 biofilms, compared to the 24-h time point, with FC ranging from −1.71 to −3.61 (Table 3). By contrast, HEMRSA-15 and ST239-16 48-h biofilms only showed differential expression of two components: pyrFE (−1.07- and −1.01-fold, respectively) and pyrF and pyrR (1.08- and 1.38-fold, respectively, Table 3). On the other hand, the high-affinity K⁺-specific transport system, upregulated in UAS391 and HEMRSA-15 24-h biofilms compared to planktonic cells, showed subsequent downregulation in 48-h biofilms (Table 2). KpdABC and kpdE showed negative FCs when comparing 48-h with 24-h biofilms for both strains. In the 72-h biofilms, none of the genes was identified as DEGs (Table 3).

Topological properties in the observed networks significantly differed from those in the corresponding random networks (p < 0.05, one-sample t-test, Supplementary Table 10), indicating that the observed patterns of gene co-expression along the time points (24, 48, and 72 h) in each strain were not stochastic, in line with the null model results that the observed transcriptional profiles significantly differed (p < 0.01) from the expected profiles (Supplementary Table 11). Positive connections dominated in all networks, indicating genes showed similar expression patterns along the time points. Gene co-expression patterns were more similar between UAS391 and HEMRSA-15, with higher network connectivity (indicated by average degree, and numbers of nodes and edges) than ST239-16 networks. In UAS391 and HEMRSA-15, there was lower network connectivity in biofilms than in planktonic networks, while an opposite pattern was observed in ST239-16 networks. Inside the modules of the UAS391 and HEMRSA-15 networks, some specific pathways such as ribosomal protein synthesis/activity, sugar metabolisms, and/or even oxidative phosphorylation dominated. However, no domination of specific pathways in the ST239-16 module was observed (Supplementary Figure 7).

To observe the impact of the ACME, which is uniquely carried by USA300, on biofilm formation, we excised the element, which resulted in a small deficit in biofilm formation capacity (OD492 = 0.619, 83%, p = 0.13), as compared to the wild type (Supplementary Figure 8A). Knocking out each gene on the ACME revealed a large biofilm decrease when opp3D was disrupted in USA300-JE2 (OD492 = 0.372; 57%; p ≤ 0.05, Supplementary Figure 8B), while other transposon disruptions only resulted in small or no biofilm defects. Since our results showed a significant upregulation of the ADI pathway in 24-h biofilms, it was also of interest that ACME harbors an additional ADI gene cluster. We observed upregulation of arcA, arcB, arcC, and arcD in UAS391 biofilms. Here, arcABCD was on average 2-fold upregulated in 48-h biofilms, compared to the 24-h time point (Table 3).

Arginine catabolic mobile element is closely associated with methicillin resistance encoding SCCmec, where transcriptomic data showed no differential expression of SCCmec genes between biofilms and planktonic phenotypes in all three strains (p > 0.05). However, lower expression of mecA and mecR1 components was observed in later biofilm time points (48 and 72 h) compared to 24 h in UAS391 and HEMRSA-15. The deletion of SCCmec, as with ACME, resulted in a small biomass decrease of 18% (OD492 = 0.429, 82%; p = 0.22, Supplementary Figure 8A). On the other hand, combined deletion of both elements resulted in a much larger biofilm defect (OD492 = 0.614, 58%; p < 0.001).