Changes in Foot Skin Microbiome of Patients with Diabetes Mellitus Using High-Throughput 16S rRNA Gene Sequencing: A Case Control Study from a Single Center.

BACKGROUND Worldwide, the treatment of complications associated with type 2 diabetes mellitus, including diabetic foot ulcer (DFU), results in an economic burden for patients and healthcare systems. This study aimed to use high-throughput 16S rRNA gene sequencing to investigate the changes in foot skin microbiome of patients with diabetes mellitus from a single center in China. MATERIAL AND METHODS Fifty-two participants were divided into 4 study groups: healthy controls (n=13); patients with short-term diabetes (<2 years; n=13); patients with intermediate-term diabetes (5-8 years; n=13); and patients with long-term diabetes (>10 years; n=13). Swabs were analyzed from the intact skin of the foot arch using high-throughput 16S ribosomal RNA sequencing. RESULTS Microbiome phylogenic diversity varied significantly between the study groups (whole tree, P<0.01; Chao1, P<0.01), but were similar within the same group. The findings were supported by non-parametric multidimensional scaling (stress=0.12) and principal component analysis (principal component 1, 8.38%; principal component 2, 5.28%). In patients with diabetes mellitus, the dominant skin microbial phyla were Firmicutes, Proteobacteria, Actinobacteria, and Bacteroidetes. CONCLUSIONS High-throughput 16S rRNA gene sequencing showed dynamic changes in the skin microbiome from the foot during the progression of diabetes mellitus. These findings support the importance of understanding the role of the skin microbiota in the pathogenesis of DFU.

Background

Diabetic foot ulcer (DFU) is a serious and highly disabling complication of diabetes. In China, the prevalence of the total diabetes and prediabetes is 9.7% and 15.5%, respectively, accounting for 92.4 million people with diabetes and 148.2 million people with prediabetes [1]. About 15–25% of diabetic patients develop DFU, 21.5% of which result in amputation [2]. DFU formation is a complex process involving multiple factors such as peripheral neuropathy, peripheral vascular diseases, foot deformity, soft tissue lesion, and persistent local pressure [3-5]. The large population, along with significant changes in diet and labor delivery modes have resulted in a larger diabetic population and increased diabetic foot problems in China. DFUs are persistent, hard to heal, and easily infected, making them the major cause of lower extremity non-traumatic amputation worldwide [6]. Physical limitations, and the economic and psychological burdens limit the participation of patients with DFU in social activities.

Infection is a major cause of persistent ulcers and non-healing wounds, which is related to host immune system dysfunction, wound condition, and the type of infectious bacteria [7,8]. Undoubtedly, microbes play an important role in causing infections [9]. Some potentially pathogenic bacteria are part of the human microbiome and can aggravate the infection by forming biofilms that inhibit the actions of antibiotics [10-12]. Moreover, microorganisms can regulate our innate and adaptive immune system. Increasing evidence suggests that dysbiosis of the skin microbiome, in the absence of an invading pathogen, leads to skin disease. The observed dysbiosis is probably the result of an adaptive response of the microbiome to several disease-induced body changes, such as changes in pH, immune system, and skin barrier integrity.

In this study, we aimed to determine whether there is an observable dynamic change in the microbiome during the course of diabetes. This study is part of our laboratory’s continuing efforts to investigate microbiota changes to further assess the diabetic state of patients. In patients with diabetes, assessing the physical condition as well as the risk of related complications based only on disease duration is insufficient. It is unclear when these complications are likely to occur. This cannot be clarified by examining only the patient’s blood glucose level or duration of the illness. Furthermore, there has no objective laboratory tests, which may be more helpful by offering some rigorous evaluation results. The patient’s original physical condition, underlying diseases, genetic factors, lifestyle, gender, and body habitats can affect the occurrence of diabetes-related complications and the risk of amputation. Our study focusing on foot skin microbiome might provide a theoretical foundation to help us comprehensively understand the development of diabetic complications and the relation between microbiota and diabetes.

Material and Methods

Ethics statement

The study protocol was approved by the Ethics Committee of Guangzhou General Hospital of Southern Theater Command, PLA (approval no. [2018] 021). All swabs were collected in the outpatient department of General Hospital of Southern Theater Command of China. Written consent was obtaining from all patients for sample collection and subsequent analyses.

