ORIGINAL RESEARCH article
Genomic and Long-Term Transcriptomic Imprints Related to the Daptomycin Mechanism of Action Occurring in Daptomycin- and Methicillin-Resistant Staphylococcus aureus Under Daptomycin Exposure
- 1Department of Biomedical and Biotechnological Sciences, University of Catania, Catania, Italy
- 2Department of Clinical and Experimental Medicine, University of Catania, Catania, Italy
Daptomycin (DAP) is one of the last-resort treatments for heterogeneous vancomycin-intermediate Staphylococcus aureus (hVISA) and vancomycin-intermediate S. aureus (VISA) infections. DAP resistance (DAP-R) is multifactorial and mainly related to cell-envelope modifications caused by single-nucleotide polymorphisms and/or modulation mechanisms of transcription emerging as result of a self-defense process in response to DAP exposure. Nevertheless, the role of these adaptations remains unclear. We aim to investigate the comparative genomics and late post-exponential growth-phase transcriptomics of two DAP-resistant/DAP-susceptible (DAPR/S) methicillin-resistant S. aureus (MRSA) clinical strain pairs to focalize the genomic and long-term transcriptomic fingerprinting and adaptations related to the DAP mechanism of action acquired in vivo under DAP pressure using Illumina whole-genome sequencing (WGS), RNA-seq, bioinformatics, and real-time qPCR validation. Comparative genomics revealed that membrane protein and transcriptional regulator coding genes emerged as shared functional coding-gene clusters harboring mutational events related to the DAP-R onset in a strain-dependent manner. Pairwise transcriptomic enrichment analysis highlighted common and strain pair-dependent Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, whereas DAPR/S double-pair cross-filtering returned 53 differentially expressed genes (DEGs). A multifactorial long-term transcriptomic-network characterized DAPR MRSA includes alterations in (i) peptidoglycan biosynthesis, cell division, and cell-membrane (CM) organization genes, as well as a cidB/lytS autolysin genes; (ii) ldh2 involved in fermentative metabolism; (iii) CM-potential perturbation genes; and (iv) oxidative and heat/cold stress response-related genes. Moreover, a D-alanyl–D-alanine decrease in cell-wall muropeptide characterized DAP/glycopeptide cross-reduced susceptibility mechanisms in DAPR MRSA. Our data provide a snapshot of DAPR MRSA genomic and long-term transcriptome signatures related to the DAP mechanism of action (MOA) evidencing that a complex network of genomic changes and transcriptomic adaptations is required to acquire DAP-R.
Methicillin-resistant Staphylococcus aureus (MRSA) remains one of the major multidrug-resistant pathogens responsible for severe infections with high mortality rates (Stefani et al., 2015; Yoon et al., 2016; Britt et al., 2017). The subset of daptomycin-resistant (DAPR) S. aureus is of particular concern for the cumulative non-reversible metabolic changes demonstrated in resistant strains and for the difficulty in the treatment of severe infections (Bayer et al., 2013; Reed et al., 2019).
The DAP mechanism of action (DAP-MOA) has not been fully elucidated yet. The model of DAP-MOA proposes that DAP binds the cytoplasmic membrane leading to its permeabilization and depolarization caused by a loss of cytoplasm potassium ions (Allen et al., 1987). This mechanism, which can account for DAP’s bactericidal effect, would correlate with several changes in the membrane components, for example, with the level of phosphatidylglycerol. The complex mechanism of DAP resistance (DAP-R) was reviewed by Miller et al. (2016). Briefly, S. aureus can adopt different strategies to resist DAP. The first is to alter the net cell-surface charge, preventing the positively charged DAP-calcium binding to the cell envelope by electrostatic repulsion, mainly due to dlt over-expression—increasing the alanylation rate of the wall teichoic acids (Yang et al., 2009; Cafiso et al., 2014)—and to mprF mutations increasing the amount of the positively charged lysyl-phosphatidylglycerol on the outer membrane, conferring a “gain-of-function” and positively increasing the cell-envelope charge (Rubio et al., 2011; Bayer et al., 2015; Ernst et al., 2018; Sabat et al., 2018; Ernst and Peschel, 2019). The second is the alteration of membrane phospholipid composition, with decreased phosphatidylglycerol amount, or changes in membrane fluidity interfering with DAP binding and oligomerization (Jones et al., 2008; Kilelee et al., 2010; Mishra et al., 2011; Zhang et al., 2014). A third mechanism is the involvement of global regulatory genes modulating cell-envelope stress and maintenance, affecting the expression of the cell-wall (CW) “stimulon” (Utaida et al., 2003; Cafiso et al., 2012). Specifically, VraSR operon (Muthaiyan et al., 2008; Mehta et al., 2012; Sabat et al., 2018; Taglialegna et al., 2019), orthologous to the Bacillus subtilis LiaSR (Jordan et al., 2006), was involved and up-regulated by vancomycin and DAP exposure, and associated with CW biosynthesis via transcription of PBP2 (penicillin-binding protein 2), tagA (wall teichoic acids-synthesis), prsA (a-chaperone), and murZ (UDP-N-acetylglucosamine-enolpyruvyl transferase) (Kuroda et al., 2003; Mwangi et al., 2007). YycFG (also named WalKR), among the two-component regulatory systems (TCRSs), was implicated in the control of the peptidoglycan biosynthesis through the regulation of LytM and AtlA (Dubrac and Msadek, 2004; Fukushima et al., 2011). A YycFG down-regulation was associated with a decreased peptidoglycan turnover, augmented cross-linking, increased glycan chain length, and resistance to lysis by Triton X-100 (Dubrac et al., 2007). YycG senses changes in membrane fluidity and responds by adjusting CW cross-linking, compensating the stresses caused by osmotic pressure (Türck and Bierbaum, 2012). In DAPR strains, yycFG can accumulate mutations (Friedman et al., 2006; Howden et al., 2011), impairing the YycFG function and altering CW homeostasis to survive the DAP action, that is, lowering peptidoglycan turnover and increasing cross-linking. Mutations in the RNA polymerase rpoB and rpoC subunits were associated with DAP-R (Friedman et al., 2006). A621E mutation in RpoB was associated with an increased expression of the dlt-operon and correlated with an increase in positive cell-surface charge, whereas A621E and A477D substitutions were linked to an increased CW biosynthesis and thickness (Cui et al., 2010; Bæk et al., 2015).
This investigation aimed to face the analysis of DAP-R mechanisms in MRSA strains on the basis of the acquired knowledge and considering new aspects. First, the “already known” DAP-R mechanisms (dlt over-expression, mprF mutations, and increased net positive cell-surface charge) do not always correlate with changes in DAP minimum inhibitory concentrations (MICs) (Mishra et al., 2014); second, during an infection, the bacteria elapse long periods in limited/arrested growth as in the late growth phases (Kolter et al., 1993); third, DAP is also active on non-dividing and not metabolically active bacteria (Mascio et al., 2007; McCall et al., 2019).
On this rationale, our study investigated genomics and, for the first time, the comparative late growth-phase transcriptomics of clinical DAPR versus DAPS isogenic MRSA by using whole-genome sequencing (WGS), RNA-seq, and bioinformatics to reveal their genomic traits as well as the long-term fingerprinting and adaptations related to DAP-MOA occurring in the acquisition of DAP-R by DAPS isolates.
These new findings could advance the knowledge in the field, revealing the genomic and long-term transcriptomic key features acquired in vivo by the DAPR versus their parental isogenic DAPS ancestor under DAP exposure and maintained after DAP removal, becoming stable traits characterizing DAPR MRSA.
Our data showed that DAP-MOA-related strain-dependent non-synonymous single-nucleotide polymorphisms (nsSNPs) in cell membrane (CM), structural/functional protein, and transcriptional regulator coding genes, combined with a complex transcriptional network, appeared as key factors characterizing DAPR MRSA versus DAPS parents.
Materials and Methods
Two S. aureus epidemiologically unrelated isogenic strain pairs of DAPS (1A, 3A) and DAPR MRSA (1C, 3B) isolated under DAP therapy, previously collected from patients hospitalized in two Italian hospitals, were investigated. The two isogenic strain pairs were provided by two different Italian hospitals (Cafiso et al., 2014); consequently, an ethical approval was not necessary as per institutional and national guidelines and regulations. Specifically, the 1A and 1C strains were isolated in Palermo by a skin infection, whereas the 3A and 3B strains were isolated by a soft tissue infection in a patient in Bergamo. Each pair included an initial pre-DAP therapy strain (1A and 3A) and its isogenic isolate after development of DAP-R during DAP administration. Each strain pair was previously assigned to the same pulsotype, multilocus sequence typing (MLST) type, SCCmec type, and agr group (Cafiso et al., 2014). The D-alanylated wall teichoic acid content, the lysinylated-phosphatidylglyderol amount, the CM fluidity, and the susceptibility to prototypic cationic host defense peptides of platelet and leukocyte origins were previously determined (Mishra et al., 2014). Finally, exponential or stationary growth-phase expression and SNPs in the “already DAP-R associated genes,” mprF and dltA, were previously reported (Table 1) (Cafiso et al., 2014; Mishra et al., 2014).
Whole-genome sequencing was performed using the Illumina MiSeq sequencing system.
DNA Extraction and Illumina Shotgun Paired-End Library Preparation
Genomic DNA was extracted using PureLink Genomic DNA Mini Kit (Invitrogen) according to the manufacturer’s protocol. The quality of the DNA was assessed using Qubit and its concentration ascertained by Picogreen (Life Technologies). Paired-end (PE) reads libraries were prepared by Nextera XT DNA Library Prep Kit (Illumina, San Diego, CA, United States) following the manufacturer’s protocol, and their quality evaluation was performed as previously published (Cafiso et al., 2019). The indexed libraries were quantified as previously published (Cafiso et al., 2019) pooled at a final concentration of 2 nM and used for Illumina MiSeq sequencing with a PE 300 (2 × 300 bp). Raw reads were processed using FastQC (v.0.11.7) to assess data quality. Cutadapter tool (v.1.16), implemented in Python (v.3.5.2), was used to remove residual PCR primer and to filter low-quality bases (Q_score < 30) and short reads (<150 bp). The filtered trimmed reads were included in the downstream analysis. The total number of PE reads was reported with the estimated coverage in Supplementary Figure S1.
De novo Genome Assembly
De novo genome assembly was performed using SPAdes software (v3.12.0), producing a contig file for each sample. Post-assembly controls and metrics were generated using Quast (v.4.6.3). Supplementary Figure S2 shows de novo Genome Assembly Report.
Genome alignments were per performed by Mauve (Multiple Genome Alignments) (v.2.4.0).
The assembled scaffolds were processed using Prokka (v.1.12) software.
DNA-Sequencing Data Accession Number
The genomic reads were deposited in the National Center for Biotechnology Information (NCBI) Genome database in the Sequence Read Archive (SRA) under study accession no. SRP166981 (BioProject: PRJNA498510).
Single-nucleotide polymorphisms calls were carried out from the PE library raw reads as previously published (Cafiso et al., 2019). Briefly, Illumina raw reads were trimmed using Trimmomatic (v.0.38), requiring a minimum base quality of 20 (Phred scale) and a minimum read length of 36 nucleotides, and aligned by “Bwa mem” (v.0.7.17). For detection of SNPs, and insertions and deletions (indels), each.bam file was sorted using Samtools (v.1.9), and duplicated reads were marked using the Picard Mark Duplicates utility. Complex variants, SNPs, and indels were detected using “Freebayes” (v.1.2.0), which required a minimum mapping quality of 30 (Phred scale) and a minimum base quality of 20. CSI phylogeny, a Center for Genomic Epidemiology service1, was used to identify the closest relationships between the strain pairs and four different reference genomes (RefGen) deposited in GenBank, that is, S. aureus NCTC8325 (CP000253.1), S. aureus USA 300 (CP000255.1), S. aureus MU50 (BA000017.4), and S. aureus ST398 (NC_017333.1).
Staphylococcus aureus NCTC8325 (CP000253.1) was selected as common RefGen for SNP mapping in agreement to the data obtained. The graphical output was a Phylogenetic tree shown in Supplementary Figure S3.
The majority of the sequenced reads was properly aligned with the corresponding reference genome, in detail, 90.98% for 1A, 91.08% for 1C, 96.45% for 3A, and 95.81% for 3B.
To select SNPs present in the DAPR strains, all SNPs were computationally filtered in both pairwise mode and double-cross filtering to find out if only nsSNPs present simultaneously in both two DAPR strains and the DAPS parents.
Whole-genome sequencing raw data were analyzed to investigate the genomic epidemiology with ResFinder (v.3.2) and Point Finder (v.3.1.0) services for the detection of the acquired antimicrobial resistance (AMR) genes and to detect the known nsSNPs related to AMR using a 98% threshold for nucleotide sequence identity and 60% for the minimum length coverage (Zankari et al., 2012). MLST was genomically checked using the MLST (v2.0) (Larsen et al., 2012), spa-typing was confirmed by spaTyper (v.1.0) (Bartels et al., 2014), the identification of acquired virulence genes was performed by VirulenceFinder (v.2.0) (Joensen et al., 2014), and the plasmid strain profile was investigated by PlasmidFinder server (v.2.1) (Carattoli et al., 2014). All tools are available on the Center for Genomic Epidemiology website1. The genomic epidemiology of strain pairs was also analyzed by Nullarbor pipeline2.
Core Genome Single-Nucleotide Polymorphisms
The core genome SNP (cgSNP) detection shared among all strains was computationally carried out by Nullarbor2.
Genomic Single-Nucleotide Polymorphism Effect Prediction
The prediction of the whole-genome SNPs (SNPome) effects was performed by SnpEff (v.4.3T), a genomic variant annotation and functional effect prediction toolbox. High (HI), low (LI), moderate (MI), and modifier impacting (MFI) effects were assigned according to the criteria previously published and in use in the tool4 (Cingolani et al., 2012). High impact: the variant is assumed to have disruptive impact in the protein, probably causing protein truncation and loss of function or triggering nonsense mediated decay. Low impact: the variant is assumed to be mostly harmless or unlikely to change protein behavior. Moderate impact: the variant is a non-disruptive variant that might change protein effectiveness. Modifier impact: the variant is a usually non-coding variant or variants affecting non-coding genes where predictions are difficult or there is no evidence of impact.
RNA-Seq Bacterial Cultures
For the construction of all RNA-seq libraries, an aliquot of an overnight culture was diluted 1:50 in 30 ml of brain heart infusion (BHI) in a sterile 50-ml flask (OD600 nm 0.05) to obtain an approximately 5 × 105 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 late post-exponential growth phase (1.5–4 × 1010 CFU/ml, ∼18 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. Growth curves of the four DAP-free MRSA strains were shown in the Supplementary Figure S4.
RNA-seq was performed using the Illumina MiSeq sequencing system, and biological replicates were performed using two different libraries, a single-end library with 50-bp reads [Short-Insert Library (SI)] and a PE read library with 150-bp reads [Tru-Seq Library (TS)] starting from two different RNA extractions according to the following protocols as a strategy to optimize the collected RNA-seq data (Pereira et al., 2017; Cafiso et al., 2019).
Tru-Seq Library RNA Extraction
RNA was extracted using the NucleoSpin RNA (Macherey-Nagel, Dueren, Germany) kit following the manufacturer’s protocol with minor modifications. Pellets were lysed by the bead-beating procedure with 50 μl of RA1 buffer. Then 3.5 μl of β-mercaptoethanol was added, and the lysates were filtered through NucleoSpin Filters. The RNA-binding condition was adjusted by adding 70% ET-OH to the lysates, and the RNA was extracted with TRIzol reagent (Invitrogen) according to the manufacturer’s protocol. 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 (Fleige and Pfaffl, 2006). 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 (Cafiso et al., 2019).
Short-Insert Library RNA Extraction
RNA was extracted with TRIzol reagent (Invitrogen) according to the manufacturer’s protocol. After ribosomal depletion, sequencing libraries were prepared 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 evaluated 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 was used for MiSeq sequencing with an A single-end stranded library with reads of 50-bp sequencing module (Cafiso et al., 2019).
Tru-Seq Library Preparation and Sequencing
The samples were processed using the Illumina MiSeq, using an A PE library with reads of 150 bp and average insert size of 350/400 bp. 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. Similarly, only trimmed reads were included in the downstream analysis (Cafiso et al., 2019).
Short-Insert Library Preparation and Sequencing
The samples were processed using the Illumina MiSeq technology with an A single-end stranded library with reads of 50 bp. 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 (Cafiso et al., 2019).
Tru-Seq and Short-Insert Library Analysis
TS and SI RNA-seq reads were annotated on S. aureus NCTC 8325 RefGen (CP000253.1), as well as transcripts assembled and quantified using Rockhopper (v.2.03) (McClure et al., 2013; Tjaden, 2015), specifically designed for bacterial gene structures and transcriptomes. 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 on the basis of 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 (McClure et al., 2013).
Finally, filtering analyses were computationally carried out for sorting, first, the DEGs in the DAPR strains versus their DAPS parental strain and, then, a selection of only those present contemporarily in both DAPR isolates showing the same expression (under-expression or over-expression) (Cafiso et al., 2019).
Determination of the Affected Pathways
The online tool DAVID (v.6.8) was selected to identify the affected pathways among the DEGs. The gene lists were uploaded as Official Gene Symbols of the S. aureus NCTC 8325 RefGen and automatically selecting the list type (Gene list) of S. aureus. The Functional Annotation Chart was visualized using the p-value threshold of 0.01 and a minimum count number of four genes. The information on the affected pathways was obtained from Kyoto Encyclopedia of Genes and Genomes (KEGG) within the analysis in DAVID, using the mentioned thresholds. The affected biological processes were obtained from the Gene Ontology (GO) Consortium within the analysis in DAVID (Cafiso et al., 2019).
RNA Functional Categories
Functional categories of DEGs were investigated by BLAST, PANTHER Classification System, GO Consortium, ExPASy (STRING), and KEGG pathway (Cafiso et al., 2019).
RNA-Sequencing Data Accession Number
RNA-seq raw reads were deposited at the Gene Expression Omnibus (GEO) database under study accession no. GSE121797 (BioProject: PRJNA498510).
Downstream Statistical Analysis
An heatmap of the entire set of DEGs and principal component analysis (PCA) on the five discovered DEG functional categories, in DAPR MRSA versus their DAPS parents, were performed by XLSTAT 2020 Excel data analysis add-on (v.2020.1.3) using default setting.
To validate RNA-seq data, real-time qPCRs were performed on the most DAPR-relevant transcripts (SAOUHSC_02317, SAOUHSC_00022, SAOUHSC_00486, SAOUHSC_02922, SAOUHSC_01806, SAOUHSC_01334, and SAOUHSC_00545) using RNAs of three experiments performed under the RNA-seq growth conditions. Real-time qPCRs were carried out following published protocols (Cafiso et al., 2012, 2014); the primer sequences are shown in Supplementary Table S1 and used with the following thermal profile: 95°C for 3 min; 35 cycles at 95°C for 10 s, 60°C for 30 s, and 72°C for 45 s; and a final cycle at 95°C for 1 min, 60°C for 30 s, and 95°C for 30 s. gyrB was the normalizer. Statistical analyses were conducted as previously described (Pfaffl et al., 2002).
L-Lactic Acid Concentration Assay
Single strains were grown under shaking at 250 rpm in standard atmospheric conditions at 37°C and harvested in the late post-exponential growth phase. The lactate conversion was measured by enzymatic lactate oxidation into pyruvate using a D-lactic acid/L-lactic acid commercial test (R-Biopharm, Germany). L-Lactic acid concentrations were calculated using the D-lactic acid/L-lactic acid ultraviolet method according to the manufacturer’s instructions (R-Biopharm), which is based on the conversion of D-lactate or L-lactate by D-lactate dehydrogenase and L-lactate dehydrogenase, respectively, to pyruvate and NADH. The L-lactic acid concentrations were expressed in g/L. Three biological and technical replicates were performed for each strain. A statistical significance was attributed in Student’s t-test with p-values < 0.05.
Core Genome Single-Nucleotide Polymorphism Phylogeny
Pairwise cgSNP distance analysis showed a very close phylogeny of the DAPR strains respect to the DAPS parents. Thirty-six cgSNP differences were found between the 1A and 1C, whereas 41 cgSNPs were found between the 3A versus 3B isolated during a period of both 4 and 6 months under DAP therapy, respectively.
The genomic epidemiology showed the membership of the 1A/C strain pair to the ST398-SCCmecIVa-t1939-agrI both harboring the plasmid pDLK1 carrying ermC, together with the pS194 harboring the str gene (Table 1) only in 1C, both conferring aminoglycoside resistance. The 3A/B strain pair was confirmed to be ST8-SCCmec-IVc-t008-agrI having the SAP063A plasmid carrying blaZ determining β-lactam resistance, and with the 3B having in addition the pDLK1 carrying ermC (Table 1).
Integrating the Center for Genomic Epidemiology and Nullarbor outputs, the strain pair resistomes (acquired AMR genes and known chromosomal SNPs conferring antimicrobial resistance) and virulomes (chromosomal and acquired virulence genes) are outlined as shown in Table 1.
Mauve Genome alignments of the single DAPR MRSA genome versus its DAPS parent showed the lack of new genetic elements associable to DAP-R. The unique additional genetic traits present in DAPR MRSA strains were the two plasmids, that is, pS194 in 1C and pDLK1 in 3B.
SNPome and Possibly Related to Daptomycin Resistance Acquisition Non-synonymous Single-Nucleotide Polymorphisms
DAPR/S strain-pair SNPome carried out versus S. aureus NCTC8325 RefGen recorded a variant rate of 64 in the DAPR/S 1A/C strain pair with the presence of 43,720 variants in 1C and 43,947 in 1A; otherwise, a variant rate of 893 was found in 3B and 592 in 3A, with the presence of 3,157 variants and 4,762 variants, respectively.
In the DAPR/S 1 A/C strain pair, according to the criteria of SnpEff tool, 0.924% mutations were putatively predicted with HI on their products; 57.118% with LI; 23.665% with MI; and 18.293% with MFI, in 1C. Among these, 23.376% were missense (M), 0.25% were nonsense (NS), and 76.374% were silent (S). In 1A, similarly, 0.978% had HI, 56.833% had LI, 23.759% had MI, and 18.431% had MFI, with 23.701% of M mutations, 0.314% of NS mutations, and 75.984% of S mutations.
No shared HI and MI acquired nsSNPs in DAP-MOA-related targets were found in both DAPR and DAPS MRSA strain pairs (Figures 1, 2); however, high stringent DAVID enrichment analysis evidenced HI and/or MI DAP-MOA-related nsSNPs in the same functional protein clusters, with recovered coding genes and nsSNPs predominantly different between the two DAPR/S MRSA strain pairs.
Figure 1. High (HI) and MI nsSNP overview of 1C versus 1A. HI, high impact; MI, moderate impact; nsSNP, non-synonymous single-nucleotide polymorphisms.
Figure 2. High and MI nsSNP overview of 3B versus 3A. HI, high impact; MI, moderate impact; nsSNP, non-synonymous single-nucleotide polymorphisms.
In details, HI and/or MI DAP-MOA-related nsSNPs were detected in two UP_SEQ_FEATURE clusters of targets putatively related to DAP-MOA, that is, membrane proteins (including transmembrane proteins, lipoproteins, and two-component histidine kinase systems and transporters) and transcriptional regulators. In addition, only in DAPR 3B versus DAPS 3A MRSA were HI DAP-MOA-related nsSNPs also recovered in the UP_SEQ_FEATURE Signal Peptide Protein cluster including gene coding proteins involved in CW turnover and teichoic acid metabolism. Furthermore, in both DAPR and DAPS MRSA strain pairs, MI DAP-MOA-related nsSNPs were also detected in a miscellaneous of gene coding proteins putatively related to DAP-MOA including proteins related to autolysis, cell division, peptidoglycan transport and cell shape, diacylglycerol metabolism, phosphatidylglycerol modifications (MprF), proteins involved the cardiolipin biosynthetic process (Cls1), and lipoteichoic acid biosynthesis (DltB) (Supplementary Table S2).
DAVID Pairwise Enrichment Analysis
Integrating data obtained from the two libraries, DAVID enrichment pairwise analysis (p-value ≤ 0.05) (Huang et al., 2009a, b) of the over-expressed DEGs evidenced the enrichment of only the ABC transporter KEGG pathway in both DAPR/S strain pairs (Supplementary Tables S3, S4). In addition, among the under-expressed DEGs, we found several enriched KEGG pathways responsible for ABC transporter, glycolysis and pyruvate metabolism, purine and propanoate metabolism, phosphotransferase system (PTS), and two-component regulatory systems. Among the TCRSs, SaeR/S (SAOUHSC_00714/SAOUHSC_00715) was found in both pairs and WalK (SAOUHSC_00021) only in the 3A/B pair (Supplementary Table S3). Diverse strain-dependent enriched over-expressed/under-expressed KEGG pathways were observed in the single 1A/C and 3A/B pair as shown in Supplementary Tables S3, S4.
Cross-Filtered Double-Pair Differentially Expressed Genes
Stringently double-pair cross-refiltering of the RNA-seq data—to define the signatures exclusively present in either DAPR MRSA strains or their DAPS parents, and with the same expression profiling in the same library (SI or TS)—highlighted the co-occurrence of 53 statistically significant DEGs: nine over-expressed and 44 under-expressed protein-coding genes were found (Figures 3, 4 and Supplementary Table S5).
Figure 3. Heatmap of DEGs related to DAP-MOA and accessory traits maintained after DAP resistance onset in DAPR MRSA. COG/GO cluster of genes is shown on the top of the figure, and the transcript ID and gene name are shown on the bottom. DAPS (1A, 3A) and DAPR MRSA (1C, 3B). DEGs, differentially expressed genes; DAP, daptomycin; COG, clusters of Orthologous Group; GO, Gene Ontology; DAPR, daptomycin resistant; DAPS, daptomycin susceptible.
Figure 4. Principal component analysis of the five discovered DEG functional categories versus the DAPR/DAPS MRSA. PCA plots show the relations between strains and genes in the different clusters. PCA, principal component analysis; DEG, differentially expressed gene; DAPR, daptomycin resistant; DAPS, daptomycin susceptible; MRSA, methicillin-resistant Staphylococcus aureus.
According to the GO numbers and Clusters of Orthologous Group (COG) groups, we categorized the characterizing DAPR DEGs in different clusters, as follows.
Cell division, cell wall, and cell membrane
Three DEGs implicated in peptidoglycan biosynthesis, two DEGs for the CW autolysis, and five DEGs involved in cell division and CM structure were found.
For peptidoglycan biosynthesis, we found under-expression in the UDP-N-acetylmuramoyl-tripeptide-D-alanyl–D-alanine ligase murF, in SAOUHSC_00751 hypothetical protein encoding gene (neighboring and consequently predicted as functionally related to the UDP-N-acetylenolpyruvoylglucosamine reductase MurB for the CW formation STRINGscore: 0.808), and in YycFG (WalKR) negative regulator, yycH. In CW autolysis, over-expression was shown in cidB autolysin and under-expression in the sensor histidine kinase lytS.
As concerns cell division and CM structure, under-expression was displayed in ftsH, incorporating the PBPs into the CM, in SAOUHSC_01919, SAOUHSC_02376, and SAOUHSC_02391 uncharacterized proteins located—according to the GO Cellular Component term (GO-CC): 0016021—in the CM, whereas over-expression was recorded only in SAOUHSC_03035 integral membrane protein.
A complex network involving differential expression was shown in metabolic genes. As regards the generation of precursor metabolites and energy, L-lactate dehydrogenase-2 encoding gene ldh2 involved in lactic fermentation was over-expressed, whereas the pyruvate kinase gene pyk involved in glycolysis and the protoporphyrinogen-oxidase SAOUHSC_01960 (porphyrin biosynthesis) were under-expressed. In the cofactor metabolism, under-expression was recorded in 2-succinyl-6-hydroxy-2,4-cyclohexadiene-1-carboxylate synthase menH.
Regarding the compound metabolism, under-expression was found in cysteine-desulfurase SAOUHSC_01727 for sulfur-compound metabolism; in the porphobilinogen-deaminase hemC related to the nitrogen-compound metabolism; in SAOUHSC_02396 Cof-like hydrolase, SAOUHSC_01055 inositol-monophosphatase, and SAOUHSC_00847 ABC transporter involved in the metabolism of phosphate-containing compounds; in SAOUHSC_00861 lipoyl-synthase lipA for the coenzyme and cofactor metabolism; and in the molybdenum cofactor biosynthesis moaA.
With regard to nucleotide metabolism, phosphopentomutase deoB, GMP reductase guaC, carbamoyl-phosphate synthase small chain carA, and uridylate kinase pyrH were under-expressed. With regard to the amino acid and lipid catabolism, our data evidenced under-expression in alanine dehydrogenase ald2 and lipase-1 lipA. In the cellular protein modification process, under-expression was observed in nrdI coding the Class Ib ribonucleotide reductase stimulatory protein.
Nucleic acid metabolism
Under-expression was found in the SAOUHSC_01693 DNA-binding protein for DNA repair, and in the SAOUHSC_01690 of the DNA polymerase III complex and the ATP-dependent DNA-elicase, pcrA. In the RNA transcription and regulation, our data showed under-expression in DNA-binding protein HU coding gene functionally related to ClpC (STRINGscore: 0.555), in transcription anti-termination nusB, in SAOUHSC_00031 tRNA-dihydrouridine synthase, and in the aspartyl/glutamyl-tRNA amidotransferase subunit B gatB (tRNA metabolism). Examining ribosome biogenesis and maturation, we found over-expression in ATPase, whereas under-expression was found in rsgA (ribosome biogenesis). In SOS response, we observed over-expression in SAOUHSC_01334 hypothetical protein functionally related to LexA repressor (STRINGscore: 0.614) (DNA-damage response).
Referring to the stress response, SAOUHSC_00304 luciferase-like monooxygenase (oxidative stress response) and SAOUHSC_00406 uncharacterized protein [previously annotated as poly(3-hydroxybutyrate) (PHB) depolymerase family protein] determining the use of PHB as a carbon source and energy storage in starvation conditions, were over-expressed. dnaJ (prevention of stress-denatured protein aggregation in response to hyperosmotic condition and heat shock) and the SAOUHSC_00204 globin domain protein, nitric oxide dioxygenase (nitrosative stress response), were under-expressed.
Under-expression was observed in SAOUHSC_00099, SAOUHSC_00137, and SAOUHSC_02700 transporters; in the SAOUHSC_00888 Na+/H+ antiporter mnhB1; and in SAOUHSC_00634 ABC metal ion transporter involved in cell adhesion. Over-expression was only found in SAOUHSC_01387 inorganic phosphate transmembrane transporter.
Results from carbohydrate transporters showed under-expression in SAOUHSC_00177 maltose ABC permease, in the SAOUHSC_02400 putative mannitol-specific PTS component, and in the SAOUHSC_00158 N-acetylmuramic acid transporter belonging to the phosphotransferase system.
Under-expression was shown in the serine-aspartate repeat-containing sdrD, a staphylococcal virulence factor.
Over-expression was observed in SAOUHSC_00826 encoding a conserved uncharacterized protein.
RNA-Seq Data Validation
RNA-seq data validation performed by real-time qPCR on the most DAPR-relevant genes confirmed the murF (SAOUHSC_ 02317), yycH (SAOUHSC_00022), ftsH (SAOUHSC_00486), ldh2 (SAOUHSC_02922), pyk (SAOUHSC_01806), SOS response protein related to LexA (SAOUHSC_01334), and sdrD (SAOUHSC_00545) expression profiles detected in RNA-seq data in both DAPR S. aureus and DAPS parents (Supplementary Figure S5).
Lactic Acid Quantification Assay
To biologically validate the bioinformatic prediction of the over-expression in L-lactate dehydrogenase-2, quantification assays of the amount of lactic acid were performed, showing a statistically significant increased concentration of L-Lactate in both DAPR MRSA and their DAPS parental strains (Figure 5).
Figure 5. L-Lactic acid quantification. L-Lactic acid production by the four strains was calculated using the L-lactic acid ultraviolet method according to the manufacturer’s instructions (D-lactic acid/L-lactic acid kit, R-Biopharm), based on conversion of L-lactate by L-lactate dehydrogenase to pyruvate and NADH. L-Lactate production differences with p-values < 0.05, obtained by Student’s t-test, were considered statistically significant. * is for the statistical significance.
In S. aureus, DAP-R is a complex and multifactorial mechanism shaped by the co-occurrence of intrinsic, acquired, and adaptive determinants that, harmonizing in concert, confer resistance. We investigated two clinical epidemiologically unrelated DAPR/S MRSA strain pairs isolated from different patients in two different Italian hospitals. The two patients failed glycopeptide therapy and were then treated with DAP (Cafiso et al., 2014).
Previous investigations on the same DAPR/S MRSA strain pairs characterized some traits related to DAP-R, including a dysregulation in the two key determinants of net positive surface charge, dltA and mprF, both during exponential and stationary growth phases; a significant increase in the D-alanylated wall teichoic acids amount correlating with DltA gain-in-function; a heightened elaboration of lysyl-phosphatidylglycerol reflecting MprF gain-in-function; an increased CM fluidity; a strain-dependent CM fatty acid perturbation due to an increase in the anteiso-branch chain species corresponding to a reduction in the major iso-branched chain and saturated fatty acids (SFAs); and a reduced susceptibility to prototypic cationic host defense peptides of platelet and leukocyte origins (Cafiso et al., 2014; Mishra et al., 2014; Boudjemaa et al., 2018).
In the present study, genome-wide CSI phylogeny and the cgSNP analysis of the same strain pairs revealed the close relation of the DAPR/S MRSA isolated from the same patient under DAP therapy administered in hospital during a period of several months. Concomitantly, genomic epidemiology confirmed their previously published molecular typing (Cafiso et al., 2014) as well as plasmids and intrinsic resistance trait content; whereas acquired AMR genes and SNP profiling highlighted the acquisition of genetic traits related to the onset of additional antimicrobial resistances. In addition, virulome analysis showed a high virulence gene toxin content of the ST8 DAPR/S 3A/B strain pair with respect to the ST398 DAPR/S 1A/C.
Comparative genomics did not reveal newly acquired mobile genetic elements in both DAPR MRSA and their DAPS parents putatively related to DAP-R onset. On the contrary, in depth, SNPomics showed HI and MI DAP-MOA-related nsSNPs in shared functional coding-gene clusters, despite the diversity of the ST398 and ST8 genomic backgrounds with respect to the S. aureus NCTC8325 RefGen (proved to be the genome most closely related to both strain pairs in study).
Comparing DAPR versus their DAPS parents on S. aureus NCTC8325 RefGen, functional protein coding-gene clusters were found to be a hot spot of nsSNPs in a pairwise shared way. However, within these clusters, the detected putative HI and MI DAP-MOA-related nsSNPs in the target genes varied between the two strain pairs, indicating a strain-dependent behavior, defining two different DAPR genomic backgrounds.
Our data, supported by previous findings on laboratory induced DAPR MRSA of specific clones HG003, USA300-TCH1516, MSSA476, MW2, and MRSA252 mutants (Coe et al., 2019), first described and defined in two clinical DAPR/S MRSA of genomic backgrounds (ST8 and ST398) strain-dependent shared and functional protein clusters as hot spots of genomic variations. In particular, membrane protein (transmembrane proteins, lipoproteins, and transporters), transcriptional regulator (Sigma-70, RpoB, RpoC, RsbU, GraX, SarR, SarU, SarX, ArlS, WalK, AgrC/A, MsrR, Msa, KdpD, SAOUHSC_02390, and SAOUHSC_00673), and cell-envelope modification (Sle1, UgtP, DltB, FmtA, LspA, Cls1, MprF, and SAOUHSC_01063) protein coding-gene clusters emerged as accumulation sites of mutational events related to the DAP-R onset, randomly occurring in the same gene. Therefore, the CM structural and the functional proteins (cell envelope, cell division, and stress response system proteins) together with several transcriptional regulators represented the cellular targets undergoing genomic gain-in changes under DAP pressure.
The late growth phase transcriptional analysis, mimicking the in vivo infection-site growth, outline for the first time the long-term transcriptomic fingerprinting and networks involved in the adaptation to become DAPR MRSA. DAVID analysis highlighted common enriched pairwise over-expressed/under-expressed KEGG pathways in transport and TCRSs, SaeS/R, and WalKR, even though not in both strain pairs for WalK. Some differences in the enriched KEGG pathways can be attributable to strain-dependent features and different genomic background. The SaeR/S regulatory system (SAOUHSC_00714/SAOUHSC_00715) controls the production of exoproteins involved in adhesion and invasion of host cells, that is, hemolysins (hla and hlb), coa, Dnase, and spa- and CW-associated proteins (emp, eap, and fnbA), whereas WalK was previously related to glycopeptide resistance (Kolter et al., 1993; Shoji et al., 2011).
Transcriptomes revealed long-term imprints closely or indirectly related to the DAP-MOA and several accessory traits corroborating the multiple transcriptomic adaptations acquired with DAP-R onset, following DAP exposure and maintained by DAPR MRSA.
In particular, we considered as closely related features the CW and CM structure and organization traits as well as the primary metabolism ones owing to their link to the DAP-MOA and indirectly related features whether involved in biological pathways not directly related to DAP-MOA as the survival mechanisms of bacteria, that is, oxidative stress response, reactive oxygen species (ROS) detoxification, and ABC transporters.
The cell-envelope organization and structure are a key trait closely related to the DAP-MOA showing a multilevel feedback. Changes in the cell-envelope organization and structure appear, in fact, related to the modulation of four different pathways—that is, the peptidoglycan biosynthesis, cytolysis, cell division, and CM structure—considered related to DAP-MOA targets.
In details, adaptations in CW and CM organization and changes in the profile and content of the membrane proteins—affecting the availability of these DAP targets (Pogliano et al., 2012)—can be speculated via a modulation of the YycFG expression by yycH and membrane-protein gene under-expression. In Bacillus subtilis, yycHI is accessory system to YycFG repressing the YycG histidine kinase function (Fukushima et al., 2008). Changes in genes modulating cell-envelope stress and maintenance, including YycFG, had been previously associated with the development of DAPR in clinical and laboratory-derived DAPR S. aureus mutants, as well as with resistance to vancomycin (Utaida et al., 2003; Gaupp et al., 2015). The YycFG was reported as a key regulator of different processes affecting CW metabolism, CM lipid homeostasis and biofilm formation, changes in membrane fluidity, and CW cross-linking compensating osmotic pressure stresses (Dubrac and Msadek, 2004; Fukushima et al., 2011; Delauné et al., 2012).
Alterations in peptidoglycan synthesis were supposed via murF and SAOUHSC_00751 (functionally related to MurB) under-expression, indicating a putative decreased addition of the C-terminal D-alanyl–D-alanine dipeptide to the CW precursor muropeptide, a critical control point of peptidoglycan synthesis in S. aureus. The biosynthesis and attachment of the D-alanyl–D-alanine dipeptide are catalyzed by the different proteins encoded by ddlA and murF, with the MurF attaching the dipeptide to the UDP-N-acetylmuramic acid (MurNAc)-tripeptide, thus completing the biosynthesis of the peptidoglycan block, the UDP-linked MurNAc-pentapeptide. The D-alanyl–D-alanine C-terminal residues are essential for reactions taking place at the CW (Sobral et al., 2006). Even though DAP is structurally related to amphomycin, and similar lipopeptides as well as peptidoglycan biosynthesis inhibitors, no experimental studies share evidence on a similar MOA (Taylor and Palmer, 2016; Gómez Casanova et al., 2017). In our opinion, this signature could represent a new putative intrinsic daptomycin/glycopeptide cross-resistance mechanism because the D-alanyl–D-alanine dipeptide also represents glycopeptide targets. This transcriptomic trait experimentally supports the association of DAP-R and glycopeptide reduced susceptibility in MRSA, as well as, for the first time, providing experimental evidence on a relationship between DAPR onset and dysfunctionality in peptidoglycan biosynthesis, likely linked to a concomitance DAP inhibitor activity as reported for similar lipopeptide antimicrobials.
In addition, PBP incorporation into the CM and the quality control of cytoplasmic and integral membrane proteins could be affected through damaged-protein degradation and the protein-folding by ftsH under-expression (Lithgow et al., 2004). Autolysis is a crucial signature of the DAPR transcriptomic pathway and is closely related to the DAP-MOA, as demonstrated by cidB over-expression and lytS under-expression. Our data showed that DAP-R could lead to an increased activity of extracellular murein hydrolases via cidB over-expression. CidAB and LytS/R represent the key of the bacterial lysis regulatory pathway (Groicher et al., 2000; Patton et al., 2006). Furthermore, the LytSR system has been hypothesized to function as a staphylococcal “voltmeter,” rapidly sensing Δψ changes and leading to adaptations for resistance to host defense cationic antimicrobial peptides (HDPs) (Yang et al., 2013). Because the bactericidal mechanisms of action of the HDPs, as well as DAP, involve disruption of Δψ associated with the CM (either primarily or secondarily), DAP—owing to its similarity to HDPs—could perturb the staphylococcal CM and alter transmembrane potential (ΔΨ) impacting on programmed cell death and autolysis. Under stress conditions, for example, the presence of antimicrobials, CidAB can collapse the proton motive force and allow access of autolysins to their substrate resulting in cellular lysis modulating programmed cell death (Patton et al., 2006).
Different metabolic adaptations are directly linked to the DAP-MOA characterize DAPR MRSA. A fermentative metabolism appears to be the main metabolic pathway as expected in late growth phase, as demonstrated by the L-lactate dehydrogenase ldh2 over-expression associated with the increase in the lactic acid concentration reflecting the L-lactate dehydrogenase activity and the concomitant pyruvate kinase pyk under-expression. Shifting toward a fermentative metabolism leads to ATP deficiency and an alteration of membrane potential (ΔΨ) as a consequence of a decrease in the loss of oxidative phosphorylation, in agreement with previous findings showing that the succinate dehydrogenase levels, TCA enzyme, were lower in a DAPR strain with respect to DAPS strains (Fischer et al., 2011). menH under-expression (also functionally related to MurF) involved in menaquinone and phylloquinone epoxide biosynthesis was observed in DAPR strains. Menaquinone is a component of the staphylococcal membrane required for correct electron transport chain (ETC) functionality; thus, the menH under-expression may also be related to DAP-R acquisition.
Carbohydrate and inorganic phosphate transmembrane transport represents another key point indirectly related to the DAP-MOA in DAPR MRSA. In detail, our data showed a predominant under-expression of various transporter coding genes, including three genes responsible for carbohydrate transport. On the contrary, over-expression of the inorganic phosphate transmembrane transporter was observed. In particular, the over-expression of the inorganic phosphate transmembrane transporters could be implicated for its STRING predicted functionality related to the walR, arlR, and ssrA master regulators involved in autolysis, biofilm formation, CW metabolism, adhesion, multidrug resistance, virulence, and the global regulation of staphylococcal virulence factors in response to environmental oxygen levels as well as repression of agr, spa, and tst transcription.
Accessory DAP-R fingerprints were also discovered in DAPR MRSA. DAP-R alters the metabolism of amino-acid, lipid, cofactors, pyrimidine, and purine as demonstrated by the under-expression of various metabolic genes. Alterations in the pyrimidine and purine metabolic pathways have been indeed previously reported in DAPR S. aureus (Cui et al., 2010). Looking at the nucleic acid metabolism, different genes required for DNA replication, transcription, and translation appeared differentially expressed comparing DAPR versus DAPS strains. The over-expression of SAOUHSC_01334—annotated as sosA—could indicate an SOS-response block. sosA appears to be LexA regulated, exhibiting an over-expression similar to that of the SOS genes (Cirz et al., 2007; Mesak et al., 2008). An SOS-response block related to hypermutator behavior is considered a key mechanism of mutational antibiotic resistance. As regards the oxidative stress response, SAOUHSC_00304 mono-oxygenase over-expression could reflect the generation of ROS (Gaupp et al., 2012). In addition, the globin domain containing protein with nitric oxide dioxygenase activity (SAOUHSC_00204) under-expression—catalyzing nitric oxide (NO) to nitrate (NO3–) conversion—could cause an accumulation of reactive oxidant species such as NO. All these data could indicate the involvement of oxidative stress response enzymes in DAP-R responsible for cell survival because of an increased activity of ROS detoxification. Other different heat and cold shock stress response genes were found differentially expressed in DAPR S. aureus. In particular, dnaJ under-expression could be related to a decreased production of the chaperone protein DnaJ, involved in the response to hyper-osmotic and heat shock stress, preventing or restoring aggregation of denatured proteins. The use of different carbon sources in stress and nutrient limitation conditions were found, which could facilitate the survival of DAPR strains in starvation and stress conditions via PHB-depolymerase family protein over-expression. As regards cell adhesion, our data showed a putative decreased adhesion ability of DAPR versus DAPS S. aureus via sdrD under-expression, although this was in contrast with previous findings (Song et al., 2013). SdrD, cell surface-associated calcium-binding protein, interacts with the extracellular matrix of higher eukaryotes.
In addition, different PCAs supported the correlation of the five distinctive functional DEG clusters with a daptomycin resistance or susceptible phenotype.
New consideration should be given for mprF and dltA expression. In the late post-exponential growth phase RNA-seq data, mprF expression trend showed increased transcripts in 1C versus 1A and decreased in 3B versus 3A, even though not statistically significant. Even dltA expression did not display statistically significant differential expression in any strain pairs; however, the expression trend displayed a decreased amount of transcripts in DAPR strains in both strain pairs (data shown in files submitted to GEO). Based on these data, we hypothesize a growth-phase/strain-dependent mprF expression in the late growth phase transcriptomes correlating with the increased lysyl-phosphatidylglycerol synthesis in 1C versus 1A strain pair and a decrease in 3B versus 3A strain pair described in our previously published data on stationary cultures (Mishra et al., 2014). On the contrary, a growth phase-dependent dltA expression was observed in the late post-exponential growth phase in both strain pairs in agreement with other findings (Yang et al., 2009). Ultimately, although the bioinformatic cutoffs used to analyze these transcriptomic data are the standard ones, they could be too stringent both to simultaneously identify all RNAomic features characterizing DAPR S. aureus and to fully reflect the complexity of this biological system.
Our data define the complex genomic and long-term transcriptomic fingerprinting and adaptations of DAPR MRSA, providing new insights into their distinctive traits focusing on targets related to DAP-MOA. Briefly, we can summarize that DAPR MRSA acquired diverse genomic and transcriptomic changes to cope with and preserve bacterial cell from DAP action. CM structural/functional proteins and transcriptional regulators emerged as the cell targets related to genomic changes gained under DAP pressure. Furthermore, DAP-MOA-related transcriptomic adaptations were found in CW and CM organization, that is, peptidoglycan biosynthesis, cell division, and CM structure, as well as in the autolytic system, in primary metabolism via a shift toward fermentation, and in CM-potential perturbation. Finally, accessory traits can also impact on DAPR such as multilevel amplified stress responses mainly including oxidative stress response determining cell survival due to an increased ROS detoxification and the ABC transporters.
Data Availability Statement
The datasets generated for this study can be found in this article/Supplementary Material. WGS raw reads were deposited at Sequence Read Archive (SRA) under study accession no. SRP166981 (BioProject: PRJNA498510). RNA-seq raw reads were deposited at the Gene Expression Omnibus database (GEO) under study accession no. GSE121797 (BioProject: PRJNA498510).
VC and SSte conceived and designed the study. VC, SStr, FL, ID, and AZ performed the transcriptomics, real-time qPCR, and bioinformatics. GP contributed to the bioinformatics analysis. All authors analyzed the data and contributed to the manuscript.
This study was supported by a research grant PRIN 2017SFBFER from MIUR Italy.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We wish to thank the Scientific Bureau of the University of Catania (Italy) for language support services.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.01893/full#supplementary-material
- ^ http://genomicepidemiology.org/
- ^ http://github.com/tseemann/nullarbor
- ^ http://david.abcc.ncifcrf.gov/
- ^ http://snpeff.sourceforge.net
Allen, N. E., Hobbs, J. N., and Alborn, W. E. Jr. (1987). Inhibition of peptidoglycan biosynthesis in gram-positive bacteria by LY146032. Antimicrob. Agents Chemother. 31, 1093–1099. doi: 10.1128/AAC.31.7.1093
Bartels, D. M., Petersen, A., Worning, P., Boye Nielsen, J., Larner-Svensson, H., and Krogh Johansen, H. (2014). Comparing whole-genome sequencing with Sanger sequencing for spa typing of methicillin-resistant Staphylococcus aureus. J. Clin. Microbiol. 52, 4305–4308. doi: 10.1128/JCM.01979-1914
Bayer, A. S., Mishra, N. N., Chen, L., Kreiswirth, B. N., Rubio, A., and Yang, S. J. (2015). Frequency and distribution of single-nucleotide polymorphisms within mprF in methicillin-resistant Staphylococcus aureus clinical isolates and their role in cross-resistance to daptomycin and host defense antimicrobial peptides. Antimicrob. Agents Chemother. 59, 4930–4937. doi: 10.1128/AAC.00970-15
Bayer, A. S., Schneider, T., and Sahl, H. G. (2013). Mechanisms of daptomycin resistance in Staphylococcus aureus: role of the cell membrane and cell wall. Ann. N.Y. Acad. Sci. 1277, 139–158. doi: 10.1111/j.1749-6632.2012.06819.x
Boudjemaa, R., Cabriel, C., Dubois-Brissonnet, F., Bourg, N., Dupuis, G., Gruss, A., et al. (2018). Impact of bacterial membrane fatty acid composition on the failure of daptomycin to kill Staphylococcus aureus. Antimicrob. Agents Chemother. 62:e00023-18. doi: 10.1128/AAC.00023-18
Britt, N. S., Patel, N., Shireman, T. I., El Atrouni, W. I., Horvat, R. T., and Steed, M. E. (2017). Relationship between vancomycin tolerance and clinical outcomes in Staphylococcus aureus bacteraemia. J. Antimicrob. Chemother. 72, 535–542. doi: 10.1093/jac/dkw453
Bæk, K. T., Thøgersen, L., Mogenssen, R. G., Mellergaard, M., Thomsen, L. E., Petersen, A., et al. (2015). Stepwise decrease in daptomycin susceptibility in clinical Staphylococcus aureus isolates associated with an initial mutation in rpoB and a compensatory inactivation of the clpX gene. Antimicrob. Agents Chemother. 59, 6983–6991. doi: 10.1128/AAC.01303-15
Cafiso, V., Bertuccio, T., Purrello, S., Campanile, F., Mammina, C., Sartor, A., et al. (2014). dltA over-expression: a strain-independent keystone of daptomycin resistance in methicillin-resistant Staphylococcus aureus. Int. J. Antimicrob. Agents 43, 26–31. doi: 10.1016/j.ijantimicag.2013.10.001
Cafiso, V., Bertuccio, T., Spina, D., Purrello, S., Campanile, F., Di Pietro, C., et al. (2012). Modulating activity of vancomycin and daptomycin on the expression of autolysis cell-wall turnover and membrane charge genes in hVISA and VISA strains. PLoS One 7:e29573. doi: 10.1371/journal.pone.0029573
Cafiso, V., Stracquadanio, S., Lo Verde, F., Gabriele, G., Mezzatesta, M. L., Caio, C., et al. (2019). Colistin resistant A. baumannii: genomic and transcriptomic traits acquired under colistin therapy. Front. Microbiol. 7:3195. doi: 10.3389/fmicb.2018.03195
Carattoli, A., Zankari, E., Garcia-Fernandez, A., Voldby Larsen, M., Lund, O., and Villa, L. (2014). In silico detection and typing of plasmids using plasmidfinder and plasmid multilocus sequence typing. Antimicrob. Agents Chemother. 58, 3895–3903. doi: 10.1128/AAC.02412-14
Cingolani, P., Platts, A., Wang le, L., Coon, M., Nguyen, T., Wang, L., et al. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 6, 80–92. doi: 10.4161/fly.19695
Cirz, R. T., Jones, M. B., Gingles, N. A., Minogue, T. D., Jarrahi, B., Peterson, S. N., et al. (2007). Complete and SOS-mediated response of Staphylococcus aureus to the antibiotic ciprofloxacin. J. Bacteriol. 189, 531–539. doi: 10.1128/JB.01464-06
Coe, K. A., Lee, W., Stone, M. C., Komazin-Meredith, G., Meredith, T. C., Grad, Y. H., et al. (2019). Multi-strain Tn-Seq reveals common daptomycin resistance determinants in Staphylococcus aureus. PLoS Pathog. 15:e1007862. doi: 10.1371/journal.ppat.1007862
Cui, L., Isii, T., Fukuda, M., Ochiai, T., Neoh, H. M., Camargo, I. L., et al. (2010). An RpoB mutation confers dual heteroresistance to daptomycin and vancomycin in Staphylococcus aureus. Antimicrob. Agents Chemother. 54, 5222–5233. doi: 10.1128/AAC.00437-10
Delauné, A., Dubrac, S., Blanchet, C., Poupel, O., Mäder, U., Hiron, A., et al. (2012). The WalKR system controls major staphylococcal virulence genes and is involved in triggering the host inflammatory response. Infect. Immun. 80, 3438–3453. doi: 10.1128/IAI.00195-12
Dubrac, S., Boneca, I. G., Poupel, O., and Msadek, T. (2007). New insights into the WalK/WalR (YycG/YycF) essential signal transduction pathway reveal a major role in controlling cell wall metabolism and biofilm formation in Staphylococcus aureus. J. Bacteriol. 189, 8257–8269. doi: 10.1128/JB.00645-07
Dubrac, S., and Msadek, T. (2004). Identification of genes controlled by the essential YycG/YycF two-component system of Staphylococcus aureus. J. Bacteriol. 186, 1175–1181. doi: 10.1128/JB.186.4.1175-1181.2004
Ernst, C. M., Slavetinsky, C. J., Kuhn, S., Hauser, J. N., Nega, M., Mishra, N. N., et al. (2018). Gain-of-function mutations in the phospholipid flippase MprF confer specific daptomycin resistance. mBio 9:e01659-18. doi: 10.1128/mBio.01659-18
Fischer, A., Yang, S. J., Bayer, A. S., Vaezzadeh, A. R., Herzig, S., Stenz, L., et al. (2011). Daptomycin resistance mechanisms in clinically derived Staphylococcus aureus strains assessed by a combined transcriptomics and proteomics approach. J. Antimicrob. Chemother. 66, 1696–1711. doi: 10.1093/jac/dkr195
Friedman, L., Alder, J. D., and Silverman, J. A. (2006). Genetic changes that correlate with reduced susceptibility to daptomycin in Staphylococcus aureus. Antimicrob. Agents Chemother. 50, 2137–2145. doi: 10.1128/AAC.00039-06
Fukushima, T., Furihata, I., Emmins, R., Daniel, R. A., Hoch, J. A., and Szurmant, H. (2011). A role for the essential YycG sensor histidine kinase in sensing cell division. Mol. Microbiol. 79, 503–522. doi: 10.1111/j.1365-2958.2010.07464.x
Fukushima, T., Szurmant, H., Kim, E. J., Perego, M., and Hoch, J. A. (2008). A sensor histidine kinase co-ordinates cell wall architecture with cell division in Bacillus subtilis. Mol. Microbiol. 69, 621–632. doi: 10.1111/j.1365-2958.2008.06308.x
Gaupp, R., Lei, S., Reed, J. M., Peisker, H., Boyle-Vavra, S., Bayer, A. S., et al. (2015). Staphylococcus aureus metabolic adaptations during the transition from a daptomycin susceptibility phenotype to a daptomycin non susceptibility phenotype. Antimicrob. Agents Chemother. 59, 4226–4238. doi: 10.1128/AAC.00160-15
Groicher, K. H., Firek, B. A., Fujimoto, D. F., and Bayles, K. W. (2000). The Staphylococcus aureus lrgAB operon modulates murein hydrolase activity and penicillin tolerance. J. Bacteriol. 182, 1794–1801. doi: 10.1128/JB.182.7.1794-1801.2000
Howden, B. P., McEvoy, C. R., Allen, D. L., Chua, K., Gao, W., Harrison, P. F., et al. (2011). Evolution of multidrug resistance during Staphylococcus aureus infection involves mutation of the essential two component regulator WalKR. PLoS Pathog. 7:e1002359. doi: 10.1371/journal.ppat.1002359
Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2009a). Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 37, 1–13. doi: 10.1093/nar/gkn923
Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2009b). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211
Joensen, K. G., Scheutz, F., Lund, O., Hasman, H., Kaas, R. S., Nielsen, E. M., et al. (2014). Real-time whole-genome sequencing for routine typing, surveillance, and outbreak detection of verotoxigenic Escherichia Coli. J. Clin. Microbiol. 52, 1501–1510. doi: 10.1128/JCM.03617-3613
Jones, T., Yeaman, M. R., Sakoulas, G., Yang, S., Proctor, R. A., Sahl, H. G., et al. (2008). Failures in clinical treatment of Staphylococcus aureus infection with daptomycin are associated with alterations in surface charge, membrane phospholipid asymmetry, and drug binding. Antimicrob. Agents Chemother. 52, 269–278. doi: 10.1128/AAC.00719-07
Jordan, S., Junker, A., Helmann, J. D., and Mascher, T. (2006). Regulation of LiaRS-dependent gene expression in Bacillus subtilis: identification of inhibitor proteins, regulator binding sites, and target genes of a conserved cell envelope stress-sensing two-component system. J. Bacteriol. 188, 5153–5166. doi: 10.1128/JB.00310-06
Kilelee, E., Pokorny, A., Yeaman, M. R., and Bayer, A. S. (2010). Lysyl-phosphatidylglycerol attenuates membrane perturbation rather than surface association of the cationic antimicrobial peptide 6W-RP-1 in a model membrane system: implications for daptomycin resistance. Antimicrob. Agents Chemother. 54, 4476–4479. doi: 10.1128/AAC.00191-10
Kuroda, M., Kuroda, H., Oshima, T., Takeuchi, F., Mori, H., and Hiramatsu, K. (2003). Two-component system VraSR positively modulates the regulation of cell-wall biosynthesis pathway in Staphylococcus aureus. Mol. Microbiol. 49, 807–821. doi: 10.1046/j.1365-2958.2003.03599.x
Larsen, M. V., Cosentino, S., Rasmussen, S., Friis, C., Hasman, H., and Lykke Marvig, R. (2012). Multilocus sequence typing of total-genome-sequenced bacteria. J. Clin. Microbiol. 50, 1355–1361. doi: 10.1128/JCM.06094-6011
Mascio, C. T., Alder, J. D., and Silverman, J. A. (2007). Bactericidal action of daptomycin against stationary-phase and nondividing Staphylococcus aureus cells. Antimicrob. Agents Chemother. 51, 4255–4260. doi: 10.1128/AAC.00824-07
McCall, I. C., Shah, N., Govindan, A., Baquero, F., and Levin, B. R. (2019). Antibiotic killing of diversely generated populations of nonreplicating bacteria. Antimicrob. Agents Chemother. 63:e02360-18. doi: 10.1128/AAC.02360-18
McClure, R., Balasubramanian, D., Sun, Y., Bobrovskyy, M., Sumby, P., Genco, C. A., et al. (2013). Computational analysis of bacterial RNA-seq data. Nucleic Acids Res. 41:e140. doi: 10.1093/nar/gkt444
Mehta, S., Cuirolo, A. X., Plata, K. B., Riosa, S., Silverman, J. A., Rubio, A., et al. (2012). VraSR two-component regulatory system contributes to mprF-mediated decreased susceptibility to daptomycin in in vivo-selected clinical strains of methicillin-resistant Staphylococcus aureus. Antimicrob. Agents Chemother. 56, 92–102. doi: 10.1128/AAC.00432-10
Mesak, L. R., Miao, V., and Davies, J. (2008). Effects of subinhibitory concentrations of antibiotics on SOS and DNA repair gene expression in Staphylococcus aureus. Antimicrob. Agents Chemother. 9, 3394–3397. doi: 10.1128/AAC.01599-07
Miller, W. R., Munita, J. M., and Arias, C. A. (2016). Mechanism of action and resistance to daptomycin in Staphylococcus aureus and Enterococci. Cold Spring Harb. Perspect. Med. 6:a026997. doi: 10.1101/cshperspect.a026997
Mishra, N. N., Bayer, A. S., Weidenmaier, C., Grau, T., Wanner, S., Stefani, S., et al. (2014). Phenotypic and genotypic characterization of daptomycin-resistant methicillin-resistant Staphylococcus aureus strains: relative roles of mprF and dlt operons. PLoS One 9:e107426. doi: 10.1371/journal.pone.0107426
Mishra, N. N., McKinnell, J., Yeaman, M. R., Rubio, A., Nast, C. C., Chen, L., et al. (2011). In vitro cross-resistance to daptomycin and host defense cationic antimicrobial peptides in clinical methicillin-resistant Staphylococcus aureus isolates. Antimicrob. Agents Chemother. 55, 4012–4018. doi: 10.1128/AAC.00223-11
Muthaiyan, A., Silverman, J. A., Jayaswal, R. K., and Wilkinson, B. J. (2008). Transcriptional profiling reveals that daptomycin induces the Staphylococcus aureus cell wall stress stimulon and genes responsive to membrane depolarization. Antimicrob. Agents Chemother. 52, 980–990. doi: 10.1128/AAC.01121-07
Mwangi, M. M., Wu, S. W., Zhou, Y., Sieradzki, K., de Lencastre, H., Richardson, P., et al. (2007). Tracking the in vivo evolution of multidrug resistance in Staphylococcus aureus by whole-genome sequencing. Proc. Natl. Acad. Sci. U.S.A. 104, 9451–9456. doi: 10.1073/pnas.0609839104
Patton, T. G., Yang, S. J., and Bayles, K. W. (2006). The role of proton motive force in expression of the Staphylococcus aureus cid and lrg operons. Mol. Microbiol. 59, 1395–1404. doi: 10.1111/j.1365-2958.2006.05034.x
Pereira, M. A., Imada, E. L., and Muniz Guedes, R. L. (2017). RNA-seq: applications and best practices,” in Applications of RNA-Seq and Omics Strategies: From Microorganisms to Human Health, eds F. Marchi, P. Cirillo, and E. C. Mateo (Norderstedt: BoD). doi: 10.5772/intechopen.69250
Pfaffl, M. W., Horgan, G. W., and Dempfle, L. (2002). Relative expression software tool (REST) for group wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 30:e36. doi: 10.1093/nar/30.9.e36
Pogliano, J., Pogliano, N., and Silverman, J. A. (2012). Daptomycin-mediated reorganization of membrane architecture causes mislocalization of essential cell division proteins. J. Bacteriol. 194, 4494–4504. doi: 10.1128/JB.00011-12
Reed, J. M., Gardner, S. G., Mishra, N. N., Bayer, A. S., and Somerville, G. A. (2019). Metabolic interventions for the prevention and treatment of daptomycin non-susceptibility in Staphylococcus aureus. J. Antimicrob. Chemother. 74, 2274–2283. doi: 10.1093/jac/dkz194
Rubio, A., Conrad, M., Haselbeck, R. J., Kedar, G. C., Brown-Driver, V., Finn, J., et al. (2011). Regulation of mprF by antisense RNA restores daptomycin susceptibility to daptomycin-resistant isolates of Staphylococcus aureus. Antimicrob. Agents Chemother. 55, 364–367. doi: 10.1128/AAC.00429-10
Sabat, A. J., Tinelli, M., Grundmann, H., Akkerboom, V., Monaco, M., Del Grosso, M., et al. (2018). Daptomycin resistant Staphylococcus aureus clinical strain with novel non-synonymous mutations in the mprF and vraS genes: a new insight into daptomycin resistance. Front. Microbiol. 9:2705. doi: 10.3389/fmicb.2018.02705
Shoji, M., Cui, L., Iizuka, R., Komoto, A., Neoh, H. M., Watanabe, Y., et al. (2011). walK and clpP mutations confer reduced vancomycin susceptibility in Staphylococcus aureus. Antimicrob. Agents Chemother. 55, 3870–3881. doi: 10.1128/AAC.01563-10
Sobral, R. G., Ludovice, A. M., de Lencastre, H., and Tomasz, A. (2006). Role of murF in cell wall biosynthesis: isolation and characterization of a murF conditional mutant of Staphylococcus aureus. J. Bacteriol. 188, 2543–2553. doi: 10.1128/JB.188.7.2543-2553.2006
Song, Y., Rubio, A., Jayaswal, R. K., Silverman, J. A., and Wilkinson, B. J. (2013). Additional routes to Staphylococcus aureus daptomycin resistance as revealed by comparative genome sequencing, transcriptional profiling, and phenotypic studies. PLoS One 8:e58469. doi: 10.1371/journal.pone.0058469
Stefani, S., Campanile, F., Santagati, M., Mezzatesta, M. L., Cafiso, V., and Pacini, G. (2015). Insights and clinical perspectives of daptomycin resistance in Staphylococcus aureus: a review of the available evidence. Int. J. Antimicrob. Agents 46, 278–289. doi: 10.1016/j.ijantimicag.2015.05.008
Taglialegna, A., Varela, M. C., Rosato, R. R., and Rosato, A. E. (2019). VraSR and virulence trait modulation during daptomycin resistance in methicillin-resistant Staphylococcus aureus infection. mSphere 4:e00557-18 doi: 10.1128/mSphere.00557-18
Utaida, S., Dunman, P. M., Macapagal, D., Murphy, E., Projan, S. J., Singh, V. K., et al. (2003). Genome-wide transcriptional profiling of the response of Staphylococcus aureus to cell-wall-active antibiotics reveals a cell-wall-stress stimulon. Microbiology 149, 2719–2732. doi: 10.1099/mic.0.26426-0
Yang, S. J., Kreiswirth, B. N., Sakoulas, G., Yeaman, M. R., Xiong, Y. Q., Sawa, A., et al. (2009). Enhanced expression of dltABCD is associated with the development of daptomycin non susceptibility in a clinical endocarditis isolate of Staphylococcus aureus. J. Infect. Dis. 200, 1916–1920. doi: 10.1086/648473
Yang, S. J., Xiong, Y. Q., Yeaman, M. R., Bayles, K. W., Abdelhady, W., and Bayer, A. S. (2013). Role of the LytSR two-component regulatory system in adaptation to cationic antimicrobial peptides in Staphylococcus aureus. Antimicrob. Agents Chemother. 57, 3875–3882. doi: 10.1128/AAC.00412-13
Yoon, Y. K., Park, D. W., Sohn, J. W., Kim, H. Y., Kim, Y. S., Lee, C. S., et al. (2016). Effects of inappropriate empirical antibiotic therapy on mortality in patients with healthcare-associated methicillin-resistant Staphylococcus aureus bacteremia: a propensity-matched analysis. BMC Infect. Dis. 15:331. doi: 10.1186/s12879-016-1650-8
Zankari, E., Hasman, H., Cosentino, S., Vestergaard, M., Rasmussen, S., Lund, O., et al. (2012). Identification of acquired antimicrobial resistance genes. J. Antimicrob. Chemother. 67, 2640–2644. doi: 10.1093/jac/dks261
Keywords: daptomycin-resistant methicillin-resistant Staphylococcus aureus, genomics, SNPome, virulome, resistome, RNA-seq, transcriptomics
Citation: Cafiso V, Stracquadanio S, Lo Verde F, De Guidi I, Zega A, Pigola G and Stefani S (2020) Genomic and Long-Term Transcriptomic Imprints Related to the Daptomycin Mechanism of Action Occurring in Daptomycin- and Methicillin-Resistant Staphylococcus aureus Under Daptomycin Exposure. Front. Microbiol. 11:1893. doi: 10.3389/fmicb.2020.01893
Received: 09 December 2019; Accepted: 20 July 2020;
Published: 14 August 2020.
Edited by:Flavia Marinelli, University of Insubria, Italy
Reviewed by:Pedro Matos Pereira, New University of Lisbon, Portugal
Iñigo Lasa, Universidad Pública de Navarra, Spain
Copyright © 2020 Cafiso, Stracquadanio, Lo Verde, De Guidi, Zega, Pigola and Stefani. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Viviana Cafiso, firstname.lastname@example.org