From August 2018 to October 2020, 44 VLBW preterm infants were signed in for this study. We excluded 10 infants who had no stool passage within the first postnatal week, one infant who missed the third stool collection time, and one who died of pneumothorax with pulmonary hypertension after second stool collection. Finally, 32 VLBW infants and 10 healthy term infants were enrolled. Ten preterm infants developed 12 episodes of invasive bacterial infections, with their fecal samples being classified as the preterm invasive group (PI). The other fecal samples from 22 preterm infants without invasive infections were designated as the preterm non-invasive group (PNI). Stools collected from term infants were classified as the healthy term group (HT). The clinical information is summarized in Table 1. The mean gestational age of the preterm infants in the PI group was younger, relative to the preterms in the PNI group (26.9 ± 2.2 vs. 28.7 ± 2, p = 0.027), but the average birth weight was similar between the two groups (1021.6 ± 283.4 vs. 1155.8 ± 161.4, p = 0.188). Infants with invasive bacterial infections also had less cumulative human milk, formula milk, and total feeding consumptions during the first week, but only formula milk was statistically significant (Table 1). The duration of empiric antibiotic use and proportion of probiotic use in the first week was similar between the PI and PNI groups; the DOTs until the third stool collection time (7 weeks) were also not significantly different between the two groups (Table 1). The duration of hospitalization was longer for the PI group, but the difference was not statistically significant between the two groups (Table 1).

The bacterial isolates of sepsis in these preterm infants are listed in Table 2. Nine had bloodstream infection, and three had central nervous system (CNS) infection. These infections were caused by Gram-negative bacilli (GNB) in three (25.0%), Gram-positive bacilli (GPB) in one (8.3%), and Gram-positive cocci (GPC) in eight (67.7%). All GPC were Staphylococcus species. Patient C06 experienced two invasive infections: one by Staphylococcus capitis at 27 days of age and the other by carbapenem-resistant Serratia liquefaciens at 48 days of age. Another patient C38 had hydrocephalus with ventriculo-peritoneal shunt (VPS) implantation and developed two distinct invasive CNS infections: one by Staphylococcus epidermidis sequence type (ST)35 at 95 days of age and the other by S. epidermidis ST59 at 116 days of age. The MLST analysis revealed that infection in patient C03 was also caused by S. epidermidis ST35, and infection in patient C28 was caused by S. epidermidis ST59. Based on the time of onset, the 12 episodes of invasive infections included 1 EOS, 8 LOS, and 3 VLOS.

A total of 20 stool samples collected from healthy term infants and 96 from VLBW infants were subject to analysis. Figure 1 shows α- and β-diversity among the three groups. According to the Shannon index, the infants with invasive infection (PI group) had significant lower α-diversity (p < 0.05) since birth and persisted to 4 and 7 weeks of age (Figure 1B). According to the Chao1 and Shannon index, the dysbiosis persisted to 7 weeks of life (Figures 1A,B). Preterm infants without invasive infections (PNI group) showed higher microbial diversity in the gut at birth, but the diversity decreased significantly in the following weeks (p < 0.001). Healthy terms (HT group) maintained higher diversity at birth and thereafter (Figures 1A,B). The analysis of β-diversity (covariate adjusted principal coordinates analysis) showed significant variations among each group (p < 0.001) (Figure 1C). In terms of the dynamic change overtime, we found that microbiota in the HT group remained stable from birth to 4 weeks of age (Figure 1C, two circles in the upper part), but that of the PI and PNI groups changed with age from birth to 7 weeks of age (Figure 1C, the representative circles moved from the left lower part to the right lower part with time).