Participant selection

Adult participants were divided into the following 4 groups (n=52, 13 in each group): group A, control group (healthy adults); group B, short-term group (diagnosed with diabetes for less than 2 years); group C, middle-term group (diagnosed with diabetes 5–8 years ago); and group D, long-term group (diagnosed with diabetes more than 10 years ago). The inclusion and exclusion criteria are listed in Table 1.

Table 1

The inclusion and exclusion criteria used in the present study.

The inclusion criteria*	Age >40 years old No history/family history of diabetes No history/family history of autoimmune diseases No foot disease (psoriasis, dermatophytosis, etc.) No antibiotic therapy within the last three months
The inclusion criteria**	Age > 40 years old, with type 2 diabetes mellitus, without complications The duration of illness met the requirement of each group No history/family history of autoimmune diseases No foot disease (psoriasis, dermatophytosis, etc.) No antibiotic therapy within the last three months
The exclusion criteria	An unstable blood glucose level, a request for hospital treatment in the last six months, or the patient has any complications Severe cardiovascular diseases Metabolic diseases Immune diseases A previous history of radiotherapy or chemotherapy

The inclusion criteria applied to group A;

the inclusion criteria applied to groups B, C, and D.

The exclusion criteria applied to groups A, B, C, and D.

Swab sample collection

Sterile rayon tipped swabs were pre-moistened with a sterile solution containing 0.15 M NaCl and 0.1% Tween 20. All swabs samples were collected from the intact skin at the left foot arch. Skin swabs were collected by firmly rubbing the moistened swab over the foot skin surface over 100 times, for a period of 30 seconds. All swabs were placed in sterile tubes after sampling and then stored in liquid nitrogen within 60 minutes.

DNA isolation and purification

DNA was extracted from all swab samples using a DNA extraction kit (MOBIO Laboratories, Carlsbad, USA). DNA concentration and purity were then measured (Thermo Fisher Scientific, USA).

16S amplicon sequencing (16S AS) and sequence analysis

The V4 region of the 16S rRNA of the DNA samples was amplified using the following primers:

515F: 5′-GTGCCAGCMGCCGCGGTAA-3′;

806R: 5′-GGACTACHVGGGTWTCTAAT-3′.

The primers were synthesized by Invitrogen Corporation (Carlsbad, CA, USA). Polymerase chain reaction (PCR) was performed on an S1000 thermal cycler (Bio-Rad, Hercules, USA) with 2×Premix Taq (25 μL; Takara Biotechnology, China), each primer (1 μL), and DNA template (3 μL). PCR thermocycling conditions were as follows: 94°C for initialization (5 minutes); 94°C for denaturation (30 cycles of 30 seconds), 52°C for annealing (30 seconds), and 72°C for extension (30 seconds); followed by a final elongation at 72°C (10 minutes).

Detection, pooling, purification, and sequencing

Agarose gel electrophoresis (1%) was used to detect the lengths and concentrations of the PCR products. The bright-main-band samples were used for further experiments. The PCR products were mixed, and then purified (EZNA Gel Extraction Kit; Omega, USA). Sequencing libraries were constructed using the NEBNext® Ultra™ DNA Library Prep Kit (New England Biolabs, USA). The library quality was assessed with the Qubit@2.0 Fluorometer (Thermo Fisher Scientific) and the Agilent 2100 Bioanalyzer. The libraries were then sequenced using the Illumina HiSeq platform (USA).

Data analysis

