Comparative transcriptomic adaptations of Staphylococcus aureus to the wound environment in non-diabetic and diabetic mice.

Infection is a major source of complications in delayed diabetic wound healing. Increased understanding of differential bacterial responses to diabetic wounds will enable us to better understand chronic wound pathogenesis. Here we create delayed-healing wounds infected with Staphylococcus aureus in non-diabetic and diabetic mice and used RNA-seq to compare bacterial gene expression profiles 3 or 7 days after infection. Analysis at day 3 demonstrated substantial transcriptomic differences between bacteria colonising non-diabetic and diabetic wound beds. Most of these transcriptional differences resolved by day 7, suggesting normalisation of many bacterial phenotypes later in the diabetic wound healing process. Lingering differentially expressed genes at day 7 were enriched for genes related to carbohydrate metabolism, which includes genes of the lac operon, and capsular polysaccharide synthesis, which includes the cap8 locus. These data encourage further research into host-pathogen interactions in wound healing and how they influence differential outcomes in the diabetic wound environment.

INTRODUCTION

Microbial infection is a critical complication of delayed wound healing; the wound bed can become infected with nosocomial pathogens or pathogens arising elsewhere in the environment, leading to slowing of wound healing, maintaining the open wound and leading to worsened patient outcomes. One of several infectious wound‐colonising microorganisms of particular concern is Staphylococcus aureus, of which strains resistant to antibiotics are commonly found in hospital settings. S. aureus is increasingly responsible for wound infections, particularly for those of chronic diabetic wounds including diabetic foot ulcers (DFUs). , , Much research has sought to uncover host organism factors underlying the healing deficiencies associated with diabetic wounds relative to their non‐diabetic counterparts. However, less research has sought to understand the differences in bacteria colonising diabetic and non‐diabetic wounds. Here we utilise use RNA‐sequencing in order to compare transcriptional profiles of a well‐characterised clinical isolate of S. aureus commonly used in research (UAMS‐1) during infection of the wound bed of full‐thickness splinted excisional wounds in non‐diabetic and diabetic mice, in order to understand the differences in transcriptional adaptations within these wounds.

METHODS

All methods are detailed in the Supporting Information Methods section (Appendix S1).

RESULTS

Generation and mapping of wound bacterial RNA samples

In order to compare S. aureus gene expression in infected wounds between non‐diabetic and diabetic mice, we utilised an infected wound model (Figure 1A). Splinted, full‐thickness excisional wounds were performed on mouse dorsa, inoculated with S. aureus strain UAMS‐1, and covered in Tegaderm wound dressing. Wounds were allowed to heal normally and infected wound beds were harvested after 3 or 7 days, RNA was isolated, and RNA‐seq analysis was performed. Since the UAMS‐1 genome has already been assembled and annotated, we used this construction to map the reads from our RNA‐seq analysis. In order to ensure that our UAMS‐1 isolate was well‐represented by this published genome construction, we also harvested RNA from in vitro planktonic cultures of UAMS‐1 and performed RNA sequencing on these samples; when mapping these reads to the published UAMS‐1 genome, we found that the mean fraction of mapped reads across these samples was >99% (Figure 1B). For wound samples, since the RNA was harvested from infected wound beds, even though we performed a low‐speed centrifugation step to deplete host cells (see Supporting Information Methods section [Appendix S1]), we suspected that there may also be mouse RNA in these samples. Therefore, we mapped all reads from each wound sample to both the UAMS‐1 genome and the M. musculus genome (mm10 construction). Mapping of wound samples demonstrated that significantly greater fractions of reads from d3 samples mapped to the UAMS‐1 genome compared with d7 samples, while no difference in fractions of reads mapping to the UAMS‐1 genome was found between wild type (WT) and db/db (db) wound beds at either time point (Figure S1A). Consistent with this finding, significantly greater fractions of reads from d7 samples mapped to the M. musculus genome compared with d3 samples, while no difference in fractions of reads mapping to the M. musculus genome was found between WT and db wound beds at either time point (Figure S1B). The most abundant transcripts represented in the reads that mapped to the M. musculus genome included S100a8, S100a9, Il1b, Cxcl2, and Ccl3 (Figure S2), suggesting that the murine RNA contained in the samples was likely derived predominantly from immune cells. As would be expected, the fraction of reads from planktonic UAMS‐1 cultures that mapped to the M. musculus genome was very low (1.41% ± 0.08% of reads, Figure S1B). Visualisation via principle component analysis (Figure S3) and heatmap (Figure S4) demonstrated dramatic differences between gene expression of wound bed and planktonic samples. This is expected, since bacteria must adapt when they encounter the complex wound environment (differentially‐expressed genes [DEGs] between wound bed and planktonic samples in are listed in Table S1, and corresponding gene ontology analysis detailed in Table S2). Taken together, these data suggest that samples harvested from wound beds following our protocol contained a mixture of RNA from S. aureus and M. musculus, and that S. aureus reads mapped to the UAMS‐1 genome could be used to further investigate differential gene expression profiles.

