This prospective observational cohort study was conducted from September 2022 to January 2025 at a tertiary care neonatal unit in Northern India. Preterm infants born between 26 and 32 weeks of gestation, with a birth weight of less than 1500 grams, and admitted to the NICU were recruited. Participants’ demographic information and clinical details were collected prospectively until discharge or death using a standardized proforma. NEC was diagnosed according to the Bell staging criteria and classified as cases. Healthy preterm infants matched for gestational and postnatal age were recruited as controls. The study was approved by the Institutional Ethics Committee (IEC) of the Postgraduate Institute of Medical Education and Research (PGIMER), Chandigarh, and was conducted in accordance with the Declaration of Helsinki.

The fecal samples from infants were collected in sterile containers during the first month of life as described in previous studies (Devarajalu et al., 2024; Devarajalu et al., 2025). Fecal samples from control subjects were collected on days of life matched to those of the NEC group. Participants were recruited after obtaining written informed consent from one of the parents. Stool samples were stored in an ultra-low temperature freezer (-80 °C) until further processing.

Genomic DNA was isolated from approximately 100 milligrams of fecal specimens using the QIAamp PowerFecal Pro DNA Kit (Qiagen Inc., Germany). The purity of the extracted DNA was assessed by 1% agarose gel electrophoresis. DNA concentration was determined using the Qubit DNA High Sensitivity (HS) assay (Invitrogen). A paired-end DNA library was prepared using the Twist EF Library Prep Kit (Illumina, Inc., USA). Briefly, the DNA was first fragmented to the desired size, then end-repaired and mono-adenylated at the 3’ end in a single enzymatic reaction. Next, adapters were ligated to the DNA fragments using a T4 DNA ligase-based reaction. Following ligation, the fragments were prepared as substrates for PCR-based indexing in the subsequent step. During PCR, barcodes were incorporated using unique primers for each sample, enabling multiplexing. All prepared libraries were assessed for fragment distribution using a 5300 Fragment Analyzer system (Agilent Technologies, Inc., CA, United States). The resulting libraries were pooled and diluted to achieve optimal loading concentrations. Finally, the pooled libraries were loaded onto a NovaSeq X Plus system (Illumina, Inc., CA, United States) to generate 150 bp paired-end reads. The average read count of the samples is 22.8 million. Sequencing statistics are shown in Supplementary Table 1 .

The adapter sequences were trimmed using the fastq-mcf tool (version 1.04.803). The trimmed reads were then aligned to the human genome (hg19) using BWA (version 0.7.12) to remove human contamination. The remaining unaligned reads were de novo assembled with MEGAHIT (version 1.2.9). The resulting assembled genome was used for open reading frame (ORF) prediction and annotation with Prodigal (version 2.6.3). Taxonomic classification was performed using datasets from the National Center for Biotechnology Information (NCBI), and functional pathways were identified from the predicted ORFs using SEED protein classification in MEGAN 6 (Huson et al., 2016).

Alpha diversity metrics, including observed counts and the Shannon index, were calculated using the Phyloseq package (McMurdie and Holmes, 2013). Differences in alpha diversity between groups were assessed using the Wilcoxon rank-sum test. Beta diversity metrics were computed with the Vegan package. A permutational analysis of variance (PERMANOVA) was performed using the adonis2 function from the pairwise package to identify variation between groups by fitting Bray-Curtis and Jaccard distance matrices separately, incorporating covariates, with 999 permutations (Oksanen et al., 2001).