Descriptive statistics [mean±standard deviation (x±SD)] were used to describe the data. Group differences were determined using one-way analysis of variance (ANOVA) with SPSS 19.0 software. The means of 2 groups were compared with the Dunnett T3 test. A value of P<0.05 was considered statistically significant. The paired-end raw reads were filtered to obtain high-quality clean reads. The barcodes and primers were then removed, leaving the paired-end clean reads. If more than 10 of the reads overlapped with the same DNA fragment, the raw clean reads were merged using FLASH v1.2.11. The effective clean tags were obtained by filtering raw tags. Next, operational taxonomic unit (OTU) clustering was performed. Sequences with similarity more than 97% were assigned to the same OTU. The singleton OTUs were removed. The chimera sequences were detected and removed (UCHIME de novo algorithm). The GreenGene database and Qiime were used to annotate taxonomic information. For further analysis, the number of effective tags was calculated and an OTU taxonomy synthesis information table was generated. To study the phylogenetic relationships among different OTUs, KRONA software (KronaTools-2.7) was used. GraPhlAn software (version 0.9.7) was used to study the genera composition and abundance. In order to study the differences in dominant species among different samples, an OTU representative sequence was selected (relative abundance >50). Multiple sequence alignments were conducted (FastTree 2.1), and information regarding relative abundance and species annotation was displayed. Abundance information of OTUs was normalized using a sequence-quantity standard. Subsequent analyses, such as alpha and beta diversity analyses were performed based on output data.

Results

Cohort characteristics

From September 2018 to November 2018, 64 participants were enrolled in the study. Twelve participants were excluded (9 patients eventually refused to participate; 3 patients were suspected to have had diabetes long before diagnosis; and 1 patient provided a contaminated sample). Detailed information about the final 52 participants in the 4 groups is shown in Table 2. The patients in group A were significantly younger than those in group D (P<0.01). The recent serum glucose levels (mmol/L) were significantly different between group A and group C (P<0.01). In addition, the disease duration was significantly different between every 2 groups (P<0.01). No significant difference in height, weight, and body mass index (BMI) was observed among the 4 groups (P>0.05) (Table 2).

Table 2

Detailed information of the patients in the 4 groups.

ITEM	Group A	Group B	Group C	Group D	F
Disease Duration	–	1.1±0.7	6.2±0.9	15±3.9	167.674**
Male/Female	7/6	7/6	7/6	6/7
Age	52.5±6.1	51.7±8.8	60.8±9.5	63.7±12.2	5.254**
Height, cm	163.2±7.9	167±5.4	167.1±9.2	161±8.0	1.834
Weight, kg	63.5±11.1	72.5±6.5	68.5±12.1	68.1±11.9	1.549
BMI	23.7±3.0	26.1±1.6	24.7±4.4	26.3±4.2	1.615
Recent serum glucose level, (mmol/L)	6.1±0.7	6.4±0.5	7.0±0.7	7.7±2.2	4.644**

Data are expressed as no. of subjects or mean±SD. SD – standard deviation; BMI – body mass index (weight, kg/(height, m×height, m)).

Dunn’s test, P<0.01.

Microbiome diversity and richness increased during diabetes development

The alpha diversity index gradually increased from group A to D (Figure 1). The diabetic groups exhibited higher diversity indices than the control group. Group D had a higher diversity index than group B for phylogenetic diversity (PD) whole tree (P<0.01, statistically significant) and the Shannon index (P>0.2, not statistically significant). The richness measure (Chao1) showed the same distribution characteristics among the 4 groups (P<0.01), along with observed OTUs (P<0.01).

Figure 1

Box plots of the alpha diversity indices of the 4 groups. (A) The observed species number showed that the numbers of genera in the diabetic groups were higher than those in the control group. The group D has a higher number of species than the group B. (B) The Chao1 estimator values were significantly different among the 4 groups (P<0.01). (C) The phylogenetic diversity (PD) whole tree index among the 4 groups (P<0.01, statistically significant) were based on operational taxonomic units (OTUs) clustered at 97% similarity for the control skin group. (D) The Shannon index showed an increasing trend from group A to group D; however, the increase was not statistically significant. A, control group; B, the short-term diabetes group; C, the middle-term diabetes group; D, the long-term diabetes group, patients with no complications.

Foot skin microbiome composition of the 4 groups was different

The compositions of the foot skin microbiome among the 4 groups were different (P<0.05). Our results showed that the composition of the samples from the same group was similar, with obviously closer evolutionary relationships (Figure 2A), especially in groups C and D. The similarity between microbiome components was higher for samples from the same group (the distance between points represents the similarity between them), but lower for samples from different groups (Figure 2B, 2C).

Figure 2

