This study was approved by the Ethical Committee on Clinical Research of La Paz University Hospital, Madrid (Spain) on June 27th, 2018 (Code: HULP PI-3199). Collection of samples and clinical information was performed after obtaining informed written consent from the parents or legal representatives of the infants participating in this study.

A total of 42 preterm infants were initially recruited at the Service of Neonatology of La Paz University Hospital, Madrid (Spain) from April to October 2018. To detect significant variations in NEFTs bacterial colonization and abundance that may be clinically relevant, we set a threshold of a 1.5 log₁₀ cfu/cm² increase in the attached bacterial population during this period. Initially, a sample size of 33 subjects was estimated, to ensure the reliability and meaningfulness of our results (G*Power version 3.1.9.2, Universität Kiel, Germany; Wilcoxon signed-rank test for matched pairs, α = 0.05, power = 0.95, and effect size = 0.6). Taking into account a 30% dropout rate, the final sample size was increased to 42 subjects. The subjects were born at a gestational age of <32 weeks and admitted to the NICU. Exclusion criteria included preterm infants with any malformations or any congenital metabolic diseases. Immediately after birth, a NEFT was inserted in each infant to provide enteral feeding by gravity eight times a day (every 3 h). On average, each NEFT is inserted for a period of 48 h. The first NEFT (NEFT-1) was removed after 2 days from birth, the second NEFT (NEFT-2) was removed after the first hospitalization week (at day 9) and the third NEFT (NEFT-3) was removed after the second hospitalization week (at day 17). In parallel, samples of meconium (Me) and feces obtained after 2 days (F1), 9 days (F2) and 17 days (F3) of life were also collected from each patient (Supplementary Figure S1).

All samples were collected by the medical and nursing staff of the Neonatology Unit. The fecal samples were collected from the diaper of each infant, using disposable sterile spatulas and sterile single-use gloves, placed in a sterile tube, and stored at −80 °C. NEFTs were removed from the infants, placed into individual sterile plastic bags and frozen at −80 °C. The collected samples were then transferred to the laboratory in containers provided of dry ice, and stored at −80 °C until processing.

DNA extraction from feces was performed using a Maxwell® RSC instrument and a Maxwell® RSC PureFood GMO and Authentication Kit (Promega Corporation, Madison, WI, USA) as described by the manufacturer with the following modifications: 0.1 g of feces were placed into a 2 ml microcentrifuge tube and mixed with 0.8 ml of CTAB buffer (Promega Corporation, Madison, WI, USA) by vortexing for 30 s. Samples were heated at 95 °C for 15 min and allowed to cool for 5 min on the bench top. The content was transferred to a MP Biomedicals Lysing Matrix-A tube (MP Biomedicals, Santa Ana, CA, USA) and mechanically lysed using the MP Biomedicals FastPrep-24 Tissue and Cell Homogenizer (MP Biomedicals, Santa Ana, CA, USA) at a vertical velocity of 6 m/s for 3 cycles of 60 s. The lysate solution was treated with 32 μl of proteinase K (100 μg/ml) and 16 μl of RNase A (4 mg/ml) for 10 min at 70 °C. After incubation, the mixture was centrifuged at 16,000 × g for 5 min, and 300 μl of the supernatant were transferred to the Maxwell RSC Cartridge. The Maxwell RSC Pure Food Protocol was run following the indications of the manufacturer and the DNA was eluted from the samples in 100 μl of the elution buffer. DNA samples were frozen at -20 °C until the metataxonomic analysis was performed.