Figures 2, 3 show the taxonomy profile and the proportion of major phyla and genera in heatmap, respectively. The phyla of gut microbiota in healthy term infants after birth (HT1) were similar to those in preterm infants after birth (PNI1 and PI1), in which Proteobacteria accounted for 40.3–48.8%, followed by Firmicutes (23.9–29.5%), Bacteroidetes (10.9–28.4%), and Actinobacteria (5.2–8.5%). At 4–7 weeks of age, however, the proportion of Proteobacteria decreased to 16.3% in HT2 but increased to 72.8 and 80.3% in PNI3 and PI3, respectively. Meanwhile, the proportion of Firmicutes increased to 42.6% in HT2 but decreased to 17.9 and 15.7% in PNI3 and PI3, respectively. The trends indicated that gut microbiota evolved toward dysbiosis in preterm infants, especially in the PI group (Shin et al., 2015). In the level of genera, the proportion of GNB (Enterobacter, Escherichia/Shigella) increased over time among preterm infants, and the magnitude was higher in the PI group than in the PNI group. The organisms with significant differences in the three groups are shown by cladogram in Figure 3. Bacteroidetes was the most significant phylum in HT1, Actinobacteria in PNI1, and Proteobacteria in PI1.

All preterm infants in our NICU had antibiotic exposure during hospitalization. The median DOT before the third stool collection of 32 preterm infants was 20.5 days (range, 4–77 days). The fecal samples collected from preterm infants whose DOT was higher than the median were classified into the DOT_>median group and the other into the DOT_<median group. Fecal samples of healthy term infants without antibiotic exposure collected at the fourth week of age (HT2) were used as control. The rates of invasive infections were not significantly different between infants of the DOT_<median group and the DOT_>median group (25.0 vs. 37.5%, p = 0.70).

The microbiome of fecal samples was analyzed to assess the impact of antibiotic use on gut microbiota. The α-diversity (Chao1 and Shannon index) was significantly lower in preterm infants (DOT_<median and DOT_>median) than in healthy term infants (HT2) (Figure 4). Analysis of β-diversity also showed significantly different microbiota composition in the gut among HT2, DOT_<median, and DOT_>median groups (p < 0.001) (Figure 4). The difference of the gut microbiota between healthy terms and preterm infants according to antibiotic exposure is shown in Figure 5 and Supplementary Figure 1. The major phyla in healthy terms were Firmicutes (42.6%), Actinobacteria (26.8%), Proteobacteria (16.3%), and Bacteroidetes (14.3%); however, the major phyla in both DOT_<median and DOT_>median were Proteobacteria (71.7 and 78.3%, respectively) and Firmicutes (18.7 and 15.7%, respectively), while the proportion of Bacteroidetes was low in both (Figure 5). In the genus level, the top five predominant genera in healthy terms were Bifidobacterium (26.5%), Veillonella (15.6%), Bacteroides (11.2%), Escherichia/Shigella (10.2%), and Streptococcus (9.7%), while in DOT_{< median} and DOT_>median, the predominant genera were Enterobacter (46.3 and 46.7%, respectively), Escherichia/Shigella (17.7 and 24.4%, respectively), and Bifidobacterium (9.5 and 5.7%, respectively). Gut microbiome in the DOT_>median group was composed of significantly more bacteria like genera Enterobacter and Pseudomonas than other groups (Supplementary Figure 1). The genus Bifidobacterium accounted for 26.5, 9.5, and 5.7% in healthy terms, DOT_<median, and DOT_>median, respectively. The genus Lactobacillus accounted for 7.4, 2.5, and <0.1% in HT2, DOT_<median, and DOT_>median, respectively. Thus, prolonged exposure to antibiotics did not reduce the proportion of bacteria in the genera Enterobacter and Pseudomonas but reduced beneficial bacteria, such as Bifidobacterium and Lactobacillus, in the gut.

