A metatranscriptomic approach to explore longitudinal tissue specimens from non-healing diabetes related foot ulcers.

Cellular mechanisms and/or microbiological interactions which contribute to chronic diabetes related foot ulcers (DRFUs) were explored using serially collected tissue specimens from chronic DRFUs and control healthy foot skin. Total RNA was isolated for next-generation sequencing. We found differentially expressed genes (DEGs) and enriched hallmark gene ontology biological processes upregulated in chronic DRFUs which primarily functioned in the host immune response including: (i) Inflammatory response; (ii) TNF signalling via NFKB; (iii) IL6 JAK-STAT3 signalling; (iv) IL2 STAT5 signalling and (v) Reactive oxygen species. A temporal analysis identified RN7SL1 signal recognition protein and IGHG4 immunoglobulin protein coding genes as being the most upregulated genes after the onset of treatment. Testing relative temporal changes between healing and non-healing DRFUs identified progressive upregulation in healed wounds of CXCR5 and MS4A1 (CD20), both canonical markers of lymphocytes (follicular B cells/follicular T helper cells and B cells, respectively). Collectively, our RNA-seq data provides insights into chronic DRFU pathogenesis.

Acute wounds in healthy people proceed through an orderly process to repair injured tissue. People with diabetes are at higher risks of impaired wound healing due to several ill‐defined immunological disturbances [1, 2] and multiple comorbidities [3]. Over the last quarter century, there has been a wealth of data and development in the field of wound repair and regeneration in context to human skin wounds [1]. Despite this, there has been slow progress in translating this data into methods for reducing the number of people who develop chronic wounds and or in drastically improving wound healing rates. There has been limited advance in identifying host molecular or microbial markers that can distinguish healing from non‐healing diabetes related foot ulcers (DRFUs), infected vs non‐infected DRFUs, or in providing clarity on why some DRFUs heal and others do not. Historically research has focussed on the cellular/biological mechanisms of wounding (normal healing and impaired healing) in vitro, ex vivo, and in animal models [2], with scant data from studies using in vivo human samples [3, 4, 5, 6, 7, 8].

More recently, the utilization of molecular techniques including quantitative polymerase chain reaction and RNA microarrays have allowed researchers to explore differentially expressed genes (DEGs) involved in the reparative process of wound healing. This has provided deeper insights into understanding which host genes are up‐regulated or down‐regulated under specific wounding conditions in various animal models, and to a lesser extent in human samples. The major limitation of this approach is that these techniques have often relied on a priori assumptions about host genes of interest. Microarray techniques explore known sets of RNA and are limited to analysis of a relatively small number of genes. To circumvent these limitations, we employed whole transcriptome sequencing (RNA‐seq) of human chronic DRFU tissue and control healthy foot skin. Identifying gene expression underlying cellular activities (e.g. life‐cycle processes, responses to environmental cues) and enrichment of biological pathways is crucial for improving our understanding of why some wounds heal and others do not. We explored host and microbial RNA‐seq libraries and performed a temporal analysis of DRFUs receiving active treatment with conservative sharp debridement and application of a topical concentrated surfactant gel (CSG).

MATERIALS AND METHODS

Participant population

Participants from a previously reported clinical study [9] had additional tissue specimens obtained for the purpose of performing total RNA‐sequencing. Individuals aged over 18 years presenting with a chronic DRFU were prospectively recruited over an 18‐month period. Participants were eligible if they had a non‐healing neuropathic or neuroischaemic DRFU in the presence of standard care. In total, 10 patients with non‐healing DRFUs were recruited over the study period.

Specimen collection

Tissue specimens were obtained via dermal curettage (4‐mm dermal ring curette, Kai Medical) from the ulcer base and adjacent to the leading edge of each ulcer after debriding and cleansing with NaCl 0.9%. Tissue specimens were collected at baseline (week 0), mid‐point (week 3) and end of treatment (week 6). All tissue specimens were immediately placed in to RNAlater (Thermo Fisher Scientific, Waltham, MA, USA) stabilization solution for 24 h at 4°C and then stored at −80°C until processed.

Study eligibility