The internal content of a NEFT was obtained following the next protocol. Collected NEFTs were externally cleaned using 70% (v/v) ethanol. The residual liquid inside the tubes was firstly flushed into a 50 ml sterile tube by pumping air using a syringe. Then, 4 ml of 0.85% (w/v) sterile saline solution were injected into the same tube to remove weakly attached microorganisms. Finally, to fully detach the biofilms developed in the inner surface of the NEFTs, those were cut into 2 cm pieces using sterilized scissors that were collected into a new Falcon test tube containing 10 ml of sterile saline and vortexed for 1 min at maximum speed. Cells were harvested from the resulting suspension by centrifugation at 31,231 × g for 10 min at 4 °C. Pellets were suspended into 1 ml of sterile saline and stored at −80 °C until use. For DNA extraction, 1 ml aliquots were centrifuged at 15,000 × g for 10 min at 4 °C and the pellet was suspended into 400 µl of TE50 buffer (10 mM Tris-HCl: 50 mM EDTA, pH 8.0). The samples were then transferred to Fast-Prep tubes and processed at a Fast-Prep® FP120 equipment (Thermo) for 3 min at maximum speed. Then, tubes were centrifuged at 15,000 × g for 1 min at 4 °C, and the supernatants were transferred to a 2 ml Eppendorf tube. A mix of enzymes was added (120 µl of lysozyme [10 mg/ml], 4 µl of mutanolysine [10 KU/ml] and 3 µl of lysostaphin [4,000 KU/ml]) and blended with the samples. The samples were incubated at 37 °C for 90 min in a thermoblock. After that, 25 µl of Proteinase K (100 μg/ml) and 500 µl of AL Buffer (QIAamp® DNA Mini Kit 250) were added to the samples and incubated again at 56 °C for 30 min. Finally, 100 µl of 3M sodium acetate and 500 µl of 96% ethanol were added to each sample and mixed. DNA purification was performed using a NucleoSpin® Gel and the PCR Clean-up (QIAGEN) protocol. DNA was eluted in 25 µl of molecular grade water (Sigma). DNA samples were frozen at −20 °C until the metataxonomic analysis was performed.

A dual-barcoded 2-step PCR was used to amplify a fragment of the V3-V4 hypervariable region of the bacterial 16S ribosomal RNA gene. Barcodes were attached to 3′and 5′ends of the amplicons to allow separation between forward and reverse sequences. The DNA concentrations of the PCR products were pooled as equimolar and run on an agarose gel. The bands were excised and purified using the commercial QIAEX II Gel Extraction Kit (QIAGEN) and quantified with PicoGreen (BMG Labtech Jena, Germany). The sequencing of the aliquots of the purified barcoded DNA amplicons was performed using the Illumina MiSeq paired-end protocol (Illumina Inc., San Diego, CA) at the facilities of the Scientific Park of Madrid (Spain). Illumina software (version 2.6.2.3) was used to demultiplex the sequences according to the manufacturer's guidelines and pipelines. Bioinformatic performance and analyses were carried out combining QIIME 2 (version 2021.10) (23) and R software (version 3.5.1, https://www.r-project.org/). For denoising, the DADA2 pipeline (24) was used following this set: the forward reads were truncated at position 294 and their first 15 nucleotides were trimmed, while the reverse reads were truncated at position 276 and their first 10 nucleotides were trimmed, to discard positions for which nucleotide median quality was Q19 or below. Taxonomy data was assigned to Amplicon Sequence Variants (ASVs) using the q2-feature-classifier (25) employing classify-sklearn naïve Bayes taxonomy classifier with the SILVA 138.1 as reference database (26). To identify, visualize and remove contaminating DNA, the decontam package version 1.2.1 (27) was applied with a negative control sample. A table of ASVs, genera and phyla sequences per sample was generated, and bacterial taxa abundances were normalized with the Total Sum Scaling (TSS) normalization (relative abundance). After this normalization step, an additional filtering process was implemented where sequences showing a relative abundance of less than 0.01% were excluded from further analysis. This threshold was strategically chosen to focus on the more significant bacterial taxa and to ensure that the study was not impacted by the presence of extremely rare sequences.

For isolates identification, 1 ml aliquot obtained after washing the internal contents of the probes was used. Decimal dilutions were made and 20 μl of each were plated in MacConkey (MCK, Oxoid) for isolation of Enterobacteriaceae. All plates were incubated 24 h at 37 °C under aerobic conditions. At least one representative colony of each morphology was selected and cultured on Brain Heart Infusion (Oxoid) agar for checking purity. Then, selected isolates were identified by Sanger sequencing of the 16S ribosomal RNA (rRNA). Purified PCR products were sequenced (STAB VIDA, Caparica, Portugal) and identified using BlastN (version 2.10) with the 16S rRNA database of NCBI.

