This study, which was conducted in the Level-III NICU at the First Affiliated Hospital of Nanjing Medical University in Nanjing, China, was approved by the Ethics Committee of the hospital (2022-SR-720) and conducted in accordance with ethical standards as laid down in the 1964 Declaration of Helsinki. All parents provided written informed consent before the study began.

Preterm infants at the gestational age of 24+⁰–31+⁶ weeks were eligible for inclusion. The exclusion criteria included major congenital anomalies, gastrointestinal tract operations, mothers who were alcohol drinkers or smokers, and those who failed to submit informed consent. We enrolled 50 very preterm infants from February 2023 to September 2023. After applying the exclusion criteria, 30 infants remained in the study. The eligible were divided into the BPD group (n = 15) and control group (n = 15), as shown in Figure 1.

The participants were at the gestational age of 24+⁰–31+⁶ weeks and had a birth weight of 520–2,050 g. The antibiotics used for preterm infants were a combination of ampicillin and cefoperazone-sulbactam, according to the treatment practices of the local hospital.

The diagnostic criteria for BPD, as proposed by the National Institute of Child Health and Human Development in 2018, applied to preterm infants with gestational age less than 32 weeks. These infants required more than 3 days of oxygen therapy to keep arterial oxygen saturation between 90 and 95% at a postmenstrual age of 36 weeks. Additionally, these criteria necessitated the presence of persistent parenchymal lung disease confirmed by imaging.

A collection of fecal samples was performed on days 1, 7, 14, 21, and 28, respectively. When the feces were discharged, the researchers collected the samples promptly using a sterile sampling spoon. The samples were subsequently stored in sterile containers at −80°C for 30 min. We collected 139 stool samples and then delivered them on dry ice to BGI Genomics in Wuhan, China, for examination. With four samples that failed to be tested because of low concentration of DNA, the data from 135 fecal samples were eligible for subsequent analyses. We reserved a portion of the fecal samples collected on day 14 for the analysis of SCFAs.

We applied high-throughput sequencing methods to investigate the 16S rDNA of the fecal bacteria. Approximately 0.1 g of feces were placed into a centrifuge tube, and an appropriate amount of lysis buffer (Qiagen, Germantown, MD, USA) was added. A suitable amount of glass beads (Thermo Fisher Scientific, Waltham, MA, USA) was then added to the sample, which was vortexed vigorously for 3–5 min to mechanically disrupt the bacterial cell walls and incubated at 56°C for a duration of 30–60 min thereafter. The samples were heated to 95°C for 10 min and centrifuged at 12,000 rpm for 10 min. The DNA pellets were washed with cold 70% ethanol (Thermo Fisher Scientific). After air drying to remove residual ethanol, the DNA pellets were dissolved in an appropriate volume of nuclease-free water. DNA concentration and purity were measured using a spectrophotometer (Thermo Fisher Scientific) by checking the A260/A280 ratio.

The variable regions V3 and V4 of the 16S rDNA were amplified using polymerase chain reaction (PCR) with the corresponding fusion primers 338F (5′ ACTCCTACGGGAGGCAGCAG 3′) and 806R (5′ GGACTACHVGGGTWTCTAAT 3′). The PCR reaction system was set up using 30 ng of DNA samples and the specified primers. Agencourt AMPure XP magnetic beads (A&D Technology, Tokyo, Japan) were applied to the purification of PCR amplification products to be dissolved in Elution Buffer (Thermo Fisher Scientific). Labels were attached to complete the library preparation. The Agilent 2100 Bioanalyzer (Agilent Technologies, Beijing, China) was used to measure the fragment range and concentration of the library. Based on the size of the inserted fragments, the qualified libraries were then sequenced on an Illumina MiSeq (Illumina, Inc., San Diego, CA, USA) platform.

After the raw data was filtered, high-quality clean reads were obtained. The primers and adapter sequences of the reads matching the primers were removed to obtain the fragments of the target region with Cutadapt (v2.6, Biotronik, Berlin, Germany). When the average value was lower than 20 over a 30 bp sliding window, the terminal sequence of reads was truncated; reads truncated to less than 75% of their original length, as well as those with ambiguous bases or low complexity, were removed.