FIGURE 1

RNA‐seq analysis on infected wound bed samples mapped to the UAMS‐1 genome. (A) Photographs of wound bed (top) directly after performing the wounding procedure and (bottom) immediately before harvest. (B) Mean fractions of processed reads that mapped or failed to map to the S. aureus UAMS‐1 genome. BL6 (wild type), db (db/db), Pl (planktonic). (C) Principle component analysis comparing wound bed‐derived sample expression profiles of wild type and db/db mice at day 3 and day 7. (D,E) Gene‐clustered heatmaps demonstrating Z scores computed across all wound bed samples mapped to the UAMS‐1 genome. Genes pictured are all of those significantly differentially expressed between wild type and db/db wound beds at (D) d3 or (E) d7. (F) Venn diagram depicting the intersecting and non‐intersecting sets of differentially expressed genes (DEGs) between wild type and db/db mice at d3 and d7

Differential gene expression between S. aureus colonising wild type and db/db mouse wounds

Next, we wished to examine the existence and nature of differential gene expression in S. aureus colonising the wound bed in wild type and db/db mice. Principle component analysis of wound bed samples revealed that gene expression profiles between wild type and db/db mice varied substantially more at day 3 than at day 7 (Figure 1C). Next, we used EdgeR to investigate significant (Benjamini–Hochberg FDR P < 0.05) S. aureus DEGs in day 3 and day 7 samples between wild type and db/db mice. At day 3, 819 genes were differentially expressed between samples from wild type and from db/db mice, most of which were no longer differentially expressed at day 7 (Figure 1D,F, Table S3). Gene ontology analysis of the genes differentially expressed at d3 revealed enrichment of cellular responses to DNA damage and stress, DNA replication and repair, protein folding, translation, and numerous metabolic pathways (Table S4). In contrast, at day 7, only 90 genes were differentially expressed (Figure 1E,F, Table S5), which were enriched almost exclusively for metabolic processes (Table S6). Of the 90 genes differentially expressed at day 7, ~74% (67/90) of these genes had also been dysregulated at day 3 (Figure 1E,F). Thus, these 67 genes (listed in Table S7) represented a conserved set of genes dysregulated in db/db wounds that did not resolve by day 7. Gene ontology analysis of this list of genes demonstrated enrichment of several carbohydrate‐related metabolic processes (Table S8). We then extracted the normalised count values of two families of genes, cap8 and lac, that were represented in this conserved list of 67 dysregulated genes. Graphing expression of the genes in the cap8 (Figure 2A) and lac (Figure 2B) families demonstrated decreased expression of lac genes and increased expression of cap8 genes in db/db wound beds compared with wild type wound beds at both time points. We also extracted expression values for numerous genes canonically associated with biofilm formation and virulence, but most of these genes did not vary significantly among experimental groups (data not shown). Exceptions include arlR and arlS, genes encoding the ArlRS two‐component system, the virulence‐associated transcriptional regulators mgrA and sigB, and genes encoding members of the bacterial ssp protease system (Figure 2C). Most of these genes demonstrated significant upregulation in db wounds at day 3 and, in the case of arlS, at day 7 as well. Taken together, these data demonstrate extensive normalisation of aberrant bacterial transcription profiles over the course of infected diabetic wound healing to better resemble expression of those infecting non‐diabetic wounds, but also suggest remnant, lingering aberrations in bacterial metabolism and in other potentially virulence‐associated factors as well.

FIGURE 2