All those isolates that were identified as belonging to species of the genera Klebsiella, Serratia and Enterobacter were submitted to genotyping profiling by Random Amplification of Polymorphic DNA [RAPD] procedure. A sample of 5 µl of genomic DNA of each isolate was used in sub-sequent PCR amplifications performed in an iCycler® thermocycler (Bio-Rad, California, USA). PCR conditions are shown in Supplementary Table S1. The primer used for the reaction was OPL5 (5′-ACGCAGGCAC-3′) (Sigma-Aldrich, Schnelldorf, Germany). Biomix Red® (Bioline Meridian, Memphis, USA) was used to perform the reaction. DNA band patterns were obtained after gel electrophoresis [2% (w/v) agarose gel] of the purified RAPD-PCR products (10 µl) obtained using the NucleoSpin® kit Gel and PCR Clean-up (Macherey-Nagel, Düren, Germany). Gels were run using a PowerPac Basic (Bio-Rad) and a Wide mini-sub cell GT (Bio-Rad) for 90 min at 45 V. Gel Red® Nucleic Acid Stain (Biotium, California, USA) (1:100) was used to visualize the DNA bands and the gels were observed in a Bio-Rad Universal Hood II transilluminator (Bio-Rad). Gel images were acquired using Gel Doc 1,000 Documentation System software (Bio-Rad). The band assignment of RAPD profiles was performed using Phoretix 1D Advanced software (Nonlinear Dynamics; Biostep, Jahnsdorf, Germany). Similarity between samples was compared using Phoretix 1D Database (Nonlinear Dynamics; Biostep, Jahnsdorf, Germany). Dendrograms were plotted using the unweighted pair group method with arithmetic averages (UPGMA) (28). Profiles were considered similar when the score values were >0.8.

Quantitative data were expressed as the median and interquartile range (IQR). The Shannon diversity index (29) was used to assess alpha diversity, which accounts for the quantity and evenness of ASVs, via the R vegan package (version: 2.5.6) (30). The Kruskal–Wallis test was performed for determining statistical differences among three or more groups, with Bonferroni p-value adjustment for pair-wise comparisons and Wilcoxon rank test was performed for determining statistical differences among two groups; a p-value < 0.05 was considered statistically significant. Values were expressed as median and quartiles 1 and 3 [Q1-Q3]. Beta diversity was studied using principal coordinates analysis (PCoA) for visually demonstrating patterns of bacterial profiles using a distance matrix filled with dissimilarity values. The Bray–Curtis and binary Jaccard indices were employed for quantitative and qualitative analyses, respectively. Permutational multivariate MANOVA (PERMANOVA) was used with 999 permutations. Heatmaps' cladograms were generated with Hclust hierarchical cluster analysis using the complete linkage method from the R's core “stats” package and the “ggplot2” package (31). Fisher's exact test was employed to evaluate correlations between NEFT and fecal metataxonomic data at a same sampling time. Kendall rank correlation coefficient was used to assess correlations between different genera within a same sample time (either NEFT or fecal samples).

The 16S rRNA-based metataxonomic analysis of the NEFT samples yielded 3,476,951 high-quality filtered sequences and the number of sequences per sample ranged between 13,908 and 83,114, with a median value of 50,432 sequences per sample. Overall, a total of 1,164 ASVs were identified.

When analyzing alpha diversity over time, no statistical differences were detected among NEFT-1 [1.71 (1.38–2.14)], NEFT-2 [1.64 (1.20–2.98)], and NEFT-3 [1.59 (1.18-2.04)] samples (p = 0.56 and p = 0.95 for Shannon and Simpson indexes, respectively) (Figures 2A,B). However, significant differences were observed in relation to beta diversity; on the one hand, the Bray-Curtis distance matrix, which considers the relative abundance of ASVs, showed that NEFT-3 samples were significantly different when compared to those obtained at previous sampling times (p = 0.003). On the other hand, binary Jaccard distance matrix analysis, which considers the presence or absence of ASVs, also revealed significant differences among the different groups of NEFT samples according to the sampling time (p = 0.003) (Figures 2C,D).