(A) Analysis of all samples by unweighted pair-group method with arithmetic mean (UPGMA). The unweighted UniFrac tree displays species classification and evolutionary relationships based on unweighted UniFrac distances. The structure of the clustering tree and the distance between the branches are shown in the left panel. The relative abundance of each sample at the phylum level is shown in the right panel. We estimated the differences between the 4 groups [control skin group (A), the short-term diabetes group (B), the middle-term diabetes group (C), and the long-term diabetes group (D)] using the UPGMA method with average linkage to interpret the distance matrix. The results showed that samples from the same group displayed obvious aggregation, indicating differences in microbiome components between different groups. (B) Non-metric multidimensional scaling (NMDS) analysis was performed. NMDS ordination is based on the unweighted UniFrac distance matrix between all samples. Each point represents a specific sample. The samples in a group are represented by a specific icon. The closer the elements are to each other in these dimensions, the more similarity they exhibit (stress=0.12). (C) Unweighted UniFrac principal coordinate analysis (PCoA) was performed. Principal component (PC)1 and PC2 in the x- and y-axis, respectively, represent two principal discrepancy components. The contribution value to the discrepancies by the component is shown in brackets (percentage). Each point represents an individual element. Each icon represents a specific group. (PC1=8.38%, PC2=5.28%).

Dynamic microbiome changed over the course of diabetes

According to the aforementioned results, species abundance and diversity showed an increasing trend from group A to D. The components of individual samples from the same group showed significantly closer evolutionary relationships and were similar to each other.

On investigating the changes in abundance (Figure 3A, 3B) at the phylum level, results showed that during the course of the disease, including disease development, disease progression, and acquiring stable disease state without complication, certain microorganisms, such as Actinobacteria and Bacteroidetes showed significant dynamic changes. Despite this, the highly abundant taxa among the 4 groups were similar. The species that dominated the foot skin was Firmicutes, followed by Proteobacteria, Actinobacteria, and Bacteroidetes. The relative abundance of bacteria was analyzed at the phylum, class, order, family, and genus levels. The variations in the abundance of representative organisms in each group are summarized. Heatmaps (Figure 4) were obtained showing the relative abundance of different phyla. The heatmap shows an overview of the numeric difference. The control group (group A) was dominated by OP3, Chlorobi, and Firmicutes. The short-term group (group B) was dominated by Parvarchaeota, Chloroflexi, Euryarchaeota, Gemmatimonadetes, OP11, OD1, Elusimicrobia, TM6, Planctomycetes, OP3, and Firmicutes. The middle-term group (group C) was dominated by Actinobacteria, Armatimonadetes, Verrucomicrobia, Cyanobacteria, Thermi, and Tenericutes. The long-term group (group D) was dominated by Nitrospirae, Fusobacteria, Thermotogae, Aquificae, SR1, Chlamydiae, Spirochaetes, WPS-2, Acidobacteria, Crenarchaeota, and Proteobacteria. The visual graph showed that at certain time points during the disease, certain microorganisms with closer evolutionary relationships appeared to have a growth advantage; however, this did not significantly alter the microbiome composition during disease progression.

Figure 3

Changes in abundance. (A) Relative abundance of species at the phylum level in individual samples. (B) Relative abundance of species at the phylum level in each group [control skin group (A), the short-term diabetes group (B), the middle-term diabetes group (C), and the long-term diabetes group (D)].

Figure 4

Heat map showing the relative abundance of different phyla in different groups, at the phylum level. Species annotations are listed longitudinally, and group information is shown transversally. The left tree is a species-related clustering tree, and the upper tree is a group-related clustering tree. The heatmap was constructed based on the differences in relative abundance of species between the different groups, with colors gradually changing from deep red to deep blue, corresponding to high and low relative abundance, respectively.

Characteristics of foot skin microbiome in normal and diabetic individuals

