Analysis of the Composition and Functions of the Microbiome in Diabetic Foot Osteomyelitis Based on 16S rRNA and Metagenome Sequencing Technology

Diabetes-related lower-extremity complications are a large and growing contributor to the disability burden worldwide (1). Diabetic foot osteomyelitis (DFO) is the most severe stage of DF. More than 20% of patients with severe diabetic foot infection (DFI) and 50–60% of patients with moderate DFI develop DFO, and the minor amputation rate of patients with DFO is 59.4% (2). The 5-year mortality for patients undergoing minor amputations is 46.2%. Long-term survival was worse in patients who underwent a major amputation with a 5-year mortality of 65.6%, which exceeds the reported 5-year overall mortality rate of cancer by 31.0%. DFO is common, complicated, and costly. In 2017, the direct cost of diabetes care was $237 billion, and more than one-third of the direct cost was attributed to the lower extremities (3). In recent years, the importance of antibiotics in DFO treatment has gradually been recognized, but the application of antibiotics is limited due to the limitations and one-sidedness of microbial identification technology (4). Fast and comprehensive analysis of the microbiome in wounds is a prerequisite for reducing the amputation rate and mortality. Although culture methods have a long history and have been widely used, shortcomings remain (e.g., time consumption and low sensitivity) (5). These deficiencies have hindered a correct understanding of the microbiome and have made it difficult to meet the needs of clinical application and the development of modern biological research. Consequently, technological updates are urgently needed in the field of microbial research. High-throughput sequencing technology does not rely on culture and can be used to examine the microbiome and its functions in wounds quickly and comprehensively, providing a new direction for the development of microbial research. The most widely used high-throughput sequencing technology for identifying bacteria is 16S rRNA sequencing. Among its advantages, 16S rRNA sequencing does not require culture and is simple, fast, low-cost, and widely applied. However, 16S rRNA sequencing provides only partial bacterial genome information. Even after full optimization, it is impossible to obtain detailed information at the species level, and it is difficult to analyze the functions of microorganisms in detail. Compared with 16S rRNA sequencing, metagenome sequencing is a revolutionary method that can examine all DNA in specimens without the use of culture or PCR amplification. Metagenome sequencing can provide sufficient and accurate species information, thus compensating for the shortcomings of the above two methods (6). However, there have been no studies detecting microorganisms in infected bone specimens based on metagenome sequencing. We therefore conducted a prospective study with the goal of better observing the microbiome and its functions in wounds of DFO based on cultivation, 16S rRNA sequencing, and metagenome sequencing.

No study has analyzed the differences and similarities between the microorganisms in wounds of DFO and posttraumatic foot osteomyelitis (PFO). Over a 1-year period, we prospectively enrolled 28 consecutive patients aged over 18 years who presented with osteomyelitis at Nanfang Hospital from September 2017 to September 2018. The patients were divided into two groups, namely, the DFO group (17 patients) and the PFO group (11 patients), according to whether or not they had diabetes. Osteomyelitis was suspected based on the patient’s clinical manifestations, laboratory tests, and imaging examination: 1) clinical manifestations included bone exposure (especially if the area was >2 cm²) or probe-to-bone test was positive and the occurrence of erythema and hardening of the toes (sausage-like toes); 2) laboratory tests included erythrocyte sedimentation rate (ESR) >70 mm/h; C-reactive protein (CRP), procalcitonin (PCT), and white blood cell (WBC) levels were increased; and microorganisms could be cultivated from specimens; 3) radiograph examination included radiographs that showed loss of cortical bone, accompanied by bone erosion or demineralization, and focal trabecular morphology loss or loss of bone marrow radio permeability; and 4) MRI examination that showed T1-weighted low-focus signals, T2-weighted high-focus signals, and high bone marrow short-term inversion recovery sequence signals (4). Patients were included if osteomyelitis was diagnosed, microorganisms could be cultivated from specimens, and the infected bone was exposed. Furthermore, the patients were included if they exhibited no obvious skin lesions on the wound surface or surroundings, were in good health, and were able to withstand debridement. The patients and their families were informed and agreed to participate in the study. The patients were ≥18 years old. Patients were excluded if they had an immune system that could not tolerate debridement or had taken immunosuppressants within 3 months before admission. The study was designed and implemented in accordance with the Declaration of Helsinki (2013), approved by the ethics committee of Nanfang Hospital (NFEC-2017–013), and registered on ClinicalTrials.gov (NCT04240964). Informed consent was obtained from all patients.

