Comparative transcriptomic analysis of Staphylococcus epidermidis associated with periprosthetic joint infection under in vivo and in vitro conditions.

Nineteen S. epidermidis isolates and sonicate fluids from 19 individuals diagnosed with PJI were studied. Clinical profiles of the subjects are shown in Table 1 . Subject ages ranged from 47 to 83 years (average of 65 years) and included 13 (68%) males. Samples included those from 13 (68%) knee and 6 (32%) hip resection arthroplasties. Implant ages ranged from 0.3 to 12.7 years (average of 2.6 years). Eleven (58%) of 19 subjects had had symptoms for ≥1 year, with eight (42%) presenting <1 year after symptom onset. Pre-operative antibiotics had been administered in 9 (47%) of 19 cases. Whole genome sequencing, followed by de novo genome assembly was performed to characterize the PJI-associated S. epidermidis isolates. Multilocus sequence typing via PubMLST found that the 19 PJI-associated isolates spanned a representative sampling of sequence types (ST), with ST-5, ST-130, and ST-83 accounting for 5 (26%), 4 (21%), and 2 (11%) of the 19 isolates in the cohort, respectively ( Table 1 ). RNA-sequencing was performed to assess differential gene expression between in vivo S. epidermidis, tested directly from sonicate fluid, and that of the matching isolates grown under in vitro laboratory conditions. In all, 1919 transcripts were found across all samples, 145 of which were DEGs as defined by adjusted p-values of ≤0.05 and log 2 fold-change values ≥2 or ≤−2 ( Fig. 1 ). Of the 145 DEGs, 42 were upregulated and 103 downregulated in in vivo compared to in vitro samples ( Fig. 1A ). Only 570 (30%) of the 1919 transcripts were annotated, leaving 1349 (70%) unannotated transcripts encoding for “hypothetical” proteins ( Fig. 1B ). Seventy-seven (53%) of the 145 DEGs were annotated. Of those, 13 (17%) were upregulated and 64 (83%) downregulated. Sixty-eight (47%) of 145 DEGs were unannotated. Of those, 29 (43%) were upregulated and 39 (57%) downregulated ( Fig. 1C ). A full breakdown of upregulated and downregulated DEGs is shown in Table S1 and Table S2 of the supplemental material , respectively. Of note, genes related to metal ion sequestration were elevated in a subset of in vivo compared to in vitro samples. Four genes within the cnt operon ( cntA, cntK, cntL, and cntM ), encoding for production of staphylopine-related proteins, were upregulated in a subset of in vivo compared to in vitro samples. cntA, cntK, cntL, and cntM were upregulated in 7, 2, 6, and 6 of the 19 in vivo isolates, respectively ( Fig. 2 ). cntA had average transcript per million (TPM) counts of 12.8 and 4.2 in vivo and in vitro, respectively ( Fig. 2A ). cntK had average TPM counts of 22.0 and 7.3 in vivo and in vitro, respectively ( Fig. 2B ). cntL had average TPM counts of 46.9 and 7.0 in vivo and in vitro, respectively ( Fig. 2C ). cntM had average TPM counts of 11.9 and 3.7 in vivo and in vitro, respectively ( Fig. 2D ). Further analysis was conducted to determine whether individual S. epidermidis strains could be differentially clustered based on whether they were harvested from in vivo or in vitro sources ( Fig. 3 ). Heatmap analysis using k-means clustering resulted in separation of samples based on harvest source, but not by strain ( Fig. 3A ). PCA was used for multivariate reduction analysis to compare in vivo and in vitro transcriptomic expression profiles, where samples with more similar transcriptomic expression profiles cluster more closely on a 2-dimensional plot. Both dimensions (Dim1 and Dim2) correspond to a percentage of total variation within the data set ( Jolliffe, 2002 ). PCA using transcriptomic expression profiling of all transcripts separated in vivo and in vitro samples, with Dim1 and Dim2 accounting for 50.4% and 21.9% of the total dataset variation, respectively ( Fig. 3B ). PCA using only the DEGs was also able to separate in vivo and in vitro samples, with Dim1 and Dim2 accounting for 47.9% and 18.7% of the total dataset variation, respectively ( Fig. 3C ). in vivo samples from subjects with and without pre-operative antibiotic treatment were compared, to assess whether treatment may have affected transcriptomic expression ( Fig. 4 ). Neither heatmap k-means clustering ( Fig. 4A ), nor PCA ( Fig. 4B ) differentiated samples based on pre-operative antibiotic treatment. Genomics may influence the propensity for a bacterial strain to exhibit a pathogenic or commensal phenotype. As such, the genomes from the 19 PJI-associated S. epidermidis isolates were compared to 21 publicly available genomes of commensal S. epidermidis isolates from NCBI GenBank ( Fig. 5 ). Overall, there were no significant differences between the genomic make-up of PJI and commensal S. epidermidis isolates. Functional comparisons were conducted using Cluster of Orthologous Groups of proteins (COG) 20 and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses ( Galperin et al., 2021 ; Kanehisa et al., 2023 ; Kanehisa, 2019 ; Kanehisa and Goto, 2000 ). While there were slight differences in enriched functional pathways with some individual isolates, there were no significant differences when comparing PJI to commensal S. epidermidis overall.

