Abstract
Background
The variant microbiome is associated with necrotizing enterocolitis (NEC). We aimed to analyze remnant formalin fixed paraffin-embedded (FFPE) intestine tissue for microbiome profiling in preterm infants with surgical NEC.
Methods
We analyzed FFPE small intestine tissues from 16 infants with NEC (8 survivors and 8 non-survivors). Extracted DNA from FFPE tissue blocks underwent 16S rRNA sequencing. We compared the microbiota profiles in survivors and non-survivors. Alpha- and beta diversity metrics were calculated using QIIME2. To assess differences in overall microbial community structure, we performed a Permutational Multivariate Analysis of Variance (PERMANOVA). The analysis was performed in MaAsLin2 R package to determine the specific microbial taxa whose relative abundances were significantly associated with survival status in a multivariable linear model.
Results
Sequencing of FFPE extracted DNA resulted in high-quality sequence reads in 16 cases.
Analysis
Analysis of microbial communities from 16 cases revealed a significant association between microbiome structure and survival status. Beta diversity analysis demonstrated distinct clustering of microbiome profiles between survivor and non-survivor groups. Alpha diversity metrics further characterized these differences: the non-survivor group exhibited a more complex and even microbiota (Shannon entropy, p = 0.02; Pielou’s evenness, p = 0.017), whereas the survivor group’s microbiome was significantly richer in observed features (p = 0.004). Notably, this association was specific to survival outcome, as overall community structure did not significantly differ when grouped by histological features of disease severity such as necrosis, inflammation, or hemorrhage. Linear mixed effect models and direct comparisons further identified numerous taxa potentially associated with survival status.
Conclusion
FFPE intestinal tissue enabled retrospective and spatially relevant microbiota assessment at the disease site. The non-survivors had complex microbiota and had distinct bacterial communities compared to survivors. Our findings suggest that the gut microbiome is a key factor related to prognosis, independent of other measures of NEC severity.
Introduction
Necrotizing enterocolitis (NEC) is the most common acute gastrointestinal illness during the neonatal period, affecting about 5%–10% of preterm infants with a birth weight of less than 1500 grams.1,2 NEC, especially surgical type, remains a leading cause of morbidity and mortality among premature neonates and leads to increased cost of care and resource utilization.3–9 The etiology of NEC is still unknown and is likely multifactorial. Many studies report the role of bacteria, viruses, and fungi of the gut microbiome in the pathogenesis of NEC and late-onset sepsis.10,11
Our previous study demonstrates that NEC-associated sepsis was present in one out of four cases in infants with medical and surgical NEC. 12 Gram-positive bacteria, gram-negative bacteria, and Candida were cultured in blood in 15.3%, 10.5%, and 2.8% of NEC cases. A recent study reported higher mortality in infants with bloodstream infections within 3 days following NEC diagnosis. 13 A recent metanalysis by Pammi et al14,15 reported that fecal microbiome from preterm infants with NEC had increased relative abundances of Proteobacteria and decreased relative abundances of Firmicutes and Bacteroidetes prior to NEC onset. Existing studies are largely limited to non-invasive stool samples, which may not be reflective of the bacteria located at anatomical sites of disease. The gut dysbiosis seen in preterm infants with NEC is associated with adverse outcomes. 16 In this study, we aim to determine and compare the microbiota profiles in survivors and non-survivors by using the historical formalin fixed paraffin-embedded (FFPE) tissue from NEC preterm infants. The FFPE intestinal tissue analysis will reflect the microbes translocating from the lumen into the lamina propria and submucosal tissues and not that in the lumen of the GI which is an important distinction perhaps from stool sample analyzing microbiome studies.
Methods
This single-center, retrospective observational cohort study was conducted at the University of Mississippi Medical Center (UMMC) at Jackson, Mississippi and the Wake Forest School of Medicine, Winston Salem, North Carolina, after the Institutional Review Board (2017-0127) and material transfer agreement approval. UMMC’s level 4 neonatal intensive care unit (NICU) is the only referral center for neonates with surgical NEC in Mississippi. From a detailed review of the electronic medical records, we identified 36 patients with surgical NEC (Bell stage III) who underwent management of NEC in the period between January 2013 and June 2018. We excluded infants with incomplete clinical data or confounding medical conditions such as congenital heart disease, intestinal atresia, or spontaneous intestinal perforation (SIP).
Clinical information
We recorded demographic characteristics, including birth weight, gestational age, sex, race/ethnicity, mode of delivery, and inborn versus out born status. We collected information regarding maternal factors, including pregnancy-induced hypertension, chorioamnionitis, and antenatal steroids.
We collected post-operative information such as post-operative ileus days (being nil per oral after surgery), antibiotic days following NEC, time to full feeds (≥120 mL/kg/day), total parenteral nutrition days, length of stay, and hospital mortality were measured. We defined mortality as death due to any reason before hospital discharge. The infants with surgical NEC received antibiotics for 10-121 and vancomycin, amikacin and +/− metronidazole were commonly antibiotics used following surgical NEC 12 as reported in our previous study. The surgeons and the surgical practices did not change during the study period.
Histopathological evaluation
Hematoxylin & eosin-stained intestinal tissue sections were evaluated by a pathologist for necrosis, inflammation, hemorrhage, and reparative changes. The tissue sections were taken from the representative most affected area of the intestine, which was determined by the trained grossing personnel blinded with the clinical course of the included infants. The severity of NEC was assessed by the depth to which these histopathological changes (coagulative necrosis, inflammation, hemorrhages, and reparative changes), were seen; grade 1 was limited to the mucosa, grade 2 changes extended to the submucosa, grade 3 to the muscularis, and grade 4 change was transmural. The resected intestine was evaluated by the depth up to which these changes were seen. A score of 0 was assigned when the exam appeared normal, 1 for 1%–25% necrosis/ inflammation, 2 when 25%–50% area involved, 3 when 50%–75% area was affected, and 4 when >75% changes were seen. 17 The reparative changes included the presence of neovascularization, fibroblasts, and the epithelial regeneration.
16S sequencing
The 16S rRNA sequencing protocol targets the V3-V4 region, which was amplified by PCR using barcoded Illumina adapter-containing primers 341F (5′- CCTACGGGNGGCWGCAG-3′) and 806R (5′- GACTACHVGGGTATCTAATCC-3′) and sequenced with 2 × 300 bp cartridges. At the start of each forward read, adapters matching the forward primer were identified over a minimum of 10 bases with an allowable error rate of 0.1. Reads beginning with the designated primer had the primer removed, and the remaining sequence was retained. Reads without a primer match were discarded. The same approach was applied to the reverse reads, using the reverse primer.
Adapters were trimmed using Cutadapt, and per-base quality scores were assessed with QIIME2. The reverse reads typically exhibited lower quality compared to the forward reads.
To ensure high-quality reads and successful merging of paired-end sequences, we applied trimming based on both quality scores and the expected amplicon length. Specifically, forward reads were trimmed at base 284 and reverse reads at base 220. These trimming positions were selected after examining per-base quality score distributions, ensuring that the retained sequences maintained high quality while also preserving sufficient overlap for accurate merging.
The V3–V4 region of the 16S rRNA gene targeted by our primers is approximately 464 base pairs in length. Given the read lengths after trimming (284 bp forward + 220 bp reverse = 504 bp total), this configuration provides approximately 40 bp of overlap between the forward and reverse reads. This overlap is essential for DADA2’s merging algorithm to confidently join read pairs and correct sequencing errors.
All trimming, denoising, and merging steps were performed using the DADA2 pipeline, which also includes quality filtering and chimera removal. Demultiplexing was performed prior to this step using the sample-specific index sequences.
Next, we constructed the feature table and feature data. The feature table details the presence and abundance of amplicon sequence variants (ASVs) across samples, indicating the number of occurrences of each ASV per sample. The feature data consists of the sequences that define each ASV. Summaries for both files were generated and explored to assess their content and distribution.
For taxonomy assignment, a Naïve-Bayes classifier was trained using the latest SILVA SSU database (version 138.1) as the reference. The reference sequences were trimmed to include only the V3-V4 region, and Eukaryotes were excluded from the reference dataset. The reference taxa and sequences were then dereplicated before being used to train the classifier. These steps were carried out using the RESCRIPt tool in QIIME2.
Using the trimmed, denoised, and demultiplexed sequences, taxonomy assignment was performed with the Naive Bayes classifier trained earlier. Reads classified as Archaea, along with those that could not be assigned at the phylum level, were discarded. The feature tables were subsequently filtered to reflect these exclusions.
Statistical analysis
Using the filtered ASV sequences, a phylogenetic tree was constructed with QIIME2. The final rooted phylogenetic tree will serve as the foundation for calculating phylogenetically informed diversity metrics. Samples were rarefied to 699,778 reads to retain all samples for further analysis while maximizing sampling depth as much as possible. Alpha- and beta diversity metrics were calculated using QIIME2, from the constructed rooted tree, the filtered feature table, and the specified sampling depth.
For each individual, Shannon’s entropy, Pielou’s evenness, Faith’s Phylogenetic Diversity, and the number of observed features were calculated as alpha diversity metrics. Group-wise comparisons between Non-survivors and Survivors were conducted using the Mann-Whitney U rank sum test.
For beta diversity, a distance metric was used to calculate the distance between all possible sample pairs. Unweighted UniFrac distances, weighted UniFrac distances and Bray-Curtis dissimilarity were employed as the beta diversity metrics. The weighted UniFrac distances between sample pairs were visualized using UMAP. To assess differences in overall microbial community structure, we performed a Permutational Multivariate Analysis of Variance (PERMANOVA). The test was conducted on Bray-Curtis dissimilarity, weighted UniFrac, and unweighted UniFrac distance matrices using the adonis2 function from the vegan R package with 999 permutations.
The bias-corrected differential abundance analysis of the microbiome between survivor and non-survivor groups was performed using ANCOM-BC, integrated within QIIME2. Enrichments with adjusted p-values below 0.01 were considered significant, and the results were collapsed at various taxonomic levels.
Additionally, the abundance levels of different microbiota in each sample were directly compared using the Mann-Whitney U rank sum test, with p-values below 0.05 deemed significant. Direct comparison of abundance levels is particularly useful for clinicians in making therapeutic decisions.
To build upon these univariate comparisons and provide a more robust statistical assessment, we utilized the MaAsLin2 R package. This analysis was performed to determine the specific microbial taxa whose relative abundances were significantly associated with survival status in a multivariable linear model. This approach allows for a more rigorous evaluation by accounting for the underlying data distribution and, if necessary, adjusting for clinical confounders. For this analysis, associations with a False Discovery Rate (FDR) corrected p-value (q-value) below 0.25 were considered statistically significant.
Statistical methods for clinical data
Demographic and clinical factors in infants with surgical NEC were compared in survivors and the non-survivors. Continuous data were summarized as median (1st quartile, 3rd quartile) with Mann-Whitney U or Kruskal-Wallis tests for differences. Categorical data were expressed as counts and percentages, with group differences assessed by Chi-squared or Fisher’s exact tests. All analyses were performed in SPSS and a p-value < 0.05 was considered statistically significant.
Results
Clinical information of the included preterm infants (death vs discharge).
Diversity analysis
Beta diversity analyses revealed that microbial community composition significantly differed by survival status. Unweighted UniFrac distances, which consider the presence or absence of phylogenetically distinct taxa, showed a significant separation between survivor groups (p = 0.001, R2 = 0.195), indicating that approximately 19.5% of the variation in community structure is explained by survival status. Weighted UniFrac distances, which account for both phylogeny and relative abundance, demonstrated an even stronger association (p = 0.001, R2 = 0.593), suggesting substantial divergence in microbial communities based on survival. Similarly, Bray-Curtis dissimilarity revealed a significant difference in community composition between survivors and non-survivors (p = 0.001, R2 = 0.497), with survival status explaining nearly half of the variance. In contrast, no significant differences in microbial composition were observed for hemorrhage group (p = 0.715), inflammation level (p = 0.362), or necrosis score (p = 0.720), indicating these clinical features did not account for substantial variation in community structure.
Alpha diversity analyses revealed significant differences in microbial diversity between survivor groups. Non-survivors exhibited greater microbial diversity based on Shannon’s entropy (p = 0.02) and higher Pielou’s evenness (p = 0.017), suggesting a more even and compositionally complex community. In contrast, survivors had a significantly higher number of observed features (p = 0.004), indicating greater richness. These findings highlight distinct patterns in community complexity and richness associated with survival status.
The clinical groups based on necrosis, inflammation, and hemorrhage did not show distinct communities based on alpha- and beta diversity.
The data are summarized in Figures 1 and 2. Alpha diversity comparison between groups, and the corresponding p-values. Community structure of survivor groups with respect to beta diversity metrics.