Thereafter, clean data were prepared for the subsequent analysis. The reads were spliced into tags based on overlapping regions. When the paired-end reads overlapped, a consensus sequence was generated using FLASH (v1.2.11; Biotronik). The tags were clustered to generate operational taxonomic units (OTUs) using USEARCH (v7.0.1090; Qiagen) according to 97% sequence similarity. Chimeras were filtered using UCHIME (v4.2.40; University of California, Santa Cruz, CA, USA).

OTU representative sequences were matched against the database for taxonomic annotation using the RDP classifier (v1.9.1; University of Michigan, Ann Arbor, MI, USA). With the database compared, OTUs were taxonomically classified at phylum, class, order, family, genus, and species levels. We selected the top 10 most abundant species, defined them as key species, and displayed their average relative abundance in each group along with the significance of the differential tests.

Ten samples were selected from the BPD group and ten from the control group. A test was performed on the SCFAs (acetate, butyrate, propionate, valerate, caproate, isobutyrate, and isovalerate) in these samples at the BGI. The feces-derived SCFAs were quantitatively detected by LC–MS/MS using a detection platform equipped with a Waters Iclass-AB Sciex 6,500 LC/MS system (Waters Corporation, Middleton, MA, USA).

A fecal sample (25 mg) was placed into a centrifuge tube, with the addition of 400 uL of methanol-acetonitrile mixture (Thermo Fisher Scientific) and magnetic beads to break it, followed by a 2-min centrifugation performed at 25,000 rpm. For the standard curve, seven types of SCFA mixtures were used for the gradient dilution. To the 20ul sample and standard curve, 60ul cold MeOH/ACN (2:1) (Aslan Chemical Co., Ltd., Shanghai, China) was added, followed by 5 min of shaking and 4 h of precipitation at −20°C. Subsequently, the sample was centrifuged at 20,000 rpm at 4°C for 15 min, and 40 μL of the supernatant was extracted from the EP tube (Eppendorf, Hamburg, Germany). For derivatization, 20 μL of 200 mM 3-NPH (Waters Corporation) and 20 μL of a 120 mM EDC-6% pyridine mixture (Thermo Fisher Scientific) were added to the EP tube. The tubes were then incubated in a metal bath and stirred at 40°C for 30 min. After incubation, the tube was cooled to room temperature on ice before being centrifuged briefly. Finally, 80 μL of 1000D internal standard (solvent: 10% ACN) (Agilent Technologies) was added to the derivatized EP tube and mixed well. To the filter plate were added 90 uL H2O and 90 uL derivatized sample before 5-min centrifugation at 3,000 rpm at 4°C. From a new 96-well plate, 90 uL filtered liquid was collected and 10 uL loaded onto the LC–MS.

We used Statistical Package for the Social Sciences (SPSS version 26; SPSS Inc., Armonk, NY, USA) to analyze the clinical data. Normally distributed numerical variables of the two groups were compared using Student’s independent sample t-test, expressed as mean and standard deviation (SD), while the Wilcoxon rank sum test was performed to compare variables that did not follow a normal distribution, expressed as median and interquartile range (IQR). Counts and percentages were used for categorical variables, and Fisher’s exact test was used to test for differences in the data. A p-value less than 0.05 was deemed statistically significant. A difference analysis was conducted to test differences in the relative abundance of bacteria. The Mann–Whitney U test was used to compare two groups, and the Kruskal–Wallis test was used to compare more than two groups. Alpha diversity was measured using Shannon and Simpson indices. Beta diversity analysis was performed to evaluate the sample differences in species complexity. Additionally, PICRUSt2 (v2.3.0, University of California) was used to predict the functional abundance of the microbial communities based on marker gene sequencing profiles. In the analysis of SCFA differences between the two groups, the ratio analysis and t-test were performed to evaluate them.

The gut microbiota we studied was during the first 28 days after birth; therefore, the clinical data related to the microbiota we collected were mostly during the first 28 days after birth. Subsequently, a comparison was made between the clinical data of the two groups, and the results indicated a lower gestational age in the BPD group than in the control group (p = 0.007; Table 1), but no significant difference in birth weight between the two groups. As shown in Table 1, the BPD group had a higher rate of endotracheal intubation (p = 0.017) and a lower Apgar score at the time point of 5 min (p = 0.008) at birth. Moreover, in the BPD group, a higher rate of invasive ventilation (p = 0.002), longer duration of invasive ventilation (p = 0.002), and longer duration of non-invasive ventilation (p < 0.001) were observed. The BPD group also received a longer administration of antibiotics (p = 0.003) during the 28 days after birth, and it took a longer time for them to initiate enteral nutrition (p = 0.002) and terminate parenteral nutrition (PN) (p < 0.001). Notably, the BPD group had more cases of maternal chorioamnionitis (4 vs. 1) and more infants diagnosed with Ureaplasma urealyticum (UU) infection than the control group (8 vs. 0; p = 0.002).

