This observational study included seven Chilean patients (six men and one woman; P1-P7) with a mean age of 61 years, all diagnosed with DFO (IDSA 3, SEWSS, moderate PEDIS). None had received antibiotics in the two weeks prior to sampling. After providing informed consent approved by the Scientific Ethics Committee of Universidad Católica del Norte and Hospital Regional de Antofagasta (HRA), patients underwent bone biopsy in a secondary-level procedure room or private clinic’s surgical pavilion, following the protocol by Féron et al. (2021) (Féron et al., 2021).

Under sterile conditions, bone samples were taken from the exposed lesion site (Table 1). In patients P3, P4 and P7 -who had more severe lesions- samples were collected in three layers (top, middle, bottom; approximately 3 x 4 mm in size). A single sample was taken from patients P1, P2, P5 and P6 due to less extensive involvement.

Each sample was divided into three subsamples: one for bacterial culture (in a sterile 0.9% saline solution and sent to the HRA microbiology lab), one for anatomical and pathological analysis (fixed in 10% formalin and processed by the pathology lab) and one for bacterial community analysis (preserved in RNALater™ stabilization solution, ThermoFisher, USA and stored at -80°C for DNA/RNA extraction).

Subsamples (single or triplicate per patient) were assessed according to each center´s protocols. Fixation occurred within 24–48 hours, followed by 72 hours EDTA-decalcification. Standard processing produced paraffin-embedded blocks. Serial sectioning (3-4 μm thick) was stained with hematoxylin and eosin, and examined under a light microscope by a single pathologist.

Histological diagnosis was categorized as:

Upon arrival at the HRA microbiology lab, subsamples were processed immediately. After confirming they met acceptance criteria, each was numbered consecutively. In a biosafety cabinet under sterile conditions, subsamples were homogenized in brain-heart infusion (BHI) broth using one of two methods: (a) large fragments were cut with bone pliers and immersed in BHI; (b) smaller samples were directly immersed in BHI within their transport containers. All samples were vortexed for 10 seconds and incubated (37°C, 12–18 hours) (Burillo et al., 2007).

Post-incubation, the homogenates were plated on blood, chocolate, and MacConkey agar using sterile swabs and loops. Blood and chocolate agar were incubated (37 ± 2°C, 5% CO₂ for 18–24 hours); MacConkey agar was incubated aerobically for the same duration. Plates without growth were incubated up to 48 hours.

All colonies were considered significant and underwent Gram staining. Identification and susceptibility testing were performed with the Vitek^®2 system, following CLSI 2016 guidelines. In rare or ambiguous phenotypes, the Kirby-Bauer method and/or Epsilon test were used (Bauer et al., 1959; Parada et al., 2016).

Anaerobic bacteria were not assessed, as this was not part of the standard procedures at HRA.

A total of 13 subsamples were processed in 2 analytical groups: One single subsample from each P1, P2, P5 and P6 was used for both gDNA and RNA extraction, while single subsamples from the top, middle and bottom layers from P3, P4 and P7 were used only for gDNA extraction. The gDNA and RNA extractions were performed using PowerSoil DNA and RNA isolation kits (MO BIO Laboratories Inc., CA, USA) following the manufacturer’s instructions. RNA was reversed-transcribed into cDNA using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems) from 1 µg/µL of total RNA per subsample, to profile the ribosomally active bacterial community.

The qualities of gDNA and cDNA were assessed using the Qubit dsDNA HS assay (Invitrogen, CA, USA). PCR amplification of the 16S rRNA gene using universal primers 27F (forward; 5′-AGAGTTTGATCATGGCTCAG-3′) and 1492R (reverse; 5′-GGTTACCTTGTTACGACTT-3′) was performed to verify template integrity. Products visualized on 2% agarose gels stained with ethidium bromide under UV light, with Escherichia coli 16S DNA as a control.

For profiling the total bacterial communities and ribosomally active fraction, the V4–V5 regions of the 16S rRNA gene were amplified from gDNA and cDNA, respectively, using primers 515F-Y (5′-GTGYCAGCMGCCGCGGTAA-3′) and 926R (5′-CCGYCAATTYMTTTRAGTTT-3′) (Parada et al., 2016). Amplicons were sequenced (2×301 bp, paired-end) on a MiSeq platform at Macrogen Inc. (Seoul, Republic of Korea), following the manufacturer’s protocols.

