A longitudinal study of the diabetic skin and wound microbiome

Study design, ethics approval, and sample collection

Ethical approval for the study was obtained from both the University of Technology Sydney Human Research Ethics Committee (approval number 2013000170), and the Western Sydney Local Health District Human Research Ethics Committee (approval number HREC2013/9/5.3(3809) AU RED LNR/13/WMEAD/294). Diabetic individuals and control subjects provided written consent for sample collection and all subsequent analyses.

Diabetic adults (n = 8) (Table 1) were selected for inclusion in the study based on medical diagnosis of type II diabetes, the presence of a chronic wound on one foot (chronic wound = present for six or more weeks) and no antibiotic therapy within the previous four weeks. Three swabs were collected for each diabetic subject every two weeks for a 10 week period using sterile rayon tipped swabs (Copan) that had been pre-moistened with a sterile solution of 0.15 M NaCl and 0.1% Tween 20. Two skin swabs were collected from intact foot skin (1) adjacent to the chronic wound (skin adjacent, SA) and (2) contralateral site to the chronic wound (skin contralateral, SC). Skin swabs were collected by firmly rubbing the moistened swab over the base of the foot skin surface for a period of 30 seconds. The whole base of the foot was used to maximise the DNA yield. Skin swab samples were taken prior to any cleaning of the skin surface that routinely took place before debridement of wound tissue. Chronic wounds were cleaned by applying gauze soaked with Prontosan wound irrigation solution (B. Braun Medical, Sheffield, UK) for ten minutes prior to sharp debridement of tissue from the top of the wound (wound debridement, WD). Wound debridement samples were only taken where debridement was deemed to be necessary for the standard wound care. Wound swabs were taken after irrigation of the wound with Prontosan to remove loose tissue, using a dry swab and the Z swab method (wound swab, WS). The Z swab method was the routine method used in the clinic at the time of sampling.

Control subjects (n = 8) (Table 1) were recruited from Sydney, Australia. The criteria for inclusion were not to have been diagnosed as diabetic, between 50–80 years of age, and without the use of antibiotics within the previous four weeks. Skin swabs were collected from the left and right feet of control subjects as described above. Samples were taken from all participants every two weeks for a 10-week period (6 time points in total). All samples were processed for DNA on the day of collection, or stored at 4°C until processing the next day. These storage conditions have been shown to adequately preserve the microbial profile of skin swab samples (Lauber et al., 2010).

Extraction of microbial DNA from skin and wound swabs and wound debridement tissue

Genomic DNA was extracted from all skin and wound samples using the BioStic DNA extraction kit (MO BIO Laboratories, Carlsbad, CA, USA). Swab heads were cut off the plastic applicator using sterile surgical scissors into the bead beating tube from the DNA extraction kit, before addition of buffer CB1. For wound debridement tissue, the tissue was directly placed into the bead beating tube. All subsequent steps were in accordance with the manufacturer’s instructions, and DNA was eluted in 50 µl of solution CB5 (10 mM Tris pH 8). The extracted DNA was quantified on a Qubit^® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) with a Qubit^® dsDNA HS Assay Kit (Life Technologies, Carlsbad, CA, USA).

Preparation of 16S rRNA gene libraries for Illumina sequencing

