Introduction

Skin ulcers in the lower leg or foot ulcers are a leading cause of impaired quality of life, hospitalization, and, in the worst cases, of amputation in diabetic patients1. The wound healing process is severely compromised in diabetes, and is often complicated by neuropathies, that reduce sensation in the limb and make the patient more insensitive to traumatic wounds. In addition, diabetic vascular dysfunction leads to stenotic limb conditions and decreases angiogenesis, which is fundamental to ulcer healing.

Diabetic ulcers are often characterized by persistent, unresolved inflammation and are highly susceptible to infection, both due to prolonged wound exposure to environmental factors and to impaired immune responses2,3. Pathogenic bacteria such as S. aureus and Pseudomonas species are commonly present1. Commensal bacteria (S. epidermis and Corynebacterium species) coexist with anaerobic bacteria such as Peptostreptococcus species, competing or cooperating with each other and modulating virulence1. Indeed, the prevalence of anaerobic bacteria has been observed in the diabetic ulcer lesion4. Moreover, the formation of an ulcer bacterial biofilm contributes to the persistence of infection and increases resistance to classical antibiotic treatment5. To limit bacterial biofilm formation within the diabetic ulcer, debridement is commonly used to remove the necrotic, devitalized tissue and better preserve the healthy skin6.

In recent years, studies on the composition of microorganisms (microbiota) stably residing in our body and in particular those of the skin (cutaneous microbiota) have multiplied thanks to the high throughput sequencing technology7. 16S ribosomal RNA gene sequencing has been shown to provide a more comprehensive assessment of the skin bacterial communities than traditional microbiological cultivation methods8.

The cutaneous microbiota plays an essential role in protecting healthy skin from pathogenic microorganisms, in shaping the responses of the skin immune system, and in the metabolism of some natural products such as antimicrobial peptides, short-chain fatty acids, and polyamines7. Lower leg and foot, the anatomical zones where diabetic ulcers mostly occur, comprise skin areas with moist (plantar), or dry (lower leg, instep, and heel of the foot) characteristics and consequently with different microbial inhabitants9. In addition, the cutaneous microbiome is dynamic during human lifespan, with physiological changes such as hormonal changes and lifestyle10.

Many skin diseases are associated with changes in the cutaneous microbiota (dysbiosis). These changes may also be associated with specific therapeutic treatments11,12,13 and may represent potential biomarkers of response to therapy14.

Some studies on the skin microbiota of individuals with type 2 diabetes mellitus have been performed both in mouse models15 and in human patients16,17,18,19,20. Significant differences have been found when comparing the skin microbiome between diabetic patients and healthy individuals. For example, S. epidermis and S. aureus have a higher prevalence in diabetic patients16. Another study showed that dynamic changes in the cutaneous microbiome occur during the onset, progression, and long-term establishment of diabetes21. The proportion and diversity of cutaneous microbiota genera reported in different literature articles are not always consistent20,22,23,24, but this dysbiosis increases the risk of developing infections in the case of skin injuries and could have an impact on wound evolution to chronicity19,23,25,26. Using a microbiota-transplanted rat model, Kunimitsu et al. showed that dysbiosis significantly affects wound dimension and the presence neutrophil and regulatory T cell infiltration27.

Despite the ever-increasing knowledge of the cutaneous microbiome in the diabetic skin, changes in the cutaneous microbiota specifically associated with diabetic ulcer healing and the influence of medications on microbiome composition have been little investigated so far20,24,28.

Therefore, we decided to address this issue by analyzing changes in the cutaneous microbiome within the ulcer bed of ten diabetic patients before and after 20 days of therapy with a fluorescein-based galenic treatment named Fluorexin. This approach allowed us to induce ulcer healing without the use of antibiotic therapy and to investigate the composition of the cutaneous microbiome that correlates with the promotion of healing. To this intent, cutaneous microbiome was also compared between the ulcer bed and the perilesional tissue prior to initiation of therapy.

Results

Patient characteristics and treatment responses