Raw data results were processed using Mothur pipeline. Rare sequences differing by 1–2 base pairs were preclusted and chimeras were removed using VSEARCH (RRID: SCR_016964) (Schloss et al., 2009). Taxonomic classification was performed, non-bacterial sequences were removed, and OTUs were calculated at a distance of 0.03 using cluster.split command. For family-level taxonomic distribution (P1–P7), OTUs with identical taxonomic classification were retained as distinct units.

Sequences (from Macrogen) were analyzed with Mothur 1.39.5 and the SILVA NGS reference database (v132: 152308 Bacteria, 3901 Archaea and 16209 Eukarya sequences; RRID: SCR_006423). Data were processed according to the MiSeq SOP protocol. Richness and diversity metrics-including Observed OTUs (Subject observation index, S_obs), Chao1, Shannon and Simpson indices-were calculated and singletons and chimeras were excluded. Relative bacterial abundances at the phylum, family and genus levels were calculated for both total (gDNA) and ribosomal active (cDNA) bacterial communities.

To generate microbial interaction networks,.shared and.taxonomy files from Mothur (RRID: SCR_011947) were converted into Cytoscape-compatible format using a modified RStudio script (RRID: SCR_000432) (Neave et al., 2014). Nodes represented samples and total OTUs per sample; edges were assigned as distances between the OTU node and sample. The network was visualized in Cytoscape v.3.5.0 (RRID: SCR_003032) (Shannon et al., 2003) using the Force Atlas Algorithm.

For final microbial network, Gephi software (RRID: SCR_017284) was used with the following Force Atlas parameters: inertia, 0.1; repulsion force, 100.0; attraction force, 100.0; autostabilization (force, 100.0; sensitivity, 0.2); gravity, 55.0; attraction distribution, 1.0; speed, 1.0. Node sizes were adjusted to highlight sample sources.

The patients included in this study were aged between 55 and 69 years, comprising six men and one woman. All patients had a history of diabetes mellitus lasting more than 10 years and, at the time of bone sampling, presented with an average glycated hemoglobin (HbA1c) level of 8.5 ± 1.1%. Renal function was impaired in all cases, with one patient diagnosed with renal failure. The average duration of OPD treatment was 12 months, except for patient P2, who received intermittent antibiotic treatment for 36 months until his death following a major amputation. This patient was the only one to die from complications related to diabetic foot over the eight-year follow-up period. All remaining patients had a survival period of more than four years. Five of them died during the COVID-19 pandemic: two due to virus-associated pneumonia and two from cardiovascular causes (stroke and cardiopulmonary arrest). One patient died after three years as a result of sepsis secondary to a diabetic foot infection, also associated with a major amputation (Table 1). While some of patients were within the high-risk age groups for COVID-19 and had significant cardiovascular risk factors, the deaths observed in the study were unrelated to the underlying DFO or the successful outpatient treatment received. These deaths did not affect the study outcomes because the metagenomic and ribosomal activity analyses were conducted on samples obtained prior to these clinical events.

In five out of the seven patients, three histological samples were obtained from different bone depth layers (P2–P4, P6, and P7). Acute osteomyelitis was diagnosed in four of these patients, based on suppuration due to neutrophilic infiltration, the presence of pyocytes or necrotic debris, osteonecrosis, and visible bacteria. These patients also exhibited osteopenia, fibrosis, or bone fragmentation. Two patients presented with chronic osteomyelitis, characterized by bone marrow fibrosis and plasma cell infiltration, without visible bacterial presence (Table 1).

A total of 19 bacterial strains were isolated from 13 bone samples collected from seven patients with histologically confirmed DFO. These strains corresponded to eight different bacterial species. Of the isolates, 95% were Gram-negative bacilli, with the most prevalent being Proteus mirabilis (26.3%; 5 isolates), Enterobacter cloacae (26.3%; 5 isolates), and Pseudomonas aeruginosa (21.1%; 4 isolates). Less frequent species included Providencia stuartii, Escherichia coli, Morganella morganii and Pseudomonas putida (each 5.3%; 1 isolate). The remaining 5% corresponded to a single Gram-positive coccus, Staphylococcus aureus (Table 2). In 38% (5/13) of the samples, only one bacterial species was identified. In 54% (7/13), two species were isolated. One sample showed no bacterial growth after 48 hours of incubation.

