Balancing the Virulence and Antimicrobial Resistance in VISA DAP-R CA-MRSA Superbug

Background: Methicillin-resistant Staphylococcus aureus (MRSA) with intermediate resistance to Vancomycin (VISA) is reported worldwide. These strains frequently emerge among hospital-associated (HA)-MRSA and rarely within community-acquired (CA)-MRSA. Here, the genomic and transcriptomic adaptations distinguishing VISA daptomycin resistant (DAP-R) CA-MRSA, which emerged in a hospitalized patient under glycopeptide treatment, were explored. Methods: Whole-genome sequencing, RNA-Seq and bioinformatics were carried out. Results: Our CA-MRSA clustered in the USA400 lineage showing additional antimicrobial resistance (AMR) versus DAP and glycopeptides. Resistomics revealed adaptations related to glycopeptide, daptomycin and rifampin resistance (mprF nsSNPS and overexpression of glycopeptide and daptomycin-resistance related genes). Similar changes were detected in virulence traits (agrA HI-nsSNPs and toxin gene underexpression), in which a decrease was observed despite the abundance of virulence-related genes. Our results predicted a balance in adaptations, decreasing the virulence and biological costs to support the co-occurrence of extensive AMR in a hypervirulent genomic background. Conclusion: Our data show that VISA DAP-R CA-MRSA shifts the potential hypervirulent behavior of CA-MRSA towards the acquisition and maintenance of extensive AMR, by a decrease in virulence and biological costs mediated by a “compensatory modulatory mutation” silencing the Agr quorum-sensing cascade.


Introduction
Hospital-and community-acquired methicillin-resistant Staphylococcus aureus (MRSA) are considered high-priority micro-organisms because the severity of their infections is often difficult to treat [1]. Exploiting a plethora of complex infection-related mechanisms of infection, S. aureus causes a wide range of serious infections, i.e., abscesses of various organs, bacteremia, infective endocarditis (IE), pneumonia, as well as osteoarticular, skin and soft tissue, pleuropulmonary and device-related infections. Complicated therapeutical management is often necessary to treat MRSA infections, even with the emergence of new Gram-positive drugs [1]. From an evolutionary viewpoint, S. aureus constantly has reprogrammed its antimicrobial resistance (AMR) mechanisms in response to increasing antimicrobial pressure and virulence mechanisms to adapt its host defense response. Antimicrobial and immune system activity are key factors in conditioning bacterial survival [2][3][4][5][6][7] and in the selection of new lineages or pathogen variants. Both mechanisms are crucial for survival under adverse conditions such as host immune system defenses, antimicrobial treatment and the need to survive in new challenging niches.
Reflecting these two evolutionary and selective drives, MRSA selected three different variants related to the modality of transmission: (i) healthcare-associated methicillin-resistant S. aureus (HA-MRSA), (ii) community-acquired S. aureus (CA-MRSA) and (iii) livestockassociated S. aureus (LA-MRSA), each of which is strongly adapted to different environments.
Among these, HA-MRSA infections are mainly reported following contact with healthcare settings. Community-associated MRSA (CA-MRSA) infections have been reported since the 1990s in individuals who had no prior hospitalization [8], whilst LA-MRSA are associated with livestock [9]. CA-MRSA is particularly virulent compared to HA-MRSA in that, on the contrary, it is highly resistant to antimicrobial treatment [10,11]. This aspect reflects different genomics of HA-MRSA versus CA-MRSA. HA-MRSA is characterized by a huge pool of antimicrobial resistance (AMR) genes and lower virulence-related genes. CA-MRSA has a lower content of AMR genes (in general restricted to β-lactams and macrolides) and a higher pool of virulence-related genes, including the Panton-Valentine Leukocidin (PVL) toxin gene. This is a β-pore-forming cytotoxin creating pores in the membranes of infected cells that can cause necrotic lesions in the skin, in the mucosa and can also determine necrotic hemorrhagic pneumonia [10,11].
Virulence and resistance factors share similar mechanisms of dissemination and coselection between species or genera. These mechanisms are mediated by horizontal gene transfer (HGT) of mobile genetic elements (MGEs), compensatory or adaptive mutations [12,13] and intrinsic factors related to the nature of micro-organisms such as biofilm-producing or intracellular bacteria [2,3,6,7], involvement of the same cellular structures such as efflux pumps [4], porins [5], cell wall alterations [6] and two-component regulatory systems [7].
MRSA antimicrobial resistance towards first-choice and last resort antimicrobials such as glycopeptides and daptomycin is one of the great challenges for the current management of MRSA infections. These strains emerge typically and frequently among HA-MRSA because of the selective pressure arising during severe antimicrobial treatment; on the other hand, they are very rare among CA-MRSA [14,15].
In our experimental plan, global integrated genomic and transcriptomic profiling was applied to better understand the adaptations of a new VISA DAP-R CA-MRSA superbug that attracted our attention for its biological characteristic of extensive antimicrobial resistance, including glycopeptide and daptomycin resistance.