After the wounds were rinsed with sterile saline and hydrogen peroxide solution, 17 bone specimens were collected in the clinic setting from patients with DFO who required debridement to manage their osteomyelitis after removal of necrotic tissue from the surface to avoid soft tissue or sinus tract contamination, which was not predictive of the presence of pathogen with sufficient accuracy (7). After the wounds were rinsed with sterile saline and iodophor solution, 11 bone specimens were routinely collected in an operating room from patients with PFO. We obtained bone specimens and divided them into two parts under sterile conditions. One part was used for routine microbial culture, and the other was used for high-throughput sequencing.

Genomic DNA was extracted using a DNA extraction kit (YiRui, ShenZhen, China) according to the manufacturer’s instructions. The extracted DNA was quantified and quality controlled by 1% agarose gel electrophoresis (JS-power 300; PeiQin, Shanghai, China). Then, we amplified the V3-V4 variable region of the 16S rRNA gene for sequencing using forward and reverse fusion primers (341F: 5′-CCTAYGGGRBGCASCAG-3′ and 806R: 5′-GGACTACNNGGGTATCTAAT-3′). PCR was performed in a total volume of 20 μL containing 4 μL of 10× PCR buffer, 2 μL of 2.5 mmol/L deoxynucleotide triphosphates, 0.8 μL of forward primer (5 μmol/L), 0.8 μL of reverse primer (5 μmol/L), 0.4 μL of FastPfu polymerase, template DNA (10 μL) and H₂O (2 μL). PCR amplification was conducted under the following conditions: initial denaturation at 95°C for 2 min, followed by 27 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 45 s. A final extension was performed at 72°C for 10 min. Then, we purified and quantified the PCR products, followed by 16S rRNA sequencing on the Illumina sequencing platform with the PE250 sequencing protocol. The data were obtained in fastq format after sequencing. Then, we separated the sample data from the irregular sequencing data according to the barcode sequences. To obtain high-quality clean reads, raw reads were filtered according to the following rules: reads containing 10% unknown nucleotides and reads containing less than 80% bases with quality (Q-value) >20 were removed. The filtered reads were then assembled into raw tags according to overlap between paired-end reads with more than 10 base pair (bp) overlap and <2% mismatch. Clean tags were obtained after filtering out short, low-quality and chimeric raw tags. Then, we clustered raw tags into operational taxonomic units (OTUs). Taxonomy and abundance were assigned to OTUs by blasting against the Ribosome Database Project database.

We also selected a subset of specimens after MicroPITA analysis for metagenome sequencing. In brief, DNA fragments of the appropriate length (∼350 bp) were obtained by sonication. Then, the fragments were A-tailed and ligated to adapters. NovaSeq sequencing systems (Illumina) were used for sequencing and library validation. The raw data obtained by sequencing had a certain proportion of low-quality data, which were removed to obtain clean data. The remaining reads were mapped to the human genome, and the matching reads were removed as contaminants. The De Bruijn graph assembly method was used to splice clean reads by single-sample assembly and mixed assembly, and contigs were obtained. Contigs less than 1,000 bp in length were filtered out. This study was based on a basic local alignment search tool (BLAST) of the nonredundant contigs to the NCBI database to obtain taxonomic annotation information for the microbiome.