A library of the V4 region of the 16S rRNA gene was prepared for Illumina sequencing from the isolated microbial DNA samples. Samples were amplified using primers based on the Caporaso et al. design (Caporaso et al., 2012), which were modified to include eight nt rather than 12 nt barcodes, and include a barcode on both the forward and reverse primer (V4_forward and V4_reverse; Table 2). Different barcoded primers were used for each sample. For skin samples, the V4 region was amplified from 500 pg template DNA; for wound samples template DNA started at 10 ng, but in some cases up to 50 ng was used where a PCR product was not obtained with lower amounts of template DNA. Each sample was subjected to 10 cycles of PCR with 0.5 µM each of V4_forward and V4_reverse barcoded primers in a 50 µl PCR reaction that contained 1 × Taq core PCR buffer (Qiagen, Venlo, Netherlands), 1 × Q solution, 250 µM dNTPs, and 1.25 U Taq DNA polymerase. Thermal cycling was carried out at 95°C for two minutes, followed by 10 cycles of 95°C for 15 s, 50°C for 30 s and 72°C for 90 s, followed by a final extension at 72°C for five minutes. Excess primer was removed via a magnetic bead clean-up using 0.8 volume of Axygen^® AxyPrep Mag beads (Corning, NY, USA) and the eluted amplicons were subjected to a further 20 cycles of PCR with 0.25 µM enrichment primers (Illumina_E_1 and Illumina_E_2; Table 2). The PCR reaction and cycling was performed as described above, except that the annealing temperature was increased to 55°C and 20 thermal cycles were performed. Following confirmation of the PCR product on a 1% agarose gel, the amplicons were purified using Axygen^® AxyPrep Mag beads (Corning, NY, USA) and quantified on a Qubit^® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) with a Qubit^® dsDNA HS Assay Kit (Life Technologies, Carlsbad, CA, USA). Equimolar (2 ng) amounts of the 16S amplicons obtained for each skin and wound sample were then pooled and the molarity of the pooled amplicons determined using a Bioanalyser High Sensitivity DNA chip (Agilent Technologies, Santa Clara, CA, USA).

Illumina sequencing and data analysis

The PCR amplicons from 264 samples (including positive and negative controls) were sequenced over two separate runs on an Illumina Miseq using 500 cycle V2 kits. Sequences were demultiplexed using phylosift (Darling et al., 2014) and read pairs merged using FLASH (Magoc & Salzberg, 2011). Sequences were quality filtered and processed into OTUs using USEARCH v 1.8.1 (Edgar, 2010) (fastq_filter command with the fastq_maxee option set to ‘2’ to remove all sequences with two or more expected errors). Further quality filtering and operational taxonomic unit (OTU) clustering was carried out in QIIME (Caporaso et al., 2010b) version 1.9.0. The split_libraries.py command was used with the –l and –L options set to 240 and 260 respectively, to remove sequences shorter than 240 and longer than 260 base pairs. Sequences were clustered into OTUs at 97% similarity using the pick_open_reference_otus.py script using default settings except that singleton OTUs were removed, and the usearch61 method was used for chimera filtering.

Taxonomy was assigned to OTUs (assign_taxonomy.py) using the UCLUST method (Edgar, 2010) against the Greengenes (DeSantis et al., 2006) database pre-clustered at 97% similarity, accessed from the QIIME website (ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_8_otus.tar.gz). Representative sequences from each OTU were aligned against the Greengenes alignment using Pynast (Caporaso et al., 2010a) (align_seqs.py), OTUs which failed alignment were filtered from the final OTU table (filter_otus_from_otu_table.py). A phylogenetic tree was built from aligned representative OTU sequences (make_phylogeny.py script) using Fasttree2 (Price, Dehal & Arkin, 2010), with the –r option set to midpoint for tree rooting.

Diabetic skin samples adjacent to wounds were found to be more similar to wound than contralateral skin samples (see Fig. S1), and were removed so as not to confound comparisons between diabetic and non-diabetic skin. To ensure more even sample sizes between the diabetic and non-diabetic groups, only the right foot samples were included from the non-diabetic group for all downstream analyses. Alpha diversity was calculated using Phyloseq (McMurdie & Holmes, 2013) for the observed number of OTUs, Chao 1 and Shannon diversity indices on data rarefied to 30,000 sequences per sample. Significance testing was carried out on alpha diversity estimates using the Wilcoxon rank sum test in R.