To study changes in the skin microbiome during the healing process of diabetic ulcers, 10 adult patients were enrolled, 9 with type II diabetes and 1 patient (pt_009) with type I diabetes, with a grade I or II lower leg or foot ulcer. For each patient, ulcers were photographed (Fig. 1A), and the following characteristics were considered: age, sex, comorbidities, glycated hemoglobin (HbA1c) levels and diabetes therapy, degree and ___location of the ulcer, body mass index (BMI), patient-reported smoking or physical activity. As shown in Tables 1, 4 patients (40%) were female, and the mean age of the patients was 71 years. BMI average value was 26.62 kg/m2, with 3 patients (30%) in the normal range, 5 patients (50%) overweight, 1 patient (10%) class I obese, and 1 patient (10%) class II obese. Only 1 in 10 patients reported smoking, and 5 patients reported light physical activity. Regarding comorbidities, hypertension and/or heart problems were highly represented (80%). HbA1c levels varied from 37 to 137 mmol/mol, with 20–48 mmol/mol considered in the normal range. Diabetic therapy was insulin in 7 patients, metformin in 2 patients, and was not assessed in 1 patient who had normal HbA1c levels.

As shown in Fig. 1B and Table 1, ulcers were located on the heel (2 patients), plantar midfoot (3 patients), mid/lower leg (2 patients), shin, malleolus, and little toe (1 patient each). Prior to the initiation of therapy (t0), stage, according to the University of Texas Staging System for Diabetic Foot Ulcers29, and Falanga’s score30 were measured, debridement was performed, and pH was measured in the ulcer bed. Out of 10 patients, 8 patients (80%) had a pH value between 5.0 and 5.5 at t0 (Table 1). The ulcers were then treated with the Fluorexin solution, a galenic compound that combines fluorescein, known for its anti-inflammatory properties; Uncaria tormentosa extract, which has both anti-inflammatory and immunomodulatory properties; allantoin, a hydrating and lenitive molecule; and Malaleuca Alternifolia extract, which has antiseptic and mild anesthetic properties. This treatment was chosen to promote ulcer healing without the application of topical classical antibiotics, which would interfere with microbiome analyses. In addition, the Fluorexin solution has the advantage of being lightly colored, unlike other traditionally used antiseptic compounds based on fluorescein derivatives (i.e. eosin), which are strongly colored and difficult to remove. Patients were instructed to continue daily application of the Fluorexin solution to the ulcer.

Table 1 Clinical characteristics of the enrolled patients.
Figure 1
figure 1

Ulcers characteristics. A. Photographs of diabetic ulcers at day 0 (T0; a, c, e), and after 20 days of treatment with Fluorexin (tT0; b, d, f). Representative images are shown of pt_001 (a, b); pt_006 (c, d); pt_007 (e, f). B. Stage, according to the University of Texas staging system for diabetic foot ulcers29, where stage A comprises pre- or post-ulcerative lesions completely epithelialized (0); superficial ulcer, not involving tendon capsule or bone (1); ulcer penetrating to tendon or capsule (2). Stage B indicates the presence of infection in addition to what reported for stage A. Falanga’s score30 was calculated for each patient ulcer: (A) granulation tissue presence in the 100% of the ulcer; (B) granulation tissue presence in the 50-100% of the ulcer and presence of fibrinous tissue; (C) granulation tissue presence in <50% of the ulcer and presence of fibrinous tissue.

After 10 days of treatment (t10), patients returned for an interim clinical evaluation. After a further 10 days (t20) of therapeutic treatment, ulcers were photographed (Fig. 1A), stage, Falanga’s score, and pH were measured, and in the case of open wounds, the ulcer was debrided, and the patient was referred for further therapeutic treatments if necessary. As shown in Fig. 1B, ulcer stage improved in 9 out of 10 patients. Six patients (60%) started with grade 2B at t0 and in 2 patients (20%) the stage decreased to 1B at t20, whereas in 4 patients (40%) the stage decreased to 2A. In the 3 patients (30%) who started with a 1B stage at t0, the stage decreased to 1A or 0 at t20. One patient (pt_010) showed a stage increase from 1A at t0 to 2A at t20, but this patient experienced an initial worsening of the ulcer lesion, with a stage increase at t10 (2B). Therefore, from t10 to t20 the patient had a clear improvement, and the ulcer was considered to be in the healing phase at t20.