The proportions of the bacterial species identical to the causative pathogens for sepsis in gut microbiota are summarized in Table 2 and Supplementary Figures 2–5, and the results of stool PCR are shown in Table 3. Bacterial species identical to the causative pathogens has been predominant in the gut microbiota of three premature infants with sepsis (C08, C28, and C44). Patient C08 developed LOS by Klebsiella pneumoniae between the second and third stool collection, and consistently, the proportion of Klebsiella genus increased from 43.9% in the first stool sample to 82.4 and 93.3% in the second and third stool samples, respectively (Table 2 and Supplementary Figure 2, Klebsiella). The pathogenic K. pneumoniae was only detected by PCR in the second and third stool samples, suggesting that K. pneumoniae being present in the first sample was likely not pathogenic (Table 3). The mother of patient C44 had chorioamnionitis by Escherichia coli (data not shown), and patient C44 developed EOS by E. coli before the first stool collection. The microbiota showed that the Escherichia genus accounted for 99.8, 1.8, and 44.9% in the first, second, and third stool samples, respectively (Table 2 and Supplementary Figure 3, Escherichia). The pathogenic E. coli was indeed detected in the first fecal sample. After antibiotic treatment, the E. coli disappeared in the second stool sample; however, it reappeared in the third stool sample (Table 3). Patient C22 had external ventricular drainage (EVD) for post-hemorrhagic hydrocephalus and developed LOS and EVD shunt infection by S. aureus. The microbiota in the genus level showed that Staphylococcus only accounted for 0.3, 0.3, and 2.4% in the three stool samples (Table 2 and Supplementary Figure 4, Staphylococcus). The PCR confirmed that pathogenic S. aureus was present in the second and third fecal samples (Table 3). Patient C28 developed LOS by S. epidermidis after third stool collection. Staphylococcus in gut microbial composition was 93.4, 1.1, and 0.1% in the three fecal samples (Table 2 and Supplementary Figure 5, Staphylococcus); however, PCR confirmed the presence of pathogenic S. epidermidis in the second stool collection (Table 3).

We designed strain-specific primers using genome sequences of the bacterial isolates and performed PCR upon the stool samples to track the dynamic change of the pathogens. The results clearly distinguished the pathogenic clone from the colonizing bacteria of the same species. Colonization of pathogenic bacteria in the gut of patients C08, C22, and C28 with LOS and of patient C44 with EOS was confirmed (Table 3). During the study period, there was an outbreak of B. cereus sepsis in NICU. B. cereus was confirmed from contaminated linen (data not shown). The present study found no B. cereus in the stool samples of patient C42, indicating that B. cereus was derived from the environment (Table 3).

PCR with strain-specific primers also examined whether the pathogenic bacteria of the same genotype in different patients were identical strains or not. MLST revealed that S. epidermidis ST35 contributed to infections in both patients C03 and C38 (first infection), and S. epidermidis ST59 to those in patients C28 and C38 (second infection). The results of PCR indicated that S. epidermidis ST59 from patients C28 and C38 (second CNS infection) who had hydrocephalus with VPS infection were identical because both were identified by the same strain-specific primers (Table 3). Notably, the two infections occurred nearly 5 months apart, and S. epidermidis ST59 was not detected in any fecal sample of patient C38. Furthermore, PCR confirmed that S. epidermidis ST35 from patients C03 and C38 (first CNS infection) belonged to different bacterial clones because strain-specific primers only recognized the bacterial isolates these primers were derived from (data not shown).

This prospective cohort study was approved by the Institutional Review Board of Chang Gung Memorial Hospital (201702153B0). Preterm infants with birth weight between 500 and 1,500 g and no congenital gastrointestinal anomaly were enrolled after obtaining consent within 72 h of age between August 2018 and October 2019. Fecal samples in the first, fourth, and seventh weeks of life were collected. All participants were followed up for 6 months. Those who had clinical symptoms of infection and bacteria yielded from sterile sites, such as blood or cerebral spinal fluid were defined as having invasive infections or sepsis (Shane et al., 2017). Clinical information about pregnancy complications, perinatal conditions, and any postnatal infections was also recorded. We also collected stool samples from healthy term infants as controls within the first and fourth weeks of age during the study period. All stool samples were analyzed with 16S rRNA gene sequencing to investigate the intestinal flora. All invasive bacterial isolates were subjected to genome sequencing. According to the genomic sequences, strain-specific primers were designed for each bacterial isolate, and PCR was done to check whether the invasive bacteria derive from the gut. In addition, the amount of antibiotics for each patient was measured by days of therapy (DOTs) and was calculated as the sum of the total number of days of all used antibiotics (Momattin et al., 2018).