Initial beta-diversity analysis was carried out in QIIME on a rarefied OTU table (30K sequences per sample) using the weighted unifrac metric, and the generate_boxplots.py script used to compare unifrac distances between groups of samples. Futher beta diversity analyses, were carried out in Phyloseq, using weighted unifrac distances calculated from an OTU table with raw counts subject to variance stabilising transformation implemented in DEseq2 (Love, Huber & Anders, 2014) as described here (McMurdie & Holmes, 2014). Weighted unifrac distances matrices were also subject to principal coordinates analysis using the Phyloseq package, and significant differences in variance between groups (diabetic and control skin) were determined with PERMANOVA (adonis function) implemented in the Vegan package (Oksanen et al., 2015) in R, using a nested model formula (∼health/subject + subject) and the weighted unifrac distance matrix.

The Wald test for differential abundance was used as implemented in the DESeq2 package in R. Multivariate correlation analysis was carried out against OTUs and wound duration and area using Pearson scores with Bonferroni correction, and p-values were determined via bootstrapping with 100 permutations (implemented in QIIME using the observation_metatdata_correlation.py command). OTU tables were filtered to remove OTUs present in less than 10% of samples for both differential abundance and correlation tests.

A Random forest learning algorithm implemented in R (Liaw & Wiener, 2002) was used to determine if diabetic status could be predicted from the foot skin microbiome. Skin samples were randomly divided into two equal subsets (restricting samples from the same participant to the same subset) for training and testing of learning algorithms. The variance stabilizing transformed OTU table was filtered to include skin samples only, and to remove OTUs observed in less than 10% of samples, and used as the input matrix for the Random forest algorithm. The Random forest fitted on the training subset was created using bootstrapping of one third of the training samples with replacement. As a general practice the rest of the samples were used as a validation set in order to decrease the risk of over-fitting associated with classification algorithms. An optimisation to minimise the out of bag error (classification error on validation data) was used to obtain the optimal number of taxonomic units accessed at each iteration of decision tree creation. Two hundred decision trees consisting of 30 OTUs evaluated at each node of the tree were created. The Random forest model was then used to predict the health status of the subjects in the test subset.

Analysis of the stability of skin microbial communities over time was carried out by comparing intrapersonal weighted unifrac distances between the diabetic and control skin samples, along with intrapersonal distances for all samples. Kruskal–Wallis tests were used to determine significant differences between groups.

Pearson’s Product Moment Correlation was used to test for correlations between wound size or duration and OTU abundance in wound samples as implemented in QIIME (observation_metadata_correlation.py). P-values were calculated using bootstrapping with 100 permutations, and Bonferroni correction for multiple testing. Kruskal-Wallis tests for OTUs that were differentially abundant in healing vs non-healing wounds were implemented in QIIME (group_significance.py). Wounds were classified as healing or non-healing based on a reduction in wound area since the last sampling time (healing) or no change or greater wound size area since the last sampling (non-healing). OTU tables were filtered to remove OTUs present in less than 10% of samples prior to testing.

Inter-visit weighted unifrac distances were compared to the overall degree of healing (1 − (final wound area/initial wound area)) using the lm function of the stats package in R.

Quality filtered sequence data has been deposited in the European Nucleotide Archive under study accession number PRJEB17696. A script containing the code used to process the data in R is provided as supplementary data, along with all the necessary input files, including OTU table and phylogenetic tree.

Cohort characteristics

The diabetic cohort (n = 8) consisted of five males and three females, with an average age of 68.9 ± 8.2 (range 58–81), average BMI of 35.4 ± 5.9 (range 27.2–47.1), and all had at least one foot ulcer which had been present for a average time of 9.1 ± 8.4 months (range 1.5–24 months). All wounds were neuropathic, with the exception of Patient 6 where the wound was ischemic. Two of the eight wounds healed during the course of sampling. Wounds were dressed with either Allevyn foam (Smith and Nephew) to promote moist wound healing, Zetuvit dressing (Hartmann) to remove excess wound exudate, Inadine antimicrobial dressing (10% povidone-iodine) (Johnson and Johnson), or Acticoat flex (antimicrobial silver coated) (Smith and Nephew), as deemed appropriate by the treating podiatrist or wound care nurse. All wounds were located on the plantar aspect of the foot. Details of the specific location of each wound, along with size and treatment over time and are provided in Table S3.