The next step was gene prediction. We first predicted the open reading frames of the assembled sequences. Then, we filtered out redundant genes to obtain nonredundant gene sets. Genes with similarity >95% were clustered together, and the longest gene in the same class after clustering was the representative gene. To obtain functional information, we compared the gene sequences were with five functional databases: Kyoto Encyclopedia of Genes and Genomes (KEGG), Nonsupervised Orthologous Groups (eggNOG), Carbohydrate-Active Enzymes (CAZy), Antibiotic Resistance Genes Database (ARDB), and Virulence Factors Database (VFDB).

Patient demographics and laboratory and clinical data were examined using χ² and t tests. Spearman correlation analysis was used to analyze the correlations between the dominant species and the clinical indicators of DFO patients. For all comparisons, the level of significance was set at 0.05.

The data generated and analyzed during the current study are available from the corresponding author upon reasonable request.

The main characteristics of the included patients are summarized in Table 1. Among these patients, 17 were diagnosed with DFO, and 11 patients were diagnosed with PFO. There were nine males and eight females in the DFO group, and the average age was 60.59 ± 9.39 years. The seven males and four females in the PFO group had an average age of 46.90 ± 13.62 years. There was no significant difference in sex composition (P = 0.580), the locations of the wounds (P = 0.157), and antibiotic use (P = 0.157) between the two groups. In terms of infection indicators, the WBC (P = 0.002) and neutrophil (N; P < 0.001), PCT (P = 0.003), CRP (P < 0.001), and ESR (P < 0.001) levels were significantly higher in the DFO group than in the PFO group (Table 1).

Microorganisms obtained from cultured infected bone specimens of DFO patients belonged to Proteobacteria and Firmicutes, corresponding to Gram-negative (G−) and Gram-positive (G+) bacteria, respectively. Twenty-three isolates from 10 genera were cultured from the infected bone specimens of DFO. The most frequently occurring genera were Streptococcus (4/23; 17.4%) and Enterococcus (4/23; 17.4%). Among the genera obtained by the culture method, G+ bacteria accounted for 47.8% (11/23), G− bacteria accounted for 52.2% (12/23), and there were no anaerobes. A total of 41.2% of the cultured infected bone specimens showed polymicrobial infections (Table 2). The 16S rRNA sequencing of infected bone specimens yielded five dominant phyla, including Proteobacteria and Firmicutes. A total of 242 genera were obtained, 18 of which were dominant: G− bacteria accounted for 61.1% (11/18), G+ bacteria accounted for 38.9% (7/18), aerobes accounted for 22.2% (4/18), facultative anaerobes accounted for 27.8% (5/18), obligate anaerobes accounted for 50.0% (8/18), and all samples showed polymicrobial infections (Table 3).

The 16S rRNA sequencing generated 1,034,689 counts, which were clustered at 97% similarity and indicated 7,555 unique OTUs. The Shannon and Simpson index values were significantly higher in the DFO group than in the PFO group (P < 0.001) (Fig. 1A and B). We found two obvious groups of microbes using Weighted UniFrac distance-based principal coordinate analysis, which showed that the microbiome of the DFO group was different from that of the PFO group (PERMANOVA [permutational multivariate analysis of variance] F value: 10.091; R²: 0.32457; P value <0.001) (Fig. 1C). A clustered heatmap was used according to the relative abundance profiles of the top 200 genera. The microbiota was more evenly distributed in the PFO group than in the DFO group (Fig. 1D). In total, 18 genera that had a relative abundance of >1% in samples from at least two subjects were authenticated. Prevotella was the most abundant genus in the DFO group, followed by Streptococcus, Anaerococcus, Staphylococcus, and Enterococcus. Halomonas was the most abundant genus in the PFO group, followed by Streptococcus, Bacteroides, Corynebacterium, Providencia, and Bradyrhizobium (Fig. 1E). The most abundant microorganisms were G− anaerobes in the DFO group and G− aerobes in the PFO group. To further characterize the impact of the DFO group and the PFO group on the circulating microbiome, a discriminative features cladogram and histogram based on an effect size cutoff of 3.6 for the logarithmic linear discriminant analysis (LDA) scores were plotted. Firmicutes (LDA score = 5.25; P = 0.001), Clostridia (LDA score = 5.03; P = 0.029), and Clostridiales (LDA score = 5.03; P = 0.029) were the top three identified microorganisms in the DFO group, followed by Prevotellaceae (LDA score = 5.02; P = 0.022) and Prevotella (LDA score = 5.02; P = 0.022). Oceanospirillales (LDA score = 5.46; P < 0.001), Halomonadaceae (LDA score = 5.46; P < 0.001), and Halomonas (LDA score = 5.46; P < 0.001) were the most widely identified microorganisms in the PFO group (Fig. 1F and G).