Sequences belonging to 71 bacterial genera were detected in the NEFT samples. Among them, the 18 most abundant genera are shown in Figure 2E. NEFT-1 samples were dominated by sequences belonging to the genus Staphylococcus, representing 62.44% (39.79-84.54) of the total bacteria detected at this sampling point, followed by Lactobacillus [0.98% (0.28–3.43)], Streptococcus [0.47% (0.12–4.53)], and Acinetobacter [0.43% (0.16–1.10)]. The relative abundances of these genera in NEFT-2 samples were similar to those found in the NEFT-1 ones. However, the abundance of Staphylococcus sequences was significantly lower (p < 0.05) in NEFT-3 samples [5.52% (0.59–14.84)], while that of some genera belonging to the Family Enterobacteriaceae, including Klebsiella, Enterobacter, and Serratia, increased in comparison to the previous samples. In addition, the relative abundance of minor genera in this type of samples gradually decreased over time.

The metataxonomic analysis of the fecal samples yielded 1,683,527 high-quality filtered sequences, ranging from 4,956 to 38,090 sequences per sample (median value: 22,278 sequences per sample). In this case, a total of 173 ASVs were identified. Alpha diversity analysis did not reveal any statistical difference among Me/F1 [1.37 (1.05–2.06)], F2 [1.34 (1.16–1.74)], and F3 [1.70 (1.35–2.01)] samples (p = 0.18 and p = 0.12 for Shannon and Simpson indexes, respectively), indicating a similar level of species diversity within these groups (Figures 3A,B). However, the Bray-Curtis distance matrix found significant differences between Me/F1 and F3 samples (p = 0.017). While the binary Jaccard distance matrix showed differences in F3 when compared to Me/F1 and F2 samples (p = 0.003 and p = 0.018, respectively). Me/F1 samples also differed statistically with F2 samples (p = 0.003) (Figures 3C,D).

Sequences corresponding to 47 different genera were detected in the fecal samples. Among them, the 18 most abundant bacterial genera in the samples of the 28 patients are represented in Figure 3E. Me/F1 samples were characterized by the presence of sequences of the genera Staphylococcus [14.44% (2.07–50.77)] and, in comparison to other sampling times, by a higher abundance of the minor genera group [0.74% (0.08-8.54)]. In F2 samples, the relative abundance of Staphylococcus spp. decreased up to 5.98% [1.54–21.20]. This downward trend was later confirmed in F3 samples in which Staphylococcus relative abundance significantly (p < 0.05) drops to 0.45% [0.06–3.29]. Parallel, an increase in the relative abundance of different genera of Enterobacteriaceae (Enterobacter, Klebsiella or Escherichia/Shigella) was observed in F3 samples (Figure 3E). Similarly, Clostridium sensu stricto relative abundance tended to increase over time from 0.45% [<0.01–2.66] in Me/F1 to 6.92% [<0.01–20.02] in F3 samples, although no significant differences were observed (p = 0.10).

The assessment of potential correlations between the NEFT and fecal sequences obtained in the same sampling time was carried out excluding NEFT-1 and Me/F1 samples because the DNA concentration in these was extremely low. Such assessment was focused on the most frequent and/or abundant genera found in this study (Staphylococcus, Streptococcus, Klebsiella, Serratia, Enterococcus, and Lactobacillus).

Strong significant associations were observed for the presence of both Staphylococcus (p = 0.001) and Enterococcus (p = 0.002) in NEFT-2 and F2 (Table 2), whilst no significant associations were identified for Streptococcus, Klebsiella, and Lactobacillus (p > 0.05). A non-significant trend for association between these two types of samples was found for Serratia (p = 0.081). When comparing NEFT-2 and F2 samples for potential interactions, significant positive correlations were observed between the relative abundances of Serratia (r = 0.48, p < 0.05), Staphylococcus (r = 0.37, p < 0.05) and Klebsiella (r = 0.33, p < 0.05) (Figure 4). Subsequently, both types of samples were analyzed independently (either NEFTs or feces) (Figure 4). In NEFT-2 samples, a significant negative correlation was found between the relative abundances of Staphylococcus and Klebsiella (r = −0.34, p < 0.05), and between those of Staphylococcus and Serratia (r = −0.14, p < 0.05). In contrast, no significant correlations were observed between these pairs of genera in F2 samples (Figure 4).

