This study involved newborn infants admitted to the Shenzhen Hospital of the Chinese Academy of Sciences from June 2021 to January 2022. Researchers collected demographic data from medical records. They categorized participants into four groups based on gestational age (GA) at birth: PreA group (34–36 weeks), PreB group (32–33 weeks), PreC group (28–31 weeks), and control group (37–42 weeks). Inclusion criteria required that newborns be under three days old, have not received antibiotics or probiotics, and that their legal guardians signed informed consent forms. Exclusion criteria included pregnant mothers who used microecologics or antibiotics within one week prior to delivery, preterm infants who underwent surgeries, and those with congenital malformations, genetic metabolic disorders, or other congenital anomalies. Maternal complications during pregnancy (e.g., preeclampsia, gestational diabetes) were recorded and excluded if they required antibiotic or immunosuppressive therapy. Furthermore, the control group had to demonstrate good health without the use of antibiotics, probiotics, or any other medications. Control group infants were exclusively vaginally delivered to minimize confounding effects of cesarean section on initial microbiota colonization. Researchers obtained written informed consent from each participant’s guardians before enrollment, adhering to the approval granted by the Ethics Committee of the Shenzhen Hospital of the Chinese Academy of Sciences (approval number: LL-KT-21234).

Stool samples were collected from postpartum individuals within three days post-delivery. Nursing personnel collected approximately 0.5–1g of fecal material using sterile spatulas, transferred into pre-weighed 2mL cryotubes, and immediately stored at -80°C. Samples were transported on dry ice to Novogene for 16S rRNA gene sequencing and metabolomics analysis to ensure integrity.

The CTAB method was employed for the extraction of total genomic DNA from the samples. Specific primers (515F-806R for 16S V4) with barcodes were utilized to amplify distinct regions of the genes (16S V3-V4). For PCR reactions, Phusion^® High-Fidelity PCR Master Mix (New England Biolabs) along with forward and reverse primers at a concentration of 2µM and approximately 10ng template DNA were used. Subsequently, TruSeq^® DNA PCR-Free Sample Preparation Kit (Illumina, USA) was employed following the manufacturer’s instructions to generate sequencing libraries with added index codes. The quality assessment of the library was performed using Qubit@2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. Finally, an Illumina NovaSeq platform was utilized to sequence the library resulting in paired-end reads of length 250 bp.