After filtering, a total of 1.2726 × 1,011 bp and 846,819,798 reads were obtained, totaling 127 GB of data (12.7 GB per sample). After assembly, 95,609 contigs were obtained, with an average length of 3,151.68 bp and N50 of 7,698 bp. In addition to the overview of the species in the DFO group and the PFO group, a heatmap containing the top 200 relative abundances of species was graphed. Compared with the species in the PFO group, the species in the DFO group showed a significantly high abundance (Fig. 2A). At the phylum level, the relative abundance of Firmicutes was significantly higher in the DFO group than in PFO group (t = 3.671; P = 0.001). The relative abundance of Proteobacteria was significantly higher in the PFO group than in the DFO group (t = −4.132; P < 0.001) (Fig. 2B). Prevotella was the most abundant genus in the DFO group, and its relative abundance was significantly higher than that in the PFO group (t = −3.817; P = 0.002). Halomonas was the most abundant genus in the PFO group, and its relative abundance was significantly higher than that in the DFO group (t = −5.074; P < 0.001) (Fig. 2C). Among the 22 dominant species, 6 belonged to Prevotella, the total relative abundance of which reached 13.6% in the DFO group. This genus had the highest proportion among all species in the DFO group. The most abundant species in the DFO group and the PFO group was Klebsiella pneumoniae, followed by Veillonella parvula. The third most abundant species in the DFO group was Prevotella intermedia, followed by Prevotella denticola, Thielavia terrestris, and Yersinia enterocolitica. The third most abundant species in the PFO group was Y. enterocolitica, followed by T. terrestris, Kluyveromyces marxianus, and P. denticola (Fig. 2D). On the basis of the LDA effect size, Prevotellaceae (LDA score = 5.49; P = 0.028) and P. intermedia (LDA score = 5.20; P = 0.009) were the two most widely identified microorganisms in the DFO group. No notable species were found in the PFO group (Fig. 2E and F).

P. intermedia (ρ = 0.31), P. denticola (ρ = 0.36), Prevotella jejuni (ρ = 0.46), and Prevotella fusca (ρ = 0.46) were all positively correlated with the duration of diabetic foot infection (DFI_d). K. pneumoniae was positively correlated with Wagner classification (ρ = 0.87) and negatively correlated with DFI_d (ρ = −0.46). B. fragilis was negatively correlated with WBC (ρ = −0.70), N (ρ = −0.70), and CRP (ρ = −0.70). Proteus vulgaris was positively correlated with WBC (ρ = 0.94), N (ρ = 0.92), and CRP (ρ = 0.68) (Fig. 3).

The microbiome of the DFO group had more functional genes than that of the PFO group (Fig. 4A). The most abundant virulence factor (VF) in the DFO group was Hsp60, followed by ClpC and ClpE (Fig. 4B). Antibiotic resistance (AR) in pathogens through the acquisition of resistance genes or mutations makes infections difficult to treat (8). The number of AR genes was much higher in the DFO group than in the PFO group (P = 0.030). The most abundant AR gene was the streptomycin resistance genes, followed by the lincosamide, macrolide, and tetracycline resistance genes (Fig. 4C). The KEGG database is a collection of large-scale molecule-level data sets (9). After blasting the filtered genes against the KEGG database, 4,626 KOs and 75,630 pathways were obtained. The most abundant pathway was signal transduction and glycolysis gluconeogenesis in the two groups (Fig. 4D). Evolutionary genealogy of genes: eggNOG is a public database of orthologous relationships (10). In the DFO group, replication- and repair-related functional genes were the most abundant, followed by translation (Fig. 4E). Genes involved in complex carbohydrate metabolism are listed in the Carbohydrate-Active Enzymes (CAZy) database (11). The abundance of CAZymes was higher in the DFO group than in the PFO group (Fig. 4F).