Eligible DRFUs included wound – grades 1 and 2 (excluding exposed deep structures or bone involvement), ischaemia – grades 0–2, infection grade 0, as per the risk stratification of Wound, Ischaemia, and foot Infection classification (WIfI) [10]. A chronic DRFU was defined as being > 6 weeks in duration and failing to respond to standard care, in addition to no observed changes in wound metrics over a lead in period of four consecutive weeks prior to enrolment. Standard of care was defined as being weekly treatments by a podiatrist performing appropriate wound bed preparation through sharp conservative debridement or curettage, wound cleansing with NaCl 0.9%, and the use of a non‐adherent absorbent wound dressing. Offloading of plantar DRFUs was through a removable cast boot (DH Offloading Walker®, Össur, Australia) and for non‐plantar DRFUs, a post‐operative shoe (Darco all‐purpose boot®, Össur, Australia).

Exclusion criteria were clinical signs and current diagnosis of an acute infection as per the International Working Group on the Diabetic Foot – diabetic foot infection guidelines [11], known diabetic foot osteomyelitis that was associated with an ulcer, or patients who had received any new topical or systemic antimicrobial therapy 2 weeks prior to enrolment. In addition to standard of care, patients received the study topical agent (concentrated SG, Plurogel®, Medline Industries Inc), which was applied to the surface of DRFUs every 2 days up to a maximum of 6 weeks. A non‐adherent absorbent dressing pad was applied over the topical CSG for exudate management.

Patient level data

Broad demographic data was collected at each timepoint on wound metrics using 3D wound imaging (eKare Inc, Fairfax, VA, USA). Data included length × width (cm), wound area (mm2), depth/volume (mm3), % reductions in wound area between time points (Tables S1 and S2). Participants were grouped into healing or non‐healing DRFUs. We defined healers as those participants having demonstrated cumulative reductions in the wound area size over the duration of the study period. Non‐healers as being participants who demonstrated no reduction in wound area over the duration of the study, or, who demonstrated an increase in the wound area size from baseline. P4 was defined as a non‐healer despite showing a 3.6% reduction in wound size from baseline to end of treatment. We decided this rate of healing over a six‐week period demonstrated a poor response to healing and the 1.2% change in wound surface area per timepoint could be attributed to an error rate in wound measurements alone.

RNA isolation and library preparation

Total RNA was isolated from tissue specimens using TRIzol reagent (Thermo Fisher Scientific). Briefly, samples were homogenized in TRIzol using a Tissue Lyser II (Qiagen, Hilden, GER) with 0.1‐ and 0.5‐mm beads for 2 × 5 min cycles at 30 Hz. Chloroform was added and the samples was centrifuged to isolate the RNA containing aqueous phase. RNA was then precipitated from the aqueous phase using Isopropanol and washed with 70% ethanol before being resuspended in nuclease free water and stored at −80°C. RNA samples were then submitted to the Next Generation Sequencing Facility at Western Sydney University for library preparation using the Zymo‐Seq RiboFree Total RNA library kit (Zymo Research, Irvine, CA, USA). Sequencing was then carried out on the Novaseq 6000 S4 platform (Illumina, CA, USA) at 2 × 150 bp to ensure an output of >200 million reads per sample.

Processing of metatranscriptomic data

Due to poor RNA yield from some DRFU tissue specimens, only a subset of the total patient population (6/10 samples) could be used for RNA‐seq analysis. RNA‐seq generated approximately 150 million (~20 SEM) reads per sample with a mean ribosomal content of 13.59%. Reads were trimmed using TrimGalore/Cutadapt [12] and aligned in paired‐end mode to GRCh38.p12 with alternative haplotypes and unlocalised contigs removed, using STAR2.5.4b [13]. Gene and isoform‐level quantification was performed using RNA‐Seq by Expectation Maximization (RSEM) [14]. Only samples that had a human library size of >30 million reads were retained for analysis, with samples exhibiting poor mapping likely being reflective of high microbial content. Following quality filtering, the mean effective library size which mapped to the human genome was approximately 61 million reads per sample, providing sufficient coverage for robust host RNA profiling. Unmapped reads were retained for microbial RNA analysis.

Annotation of metatranscriptome datasets