The control cohort (n = 8) consisted of 2 males and 6 females, with an average age of 62.8 ± 13.4 (range 50–81), average BMI of 28.0 ± 6.6 (range 20.4–37.9), and did not have wounds present on the feet.

Sample processing, 16S PCR and sequencing

A total of 242 samples were collected from the diabetic and control cohorts, including 170 skin swabs (85 diabetic and 85 control), 40 wound swabs and 32 wound debridement samples. Full details for samples collected at each time point for each participant can be found in Table S1 (diabetic participants) and Table S2 (control participants).

DNA yields obtained from diabetic skin swabs varied from 0.51 to 600 ng, with a median of 8.5 ng. Three skin samples did not yield enough DNA to be measured by the Qubit assay, however 16S rRNA gene PCR products were still obtained. Control skin sample DNA yields ranged from 0.5 to 41.7 ng (median 5.55), with 20 samples falling below the detection limit of the Qubit assay (<5 pg/µl). Of these 20 samples, PCR products were obtained for all but 3. DNA yields from wound swab samples ranged from 15 ng to 5.6 µg (median 760 ng) and for wound debridement samples ranged from 170 ng to 5.8 µg (median 1.2 µg). One wound swab sample did not yield enough DNA to be detected. Negative control swabs (n =4) did not yield enough DNA to be detected, and also did not yield detectable PCR products.

PCR products from the V4 region of the 16S rRNA gene were obtained for 257 of the 273 samples collected. Repeated attempts were made with increased amounts of template for those samples that did not initially yield a PCR product, however no PCR product was obtained (detailed in Tables S1 and S2). Amplicons from the remaining 257 samples were pooled and paired-end sequenced over two separate MiSeq runs with V2-500 cycle kits. Sample from four diabetic and four control subjects were sequenced in each run (Table S4). A median coverage of 73,599 sequences per sample was obtained (minimum 1,683, maximum 297,817). Negative controls (two blank swab and 2 no DNA PCR controls) had between 1,508 and 27,840 sequences assigned. The final sequencing coverage obtained for each sample can be found in Table S4 . Because negative control samples contained taxa that are similar to those found on skin (e.g., Staphylococcus, Corynebacterium and Acinetobacter) specific taxa were not removed from the data, rather samples with less than 30,000 sequences (n = 5) were removed from the analysis, based on the highest level of sequencing reads obtained from negative controls. A PERMANOVA test was run on a weighted unfrac distance matrix generated from variance stabilising transformed counts to assess the amount of variance attributable to the two different sequencing runs, (run + subject). Sequencing run was a significant factor accounting for 3.0% of the variance (p < 0.001), while inter-individual differences accounted for 34.5% (p < 0.001).

The microbiome of diabetic skin is less diverse than control skin

Diversity in all three groups was significantly different for observed richness, Chao1 and Shannon diversity indices (likelihood ratio test, p < 0.01). Diabetic skin was significantly less diverse than control skin for richness and Chao1 indices (Wilcoxon rank sum test, p < 0.01) (Fig. 1). Control skin had a median of 998.5 observed OTUs, compared to 435 for diabetic skin. Wounds were also significantly less diverse than diabetic skin with a median of 145 observed OTUs.

The skin microbiome is significantly different between diabetic and control subjects

Skin microbial communities overall were significantly different between diabetic and control skin (Fig. 2). A clear distinction can be observed between the sample types, and this was confirmed by a PERMANOVA test (∼health/subject + health), where health (diabetes vs control) was a significant factor accounting for 11.7% of the variance (R² = 0.117, p = 0.001). Subject (inter-individual differences) was the most significant factor accounting for 34.6% of the observed variance (R² = 0.346, p = 0.001).

