Distinct mucosal microbial communities in infants with surgical necrotizing enterocolitis correlate with age and antibiotic exposure
Distinct mucosal microbial communities in infants with surgical necrotizing enterocolitis correlate with age and antibiotic exposure
Joann Romano-Keeler 0 4 5 6
Meghan H. Shilts 2 4 5 6
Andrey Tovchigrechko 2 4 5 6
Chunlin Wang 4 5 6
Robert M. Brucker 1 4 5 6
Daniel J. Moore 0 3 4 5 6
Christopher Fonnesbeck 4 5 6
Shufang Meng 0 4 5 6
Hernan Correa 3 4 5 6
Harold N. Lovvorn 4 5 6
Yi-Wei Tang 4 5 6
Lora Hooper 4 5 6
Seth R. Bordenstein 1 3 4 5 6
Suman R. DasID 2 4 5 6
Jo? rn-Hendrik WeitkampID 0 4 5 6
0 Department of Pediatrics, Vanderbilt University , Nashville, Tennessee , United States of America
1 Department of Biological Sciences, Vanderbilt University , Nashville, Tennessee , United States of America
2 Department of Medicine, Vanderbilt University Medical Center, Nashville, Tennessee, United States of America , 3 Research Bioinformatics, Medimmune, Gaithersburg, Maryland , Tennessee, United States of America, 4 Genome Technology Center, Stanford University , Palo Alto, California , United States of America
3 Department of Pathology, Microbiology & Immunology, Vanderbilt University , Nashville , Tennessee, United States of America, 7 Vanderbilt Institute for Infection, Immunology and Inflammation, Vanderbilt University Medical University , Nashville , Tennessee, United States of America, 8 Department of Biostatistics, Vanderbilt University , Nashville , Tennessee, United States of America, 9 Department of Pediatric Surgery, Vanderbilt University , Nashville , Tennessee, United States of America, 10 Department of Laboratory Medicine, Memorial Sloan Kettering Cancer Center , New York , New York, United States of America, 11 Department of Immunology, The University of Texas Southwestern Medical Center , Dallas, Texas , United States of America
4 Editor: Efrem Lim, Arizona State University , UNITED STATES
5 Current address: The Rowland Institute at Harvard, Harvard University , Cambridge, Massachusetts , United States of America
6 American Academy of Pediatrics Marshall Klaus Perinatal Research Award (to J.R.K.) and the Eunice Kennedy Shriver National Institute Of Child Health & Human Development (NICHD) [T32HD068256 (to J.R.K.) , K08HD061607 (to J.H
Necrotizing enterocolitis (NEC) is the most common surgical emergency in preterm infants,
and pathogenesis associates with changes in the fecal microbiome. As fecal samples
incompletely represent microbial communities in intestinal mucosa, we sought to determine
the NEC tissue-specific microbiome and assess its contribution to pathogenesis.
We amplified and sequenced the V1-V3 hypervariable region of the bacterial 16S rRNA
gene extracted from intestinal tissue and corresponding fecal samples from 12 surgical
patients with NEC and 14 surgical patients without NEC. Low quality and non-bacterial
sequences were removed, and taxonomic assignment was made with the Ribosomal
Database Project. Operational taxonomic units were clustered at 97%. We tested for differences
between NEC and non-NEC samples in microbiome alpha- and beta-diversity and
differential abundance of specific taxa between NEC and non-NEC samples. Additional analyses
were performed to assess the contribution of other demographic and environmental
confounding factors on the infant tissue and fecal microbiome.
W.)]; the National Institute of Diabetes and
Digestive and Kidney Diseases (NIDDK)
[K08DK090146 (to D.J.M.)], the National Institute
of Health (NIH) Division of Loan Repayment
National Institute of Diabetes and Digestive and
Kidney Diseases Award (to J.R.K.), the National
Science Foundation [DEB-1046149] and The
Vanderbilt Microbiome Initiative (to S.R.B.), and
U19AI095227 and P30 AI110527, the National
Institute of Allergy and Infectious Diseases (NIAID)
(to S.R.D.). The content is solely the responsibility
of the authors and does not necessarily represent
the official views of the NICHD, NIDDK, NIAID or
the National Institutes of Health (NIH). This project
was also funded by the Vanderbilt Digestive
Disease Research Center [P30DK058404],
Vanderbilt Diabetes Center [P30DK20593], and the
Vanderbilt CTSA Grant UL1 RR024975-01 from the
National Center for Research Resources (NCRR/
NIH). The funders had no role in study design, data
collection and analysis, decision to publish, or
preparation of the manuscript.
Competing interests: The authors have declared
that no competing interests exist.
The fecal and tissue microbial communities were different. NEC was associated with a
distinct microbiome, which was characterized by low diversity, higher abundances of
Staphylococcus and Clostridium_sensu_stricto, and lower abundances of Actinomyces and
Corynebacterium. Infant age and vancomycin exposure correlated with shifts in the tissue
The observed low diversity in NEC tissues suggests that NEC is associated with a bacterial
bloom and a distinct mucosal bacterial community. The exact bacterial species that
constitute the bloom varied by infant and were strongly influenced by age and exposure to
Necrotizing enterocolitis (NEC) is a common and frequently fatal intestinal complication in
premature infants [
]. Experiments in germ-free animals and toll-like receptor targeted
knock out mice strongly suggest a bacterial antigen is critical for the initiation of intestinal
inflammation and NEC development [
]. Bacterial DNA is present in larger quantities in
acute human NEC specimens compared to samples collected after NEC has clinically resolved
. A number of different gram-positive and gram-negative bacteria as well as viruses have
been associated with NEC [
]. Indeed, microbial community studies using 16S rRNA gene
sequencing of the fecal microbiome demonstrate a reduction in microbial community
diversity with a shift towards potentially pathogenic subgroups [
We previously detected significant differences in the microbiome between surgical tissue
and parallel collected fecal samples in preterm infants without NEC [
]. We hypothesized the
existence of a specific microbial profile at the site of injury in the small intestinal mucosa of
premature infants with NEC that has not been previously recognized in fecal microbiome
studies. Hence, we sought to interrogate differences in the tissue-level and fecal microbiomes
in infants with and without NEC to determine bacterial communities at the site of injury and
their representation in feces. As intestinal tissue cannot ethically be collected from healthy
infants, we included infants with intestinal diseases other than NEC in this study for
comparison. We detected a statistically significant increase in the abundance of Staphylococcus and
Clostridium_sensu_stricto in NEC compared to non-NEC tissue samples when controlling for
age and history of antibiotic exposure.
Materials and methods
This study was approved by the Vanderbilt University Institutional Review Board (protocol
number 090161). All infants hospitalized at the Monroe Carell Jr. Children?s Hospital at
Vanderbilt were eligible for the study if they underwent intestinal resection at <180 days of age.
We obtained written informed consent from parents, the next of kin, caretakers, or guardians
on behalf of the minors/children enrolled in the study to permit collection of tissue and
metadata from the medical records including gestational age, birth weight, race, sex, mode of
2 / 18
delivery, maternal or fetal indications for delivery, antibiotic exposure, enteral feeding
regimens, diagnoses and type of surgical resection.
Tissue collected at the time of surgery was gently rinsed with sterile saline solution, and
immediately cryopreserved in sterile containers [
]. Fecal material was collected by either taking
the patient?s first post-operative stool or by scraping surgical tissue; samples were immediately
cryopreserved (Table 1). The clinical and intraoperative diagnosis of NEC was confirmed by a
pediatric pathologist after histologic examination of the resected specimen and by review of
the operative and surgical pathology reports.
DNA extraction and amplification of 16S rRNA gene
We extracted DNA from fresh NEC and non-NEC surgical tissue and corresponding fecal
samples as previously described [
]. Briefly, we extracted DNA from 15?25 mg of intestinal
tissue and 180?200 mg of feces and amplified the V1-V3 hypervariable region of bacterial 16S
rRNA with previously validated primers: 5F (5?-TGGAGAGTTTGATCCTGGCTCAG-3?)
and 532R (5?-TACCGCGGCTGCTGGCAC-3?) [
]. PCR was conducted as described
] and barcoded amplicons were gel purified (Qiagen), quantified, and pooled prior to
sequencing on a 454 FLX Titanium sequencer. Sequencing negative controls?template-free
sterile water, processed with the same DNA extraction and PCR amplification kits as the real
samples?were sequenced on the same run [
Pyrosequencing and data analysis
Sequences generated from the pyrosequencing of barcoded 16S rRNA gene PCR amplicons
were analysed using mothur (http://mothur.org) [
] by following the 454 SOP as of 13 March
2017. Sequences were aligned to the SILVA database release 123 [
] and taxonomically
classified with the Ribosomal Database Project (RDP) classifier 11 [
]. Chimeric sequences as
detected by UCHIME were removed [
]. OTUs were clustered at 97% similarity. Prior to
statistical analysis, samples with <400 reads were discarded (N = 3).
Phylogenetic Investigation of Communities by Reconstruction of Unobserved States
(PICRUSt) was used to predict metagenomic and functional composition of the samples from
16S rRNA sequences [
]. Prior to PICRUSt analyses, closed reference OTUs were picked
against the GreenGenes database 13_5 [
] using uclust [
] in QIIME 1.9.1 [
assignments were made using the RDP Classifier 2.2 . Functions of genes were assigned
using the KEGG Orthology database [
Statistical analysis was performed in R using MGSAT (https://bitbucket.org/andreyto/
mgsat), which wraps a number of R packages, including vegan [
] to perform alpha- and
beta- diversity analyses and DESeq2 [
], GeneSelector [
], and stabsel [
] for testing
taxonomic associations with metadata. When testing taxonomic associations with metadata, we
report the q-values computed with the Benjamini & Hochberg false discovery rate method to
adjust for multiple comparisons [
For diversity and richness estimates, full count matrices as produced by the mothur
annotation were used [
]. To compare microbial alpha diversity estimates between groupings, we
estimated Shannon-Wiener (H?) and Simpson?s diversity indices; to compare microbial
richness estimates, we estimated observed OTUs and calculated S. chao1 estimates. Counts were
rarefied to the lowest library size of all the samples (number of reads per sample = 445), and
then abundance-based and incidence-based alpha diversity indices and richness estimates
were computed. This was repeated multiple times (n = 400), and the results were averaged.
3 / 18
A, G, V, MZ
V, G, P
V, G, M
A, G, V, C,
A, G, V, MZ
A, G, V
A, G, MZ
A, G, V, P
A, G, MZ
A, G, MZ
A, G, CL, V
A, G, CL, V
A, G, C, CP,
CT, M, V
A, G, C, MZ,
A, G, MZ, P,
A, G, M, V,
0, less than 24 hours
? A, ampicillin; G, gentamicin; V, vancomycin; C, cefotaxime; CL, clindamycin; CP, cefepime; CT, ceftriaxone; M, meropenem; P piperacillin-tazobactam; MZ,
metronidazole; T, tobramycin
? Feces adherent to collected mucosa; all other fecal samples collected at patient?s first post-operative stool
NEC, necrotizing enterocolitis; NPO, nil per mouth; EBM, expressed breast milk; DBM, donor breast milk
Incidence-based estimates were computed on pools of observations split by the relevant
metadata attribute, and in each repetition, observations were also stratified to balance the number
of observations at each level of the metadata attribute. Inverted Simpson and Shannon
diversity indices were converted into corresponding Hill numbers [
]. Linear models were fit to
test for associations between abundance-based richness and diversity estimates and metadata
We applied the PermANOVA (permutation-based analysis of variance) [
] test of
statistical significance (as implemented in the Adonis function of the R vegan package) [
] on the
association between the abundance profile dissimilarities and the metadata variables. We used
4 / 18
the Bray-Curtis dissimilarity index [
] and 4000 permutations. The counts were normalized
to simple proportions within each observation.
When differential abundance analysis was performed, in order to remove the likely
noninformative features and to reduce the associated penalty from the multiple testing correction
applied after univariate tests, we used unbiased metadata-independent filtering at each
taxonomy level by eliminating all taxa that were detected with a mean proportional abundance of
less than 0.0005. The absolute counts from the removed features were aggregated into a
category ?other,? which was taken into an account when computing simple proportions during
data normalization, but were otherwise discarded. When testing taxonomic associations with
metadata, for each feature, we also obtained, from the same test done on the full dataset, the
pvalue computed using the test implementation from R exactRankTests package [
qvalue computed with the Benjamini & Hochberg false discovery rate method in the package
function p.adjust [
], and several types of the effect size such as common language effect size
and rank biserial correlation [
]. To evaluate the influence of confounders, models were built
in DESeq2 with pre-selected covariates added in.
Stabsel is a stability selection approach implemented in the R package stabs [
feature selection method implements a stability selection procedure described in [
] with the
improved error bounds described in [
]. Elastic net (from R package glmnet ) was used as
the base feature selection method that was wrapped by the stability protocol. For groupings
with two factor levels, a binomial family model was built with the grouping as a response and
the matrix of the abundance values as predictors. The mixing parameter ? of glmnet was
selected based on a 15-fold cross-validation minimizing deviance on the full dataset. The
predictors were first normalized to simple proportions within each multivariate observation,
transformed with the inverse hyperbolic sign log?x ? ?x2 ? 1??, and then standardized to
zero means and unit variances. With its multivariate base feature selection method, this
protocol can potentially detect those correlated groups of biologically relevant features that will be
missed by the univariate methods. The ranking of taxa and their probability of being selected
into the model were reported, as well as the probability cutoff corresponding to the per-family
error rate (PFER) that is controlled by this method. Our PFER cutoff was set to 0.05, and the
target number of features selected by the base classifier was set to ?0:8 p? where p is the
total number of features (39). In our experience with omics datasets, the PFER control in this
method is fairly conservative, and we typically look at the ranking of features as opposed to
only concentrating on features that pass the PFER cutoff.
All sequences reported in this paper have been deposited into the NCBI sequence short read
archive (accession no. SRR7993700-SRR7993745).
Demographic and antimicrobial exposure characteristics were similar
between NEC and non-NEC infants
We collected and analyzed fresh surgical tissue and corresponding fecal samples from 10
patients with NEC and 14 patients without NEC; in total, 44 samples were analyzed (fecal
N = 21; tissue N = 23) (Table 1). Surgical samples included patients with spontaneous
intestinal perforations, ileal and jejunal atresias, midgut volvulus, and mesenteric ischemic bowel
injuries. Mean gestational age, birth weight and postnatal age were 29 weeks (range 25?33
weeks), 1,274 grams (range 440?2,101 grams), and 17 days (range 5?46 days) for NEC infants
5 / 18
and 30 weeks (range 24?39 weeks), 1,662 grams (range 650?3,454 grams), and 31 days (0?132
days) for non-NEC patients, respectively (all t-tests p>0.05). Female infants represented 60%
and 50% of the study population in the NEC and non-NEC groups, respectively. Except for
two colon samples among the non-NEC group and two colon samples within the NEC group,
all analyzed tissues were from the ileum or jejunum. For the non-NEC group, one fecal sample
(C5) was adherent to the mucosa when collected, for the NEC group there was one (N27). All
but one infant from the non-NEC group had perinatal antibiotic exposure. Mean number of
antibiotic exposure days prior to surgery were less in the NEC group (5 days, range 1?22)
compared to the non-NEC group (17 days, range 0?131) but means were not statistically different
(t-test with Welch?s correction, p = 0.180). Both the NEC and non-NEC groups contained
infants receiving breast milk, infant formula, or no enteral nutrition prior to sample collection.
Of non-NEC infants, 36% were delivered via C-section compared to 50% of infants in the
Microbial diversity was reduced in NEC samples compared to non-NEC
After quality filtering and removal of chimeras and non-bacterial sequences, barcoded 16S
rRNA amplicons generated a total of 59,778 sequences for fecal and 72,791 sequences for tissue
samples. The mean (range) number of fecal sample sequences was 4,719 (697?15,319) for NEC
and 2,510 (589?6,530) for non-NEC subjects, and the mean (range) number of tissue sample
sequences was 2,322 (445?6,066) sequences for NEC and 2,799 (634?7,906) for non-NEC
Prior to estimating microbial alpha-diversity or richness, samples were rarefied to the
lowest library size of all the samples (445 reads per sample). When testing for a non-zero
coefficient of a normal linear model that used NEC/non-NEC group membership as predictor of
richness, microbial richness and diversity were lower in NEC samples compared to non-NEC
samples (Fig 1A). In tissue samples, when comparing microbial richness or diversity in tissue
from NEC and non-NEC subjects, p-values for all tested richness estimates were <0.05 and
there was a trend towards lower alpha diversity estimates in tissue from infants with NEC
compared to those without NEC (p-values: N1 = 0.081, N2 = 0.168) (Fig 1A). Both microbial
richness (observed OTU counts (S.obs) p-value = 0.046, S.Chao1 p-value = 0.065) and alpha
diversity (N1 p-value = 0.075, N2 p-value = 0.078) were at or near significantly lower in stool
from infants with NEC compared to those without (Fig 1A).
NEC and non-NEC samples exhibited distinct microbial profiles
Prior to estimating beta diversity, samples were rarefied to the lowest library size (445 reads/
sample). Principal Component Analysis (PCoA) using pairwise Bray-Curtis dissimilarities
demonstrates distinct microbial genus composition of tissue samples from NEC versus
nonNEC patients (Adonis test p-value = 0.0003) (Fig 1B). The microbial communities isolated
from NEC and non-NEC fecal samples were also significantly dissimilar (Fig 2, Bray-Curtis
dissimilarities calculated at the genus level, Adonis test p-value = 0.003). In contrast to the
more uniform pattern in non-NEC tissues, microbial composition in NEC tissue clustered in
separate coordinates indicating discrete colonization types.
Specific taxa associated with the differential microbial profiles of NEC and non-NEC
samples. Fig 3 shows a heatmap of the top 30 most abundant genera found in tissue samples across
the bottom, with sample clustering on the left and each individual sample marked on the right
with both infant age in days at time of collection and whether the sample was from an infant
with or without NEC. NEC and non-NEC samples generally formed two distinct clusters.
6 / 18
Fig 1. A) Boxplots of tissue microbial diversity and richness in infants with and without necrotizing enterocolitis (NEC) at the operational taxonomic unit (OTU) level
for all samples, stool alone, and tissue alone. After rarefaction to the lowest library size of all the samples (445 reads per sample), ? diversity and richness estimates were
calculated per each sample. This process was repeated 400 times and results were averaged. The Shannon and inverse Simpson indices were calculated to estimate
abundance-based OTU diversity, while the Chao1 estimator and observed taxa counts were calculated to estimate abundance-based OTU richness. Displayed p-values
were obtained after testing for a non-zero coefficient of a normal linear model that used NEC/non-NEC group membership as predictor of richness or diversity. All
tested richness and diversity indices for both tissue and stool samples were at or near significantly lower in NEC compared to non-NEC samples. B) Principal
coordinates analysis (PCoA) plot of tissue samples, labelled by NEC status. Bray-Curtis dissimilarities between samples were calculated at the genus level after
normalizing read counts to simple proportions and after rarefaction to the lowest library size (445 reads per sample). The centroids between the NEC and non-NEC
samples were significantly dissimilar (Adonis PerMANOVA p-value = 0.0002).
Bacterial genus level assignments for tissue and fecal samples comparing NEC with non-NEC
patients are depicted in Fig 4; NEC tissue samples were more likely to be dominated by a single
genus (Fig 1A), including Staphylococcus, Clostridium, Escherichia, or Bacteroides than
nonNEC samples. Stool and tissue communities were significantly dissimilar (Bray-Curtis
dissimilarities Adonis test p-value = 0.0005), with tissue and stool communities from the same infant
sharing little overlap (Fig 4).
When tissue samples were analyzed alone, 15 taxa at the genus level had differential
abundances in NEC compared to non-NEC samples with DESeq2 test q-values < 0.1 (Fig 5).
Staphylococcus was ranked first in the DESeq2 model as most significantly different between NEC
and non-NEC samples. Clostridium_sensu_stricto was near significantly more abundant in
NEC tissue compared to non-NEC tissue (Table 2). Both groups were the only two genera
identified as being significantly or near significantly more abundant in NEC compared to
non-NEC samples. Clostridium_sensu_stricto was significantly more abundant in NEC than
non-NEC tissues when the GeneSelector test using Wilcoxon test rankings was applied
(qvalue = 0.021). Clostridium_sensu_stricto abundance being higher in NEC infants appears to
be due mostly to a single infant, N27, who had nearly 100% Clostridium_sensu_stricto
abundance; the relative abundance of this genus was low for the remainder of samples (Fig 4).
7 / 18
Fig 2. Principal coordinates analysis (PCoA) plots of stool samples, labelled by necrotizing enterocolitis (NEC) status. Bray-Curtis
dissimilarities between samples were calculated at the genus level after normalizing read counts to simple proportions. The centroids
between the NEC and non-NEC samples were significantly dissimilar (Adonis PerMANOVA p-value = 0.003).
When fecal samples were analyzed alone, while Clostridium_sensu_stricto did not differ in
abundance between NEC and non-NEC samples, Staphylococcus as a genus was more
abundant during NEC (Table 2), consistent with findings from a recent study describing fecal
microbiome samples from NEC patients . A single Staphylococcus OTU, identified as
OTU0004, was dominated in NEC fecal samples (Table 2) compared to non-NEC samples.
This same Staphylococcus OTU0004 was also found to be significantly more abundant in tissue
samples in infants with NEC (Table 2) compared to those without NEC. Due to the limited
read lengths obtained, this OTU could not confidently be classified below the genus level.
NEC-associated changes in the microbiome were stronger than the
influence of other measured potential confounders
Although infants in the NEC and non-NEC groups were similar demographically and had
similar environmental exposures in aggregate (all p-values >0.05), we conducted additional
analyses to assess the effect of potential gut microbiome confounders.
Mode of delivery did not have a significant correlation with infant microbiome richness,
alpha or beta diversity, or abundance of specific taxa in either stool or tissue samples. Infant
sex did not have a significant correlation with infant microbiome richness or alpha or beta
diversity in either tissue or fecal samples; however, tissue from males had lower abundance of
8 / 18
Fig 3. Heatmap of microbial abundance profiles of infant gut tissue at the genus level. The 30 most abundant genera are shown. Infant age, in days,
and necrotizing enterocolitis (NEC) status are labelled for each sample. Clustering due to NEC status can be observed.
Staphylococcus (Table 2). Age of the infant at time of sampling was found to correlate with
trends in the gut microbiome: overall, the microbial communities significantly differed
between age groups (pairwise Bray-Curtis dissimilarities were calculated between each sample
and infants were quartered into age groups of as even size as possible and the PermANOVA
Adonis test was performed on these groupings; p-value = 0.011). We observed an association
between microbial richness, diversity and infant age and this was strongly correlated with
NEC status (Fig 2).
Prior to sampling, all but one infant had been exposed to antibiotics (Table 1). Of infants
who had received antibiotics, all had received at least two different antibiotics and at least one
broad-spectrum antibiotic. Half of all infants in this study (12/24) were treated with
vancomycin. Staphylococcus abundance was higher in tissue taken from infants who had received
vancomycin (Table 2). In contrast, Staphylococcus was not significantly more abundant in stool
taken from infants with vancomycin exposure. Vancomycin exposure was not significantly
associated with differential abundance of any other taxa or with alpha diversity or richness in
either fecal or tissue samples.
To further assess the influence of confounding variables on the effect of NEC on the
microbiome, we built models in DESeq2 to explicitly account for a priori selected covariates that
may affect the gut microbiome (delivery mode, infant sex, infant age, diet, tissue type, exposure
to vancomycin). Regardless of covariates added, the DESeq2 calculated log2 fold change in
Staphylococcus abundance between infants with and without NEC directionality did not
change (i.e., Staphylococcus abundance was always higher in NEC infants). Infant age and
9 / 18
Fig 4. Stacked bar graph showing genus level taxonomic composition for each individual sample, expressed as a proportion of reads. Infant
sample ID is on the x-axis. The top 15 genera with the highest average relative abundance are shown. Samples are stratified by necrotizing
enterocolitis (NEC) status, and whether sample was tissue or stool. Tissue and stool samples from the same infant had dissimilar microbial
exposure to vancomycin had the strongest effect on the association between NEC on
Staphylococcus abundance: young infants with NEC who had been exposed to vancomycin generally
had high Staphylococcus tissue abundance (Fig 6A).
For ethical reasons, intestinal tissue cannot be collected from healthy infants; therefore, all
the non-NEC infants in this study were in the hospital for ailments unrelated to NEC. There
were two disparate groups of non-NEC infants recruited for this study: those who were very
young (<1 week old) and those who were older (40+ days). To assess the impact of these two
10 / 18
Fig 5. Comparison of the abundance of tissue bacterial genera between infants with and without necrotizing enterocolitis (NEC). Only bacterial
genera that were significantly different between groups after adjusting for multiple comparisons using the DESeq2 package (see text for details) are
indicated by an asterisk. Log2 fold change and log2 fold change standard error of tissue bacterial genera according to NEC status as calculated with the
DESeq2 analysis. A log2 fold change of >0 (pink bars) indicates that abundance was detected to be higher during NEC, while a log2 fold change <0
(blue bars) indicates that abundance was detected to be higher in infants without NEC.
disparate non-NEC groups, we repeated our main analyses, after separating out the samples
into three groups: 1) NEC 2) ?young? non-NEC and 3) ?old? non-NEC. A PCoA plot
constructed using pairwise Bray-Curtis dissimilarities at the genus level revealed that the two
nonNEC groups clustered together, separately from samples from infants with NEC (S1 Fig).
Additionally, microbial richness was lowest in infants who had NEC and was similar in the
two non-NEC groups, regardless of infant age (S2 Fig). Staphylococcus abundance was highest
in the NEC group, followed by the youngest non-NEC group, and lowest in the older
nonNEC group suggesting intestinal Staphylococcus colonization at early age.
PICRUSt analyses reveals different functional profiles in NEC compared to
Overall, the predicted functional profile in NEC and non-NEC samples were generally distinct
in stool (Fig 7A) and tissue samples (Fig 7B). Processes and pathways related to signatures of
infectious diseases, i.e., bacterial toxins (base mean = 661.6, log2 fold change = 1.049,
qvalue = 2.045e-05) and Staphylococcus aureus infection (base mean = 831.7, log2 fold
change = 0.724, q-value = 1.209e-02) were enriched in NEC tissue samples compared to
11 / 18
Base means are calculated in the DESeq2 package for each taxon after normalizing read counts for each sample to account for differences in sequencing depth.
? Log2 fold changes are calculated by the DESeq2 package and indicate the magnitude of the difference in abundance for each comparison. For categorical tests, positive
values indicate that DESeq2 estimated the taxon was more abundant in the first tested group while negative values indicate that DESeq2 estimated the taxon was more
abundant in the second group. When age in days was used as the group to test, a positive value indicates that DESeq2 found that taxon increased in abundance with each
? Reported q-values are the result of a Wald test with the Benjamini and Hochberg correction for multiple comparisons.
NEC, necrotizing enterocolit
Only a few studies have interrogated the tissue-level intestinal microbiome in NEC, despite the
relative proximal location of intestinal injury and previous reports on the existence of a
sitespecific intestinal microbiome [
]. Here, we report a tissue-specific overrepresentation
of Firmicutes, specifically Staphylococcus sp. and Clostridium sp. in NEC. We are aware of only
two other reports on the NEC tissue-level microbiome in humans: a study from Denmark
performed a retrospective analysis of formalin-fixed and paraffin-embedded tissue specimens
using fluorescent in situ hybridization with bacterial rRNA-targeting oligonucleotide probes
. They detected Proteobacteria (49.0%), Firmicutes (30.4%), Actinobacteria (17.1%) and
Bacteroidetes (3.6%) in tissue samples. More recently Brower-Sinning et al. applied 16S rRNA
technology to compare the microbiome of 16 cryopreserved NEC samples and 10 controls
. Except for a higher bacterial load in NEC tissues, no statistically significant distinction
was found between the composition of NEC and non-NEC microbial communities. The
different results in our study may be explained by the fact that in the work by Brower-Sinning et al.
all but one control patient were former NEC patients. In contrast, we included samples from
infants with no history of NEC.
While we observed that the infant gut microbiome was significantly dissimilar in infants
with NEC compared to those without NEC, we conducted a number of analyses to test the
influence of potential confounders. After adding multiple covariates to the model in DESeq2
suggests that a combination of variables is likely to influence the infant tissue microbiome, for
example age, vancomycin exposure, and NEC were found to correlate with Staphylococcus
abundance. We observed that very young infants with NEC who had been exposed to
vancomycin were most likely to have high Staphylococcus abundance in their gut tissue. We do not
know how to explain this unexpected finding except by the fact that vancomycin does not
penetrate tissue very well. Delivery by C-section has been associated with colonization of the
neonate with Staphylococcus . Therefore, we were surprised by our finding that mode of
12 / 18
Fig 6. Age of the infant in days at time of sample collection is plotted on the x-axis and Staphylococcus relative abundance is plotted
on the y-axis. Tissue Staphylococcus abundance is highest in infants with early necrotizing enterocolitis (NEC) resection, who received
delivery did not correlate with specific taxa in our dataset. However, three out of four samples
with high abundance of Staphylococcus were from C-section-delivered infants indicating that
our sample size may have been insufficient to detect a statistical significance.
One unique aspect of our study is the direct comparison between tissue and fecal
samples. This allows for an additional level of quality control as each patient is his/her own
control and results between fecal and tissue samples were distinct in both non-NEC and
NEC patients. Consistent with previous studies in preterm infants [
], we confirmed
the dominant phyla as Proteobacteria and Firmicutes, with a smaller contribution (<20%)
from Bacteroidetes and Actinobacteria. Several fecal microbiome studies reported a bloom of
13 / 18
Fig 7. Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) was used to predict metagenomic and functional
composition of the samples from 716S rRNA sequences. Heatmaps of normalized counts of microbial function pathways detected with the PICRUSt pipeline for
infant A) stool and B) tissue are shown. Infant age, in days, and necrotizing enterocolitis (NEC) status are labelled for each sample. Both sample types display
clustering associated with NEC status.
?-Proteobacteria with a concomitant decrease in Firmicutes in NEC patients [
]. This shift
in microbial communities in NEC patients appears to start 1?2 weeks prior to diagnosis and
has been associated with metabolic changes . While our data do not replicate this shift in
Proteobacteria in fecal samples, possibly as we measured the gut microbiome during rather
than prior to NEC diagnosis, we confirmed the previously reported reduced microbial
diversity and loss of Actinobacteria in NEC patients, especially patients with severe (surgical)
Given numerous previous reports on the dominance of Proteobacteria in NEC [
were surprised to find the high prevalence of Firmicutes and specifically Staphylococcus sp. in
NEC tissue. However, different forms of dysbiosis have been reported in NEC [
including recently an association between Clostridium and Staphylococcus with NEC in European
preterm infants . Importantly, NEC dysbiosis with Firmicutes including Staphylococcus
has been associated with earlier disease and higher mortality . Our study included only
infants with surgical NEC, the group of patients with highest mortality . When comparing
NEC patients with heavy versus light Staphylococcus abundance, NEC patients with high
abundance required surgical resection significantly earlier. Staphylococcus is the major colonizing
organism of the infant gut shortly after birth [
]. In preterm neonates, culture-based
studies detected Staphylococcus in 50% of meconium and 100% of fecal samples from the first
week post-partum . Staphylococcus sp. are frequently cultured from meconium and have
been associated with increased risk for NEC [
Our study has limitations. While we collected tissue and fecal samples prospectively,
technical and ethical limitations do not allow for tissue sampling prior to surgical resection.
Therefore, we cannot perform time series experiments to evaluate the dynamic microbiome changes
in NEC tissue. Similarly, since it is currently not possible to sample intestinal tissue from
normal infants, we lack a healthy control cohort in which to characterize the standard infant tissue
microbiome. While we attempted to match for important variables such as mode of delivery,
antibiotic exposure and type of feeding, given the nature of this human study that explores
both tissue and stool of a surgical emergency in a very vulnerable population, we were not able
to control for all possible microbiome confounders. In addition, the lack of shotgun
14 / 18
metagenomic sequencing prohibits further classification of the bacteria, especially those of
important genera e.g., Staphylococcus and Clostridium identified in this study. However, based
on the recent findings by Roze? et al, we speculate that the majority of Staphylococcus and
Clostridium species would be S. aureus and C. neonatale . Future studies implementing whole
genome sequencing will be necessary to address strain identification and implications for
derangements in metabolic function associated with the distinct microbial community
structure we detected. An additional future aim could be to measure Staphylococcus-specific
endotoxin production in stool samples, especially as our PICRUSt data suggests there was an
increase in bacterial toxin pathways in NEC compared to non-NEC tissue samples.
To the best of our knowledge, we define here for the first time corresponding fecal and
tissuelevel microbial communities comparing NEC patients with patients without a history of NEC
and confirm age and antimicrobial exposure as defining factors.
S1 Fig. Principal coordinates analysis (PCoA) plots of tissue samples, labelled by
necrotizing enterocolitis (NEC) status and whether the infant was from the young or old non-NEC
group. Bray-Curtis dissimilarities between samples were calculated at the OTU level after
normalizing read counts to simple proportions. NEC and non-NEC samples are observed to
cluster separately, while both the young and old non-NEC samples clustered together.
S2 Fig. Microbial richness was estimated using two indices, the Chao estimator (S.chao1)
and estimated number of operational taxonomic units (OTUs) (S. obs) and each of these
indices was plotted as a function of infant age in days. Two age disparate groups of
non-necrotizing enterocolitis (NEC) infants were included in the analysis; however, from this figure it
can be observed that NEC/non-NEC status had a much stronger effect on microbial richness
than infant age.
Conceptualization: Joann Romano-Keeler, Yi-Wei Tang, Jo?rn-Hendrik Weitkamp.
Data curation: Meghan H. Shilts, Andrey Tovchigrechko, Chunlin Wang, Suman R. Das,
Formal analysis: Joann Romano-Keeler, Meghan H. Shilts, Chunlin Wang, Daniel J. Moore,
Christopher Fonnesbeck, Shufang Meng, Suman R. Das.
Funding acquisition: Joann Romano-Keeler, Jo?rn-Hendrik Weitkamp.
Investigation: Joann Romano-Keeler, Jo?rn-Hendrik Weitkamp.
Methodology: Joann Romano-Keeler, Andrey Tovchigrechko, Chunlin Wang, Robert M.
Brucker, Daniel J. Moore, Christopher Fonnesbeck, Hernan Correa, Yi-Wei Tang, Lora
Hooper, Seth R. Bordenstein, Suman R. Das, Jo?rn-Hendrik Weitkamp.
Project administration: Joann Romano-Keeler, Jo?rn-Hendrik Weitkamp.
Resources: Daniel J. Moore, Hernan Correa, Harold N. Lovvorn, III, Yi-Wei Tang, Lora
Hooper, Seth R. Bordenstein, Suman R. Das, Jo?rn-Hendrik Weitkamp.
15 / 18
Supervision: Jo?rn-Hendrik Weitkamp.
Validation: Suman R. Das.
Visualization: Meghan H. Shilts, Christopher Fonnesbeck.
Writing ? original draft: Joann Romano-Keeler, Meghan H. Shilts, Jo?rn-Hendrik Weitkamp.
Writing ? review & editing: Meghan H. Shilts, Chunlin Wang, Robert M. Brucker, Daniel J.
Moore, Christopher Fonnesbeck, Hernan Correa, Harold N. Lovvorn, III, Yi-Wei Tang,
Lora Hooper, Seth R. Bordenstein, Suman R. Das, Jo?rn-Hendrik Weitkamp.
16 / 18
17 / 18
1. Gordon P , Christensen R , Weitkamp JH , Maheshwari A . Mapping the New World of Necrotizing Enterocolitis (NEC): Review and Opinion . EJ Neonatol Res . 2 ( 4 ): 145 - 72 . PMID: 23730536
Weitkamp JH . More than a gut feeling: predicting surgical necrotising enterocolitis . Gut . 2014 ; 63 ( 8 ): 1205 - 6 . https://doi.org/10.1136/gutjnl-2013 -305928 PMID: 24064006
3. Musemeche CA , Kosloske AM , Bartow SA , Umland ET . Comparative effects of ischemia, bacteria, and substrate on the pathogenesis of intestinal necrosis . Journal of pediatric surgery . 1986 ; 21 ( 6 ): 536 - 8 . PMID: 3723307
4. Rozenfeld RA , Liu X , DePlaen I , Hsueh W. Role of gut flora on intestinal group II phospholipase A2 activity and intestinal injury in shock . American journal of physiology Gastrointestinal and liver physiology . 2001 ; 281 ( 4 ): G957 - 63 . https://doi.org/10.1152/ajpgi. 2001 . 281 .4. G957 PMID : 11557516
5. Jilling T , Simon D , Lu J , Meng FJ , Li D , Schy R , et al. The Roles of Bacteria and TLR4 in Rat and Murine Models of Necrotizing Enterocolitis1 . J Immunol . 2006 ; 177 ( 5 ): 3273 - 82 . PMID: 16920968
6. Leaphart CL , Cavallo J , Gribar SC , Cetin S , Li J , Branca MF , et al. A critical role for TLR4 in the pathogenesis of necrotizing enterocolitis by modulating intestinal injury and repair . J Immunol . 2007 ; 179 ( 7 ): 4808 - 20 . PMID: 17878380
7. Bucher BT , McDuffie LA , Shaikh N , Tarr PI , Warner BB , Hamvas A , et al. Bacterial DNA Content in the Intestinal Wall from Infants with Necrotizing Enterocolitis . Journal of pediatric surgery . 2011 ; 46 ( 6 ): 1029 - 33 . https://doi.org/10.1016/j.jpedsurg. 2011 . 03 .026 PMID: 21683193
8. Coggins SA , Wynn JL , Weitkamp JH . Infectious causes of necrotizing enterocolitis . Clinics in perinatology. 2015 ; 42 ( 1 ): 133 - 54 , ix. https://doi.org/10.1016/j.clp. 2014 . 10 .012 PMID: 25678001
9. Claud EC , Keegan KP , Brulc JM , Lu L , Bartels D , Glass E , et al. Bacterial community structure and functional contributions to emergence of health or necrotizing enterocolitis in preterm infants . Microbiome . 2013 ; 1 ( 1 ): 20 . https://doi.org/10.1186/2049-2618-1-20 PMID: 24450928
10. Mai V , Young CM , Ukhanova M , Wang X , Sun Y , Casella G , et al. Fecal microbiota in premature infants prior to necrotizing enterocolitis . PloS one . 2011 ; 6 ( 6 ):e20647. https://doi.org/10.1371/journal.pone. 0020647 PMID: 21674011
11. Sim K , Shaw AG , Randell P , Cox MJ , McClure ZE , Li MS , et al. Dysbiosis anticipating necrotizing enterocolitis in very premature infants. Clinical infectious diseases: an official publication of the Infectious Diseases Society of America . 2015 ; 60 ( 3 ): 389 - 97 .
Wang Y , Hoenig JD , Malin KJ , Qamar S , Petrof EO , Sun J , et al. 16S rRNA gene-based analysis of fecal microbiota from preterm infants with and without necrotizing enterocolitis . ISME J . 2009 ; 3 ( 8 ): 944 - 54 . https://doi.org/10.1038/ismej. 2009 .37 PMID: 19369970
13. Romano-Keeler J , Moore DJ , Wang C , Brucker RM , Fonnesbeck C , Slaughter JC , et al. Early life establishment of site-specific microbial communities in the gut . Gut microbes . 2014 ; 5 ( 2 ): 192 - 201 . https:// doi.org/10.4161/gmic.28442 PMID: 24637795
14. Tang YW , Ellis NM , Hopkins MK , Smith DH , Dodge DE , Persing DH . Comparison of phenotypic and genotypic techniques for identification of unusual aerobic pathogenic gram-negative bacilli . Journal of clinical microbiology . 1998 ; 36 ( 12 ): 3674 - 9 . PMID: 9817894
15. Salter SJ , Cox MJ , Turek EM , Calus ST , Cookson WO , Moffatt MF , et al. Reagent and laboratory contamination can critically impact sequence-based microbiome analyses . BMC biology . 2014 ; 12 : 87 . https://doi.org/10.1186/s12915-014-0087-z PMID: 25387460
16. Schloss PD , Gevers D , Westcott SL . Reducing the effects of PCR amplification and sequencing artifacts on 16S rRNA-based studies . PloS one . 2011 ; 6 ( 12 ):e27310. https://doi.org/10.1371/journal.pone. 0027310 PMID: 22194782
17. Pruesse E , Quast C , Knittel K , Fuchs BM , Ludwig W , Peplies J , et al. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB . Nucleic acids research . 2007 ; 35 ( 21 ): 7188 - 96 . https://doi.org/10.1093/nar/gkm864 PMID: 17947321
18. Cole JR , Wang Q , Cardenas E , Fish J , Chai B , Farris RJ , et al. The Ribosomal Database Project: improved alignments and new tools for rRNA analysis . Nucleic acids research . 2009 ; 37 (Database issue): D141 - 5 . https://doi.org/10.1093/nar/gkn879 PMID: 19004872
19. Edgar RC , Haas BJ , Clemente JC , Quince C , Knight R. UCHIME improves sensitivity and speed of chimera detection . Bioinformatics . 272011. p. 2194 - 200 .
20. Langille MG , Zaneveld J , Caporaso JG , McDonald D , Knights D , Reyes JA , et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences . Nature biotechnology . 2013 ; 31 ( 9 ): 814 - 21 . https://doi.org/10.1038/nbt.2676 PMID: 23975157
21. DeSantis TZ , Hugenholtz P , Larsen N , Rojas M , Brodie EL , Keller K , et al. Greengenes, a chimerachecked 16S rRNA gene database and workbench compatible with ARB . Applied and environmental microbiology . 2006 ; 72 ( 7 ): 5069 - 72 . https://doi.org/10.1128/AEM.03006-05 PMID: 16820507
22. Edgar RC . Search and clustering orders of magnitude faster than BLAST . Bioinformatics . 2010 ; 26 ( 19 ): 2460 - 1 . https://doi.org/10.1093/bioinformatics/btq461 PMID: 20709691
23. Caporaso JG , Kuczynski J , Stombaugh J , Bittinger K , Bushman FD , Costello EK , et al. QIIME allows analysis of high-throughput community sequencing data . Nat Methods. 7. United States2010 . p. 335 - 6 . https://doi.org/10.1038/nmeth.f.303 PMID: 20383131
Wang Q , Garrity GM , Tiedje JM , Cole JR . Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy . Applied and environmental microbiology . 2007 ; 73 ( 16 ): 5261 - 7 . https://doi.org/10.1128/AEM.00062-07 PMID: 17586664
25. Kanehisa M , Institute for Chemical Research KU, Uji, Kyoto 611 -0011, Japan, Sato Y , Healthcare Solutions Department FKSL , Hakata-ku, Fukuoka 812 -0007, Japan, Kawashima M , Healthcare Solutions Department FKSL , Hakata-ku, Fukuoka 812 -0007, Japan , et al. KEGG as a reference resource for gene and protein annotation . Nucleic acids research . 2017 ; 44 ( D1 ).
26. Oksanen J , Blanchet FG , Kindt R , Legendre P , Minchin PR , O'Hara RB , et al. vegan: Community Ecology Package. R package version 2 . 0 - 10 . 2014 .
27. Love MI , Huber W , Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 . Genome Biology . 2014 ; 15 ( 12 ).
28. Boulesteix AL , Slawski M. Stability and aggregation of ranked gene lists . Briefings in bioinformatics. 2009 ; 10 ( 5 ): 556 - 68 . https://doi.org/10.1093/bib/bbp034 PMID: 19679825
29. Hofner B , Hothorn T. stabs: Stability Selection with Error Control 2014 [updated 2014 . Available from: http://cran.r-project.org/package=stabs.
30. Shilts MH , Rosas-Salazar C , Tovchigrechko A , Larkin EK , Torralba M , Akopov A , et al. Minimally invasive sampling method identifies differences in taxonomic richness of nasal microbiomes in young infants associated with mode of delivery . Microbial ecology . 2015 .
31. Schloss PD , Westcott SL , Ryabin T , Hall JR , Hartmann M , Hollister EB , et al. Introducing mothur: opensource, platform-independent, community-supported software for describing and comparing microbial communities . Applied and environmental microbiology . 2009 ; 75 ( 23 ): 7537 - 41 . https://doi.org/10.1128/ AEM.01541-09 PMID: 19801464
32. Hill MO . Diversity and evenness: a unifying notation and its consequences . Ecology . 1973 ; 54 ( 2 ): 427 - 32 .
33. Anderson MJ . A new method for non-parametric multivariate analysis of variance . Austral Ecol . 2001 ; 26 ( 1 ): 32 - 46 .
34. Oksanen J , Blanchet FG , Kindt R , Legendre P , Minchin PR , O'Hara RB , et al. vegan: Community Ecology Package 2014 [updated 2014 . Available from: http://cran.r-project.org/package=vegan.
35. Bray JR , Curtis JT . An ordination of the upland forest communities of southern Wisconsin . Ecol Monogr . 1957 ; 27 ( 4 ): 325 - 49 .
36. Hothorn T , Hornik K. exactRankTests: Exact Distributions for Rank and Permutation Tests . 2013 .
37. Benjamini Y , Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing . J R Stat Soc Series B Stat Methodol . 1995 ; 57 ( 1 ): 289 - 300 .
38. Grissom RJ , Kim JJ . Effect Sizes for Research: Univariate and Multivariate Applications, Second Edition . London: Routledge; 2012 2012 /04/23/. 453 p.
39. Meinshausen N , Bu?hlmann P. Stability selection. J R Stat Soc Series B Stat Methodol . 2010 ; 72 ( 4 ): 417 - 73 .
40. Shah RD , Samworth RJ . Variable selection with error control: another look at stability selection . J R Stat Soc Series B Stat Methodol . 2013 ; 75 ( 1 ): 55 - 80 .
Friedman J , Hastie T , Tibshirani R . Regularization paths for generalized linear models via coordinate descent . J Stat Softw . 2010 ; 33 ( 1 ): 1 - 22 . PMID: 20808728
Am J Clin Nutr . 2017 ; 106 ( 3 ): 821 - 30 . https://doi.org/10.3945/ajcn.117.152967 PMID: 28659297
Eckburg PB , Bik EM , Bernstein CN , Purdom E , Dethlefsen L , Sargent M , et al. Diversity of the human intestinal microbial flora . Science . 2005 ; 308 ( 5728 ): 1635 - 8 . https://doi.org/10.1126/science.1110591 PMID: 15831718
Gevers D , Kugathasan S , Denson LA , Vazquez-Baeza Y , Van Treuren W , Ren B , et al. The treatmentnaive microbiome in new-onset Crohn's disease . Cell Host Microbe . 2014 ; 15 ( 3 ): 382 - 92 . https://doi.
org/10 .1016/j.chom. 2014 . 02 .005 PMID: 24629344
Smith B , Bod e? S, Petersen BL , Jensen TK , Pipper C , Kloppenborg J , et al. Community analysis of bacteria colonizing intestinal tissue of neonates with necrotizing enterocolitis . BMC Microbiol . 112011 . p.
Brower-Sinning R , Zhong D , Good M , Firek B , Baker R , Sodhi CP , et al. Mucosa-associated bacterial diversity in necrotizing enterocolitis . PloS one . 2014 ; 9 ( 9 ):e105046. https://doi.org/10.1371/journal.
pone.0105046 PMID: 25203729
Proc Natl Acad Sci U S A . 2010 ; 107 ( 26 ): 11971 - 5 . https://doi.org/10.1073/pnas.1002601107 PMID: 20566857
Morrow AL , Lagomarcino AJ , Schibler KR , Taft DH , Yu Z , Wang B , et al. Early microbial and metabolomic signatures predict later onset of necrotizing enterocolitis in preterm infants . Microbiome . 2013 ; 1 ( 1 ): 13 . https://doi.org/10.1186/2049-2618-1-13 PMID: 24450576
Koenig JE , Spor A , Scalfone N , Fricker AD , Stombaugh J , Knight R , et al. Succession of microbial consortia in the developing infant gut microbiome . Proc Natl Acad Sci U S A . 2011 ; 108 Suppl 1 : 4578 - 85 .
McMurtry VE , Gupta RW , Tran L , Blanchard EEt , Penn D , Taylor CM , et al. Bacterial diversity and Clostridia abundance decrease with increasing severity of necrotizing enterocolitis . Microbiome . 2015 ; 3 : 11 .
https://doi.org/10.1186/s40168-015 -0075-8 PMID: 25810906
Guthrie SO , Gordon PV , Thomas V , Thorp JA , Peabody J , Clark RH . Necrotizing enterocolitis among neonates in the United States . Journal of perinatology: official journal of the California Perinatal Association . 2003 ; 23 ( 4 ): 278 - 85 .
52. Jacquot A , Neveu D , Aujoulat F , Mercier G , Marchandin H , Jumas-Bilak E , et al. Dynamics and clinical evolution of bacterial gut microflora in extremely premature patients . J Pediatr . 2011 ; 158 ( 3 ): 390 - 6 . https://doi.org/10.1016/j.jpeds. 2010 . 09 .007 PMID: 20961563
Subramanian S , Huq S , Yatsunenko T , Haque R , Mahfuz M , Alam MA , et al. Persistent gut microbiota immaturity in malnourished Bangladeshi children . Nature . 2014 ; 510 ( 7505 ): 417 - 21 . https://doi.org/10.
1038 /nature13421 PMID: 24896187
Moles L , Gomez M , Heilig H , Bustos G , Fuentes S , de Vos W , et al. Bacterial diversity in meconium of preterm neonates and evolution of their fecal microbiota during the first month of life . PloS one . 2013 ; 8 ( 6 ):e66986. https://doi.org/10.1371/journal.pone. 0066986 PMID: 23840569
Stewart CJ , Marrs EC , Magorrian S , Nelson A , Lanyon C , Perry JD , et al. The preterm gut microbiota: changes associated with necrotizing enterocolitis and infection . Acta paediatrica (Oslo , Norway: 1992 ).