For host analysis, a fractional factorial design was fitted using EdgeR [15] to all 18 DRFU samples with timepoint as an ordered factor, wound healing as a binary factor and reduction in wound area as a covariate and proxy for patient ID. All DEGs were calculated using quasi‐likelihood F‐tests from the same model as, for the pairwise timepoint comparisons, a set of post‐hoc contrasts, and as a linear interaction for the second order effect between ordered timepoint and wound healing. DEGs were defined as having a FWER <0.05 and absolute log2 fold change >2.

Microbial analysis of RNA data was completed using the SqueezeMeta pipeline (v1.3) [16] utilizing the co‐assembly option with no binning and singletons included. Briefly, paired end reads were assembled using RNAspades prior to taxonomic and functional annotation using the DIAMOND sequencing aligner to the GenBank and KEGG databases, respectively. Reads of individual samples were then be mapped to assembled contigs for the estimation of taxonomic and functional abundances using Bowtie2 [17]. Relative RNA plots for each DRFU phenotype were then generated using R, based on the raw read counts mapping to each taxon. Raw read counts varied across samples, ranging from 3 to 13 million reads available for downstream microbial RNA profiling. Taxonomic and functional analysis was completed in parallel using a Bioconductor software package for examining differential expression of replicated count data (EdgeR) and the SQMtools package. Read counts were normalized (LogCPM), and a PCA was performed to demonstrate any variation among datasets based on their taxonomic and functional profiles. Data was subset using the subsetFun() function within the SQMtools package.

Statistical analysis

The R Statistical Package (R Core Development Team, 2017) was used to generate all figures and compute statistical analysis.

RESULTS

Differential expression analysis of DRFU biopsies

Ten patients were recruited from a high‐risk foot service as part of a previously published study on chronic DRFUs [9]. Punch biopsies from foot ulcer tissue were obtained prior to treatment (baseline), after 3 weeks of treatment (midpoint) and after 6 weeks of treatment (study endpoint). Demographic and clinical metadata were collected, in addition to wound metrics obtained at each timepoint (Tables S1 and S2). Biopsies were also collected of healthy plantar foot skin in two persons with diabetes and were used as controls. Total RNA sequencing was performed which produced a median effective library size range of 21–80 million read pairs per sample for host differential expression analysis. Reads which failed to map to a human reference genome were retained for microbial analysis which yielded library sizes of 3–13 million read pairs. Four patient datasets were omitted from the final analysis due to missing sample timepoints; in these cases, low RNA yields resulted in library failure.

Microbial taxonomic and functional activity within chronic DRFUs were explored through Principal Coordinate Analysis (PCoA) which showed that midpoint and endpoint samples generally cluster closer together compared to baseline, both in the context of microbial taxonomy and function (Figs S1 and S2). Initially, we aimed to obtain a synoptic overview of all human RNA samples to examine their dissimilarity based on treatment time (baseline, midpoint, and endpoint) relative to controls. Normalized read counts (log counts per million – CPM) were used to perform principal component analysis (PCA) (Fig. S3). This identified that control healthy foot skin is highly dissimilar to chronic DRFUs. We then sought to identify DEGs unique to chronic DRFUs (baseline samples) relative to controls. An unpaired comparison model was performed with EdgeR [10] on normalized read counts, which identified 306 DEGs with a LogFC > 2 and FWER < 0.05 (Fig. S4), with 72 genes upregulated and 234 genes downregulated (Data S1).

Chronic DRFU microbiome is heterogenous

Relative microbial RNA transcripts (%) for DRFU samples can be seen in Fig. 1. In baseline chronic DRFU samples, Streptococcus, Staphylococcus, Streptomyces, and Corynebacterium demonstrated the highest relative microbial activity. When performing a temporal analysis of chronic DRFU samples, we commonly observed small shifts in the active communities from baseline (the onset of treatment) to midpoint and from mid‐point to end of study (Fig. 1). These were from within‐sample taxa that either increased or decreased in relative RNA transcripts. The most striking change in relative RNA was observed in the case of P5 (non‐healing DRFU), who demonstrated an increase in the proportion of Streptococcus RNA by end of treatment.

Fig. 1

A Bar chart showing the relative microbial RNA transcripts of the top 15 bacteria which were stratified at the genus level. *(B = baseline [week 0], M = mid‐point of treatment [week 3], E = end of treatment [week 6]).

To determine the functional activity of microbial communities in chronic DRFUs, EdgeR was utilized for DEG analysis of all KEGG annotated transcripts. However, no DEGs were identified between chronic DRFU treatment timepoints.