As for the Falanga’s score, the improvement was seen in 3 out of 10 patients, while 5 maintained the same score and 2 showed worsening, including pt_010 (Fig. 1B). This patient also experienced a worsening of the Falanga’s score at t10 (from A to D) and a subsequent improvement (from D to C). On the other hand, pt_009 showed an improvement in the score (from 2B to 2A), but at t20 a small point of the ulcer bed showed a clear worsening of the granulation tissue. Therefore, the Falanga’s score was reported to change from A to B.

Microbiome composition in the ulcers of diabetic patients at t0 versus t20

DNA was extracted from swabs taken in the ulcer bed at t0, and in the ulcer bed or on the reconstituted tissue at t20, depending on the clinical situation of the ulcer. Next Generation Sequencing (NGS) of the 16 S rRNA was then performed, and the data obtained were subjected to bioinformatic analysis. In each run, > 1.2 giga-bases were sequenced and a final library of 5 mega-bases was obtained in both runs, 5.068550 and 5.405278, respectively. Average counts per sample were 437,200 (range: 235200–719400) with 1760 operational taxonomic units (OTUs) with at least 2 counts. As shown in Fig. 2 and Supplementary Table 1, analysis of the microbiome composition at t0 and t20, considering the 20 most representative bacterial genera, indicated that bacteria colonizing the diabetic ulcers before therapy were mainly of the genus Staphylococcus (35.4%), which was reduced by therapy at t20 (24.9%). Corynebacterium was the second most represented genus and increased with therapy from t0 (12.6%) to t20 (24.5%). The genera Anaerococcus (8.2%) and Methylobacterium (1%) were stably present at t0 and t20, whereas the genera Pseudomonas (from 5.4 to 9.1%), Propionibacterium (from 3.3 to 6.7%), Streptococcus (from 1 to 2.9%), Arthrobacter (from 0.4 to 1.6%), and Sphingomonas (from 0.2 to 1.2%) increased during the 20-day Fluorexin treatment.

Figure 2
figure 2

Ulcer bed microbiome mean data at t0 and t20. (A) Stacked bar plot of the most abundant 20 genera found in the lesional areas before (t0) and after the 20-day treatment (t20) with Fluorexin. Text layout in the plot was edited with GNU Image Manipulation Program (GIMP v2.10.34) (B) Beta diversity in lesional areas before (t0) and after the 20-day treatment (t20). Principal coordinate analysis (PCoA) plots were performed by the Bray-Curtis distance measures. Statistical significance was assessed by PERMANOVA.

An important reduction in Enterococcus (from 6.5 to 1%), Campylobacter (2.4%-0.6%), and Alcaligenes (1.5%-0.3%) genera at t20 was observed.

The category named as “others” included the bacterial genera after the 20 most represented ones.

The differences were not significant due to the high variability among samples and to the small cohort of patients (Supplementary Table 2).

The relative abundance of some of the species found in the intralesional samples at t0 and t20 revealed that the Gram-positive aerobe S. aureus was the most abundant species at t0 and tended to be reduced by therapy (Table 2). S. pettenkoferi, S. warneri, S. epidermidis, and S. simulans, were also present. As shown in Table 2, an increase in Gram-positive P. acnes and in Gram-negative P. aeruginosa was observed with this therapeutic approach. As for Streptococcus, and Corynebacterium genera, the most represented species were S. intermedius and C. striatum, C. tuberculostearicum, and C. simulans, respectively.

Only one patient (pt_002) was positive at t0 for the Gram-negative anaerobe P. mirabilis, which was importantly reduced after therapy. The relative abundance of A. faecalis was also reduced at t20.