From 16S rRNA gene sequencing of both genomic DNA (gDNA) and complementary DNA (cDNA) derived from patient bone samples, a total of 2,073,017 reads and 3,412 operational taxonomic units (OTUs) were identified. Sample coverage ranged from 0.98 to 0.99, ensuring reliable detection of taxa and accurate estimates of richness and diversity (Supplementary Table S1).

Richness scores varied significantly among samples and sequencing type. For the total bacterial community (gDNA), the observed OTUs (S_obs index) ranged from 62 to 515. Higher richness was observed in samples P3, P4, P5, P6, and P7 (282–515 OTUs), while P1 and P2 had lower values (22–62 OTUs). The Chao1 estimator showed similar trends, with values ranging from 1351.67 to 30,709.15 for P3–P7 and from 77 to 130 for P1 and P2. Ribosomal transcript–based bacterial communities (cDNA) followed the same pattern (Figure 1).

Shannon and Simpson diversity indices revealed that the highest bacterial diversity was observed in all P4 samples (Shannon: 2.03–5.67), followed by samples P1, P2, P5, P6, and all P7 samples (1.01–6.12). P3 samples exhibited the lowest diversity (Shannon index from 1.25 × 10^-7 to 1.13). Interestingly, P1 and P2 showed higher ribosomal activity (Simpson index: 2.97–15.59) compared to P5 and P6 (1.22–3.4) (Figure 1).

The identified OTUs were classified into 12 phyla, 11 classes, 91 families, and 118 bacterial genera. At the phylum level, the most abundant taxa across all samples were Proteobacteria (58.4%), Bacteroidetes (29.8%), Firmicutes (10.3%), and Actinobacteria (≥0.1%) (Figure 2A). At the family level (Figure 2B), the top ten most abundant families were Pseudomonadaceae (27.3%), Enterobacteriaceae (26.9%), Prevotellaceae (19.5%), Bacteroidaceae (9.6%), Family_XI (3.0%), Moraxellaceae (4.5%), Enterococcaceae (2.8%), Peptostreptococcaceae (2.3%), Staphylococcaceae (0.8%) and Corynebacteriaceae (0.6%).

For the total ribosomally active bacterial community, the dominant phyla were similar: Proteobacteria (43.5%), Firmicutes (20.7%), Bacteroidetes (19.8%), Actinobacteria (13.6%), and Chloroflexi (9.4%) (Figure 3A). The most abundant families included Enterobacteriaceae (18.5%), Bacteroidaceae (11.8%), Enterococcaceae (9.6%), Propionibacteriaceae (8.1%), Burkholderiaceae (4.6%), and others (Figure 3B).

The total bacterial composition for each patient was determined at the taxonomic family level, 16S rDNA-based analysis of DFO samples revealed a diverse bacterial community representing the total microbiota, including viable, nonviable and dormant bacteria. Enterobacteriaceae was the only family detected in all samples, with highly variable relative abundances (1.6%–91.1%). The dominant family varied by patient, for example: Prevotellaceae dominated P6 (61.8%) and all P7 depths (60.3%), Enterobacteriaceae dominated P1 (52.1%), P2 (91.1%), P4 (36%), and P5 (83.8%), while Pseudomonadaceae dominated all depths of P3 (94.8%). On the other hand, one patient showed a more balanced microbial profile, with several low-abundance families. The bacterial family composition and its variability across samples are illustrated in Figure 4.

At the genus level, on average, 37 ± 15 bacterial genera were identified per sample. Patient-specific dominant genera were observed, with Proteus (P3), Pseudomonas (P1–P2), and Prevotella 7 (P7) predominating. The most prevalent Gram-negative aerobic genera were Pseudomonas, Acinetobacter, Sphingomonas and Alcaligenes. Notably, anaerobic bacteria predominated, accounting for 10/18 genera with the highest RA%, with Anaerococcus being the only genus detected in all samples, despite its low relative abundance (8.7%–0.3%). Staphylococcus (~3.5%), a facultative anaerobe, was detected only in patient P4. Overall, the coexistence of anaerobic, aerobic and facultative taxa underscores the polymicrobial nature of DFO (Figure 5A). Within the ribosomally active bacteria showing the highest RA%, Cutibacterium was present in all samples, although at low relative abundance (27% to 1.6 × 10^-2%), whereas Morganella (90%) and Bacteroides (64%) were exclusively detected in P5 and P6, respectively. Consistent with metagenomic analysis, anaerobic genera (6/18) and facultative anaerobes (8/18) predominated over aerobic genera (2/18). Several genera, including Bacteroides, Proteus, Finegoldia, Bilophila, Acinetobacter, Peptoniphilus and Corynebacterium, were detected in both 16S rDNA- and 16S rRNA-derived datasets, indicating their presence in the total community as well as their association with the ribosomally active fraction (Figure 5B).