Signatures for mechanical stress and microbial presence provide stimuli to the host immune system

Analysis of DEGs enriched within baseline chronic DRFUs relative to control skin identified several genes involved in cell stress, injury, and signalling. This included upregulation of signalling molecules associated with mechanical stress; ITGB3 (5‐LogFC), ITGA5 (4.7‐LogFC) and pathogen associated molecular patterns (PAMPS); S100A4 (4.4‐LogFC), FPR1 (5.2‐LogFC), FPR2 (7.7‐LogFC) (Data S1).

Early adaptive immune responses, inflammation and fibrosis are predominant features in chronic DRFUs compared to healthy skin

Pro‐inflammatory mediators were a defining feature of chronic DRFUs (Fig. 2A,B). A large subset of enriched genes associated with the immune response were observed including CXCL5, PROK2, PTPRN, OSM, PADI4, CSF3, IL1B, and CXCL3, which were significantly upregulated with a LogFC > 10. In baseline chronic DRFU samples, gene set enrichment analysis (GSEA) on Gene Ontology (GO) [11] identified enrichment of immune biological processes including Inflammatory response, TNF signalling via NFKB, IL6 JAK‐STAT3 signalling, IL2 STAT5 signalling and reactive oxygen species pathway (Fig. 2B). We further observed DEGs potentially involved in fibrosis via dysregulation of injury‐triggered Epithelial‐to‐Mesenchymal Transition (EMT) including up‐regulation of CXCL1 (7.5‐LogFC), CTGF (4.9‐LogFC), ADAM12 (4.1‐LogFC), TNFRSF12A (4‐LogFC), ITGB 3 (5‐LogFC) and ITGA 5 (4.7‐LogFC). Enrichment of Hallmark gene sets also identified EMT as the most enriched biological process, in addition to enrichment of Transforming Growth Factor Beta (TGF‐β), a principal regulator of EMT (Data S2).

Fig. 2

RNA‐seq identifies unique transcriptome profiles of DRFUs relative to healthy control tissue. A heatmap identifying the top 30 DEGs between control and chronic wound (baseline) tissue samples with clustering performed using ward. Colour and intensity of genes expressed is based on LogCPM counts. B top 20 pathways identified to be significantly enriched within baseline samples based on GSEA using gene ontology (FWER < 0.05).

In baseline chronic DRFUs, 234 genes were downregulated with the majority of these being associated with skin development and homeostasis, including KRT2 (16.1‐LogFC), KPRP (16.1‐LogFC), DSC1 (15.7‐LogFC), KRT9 (14.6‐LogFC), KRT1 (14.2‐LogFC), SERPINA12 (13.9‐LogFC), KRT10 (11.6‐LogFC), and FOXN1 (10.2‐LogFC). Additionally, GO gene set enrichment analysis identified the enrichment of pathways associated with Skin development, Epidermis development, Epidermal cell differentiation, Keratinocyte differentiation and Keratinization (Data S3).

Treatment of chronic DRFUs with conservative sharp debridement and a topical concentrated surfactant gel (CSG) alter the host transcriptome

We next examined longitudinal changes to the host metatranscriptome resulting from treatment of DRFUs with conservative sharp debridement and a topical CSG. We identified several upregulated DEGs at midpoint compared to baseline across the patient cohort, the most significant being RN7SL1 (16.7‐LogFC), IGHG4 (14.5‐LogFC), ATP6V1G2‐DDX39B (15.2‐LogFC) and CXCL13 (6.7‐LogFC) (Fig. 3A and Data S4).

Fig. 3

Temporal DEG analysis of all study timepoints including a baseline vs midpoint, B midpoint vs endpoint and C baseline vs endpoint. For each comparison, the top 4 DEG LogCPMs are shown and distinguished by treatment stage. DEGs were considered significant if LogFC>2 and FWER<0.05.