The raw tags were subjected to specific filtering conditions in order to ensure the accuracy of subsequent analysis (Bokulich et al., 2013) utilizing the QIIME software (V1.9.1, http://qiime.org/scripts/split_libraries_fastq.html) (Caporaso et al., 2010). Quality control measures were implemented during this process. Subsequently, the tags were compared against the Silva database (https://www.arb-silva.de/) using usearch software (https://github.com/torognes/vsearch/) (Rognes et al., 2016) for identification and elimination of any chimera sequences (Haas et al., 2011). Sequence analysis was performed using Uparse software (Uparse v7.0.100, http://drive5.com/uparse/) (Edgar, 2013). The filtered sequences were then clustered into operational taxonomic units (OTUs) with a similarity threshold of ≥97%.A representative sequence was selected for each OTU and further annotated accordingly. Alpha diversity analyses were conducted based on this normalized dataset to assess species richness and evenness in a sample using abundance indices (Chao1 and ACE), as well as diversity indices (Shannon and Simpson).Beta diversity, measured by weighted and unweighted unifrac distances, was calculated using QIIME software (Version 1.9.1). Cluster analysis was performed via principal component analysis (PCA) and Principal Coordinate Analysis (PCoA).

The Vanquish UHPLC system (ThermoFisher, Germany) was coupled with an Orbitrap Q ExactiveTM HF-X mass spectrometer (Thermo Fisher, Germany) for sample analysis at Novogene Co., Ltd (Beijing, China). A 100 mg sample was taken and placed in an EP tube, followed by the addition of 500µL of an 80% methanol aqueous solution. The mixture was vortexed and then allowed to sit in an ice bath for 5 minutes. Subsequently, it was centrifuged at 15,000g for 20 minutes at 4°C. A specific volume of the supernatant was diluted with mass spectrometry-grade water to achieve a methanol concentration of 53%. The solution was then centrifuged again at 15,000g for 20 minutes at 4°C, and the supernatant was collected for analysis via LC-MS. Hypesil Gold column (100×2.1 mm, 1.9μm) was utilized for the samples using a linear gradient over a duration of 17 minutes at a flow rate of 0.2 mL/min. In positive polarity mode, eluent A consisted of 0.1% formic acid in water and eluent B comprised methanol. For negative polarity mode, eluent A contained ammonium acetate at a concentration of 5 mM with pH adjusted to 9.0 and eluent B remained as methanol. The solvent gradient was set as follows: 2% B, 1.5 min; 2-85% B, 3 min; 100% B, 10 min;100-2% B, 10.1 min;2% B, 12 min.Q ExactiveTM HF mass spectrometer operated under positive/negative polarity mode with spray voltage of 3.2 kV, capillary temperature of 320°C, sheath gas flow rate of 40 arb and aux gas flow rate of 10 arb, Funnel RF level of 40, Aux gas heater temperature of 350°C.

The UHPLC-MS/MS generated mass spectrometry raw data underwent processing using Compound Discoverer 3.1 (CD3.1, ThermoFisher) to align peaks, select metabolites, and quantify them. Statistical analyses were performed using the statistical software R (version R-3.4.3), Python (version 2.7.6) and CentOS (release 6.6). When the data were not normally distributed, attempts were made to normalize them using the area normalization method. The identification of these metabolites was accomplished by referencing the KEGG database (https://www.genome.jp/kegg/pathway.html), HMDB database (https://hmdb.ca/metabolites), and LIPIDMaps database (http://www.lipidmaps.org/). MetaX was employed for conducting principal components analysis (PCA) and partial least squares discriminant analysis (PLS-DA) (Wen et al., 2017).

Differential metabolites were identified as those exhibiting VIP>1, P value< 0.05, and fold change≥2 or FC ≤ 0.5. Volcano plots were generated using ggplot2 in the R language based on the log2 (FoldChange) and -log10 (p-value) of these metabolites to identify metabolites of interest. Pearson’s method was used in Rlanguage to analyze the correlation between these differential metabolites. The KEGG database was utilized for investigating functions and metabolic pathways associated with these metabolites, while metabolic pathway enrichment analysis was conducted by assessing ratios x/n>y/N for pathway enrichment and considering a metabolic pathway statistically significant if its P value <0.05.

The statistical analyses were conducted using SPSS 24.0 software (SPSS, Inc., Chicago, IL, USA) and R software (Version 2.15.3). Demographic characteristics and delivery data were compared between groups using Student’s t-test, Wilcoxon rank-sum test or Mann-Whitney U test. Differences in categorical variables were determined using the χ² or Fisher’s exact test. A principal components analysis was employed to compare samples among groups and identify statistically significant species at each classification level. The comparison of four groups for alpha diversity was assessed by the Wilcoxon rank-sum test. Unique taxa or functions that distinguished the groups as statistically significant biomarkers were identified through linear discriminant analysis effect size (LEfSE) with an LDA cutoff of 4. Associations between gut microbiome composition and metabolites were demonstrated using Spearman correlation coefficient. A significance level of P < 0.05 was deemed appropriate.

A total of 48 subjects meeting the inclusion criteria were enrolled at Shenzhen Hospital of the Chinese Academy of Sciences between June 2021 and January 2022.Stool samples were collected from neonates with a median birth weight of 2.36 kg (range: 1.1-3.7 kg) and GA ranging from 28 to 41 weeks ( Supplementary Table 1 ). Supplementary Table 1 provides an overview of the baseline characteristics for each group. The mean GA of the PreA, PreB, PreC and control groups were 35.89, 33.00, 30.12 and 39.64 weeks, respectively. The 48 newborn infants were divided into the control group (n=15), the PreA group (n =13), the PreB group (n=11) and the PreC group (n=9).

The alpha diversity, which quantifies the species richness and taxonomic diversity within a sample, was evaluated using multiple indices (ACE index, Chao index, Simpson index, and Shannon index).Our analysis revealed significant differences in ACE and Chao indices between the PreA group and control group (P=0.001 for Figure 1A ; P=0.000 for Figure 1B ), indicating that gut microbiota richness was higher in the PreA group than controls. Furthermore, the Wilcoxon rank-sum test demonstrated that both ACE and Chao indices were significantly elevated in the PreC group compared to controls (P=0.046 for Figure 1A ; P=0.023 for Figure 1B ). The Shannon index of PreB samples was lower than that of PreC samples (P=0.042 for Figure 1D ), while the Simpson index of PreB samples was lower than all other groups ( Figure 1C ). Overall, our findings suggest a tendency towards greater alpha diversity of fecal microbiota in the PreA group relative to controls but no significant difference between the PreA and PreC groups based on these four indices.

The PCoA analysis revealed a distinct separation among the samples into four groups, with a significant distinction observed between the PreC group and the other three groups ( Figure 1E ). NMDS plots also indicated a significant difference in bacterial communities between the PreC group and the remaining groups (stress=0.14, Figure 1F ). ANOSIM results demonstrated a significant dissimilarity in beta diversity between the PreB group and PreC group (R=0.118, P<0.05).Consistent with these findings, Amova analysis further confirmed an extremely significant variation in gut microbiota composition between the PreB/control group and PreC/control group (P<0.05). Our data suggested that as GA increases, preterm infants exhibit a gradual convergence of gut microbiota towards term infants. Notably, preterm infants born before 32 weeks displayed significantly different gut microbiota diversity compared to other gestational ages, which may contribute to subsequent digestive complications. Therefore, future research should prioritize investigating preterm infants born before 32 weeks of gestation.

The most predominant phyla identified in the four groups were Proteobacteria and Firmicutes ( Figure 2A ). The control group exhibited the lowest relative abundance of Proteobacteria, while Firmicutes were found to be the most abundant. Additionally, it was observed that the PreC group tended to have a higher proportion of Bacteroidota and Actinobacteria, as well as a lower proportion of Firmicutes compared to the other three groups; however, no statistically significant differences were detected (P>0.05). Furthermore, an increase in Desulfobacterota and a decrease in Spirochaetota at the phylum level were observed in the control group when compared to the PreB group (P<0.05, Figure 2D ).

At the family level, a significant decrease in Peptostreptococcaceae was observed in both the PreB and PreC groups compared to the control group (P<0.01, Figure 2B and Figure 2E ). The analysis of flora composition at the family level revealed that Staphylococcaceae was more abundant in the PreA group than in the other three groups. Conversely, Moraxellaceae and Pseudomonadaceae showed decreased abundances compared to the other three groups; however, these differences were not statistically significant ( Figure 2E ). Notably, Enterococcaceae was exclusively present in the PreB group while Enterobacteriaceae and Prevotellaceae were found only in the PreC group. Burkholderiaceae and Pseudomonadaceae were enriched in the control group. As gestational age increases, Peptostreptococcaceae gradually emerged, while Prevotellaceae and Enterobacteriaceae progressively diminished, Moraxellaceae remained constant.

The most prevalent genera identified across all four groups were Acinetobacter, Enterococcus, Staphylococcus and Paeniclostridium ( Figure 2C ). In comparison to the other three groups, the control group demonstrated an elevation in Paeniclostridium, Megamonas and Prevotella at the genus level while Buchnera decreased ( Figure 2F ). Furthermore, a higher abundance of Escherichia-Shigella and Ureaplasma was observed in the PreC group ( Figure 2F ). The PreAgroup exhibited significantly increased relative abundances of Bifidobacterium and Streptococcus than the control group (P<0.05, Figure 2F ). As gestational age increases, Megamonas and Prevotella gradually emerged, while Escherichia-Shigella and Ureaplasma progressively diminished.

The dominant floras were identified by analyzing all the samples using LEfSe (LDA>4). Compared to the other three groups, the control group showed a significant increase in the abundance of Peptostreptococcaceae (family) and Peptostreptococcales-Tissierellales (order) ( Figure 2G ). The LEfSe test indicated that Methylobacterium_brachiatum(species) had a significantly higher relative abundance in the PreA group compared to the control group ( Supplementary Figure 1 ). Among the dominant fecal floras, it was observed that Peptostreptococcaceae (family), Peptostreptococcales-Tissierellales (order), and Vibrio (genus) were more prevalent in the control group than in PreB group ( Supplementary Figure 1 ). Additionally, through LEfSe analysis, it was noted that Peptostreptococcaceae (family) showed greater abundance in the control group while Methylobacterium_brachiatum (species) displayed higher prevalence in the PreC group ( Supplementary Figure 1 ).

The heatmap analysis revealed disruptions in KEGG pathways at different levels (level 1, level 2, level 3, and KO) in preterm infants compared to the control group ( Figure 3 ). The analysis of KEGG pathways indicated abnormalities in various aspects such as organismal systems energy, genetic information processing, cellular processes, metabolism, and human diseases in preterm infants. For instance, both organismal systems energy and genetic information processing were found to be higher in the PreC group than in the control group ( Figure 3A ). Notably, glycan biosynthesis and metabolism as well as transport and catabolism exhibited significant enrichment among the KEGG level 2 pathways of the PreC group; conversely, cancers and neurodegenerative diseases showed significant enrichment in the PreB group ( Figure 3B ). Our study demonstrated the enrichment for transporters, ABC-transporters, General_function_prediction_only and Transcription_factors within the control group’s KEGG level 3 pathways when compared to preterm infants ( Figure 3C ). Additionally, we observed a significantly higher abundance of unique genes such as RNA polymerase sigma-70 factor (K03088) and branched-chain amino acid transport system substrate-binding protein (K01999) in the PreC group compared to controls; however phosphoglycolate phosphatase (K01091) was lower in abundance with statistical significance (P<0.05, Figure 3D ).

The differential untargeted metabolomic substances in each group were analyzed using principal component analysis (PCA). Notably, a clear separation was observed between the PreC group and the other three groups, as illustrated in Figure 4A . The heatmaps exhibited distinct differences in intestinal metabolites among these groups ( Figure 4B ). Volcano maps revealed distinct screening results for positive ions, leading to the identification of 144 metabolites that exhibited differential abundance between the PreA group and control group. Among these metabolites, 18 were upregulated while 126 were downregulated in the PreA group compared to controls ( Figure 4C ). Specifically, a total of 229 differential metabolites were identified in the PreB group with 31 being upregulated and 198 downregulated relative to controls. Moreover, there were a total of 225 differential metabolites detected in the PreC group with 30 being upregulated and 195 downregulated when compared to controls ( Figure 4D ). The metabolites in PreB and PreC exhibited more pronounced disparities in their metabolic profiles as compared to term infants. The disparity in metabolite levels increases with the discrepancy in GA. Preterm infants displayed diverse degrees of metabolic downregulation compared to full-term infants.

The identified differential metabolites were classified and annotated using the KEGG database to gain further insights into the functions performed by these specific metabolites. Our findings indicated that these differential metabolites were primarily involved in cellular processes, environmental information processing, metabolism, and genetic information processing ( Figure 5A ). Notably, significant differences were observed in the metabolite profiles between preterm infants groups compared to the control group ( Figure 5 ). Through KEGG enrichment analysis, we discovered that pyrimidine metabolism was a distinct pathway associated with differential abundance of metabolites between PreA group and control group ( Figure 5B ). Similarly, only one main differential metabolic pathway was found between PreB group and control group-arachidonic acid metabolism ( Figure 5C ). The bubble diagram revealed significant enrichment of differential metabolites in pyrimidine metabolism and beta-alanine metabolism pathways when comparing PreC group to control group ( Figure 5D ). As gestational age increases, the metabolism of arachidonic acid progressively intensified, while pyrimidine metabolism gradually diminished, and the metabolism of secondary bile acids remained stable.

To investigate the relationship between fecal metabolites and gut microbiota in different groups, we performed a Spearman correlation analysis at the genus level. The results revealed significant correlations between Zoogloea, Luteolibacter, Gallicola and Pseudoxanthomonas with various metabolites in the PreA vs control group ( Figure 6A ). Moreover, Marinospirillum, Romboutsia, Propionicicella, Taibaiella and Pseudoalteromonas were found to exhibit significant correlations with different metabolites in the PreB vs control group. Furthermore, Psychrobacter, Vagococcus, Anaerostipes, Marinospirillum and Lachnospira demonstrated significant associations with distinct metabolites in the PreC vs control group ( Figure 6C ). There was a strong correlation (R=0.79) observed between Anaerostipes and N-Phenylacetylglutamine as depicted by Figure 6C .

The levels of control-enriched N-Phenylacetylglutamine, 3-Hydroxy-3-methylbutanoic acid, and cis,cis-Muconic acid were positively correlated with the abundance of Agathobacter, as depicted in Figure 6C .Furthermore, Romboutsia exhibited a positive association (R=0.585) with the level of control-enriched Argininosuccinic acid ( Figure 6B ). In Figure 5A , it can be observed that there is a positive correlation between Acidaminobacter abundance and control-enriched Nicotinamide levels (R=0.43). These findings suggested a close connection between the distinctive metabolites found in preterm infants and changes in gut microbiota composition. For instance, Luteolibacter displayed a positive correlation with Nicotinamide but demonstrated a negative correlation with RNK and Deoxyribose 5-Phosphate in Figure 6A .