Table 2 Relative species abundances of the 7 most common genera found in intralesional ulcer samples at t0 and t20. Records are sorted alphabetically by genus, and species are reported if their mean abundance at t0 was greater than 0.2% and the absolute difference (t20-t0) was also greater than 0.1%. For each time point, microbiome data from the 10 enrolled patients were collected and mean, 5th and 95th percentiles were calculated.

Microbiome changes over the course of therapy for each patient

Examining the differences in microbiome composition for each of the 10 patients enrolled in the study, we found heterogeneous relative abundance of genera among the patients (Fig. 3).

Two patients (pt_001 and pt_009) manifested a high abundance of Staphylococcus at t0, and the amount did not change with therapy for pt_001, whereas it was reduced for pt_009 at t20. Corynebacterium genus was present in 7 out of 10 patients at t0 and increased at t20 with the only exception of patient pt_007, where it was slightly reduced. The genus Serratia was present in only one patient (pt_008). Interestingly, microbiological culture analysis performed at t10 for this patient revealed the presence of S. marcescens, further confirming the NGS analysis. Therefore, the reduction in Serratia genus we observed when comparing t0 and t20, would likely have occurred between t10 and t20.

Figure 3
figure 3

Stacked bar plot of the most abundant 20 genera found in lesional areas before (t0) and after the 20-days treatment (t20) for each of the 10 patients enrolled in the study. Text layout in the plot was edited with GIMP.

Comparison of intralesional versus perilesional microbiome composition at t0

With the exception of pt_007, DNA was also extracted from swabs taken into the ulcer bed (intralesional) and in the perilesional tissue of the ulcer at t0. For the perilesional tissue, a sample was taken by the clinician at t0 in the uncompromised area near the ulcer. Sampling was not repeated at t20 because the margins between closed wound tissue and the surrounding healthy tissue were not always clearly visible, and the microbiome composition in the perilesional area could have been influenced by an abundant application of Fluorexin by the patient which could have involved this area as well. 16S rRNA NGS was then performed, and the data obtained were subjected to bioinformatic analysis.

As shown in Fig. 4 and Supplementary Table 3, no major differences were observed in the microbiome profile. The genus Staphylococcus was slightly more abundant in the perilesional area than in the intralesional area (42% versus 37%), whereas the genus Finegoldia was more abundant in the intralesional area than in the perilesional area (5.4% versus 2.4%). Corynebacterium, Anaerococcus, Enterococcus, Pseudomonas, and other genera listed in Fig. 4 and Supplementary Table 3, were equally present in the intralesional and perilesional zones at t0.

The differences were not significant due to the high variability among samples and to the small cohort of patients (Supplementary Table 4).

Figure 4
figure 4

Microbiome mean data of intra and perilesional tissue samples at t0. A. Stacked bar plot of the most abundant 20 genera found in the intralesional and perilesional areas before Fluorexin treatment (t0). B. Beta diversity in the intralesional and perilesional areas at t0. Principal coordinate analysis (PCoA) plots were performed by the Bray-Curtis distance measures. Statistical significance was assessed by PERMANOVA.

When the analysis was performed on the individual patients (Fig. 5), we found that one patient (pt_001) had a predominance of Staphylococcus genus, almost identical in both zones. In three other patients (pt_005, 009, and 010), Staphylococcus genus was highly present in both zones, and in two patients (pt_002 and pt_003) predominantly in the perilesional area. One patient (pt_006) has a prominent presence of Enterococcus genus in both intra and perilesional tissues, whereas Corynebacterium was highly represented in 6 patients, with variation between intra- and perilesional areas among the patients.

Figure 5
figure 5

Stacked bar plot of the most abundant 20 genera found before treatment (t0) in 9 patients, evaluating the microbiome community in the ulcer bed (intralesional), and in the perilesional tissue. Perilesional swab was not available for pt_007 that was excluded. Text layout in the plot was edited with GIMP.

Microbiome differences related to Falanga’s score or ulcer ___location