A strong significant association was also observed for the presence of Staphylococcus (p = 0.001) in NEFT-3 and F3 samples, while the association regarding Serratia sequences reached a statistical significance, too (p = 0.015) (Table 2). No significant associations were found for Enterococcus (p = 0.117), Lactobacillus (p = 0.128) and Streptococcus (p = 0.346). In this case, a non-significant trend for association between these two types of samples was found for Klebsiella (p = 0.072). A comparison of NEFT-3 and F3 samples for potential interactions, revealed a strong correlation in the relative abundances of Serratia in both sites (r = 0.65, p < 0.05). Positive correlations were also observed between the relative abundances of Klebsiella (r = 0.60, p < 0.05), Staphylococcus (r = 0.37, p < 0.05), and Lactobacillus (r = 0.49, p < 0.05), in both types of samples (Figure 5). Regarding NEFT-3 samples, no significant correlations were observed among the different genera detected in the same sample type. In contrast, significant positive correlations were found between the relative abundances of Lactobacillus and Serratia (r = 0.41, p < 0.05) and, also, between those of Streptococcus and Serratia (r = 0.24, p < 0.05) (Figure 5).

To investigate if bacterial transference from NEFT to feces and/or vice versa might have occurred, the isolates belonging to the genera Klebsiella, Serratia and Enterobacter were submitted to RAPD genotyping. The results are shown in Tables 3, 4 and in Supplementary Figure S2. Globally, 7 patients (28%) shared at least one RAPD-PCR profile in fecal and NEFT samples. The 11% of the profiles were found at least once simultaneously in NEFT and fecal samples from the same patient.

In the case of Klebsiella isolates (n = 67) RAPD-PCR revealed the existence of 9, 4 and 11 different RAPD profiles for Klebsiella pneumoniae (37 isolates), Klebsiella oxytoca (6 isolates) and Klebsiella michiganensis (24 isolates), respectively (Table 3 and Supplementary Table S1). Overall, 9 out of the 24 RAPD profiles (37.5%) were ascribed to a single isolate. RAPD profiles 2, 4 and 6 of K. pneumoniae were shared by 11 (∼ 30%), 7 (∼ 20%) and 5 (∼14%) isolates, respectively. Isolates ascribed to profiles 2 and 6 of K. pneumoniae were recovered from fecal samples of different patients who did not share the same hospitalization box. A similar pattern was observed for profile 9 of K. pneumoniae, profile 4 of K. oxytoca and profiles 1, 3, 6 of K. michiganensis. Similarly, some of the Klebsiella profiles were shared by isolates recovered from NEFTs of patients who were hospitalized in different boxes (profile 4 of K. pneumoniae and profile 2 of K. michiganensis). Profiles 3 and 5 of K. pneumoniae were found in isolates recovered from NEFT and feces samples, belonging to patients sharing the same box. This pattern was also observed for profiles 1 and 4 of K. michiganensis. Overall, 15 out of the 24 Klebsiella profiles (62.5%) were detected at the same sampling time in fecal samples of different patients, and 7 out of 24 (29.2%) persisted over time in the same patient.

The Serratia isolates were identified as Serratia nematodiphilia (24 isolates), Serratia marcescens (11 isolates) and Serratia liquifaciens (3 isolates). These isolates were classified in 8, 3 and 2 different RAPD profiles, respectively (Table 4 and Supplementary Figure S2). Overall, 4 of the 13 genotypes corresponded to a single isolate. Profile 3 of S. marcescens and profile 2 of S. liquifaciens were detected simultaneously in NEFT and fecal samples of the same patient while profiles 1, 2, 3, and 8 of S. nematodiphilia and profile 2 of S. marcescens were detected simultaneously in NEFT and fecal samples of different patients. Seven RAPD profiles (53.8%) were found in samples from different patients at the same sampling time, although only 2 (15.4%) of the Serratia profiles persisted in the same patient over time. Overall, 46.2% of the Serratia profiles were found in patients who did not share hospitalization boxes (Table 4).

The 19 isolates of the genus Enterobacter were identified as Enterobacter hormaechei. Only 3 different RAPD profiles were detected among these isolates and most of them were ascribed to profile 3 (89.4%) (Table 4 and Supplementary Figure S2). Enterobacter isolates sharing this profile were recovered from both NEFT and fecal samples and it was simultaneously detected in almost the totality of patients (5 of 6) either sharing or not the same hospitalization box. In some of these patients (2), this profile was detected over time both in NEFT-2 and NEFT-3, and F2 and F3.