Differential bacterial abundance between variables was analyzed using Linear Discriminant Analysis (LDA) Effect Size (LEfSe) with an LDA threshold of 4, implemented through the MicrobiomeMarker package (Cao et al., 2022). Additionally, a linear regression model employing the LinDA and MaAsLin2 packages was used to identify taxa and functional pathways associated with necrotizing enterocolitis (NEC) (Mallick et al., 2021; Zhou et al., 2022). Multiple testing corrections, Benjamini-Hochberg (BH) was applied, with a significance threshold of q < 0.05. Welch’s t-test was conducted for differential functional analysis using STAMP software (Parks et al., 2014). Correlation and network analyses were performed using the Network Construction and Comparison for Microbiome Data (NetCoMi) package (Peschel et al., 2021). All statistical analyses were conducted using R software (version 4.4.2), and figures were generated using the ggplot2, pheatmap, and MicrobiomeStat packages.

The baseline demographic information is presented in Table 1 . The NEC and control groups were matched for gestational and postnatal ages. One sample (NEC) was excluded from downstream analysis due to low sequencing read quality.

Microbial richness, assessed using alpha diversity, showed no significant differences in observed operational taxonomic units (OTUs) and Shannon index between infants with NEC and the control group ( Figure 1 ). Alpha diversity indices did not differ significantly when analyzed against other covariates, including mode of delivery (vaginal vs. LSCS), sex (male vs. female), birth weight (<1000 g vs. >1000 g), and gestational age (<28 weeks vs. 28–32 weeks). Furthermore, no statistically significant differences or effect size variations were observed between NEC and control groups for beta diversity, as measured by both Bray-Curtis and Jaccard distance metrics ( Figure 2 ). PERMANOVA analysis revealed that mode of delivery accounted for a greater variance in beta diversity than other covariates within our cohort (Table 2).

The taxonomic composition at the phylum level indicated that Pseudomonadota, Bacillota, and Actinomycetota were the predominant phyla during the first week of life. An increased abundance of Pseudomonadota was observed in preterm infants with NEC, accounting for 47% of the total phyla ( Figure 3a ). The two major families identified in the gut of preterm infants during the first week were Enterobacteriaceae and Enterococcaceae, which together comprised 60% of the total taxonomic families ( Figure 3b ). The abundance of Enterobacteriaceae was significantly higher in the gut of infants with NEC than without NEC, as demonstrated by linear regression analysis ( Table 3 ). There was an increased prevalence of the genera Escherichia and Klebsiella and a decreased prevalence of Enterococcus in preterm infants with NEC. Linear regression analysis indicated a significant increase in the abundance of the genus Klebsiella in NEC cases. Heatmap analysis revealed that either Escherichia or Klebsiella was present in 82% of NEC cases compared to 33.3% of controls ( Figure 4 ). Krona plots showed that the abundance of Klebsiella pneumoniae in the NEC and control groups was 24% and 8% of the total microbial OTUs, respectively, while Escherichia coli accounted for 21% in the NEC cases and 13% in controls ( Figure 5 ).

Differential abundance analysis using LEfSe revealed that Klebsiella pneumoniae and Klebsiella quasipneumoniae were more prevalent in infants with NEC, whereas unclassified Acinetobacter species were more abundant in the non-NEC group ( Figure 6 ).

The relationships between microbial communities at the genus level were evaluated using Pearson’s correlation analysis. A strong positive correlation was observed between Finegoldia and Cutibacterium as well as between Paraburkholderia and Clostridioides. Weak positive correlations were found between Shigella and Klebsiella, and between Escherichia and Klebsiella. Additionally, weak negative correlations were noted between Bifidobacterium and Enterococcus and between Enterobacter and Acinetobacter ( Supplementary Figure 1 ).

Functional analysis using SEED classification revealed that pathways associated with type IV secretion systems, conjugative transfer, Enterobacterial common antigen (LPS O-antigen), quorum sensing in Yersinia, the L-rhamnose pathway, and iron transport systems—including ABC transporters and the Shikimate kinase SK3 cluster—were significantly enriched in infants with NEC compared to controls ( Figures 7 , 8 ). These pathways were significant in both linear regression analysis and Welch’s t-test using the STAMP tool ( Supplementary Figure 2 ). Additionally, genes involved in folate biosynthesis and lactose utilization were reduced in the NEC group compared to the control group.