We analyzed the microbiome data by comparing the ulcer Falanga’s scores and dividing them into three groups: patients who had an amelioration of the Falanga’s score during the 20-day therapy (better), patients who maintained the same Falanga’s score (unchanged), and two patients who had a worsening of the Falanga’s score at t20. Alpha diversity was measured using original counts and applying Chao1, Shannon and Simpson diversity metrics. While the last ones account for both richness (count) and evenness (distribution) of the microbial community, Chao1 considers observed OTUs with account for unobserved species based on low abundance OTUs (richness only). Welch t-tests for post-hoc comparison were performed and showed no significant differences when comparing t0 and t20 in each group (Supplementary Fig. 1).

To analyze the NGS microbiome data in relation to ulcer ___location, we considered only the three body sites where two or three ulcers were localized. In particular, microbiome relative abundance in ulcers localized to the foot plantar (pt_002, 003, and 004), in the lower leg (pt_006, 007), and in the heel (pt_001, 008) was examined. As shown in Fig. 6A and B, respectively, the relative abundance at t0 (Fig. 6A) was different in the three groups, whereas at t20 (Fig. 6B) the differences between the plantar and lower leg ulcers were significant (FDR = 0.0053). In fact, the homogeneity between the plantar and the lower leg ulcers was low at t0 (Fig. 6C), whereas the relative abundance of microbiome genera was highly similar among the two lower leg samples and the three plantar samples, respectively, at t20 (Fig. 6D).

Figure 6
figure 6

Microbiome data on ulcer body site. A. Stacked bar plot of the most abundant 20 genera found in the intralesional samples of patients with ulcer located at plantar, lower leg or heel before treatment initiation. B. Stacked bar plot of the most abundant 20 genera found after 20 days of treatment (t20) in intralesional samples collected from patients with ulcer located at plantar, lower leg or heel. C. Top 20 taxa at genus level of the 5 patient samples from plantar or lower leg ulcers at t0. The table below summarizes the result of post-hoc pairwise comparison (multi-group only) performing the regular Welch t-tests. The multi-testing adjustment is based on Benjamini-Hochberg procedure (FDR). No statistically significant difference was observed in the alpha diversity between groups (FDR always >0.05). D. Top 20 taxa at genus level of the 5 patient samples from plantar or lower leg ulcers at t20. The table below summarizes the result of post-hoc pairwise comparison (multi-group only) performing the regular Welch t-tests. Statistically significant difference was observed in the alpha diversity between the plantar and the lower leg group (FDR = 0.005).

Discussion

Studies of the cutaneous microbiome have the potential to provide valuable information regarding the presence of microbial species during the course of a disease, and in response to therapeutics, although these aspects have mainly been studied in skin pathologies such as psoriasis and atopic dermatitis11,12,13. Therefore, our study was designed to address the changes in microbiome composition during the healing process of a diabetic ulcer. To this end, we selected ten diabetic patients with grade I or II foot or lower leg ulcers. We homogenously induced ulcer healing by the application of a galenic treatment and without the administration of topical or systemic antibiotics which would have interfered with the microbiome analysis. NGS analysis was then performed at t0 and after only 20 days of treatment to detect early changes.

Altogether, our study demonstrated a healing benefit for stage I and II ulcers with the application of a simple anti-inflammatory and anti-septic galenic solution. The therapeutic action of Fluorexin may be both a reduction in local inflammation and an inhibition for biofilm formation, since both inflammation and biofilm have been correlated with altered wound closure19.

Previously, in the comprehensive work by Kalan et al., microbiome analysis was performed in 100 neuropathic diabetic plantar ulcers of different grades, by tracking microbiome changes for 26 weeks28. During this period, patients underwent different treatments and microbiome composition was retrospectively examined based on ulcer outcome (healed, unhealed, foot amputation). Interestingly, the abundance of genera in their dataset was quite similar to what we found in our smaller cohort, with a predominance of Staphylococcus, followed by Corynebacterium, Pseudomonas, and Streptococcus. This pattern is also consistent with that found in other studies23,25. Unlike Kalan et al., here we evaluated ulcer healing by reporting the differences in stage and Falanga’s score after only 20 days and found that the healing process induced a decrease of Staphylococcus, and an increase of Corynebacterium, Propionibacterium, Streptococcus, and Pseudomonas genera. This was also true for one patient who experienced an initial worsening of the ulcer grade and a delayed healing process, whereas in one patient Staphylococcus genus was equally present at t0 and t20 and in the intra- and perilesional tissue at t0. From a clinical point of view, this ulcer was reported to be still infected at t20 although an improvement was observed (stage change from 2B to 1B). Kalan et al. pointed out that the relative abundance of different S. aureus strains in the ulcer bed could be related with healing28. However, our sequencing strategy does not have sufficient resolution to allow this investigation.