Comparison of microbial abundance between survival groups
Phylum level
At phylum level, the non-survivor infants had significant higher proportion of Actinobacteriota and Firmicutes phyla (p < 0.01) than survivors. However, the survivors had significant higher proportion of Proteobacteria (p < 0.001) than the non-survivors (see Figure 3). Phylum level bacteria abundance between two groups with p < 0.05.
Class level
At class level, the non-survivor infants had significantly higher proportion of Actinobacteria (p<0.001) alpha proteobacteria (p < 0.05) and bacilli (p < 0.01) than survivors. However, the survivors had significantly higher proportion of Gammaproteobacteria (p < 0.05) than the non-survivors (see Figure 4). Class level bacteria abundance between two groups with p < 0.05.
Order level
At order level, the non-survivor infants had significant higher proportion of Corynebacteriales (p < 0.001), Propionibacteriales (p < 0.001), Staphylococcales (p < 0.01). However, the survivors had higher proportion of Burkholderiales (p < 0.001), Xanthomonadales (p < 0.001), Pseudomonadales (p < 0.001) and Sphingobacteriales (p < 0.05) than the non-survivors (see Figure 5). Order level bacteria abundance between two groups with p < 0.05.
Family level
At family level, the non-survivor infants had significant higher proportion of Dermacoccaceae (p < 0.05), Enterobacteriaceae (p < 0.01), Xanthobacteraceae (p < 0.05), Staphylococcaceae (p < 0.01), Corynebacteriaceae (p < 0.01), Propionibacteriaceae (p < 0.001)than the survivors (p < 0.05). The survivors had the higher proportion of Comamonadaceae (p < 0.001), Xanthomonadaceae (p < 0.001), Pseudomonadaceae (p < 0.001), and Rhizobiaceae (p < 0.05) than the non-survivors (see Figure 6). Family level bacteria abundance between two groups with p < 0.05.
Genus level
At genus level, the non-survivor infants had higher proportion of Dermacoccus, Escherichia-Shigella, Rhodopseudomonas, Staphylococcus, and Cutibacterium (p < 0.05) than survivors. However, the survivors had higher proportion of Stenotrophomonas, Delftia, and Phyllobacteriumthan than non-survivors (see Figure 7). Genus level bacteria abundance between two groups with p < 0.05.
Multivariable association analysis of microbial taxa
To build upon the initial univariate comparisons from previous section and to more robustly quantify these findings, we next employed a more advanced linear modeling approach. As recommended for its statistical power and flexibility, this method allows for the calculation of precise effect sizes and can account for potential confounding variables. For this, we utilized the MaAsLin2 tool to model the direct linear association between the relative abundance of each microbial taxon and the survival status outcome.
Prior to the final analysis, we explicitly evaluated the potential confounding effect of other available metadata, including necrosis, hemorrhage, and inflammation scores. A preliminary statistical analysis revealed no significant association between any of these variables and the survival status outcome in our cohort (all p > 0.5). Given these results, and to prioritize statistical power and avoid the significant risk of overfitting the model with our sample size, we proceeded with a simple, unadjusted model. This final model therefore assesses the direct association between each taxon and survival status, providing a clear effect size (coefficient) and a False Discovery Rate (FDR) for each potential biomarker. A positive coefficient of (e.g., 6.41) indicates that a higher relative abundance of this microbe is associated with the survivor group. The magnitude of the coefficient represents the estimated difference in mean abundance.
Notably, the taxa identified as significant through this robust linear modeling approach were in strong agreement with the initial findings from the direct comparison tests detailed in previous section. This consistency across different statistical methods strengthens the evidence for these specific microbes as potential biomarkers of survival status. A comprehensive set of plots illustrating all significant associations from the phylum to the genus level is provided in Figures 8–12. Differentially abundant phyla associated with survival status. Boxplots show the relative abundance of microbial phyla found to be significantly associated with survival. Associations were determined using linear models (MaAsLin2), with significance defined as a False Discovery Rate (FDR) < 0.25. The calculated coefficient (effect size) and FDR are displayed for each significant taxon. Differentially abundant classes associated with survival status (MaAsLin2). Differentially abundant orders associated with survival status (MaAsLin2). Differentially abundant genera associated with survival status (MaAsLin2). Differentially abundant families associated with survival status (MaAsLin2).