The microbiome richness, diversity, and composition showed significant alterations. The abundance of several microorganisms was markedly different among the different groups; therefore, they might have a strong influence on the composition of each group. LEfSe was used to identify the different taxa among the groups (Figure 5). A cladogram was generated to represent the evolutionary position of every taxa in disease groups (from B to D) from the phylum level to the genus level, and show the most differentially abundant taxa between groups. Small circles shaded with different colors in the diagram represent abundances of those taxa in the respective groups (red, B; green, C; blue, D). The diameter is proportional to the abundance. Taxa that have no significant difference in abundance between groups are uniformly colored in yellow. Linear discriminant analysis (LDA) was performed. Only taxa meeting the LDA threshold of >1.0 and P<0.05 are shown. No taxa in group A reached the minimum LDA score. A total of 77 differential genera in group D exhibited an absolute LDA score >1.0. The taxa with high LDA scores included Xanthomonas, Axonopodis, Acidobacteria, Melaninogenica, Ellin6512, DA052, Pasteurellales, Pasteurellaceae, and Acidobacteriales. A total of 19 taxa in group C and 12 taxa in group B were identified to have an LDA score >1.0. In group B, the taxa with high LDA scores included GCA004, Aminovorans, Phyllobacteriaceae, Nasimurium, and Nitrosa. In group C, the main taxa with high LDA scores included Jeotgalicoccus, Isosphaeraceae, Delftia, Thermomicrobia, Alaskensis, Varibaculum, and Cloacibacterium. Overall, the short-term, middle-term, and long-term groups displayed significant changes in the foot skin microbial community structures, compared with the control group, and these changes were positively correlated with the duration of illness.

Figure 5

LEfSe. (A) A cladogram representing the evolutionary position of every taxa in the disease groups from the phylum level to the genus level [groups B to D: the short-term diabetes group (B), the middle-term diabetes group (C), and the long-term diabetes group (D); linear discriminant analysis (LDA) score >1.0]. (B) Histogram of the LDA scores for differentially abundant genera between groups. Only taxa meeting an LDA threshold of >1 and P<0.05 are shown.

Discussion

The high rate of DFU recurrence is a clinical challenge. Proper interventions, such as external application of Tangbi Waixi decoction to the lower extremities, custom therapeutic footwear, and off-loading in high-risk patients can reduce the risk of DFU and relieve clinical symptoms. Moreover, nursing care is a major financial burden associated with the treatment of DFUs. However, the cost of treating an ulcer is generally much lesser than the cost of amputation. It is strongly believed that early intervention, such as regular foot thermometry, foot evaluation, and foot care can reduce the incidence of DFUs, avoid amputations, or at least extend the number of “ulcer-free” days in patients [13-17]. However, in the case of a very large population with limited medical resources, it is unrealistic and expensive to implement preventive interventions for all patients. An early warning system is needed to identify patient groups that require priority intervention.

This study showed that there are dynamic changes in the skin microbiome during the onset, progression, and deterioration of diabetes. The microbiome compositions of different groups were different. However, the differences were not accounted for by changes in the most abundant members of the skin microbiota. The most abundant microbes in the foot skin included Firmicutes, Proteobacteria, Actinobacteria, and Bacteroidetes, and the changes in their abundance among the different groups were not significant. The abundance of Proteobacteria was slightly reduced at disease onset (group B) but increased as the disease progressed (long-term stable disease state; group D). Actinobacteria abundance showed the opposite change pattern. Further, our results showed that the diversity of low abundance microbes increased with disease progression, which is contrary to the findings of a previous study, which indicated that the diversity of microbiome in DFU skin was significantly lower than that in control skin [18]. In our study, the 3 diabetes groups showed higher microbial diversity than the control group, and the long-term diabetes group exhibited higher diversity than the short- and middle-term groups (Figure 1). Our results showed that samples from the same group had similar microbiome communities, with similar compositions and obviously closer evolutionary relationships. The relative abundance of each sample at the phylum level showed that the microbiome structure in samples from the same group showed obvious aggregation (Figure 2A; left panel). At the phylum level, changes in abundances were observed. The abundance of Actinobacteria showed an increasing trend at the beginning of the disease (from group A to group C), followed by a decreasing trend (from group C to group D). Similar trends were observed for Bacteroidetes (Figure 3). The relative abundance of different phyla showed numeric differences between several main microorganisms. We found that the microbiome of group A shared higher similarity with group C than groups B and D, as indicated by the intuitive information from the heatmap. Group B showed an increase in the relative abundance of Parvarchaeota, Chloroflexi, Euryarchaeota, Gemmatimonadetes, OP11, OD1, Elusimicrobia, TM6, Planctomycetes, OP3, and Firmicutes, while group D was dominated by Nitrospirae, Fusobacteria, Thermotogae, Aquificae, SR1, Chlamydiae, Spirochaetes, WPS-2, Acidobacteria, Crenarchaeota, and Proteobacteria. The relatively abundant microorganisms within the same group shared closer evolutionary relationships (Figure 4). Even though the kind of microorganism having numeric increasing in group B and D are different, they are all minority in health group (group A). The most abundant microorganisms in group A showed decreasing trends in the disease groups (B, C, and D), including Chlorobi, OP3, and Firmicutes. The clustering tree (Figure 4) demonstrated that different disease groups (B, C, and D) showed different dominant species (phylum level), which have significantly closer evolutionary relationships. In order to determine the main contributing microorganisms, LEfSe was performed to identify the most differentially abundant taxons and a cladogram was generated to determine the evolutionary relationships (phylum level) between different taxa. Highest numbers of taxa in group D were identified to have an LDA score >1.0, while no taxa in group A achieved this score (LDA score >1.0 and P<0.05). In groups C and B, the number of taxa that reached the minimum LDA score was similar, but lower than that in group D. These results partially support the notion that changes in microbial communities may play a positive role in the self-healing process.