We could not even find a microbiome pattern that could explain the worst Falanga’s score registered at t20 in two patients (pt_009 and pt_0010).

Kalan et al.28 also highlighted that diabetic wounds with A. faecalis biofilms showed an improved initial healing, although no association was found between the presence of A. faecalis and ulcer outcome. Indeed, a recent paper by White et al. showed that treatment of wounds with A. faecalis in a diabetic mouse model balanced matrix metalloproteinase expression and favored wound epithelialization31. In our patients, we found that A. faecalis was reduced in the ulcer during the healing phase, consistent with a role only at the beginning of the healing process. It is likely that some components of the microbiome are important in promoting wound healing, but they should be removed before complete tissue regeneration of the tissue is achieved.

Large individual differences did not allow us to extrapolate a microbiome signature that would be prognostic for diabetic ulcer healing. However, our study showed that the healing process was mainly associated with a reduction of S. aureus, and an increase of Corynebacterium and Propionibacterium genera. A future analysis on a larger cohort of patients will be important to reduce the impact of individual differences and to support the conclusions we obtained.

Comparing the microbiome in the ulcer bed and in the perilesional tissue of the same patient at t0, we did not observe major differences in the microbiome constituents, supporting the idea that the microbiota within the ulcer could be derived from the dissemination of commensal bacteria living in the perilesional skin. Microbiome studies have previously been performed in hard-to-heal pressure ulcers, and a correlation between healing status and the microbial similarity of the wound and peri-wound tissue was found32. No such a correlation was found in our patients. In fact, the alpha diversity analysis based on the improvement or not of the Falanga’s score registered after the 20-day treatment with Fluorexin, did not show any significant data.

In a previous paper from Gardiner et al., the microbiome present in foot skin and wound of eight diabetic patients and in the foot skin of eight healthy controls was compared24. They found a significant reduction in the diversity of the skin microbiome in the diabetic wound compared to intact diabetic skin. In our study, a similar diversity was found between intralesional and perilesional tissue, probably due to the proximity of the two areas in our case, whereas the previous work considered the other foot without the ulcer lesion.

Interestingly, we found that at t20 the microbiome present in the plantar ulcers was significantly different from that present in the lower leg. The significance of our data was probably due to the high uniformity of microbiome genera present in the three plantar and in the two leg samples, respectively. These data strongly suggest that, with the restoration of tissue integrity, the microbiome present in the two area is likely to be different, due to diverse skin characteristics between the plantar foot and the lower leg, such as moist or dry skin, presence of sebaceous glands, different skin temperature9,26. Moving towards healing, we found that at t20 foot plantar ulcers showed a preponderance of Corynebacterium and Anaerococcus, whereas lower leg ulcers showed a similar presence of Propionibacterium, and Staphylococcus.

To optimally treat a diabetic ulcer, it is essential to elucidate which bacteria in the complex ulcer microbiome are pathogenic, and to distinguish colonization from infection. A better understanding of both host-microbiota interactions and cross-relationships between commensal and pathogenic bacteria would be highly useful in selecting the most appropriate therapeutic protocol for each patient. This aspect is extremely important in the management of a diabetic ulcer. In fact, in an era of increasing antibiotic-resistant bacterial species, it is intriguing to consider that the resolution of the infections that often accompany the onset of diabetic ulcers could be achieved directly by acting on the relative abundance of specific skin commensal microbiota, able to compete with and control the growth of pathogenic species. Although obtained from a small cohort of adult patients, our data provide further knowledge into the changes in the diabetic ulcer microbiome during the healing process, representing an additional step towards the design of innovative therapeutic protocols based on microbiome manipulation.