For over 150 years, clinicians have defined causative organisms based on the results of culture-based experiments (12). Almost every European and North American study reports Staphylococcus aureus as the most common microorganism cultured from DFO, followed by Staphylococcus epidermidis (12–15). These findings are not surprising; culture-based methods select only species that flourish under the typical nutritional conditions in microbiological laboratories, not the most abundant species. Therefore, Staphylococcus grows more easily than much-needed bacteria, such as anaerobes, especially in media that are commonly used in clinical microbiological laboratories within short periods of time (16). Aerobic culture was used in this study, and no anaerobes were cultured. In fact, because of the difficulty in maintaining the anaerobic environment, many clinical microbiology laboratories have not carried out anaerobic culture (17). The culture method is not sensitive enough to identify species and lacks the judgment of microbial abundance. It cannot reflect the composition of microorganisms in the wound in a timely and comprehensive manner, which will delay the treatment. Over the past decade, the accuracy of culture results has been doubted in the field of molecular microbiology (18), with a larger number of bacterial communities being identified through the development of new molecular technology (19).

The most widely used high-throughput sequencing technique for bacterial identification is 16S rRNA sequencing. The sensitivity of high-throughput sequencing detection methods to microorganisms is significantly higher than that of culture-based methods (20). Furthermore, the number of microorganisms cultured per sample was small. In this study, 16S rRNA sequencing showed that all specimens were polymicrobial, but only 41.18% of the specimens were polymicrobial by culture method. The number of bacteria obtained by sequencing is much higher than that obtained by culture method. Price et al. (21) utilized anaerobic culture, aerobic culture methods, and 16S rRNA sequencing and showed that the number of bacteria obtained by sequencing was four times than that of the culture method, which further verified the outstanding advantages of 16S rRNA sequencing in demonstrating bacterial diversity. It has been also reported that osteomyelitis is caused not by single pathogens but by a variety of bacteria, including aerobes and anaerobes, many of which cannot be cultured by traditional culture methods (22). Coincidentally, Brook et al. (23) showed a significant increase in the mortality of mice when anaerobes were mixed with aerobes. Anaerobes play an important role in the development of DFO and may cooperate with aerobes to promote progression of disease (24). In this study, the highest relative abundance in DFO was Prevotella based on the 16S rRNA sequencing, while that was Halomonas in PFO. Compared with the DFO group, patients in the PFO group seemed to be much healthier. The possible reason is that patients in the PFO group were younger and had fewer systemic diseases, the abundance of VFs in the PFO group is also lower. In general, G− bacteria and anaerobes are dominant in DFO, while G− bacteria and aerobes are dominant in PFO, providing reference for clinical empirical selection of antibiotics.

However, the selective amplification of the target gene fragment by 16S rRNA sequencing results in limited information. Even after continuous optimization, bacteria can only be identified to the genus level, and specific information at the species level cannot be obtained. The emergence of metagenome sequencing has solved this problem. Compared with 16S rRNA sequencing, metagenome sequencing is a technique that involves sequencing and analyzing the entire genome of a specimen. More detailed information on the taxonomy and genes of microorganisms can be obtained with this method. The sequencing information cannot only identify bacteria at the species level and even the strain level but also obtain functional information of the microbiome. This study first utilized metagenome sequencing to analyze the microbiome of infected bone specimens, which confirmed the feasibility of this technique in identification of the microbiome in bone specimens. That lays a solid foundation for promoting metagenome sequencing from the research field to clinical application as soon as possible. Metagenome sequencing identified a greater variety of microorganisms than culture-based methods. The majority of microbes identified by metagenome sequencing were not identified by culture, but almost all bacteria identified by culture were also authenticated by metagenome sequencing. In this study, 832 species were obtained from the metagenome sequencing of 10 infected bone specimens, including 22 dominant species. Among them, six species belong to Prevotella, which is the highest proportion of all species in the DFO group. This result was consistent with the results of 16S rRNA sequencing, thus verifying the accuracy of metagenome sequencing.