Interestingly, we noted that signal recognition particle (RN7SL1) expression was limited exclusively to the midpoint of treatments, with the exceptions being P5 and P9 (Fig. 3A,B). IGHG4 was also notable, as its expression was higher in individuals who were healing by end of treatment (P1, P9, P10) while predominantly absent in patients who were not healing (P4, P5, P8). DEG analysis between endpoint and baseline‐midpoint identified DE genes associated with skin development, the most significant belonging to the Keratinocyte family (KRT6A; 10.3‐ LogFC, KRT6C; 10.8‐ LogFC, KRT14; 8.5‐LogFC, KRT6B; 6.8‐LogFC, KRT16; 6.7‐LogFC, KRT5; 7.2‐LogFC, KRT17; 6.6‐LogFC) (Fig. 3C and [Link], [Link]).

GSEA analysis was also performed to identify cellular processes which were enriched as treatment progressed. When comparing baseline to midpoint, we identified that the GO humoral immune response pathway was significantly enriched within the midpoint samples (NES = 1.92, FWER < 0.05). When comparing the endpoint to baseline‐midpoint, we again observed that pathways associated with skin homeostasis were significantly enriched within the endpoint samples, in addition to pathways associated with the humoral immune response (Figs S5–S7). Hallmark GSEA also revealed several pathways associated with the immune response were enriched within endpoint samples, including (1) Interferon Gamma Response, (2) P53 Pathway, (3) TNFA signalling via NFKB and (4) Interferon alpha response (Data S7).

B cell activity differentiates healing DRFUs compared to non‐healing

We next aimed to identify DEGs and cellular pathways between healing and non‐healing DRFUs in response to treatment. We tested the interaction effect between the binary healing response and all treatment stages as an ordered factor. The top 30 DEGs between healing and non‐healing can be seen in Fig. 4A. The most significant DEGs showing progressive upregulation in healing DRFUs included CXCR5 (59‐LogFC) and MS4A1 (51.6‐LogFC), which are canonical markers of follicular helper T cells, follicular B cells and memory B cells. Functionally, these genes have essential roles in B cell migration and the anti‐pathogen antibody response. Immunoglobulin antibodies IGHA1 (6.4‐LogFC) and IGKV3 (9.5‐LogFC) were also observed as being enriched in healing DRFUs. We further identified several long non‐coding RNAs AC244154.1 (46.6‐ LogFC), AC005943.1 (13.6‐LogFC), AC106886.5 (6.7‐LogFC), AC068234.1 (11‐LogFC), AL163636.2 (42.2‐LogFC). The most significant downregulated genes included RNVU1‐7 (9.6‐LogFC), KRT6C (5.8‐LogFC) and Pl15 (4.6‐LogFC) (Fig. 4B).

Fig. 4

Linear interaction effect of healing/non‐healing and timepoint. A heatmap of top 30 DEGs. B DEGs between healing and non‐healing patients. Blue/labelled points represent the top 20 DEGs based on FWER values, red points represent all other DEGs with FWER <0.05 and grey points represent non‐significant DEGs. C pathways identified to be associated with non‐healing patients as identified through GSEA utilizing GO and Hallmark gene sets (FWER<0.05).

GSEA did not identify any significant pathways associated with healing DRFUs, however several pathways were enriched in non‐healing DRFUs (Fig. 4C). GSEA using GO ontology identified pathways associated with mitochondrial function and oxidative phosphorylation, with the most significantly enriched including (i) the mitochondrial protein complex; (ii) oxidative phosphorylation; (iii) inner mitochondrial membrane protein complex; (iv) positive regulation of ligase activity and the (v) anaphase promoting complex dependent catabolic process. Additionally, GSEA using Hallmark Gene sets identified several enriched pathways associated with the immune response, including (i) IL6 JAK STAT3 signalling; (ii) interferon gamma response and (iii) interferon alpha response.

DISCUSSION

Chronic DRFUs remain a significant comorbidity of diabetes worldwide and a major contributing factor to lower limb amputations and reduced quality of life. Numerous studies have identified persistent inflammation as a key feature of chronic DRFUs [18, 19, 20], with the causative mechanisms often being attributed to tissue injury by increased plantar pressures/shear [21] and/or by microbial involvement [22]. Recent studies have shown increasing evidence of the presence of biofilms within human tissue and their tolerance to therapeutics and the host response [23]. We have previously utilized microscopy and molecular approaches to reveal that biofilms are ubiquitous in DRFUs but have yet to define how they contribute to persistent infection and delayed wound healing [24]. Previous evidence has proposed that biofilms provide constant stimuli to the immune system, perpetuating a state of chronic inflammation [25, 26].