Results
A comparative genomic and transcriptomic approach was exploited to uncover the gained adaptations by VISA DAP-R CA-MRSA, acquired under the selective pressure of teicoplanin treatment, versus its VSSA DAP-S CA-MRSA isogenic parents.
First, the genomic lineage and background characterizing VISA DAP-R CA-MRSA, as well as the genomic markers associated with the take-over of extensive antimicrobial resistance (VISA and DAP-R), were determined. In-depth genomic characterization (phylogeny, genomic epidemiology, molecular typing, whole-genome non-synonymous single nucleotide polymorphisms) was conducted by whole-genome sequencing (WGS) and bioinformatics, focusing mainly on the traits impacting AMR and Virulence.
Similarly, a pool of the statistically significant differentially expressed genes (DEGs), their expression trend (over-or under-expression) and the main enriched dysregulated KEGG pathways were determined by RNA-seq and bioinformatics.
Finally, the phenotypical biofilm production ability (together with already known delta-hemolysis production ability) and growth-curve kinetics were evaluated to assess the production of prototype (adhesion and toxigenicity, respectively) virulence factors and the biological AMR and virulence fitness costs.

Biological Fitness Costs of AMR and Virulence
AMR and virulence biological fitness costs were evaluated in terms of variation in generation time and lag growth-phase duration in 1-R vs 1-S, as previously published [16]. A longer generation time was found in 1-R (60 ) than in 1-S (30 ). In addition, growth kinetic experiments showed that the lag growth-phase of VISA DAP-R CA-MRSA 1-R was greater (2.00 h) than VSSA DAP-S CA-MRSA 1-S (1.0 h), as well as more similar to MW2 CA-MRSA (2.30 h) and S. aureus ATCC 29213 (2.25 h) (Table 1).

Biofilm Production
Biofilm production assays showed that 1-S and 1-R MRSA strains were not biofilm producers ( Table 2).

Phylogeny
CSI phylogeny clearly showed the clusterization of 1-S and 1-R MRSA strains in the USA-400 CA-MRSA phylogenetic lineage. The phylogenetic tree showed a close relationship between the 1-S/R MRSA strain pair and the ST-1 MRSA cluster. This highlighted the strong phylogenetic and genomic relationship with the MW2 CA-MRSA, notoriously a USA400, agr-III, CA-MRSA prototype (Figure 1).
Virulome analysis found the major staphylococcal virulence factors including 22 adherence related genes, 22 immune evasion related genes, 8 secretion system coding genes, 14 exoenzyme encoding genes and 26 toxin coding genes, including lukS/F encoding the PVL in both strains (Table 3).