It was reported that elevated CRP values were highly predictive of amputations and may also be associated with G− bacterial infection (25). The species most positively correlated with CRP in this study was P. vulgaris, suggesting that the appearance of P. vulgaris may be related to poor prognosis of patients, such as amputation. Also, the species with the greatest negative correlation with PCT and CRP value was B. fragilis, suggesting that B. fragilis may exist in wounds of patients with DFO in the early stage of infection, so it is necessary to select sensitive antibiotics for treatment at that time. Gardner et al. (16) analyzed the correlation between three dimensions of microorganisms in DFO wounds and six clinical variables and found that the duration of ulcers was positively correlated with bacterial diversity and the abundance of P. vulgaris and negatively correlated with the abundance of Staphylococcus. Ulcer depth was positively correlated with the abundance of anaerobes and negatively correlated with the abundance of Staphylococcus. The results of this study show that bacteria belonging to Prevotella are positively correlated with the duration of DFI, suggesting that the longer the duration of DFI, the greater the chance of emergence of G− anaerobes, which indicates that treatment should be strengthened accordingly.

The rise in morbidity and mortality associated with drug-resistant microbial infections is a major global threat facing humans today. Improved comprehension of AR will provide a new perspective for redirecting antibiotic therapy and reducing bacterial resistance, which could minimize complications arising from treatment with broad-spectrum antibiotics and facilitate appropriate treatment. To avoid prolonging the patient’s condition due to poor efficacy and increased levels of drug-resistant bacteria, streptomycin, lincosamide, and macrolide are not used as a first-line empirical anti-infective treatment for osteomyelitis. However, this study merely defined a range of possible microbial drug resistance in wounds. Because of the complexity and diversity of bacterial drug resistance mechanisms, the evaluation of drug resistance is far from simple. It is not possible to determine the resistance of certain bacteria only by identifying drug resistance genes, further proof must be obtained through drug sensitivity experiments and transcriptomics.

This is the first study of the DFO microbiome conducted by utilizing metagenome sequencing. The total DNA was analyzed with sufficient depth in the samples to provide a high level of detail, pinpoint bacteria at the species level and provide information regarding their functions. Moreover, the biases introduced by assembly, database construction, and reference database annotation in metagenome sequencing are more easily understood than those in 16S rRNA sequencing. However, the ability of metagenome sequencing to provide information regarding the actual functions of the microbiome is limited because this method cannot distinguish between expressed and nonexpressed genes. To overcome this limitation, metagenome sequencing should be combined with other molecular approaches, such as transcriptomics and proteomics, to identify biological characteristics that control the expression of metabolic activity in microbial communities (26).

DFO is a major stress to health care systems and leads to significant morbidity and mortality. The economic and social burden of DFO is as the global incidence of diabetes increases. Therefore, new treatment and management methods are urgently needed. Here, we examined specimens of DFO with metagenome sequencing to determine the microbial taxonomy and its functions. We combined this analysis with clinical data to demonstrate the role of the wound microbiome in host response and wound healing. These insights may lead to improved management and treatment based on the wound microbiome. Further studies with larger populations are needed to fully understand the DFO microbiome.

Funding. This work was supported by the National Natural Science Foundation of China for Young Scholars (81600648).

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. M.Z. and Y. Cai collected specimens, analyzed data, and wrote the manuscript. P.H. collected specimens and researched data. Y. Cao, B.Z., and X.W. analyzed data. X.L., N.J., and Q.L. collected specimens. X.F. and H.Z. contributed to the discussion. Y.X. and F.G. reviewed the manuscript. F.G. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.