The six individuals with chronic DRFUs in this study all had previously confirmed microbial aggregates/biofilms within DRFU tissue using scanning electron microscopy (SEM) and peptide nucleic acid fluorescence in situ hybridisation (PNA‐FISH) [9]. All individuals were provided with offloading devices (removable CAM walkers), but compliance with activity reduction or with wearing the offloading device was not quantified. We provide a medical illustration as a simplistic overview of our findings (Fig. 5) but also acknowledge this may not be generalisable to all chronic DRFUs. In this study, we identified differentially expressed genes in chronic DRFUs relative to healthy skin which broadly functioned in response to tissue injury‐damage, and pathogen stimuli/recognition. We further identified both DEGs and enrichment of pathways involved in early adaptive immune and inflammatory responses and fibrosis.

Fig. 5

Medical illustration summarizing the key events inferred from gene expression in six chronic DRFUs and the differentiation between healing and non‐healing ulcers.

We observed enrichment of integrins (ITGB3 and ITGA5) which can function as mechanosensors in detecting mechanical stress in cells [27]. We also observed enriched DEGs associated with recognizing stress signals that include damage/pathogen associated molecular patterns (DAMPS/PAMPS), including S100A4, FPR and FPR2. The pattern recognition receptors FPR1 and FPR2 are expressed on several immune cells including neutrophils and macrophages, as well as non‐immune cells including endothelial cells and keratinocytes [28]. Formyl peptide receptors (FPR) can bind host derived DAMPs and microbial derived PAMPs, due to the similarities between mitochondrial and bacterial N‐formyl peptides [29]. FPR2 also binds host derived cathelicidin (LL37) from granules as well as several microbial derived ligands [30, 31]. Activation of FRPs by DAMPs/PAMPs results in increased neutrophil migration, phagocytosis, and the release of superoxide at the site of injury which all contribute to tissue damage [32, 33].

Consistent with a response to mechanical stress and microbial presence, we observed enrichment of pro‐inflammatory cytokines, PROK2, PTPRN, OSM, CSF3, G‐CSF, IL1B, and a matrix metallopeptidase, MMP9. Several receptors and chemokine ligands were also enriched, including CXCL1, CXCL3, CXCL5, CCR1, CCL18 and CSF3R. Enrichment of IL1B and OSM has been reported as critical to the pathogenesis of hypertensive leg ulcers, potentially by catalysing the destruction of reconstructed epidermal tissue [34]. GSEA of enriched genes within DFRUs identified the enrichment of pathways associated with the immune and inflammatory responses. In total, GO enrichment identified 98 significantly enriched pathways in baseline chronic DRFUs, with the majority involved in cell chemotaxis and the inflammatory response (Fig. S8). The most significantly enriched pathways included GO Leukocyte Migration, GO Chemokine Mediated Signalling and GO Response to Interferon Gamma. Collectively, these results are indicative of mass neutrophil and macrophage recruitment within DRFU tissue [35, 36, 37, 38].

We also observed enrichment of PADI4 and CSF3, suggesting the formation of neutrophil extracellular traps and/or neutrophils undergoing NETosis to produce oxidative bursts within DRFUs [39]. The formation of NETs/NETosis has been previously shown through in vitro models to be driven by the presence of biofilm [40, 41, 42]. Recently, Papayannopoulos [43] demonstrated that P. aeruginosa biofilm formation induces NET release by infiltrating neutrophils, and in turn, NETs amplified biofilm formation. Several studies have also identified NETs as potential markers of impaired wound healing within DRFUs [44, 45]. Studies in humans and murine models have shown that hyperglycaemia predisposes neutrophils to release NETs when stimulated with bacterial lipopolysaccharides (LPS) [46]. In addition, protein expression of PADI4 was found to be increased fourfold in neutrophils from individuals with diabetes when compared to those from healthy controls, suggesting that increased PADI4 may favour chromatin decondensation [47].