It should be considered that the patient microbiome also includes fungi and viruses, the role of which has not been thoroughly investigated. Therefore, further studies are needed to gain a complete understanding of the microbiota that would be more beneficial for proper wound healing and control of pathogenic dysbiosis.

Methods

Patients and clinical assessment

This study was conducted according to the Good Clinical Practice Guidelines and the Declaration of Helsinki. The study was approved by the Ethical Committee of IDI-IRCCS (protocol n. 621/2021). All patients enrolled in the study provided written informed consent. Inclusion criteria were: adult diabetic patients of both sexes, suffering from grade I or II diabetic ulcer, stage A or B according to the University of Texas staging system for diabetic foot ulcers29, who had not received antibiotic therapy, either systemic or topical, in the 3 weeks preceding the study. The study enrolled 11 patients; one patient left due to pathological complications. The study did not involve any procedures on the patients different from those used in normal clinical practice.

DNA extraction

DNA extraction was performed using the QIAamp DNA Microbiome Kit (QIAGEN, Germany) according to the manufacturer’s manual. The kit provides workflow for selective isolation of bacterial DNA from samples that are intrinsically rich in host DNA.

NGS analysis

DNA isolated from 29 human samples was used to generate libraries for templating and sequencing. The samples were as follows: 10 intralesional at t0, 9 perilesional at t0, and 10 intralesional at t20. Microbial DNA was processed for NGS sequencing using the Ion AmpliSeq Microbiome Health Research Kit (Thermo Fisher Scientific, U.S.A.). The assay consists of two pools of oligonucleotide primers and associated reagents to generate amplicon libraries for NGS on Ion Torrent platforms, targeting multiple hypervariable regions in the 16S rRNA gene, as well as 73 species-specific sequences. Two pools of 16 samples were analyzed on two consecutive runs on Ion Torrent S5 using 520 chip cartridges.

Bioinformatic and statistical analysis

Sequencing data were collected, processed, and annotated using the latest available version of Ion Report (5.20.4.0) and the Metagenomics 16 S w.1.1 workflow. Curated MicroSEQ(R) 16S Reference Library v2013.1 and Curated Greengenes v13.5 were used as reference. Bioinformatics analysis on 16S dataset was performed using the method described by Dhariwal et al33. The R package version of Marker Data Profiling (MDP) module was used to collect, inspect, evaluate, and represent NGS dataset. NGS raw data (fastq files) have been deposited at SRA database with BioProject accession code: PRJNA1120631. Plots were generated with RStudio software (v2024.04.1 Build 748) embedded with the latest available packages of R v.4.3.3 (https://www.R-project.org/). Alpha diversity was measured using original counts and applying Chao1, Shannon and Simpson diversity metrics. Post-hoc pairwise comparison (multi-group only) was performed using Welch t-tests for each pair and multi-testing adjustments were based on Benjamini-Hochberg method (FDR: False Discovery Ratio). Beta diversity profiling has been evaluated through non-metric multidimensional scaling (NMDS) ordination of variance stabilized counts of taxa using Bray-Curtis dissimilarity. Permutation analysis of variance (PERMANOVA) and corresponding R-squared and p-values were calculated as reported by Nardelli et al33. We used Linear Discriminant Analysis Effect Size (LEfSe), an algorithm for high-dimensional biomarker discovery. This bioinformatics tool aims at identifying taxa that differentiate between the sample groups using a non-parametric factorial Kruskal-Wallis (KW) sum-rank test identifying features with significant differential abundance relative to the class of interest. Finally, Linear Discriminant Analysis (LEfSe) estimates the effect size of each differentially abundant feature and ranks these features accordingly. Plots were all generated with RStudio software (v2024.04.1 Build 748) embedded with the latest available packages of R v.4.3.3 (https://www.R-project.org/).