Relative expression of selected genes differentially expressed between wild type (BL6) and (db/db) wound beds. Box and whisker plots of log2(CPM) are displayed for (A) lac (B) cap8 (C) arlRS, mgrA, sigB, and sspABC family genes. n = 5 samples per group. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001. P values visualised for statistical comparisons are derived from the differential gene expression analysis by edgeR, utilising trimmed mean of M‐values (TMM) normalisation and the Benjamini–Hochberg P value correction for multiple comparisons

DISCUSSION

We set out to better understand the differences in bacterial transcriptomes of S. aureus colonising non‐diabetic and diabetic wounds, and how these transcriptional differences changed over the course of healing. Much research on infected diabetic wounds has sought to characterise aberrancies and deficiencies in the host response to infection, , , , which is doubtlessly critical to furthering our understanding of the basic science of chronic wound healing and of therapeutic modalities to promote full closure in these wounds. Less explored, however, are the ways in which the diabetic wound environment affects phenotypes and gene expression of infecting bacteria. Here, we utilised a splinted excisional wound model in wild type and db/db mice infected with S. aureus in order to simulate the delayed, infected wound environment. RNA isolation and transcriptome analysis by RNA‐seq enabled us to compare gene expression profiles of bacteria infecting diabetic and non‐diabetic wounds harvested at different time points. Large transcriptomic differences were visible between S. aureus in WT and db/db mice at day 3 post‐infection (Figure 1C,D). Most of these differences had resolved by day 7 (Figure 1E,F), suggesting normalisation of the db/db wound environment, and thus the transcriptional profiles of the infecting bacteria, at this later time point. Notably, the remaining genes dysregulated at d7 were enriched in several metabolic pathways including carbohydrate metabolism (Tables S6 and S8), suggesting that remnant differences in bacterial metabolism might contribute to diabetic infected wound phenotypes as wound healing progresses, even as some other phenotypes have become more regular. In particular, S. aureus infecting diabetic wounds demonstrates decreased expression of genes encoded by the lac operon (Figure 2A). This is consistent with diabetic hyperglycaemia, as glucose is a known negative regulator of the lac operon, as well as with the critical roles of metabolism in the regulation of S. aureus virulence. , Diabetic wound samples also demonstrated consistent upregulation of the genes of the Cap8 locus (Figure 2B), which are responsible for capsular polysaccharide synthesis, a critical determinant of S. aureus virulence and propensity for immune evasion, , suggesting that this may be a factor contributing to heightened bacterial virulence and aberrant inflammatory responses in infected diabetic wounds.

It is critical to interpret our findings in the context of their limitations. As all of the mice used here were male, we are unable to speculate on potential sex‐specific effects relevant to the diabetic wound infection paradigm. Though the splinted murine wound model heals through granulation and re‐epithelialization in a manner similar to human wounds, mouse skin still poorly simulates human skin by several metrics; compared with human skin, healthy mouse skin contains thin epidermis, lacks rete ridges, and adheres loosely to underlying tissues. Furthermore, though S. aureus is demonstrably able to scavenge metabolically critical iron from murine haemoglobin, it does so with far less efficacy than from human haemoglobin, further limiting the simulative accuracy of this model as a proxy for human S. aureus‐infected wound healing. Additionally, since we utilised the wound bed samples in this project solely for isolation of RNA for transcriptomic analysis, we were unable to characterise wound healing in the samples from which the RNA was derived, preventing us from correlating gene expression profiles with phenotypic effects or healing outcomes in host wounds. Previously, our lab has demonstrated that S. aureus‐infected wounds in a polygenic diabetic mouse model (TallyHo mice) demonstrated decelerated healing and increased bacterial burden compared with non‐diabetic (SWJ) mice, suggesting that db/db mouse infected wounds in the experiments described in this manuscript likely also healed more slowly and with potentially greater bacterial burden, though this will need to be confirmed or refuted with follow‐up experiments. Though potential differences in bacterial burden among groups are not technically responsible for differential gene expression profiles observed here, as statistical calculations for differential gene expression involve per‐sample normalisation based on the proportions of reads mapping to the S. aureus genome as well as library size, differences in bacterial densities among experimental groups may still contribute to phenotypic changes in situ, the transcriptomic signatures of which may then be detected at experimental endpoints by RNA‐seq. Another limitation of the present study is the use of wild type C57BL/6 mice as control animals, rather than the more appropriate heterozygote (db/+) mice frequently used as non‐diabetic control animals in experiments utilising db/db mice. Previously, it has been demonstrated that wounds of db/db mice infected with MRSA heal significantly more slowly compared with heterozygous (db/+) control mice. This suggests that the aberrations that we observed in db/db mouse wound healing relative to wild type C57BL/6 mice were unlikely to be due to the use of wild type rather than heterozygous control mice, though future experiments will be needed to confirm this as well.