Epithelial‐to‐Mesenchymal transition is necessary for physiological tissue repair, though continued transition from epithelial cells to myofibroblasts may result in pathological fibrosis [42]. Several DEGs from baseline chronic DRFUs may be associated with pathological EMT in our dataset. We observed down‐regulation of epithelial cadherin CDHR1 (9.9‐LogFC) and desmocollin DSC1 (15.7‐LogFC). Both are essential for the mechanical integrity of tissue and their depletion is one of the main initiation events of EMT [47]. Cell progression towards a mesenchymal phenotype is mediated through expression of proteins including (but not limited to) Integrins. We previously identified upregulation of integrins ITGB 3 (5‐LogFC) and ITGA 5 (4.7‐LogFC) as mechanosensors of stress, however Integrin signalling has also been implicated in EMT [48]. We further identified upregulation of connective tissue growth factor (CTGF; 4.9‐LogFC) which is responsible for modulating signalling pathways including myofibroblast activation and extracellular matrix deposition and remodelling [49]. Continued upregulation of CTGF as noted in baseline chronic DRFUs, could indicate pathologically prolonged myofibroblast activity and resulting fibrogenesis. Hallmark gene sets further identified Epithelial‐to‐Mesenchymal Transition as the most enriched biological process (NES = 2.47, FWER < 0.000), with supportive enrichment of TGF‐β (NES = 1.70, FWER = 0.005), a principal regulator of EMT.

We conducted a temporal analysis of the host and microbial transcriptomes of chronic DRFUs following treatment with conservative sharp debridement and a topical CSG, to identify host/microbial factors associated with wound healing. DEG analysis identified significant enrichment of RN7SL1 and IGHG4 (immunoglobulin heavy constant gamma 4) within midpoint samples relative to baseline. RN7SL1 is the RNA component of the signal recognition particles and functions primarily in facilitating protein secretion [50]. Upregulation of RN7SL1 may play a role in negatively regulating the immune response through the suppression of p53 translation, which is supported by the enrichment of the p53 signalling pathway in endpoint samples relative to midpoint, in conjunction with the decreased expression of RN7SL1 [51] (Data S3). In response to stress signals, p53 is activated in a specific manner by post‐translational modifications that induce cell senescence or cellular apoptosis [52]. We observed that RN7SL1 upregulation is predominantly limited to midpoint samples except for P9 which showed continual upregulation to the end of treatment. Potential suppression of p53 through RN7SL1 and enrichment of IGHG4 may indicate that treatment with conservative sharp debridement and a topical CSG were successful in removing or destabilizing microbial communities which is supported by our previous work [9]. The downstream effects of this may be altered pathways supressing cellular senescence and apoptosis and favouring cell activity/proliferation.

A major aim of the study was to ascertain if healing DRFUs followed a longitudinal gene expression pattern distinct from non‐healing DRFUs. Analysis identified CXCR5 as the most significant progressively upregulated DEG in healing DRFUs. We also detected upregulation of CXCL13 and its receptor CXCR5 in midpoint samples, confirming that B lymphocyte activity was induced at this time [53] and continued in those individuals with healing DRFUs. CXCR5 is the only known receptor for CXCL13 and is expressed by follicular B cells, helper T cells (Tfh) and T follicular regulatory cells (Tfr) [54]. Tfh cells are primarily responsible for germinal centre (GC) formation and instigate the differentiation of B cells to permablasts and plasma cells for the secretion of antibodies. Tfr cells regulate Tfh cell proliferation and promote selection of high affinity B cells [55]. CXCL13 and its receptor CXCR5 are pivotal in lymphoid neogenesis, and we postulate that enrichment of CXCR5 may support secondary lymphoid tissue formation and B cell homing within chronic DRFUs [56, 57].

MS4A1 encoding for CD20 was found to be significantly enriched within healing patients and is expressed on the surface of B lymphocytes. CD20 is required for T‐independent humoral immunity by acting as a calcium channel for B cell activation, making it a canonical marker for memory B cells [58, 59, 60]. Several genes (IGHA1 and IGV3‐20) associated with heavy and light chain immunoglobulins respectively, were also enriched within healing patients. Heavy and light chain immunoglobulins play an essential role in the activation of immature IgM expressing B cells from pro‐B cells [60]. IgM expression plays a pivotal role in B cell differentiation/maturation within the GC as well as responses to neutrophil derived cytokines in the early immune response [60]. Combined, this data suggests that B cell maturation correlates with healing DRFUs, supporting previous findings by Sirbulescu et al [61]. It could further suggest that in non‐healing DRFUs, there is a potential dysfunction in B cell development that maybe directly associated to the diabetes [62].