As indicated in Table 2, the BPD group had more periventricular leukomalacia (PVL), intraventricular hemorrhage (IVH), and hemodynamically significant patent ductus arteriosus (hsPDA) than the control group, although there were no statistically significant differences. Greater retinopathy of prematurity (ROP) (p = 0.042) was observed in the BPD group, along with a higher number of nosocomial infections.

As shown in Figure 2A, the analysis of species accumulation produced a flat curve at the terminal, which indicated adequate sampling quantity. Through analysis of OTUs in the two groups, we observed that the OTUs in BPD group were more abundant on postnatal day 1 but subsequently showed a declining trend, falling below the level of the control after one week of birth. In contrast, the OTU count in the control remained relatively stable (Figure 2B).

At 28 days after birth, the dominant phyla in the BPD group were Bacillota, Pseudomonadota, Mycoplasmatota, Actinomycetota, and Bacteroidota, whereas those in the control group were Bacillota, Pseudomonadota, Actinomycetota, and Bacteroidota (Figure 3A). The BPD group had a higher abundance of Mycoplasmatota (p = 0.013) and a lower abundance of Bacillota than the control group (p = 0.029; Figure 3B). In the 28 days after birth, the dominant genera were Enterobacter, Streptococcus, Enterococcus, Staphylococcus, and Ureaplasma in the BPD group, whereas those in the control group were Streptococcus, Enterococcus, Enterobacter, Escherichia, and Veillonella (Figure 3C). A significant increase in the relative abundance of Enterobacter (p = 0.020) and Ureaplasma (p = 0.010) was observed in the BPD group compared to that in the control group (Figure 3D).

From the intestinal microbiota compared at different time intervals in the BPD group, we found that Bacteroidota showed a continuous declining trend from birth to day 21 and a significant decline from day 14 to day 21 (p = 0.015; Figures 4A,B). Mycoplasmatota was the dominant phylum on day 1 in the BPD group, followed by a rapid decrease, but without a statistical difference. As shown in Figures 4C,D, a transient increase was observed in Veillonella from days 1 to 7, followed by a significant decline from days 7 to 28 (p = 0.023) in the BPD group, the main species being Veillonella dispar. From days 1 to 7, the abundance of Enterobacter showed an increasing trend and then remained high, the main species being Enterobacter hormaechei. In contrast, in the case of Veillonella, there was no significant difference between the five time points, and an increasing trend was observed from days 14 to 28 in the control group (p = 0.051). Enterobacter decreased significantly from days 1 to 7 in the control (p = 0.030), followed by a further decline from days 7 to 14 (p = 0.042). Additionally, the Shannon index decreased with time, but there was an upward trend in the beta diversity in the BPD group (Figures 4E,F).

According to comparisons between the two groups, the BPD group showed a significant increase in Mycoplasmatota (p = 0.027) and a significant decrease in Bacteroidota (p = 0.011) at the phylum level on day 1 (Figure 5A). A significant decrease was also observed in Negativicutes, Betaproteobacteria, Alphaproteobacteria, and Flavobacteriia in the BPD group. The BPD group showed a significant increase in Ureaplasma (p = 0.037; Figure 5B), where the main species was Ureaplasma urealyticum (p = 0.037). The BPD group showed a significant decrease in the Shannon index (p = 0.035) and a significant increase in beta diversity (p = 0.003; Figures 5C,D). The relative abundance of Bifidobacteria and Lactobacillus was extremely low in both groups. On day 1, the abundance of Lactobacillus was significantly lower in the BPD group than in the control group (p = 0.03). Linear discriminant analysis (LDA) was performed; the preset LDA value was 2.0, and an absolute LDA score greater than this value indicated a statistically significant difference. Linear discriminant analysis effect size (LEfSe) showed that Mycoplasmatales, Mycoplasmataceae, and Ureaplasma were the most different species in the BPD group on day 1, whereas Lactobacillales, Stenotrophobacter, and Bacteroidota were the most different species in the control (Figure 6).