Discussion
HA-MRSA and CA-MRSA emerged as different variants of clear adaptive events that occurred within the same species. The genomic composition of the two MRSA lineages clearly indicated that evolution led to highly balanced genomes with HA-MRSA strongly adapted to survive in extreme antimicrobial selective pressure whilst CA-MRSA overcame the host-response [18][19][20][21].
In this scenario, our data describe, for the first time, a new CA-MRSA superbug that emerged as a result of a convergent cross-selection mechanism that makes the CA-MRSA superbug able to decrease its hypervirulence, supporting and maintaining extensive antimicrobial resistance.
These new genomically related USA-400, ST-1, spatype-t127, agr-III, SCCmecIVa, PVL-positive, glycopeptide and daptomycin-resistant CA-MRSA were characterized by the co-occurrence of an extensively antimicrobial resistance profile and decreased virulence despite their great virulence potential.
New insights arose on the layout of genomic and transcriptomic adaptations occurring in this exceptional circumstance in which MRSA harbors simultaneous adaptations related to and impacting antimicrobial resistance and virulence.
Two convergent key motifs emerged as key points of its success: (i) the selection of a compensatory "Silencing regulatory" mutation in the agr-locus, strategic quorum sensing regulator pathway, involved in virulence and antibiotic resistance; (ii) dysregulation in essential primary metabolic pathways associated with a prevalent underexpression trend in DEGs. These co-events determine and allow a balancing of the antimicrobial resistance and virulence adaptations that lead to the setup of increased antimicrobial resistance by decreasing virulence and restoring intrinsic biological costs.