Despite these limitations, we believe that our analysis begins to shed light on paradigms in support of differential pathogen adaptation to the wound environment that are worthy of further investigation. Ongoing work seeks to validate our findings and uncover additional dysregulated gene expression paradigms in experimental models more analogous to human wound healing. Future work will also aim at developing a more complete understanding of the factors characteristic of the diabetic wound environment that drive phenotypic changes in the bacteria that colonise them, as well as the consequences of persistent adaptations to the diabetic wound environment on healing.

AUTHOR CONTRIBUTIONS

Kai Leung, Thomas Mustoe, Seok Jong Hong, and Robert Galiano conceived of the project and designed experiments. Kai Leung, Thomas Mustoe, Seok Jong Hong, and Robert Galiano provided oversight and expert input. David Dolivo, Steven Lanier, and Seok Jong Hong performed experiments and analysed data. David Dolivo drafted the manuscript and designed the figures. All authors approved of the final version of the manuscript.

CONFLICT OF INTEREST

All authors have no conflicts to report.

Appendix S1 Supporting Information Methods.

Click here for additional data file.

Figure S1 UAMS‐1 and M. musculus mapping statistics. Box and whisker plots representing the percentages of processed reads that mapped 0, 1, or >1 times to the S. aureus UAMS‐1 genome (A) or the M. musculus mm10 genome (B). n = 5 samples per group. Error bars = standard deviation. Statistical analysis for each group was performed using a one‐way ANOVA and Tukey's post‐hoc pairwise comparison. *P < 0.05, **P < 0.01, ***P < 0.001. For clarity, statistical comparisons to planktonic cultures are omitted due to triviality. BL6 (wild type), db (db/db), Pl (planktonic).

Click here for additional data file.

Figure S2 Expression of contaminating murine inflammatory genes. Reads for all wound bed samples (n = 20) were mapped to the M. musculus genome (mm10 construction). Aligned reads were extracted and transcripts per million (TPM) values for each gene for each sample were estimated using Stringtie. Violin plots were constructed to visualise the expression distributions of selected murine inflammatory genes expressed relatively highly across all wound bed samples (n = 20). Inner dotted lines denote medians and outer dotted lines denote quartile demarcations.

Click here for additional data file.

Figure S3 Principle component analysis comparing wound bed and planktonic samples. Principal component analysis comparing gene expression profiles across all samples.

Click here for additional data file.

Figure S4 Gene expression heatmap visualising wound bed and planktonic sample comparisons. Gene‐clustered heatmap demonstrating Z scores calculated among all samples mapped to the UAMS‐1 genome. Genes pictured are all of those significantly differentially expressed between planktonic (n = 5) and wound bed (n = 20 samples, 4 groups of 5 samples each). Significantly differentially expressed genes are defined as FDR P < 0.05 as computed by the Benjamini and Hochberg method.

Click here for additional data file.

Table S1 DEGs of wound bed samples relative to planktonic samples.

Click here for additional data file.

Table S2 Gene ontology analysis (biological process) of DEGs of wound bed samples relative to planktonic samples.

Click here for additional data file.

Table S3 DEGs of db/db versus wild type wound bed samples at day 3.

Click here for additional data file.

Table S4 Gene ontology analysis (biological process) of DEGs of db/db versus wild type wound bed samples at day 3.

Click here for additional data file.

Table S5 DEGs of db versus BL/6 wound bed samples at day 7.

Click here for additional data file.

Table S6 Gene ontology analysis (biological process) of DEGs of db versus BL/6 wound bed samples at day 7.

Click here for additional data file.

Table S7 Shared day 3 and day 7 DEGs of db/db versus wild type wound bed samples.

Click here for additional data file.

Table S8 Gene ontology analysis (biological process) of shared day 3 and day 7 DEGs of db/db versus wild type wound bed samples.

Click here for additional data file.