On day 7, no significant difference was observed in all levels of key species between the two groups, nor in the Shannon index (p = 0.943) and beta diversity (p = 0.968). As indicated by LEfSe, Lactobacillales was the most different species in the control but had a very low relative abundance (Figure 7).

On day 14, no significant difference was observed in all levels of key species between the two groups, nor in the Shannon index (p = 1.000) and beta diversity (p = 0.106). LEfSe showed that Oscillospiraceae, Lachnospiraceae, and Lactobacillaceae were the most different species in the BPD group, which also had very low relative abundance (Figure 8). On the same day, ten samples were selected from the BPD group, and ten from the control group and a test was performed on seven kinds of SCFAs within the feces, with no difference in SCFAs content between the two groups, as indicated by the PCA analysis (Figure 9A).

On day 21, a lower abundance of Bifidobacterium animalis (p = 0.021) and a higher abundance of Enterococcus faecalis (p = 0.010) were observed in the BPD group. However, no difference was observed in the Shannon index (p = 0.512) or beta diversity (p = 0.096) between the BPD and control groups. LEfSe showed that Macrococcus, Ligilactobacillus, and Succiniclasticum were the most different species in the BPD group but had very low relative abundance (Figure 10).

We simultaneously compared and analyzed the differences in the clinical factors between the two groups on days 7, 14, and 21. Except for a lower gestational age (p = 0.007), a higher rate of endotracheal intubation (p = 0.017), and a lower Apgar score at the time point of 5 min (p = 0.008) at birth in the BPD group; however, there was no significant difference in the duration of antibiotics, feeding mode, duration of PN, and duration of non-invasive ventilation on days 7 and 14, but we observed a longer duration of invasive ventilation on days 7 and 14 (p = 0.047; p = 0.032). We also observed a longer duration of invasive ventilation (p = 0.016) and longer administration of antibiotics (p = 0.013) on day 21.

On day 28, a significant difference was observed in diversity between the two groups. The BPD group exhibited a lower Shannon index, but the difference was not statistically significant (p = 0.060), and there was a significant increase in beta diversity (p = 0.007; Figures 9B,C). The BPD group had a lower abundance of Fusobacteriota, Clostridia, Negativicutes, and Eubacteriales than the control group. We also observed a significant reduction in Veillonella in the BPD group (p = 0.014), with the main species being Veillonella dispar (p = 0.023; Figure 9D). We acknowledge that the relatively small sample size (n = 15 per group) limited the generalizability of our findings. As indicated in Figure 11 (LEfSe), among the different species of the two groups, Veillonella was found to be the most different species. The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that the function abundance of the microbial community had significant differences between the two groups on day 28, when the BPD group had lower function abundance of “glycan biosynthesis and metabolism” (p = 0.028) and higher function abundance of “bacterial infectious diseases” than the control (p = 0.018).

In the analysis of the factors affecting microbiota, we did not observe any statistical difference in the respective rates of breast and mixed feeding and cesarean section between the two groups during the period of antibiotic exposure. When we stratified the preterm infants according to the duration of antibiotic use, the period of antibiotic exposure was 1–27 days before day 28 in the BPD group, which was divided into three subgroups: group A of ≤7 days (n = 4), group B of 8–14 days (n = 5), and group C of >14 days (n = 6). The period of antibiotic exposure was 6–14 days in the control group, which was then divided into two subgroups: group A′ of ≤7 days (n = 8) and group B′ of 8–14 days (n = 7) (Table 3).

Subsequently, we made a comparison of groups A and A′ and observed no significant differences on day 28. We found no significant differences in gestational age (p = 0.052), birth weight, feeding mode, delivery mode, or duration of parental nutrition between the two subgroups. However, we observed significant differences in the duration of invasive ventilation (p = 0.029) and non-invasive ventilation (p = 0.013). When we did the same with groups B and B′ on day 28, intriguingly, we found that the former had a significant decrease in Veillonella and Veillonella dispar (p = 0.009 vs. p = 0.006). Moreover, group B had a lower abundance of Clostridia, Negativicutes, and Eubacteriales. In comparison, the clinical data between groups B and B′ were consistent with those between groups A and A′.