Discussion
HA-MRSA and CA-MRSA emerged as different variants of clear adaptive events that occurred within the same species. The genomic composition of the two MRSA lineages clearly indicated that evolution led to highly balanced genomes with HA-MRSA strongly adapted to survive in extreme antimicrobial selective pressure whilst CA-MRSA overcame the host-response [18][19][20][21].
In this scenario, our data describe, for the first time, a new CA-MRSA superbug that emerged as a result of a convergent cross-selection mechanism that makes the CA-MRSA superbug able to decrease its hypervirulence, supporting and maintaining extensive antimicrobial resistance.
These new genomically related USA-400, ST-1, spatype-t127, agr-III, SCCmecIVa, PVLpositive, glycopeptide and daptomycin-resistant CA-MRSA were characterized by the co-occurrence of an extensively antimicrobial resistance profile and decreased virulence despite their great virulence potential.
New insights arose on the layout of genomic and transcriptomic adaptations occurring in this exceptional circumstance in which MRSA harbors simultaneous adaptations related to and impacting antimicrobial resistance and virulence.
Two convergent key motifs emerged as key points of its success: (i) the selection of a compensatory "Silencing regulatory" mutation in the agr-locus, strategic quorum sensing regulator pathway, involved in virulence and antibiotic resistance; (ii) dysregulation in essential primary metabolic pathways associated with a prevalent underexpression trend in DEGs. These co-events determine and allow a balancing of the antimicrobial resistance and virulence adaptations that lead to the setup of increased antimicrobial resistance by decreasing virulence and restoring intrinsic biological costs.
Focusing on antimicrobial-resistance adaptive traits, our findings showed features related to the acquired resistance towards glycopeptides and daptomycin. Resistomics highlighted that MI-nsSNPs in mprF are associated with glycopeptide and daptomycin resistance. In addition, AMR-related transcriptomics revealed overexpression in tcaA, associated with reduced glycopeptide resistance, as well as in dltA and mprF related to daptomycin-glycopeptide resistance. Looking uniquely at glycopeptide resistance, our data indicate that tcaA can have an impact on the level of teicoplanin resistance in terms of MIC values rather than on teicoplanin resistance phenotype acquisition. This observation could seem in contrast with previous findings reporting that tcaA inactivation and tcaRAB deletion determine an increase in teicoplanin resistance [22,23]. These data speculatively implicate that if tcaRAB inactivation increased teicoplanin resistance, on the contrary, an overexpression should decrease teicoplanin resistance. Our data show tcaA overexpression in the teicoplanin-resistant 1-R (TEC MIC 32 mg/L) versus its teicoplanin-susceptible parent 1-S, indicating a positive association with increased teicoplanin resistance. Therefore, the role of tcaA genes in teicoplanin resistance is still not fully clarified.
In S. aureus, daptomycin resistance and reduced susceptibility to vancomycin are often associated with multifactorial mechanisms [24][25][26][27]. Our data confirm dltA and mprF overexpression as the basis of a mechanism of electrostatic repulsion responsible for daptomycin resistance in MRSA, as previously published [24][25][26][27]. dltA and mprf determine an alteration of the surface envelope by alanylation and lysinylation teichoic acids as well as in the membrane phosphatidylglycerol, responsible for increased positive net charge that blocks DAP docking to its target [24,26]. From these new RNA-seq experiments, dltA and mprF overexpression were clearly, statistically significant and differentially expressed in 1-R with respect to 1-S. These data contrasted with our previous findings, likely due to the different methodology used [17].
Furthermore, in VISA DAP-R CA-MRSA, other traits associated, as expected, with other AMR-resistance traits were found. These were acquired resistant SNPs related to rifampin resistance (rpoB), and dysregulation in murF related to D-cycloserine resistance, in norA and norB implicated in FQ-resistance, and tet38 involved in tetracycline-resistance.
Focusing on virulence, virulomics revealed the characteristic huge virulence potential of VISA DAP-R CA-MRSA due to a great pool of virulence-related genes conferring potentially extraordinary abilities of adherence, anti-phagocytosis, immune system evasion and toxigenicity in agreement with previously published data [28].
However, new consideration arises from our data regarding the presence of a stop codon compensatory modulatory mutation in agrA. This allows the CA-MRSA genomic background to shift towards a wide AMR phenotype. The agrA compensatory regulatory mutation, associated or not to the dysregulation in the master regulators (vraRS, sigB, saeR, rot overexpression and walKR, srrA, agrBCA, sarS, graR underexpression), confers the ability to drastically decrease CA-MRSA toxigenicity to adapt the strain to support the antimicrobial resistance burden needed for maintaining extensive AMR. A stop codon, as a "compensatory regulatory mutation" leading to a 170 AA residue truncated AgrA, acts with a "double smart strategy". The former is as an "agr quorum-sensing silencer", the latter is a "virulence-modulator". It is well known that agr-locus positively regulates the production of numerous toxin-coding genes by antisense regulatory RNA-III and negative adhesion [29]. AgrA is the histidine kinase receptor of the agr-system involved in the regulation of S. aureus virulence in response to bacterial cell density. Truncated AgrA is crucial and strategic to block the cascade leading to the activation/repression of the regulatory pathways related both to toxigenicity and antimicrobial resistance. Virulencerelated transcriptomic data confirmed that AgrA non-functionality is concomitant with the underexpression of the agr-locus, hemolysin coding genes (hla,hld,hlgB/C) and enterotoxins (seh, MW1552), and is further supported by the lack of delta-hemolysin production ability.
Moreover, these data confirmed other results that clearly demonstrate that VISA or DAP-R MRSA is typically associated with the decreased agr-functionality frequently showed in agr-II MRSA with increased and extensive AMR [17,[24][25][26][27].
Focusing on biological fitness costs, resistant bacteria pay biological fitness costs for their changes that can diminish their growth rate; our data could appear in contrast with this observation as our VISA DAP-R CA-MRSA showed different strategies to balance the AMR fitness costs [16,30]. VISA DAP-R CA-MRSA faces AMR fitness costs by increasing the generation time compensated with an extension of the lag growth phase that restores the growth kinetics of more susceptible S. aureus. Additionally, two strategies of transcriptomic adaptations were evidenced. VISA DAP-R CA-MRSA was able to balance the biological costs shifting the transcriptional rate as follows: (i) an underexpression trend in differentially expressed genes, including numerous accessory virulence-related genes, i.e., toxin genes (as discussed above), adhesin genes (sdrC, sdrD, eap/map), type VII secretion system genes (esaA), membrane-acting toxin and superantigen genes (spa), exoenzyme genes (lip, splB), immune-modulatory genes coding for the capsular antigens (cap8C/G/P), immunodominant staphylococcal antigen B (isaB) and biofilm production (icaA), confirmed by the lack of biofilm production; (ii) dysregulation of several essential primary metabolic pathways (glycolysis, gluconeogenesis, microbial metabolism in diverse environments, mismatch repair, homologous recombination and DNA replication) and master transcriptional regulators (overexpression in vraRS, sigB, saeR and rot and underexpression in walKR, srrA, agrBCA, sarS and graR).

Conclusions
Our data reveal that convergent evolution, exerted by environmental selective pressure, selects MRSA hybrid variants with extensive antimicrobial resistance towards the last resort antimicrobials, i.e., glycopeptides and daptomycin, within "decreased hypervirulent" CA-MRSA by a balanced pool of adaptations supporting and maintaining the burden of enhanced antimicrobial resistance.
This mechanism is very interesting and to be taken into consideration as it could determine the selection of new hybrid variants that are "extremely dangerous" in the potential severity of infections and in the therapeutical options that could dramatically complicate the future scenario of bacterial infection management.

Bacterial Strains
One S. aureus isogenic strain-pair (1-S and 1-R) of ST-1, agr-III, delta-hemolysin negative, was recovered from a patient hospitalized in an Italian hospital, as previously described [17]. Briefly, a 69-year-old female patient was admitted on 30th of March 2012 to the coronary care unit (CCU) with fever and cardiac decompensation due to mitralic and aortic valve insufficiency. Echocardiography revealed infective endocarditis (IE) on the native mitral and biological aortic valves, and the 1-S MRSA was isolated from blood cultures. Treatment with intravenous (IV) teicoplanin plus gentamycin was begun with initial defervescence. On the 11th day of treatment, the patient was again febrile, and blood cultures yielded a second MRSA (1-R strain) with daptomycin and vancomycin resistance. Vancomycin and gentamycin were then suspended, and IV quinupristin/dalfopristin (Q/D) was prescribed. However, because of the initial unavailability of this drug, IV linezolid was administered for 10 days. Later, Q/D was again available and administered alone for 4 more weeks. During the latter treatment, blood cultures were negative, and the patient was afebrile. However, at the end of this antibiotic treatment, the patient developed a central venous catheter (CVC)-related infection due to Stenotrophomonas maltophilia and had fatal septic shock, as previously described [17]. In detail, the first (1-S) MRSA strain was susceptible to vancomycin, teicoplanin, daptomycin, linezolid, trimethoprim/sulfamethoxazole and rifampin but resistant to cefoxitin (CFX), ampicillin, amoxicillin/clavulanate, fluoroquinolones, erythromycin, clindamycin and tetracycline. On the contrary, the second isogenic MRSA (1-R), isolated under teicoplanin therapy, was characterized by additional resistance to glycopeptides and daptomycin according to the EUCAST guidelines 2021 [31]. The antimicrobial susceptibility and molecular characterization were in part, previously characterized [17] and reported in Tables 2 and 3.

Biological Costs Determined by the Growth Kinetics Test
The 1-S and 1-R MRSA strains were tested for maximum growth rate and length of the lag growth phase. Both MRSA strains were inoculated in Mannitol Salt Agar (MSA) (Oxoid) plates and incubated overnight at 37 • C. One colony of each isolate was resuspended in phosphatebuffered saline and diluted to an optical density of 0.05 at 600 nm. In total, 1 mL of each suspension was then added to 50 mL of BHI broth. Optical density readings at 600 nm were taken every 30 min over 6 h and plotted against time. The generation time was considered the time it takes for bacterial population duplication, and the lag phase duration was at the beginning of the maximum growth rate. Generation time and lag growth phase duration for 1-S and 1-R were compared by the independent samples t-test, and a p-value ≤ 0.05 was considered the cut-off of statistical significance. This experiment was repeated three times, and CA-MRSA MW2 and S. aureus ATCC29213 were used as controls [32].

Biofilm Production
1-S and 1-R MRSA strains were tested for their ability to produce biofilm by a spectrophotometric quantitative assay. Each strain was grown in Tryptone Soy Broth (Oxoid, Basingstoke, UK), with the addition of 0.25% glucose (TSBG), and biofilm production assays were performed in microtiter plates as previously described [33].

Whole-Genome Sequencing
Genomic DNA was extracted using the PureLink Genomic DNA Mini Kit (Invitrogen, Waltham, MA, USA) following the manufacturer's protocol. DNA quality was evaluated by Qubit, and its concentration was determined by Picogreen (Life Technologies, Carlsbad, CA, USA). Whole-genome sequencing (WGS) was performed using the Illumina Mi-Seq sequencing system, using both a paired-end library with 150 bp reads (400 bp average insert size) and a mate-pair library with 250 bp reads (8 kb average insert size). After the sequence data generation, raw reads were processed by Fast QC (v0.11.7) to assess data quality, and the Trimmomatic tool (v0.38) was used to remove sequencing adapters for paired end reads to filter low-quality bases (Q_score < 30) and short reads (<150 bp) as well as for mate pair reads to the process by requiring a minimum base quality of 20 (Phred scale) and a minimum read length of 100 nucleotides to filter out sequences composed only by Ns and to improve the per base score of the mate pair reads. The total number of PE and MP reads is reported with the estimated coverage in Supplementary Table S1. The trimmed reads were used for downstream analysis [27,34].

De-Novo Genome Assembly
De novo genome assembly was performed using SPAdes software (v3.12.0). Reads were initially normalized with khmer 1.3, and then they were error-corrected using the Bayesian Hammer utility of SPAdes. Finally, assembly was performed using the recommended parameters for such Illumina data. The SPAdes software produced a contigs file for each sample, post-assembly controls and metrics were evaluated using Quast software (v4.6.3) and are presented in Supplementary Table S2 [27,34].

Gene Annotation
The assembled contigs were processed using Prokka software (1.14.6) to predict genes and annotate those sequences using a core set of conserved prokaryotic genes [27,34].

Single Nucleotide Variants (SNVs)
For SNVs, genomic re-sequencing was performed from the paired-end library raw reads as already published [27,34]. Briefly, Illumina raw reads were trimmed by the Trimmomatic tool (v0.38), requiring a minimum base quality of 20 (Phred scale) and a minimum read length of 36 nucleotides. Only trimmed reads were included in the downstream analysis. Each sample was aligned by BWA v. 0.7.5 on S. aureus MW2 (BA000033.2) and used as a reference genome. Each .bam file was sorted by Samtools (v.0.1. 19), and duplicate reads were marked using the Picard Mark Duplicates utility. Complex variants, SNVs and indels were detected by "Freebayes" (v.0.9.14), requiring a minimum mapping quality of 25 (Phred scale) and a minimum base quality of 30.
Sequenced reads were properly aligned with the reference genome, 97.28% for 1-S and 97.29% for 1-R.
To select only SNVs present in the VISA DAP-R-CA-MRSA, wg nsSNVs were computationally filtered versus those present in the VSSA DAP-S CA-MRSA parent. All non-synonymous SNPs present in COL-R isolates were confirmed by Sanger sequencing.

Whole-Genome Single Nucleotide Polymorphisms (wgSNPs) Effect Prediction
wgSNPs effect prediction was evaluated using snpEff (v.4.3T). High (HI), low (LI), moderate (MI) or modifier impact (MFI) was assigned according to the criteria previously published [35]. In detail: high impact: the variant is assumed to have disruptive impact, probably causing protein truncation, loss of function or triggering nonsense-mediated decay; low impact: assumed to be mostly harmless or unlikely to change protein behavior; moderate impact: a non-disruptive variant that might change protein effectiveness; modifier impact: usually non-coding variants or variants affecting non-coding genes where predictions are difficult, or there is no evidence of impact [34]. Only high (HI) and moderate (MI) effect wgSNPs were considered and, consequently, described.

Phylogeny and Genomic Epidemiology
The whole-genome sequencing raw data were evaluated by free tools of the Center for Genomic Epidemiology (CGE, http://www.genomicepidemiology.org/, accessed on 1 July 2022) to investigate the genetic and molecular features of the strain-pair. In detail, spaTyper (v1.0) was used to determine the staphylococcal protein A (spa) type of each strain, SCCmecFinder (v1.2) identified the staphylococcal cassette chromosome mec (SCCmec), PlasmidFinder (v2.0) was used for plasmid search, PHAge Search Tool (PHAST) was used considering only the prophage regions detected with a completeness score >90 to detect the prophages, VirulenceFinder (v2.0) was used to identify the virulence factors to define the Virulome and ResFinder (v3.2) was used for the detection of the acquired antimicrobial resistance genes [36][37][38][39][40][41][42]. Analyses were performed with default setting parameters. An aliquot of an overnight culture was diluted 1:50 in 30 mL of brain heart infusion (BHI) in a sterile 50 mL flask (OD 600 nm 0.05) to obtain approximately 5 × 10 5 CFU/mL inoculum for each strain. Cells were grown under shaking at 250 rpm with normal atmospheric conditions at 37 • C and harvested in the exponential growth phase (OD 600 0.5, 2 × 10 8 CFU/mL ∼ 3-4 h). RNA extraction started immediately after cell harvesting to maintain RNA integrity. The cell density was determined by colony counting after plating onto Mueller-Hinton (MH) agar and incubation.

RNA-Seq Libraries
RNA-seq was carried out using the Illumina Mi-seq sequencing platform. To improve RNA-seq data, two replicates using two different libraries were conducted, a Single-End Library with 50 bp reads (SI, Short-Insert Library) and a Paired-end Read Library with 150 bp reads (TS, Tru-Seq Library) and an average insert size of 350/400 bp.

RNA-Extraction
Specific RNA extractions for the Tru-Seq Library and Short-Insert Library preparation were performed according to the specific protocols, as a strategy to optimize the collected RNA-seq data, as previously published [27,34].

Tru-Seq Library Preparation
The total RNA quality was verified using a 2200 TapeStation RNA Screen Tape device (Agilent, Santa Clara, CA, United States), and its concentration was ascertained using an ND-1000 spectrophotometer (NanoDrop, Wilmington, DE, United States). The Agilent TapeStation 2200 system, an automated instrument for nucleic acid gel electrophoresis, assigned RNA integrity number (RIN) values ranging from 1 to 10, with 10 being the highest quality. Only samples with preserved 16S and 23S peaks and RIN values > 8 were used for the library's construction. The RIN values > 8 indicate intact and high-quality RNA samples for downstream applications, as previously published [27,34]. Ribosomal RNA was removed using the Bacteria Ribo-Zero rRNA Removal Kit from 2 µg of RNA. The depleted RNA was used for the Illumina Truseq RNA stranded kit without PolyA enrichment. The obtained libraries were evaluated with high-sensitivity D1000 screen Tape (Agilent Tape Station 2200), and the indexed libraries quantified with the ABI9700 qPCR instrument using the KAPA Library Quantification Kit in triplicates was according to the manufacturer's protocol (Kapa Biosystems, Woburn, MA, United States). From the pooled library, 2 nm final concentrations were used for sequencing with a 150 PE read sequencing module [27,34].

Short-Insert Library Preparation
After ribosomal depletion, sequencing libraries were created using the Illumina mRNA-seq sample preparation kit following the supplier's instructions, except that total RNA was not fragmented, and double-stranded cDNA was size-selected (100-400 bp) to maximize the recovery of small-size RNA. The prepared libraries were valued with high-sensitivity D1000 screen Tape (Agilent Tape Station 2200), as described for the TS Library. The indexed libraries were quantified in triplicate with the ABI7900 qPCR instrument using the KAPA Library Quantification Kit, according to the manufacturer's protocol (Kapa Biosystems, Woburn, MA, United States). From the pooled library, 5 µL at a final concentration of 4 nM were utilized for MiSeq sequencing with an A single-end stranded library with reads of a 50-bp sequencing module [27,34].

Tru-Seq Library Raw Read Post-Processing
After sequence data generation, raw reads were processed using FastQC (v.0.11.2) to assess data quality. The sequenced reads were then trimmed using Trimmomatic (v.0.33.2) to remove only sequencing adapters for PE reads. A minimum base quality of 15 (Phred scale) over a four-base sliding window was required. Only sequences with a length above 36 nucleotides were included in the downstream analysis, and likewise, only trimmed reads were included in the downstream analysis [27,34].

Short-Insert Library Raw Read Post-Processing
After sequence data generation, raw reads were processed using FastQC (v.0.11.2) to assess data quality. Reads were then trimmed using Trimmomatic (v.0.33.2) to remove sequencing adapters for single-end reads, requiring a minimum base quality of 15 (Phred scale) and a minimum read length of 15 nucleotides. Only trimmed reads were included in the downstream analysis [27,34].

Tru-Seq and Short-Insert Read Analysis
TS and SI RNA-seq reads were annotated on S. aureus MW2 (BA000033.2) used as RefGen, as well as transcripts assembled and quantified using Rockhopper (v.2.03) [27,34]. Analyses were run using default parameter settings with verbose output to obtain expression data. Rockhopper normalizes read counts for each sample using the upper quartile gene expression level. Starting from the p-values calculated according to the Anders and Huber approach, differentially expressed genes (DEGs) were assigned by computing q-values ≤ 0.01 based on the Benjamini-Hochberg correction with a false discovery rate of <1%. In addition, Rockhopper is a versatile tool using biological replicates when available and surrogate replicates when biological replicates for two different conditions are unavailable, considering the two conditions under investigation as surrogate replicates for each other [27,34].

DAVID Enrichment Analysis
The online tool DAVID (v.6.8) (http://david.abcc.ncifcrf.gov/ accessed on 1 July 2022) was used to detect affected pathways among DEGs. The gene lists of the strain-pair, grouped according to over-and underexpressed genes, were uploaded as Official Gene Symbols of the S. aureus MW2 reference genome, automatically selecting the list type (Gene list) of S. aureus. The Functional Annotation Chart was obtained using an EASE score threshold ≤0.5 and a minimum count number of four genes. The DAVID Functional Categories were investigated by the KEGG pathway and PANTHER Classification System [43] and grouped in annotation clusters refined of the same genes.

Real-Time qPCR validation
To validate RNA-seq data, real-time qPCRs for a set of characterizing transcripts, murF, dltA, mprF, hld, hla, spa, agrA, icaA, and sdrD, were carried out in the same RNA-seq growth phase. Primer lists are provided in Tab.S3; gyrB was used as a normalizer gene. Real-time qPCRs and statistical analyses were carried out as previously published [31,32].

DNA-Sequencing and RNA-Sequencing Data Accession Number
The genomic and transcriptomic reads were deposited in the National Center for Biotechnology Information (NCBI) Genome database in the Sequence Read Archive (SRA) under the BioProject PRJNA860577 with the Biosample Genomic Paired End and RNA-seq raw sequences n • SAMN29849119, SAMN29849120 and with Biosample Genomic mate-pair draft sequences n • SAMN30428603, SAMN30428604.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/antibiotics11091159/s1, Table S1: Total number of paired end and mate pair reads and Estimated Coverage Generated by the Trimming after the sequencing process, Table S2A,B: Summary  Table of Post-Assembly Statistics and Metrics generated using the Quast software, Table S3. Primer set used in real time qPCR Validation.
Author Contributions: All authors contributed to the manuscript. V.C. conceived and designed the study. V.C., R.S., A.Z., E.A. and F.L.V. performed the genomics, transcriptomics and bioinformatics analysis. G.P. contributed to the bioinformatics analysis. S.S. critically revised the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This study was supported by a research grant, PRIN2017SFBFER, from MIUR, Italy. "Starting Grant" of the "Piano di incentivi per la ricerca di Ateneo" Prot. n. 253556 del 16 marzo 2021, University of Catania.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Genomic and transcriptomic reads can be found in National Center for Biotechnology Information (NCBI) Genome database in the Sequence Read Archive (SRA) under the BioProject PRJNA860577 with the Biosample Genomic Paired End and RNA-seq raw sequences n • SAMN29849119, SAMN29849120 and with Biosample Genomic mate-pair draft sequences n • SAMN30428603, SAMN30428604.