To understand the spatial distribution of the microbiota in DFO, a gDNA-based metagenomic analysis was performed across bone depth layers (top, middle and bottom) for patients P3, P4, and P7. Six families were found exclusively in the superficial layer (Top), with Methylophilaceae being the most represented (>1%). The middle layer had the highest diversity, with 19 exclusive families, notably Fusobacteriaceae (>1%). In the bottom layer, 10 unique families were identified, including Micromonosporaceae as the most abundant (>1%). Several families were shared between layers: six between middle and bottom, and five between top and middle. Dominant families present across all layers included Enterobacteriaceae (52.5%), Family_XI (5.96%), Moraxellaceae (5.01%), Peptostreptococcaceae (4.66%), Enterococcaceae (3.95%) and Corynebacteriaceae (52.09%) (Supplementary Figure S1).

Samples from P3 (all depths) and P5 showed the most homogeneous microbial composition and strongest interconnections. P3 samples were closely related to P4 and P6, while P7 samples clustered more closely with P1 and P2, though sequence numbers in these were lower (Supplementary Figure S2).

At the genus level, samples from P3 were predominantly composed of Pseudomonas across all layers (>90%). In P4, Pseudomonas (~19.1%) was detected together with other unclassified Enterobacteriaceae (~34.2%), while Prevotella accounted for more than 60% of the community in P7 (Figure 5A). Consistently, as bone depth increased, the relative coverage of the most prevalent anaerobic genera also increased.

Using cDNA sequencing of samples P1, P2, P5, and P6, ribosomally active bacterial at phylum and family levels were identified. The dominant phyla were: Proteobacteria (43.5%), Firmicutes (20.7%), Bacteroidetes (19.8%), Actinobacteria (13.6%), and Chloroflexi (9.4%) (Figure XA). The most abundant families included Enterobacteriaceae (18.5%), Bacteroidaceae (11.8%), Enterococcaceae (9.6%), Propionibacteriaceae (8.1%), Burkholderiaceae (4.6%), and others (Figure XB). At family level, Enterobacteriaceae (36.5%) and Enterococcaceae (28.3%) were the most ribosomally active, followed by Pseudomonadaceae (4.2%), Corynebacteriaceae (4.0%), Family_XI (2.9%), Lachnospiraceae (2.3%) and Bacteroidaceae (1.3%). Notably, Bacteroidaceae was the only family present in both total and active communities across all samples (Figure 6).

Comparisons between layer-specific and non-layered samples revealed that ribosomally active taxa were detected in all bone layers, although their distribution varied (Supplementary Figure S1). Some families were exclusive to the superficial/top (2), middle (3) or bottom (7) layers, while others were shared across two or all three layers. Overall, 67.8% of the taxa identified by 16S rRNA overlapped with those detected by 16S rDNA. The active-to-total ratio by layer was: 21/37 for superficial (top), 23/56 for middle, and 26/54 for bottom samples. Across all layers, the most abundant ribosomally active families were Enterobacteriaceae, Enterococcaceae, Corynebacteriaceae, Pseudomonadaceae, Family_XI, Lachnospiraceae, and Bacteroidaceae, suggesting their potential involvement in DFO-associated microbial activity.

Differences between DNA- and RNA-based profiles indicate that not all taxa detected at the DNA level are ribosomally active, while some taxa with low genomic abundance show increased ribosomal activity. Finally, most detected genera correspond to known cultivable bacteria, although several are fastidious or require strict anaerobic conditions; only one taxon (AKIW781_ge) represents an uncultured bacterial lineage identified exclusively by 16S rRNA gene sequencing.