It has been suggested that a reduction in the diversity and stability of diabetes skin microbiome is driven by various environmental changes, such as increased pH or exposure to antibiotics [19]. The microbiome has a self-regulatory function via regulation of the immune system as well as secretion of microbes and immune cells [20-26]. The markedly closer evolutionary relationships among the microbiome of the same group supported the notion that the microbial communities in the same group may need to alter their structure to adapt to the diabetes-associated skin environmental changes. If the changes in microbiome richness, diversity, and abundance according to the disease state show a noticeable pattern, it could contribute to the diagnosis and treatment of diabetes. It has been proven that microbial dysbiosis may lead to diseases without evident pathogenic invasion [27]. In certain body sites, increased microbiome diversity is associated with disease state [28]. However, changes in the microbiome may not be a causative factor of disease but might represent an adaptive survival response to the changing environment. Thus, if the disease determinants precede microbiome dysbiosis, the dysbiosis is probably an adaptive response to the changing conditions. Dysbiosis can be beneficial before the regulation is out of control.

In this study, we showed that the microbiome of the foot skin is not stable in patients with diabetes over time. The longer the disease duration, the higher was the microbiome richness, diversity, and abundance. In the past, removal of most microorganisms was considered necessary [29,30]. However, our results showed that changes in the microbiome had a positive effect on wound healing, which is contrary to previous findings, which showed that increased stability delayed healing [31,32]. A double-blind multicenter randomized controlled trial proposed that topical therapy should be used early for moderate to severe disease, whereas at later stages or for mild disease, antibiotics are unnecessary. Antibiotic therapy cannot always be the first treatment choice [31]. Hence, a significant alteration in the microbiome may not be as deleterious as previously thought. During the process of diabetes onset, development, and deterioration, not only does the body undergo constant changes, but the skin microbiome is also altered, which might be an adaptive adjustment. This finding may contribute to the diagnosis and treatment of diabetes.

However, the study was conducted at a single center, with a small number of patients in each group, and was conducted by clinicians who managed the patients; those retired and self-employed individuals may not have a strong awareness of physical examination and some patients may be diagnosed later than the actual date of disease occurrence. These factors may have introduced bias into the study. The possibility of seasonal variation in the skin microbiome and the effects of patient age and gender were not included in this study. This study was conducted to analyze the microbiome of the skin of the foot, due to the importance of DFU in this patient population, but comparative studies from the skin at other sites in patients with diabetes may identify important associations and could be undertaken in future. Also, the method used to investigate changes in the skin microbiome, high-throughput 16S small subunit ribosomal RNA (rRNA) gene sequencing, has several limitations that include the possibility of sequence similarity between some bacterial species, which may require alternative gene targets for species identification. Future studies may be conducted to compare the use of high-throughput 16S rRNA gene sequencing findings with other methods, including metagenomics.

Conclusions

Alterations in the composition of microbiome are observed during disease progression. While it is known that microorganisms can aggravate wound infection, resulting in a non-healing wound, it is important to highlight that the skin microbiome represents a major part of the skin and influences skin health. With further research, we can better understand non-healing wounds and other skin conditions associated with diabetes.