Abundant taxa from skin are similar between persons with diabetes and healthy controls

Despite the clear distinction between diabetic and control skin in the PCoA plot above, the most abundant taxa from both groups were similar. Foot skin communities from diabetic skin were dominated by the genera Staphylococcus, followed by Acinetobacter and Corynebacterium, then unclassified Enterobacteriacea. Control skin was dominated by the genera Staphylococcus, followed by Acinetobacter, Kocuria, Corynebacterium and Micrococcus, (Fig. 3).

To determine which OTUs were contributing to the significant difference detected in the PERMANOVA analysis, the Wald test as implemented in the DESeq2 package (Love, Huber & Anders, 2014), was carried out. Sixty-nine OTUs were identified as significantly different in abundance (adjusted p < 0.05), all with an average abundance of less than 1%. A full list of the results can be found in Table S5. Similar results were found when re-running the analysis at the Genera level, with 24 genera identified as significantly different, but all at an average relative abundance of less then 1% (Table S6).

The foot skin microbiome may predict diabetic status

Despite only low abundance OTUs showing significant differences between diabetic and non-diabetic skin, a Random Forrest classifier was able to predict diabetic status from the foot skin microbiome. The model achieved an overall accuracy of 85.0%, with a sensitivity of 79.2%, and specificity of 93.8%. The negative predictive value (75.0%) was lower than the positive predictive value (95.0%). The classifier’s Gini index provided a list of 106 OTUs that were important in the classification task (Table S7); the majority were low abundance OTUs (103 OTUs <  1% average relative abundance), and the majority of these were more abundant in control than diabetic skin (75 OTUs).

Stability of the diabetic skin microbiome over time

Longitudinal analysis of the skin microbiome over time showed a trend of lower stability for diabetic skin than non-diabetic skin (Fig. 4), however this difference did not reach significance (p = 0.09), while both control and diabetic skin intrapersonal differences over time were significantly smaller (i.e., more stable) than inter-individual differences (p < 0.05).

Microbiota of chronic diabetic wounds overlap with skin and differ between patient

Wound swab and debridement samples were similar in taxonomic composition, and the top ten OTUs from all wounds per patient are shown in Fig. 5. The most abundant OTU detected in wounds was also the most abundant OTU found on skin, Staphylococcus sp. (OTU 1084865), and was present in the wounds of all eight patients. Other skin associated OTUs found in wounds included Corynebacterium (OTU 1011712), which was in the top 10 OTUs in six out of eight patient’s wounds.

The Wald test for differential abundance between diabetic skin and wounds identified four OTUs that were significantly more abundant across all wounds (two classified as Enterobactericaeae, one as Serratia and one as Finegoldia). The complete list of results can be found in Table S8.

The top 10 OTUs in wounds per patient over time are shown in Fig. 6. Of the eight wounds, six are dominated by the most abundant skin OTU, at the majority of time points measured (Staphylococcus OTU 1084865). Only Patients 6 and 10 showed wound profiles dominated by non-skin associated taxa across the time period surveyed. No significant correlations were found between any abundant OTUs (average abundance >  1%) and wound duration or healing status. No significant correlation was found between the overall degree of wound healing, and inter-visit weighted unifrac distances in individual wounds (Fig. S2, p = 0.29). However, some interesting observations were made that correlated to clinical events. For example, the wound of Patient 6 had been present for 24 months at the start of the study. It was dominated by Enterobacteriacaea and showed little healing until time point 3, which coincided with an angioplasty procedure to improve blood flow to the foot. This was followed by resolution of the wound within two weeks. When Patient 7 presented to the clinic, the wound had been present for 12 months, and was dominated by an OTU from the Neisseriaceae family. Following the standard treatment of debridement and wound dressing, rapid healing was observed, as well as a shift to a community dominated by the most abundant skin OTU.