Discussion
Our data demonstrate that the FFPE intestinal tissue enabled retrospective and spatially relevant microbiota assessment at the disease site which captured the mucosal niche microbes and the tissue microbes in preterm infants with surgical NEC. There may be some loss of mucosal surface microbes due to tissue fixation and so the 16s is more reflective of the tissue compartment. In our NEC cohort, we noticed the bacterial communities belonging to the Proteobacteria, Acinetobacter, Firmicutes groups as reported in published studies. 14 However, any bacterial communities’ differences in the survivors and the non-survivors have not been reported before. Therefore, our novel data highlights the microbial differences in the survivors and the non-survivors.
The non-survivors had complex microbiota and had distinct bacterial communities compared to survivors. The increased diversity in non-survivors may reflect tissue invasion, not beneficial diversity. However, the clinical groups based on intestinal histopathology such as necrosis, inflammation, and hemorrhage did not show distinct communities. The non-survivor infants have a more diverse and even microbiome which may seem counterintuitive. This could reflect a failure of the maternal or infant immune system shaping the microbes surviving in the intestinal tissue.
In our cohort at the phylum level, the non-survivor infants had significant relative abundance of Actinobacteriota and Firmicutes phyla than the survivors. A report by Stewart et al 18 reported the profiling the microbiota at the site of NEC disease using the FFPE intestinal tissue and found the relative abundance of Proteobacteria at the phylum levels like our study. However, they did not report any microbiota changes in survivors and the non-survivors.
A recent metanalysis by Pammi et al 14 of 14 studies reported that fecal microbiome from preterm infants with NEC had increased relative abundances of Proteobacteria and decreased relative abundances of Firmicutes and Bacteroidetes prior to NEC onset. Alpha- or beta diversity indices in preterm infants with NEC were not consistently different from controls, but the authors found differences in taxonomic profiles related to antibiotic exposure, formula feeding, and mode of delivery. Exploring heterogeneity revealed differences in microbial profiles by study and the target region of the 16S rRNA gene (V1–V3 or V3–V5). However, the analysis did not look at the microbiome differences in the survivors and the non-survivors.
Warner et al 19 in a prospective study done in stool samples (166 infants, 3586 stools, 46 cases of necrotizing enterocolitis), reported the increased proportions of Gammaproteobacteria (p = 0·0011) and lower proportions of both Negativicutes (p = 0·0013) and the combined Clostridia-Negativicutes class (p = 0·0051) in infants who went on to develop necrotizing enterocolitis compared with controls. These associations were strongest in both the primary cohort and the overall cohort for infants born at less than 27 weeks’ gestation.
In a prospective study by Yu, He al, reported that NEC patients had increased Proteobacteria and decreased Firmicutes and Bacteroidetes compared to fecal control samples, and the level of butyric acid in the NEC group was lower than the control group. Fecal microbial transplant (FMT) in germ-free mice with samples from NEC patients achieved a higher histological injury score when compared to mice that received FMT with control samples. Alterations in microbiota and butyrate levels were maintained in mice following FMT. The ratio of Treg/CD4+T (Thelper) cells was reduced in both NEC patients and mice modeling NEC following FMT.
Lindberg et al 20 reported low levels of diversity in samples obtained from all preterm infants and antibiotic exposure further decreased diversity among both NEC cases and controls. Fecal microbial composition differed between NEC cases and controls, with a greater abundance of Proteobacteria and bacteria belonging to the class Gammaproteobacteria among NEC infants. Control infants demonstrated a greater abundance of bacteria belonging to the phylum Firmicutes.
A report by Wang et al 21 found that NEC is associated with severe lack of microbiota diversity that may accentuate the impact of single dominant microorganisms favored by empiric and widespread use of antibiotics. In their cohort, fecal samples from 20 preterm infants, 10 with NEC and 10 matched controls (including 4 twin pairs) showed infants with NEC had even less diversity, an increase in abundance of Gammaproteobacteria, a decrease in other bacteria species, and had received a higher mean number of previous days of antibiotics. 21
Yu et al 22 investigated the impact of early preterm infant microbiota on host gut development using a gnotobiotic mouse model. The authors did the histological assessment of intestinal development and the differentiation of four epithelial cell lineages (enterocytes, goblet cells, Paneth cells, and enteroendocrine cells) and tight junction (TJ) formation was examined. Using weight gain as a surrogate marker for health, they found that early microbiota from a preterm infant with normal weight gain (MPI-H) induced increased villus height and crypt depth, increased cell proliferation, increased numbers of goblet cells and Paneth cells, and enhanced TJs compared with the changes induced by early microbiota from a poor weight gain preterm infant (MPI-L). Laser capture microdissection (LCM) plus qRT-PCR further revealed, in MPI-H mice, a higher expression of stem cell marker Lgr5 and Paneth cell markers Lyz1 and Cryptdin5 in crypt populations, along with higher expression of the goblet cell and mature enterocyte marker Muc3 in villus populations. In contrast, MPI-L microbiota failed to induce the changes and presented intestinal characteristics comparable to a germ-free host. The study concluded that microbial communities have differential effects on intestinal development and the future studies to identify pioneer settlers in neonatal microbial communities necessary to induce maturation may provide new insights for preterm infant microbial ecosystem therapeutics.
Our study differs as the microbes are sourced from tissue (not lumen or surface) and thus may reflect the microbes translocating from the lumen into the lamina propria and submucosal tissues and not that in the lumen of the GI. This is an important distinction perhaps from other studies and may point to the barrier function and immune system as key differences in the creation of a different tissue microbiome in survivors and the non-survivors in infants with surgical NEC. In our study, microbiome assessment by 16s is more reflective of the tissue compartment than the mucosal surface due to mucosal loss due to tissue fixation and handling.
Our study is limited by single-center experience and a small sample size which affects the generalizability of our results. Our cohort did not have the stool samples to compare with the intestinal tissue results. Most infants with NEC are exposed to antibiotics before surgery, that could possibly alter the microbiome, and FFPE samples are prone to contamination. We did not take controls from the other neonatal gastrointestinal conditions such as atresia’s and the malrotations as the they have different age of presentation, pathology, and outcomes. However, our data highlights the complex and diverse microbial communities in the survivors and the non-survivors in the preterm infants with surgical NEC. Our data supports the idea that intestinal FFPE blocks can be used for profiling the microbiota at the disease site. In future, larger multicentric prospective studies are needed to assess the impact of complex microbial communities on the clinical outcomes. There is a need to develop the microbiome targeted therapeutic intervention and approaches to have improved outcomes in preterm infants with surgical NEC.
Impact
(1) The non-survivors had complex microbiota and had distinct bacterial communities compared to survivors. The increased diversity in non-survivors may reflect tissue invasion, not beneficial diversity. (2) The FFPE intestinal tissue enables retrospective and spatially relevant microbiota assessment in surgical NEC infants. (3) There is a need to develop the microbiome targeted therapeutic intervention and approaches to have improved outcomes in preterm infants with surgical NEC.
Footnotes
Acknowledgment
The Mississippi Center for Clinical and Translational Research and the Department of Pediatrics at the wake forest school of medicine for supporting the NEC research.
Author contribution
PMG designed the study; PMG, PZ, LL, WZ, KK, TH, DS, NV, GH, DJS, and PP analyzed the collected data. PMG, PZ, and all the authors wrote the article and approved the manuscript.
Funding
The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Parvesh Garg is partially supported by the NIGMS of the NIH under Award Number U54GM115428. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.
Declaration of conflicting interests
The authors declared the following potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Data Availability Statement
All data generated and analyzed during this study are included in this article and its supplementary information files.