Stool samples were collected from diapers using the Longsee Fecalpro Kit (Longsee Medical Technology Co., China). Genomic DNA extraction was performed using the QIAamp PowerFecal DNA Kit (Qiagen, United States) according to the instructions of the manufacturer. About 250 μg of the sample was added to 750 μl of PowerBead Solution and 60 μl of Solution C1, and the tube was heated at 65°C for 10 min. Then, the mixture was vortexed using a PowerLyser Homogenizer for 10 min at 1,000 RPM. The DNA was extracted and washed, eluted using nuclease-free water according to the QIAamp PowerFecal DNA Kit handbook. The extracted DNA was stored at -80°C. The concentrations and qualities of the purified DNA were determined by Qubit 4 Fluorometer with high-sensitivity dsDNA assay (Thermo Fisher Scientific).

The genomes of the cultured bacteria from patients with invasive infections were sequenced with Illumina Miseq platform (Illumina, San Diego, CA, United States). The short reads generated were de novo assembled into contigs using SPAdes version 3.11.1 (Bankevich et al., 2012). The assembled genome sequences have been submitted to the NCBI GenBank database under the accession number PRJNA674911.¹ The multilocus sequencing typing (MLST) was performed by the online database BacWGSTdb (Ruan and Feng, 2016). For each isolate, the genome was compared against the representative genomes in BacWGSTdb through the Cross_match program. The derived sequence-type specific fragments were used for designing specific primers for subsequent PCR confirmation.

The MLST was performed by using the protocol in the online database BacWGSTdb (Ruan and Feng, 2016). For each isolate, the genome was compared against the representative genomes in BacWGSTdb through the Cross_match program. The derived sequence-type specific fragments were used to design specific primers for subsequent PCR confirmation (Supplementary Table 1). For each strain, two pairs of primers targeting two genes were designed. To prove the presence of cultured pathogenic bacterial strains in feces, two unique primers to recognize the specific sequence of pathogenic bacterial strains were designed for each cultured pathogenic bacterium. Each fecal sample was processed for PCR to confirm the presence of the causative bacteria in the stool.

The procedures for primer designing include the following: (1) choose a DNA fragment that does not come from a general plasmid region, determined by means of NCBI BLASTX with a database of non-redundant protein sequence and with a specific organism at https://blast.ncbi.nlm.nih.gov/Blast.cgi; (2) select a sequence (150–1,500 bps) of the chosen DNA fragment that has a unique single nucleotide polymorphism aligned with other sequences of sources at both ends (Ruan and Feng, 2016). The two end sequences are designed as a set of PCR primers with a size of 20--25 bps and 40--60% GC content by means of Nucleotide BLAST² or Primer 3.³ Primers used in this study are listed in Supplementary Table 1. Each fecal sample was processed for PCR with unique primers to find pathogenic bacteria strains. Total RNA was reversely transcribed, and the subsequent cDNA samples were subjected to real-time PCR using SYBR Green (Vazyme, China) and the primer sets.

The bacterial 16S rRNA gene amplicon was constructed based on the hypervariable V3–V4 region (Klindworth et al., 2013) and sequenced using the MiSeq System (Illumina, United States). The 16S rRNA gene amplicon library sequencing reads were initially demultiplexed using MiSeq Reporter v2.6 according to sample barcodes. Raw reads were performed using Usearch (V.11)⁴ for pairing reads, quality filtering, and operational taxonomic unit (OTU) clustering. The RDP training set (v16) database was used as a reference set for final taxonomic assignments performed using the SINTAX classifier. The α-diversity (e.g., Chao index and Shannon index) and β-diversity measurements were calculated. Finally, we used linear discriminant analysis (LDA) of effect size (LEfSe) to determine the taxonomy that most likely explains the differences between groups (Segata et al., 2011).

The clinical characteristics were analyzed using Student’s t-test and chi-square test. The relative abundance of bacteria and alpha diversity indices were compared and displayed by GraphPad Prism 8 (GraphPad Software, Inc., La Jolla, CA, United States). Between-group comparisons were performed with a random forest for two groups. We performed multivariate logistic regression analysis on the bacteria markers and used the statistical logistic forward method to produce a panel to identify case patients. Statistical analyses were performed using the IBM SPSS Statistics software, version 25.0 (IBM Inc., Chicago, IL, United States). The tests were two-sided, and a p < 0.05 was considered statistically significant.