Sonicate fluids from 19 individuals undergoing total knee (n=13) or hip (n=6) implant removal due to S. epidermidis -associated PJI were studied. Samples were collected between 2008 and 2017 under Mayo Clinic Institutional Review Board 09–000808. Sonicate fluid samples were derived from explanted devices processed following an established laboratory protocol involving vortexing and sonication steps ( Piper et al., 2009 ; Trampuz et al., 2007 ). After processing, sonicate fluid samples were frozen at −80°C, without an RNA stabilizer, for later analysis. Infection status was determined by established diagnostic criteria ( Parvizi et al., 2018 ). Bacterial cultures of each sonicate fluid sample had yielded ≥100 colony forming units S. epidermidis /10 mL fluid, without other microorganisms, as identified by Mayo Clinic Clinical Microbiology Laboratory standard protocols. Further experimentation was performed on two different growth sources, as was previously described for S. aureus and Streptococcus agalactiae samples ( Masters et al., 2021 ; Cho et al., 2021 ). In vivo samples were sonicate fluids tested directly and in vitro samples were the original S. epidermidis isolates from sonicate fluid (frozen at −80°C) regrown in a laboratory setting by inoculation of overnight isolate subculture from the frozen stocks in 50 mL sterile tryptic soy broth at 37°C with shaking at 200 rpm until an OD 600 of 0.55–0.70 was reached via spectrophotometry. Bacterial cells were centrifuged at 5000 rpm for 5 minutes and washed twice with RNase-free water. Cell pellets were snap-frozen at −80°C and stored overnight before RNA extraction. Whole genome sequencing was performed on cultured S. epidermidis isolates to provide genomic context for downstream RNA-sequencing analyses. Genomic DNA of the 19 S. epidermidis isolates was extracted using previously described methods ( Masters et al., 2021 ). Briefly, DNA was extracted using the Zymo Research Quick-DNA Fungal/Bacterial Miniprep kit (Zymo Research, Irvine, CA, US) and quantified using a Qubit 2.0 Fluorometer (Thermo Fisher Scientific). Approximately 0.5 mg of DNA was prepared for DNA libraries using a Nextera XT PE kit (Illumina, San Diego, CA, US). Whole genome sequencing was performed on an Illumina HiSeq 4000 with a 2 ×150-bp setting and 60 sample libraries multiplexed per flow cell. SKESA v2.4.029 was used to assemble de novo bacterial genomes from raw reads ( Souvorov et al., 2018 ). Prokka v1.14.5 was used to annotate the assemblies, with S. epidermidis ATCC 12228 used as the reference genome ( Seemann, 2014 ). Multilocus sequence typing (MLST) on the 19 S. epidermidis isolates was conducted using a public S. epidermidis PubMLST scheme based on sequence variation in seven housekeeping loci, including arc, aroE, gtr , mutS, pyrR, tpiA, and yqiL ( Jolley et al., 2018 ). A pangenome was constructed using annotated fragmented de novo assemblies to identify core and accessory genes ( Page et al., 2015 ). The pangenome served as a common reference for transcript quantification and differential expression analysis between in vivo and in vitro samples. Total RNA from in vivo and in vitro samples was isolated as previously described ( Masters et al., 2021 ). Briefly, cell suspensions were mechanically homogenized, and genomic DNA and bacterial ribosomal RNA (rRNA) removed using a DNase I kit (QIAGEN, Germantown, MD, US) and a Ribo-Zero rRNA removal kit (Bacteria) (Illumina, San Diego, CA, US), respectively. Total sample RNA was isolated using the miR-Neasy Serum/Plasma kit (QIAGEN). rRNA-depleted RNA was purified using the RNeasy MinElute Cleanup kit (QIAGEN) and eluted in a final volume of 30 μL RNase-free water. RNA quantity and quality were evaluated using a Qubit 2.0 Fluorometer coupled with a Qubit RNA high-sensitive assay kit (Thermo Fisher Scientific, Waltham, MA, US) and an Agilent 4200 TapeStation system (Agilent, Santa Clara, CA, US). An Ovation SoLo RNA-Seq System (NuGEN Technologies, San Carlos, CA, US) was used to prepare cDNA libraries, as previously described ( Masters et al., 2018 ). External RNA Control Consortium RNA Spike-In Mix 1 (Thermo Fisher Scientific) was added to each total RNA sample to measure variability in the library generation process ( Jiang et al., 2011 ). Human rRNA sequences were depleted using the SoLo Any-Deplete probe mix (NuGEN Technologies). Resulting cDNA libraries were sequenced on an Illumina HiSeq 4000 platform with 10 samples multiplexed per lane, producing 100-nucleotide, paired-end reads. Human non-rRNA was not depleted before library preparation allowing for later human transcriptomic studies, so raw sequencing reads from the sonicate fluids contained both human and microbial RNA. The S. epidermidis transcriptome was analyzed using a previously described taxonomic assignment pipeline ( Masters et al., 2018 ; Thoendel et al., 2016 ), with modifications. Adapters were trimmed using Atropos 1.1.1936 ( Didion et al., 2017 ) and human reads removed using BioBloom Tools v2.1.1 ( Chu et al., 2014 ). Microbial classification was conducted using the Livermore Metagenomics Analysis Toolkit 1.2.6 with k-mer identifiers and the kMLþH.noprune.4–14.2025.db database ( Ames et al., 2013 ), and the taxonomic classifier Kaiju ( Menzel et al., 2016 ). S. epidermidis transcript quantification was performed for each sample with kallisto v0.46.1 ( Bray et al., 2016 ). A pangenome of assembled isolates and commensal sequences was constructed using Roary V3.13.0, with the gene IDs of each assembly mapping to a gene cluster ID in the pangenome to compare differences in expression ( Page et al., 2015 ). Differentially expressed genes (DEGs) were determined using a combination of txImport v1.28.0 ( Soneson et al., 2015 ) to merge transcript counts across samples using pangenome gene cluster IDs from Roary ( Page et al., 2015 ), and LinDA v1.1 ( Zhou et al., 2022 ) for differential analysis, addressing the sparsity of the clinical sample data, differences in sequencing depth between environments, and within-strain correlations. DEG statistical significance was determined as transcripts with adjusted p-values ≤0.05 and log 2 fold-change values ≥2 or ≤−2 ( Zhou et al., 2022 ). Genomes from the 19 PJI-associated S. epidermidis isolates were compared to 20 publically available genomes of known commensal S. epidermidis isolates from the National Center for Biotechnology Information (NCBI) GenBank ( Sayers et al., 2022 ). A contig database was created using Anvi’o v7.1, with functional prediction performed using the anvi-run-ncbi-cogs and anvi-run-kegg-kofams programs ( Eren et al., 2015 , 2021 ; Shaiber et al., 2020 ). Functional enrichment tests were performed using the Anvi’o programs anvi-compute-functional-enrichment-in-pan and anvi-compute-metabolic-enrichment for COG20 and KEGG annotations, respectively ( Galperin et al., 2021 ; Kanehisa et al., 2023 ; Kanehisa, 2019 ; Kanehisa and Goto, 2000 ). Functions with adjusted p-values ≤0.05 were considered to be enriched. Data was analyzed, organized, and graphed in RStudio v2022.12.0+353 ( Team" RC, 2020 ) using R-packages, including “factoextra” for principal component analysis (PCA) ( Kassambara and Mundt, 2020 ), and “ComplexHeatmap” for heatmap creation and k-means clustering ( Gu and Schlesner, 2016 ). Graphpad Prism 9 v9.5.1 (San Diego, CA) was used for volcano plot and bar graph creation.

Studies to investigate the transcriptomic response of S. epidermidis under various laboratory and clinical conditions have been previously conducted ( Liu et al., 2020 ; Carvalhais et al., 2014 , 2015 ; Franca et al., 2014 ; Chen et al., 2022 ; Franca et al., 2016 ; Moran et al., 2017 ; Teichmann et al., 2022 ), though alterations in gene expression during PJI were largely unexplored prior to this study. Here, transcriptomic profiles of paired PJI-associated S. epidermidis samples from in vivo sources – that is, those in from sonicate fluid –, and in vitro sources – that is, those grown in laboratory conditions –, were compared using RNA-sequencing to determine how alterations in gene expression may allow these bacteria to thrive under pathogenic in vivo and in vitro states. Transcriptomic profiling showed differential gene expression when comparing in vivo and in vitro samples, with 42 identified upregulated and 103 identified downregulated genes in in vivo samples, out of the total 1919 genes analyzed. Moreover, k-means clustering and PCA accurately differentially separated samples based on their in vivo or in vitro growth sources, albeit not based on their strain pairing. Even though the paired strains are genetically identical, samples have more similar transcriptomic profiles when compared to those from the same in vivo or in vitro source, than when compared their counterparts from a different source. When further analyzing the DEGs between in vivo and in vitro samples, an apparent mechanistic pattern emerged. Metal sequestration-related genes, particularly those of the staphylopine-related cnt operon, were upregulated in a subset of in vivo compared to in vitro isolates. Staphylopine is a metallophore produced by a variety of human pathogenic bacteria. Its role in the sequestration and acquisition of transition metals, such as Zn 2+ , Ni 2+ , and Co 2+ , is important for metal ion homeostasis and induction of virulence factors ( Song et al., 2018 ; Ghssein and Ezzeddine, 2022 ). Utilization of staphylopine, calprotectin, lipocalin, and other sequestration systems is part of a competitive process, also known as nutritional immunity, between host and pathogen, for transition metals within the infective microenvironment ( Price and Boyd, 2020 ; Cassat and Skaar, 2012 ; Hood and Skaar, 2012 ; Núñez et al., 2018 ; Zackular et al., 2015 ; Palmer and Skaar, 2016 ). Interestingly, upregulation of inorganic ion transport and metabolism genes, including cntA and cntL, has recently been described as crucial to the multistress response, such as during infection and colonization, in S. epidermidis ( Spoto et al., 2022 ). The results found here show that cnt operon genes cntA, cntK, cntL, and cntM are upregulated in a subset of PJI isolates. Proteins associated with cntK, cntL, and cntM are needed to synthesize the staphylopine molecule, while the protein associated with cntA is important for recognition and initiation of staphylopine/metal complex import into the cell ( Song et al., 2018 ). Upregulation of transition metal-related gene sets, particularly those involving Ni 2+ and Co 2+ , in in vivo compared to in vitro samples recapitulates results previously described in S. aureus and S. agalactiae samples, using similar methods as described here ( Masters et al., 2021 ; Cho et al., 2021 ). These findings suggest a potential role for acquisition of transition metals, such as Ni 2+ and Co 2+ , in the pathogenesis of S. epidermidis, S. aureus, and S. agalactiae -associated PJIs. Upregulation of iron-acquisition genes, such as those encoding for staphyloferrin A and B or the Isd heme consumption system, has previously been described in S. aureus PJI ( Masters et al., 2021 ). Of note, expression of these genes is entirely absent in the S. epidermidis in both the in vivo and the in vitro isolates analyzed here. This discrepancy between S. aureus and S. epidermidis analyses could represent differences in species-specific metal ion response during PJI, though study limitations, including sampling and sequencing protocol differences, may also play a role. Nonetheless, additional investigation into the role of the iron-acquisition response in both S. aureus and S. epidermidis -associated PJI is necessary to make any definitive conclusions. A major limitation of RNA-sequencing of S. epidermidis is the relative lack of annotated genes available. In this analysis, over 70% of the transcripts and almost 50% of the DEGs were unannotated. Due to these unannotated transcripts, actual functional differences between S. epidermidis from in vivo and in vitro sources remains largely undescribed. Continued investigations to annotate the S. epidermidis genome will lead to a more robust understanding of the how S. epidermidis adapts to the host environment during infection. In 2012, Conlan et al. described the presence of the formate dehydrogenase gene, fdh, as a discriminatory marker between skin commensal and nosocomial S. epidermidis ( Conlan et al., 2012 ), though another study found no association of suggested commensalism markers ( fdh and ACME-arcA ), or biofilm-related genes ( icaA, aap, bhp, embp ) with S. epidermidis isolated from PJI compared to those from nares ( Mansson et al., 2021 ). Additionally, it has been proposed that there is an evolutionarily diverse population of commensal S. epidermidis present on the skin ( Zhou et al., 2020 ), some of which may be genetically predisposed to being PJI-associated pathogens when introduced to the joint microenvironment. To determine whether genomic differences may be related to pathogenic or non-pathogenic phenotypes, the study cohort of PJI-associated isolates was compared to a cohort of known commensal S. epidermidis isolates. Genomic comparisons found that there was no significant difference in the genomic make-up of known pathogenic and commensal S. epidermidis strains. Furthermore, genomic COG20 and KEGG pathway enrichment analysis found no correlative differences in infectious and non-infectious S. epidermidis phenotypes. Further, expression of commensalism makers ( fdh and ACME-arcA ) and biofilm-related genes ( icaA, aap, bhp, embp ) was not found via RNA-sequencing, though fdhD, a gene encoding a sulfur-carrier protein required for formate dehydrogenase activity ( Thome et al., 2012 ; Arnoux et al., 2015 ), was present in the genome of 13/19 (68%) the PJI and 18/21 (86%) of commensal isolates, ACME-arcA was found in 19/19 (100%) of PJI and 21/21 (100%) of commensal isolates genomes, and icaA in 6/19 (32%) of PJI and 6/21 (29%) of commensal isolate genomes. aap and embp were not found in the genome of any isolate. Together, these results suggest that functional differences in S. epidermidis isolates found in PJI and commensal niches are primarily driven by differences in the surrounding microenvironment rather than genomic differences in the microbe. For example, the presence and/or abundance of transition metal ions may play a role in altering bacterial function during PJI, as evidenced by the elevation of staphylopine-associated transcripts in vivo. Further investigation of additional environmental characteristics is warranted. This study has several limitations. First, the sample size is relatively small. Only 19 pairs of sonicate and isolate samples were investigated, 18 (95%) of which were from Caucasian subjects. Whether race plays a role in transcriptomic differences of PJI-associated S. epidermidis is unknown. Next, RNA was isolated from clinical samples multiple years after initial harvest and samples were not stored with an RNA protector. As such, long-term storage may reduce RNA quality and impact transcriptomic profiling. Also, the process of sonication itself may impact the transcriptomic profile. Additionally, S. epidermidis ATCC 12228 was used as the sole reference genome and may not include all possible S. epidermidis genes. Therefore, it is unknown whether lack of expression of genes not found in this study is biologically accurate, or a limitation of the transcriptomic methods utilized herein. Finally, robust functional and integrated analyses were limited by the relative lack of available transcript annotation, as previously described. In conclusion, this study provides insights into differences in gene expression of S. epidermidis when grown in pathogenic and in vitro settings. Further investigation of these differences may contribute to a deeper understanding of S. epidermidis as a pathogen, potentially leading to the development of novel therapeutic targets for S. epidermidis -associated infections.

Staphylococcus epidermidis is a coagulase negative staphylococcal species that makes up a major proportion of the normal human microbiota, residing on the skin and mucous membranes, such as the nares, the gastrointestinal tract, and the lower reproductive tract ( Sabaté Brescó et al., 2017 ; Leonel et al., 2019 ; Severn and Horswill, 2023 ; Otto, 2009 ). As a pathogen, S. epidermidis is a cause of a wide range of diseases including some serious life-threatening infections, such as neonatal sepsis ( Otto, 2009 ; Post et al., 2017 ; Joubert et al., 2022 ; Greenfield et al., 2021 ). S. epidermidis is also associated with medical implant-associated biofilm-related infections, where S. epidermidis has been reported in approximately 25–40% of periprosthetic joint infection (PJI) cases and is the most frequent cause of central venous catheter and catheter-related bloodstream infections ( Sabaté Brescó et al., 2017 ; Tai et al., 2022 ; Patel, 2023 ; Fisher and Patel, 2023 ; Tande and Patel, 2014 ; Del Pozo and Patel, 2009 ; Santarpia et al., 2016 ; Wu et al., 2017 ). S. epidermidis infections are often difficult to treat due to their robust biofilm growth and risk of recidivism after treatment has concluded, leading to the need for prolonged antimicrobial therapies ( Morgenstern et al., 2016 ; Romano et al., 2011 ). Previously, S. epidermidis was thought to possess a limited virulence factor repertoire compared to Staphylococcus aureus, a related and more well-studied pathogenic Staphylococcus species, accounting for it rarity in causing invasive infections in otherwise healthy individuals ( Gill et al., 2005 ). Recently, genomic sequencing has identified an increasing number of S. epidermidis virulence factors potentially important in PJI and other infections ( Severn and Horswill, 2023 ; Franca et al., 2021 ; Cheung et al., 2014 ; Hanke and Kielian, 2012 ; Perez and Patel, 2018 ; Vuong and Otto, 2002 ). The ability of S. epidermidis to cause persistent medical device-related infections is largely attributed by its ubiquitous presence on human skin and capacity to quickly adhere to and form robust biofilms on devices. Despite the high prevalence of S. epidermidis in device-related infections, including those associated with PJI, studies elucidating specific molecular mechanisms employed by this bacterium in difficult-to-treat infections remain limited and underrepresented in the literature. Investigation of S. epidermidis pathogenicity is especially lagging behind that of S. aureus. In this study, S. epidermidis from sonicate fluid (hereafter labeled “ in vivo ”), a clinical sample-type derived from PJI subjects, and paired laboratory-grown (hereafter labeled “ in vitro ”) S. epidermidis isolates were compared for gene expression profiles using a high throughput RNA-sequencing approach. Genotypes of PJI-associated isolates were also compared to a dataset of publicly available commensal S. epidermidis isolates to elucidate potential genomic contributors to pathogenic versus commensal behavior. Results from this study provide insights that enhance understanding of pathogenic mechanisms of this bacterium.