GSEA analysis identified several GO and Hallmark pathways correlating with non‐healing DRFUs. These were predominantly associated with mitochondrial activity, oxidative phosphorylation, and the immune response (Fig. 4C). Oxidative phosphorylation is the primary mechanism by which ATP is generated within cells and occurs within mitochondria. A common by‐product of oxidative phosphorylation which is essential for the immune response and wound repair is the production of ROS to clear invading bacteria [63]. The downstream effects of prolonged ROS production are increased oxidative stress within host tissues and resultant damage to cell DNA, proteins, and lipids [62, 64]. Our data is supportive of previous work in animal models suggesting that prolonged oxidative stress is a major contributor to DRFU chronicity and a key indicator for non‐healing [65].

Study limitations

In this study we did not use PCR to validate DEGs of interest. Our rationale against doing this is provided in a recent review article by Coenye [66] who discusses the argument of how reliable RNA‐seq is in identifying DEGs, and if qPCR is needed to validate such expression differences. Coenye, provides a balanced argument detailing the historic validation of genome‐scale expression studies which stems from prior work with microarrays, and concerns about reproducibility and bias. RNA‐seq does not suffer from the same issues with several studies specifically addressing the correlation between results obtained with RNA‐seq and qPCR. Additionally, we did not identify any microbial differentially expressed genes between DRFU samples. This is likely due to the small sample size combined with the lower coverage of microbial reads in tissue. Ideally, future studies should include a larger sample size that are sequenced to a significantly higher depth of greater than 500 million reads per sample. However, we do note that costs associated with achieving these sequencing depths are high and will likely limit the number of samples that can be sequenced.

Collectively, our RNA‐seq data provides important insights into chronic DRFU pathogenesis. Larger datasets that included social/environmental determinants, clinical and laboratory data would be valuable to increase the generalisability of our data.

The authors would like to acknowledge the extended team involved in patient care at the Liverpool Hospital High Risk Foot Service. The authors would also like to thank the infrastructure support provided by South Western Sydney LHD, Ingham Institute of Applied Medical Research, the Next Generation Sequencing Facility at Western Sydney University and the Advanced Materials Characterization Facility (AMCF) at Western Sydney University. Open access publishing facilitated by The University of Sydney, as part of the Wiley ‐ The University of Sydney agreement via the Council of Australian University Librarians. [Correction added on 18 May 2022, after first online publication: CAUL funding statement has been added.]

Conflict of interest

MM was provided with an investigator‐initiated grant by Medline Industries to cover the costs of RNA‐seq. All other authors have nothing to declare.

Figure S1. PCA of microbial taxonomic transcripts across all samples.

Figure S2. PCA of microbial functional transcripts across all samples.

Figure S3. PCA plot of normalised counts (LogCPM) for all samples.

Figure S4. MD plot and Volcano plot of Differentially Expressed Genes (DEGs). DEGs between healthy control tissue and baseline chronic DRFU tissue (LogFC >2, FWER<0.05).

Figure S5. Treemap of GO pathways enriched within midpoint samples relative to baseline.

Figure S6. Treemap of GO pathways enriched within endpoint samples relative to midpoint samples.

Figure S7. Treemap of GO pathways enriched within endpoint samples relative to baseline samples.

Figure S8. Treeplot of GO pathways enriched within baseline samples relative to healthy controls.

Table S1. Patient demographics and clinical data.

Table S2. Wound metrics obtained using a 3D camera to track wound area, length x width and depth.

Click here for additional data file.

Data S1 Baseline Chronic DRFU vs Control DEGs.

Click here for additional data file.

Data S2 Baeline Chronic DRFU Hallmark GSEA.

Click here for additional data file.

Data S3 Baseline Chronic DRFU GO pathways.

Click here for additional data file.

Data S4 Baseline vs Mid‐point of treatment DEGs.

Click here for additional data file.

Data S5 Mid‐point vs End of treatment DEGs.

Click here for additional data file.

Data S6 Baseline vs End of treatment DEGs.

Click here for additional data file.

Data S7 End of treatment to baseline‐midpoint Hallmark GSEA.

Click here for additional data file.

Data S8 End of treatment to baseline‐midpoint GO pathways.

Click here for additional data file.